  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 $\varepsilon$. The most used schemes for simulating these dynamics
are the Euler integrator in $\mathbb{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 $\varepsilon$. We
propose in this paper a new consistent method with an accuracy independent of
$\varepsilon$ for solving penalized dynamics on a manifold of any dimension.
Moreover, this method converges to the constrained Euler scheme when
$\varepsilon$ 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.
