RATTLE constraint for velocities

[chemicalfiend]

I implemented the RATTLE constraint which constrains velocities along constraints to be 0. While the test works for one usage of apply_constraints works, it is not converging for multiple steps (refer the test here : https://github.com/chemicalfiend/md-constraints/blob/main/test/rattle-test.jl ; the first 2 tests pass, but the simulate! command keeps running forever)

[codecov[bot]]

Codecov Report

Merging #97 (d389340) into master (c684c87) will increase coverage by 0.08%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #97      +/-   ##
==========================================
+ Coverage   83.87%   83.96%   +0.08%     
==========================================
  Files          32       32              
  Lines        3795     3816      +21     
==========================================
+ Hits         3183     3204      +21     
  Misses        612      612              

Impacted Files	Coverage Δ
src/simulators.jl	95.34% <ø> (ø)
src/constraints.jl	96.49% <100.00%> (+2.04%)	:arrow_up:
src/types.jl	73.36% <0.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[jgreener64]

Sounds like maximum(vpset) never gets low enough in the case where it fails. If you print out maximum(vpset) each iteration what does it give you? 10^-15 may be too strict a requirement.

Also, do you have a link to the implementation you are basing this on?

[chemicalfiend]

1e-15 was too strict a requirement, indeed. I think because the units of v . r is nm^2 / ps we can only expect it to converge to about an order of O(10^-3) . I tried printing out maximum(vpset) and no matter what I did it seemed to only hit roughly 1e-2 nm^2/ps for each of the the test cases I was working on, so I relaxed the tolerance bit by bit until it worked. I'm concerned it might be a little too relaxed now...

As for references, this has the general implementation for pairs; I used the iteration scheme prescribed here ; I also referred the paper here directly .

[jgreener64]

Thanks for the links. Yes I think that is too relaxed a requirement.

Is the line r01 = vector(old_coords[i1], old_coords[i0], sys.boundary) and the corresponding line after correct? The vector function arguments would be i0 then i1 if you wanted the vector from i0 to i1. Here they are the other way around. I notice it is the same order as for SHAKE but there you are mainly taking norms or doing other operations that would cancel out an error.