Negative rotational kinetic energy by gmx traj

GROMACS version:2018.4
GROMACS modification: No
Here post your question

Hi All,

I have a quick question about the rotational kinetic energy calculated by gmx traj -ekr.
According to the manual, the gmx traj uses center of mass and velocity data to estimate the rotational kinetic energy of a target group.
The output data by -ekr can contain negative values as the example shown below.
Since kinetic energy cannot be negative, is the output by -ekr is the changes in rotational kinetic energy? If this is the case, what would be the reference value, starting point (t=0)?

Thanks a bunch in advance,
-Eric

Example output:

@ title “Center of mass rotation”
@ xaxis label “Time (ps)”
@ yaxis label “Energy (kJ mol\S-1\N)”
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend “OCT”
0 -1.72973e-15
1 4.38805e-07
2 8.05563e-06
3 -0
4 6.69746e-06
5 8.35339e-06
6 4.2699e-05
7 0
8 9.4248e-06
9 1.23636e-05
10 9.10592e-06
11 4.66338e-05
12 -8.67984e-33
13 0
14 0
15 0
16 4.72635e-05
17 4.78465e-05
18 7.50094e-05
19 0.000143586
20 0.000206293

The negative values are all extremely small to the point that they should be interpreted as zero.

Hi Justin,

Thank you so much for such prompt response!

I only posted partial data in my original post.
The rotational kinetic energy calculated from my systems shows similar scale of valule on both positive and negative ranges, as shown below.
Imgur

I wonder if the rotational energy calculated by gmx traj -ekr is actually the changes in kinetic energy.

Many thanks.

-Eric

Here is another example.
Kinetic energy from the same systems was calculated but using smaller time steps.

Imgur

All these simulations applied constrained potentials to pull a single organic molecule into an organic surfaces along the reaction coordinate.

Each system uses NVT ensemble and contains an organic surface (frozen) and 1 octane (non-polar) / octanethiol (polar) molecule.

Check the ekrot() function in src/gromacs/gmxana/gmx_traj.cpp for the math that’s being used. I haven’t had time to go all the way through it, but that’s where you’ll find what the code is doing. It’s not computing a difference, e.g. over time.

Thank you so much for your suggestions.

Perhaps because of build on the HPC, I couldn’t find the file for “gmx_traj.cpp”.
I tried different version and even different HPC clusters.
I probably missed something.

The closest one I can find is the “gmx-traj.1
which contains the following entry of ‘ekr


.sp
Options \fB-ekt\fP and \fB-ekr\fP plot the translational and
rotational kinetic energy of each group,
provided velocities are present in the trajectory file.
This implies \fB-com\fP&.
.sp

For your reference, below is my command input and output for example.

which gmx_sp
/usr/common/software/gromacs/2018.4.knl/bin/gmx_sp

find . -name '*traj*'
./include/gromacs/trajectory
./include/gromacs/trajectory/trajectoryframe.h
./include/gromacs/trajectoryanalysis.h
./include/gromacs/trajectoryanalysis
./share/man/man1/gmx-traj.1
./share/man/man1/gmx-trajectory.1
./share/man/man1/gmx-nmtraj.1

find . -name '*cpp*'
./share/gromacs/template/template.cpp

My suggestion was to inspect the source code. You won’t find the details of the calculations in an installed set of binaries and man pages.

Thank you so much for pointing this out.
I just find it.
Based on the script below, it seems that -ekr is calculated by
Ekr = 0.5 * acm ^2 * TCM

Not exactly sure the how the rotational velocity acm and ‘moment of inertia ?’ TCM are derived from the mass, position, and velocity data.
If that was the case, rotational kinetic energy shall never be negative.

Thanks again.

  • static real ekrot(rvec x[], rvec v[], real mass[], int isize, int index[])

  • {

  • real          TCM[5][5], L[5][5];
    
  • double        tm, m0, lxx, lxy, lxz, lyy, lyz, lzz, ekrot;
    
  • rvec          a0, ocm;
    
  • dvec          dx, b0;
    
  • dvec          xcm, vcm, acm;
    
  • int           i, j, m, n;
    
  • clear_dvec(xcm);
    
  • clear_dvec(vcm);
    
  • clear_dvec(acm);
    
  • tm = 0.0;
    
  • for (i = 0; i < isize; i++)
    
  • {
    
  •     j   = index[i];
    
  •     m0  = mass[j];
    
  •     tm += m0;
    
  •     cprod(x[j], v[j], a0);
    
  •     for (m = 0; (m < DIM); m++)
    
  •     {
    
  •         xcm[m] += m0*x[j][m]; /* c.o.m. position */
    
  •         vcm[m] += m0*v[j][m]; /* c.o.m. velocity */
    
  •         acm[m] += m0*a0[m];   /* rotational velocity around c.o.m. */
    
  •     }
    
  • }
    
  • dcprod(xcm, vcm, b0);
    
  • for (m = 0; (m < DIM); m++)
    
  • {
    
  •     xcm[m] /= tm;
    
  •     vcm[m] /= tm;
    
  •     acm[m] -= b0[m]/tm;
    
  • }
    
  • lxx = lxy = lxz = lyy = lyz = lzz = 0.0;
    
  • for (i = 0; i < isize; i++)
    
  • {
    
  •     j  = index[i];
    
  •     m0 = mass[j];
    
  •     for (m = 0; m < DIM; m++)
    
  •     {
    
  •         dx[m] = x[j][m] - xcm[m];
    
  •     }
    
  •     lxx += dx[XX]*dx[XX]*m0;
    
  •     lxy += dx[XX]*dx[YY]*m0;
    
  •     lxz += dx[XX]*dx[ZZ]*m0;
    
  •     lyy += dx[YY]*dx[YY]*m0;
    
  •     lyz += dx[YY]*dx[ZZ]*m0;
    
  •     lzz += dx[ZZ]*dx[ZZ]*m0;
    
  • }
    
  • L[XX][XX] =  lyy + lzz;
    
  • L[YY][XX] = -lxy;
    
  • L[ZZ][XX] = -lxz;
    
  • L[XX][YY] = -lxy;
    
  • L[YY][YY] =  lxx + lzz;
    
  • L[ZZ][YY] = -lyz;
    
  • L[XX][ZZ] = -lxz;
    
  • L[YY][ZZ] = -lyz;
    
  • L[ZZ][ZZ] =  lxx + lyy;
    
  • m_inv_gen(&L[0][0], DIM, &TCM[0][0]);
    
  • /* Compute omega (hoeksnelheid) */
    
  • clear_rvec(ocm);
    
  • ekrot = 0;
    
  • for (m = 0; m < DIM; m++)
    
  • {
    
  •     for (n = 0; n < DIM; n++)
    
  •     {
    
  •         ocm[m] += TCM[m][n]*acm[n];
    
  •     }
    
  •     ekrot += 0.5*ocm[m]*acm[m];
    
  • }
    
  • return ekrot;
    
  • }

Hi Er1ck,
Did you solve your problem?
I meet the same problem.
I tried to extract the rotational kinetic energy of a single water molecule in protein solutions.
Surprisingly, the rotational kinetic energy at some moment is negative, as shown below.

Command line:
gmx traj -s md_water1.tpr -f md0_100ps.trr -n water1.ndx -ov water1_V.xvg -ekr water1_EKR.xvg -ekt water1_EKT.xvg
gmx traj is part of G R O M A C S:

GROwing Monsters And Cloning Shrimps

@ title “Center of mass rotation”
@ xaxis label "Time (ps)"
@ yaxis label "Energy (kJ mol\S-1\N)"
@TYPE xy
@ view 0.15, 0.15, 0.75, 0.85
@ legend on
@ legend box on
@ legend loctype view
@ legend 0.78, 0.8
@ legend length 2
@ s0 legend “r_130”
0 2.52655
0.8 5.94535
1.6 0.908079
2.4 5.10452
3.2 23.8468
4 -6.76135
4.8 0.725384
5.6 3.57106
6.4 -2.08162
7.2 0.285385
8 5.13313
8.8 6.86377
9.6 -0.0344337
10.4 1.16367
11.2 0.318203
12 -0.0536116
12.8 2.48021
13.6 -6.42674
14.4 -4.91415
15.2 -3.47536
16 1.94549
16.8 -257.257
17.6 53.4847
18.4 1.48988
19.2 4.78093
20 -0.77453

I also checked the translational kinetic energy and it is just fine.
Do you have any idea about that?
Many thanks.

Hi @gongrehk

Thanks for asking.
No, I haven’t. I pivoted my project and had stopped pursuing this direction further.
I did try to look into the source code, but to no avail.

To obtain the rotational kinetic energy, perhaps you can try different version of GROMACS or use a different simulation package. Sorry I wish I could be of more help.

Hi @Er1ck

Thanks very much for your advice.
Perhaps I should try with a latest version of GROMACS first.