c-rescale doesn't work in NPH
Summary
C-rescale uses the reference temperature in its stochastic term. When using MD and no temperature coupling algorithm, grompp
silently sets the reference temperature to 0, independently of the mdp value of ref-t
. C-rescale then, just as silently, uses this reference temperature in its stochastic term.
This leads to the stochastic term cancelling out, and c-rescale essentially retrieving Berendsen barostat behavior. Under NPH, c-rescale is therefore having the same non-physical behavior as Berendsen, without any notice to the user.
Edit: As mentioned by @GiovanniBussi below, c-rescale doesn't become identical to Berendsen, as c-rescale still does velocity scaling which is not present in Berendsen.
This affects all versions of GROMACS in which c-rescale is present.
Possible fixes
I think we can
- forbid the combination (i.e. throw an error if c-rescale barostat is chosen with ref-t == 0), which would mean that c-rescale can't be used for NPH calculations.
- add an additional mdp value for the barostat reference temperature. That would likely require some sanity checks when ref-t != 0 (what happens if a user choses different temperatures for thermostat and barostat??).
- change the behavior of grompp, and avoid it setting the reference temperature to 0. I think that would be the cleanest solution, but would require careful checking that no other code is depending on this value being zeroed in NVE / NPH calculations.
Happy to hear your thoughts, @GiovanniBussi and @berkhess.