Systematic slight discrepancy of the temperature value in NPT simulation

User discussions
1

GROMACS version: 2020.3
GROMACS modification: No

Using Gromacs 2020.3 I noticed this anomaly in an NPT simulation of an aqueous solution (OPLS-AA force field; box size about 9x9x9 nm^3); mdrun runs on a Linux cluster using CPU + GPU). The average value of the temperature (controlled by the v-rescale thermostat) is slightly different from the value of the mdp directive (for example 302.77 K instead of 303 K, respectively). This systematic error is reproducible.
I found this behavior rather strange, since in past simulations the same thermostat setting gave the desired temperature value (in K) within the second or even the third decimal place.

These are the mdp directives I used:

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
Tcoupl = v-rescale
nsttcouple = -1
nh-chain-length = 10
print-nose-hoover-chain-variables = no
; Groups to couple separately
tc_grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.1
ref_t = 303
; pressure coupling
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
nstpcouple = -1
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p = 2
compressibility = 4.4e-5
ref-p = 0.987
; Scaling of reference coordinates, No, All or COM
refcoord-scaling = all

and these the output from gmx energy on four subsequent trajectory intervals each of 100 ns

Energy Average Err.Est. RMSD Tot-Drift

Temperature 302.77 0.0014 1.14204 0.000128223 (K)

Temperature 302.768 0.0039 1.14271 0.00180873 (K)

Temperature 302.766 0.0025 1.14083 -0.011248 (K)

Temperature 302.764 0.0013 1.14201 -0.00211153 (K)

Obviously for my purposes to have 303 or 302.77 K does not make difference, however I would like to be sure that this is not the clue of possible problems in the simulation setup.

Did anyone detect a similar behaviour?

Thanks in advance for any comment or suggestion

Ester

2

Did you use exactly the same mdp options in both cases?
Which integrator did you use?
Could you compare the drift in the conserved energy quantity?

3

Thank you for the answer.

Yes, all mdp options are the same and we used the leap-frog algorithm.
However, I realized that there is a difference in the algorithm used to model rigid water: in one case it is lincs, in the other one it is settles. With lincs the average temperature is slightly lower than the ref-t value.
E.g., for ref-t=293 in an NPT run of 30 ns we obtain
2.92780e+02 K (lincs)
2.92996e+02 K (settles)

4

That explains it then. SETTLE is exact, up to rounding errors, whereas LINCS is only exact in the limit of high order and iterations, but then it also has more rounding errors.

5

Thank you for these explanations.