RATTLE Implementation

[ejmeitz]

Implemented iterative scheme for RATTLE....completely untested but I have a many clarifying questions/thoughts.

• Its unclear to me if these implementations work with all integrators. The derivations all use VV. I don't think other MD codes limit you to this but I don't really see why these methods are generally applicable (they definitely shouldn't work for NVT & NPT out of the box)
  • Related to this I'm not sure if RATTLE should be called after half-step updates or before or on both. I more or less get the math but not how these constraints should be integrated into a real MD code.
• I think the whole apply_constraints! interface will need to change as RATTLE and SHAKE have to be applied at different points in the integration process. That said not sure how to fix that at the moment in a pretty way.
• What does your SHAKE implementation follow? I used this for RATTLE which has the SHAKE steps in it also. I've never seen that quadratic formula to find g in any paper. Granted every paper I read does something totally different.
• LAMMPS docs have me a little confused. They say they solve RATTLE analytically (and the source matches that statement) but they also said they follow the Andersen paper which solves it iteratively. I think analytically would be better but not really clear how that is done.

[codecov]

Codecov Report

Attention: Patch coverage is 92.45283% with 16 lines in your changes are missing coverage. Please review.

Project coverage is 72.44%. Comparing base (744f11b) to head (f15637d).
Report is 2 commits behind head on master.

Files	Patch %	Lines
src/types.jl	47.82%	12 Missing ⚠️
src/constraints/constraints.jl	97.18%	2 Missing ⚠️
src/constraints/shake.jl	97.46%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #132      +/-   ##
==========================================
+ Coverage   72.18%   72.44%   +0.26%     
==========================================
  Files          35       36       +1     
  Lines        5278     5444     +166     
==========================================
+ Hits         3810     3944     +134     
- Misses       1468     1500      +32     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jgreener64]

I more or less get the math but not how these constraints should be integrated into a real MD code.

I agree this is not clear and much of implementing this correctly will be calling the constraint functions at the correct point. See for instance OpenMM, which has examples of where to call the constraints for Verlet and Langevin dynamics. There is also a Julia RATTLE implementation here but I don't know how good it is.

I think the whole apply_constraints! interface will need to change as RATTLE and SHAKE have to be applied at different points in the integration process.

Yes I think we might want to switch to constrain_coords! and constrain_velocities! and call them at the appropriate points. Whether this transfers to other constraint algorithms remains to be seen, but it should be enough for SHAKE/RATTLE.

What does your SHAKE implementation follow?

I think it follows LAMMPS, you might get more from the original PRs (#82 and #84).

LAMMPS docs have me a little confused.

Look in the source code if possible?

[ejmeitz]

Just some notes from LAMMPS SHAKE source:

• Unconstrainted position update "assumes NVE update, seems to be accurate enough for NVT,NPT,NPH as well" Line 1685
• Shake bond constraint line 1774 (our existing code is basically a copy of this)
• 2 & 3 bond clusters lines 1877 & 2051
• 2 bonds 1 angle line 2303

[ejmeitz]

Useful threads from LAMMPS

• https://matsci.org/t/fix-shake-a-hybrid-method-between-shake-and-rattle/20507/25
• https://matsci.org/t/fix-shake-vs-fix-rigid/49193

Things to remember:

• Only works when positions and velocities use the VV update step. So something like NoseHoover is still safe, with NPT the order matters a lot
• Dont use during minimization
• No 180 deg molecules

[ejmeitz]

• The replica system constructor will not parse the constraints correctly as it does not actually call the System constructor as a function rather it specifies all of the types explicitly.

[ejmeitz]

Anyone have any clue why when I simulate hydrogen gas the rotational DoF suddenly freeze? Basically the first 200 timesteps look normal then rotations completely stop and the molecules just move around the box constrained interacting as if they cannot rotate. When I simulate the same thing in LAMMPS the H2 molecules spin like crazy.

I'm not really convinced this is an artifact of my SHAKE implementation but I don't know what else could cause it.

[jgreener64]

Could you give code to reproduce the issue?

It might also be worth comparing to the same simulation without constraints, i.e. with a strong harmonic bond and a small time step. If that simulation doesn't show the issue then it is probably something to do with the constraints.

[ejmeitz]

This is the input I am using. I uploaded the xyz file for this with 1 molecule (where you can also see what I am talking about). In the 1 molecule case the two atoms do not interact with each other so this rotation is just from the initial velocity. It might look like the molecule does not translate from the initial velocity but I think they do and its just small. The weirder thing to me is that the particles stop rotating and moving. After looking at all the output I noticed that the velocity vectors eventually end up pointing towards the COM of the molecule effectively canceling each other out. Why this is the case I have no idea. I tested the exact same system in LAMMPS with SHAKE and this does not happen. I plan to test the original implementation of SHAKE in Molly (I have only changed how it is initialized not the algorithm itself) to see if this phenomena occurs as well and also the harmonic bonds.

Edit: Tried with just Verlet algorithm as well and the same issue pops up.
Edit: Tried with version of SHAKE in Molly atm and the same thing happens. In the code from the tests the particles are much closer so they interact and keep rotating, but if I change it to have the parameters below then all rotationg stops because the gas is so diffuse. I'm not convinced this is the correct behavior.

pos.zip

r_cut = 8.5u"Å"
temp = 300.0u"K"

#Initialize atoms (0.0252259036145 molecules / nm^3 typical for H2 gas at STP)
n_atoms_half = 200
atom_mass = 1.0u"u"
atoms = [Atom(index = i, mass=atom_mass, σ=2.928u"Å", ϵ=0.074u"kcal* mol^-1") for i in 1:n_atoms_half]
max_coord = 200.0u"Å" # nm
coords = [max_coord .* rand(SVector{3}) for i in 1:n_atoms_half]

#Add bonded atoms
bond_length = 0.74u"Å" #hydrogen bond length
constraints = []
for j in range(1, n_atoms_half)
    push!(atoms, Atom(index = j + n_atoms_half, mass = atom_mass, σ=2.928u"Å", ϵ=0.074u"kcal* mol^-1"))
    push!(coords, coords[j] .+ SVector(bond_length,0.0u"Å",0.0u"Å"))
    push!(constraints, DistanceConstraint(SVector(j, j+n_atoms_half), bond_length))
end

constraint_algorithm = SHAKE(similar(coords))

velocities = [random_velocity(atom_mass, temp) for i in 1:length(atoms)]
# velocities = [[0.0u"nm/ps", 0.0u"nm/ps", 0.0u"nm/ps"] for i in 1:length(atoms)]
boundary = CubicBoundary(max_coord)

neighbor_finder = DistanceNeighborFinder(eligible = trues(length(atoms),length(atoms)), dist_cutoff = 1.5*r_cut)

sys = System(
        atoms = atoms,
        coords = coords,
        boundary = boundary,
        velocities=velocities,
        pairwise_inters=(
            LennardJones(
                cutoff = ShiftedPotentialCutoff(r_cut),
                use_neighbors = true,
                energy_units = u"kcal * mol^-1",
                force_units = u"kcal * mol^-1 * Å^-1"
                ),
            ),
        neighbor_finder = neighbor_finder,
        constraints = constraints,
        constraint_algorithm = constraint_algorithm,
        loggers=( 
            tot_eng = TotalEnergyLogger(1),
            coords_out = CoordinateLogger(1),
            vels_out = VelocityLogger(1)
            ),
        energy_units = u"kcal * mol^-1",
        force_units = u"kcal * mol^-1 * Å^-1"
        )


simulator = VelocityVerlet(dt = 0.002u"ps")
simulate!(sys, simulator, 100_000, n_threads = 1)

[ejmeitz]

I've found the issue, our original implementation of SHAKE was only modifying the positions when it needed to be modifying the forces and it was doing so in the incorrect place anyways.

There is still the units issue I messaged you on Slack Joe, but for documentation I will explain here as well. In the working version I modify the accelerations which have some bizarre units kcal Å^-1 mol^-1 u^-1 but the acceleration I calculate inside of SHAKE has units of Å/ps^2. Converting from one to the other just requires a factor of 418.6 and to realize that 1 u is 1 g/mol. I did this conversion manually to check that everything works (it appears to) but Unitful is unable to handle this conversion. Likely because the original acceleration could not cancel the molar part.

[jgreener64]

I am having a look at this but it is a large and complex change and I am not so familiar with the algorithms. There are a couple of thoughts I have so far:

• The approach of modifying the accelerations rather than the coordinates/velocities may be limiting. I guess this is how LAMMPS does it but it probably means more maths to get it working for other simulators. Does this make it harder to support things like half-step velocities and user-defined simulators? We would want to support Langevin with this PR, for example. An approach that applies constraints to coordinates or velocities with two functions would allow those functions to be called at the end of the integration step, or at least in a more understandable way. It might be worth seeing what other software does, for example the OpenMM custom integrator interface which is compatible with constraints.

• The docs describe compatibility with SHAKE and RATTLE per simulator but I'm not sure this is the way it should be presented. As I understand it SHAKE is sufficient when velocities are not directly used, and if velocities are used then both SHAKE and RATTLE are required. So the user doesn't need to know which simulator works with which, but rather whether that simulator supports correct constraints or not. From my understanding all do here all apart from Langevin and LangevinSplitting (MetropolisMonteCarlo and SteepestDescentMinimizer wouldn't tend to use constraints).

• It would be nice to see the performance of constraints. A test of the implementation and of performance would be to constrain all bonds to hydrogen in a simulation of a protein similar to that in the docs.

Let me know what you think. There's a lot of good work in here but I want to make sure it's a flexible and accurate approach to constraints. We shouldn't require too much caution on behalf of the user when using constraints, which basically all biomolecular simulations will do.

[ejmeitz]

Tests passed locally, just to summarize this PR.

New Features:

• Rattle implemented for VelocityVerlet
• Velocity corrections added to all integrators that use SHAKE
• Unified RATTLE/SHAKE interface into one class
• Setup constraint framework so that it will be easy to add more complicated algorithms
• A function to calculate the # of DoF lost by using SHAKE, this correction is applied when calculating tempearture

Things to double check before merging:

• Energy conservation with constraints, maybe do a parity check with OpenMM for Langevin
• Check DoF function, I am not sure it is generic

What I was unable to achieve:

• Constraints for NoseHoover, Overdamped Langevin, or LangevinSplitting. I think each of these should be possible to implement they just are not now.
• Angle constraints (should be a relatively easy new PR)
• Faster version of SHAKE that can be parallelized (also left to another PR)

[jgreener64]

Nice. I'll try and review this properly over the next few days.

From an initial pass, could you remove all commented-out code and anything that is just setting up for later.

Also the field in System can remain as constraints, and I don't think the eligible matrix should be updated automatically for constraints (we leave it to the user to setup elsewhere for clarity).

[ejmeitz]

Why not do the eligible matrix automatically? It is something that HAS to be done when you use constraints. I don't really see a scenario where having that be automatic would break things?

[jgreener64]

For the standard molecular case pairwise interactions are excluded for constrained atoms, but there are more general simulation cases where people may want other behaviour.

[ejmeitz]

Can I add a flag for this? I think by default we should do this, it is what I would expect the code to do and also a far more common use case imo. (then again Ive never simulated macromolecules).

I guess I dont see how even with many body interactions you wouldnt want this behavior. Arent intra molecular forces usually off?

[jgreener64]

This is starting to take shape, well done. I left lots of comments. You will also need to merge/rebase recent changes.

I think for now we should not change the neighbour list automatically. For biomolecular simulation it should be done during system setup from a file, as it is for bonds currently. That can be a future change. The user would specify which type of constraints they want, e.g. all bonds to H, then these would be added and the bonds removed.

[ejmeitz]

I removed the automatic update to the eligible list, but I've left the function to disable the interactions given a list of constraints. I also updated the docs code to reflect that this should be done with the code I put there.

[ejmeitz]

• Note to self, add test for Overdamped Langevin

[jgreener64]

Okay this looks good, just some minor changes before merge.

[jgreener64]

Nice one, thanks 🎉