What is the difference in using fix shake vs. fix rigid? They seem to have the same goals and just achieve them through different means. Is that a fair characterization or is there some reason to use one over the other? fix rigid seems to be objectively better and faster.
Both fixes have fundamentally different approaches.
With fix rigid you take a group of points/particles and move them as one (rigid) object. That means the fix determines the center of mass and moments of inertia for each rigid object and then gets the center of mass force and the torque around the center of mass from the pair atom/particle forces and torques.
With that, it does time integration for the center of mass and the rotation around the center of mass of each rigid object.
With fix shake (or fix rattle), atoms are moved just like they were individual atoms and then the fix solves the constraint equations iteratively. As a result the position and velocities updates are changed so that total kinetic energy and momentum are conserved by also the constraints fulfilled.
Fix shake is generally far superior because it allows a much larger timestep (up to 2fs for biomolecular systems), but it has many restrictions (it works for constraining individual hydrogen bonds and water molecules (by turning the angle constraint into a bond constraint between the two hydrogen atom).
Fix rigid in contrast requires a much shorter timestep and the integration of motions are far more sensitive and thus require more fine grained discretization, i.e. a shorter timestep (0.5fs or less for water). Fix rigid does not scale well in parallel since it requires multiple collective communications to collect and update the rigid body information across all processors. Fix rigid/small tries to address this, but only works for small molecules that are suitably shorter than the communication cutoff so that all constituent atoms are present on the individual subdomain and thus no global communication would be required.
Fix rigid is also capable to handling extended particles (spheres, ellipsoids), pair styles with dipoles and similar that result in torques and even overlapping extended particles (in that case, however, the moments of inertia have to be computed externally and provided as input).
So both have their use and applications where one is preferred over the other.
6 Likes
Awesome makes perfect sense.
One last clarifying question. If I were to use fix shake or fix rattle to constrain some bonds would those atoms then be updated using an NVE update even if I have an NVT fix acting on the entire system? I’ve tried to look through the source code but get a bit lost in the structure of LAMMPS outside of the RATTLE/SHAKE files. In those files the update is definitively NVE, but I am also allowed to use fix shake/rattle with NVT or NPT (e.g. in this example). Was hoping for some clarification for how other ensembles are handled by these algorithms. Is is just that the constrained atoms get an NVE update while the remaining particles use the other integrator or is it more complicated? I feel it is the latter as something like NoseHoover is basically just Velocity Verlet with extra degrees of freedom added but the literature nor documentation really touch on this point.
I understand you might not be an expert in this, just trying to understand the code I’m using, thanks!
This is a misleading statement. The position and velocity updates in LAMMPS always follow the Velocity Verlet scheme. The difference between fix nve and fix nvt is that fix nve does only do the updates, while fix nvt will couple a heat bath (through fictitious oscillators) to the updates according to the Nose-Hoover algorithm which modifies the forces used in the velocity update (see \zeta e.g. Nosé-Hoover thermostat page on SklogWiki - a wiki for statistical mechanics and thermodynamics).
Fix shake applies a correction to the position updates that preserves the momentum of the affected aggregates while maintaining the constraints. So, you first have an unconstrained update of the positions, and then fix shake iteratively solves the constraint equations. Thus for a water molecule, the center of mass momentum and the rotational/angular moments are preserved relative to the unconstrained update.
This situation with fix shake is a bit complicated, because the SHAKE algorithm is technically only derived for a regular Verlet time integration (where velocities are implicit). It thus has to recast the Velocity Verlet updates to that. This works ok in the middle of a run, but can be a problem in the zeroeth step (where LAMMPS recovers the forces from the positions, so it can do the first half velocity update). The RATTLE algorithm is an extension to SHAKE for Velocity Verlet.
Yes, in principle. But careful with fix rattle. It has constraints in which operation has to come first or later than cannot be satisfied in combination with fix npt.
The literature is quite accurate and detailed, but not an easy read. What LAMMPS primarily implements by default is the MTK paper, as far as I recall. This can fall back to other, older formulations when certain parameters are changed, e.g. the number of oscillators in the Nose-Hoover chain(s).
The documentation is rather old since the nvt/npt fixes exist in LAMMPS for a very long time and the documentation for those was written at a time, when it was very laborious to add mathematical expressions to it (Now it is much easier, after we converted to using Sphinx and MathJAX to build the documentation). So people had a tendency to avoid those.
In general, you have to understand that the documentation is not supposed to be an explanation of the algorithm (for that you have references to the original literature) or the details of its implementation (for that you have the source code). It is is more like an owner’s manual of a car: it shows you all the knobs and lights and handles and lists what they do, but it will not teach your how to drive.
With the Nose-Hoover code you are looking at one of the more complex and hard to figure out pieces of code in LAMMPS. This is because there is a lot of shared code between fix npt, nph, and nvt, so fix_nh.cpp/fix_nh.h was written as a base class that can do either depending on a few flags. Later, some additional features were added and now the code is rather convoluted.
3 Likes