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!