Calculation of kinetic energy from velocity and mass information

User discussions
1

GROMACS version: 2018.1
GROMACS modification: No
Hello dear GROMACS users and developers,
I fall into trouble during some consistency check of my molecular dynamics trajectories. I found out I am not able to reproduce the proper value of kinetic energy of the system. The value I get from gmx energy tool (~6952 kJ/mol) and the value I calculate from mass (tpr file) and velocity (trr file) data differs significantly. Thus I would rather start with a snippet of a python code (utilizing MDAnalysis package) I am using for the data calculation rather then describing the details of the system as I must be missing something elementary.

import MDAnalysis as MD
import numpy as np
mass_conversion = 1.660538921*1e-27  # convert atomic mass to kg
speed_conversion = 1e3  # convert nm/ps to m/s
N_a = 6.02214076*1e23  # Avogadro's number
topo_file = MD.topology.TPRParser.TPRParser('diala_water.tpr')
trajectory = MD.coordinates.TRR.TRRReader('diala_water.trr')
topo = topo_file.parse()
vel_arr = trajectory.ts.velocities
mass_arr = topo.masses.values
nof_particles = topo.n_atoms
nof_residues = topo.n_residues
sum = 0
for i in range(nof_particles):
    sum += mass_arr[i] * np.linalg.norm(vel_arr[i, :])**2
kin = sum / 2 * mass_conversion * speed_conversion**2 * N_a / nof_residues
print('Kinetic energy %.2f kJ/mol' % (kin/1000))

If i didn’t make a mistake in a code it is suppose to make a sum of “m_i*|v_i|^2” term over all atoms in the simulation box. Where m_i is mass and |v_i| the norm of velocity of particle i. Then take a half of this sum and multiply it by conversion factor due to mass and velocity units. Last step taken was to evalute the energy per mol. Though this was a source of a little confusion for me as I am not sure whether the value given by gmx energy tool is given as of mol of atom or molecules I tried to recalculate the value for both possibilities. And I get the values 247.58, 747.56 and 654117.50 kJ/mol in case of mol of atoms, mol of residues and mol of simulation boxes. Can you please give me a hint how to get the proper value? Thanks in advance! Jiri

2

The mol is unit is Avogadro’s number, this has nothing to do with the number of molecules in the system. If the input in is GROMACS MD units and you don’t do any unitconversion, you should get the right answer right away. If you want to convert you should remove any number of atoms/residues/mols and your net conversion factor should be 1.

3

Thank you for your reply, hess! That solves it. There were three discrepancies.

  1. TRRReader from MDAnalysis package really did some unit conversion for velocities by default. So it must be forbidden explicitly.
    trajectory = MD.coordinates.TRR.TRRReader(‘diala_water.trr’, convert_units=False)
  2. step zero in trajectory file corresponds to second line in energy xvg file not the first one.
  3. You are right with the conversion factor of 1 for mol recalculation. Thus the correct form is this:
    kin = sum / 2 * mass_conversion * speed_conversion**2 * N_a / 1
    Thank you again I appreciate your help! However you did not convinced me that the value is the not the kinetic energy per mol of simulation boxes :)