Calculating temperature from the velocity information of trr file

GROMACS version: version 2021.1
GROMACS modification: No

Hello. I am a beginner of GROMACS.
I have a question about how to calculate the temperature from the velocity information in the trr file.

I would like to analyze the energy and temperature of each type of molecule in a mixed system in which multiple types of molecules exist.

I thought this verification could be achieved by the calculation with equation (8) described in the link below (Molecular Dynamics — GROMACS 2021.1 documentation. The number of degrees of freedom is Ndf = 3N-3, and the mass of an atom is calculated by multiplying the atomic mass unit (1.66054 × 10 ^ -27kg) by the atomic number. Also, since the atomic velocity information after simulation could be obtained by generating a trr file with the gmx trjconv command, I considered performing this operation for each molecular species.

To verify the effectiveness of this method, I first tried to calculate the temperature of a system containing only water.

First, 884 molecules of water are added to a cubic cell with a side of 3 nm and energy minimization was conducted. Second, equilibration was performed for 500 picoseconds under the nvt condition (reference temperature: 300 K) and the nve condition respectively. After Equilibration was done, finally a production simulation of 500 picoseconds was performed under nve conditions.

After the calculation, the average temperature output by the gmx energy command was 297.3 K. On the other hand, the calculated average temperature based on the above equation (8) using the velocity information of trr file was 196.2 K. This deviates by 100 K or more from the data output by the gmx energy command.

I would be grateful if you could give me some advice on why such a divergence occurs and possible causes. The calculation conditions used for the nve simulation are as follows.

; Run parameters
integrator = md-vv-avek ;
nsteps = 250000 ; 0.002 * 250000 = 500 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ;
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 2 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
verlet-buffer-tolerance = 0.0000033 ; Changed from default (0.005) according to notes
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; largely irrelevant with Verlet
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.16 ; grid spacing for FFT
; Temperature coupling is off
tcoupl = no ; nose-hoover thermostat
;tc-grps = SYSTEM ; modified from default condition
;tau_t = 0.1 ; time constant, in ps
;ref_t = 300 ; reference temperature, one for each group, in K
nsttcouple = 1
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
nstpcouple = 1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; assign velocities from Maxwell distribution
;gen_temp = 300 ; temperature for Maxwell distribution
;gen_seed = -1 ; generate a random seed

Thank you in advance.

Hi,
maybe the difference you observed is due to how you calculate the Ndf. Constrains have also to considered when you evaluate Ndf (see eq 9 in documentation)
Ndf = 3N−Nc−Ncom
where Nc is the number of constraints imposed on the system. N = number of particles.
\Alessandra

Hi Alessandra,

Thank you for your prompt reply.

My simulation uses LINCS as the Constraint algorithm, how do I count the number of constraints given to the system?

Also, I think Ncom is the number of degrees of freedom, and in a non-vacuum system like my case, it should be set 3 considered to the center-of-mass motion in the x, y, z directions. Is this recognition correct?

I’m sorry for the basic question, but I would appreciate it if you could teach me.

Best regards,
Sakemi

Hi,

you find the information in the log file. But I understood that you have a system only with water molecules (as a test case). The standard water model are rigid. Thus number of contraints = n.water molecules x 3

Yes , see eq 9 here Molecular Dynamics — GROMACS 2021.3 documentation

regards
\Alessandra

Hi Alessandra,

Thank you for your teaching.

As you taught me, when I calculated with the number of constraints set to 3, I was able to calculate a reasonable temperature. I’m really thankful to you. I understand that these three constraints include two OH bond lengths and one HOH bond angle. However, when I checked the log file, it said that the number of constraints was 2. Does the number of constraints in the log file not include bond angles?

I’m sorry many times, but I would appreciate it if you could teach me.

Best regards,

Sakemi

Hi,

My guess is that md.log reports only the Nc based on the mdp file options. Do you still refer to an only-water simulation? Is the mdp file used the one above?

\Alessandra

Hi Alessandra,

No, now I trying to calculate the temperature of the cell which filled with naphthalene, but number of constraints is 3640 according to the generated log file, and manual calculated temperature is quit different from gmx energy calculated value.

The part of generated log file is shown below.

Input Parameters:
integrator = md-vv
tinit = 0
dt = 0.002
nsteps = 1500000
init-step = 0
simulation-part = 1
mts = false
comm-mode = Linear
nstcomm = 100
bd-fric = 0
ld-seed = -1142424129
emtol = 10
emstep = 0.01
niter = 20
fcstep = 0
nstcgsteep = 1000
nbfgscorr = 10
rtpi = 0.05
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstcalcenergy = 100
nstenergy = 500
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
periodic-molecules = false
verlet-buffer-tolerance = 0.005
rlist = 1
coulombtype = PME
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1
epsilon-r = 1
epsilon-rf = inf
vdw-type = Cut-off
vdw-modifier = Potential-shift
rvdw-switch = 0
rvdw = 1
DispCorr = EnerPres
table-extension = 1
fourierspacing = 0.16
fourier-nx = 32
fourier-ny = 32
fourier-nz = 32
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 0
epsilon-surface = 0
tcoupl = Nose-Hoover
nsttcouple = 1
nh-chain-length = 10
print-nose-hoover-chain-variables = false
pcoupl = No
pcoupltype = Isotropic
nstpcouple = 1
tau-p = 1
compressibility (3x3):
compressibility[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p (3x3):
ref-p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
refcoord-scaling = No
posres-com (3):
posres-com[0]= 0.00000e+00
posres-com[1]= 0.00000e+00
posres-com[2]= 0.00000e+00
posres-comB (3):
posres-comB[0]= 0.00000e+00
posres-comB[1]= 0.00000e+00
posres-comB[2]= 0.00000e+00
QMMM = false
qm-opts:
ngQM = 0
constraint-algorithm = Lincs
continuation = false
Shake-SOR = false
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 2
lincs-warnangle = 30
nwall = 0
wall-type = 9-3
wall-r-linpot = -1
wall-atomtype[0] = -1
wall-atomtype[1] = -1
wall-density[0] = 0
wall-density[1] = 0
wall-ewald-zfac = 3
pull = false
awh = false
rotation = false
interactiveMD = false
disre = No
disre-weighting = Conservative
disre-mixed = false
dr-fc = 1000
dr-tau = 0
nstdisreout = 100
orire-fc = 0
orire-tau = 0
nstorireout = 100
free-energy = no
cos-acceleration = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
simulated-tempering = false
swapcoords = no
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
applied-forces:
electric-field:
x:
E0 = 0
omega = 0
t0 = 0
sigma = 0
y:
E0 = 0
omega = 0
t0 = 0
sigma = 0
z:
E0 = 0
omega = 0
t0 = 0
sigma = 0
density-guided-simulation:
active = false
group = protein
similarity-measure = inner-product
atom-spreading-weight = unity
force-constant = 1e+09
gaussian-transform-spreading-width = 0.2
gaussian-transform-spreading-range-in-multiples-of-width = 4
reference-density-filename = reference.mrc
nst = 1
normalize-densities = true
adaptive-force-scaling = false
adaptive-force-scaling-time-constant = 4
shift-vector =
transformation-matrix =
grpopts:
nrdf: 20927
ref-t: 273
tau-t: 0.4
annealing: No
annealing-npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp-flags[ 0]: 0

Changing nstlist from 10 to 100, rlist from 1 to 1.038

Using 1 MPI thread
Using 6 OpenMP threads

Pinning threads with an auto-selected logical core stride of 1
System total charge: -0.000
Will do PME sum in reciprocal space for electrostatic interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Using a Gaussian width (1/beta) of 0.320163 nm for Ewald
Potential shift: LJ r^-12: -1.000e+00 r^-6: -1.000e+00, Ewald -1.000e-05
Initialized non-bonded Ewald tables, spacing: 9.33e-04 size: 1073

Generated table with 1019 data points for 1-4 COUL.
Tabscale = 500 points/nm
Generated table with 1019 data points for 1-4 LJ6.
Tabscale = 500 points/nm
Generated table with 1019 data points for 1-4 LJ12.
Tabscale = 500 points/nm
Long Range LJ corr.: 9.1706e-04

Using SIMD 4x8 nonbonded short-range kernels

Using a dual 4x8 pair-list setup updated with dynamic pruning:
outer list: updated every 100 steps, buffer 0.038 nm, rlist 1.038 nm
inner list: updated every 61 steps, buffer 0.001 nm, rlist 1.001 nm
At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be:
outer list: updated every 100 steps, buffer 0.207 nm, rlist 1.207 nm
inner list: updated every 61 steps, buffer 0.139 nm, rlist 1.139 nm

Using Lorentz-Berthelot Lennard-Jones combination rule
Removing pbc first time

Initializing LINear Constraint Solver

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
B. Hess and H. Bekker and H. J. C. Berendsen and J. G. E. M. Fraaije
LINCS: A Linear Constraint Solver for molecular simulations
J. Comp. Chem. 18 (1997) pp. 1463-1472
-------- -------- — Thank You — -------- --------

The number of constraints is 3640
There are: 8190 Atoms

Constraining the starting coordinates (step 0)
Center of mass motion removal mode is Linear
We have the following groups for center of mass motion removal:
0: rest
RMS relative constraint deviation after constraining: 6.26e-07
Initial temperature: 273.105 K

Started mdrun on rank 0 Thu Sep 9 21:17:28 2021

Best regards,

Sakemi