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.