In molecular dynamics, penalized overdamped Langevin dynamics are used to model the motion of a set of particles that follow constraints up to a parameter ε. The most used schemes for simulating these dynamics are the Euler integrator in R d and the constrained Euler integrator. Both have weak order one of accuracy, but work properly only in specific regimes depending on the size of the parameter ε. We propose in this paper a new consistent method with an accuracy independent of ε for solving penalized dynamics on a manifold of any dimension. Moreover, this method converges to the constrained Euler scheme when ε goes to zero. The numerical experiments confirm the theoretical findings, in the context of weak convergence and for the invariant measure, on a torus and on the orthogonal group in high dimension and high codimension.