Dear LAMMPS users,
I am running simulations of rigid bodies using Langevin dynamics, and I need to probe the high-friction limit (overdamped dynamics – sometimes also called “Brownian” dynamics). In practice, I want the acceleration (inertia) term to be negligible compared to the other terms in the Langevin equation of motion.
I am currently using fix/rigd/nve + molecule + langevin to simulate my system.
However, the Langevin thermostat as implemented in LAMMPS (Schneider and Stoll, Phys Rev B, 17, 1302, 1978), is expected to break down in the high-friction limit, as discussed in Bussi and Parrinello, Phys. Rev. E 75, 056707, 2007.
I thought therefore that I could use the Langevin thermostat of Bussi and Parrinello, which is implemented in LAMMPS (fix temp/csld).
Unfortunately, unless I’m mistaken, it seems that this thermostat cannot be used in conjunction with fix/rigid, which only supports the “normal” Langevin thermostat and the Nose-Hoover thermostat.
Is there a way to use the Bussi-Parrinello thermostat in conjunction with fix/rigid?
If this is not possible, what is the correct way to implement overdamped Langevin dynamics in LAMMPS?
Best regards,
Valerio
Hi Valerio, indeed there are several Langevin integration schemes out there, and although the more modern ones are designed to handle the high-friction limit correctly, not all of them have been implemented in LAMMPS in a way that can handle transparently both fix nve (i.e. atoms) and fix rigid (i.e. whole bodies).
The fix langevin command also provides the Grømbech-Jensen and Farago integration scheme (Mol Phys, 2013), which was developed without the the low-friction approximation of Schneider and Stoll:
https://doi.org/10.1080/00268976.2012.760055
Besides the very simple examples included in that paper, we also found it to be suitable for high-friction simulations on a realistic system (polymer chains):
https://doi.org/10.1080/00268976.2019.1649493
If you choose to use the G-JF scheme, it may be appropriate to use the half-step velocity estimator via the “gjf vhalf” keyword:
https://lammps.sandia.gov/doc/fix_langevin.html
https://doi.org/10.1080/00268976.2019.1570369
This alternative formulation of the velocity does not change the trajectory in coordinate space (which will follow the canonical ensemble), but reports more correct kinetic energies (which would be otherwise systematically biased as a result of the discrete time-stepping).
Giacomo
Dear Giacomo,
Thank you very much, this is really useful.
I didn’t know this version of the Langevin thermostat had been implemented in LAMMPS.
Also, for some reason, I thought that to use the Langevin thermostat with rigid bodies one had to run a command in the form
fix 1 all rigid/nve/small molecule langevin {mytemp} {mytemp} {damp} {lseed}
I see now that another possibility is to run, e.g.,
fix fix_rigid all rigid/nve/small molecule
fix fix_gjf all langevin {mytemp} {mytemp} {damp} {lseed} (gjf no)
where “gjf no” (or also gjf half etc.) can be included if one wishes to use the G-JF integration scheme.
May I ask, what’s the difference between the first and the second way of running a Langevin thermostat with rigid bodies?
Best regards,
Valerio
Hi Valerio, about your last question on what is the difference between using the Langevin scheme integrated within fix rigid vs. the one provided by fix langevin. One thing that stands out is that in the latter you only get what is implemented, i.e. Schneider-Stoll with a flatly-distributed noise term (as opposed to G-JF with Gaussian noise). So in general for high friction, fix langevin may be more suitable for coordinate sampling. Also, the half step formulation of the velocity in G-JF is demonstrated to be more accurate for atoms, but I don’t know if it has even been tested with rigid bodies (with the due corrections for the constrained degrees of freedom). (At all times, do keep in mind that choosing full- or half-step velocity estimators doesn’t change the coordinate trajectory).
As for the fix rigid side, I completely certain from only the user doc, but if you look in fix_rigid*.cpp at the comment of apply_langevin_thermostat(), it says that in contrast to fix langevin, the Langevin scheme embedded in fix rigid will also account for the stochastic forces in the center-of-mass forces and torques.
I am copying the people who made the most recent relevant changes to the fix langevin and fix rigid code, who may want to pitch in.
Giacomo
The langevin keyword in the fix rigid commands is referring
to Langevin terms being added to the translational
and rotatational DOFs of the entire bodies. That has
only been implemented for the vanilla (original) formulation
of Langevin thermosttatting. The more sophisticated
variants (gjf, csld, etc) know nothing about rigid bodies.
They operate on point particles only, just 3 DOFs.
If you use them in a rigid body model, you should normally
just apply them to the point-particle solvent (or polymer
chain atom, etc). If you also apply them to the rigid
body atoms, those forces will be added, but they
will just become part of the total force and torque on
the bodies. Not really the same as thermostatting.
Steve