Shake and langevin

LAMMPS LAMMPS General Discussion
1

Hello lamps user

I was wandering what is the correct order of fixes shake and langevin when used in a molecular dynamics simulations.

I understand that langevin adds forces to individual atoms during the simulations

And in the shake documentation it is said that shake should appear after fixes that add forces.

I am asking because the documentation is very specific. But there is no mention to langevin in the shake documentation and also not mention of shake in the langevin documentation. I would have expected the langevin documentation to metion how it should be used with constraints.

Thank for any clarification on this matters

Matías.

2

There are a lot of fixes that add or modify forces and there are new fixes added all the time. It would be a nightmare to add those all to the fix shake documentation and vice versa.

The LAMMPS way to handle this kind of situation can be seen by looking at the Shake::init() function where it is checked that any fix changing the box dimensions must come after fix shake. Fixes that change box dimensions have a member variable “box_change” set to 1 (default is 0). So the equivalent procedure would require to add a “force_change” member to the Fix class in fix.h and initialized it to 0 in fix.cpp and then add “force_change = 1;” to the constructor of any fix that changes forces and finally a check testing for the ordering to fix_shake.cpp.

Before spending time on implementing this and posting a pull request, it may be worth discussing with @sjplimp, and the best chance for that would be submitting a feature request issue on Issues · lammps/lammps · GitHub and assigning it to @sjplimp or just mention him (and me).

3

This is another opportunity to quote George Box:

It is inappropriate to be concerned about mice when there are tigers abroad.

I am sure there is a “correct” order in which thermostatting and rigid constraints should be applied. But both this order and the order’s effect will vary tremendously from system to system, especially on its phase and response characteristics (which are often the very properties to be measured, so what then?). This paper lists forty-two different orderings for the various manipulations in an NVT-SHAKE scheme.

So I would argue that it is impractical for LAMMPS to try and police these relative orderings, other than simply offering a common convention as described in the fix shake documentation. We could simply make additions (in bold) to the current text:

If you define fixes (e.g. fix efield, and including thermostats like fix nvt and fix langevin) that add additional force to the atoms after fix shake operates, then this fix will not take them into account and the time integration will typically not satisfy the SHAKE constraints.

and then just count on scientists to do their best. The impact of force field approximations, floating point truncation, PPPM inaccuracy, and others is likely to far outweigh inaccuracies from thermostat-constraint non-commutativity. I also suspect that most choices of timestep are either conservative enough that the effects will be limited, or large enough that the entire simulation destabilizes, and thus one would have to be unlucky (and careless) to land a timestep that is both small enough to return a largely stable simulation and large enough that flipping thermostat and constraints would give different answers.

Finally – these days I am seeing so many peer-reviewed papers sans their input scripts, and it drives me to despair. If I can’t see the script that was used to generate data, it’s really immaterial whether fix langevin came before or after fix shake, because neither choice is testable or reproducible relative to the paper data. It’s a little like running OSHA checks on the sturdiness of the deck chairs on the Titanic. The greatest need now is for strong community awareness and support for the most basic of standards – no peer review acceptance without at least a pretense of putting a script online – especially when it’s not long before someone at Google starts drafting their next Nature paper on “autonomous high-throughput non-equilibrium molecular dynamics to discover novel rheological regimes”, where they proudly tour a menagerie of LAMMPS simulations that meticulously sample complicated chemicals modelled in thoroughly unphysical ensembles. I’m personally too upset by this sort of thing to care so much any more about the finer points of O(t^2) errors from non-commutativity of integrator substeps.