Inaccurate thermostatting when simulating MD water

GROMACS version:2024.6
GROMACS modification: No

I use the TIP4P/Ice water model frequently, and I noticed in passing that the temperature output by gmx energy does not quite average out to the ref-t value when averaged over time. Example graph (block averaging done by taking 100ps wide blocks):

The block-averaged <T> came out to be 297K, rather than the ref-t = 300K, and the block averaged s.e.m. was <1K. I’m not confident that this is correct behaviour, even if the measured temperature is 1% off.

What I found most curious while further investigating myself is that this does not happen for TIP4P-Ew water (using the itps from the gromacs share directory):

Everything else was kept identical between the two runs, including the mdp and the starting gro file. Only the topol.top file was modified. I am using the V-rescale thermostat with a 1 ps time constant. I don’t think this is likely to be the problem, else it should’ve shown similar incorrect behaviour for both models.

This almost certainly appears to be a topology problem at this point, so I’ll paste my topology files here and hope someone can spot an error that I missed out. I have cross-checked the numerical values of the LJ parameters, the charges, and the bond lengths numerous times, including a few moments ago. I’m at a loss as to where the bug is.

File content that may be relevant:

npt.mdp:

title                   = 
define                  = 
; Run parameters
integrator              = md          ; leap-frog integrator
nsteps                  = 5000000     ; 5,000,000 = 10 ns
dt                      = 0.002       ; 2 fs

; Output control
nstxout                 = 500        ; save coordinates every 1 ps
nstvout                 = 500        ; save velocities every  1 ps
nstenergy               = 500        ; save energies every 1 ps
nstlog                  = 500        ; update log file every 1 ps

; Center of mass motion removal
comm-mode               = Linear       ; remove COM linear momentum

; Bond parameters
continuation            = no        
constraint-algorithm    = lincs   
constraints             = h-bonds       
lincs-iter              = 1          ; accuracy of LINCS
lincs-order             = 4          ; also related to accuracy

; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
verlet-buffer-tolerance = 0.001    ; GROMACS: To ensure the error is below 0.1%, decrease verlet-buffer-tolerance to
ns-type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; largely irrelevant with Verlet scheme
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme

; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme-order               = 4         ; cubic interpolation
fourierspacing          = 0.12      ; grid spacing for FFT

; Velocity generation
gen-vel                 = yes
gen-temp                = 300
gen-seed                = 666

; Temperature coupling is on
tcoupl                  = V-rescale         ; stochastic velocity rescaling
tc-grps                 = System
tau-t                   = 1.0               ; time constant, in ps
ref-t                   = 300               ; reference temperature, one for each group, in K

; Pressure coupling is on
nstpcouple              = 25               ; update interval for pressure coupling
pcoupl                  = C-rescale                             
pcoupltype              = isotropic                           
tau-p                   = 1.0               ; time constant, in ps
ref-p                   = 1.0
compressibility         = 4.5e-5
refcoord-scaling        = com

; Periodic boundary conditions
pbc                     = xyz           ; 3-D PBC

topol.top:

[ defaults ]
; nbfunc      comb-rule       gen-pairs       fudgeLJ fudgeQQ
  1             2               yes              1.0     1.0
 
[ atomtypes ]
; name      mass      charge    ptype   sigma         epsilon
  IW     0             0.000       D    0.0           0.0
  OWT4   15.99940      0.000       A    0.31668       0.88211
  HW     1.00800       0.000       A    0.0           0.0
 
#include "ff_liq.itp"
 
[ system ]
 System of TIP4P/Ice liquid water
 
[ molecules ]
 SOL  3000

.itp:

; TIP4P/ice model by Abascal, Sanz, Fernandez, Vega 2005 (https://pubs.aip.org/aip/jcp/article/122/23/234511/901130/A-potential-model-for-the-study-of-ices-and)

[ moleculetype ]
 ; name nrexcl
 SOL  2
 
[ atoms ]
 ; nr    type resnr residu atom cgnr  charge
   1     OWT4  1     SOL  OW   1     0         15.994
   2     HW    1     SOL  HW1  1     0.5897    1.008
   3     HW    1     SOL  HW2  1     0.5897    1.008
   4     IW    1     SOL  MW   1    -1.1794    0.0
 
[ constraints ]
 ;i j funct doh  dhh
 1       2       1       0.09572
 1       3       1       0.09572
 2       3       1       0.15139
 
[ exclusions ]
 1       2       3       4
 2       1       3       4
 3       1       2       4
 4       1       2       3
 
 ; The position of the dummy is computed as follows:
 ;
 ;               O
 ;
 ;                D
 ;
 ;       H               H
 ;
 ; const = distance (OD) / [ cos (angle(DOH))    * distance (OH) ]
 ;         0.015 nm      / [ cos (52.26 deg)     * 0.09572 nm    ]
 ; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1)
 
[ dummies3 ]
 ; Dummy from                    funct   a               b
 4       1       2       3       1       0.13458         0.13458  

Ref for TIP4P/Ice model:
J. Chem. Phys. 122, 234511 (2005)

TIA for any assistance!

Update:

It appears to be the same issue as: Systematic slight discrepancy of the temperature value in NPT simulation. Modifying my .itp files with [ settles ] instead of [ constraints ], appropriately redefined, appears to have fixed the problem. I’m still curious what other effects it might have had in the overall range of system properties when using LINCS versus SETTLE for rigid water.