Reading velocities and forces from trr file

[gmonet]

Hi all, I'm using pytraj to read trr trajectory files generated with GROMACS :

import pytraj
print(pytraj.__version__)
#2.0.5
import MDAnalysis as mda
print(mda.__version__)
#1.0.0
path_top= 'init.pdb'
path_traj ='traj.trr'
u = mda.Universe(path_top, path_traj)
traj = pytraj.TrajectoryIterator(path_traj , path_top)

diff = u.trajectory[0].positions/traj[0].coordinates
print(np.nanmin(diff), np.nanmax(diff)) 
# 0.99999994041066 1.0000000595886653 # Same coordinates

%timeit u.trajectory[0].positions
# 2.94 ms ± 3.95 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit traj[0].coordinates
# 194 µs ± 1.64 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

I don't know why pytraj is much faster than MDAnalysis but that's why I use pytraj so far. My problem appears when I want to read velocities and forces:

diff = u.trajectory[0].velocities/traj[0].velocity
print(np.nanmin(diff), np.nanmax(diff))
# 20.454998780973654 20.455001219154948

diff = u.trajectory[0].forces/traj[0].force
print(np.nanmin(diff), np.nanmax(diff))
# 4.183999900247322 4.184000299235091

As you can see, the speeds and forces read by pytraj and MDAnalysis are different. After verification, the velocities given by MDAnalysis are correct. I may have missed a unit conversion but I did not find any information about this on the pytraj site.

[hainm]

hi @drroe, can you please comment in this case. Thanks.

[hainm]

Here are the units in MDAnalysis for TRR: https://github.com/MDAnalysis/mdanalysis/blob/e1c90ae4e3151423bc25ac9a6e111232c6022491/package/MDAnalysis/coordinates/TRR.py#L61-L62

[hainm]

diff = u.trajectory[0].forces/traj[0].force
print(np.nanmin(diff), np.nanmax(diff))
# 4.183999900247322 4.184000299235091

From my google-ful: 1 kcal = 4.184 kJ so I guess cpptraj is following AMBER suite to use kcal unit.

[hainm]

diff = u.trajectory[0].velocities/traj[0].velocity
print(np.nanmin(diff), np.nanmax(diff))
# 20.454998780973654 20.455001219154948

Per David A. Case in this post: https://structbio.vanderbilt.edu/archives/amber-archive/2008/4377.php

[The unit of length in Amber is Angstroms; the unit of time is 1/20.455 ps.
Velocities have units of length/time.]

So I guess pytraj (cpptraj under the hood) reads the v/f correctly. Cheers.

[gmonet]

Thanks a lot hainm ! The documentation of cpptraj is a bit misleading about the velocity's unit. It might be worth mentioning the units used in the pytraj documentation as well as in the MDAnalysis documentation. Best.