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.