Bond Constraints with SHAKE and V-rescale thermostat

User discussions
1

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?

Thanks,

Michael DeLyser

2

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.