The need for accurcy in constrained molecular dynamics simulations

In constrained molecular dynamics a system of differential algebraic equations
$$\begin{align} \dot{q} &= v, \\ M\dot{v} & = f(q) - G(q)^T \lambda, \\ g(q) &= 0 \end{align} $$
is solved. Here $(q,v)$ is the position and velocity of all particles, $M$ is a nonsingular diagonal mass matrix, $f$ represents the real forces acting on the particles, $g$ is the (holonomic) constraint function, $G$ is the Jacobian of the constraint function and $\lambda$ is a vector of Lagrange multipliers. In computational chemistry, the standard algorithm for this problem is the celebrated SHAKE algorithm:
$$\begin{align}
 v\_{k+1/2} &= v\_{k-1/2} + h M^{-1} \left( f( q\_k) - G( q\_k)^T \lambda\_k \right) \, , \\
 q\_{k+1} &= q\_k + h v\_{k + 1/2} \, , \\
 g(q\_{k+1}) &= 0
\end{align}
$$
The final equation is a nonlinear equation for the vector $\lambda\_k$ of Lagrange multipliers. It is typically solved using the nonlinear Gauss-Seidel method which converges linearly under certain mild conditions. Over the years, there have been several attempt to replace the Gauss-Seidel with quasi-Newton methods that require fewer, but more expensive iterations.


This question concerns the need for accuracy when solving the nonlinear equation
$$
g(q\_{k+1}(\lambda\_k)) = 0$$
with respect to the Lagrange multiplier $\lambda\_k$.


The SHAKE algorithm has several nice properties that all hinge on the constraint equation being solved exactly. These properties are lost when we solve the constraint equation with a tolerance that is greater than zero, and quality of the simulation tends to degrade as we increase the tolerance. In the absence of detailed knowledge of the constraint function $g$, the only really "safe" value for the tolerance is the unit roundoff used by the computer. In particular, we lose the second order accuracy of the method and we loose the ability of compute reliable error estimates using Richardson extrapolation when the tolerance is not tiny.


The SHAKE algorithm has excellent conservation of energy when we use a tiny tolerance, but the total mechanical energy can increase or decrease when the tolerance is not tiny. Moreover, if the tolerance is too large, then it is possible to venture so far from the constraint surface, that the constraint solver is unable to find a solution and the simulation fails. The statements in this paragraph can easily be demonstrated in the case a simple pendulum, i.e., a single constraint.


Mathematically, errors associated with the Lagrange multiplier are equivalent to an external force field. Hence, if the objective is to simulate an isolated system, then we can hardly accept anything but the smallest possible tolerance that can be delivered by the hardware.


I would like to know a few numbers that can be reliably computed from constrained simulations that use a large tolerance when solving the constraint equations. I am primarily interested in numbers that have real life applications such as the development of medicine and vaccines.


Similarly, I would like to know a few numbers that cannot be trusted unless the underlying simulation are executed using a tiny tolerance. Again, I am primarily interested in numbers that have real life applications.


I have a fair understanding of differential equations, numerical analysis, numerical linear algebra and parallel computing, i.e., the skills needed to produce numbers reliably and quickly. On the other hand, I know very little about chemistry, so I may need more basic information than is customary here.


