Fix RATTLE constraints for MD

The current implementation of RATTLE for rigid bonds (FixBondLengths class in constraints.py) uses the old atoms positions for calculating the undetermined multipliers to apply the constraints. This means that in a MD simulation where one uses RATTLE to fix certain bonds, the constraints at each step will be applied based on the positions at the previous step, which could lead to propagation of errors on the bondlengths. For this reason the convergence criteria should be kept as low as 1e-10 - 1-13 to avoid significant drifts in the simulation box. Following the original formulation of RATTLE for MD simulations with the Velocity Verlet integrator (https://doi.org/10.1016/0021-9991(83)90014-1), I have implemented a new version of FixBondLengths, which I named FixBondLengthsMD, that calculates the undetermined multipliers based on a fixed bondlength value. The value of the bondlength one wants to fix has to be specified together with the pair of atoms participating in the bond:r01 = 0.95 c = FixBondLengthMD(0, 1, r01) atoms.set_constraint(c)I have tested the performance of the new implementation vs. the old (current) implementation on a box of 112 TIP4P water molecules previously equilibrated at 300 K, employing Velocity Verlet for the propagation. In both cases I used a tolerance of 1e-5 to keep the water moleculs "rigid". Below I plotted the bondlengths inside a single water molecules in the box vs. simulation time:
Using the old implementation with this convergence criteria causes the bondlengths and energy to drift significantly until the box is blown apart:
The new implementation, on the contrary, seems to give results that are stable enough on a 200 ps time scale. I have already modified the rattle.py test to exercise the new implementation.