Bond Constraints with SHAKE and V-rescale thermostat

GROMACS version: 2019.4
GROMACS modification: No

I am running a large CG simulation with a martini-like force field. Some of the bonds are constrained to a set distance, others are allowed to fluctuate. I am running in the NPT ensemble at 1 bar (Parrinello-Rahman barostat) and 300K. Time step is 10 fs, thermostat coupling constant is 2 ps, barostat coupling constant is 5 ps.

When I use the v-rescale thermostat and SHAKE to constrain the bonds, the average temperature comes out to around 275K. Further, when I use gmx energy to output the temperature, the output is:

Energy Average Err.Est. RMSD Tot-Drift
Temperature 5519.72 1600 31223 11013.5 (K)

yet the temperature.xvg file created by gmx energy gives values around 275K for the temperature at each output step.

I tried using Nose-Hoover & SHAKE and I get the right temperature. I tried using v-rescale & LINCS and I get the right temperature. I tried changing the constrained bonds to harmonic bonds and I get the right temperature.

Is/was there a known issue when combining the SHAKE algorithm with the v-rescale thermostat?


Michael DeLyser

My guess is that your SHAKE parameters are not sufficiently accurate, or maybe cannot be sufficiently accurate, and the constraint errors pump energy into the system faster than the v-rescale thermostat can pump out with a coupling time of 2 fs. If my guess is correct, you should not use SHAKE, or at least not with the SHAKE parameters you are using.