Negative rotational kinetic energy by gmx traj

Er1ck · June 13, 2020, 1:44am

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
…

jalemkul · June 15, 2020, 1:38am

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

Er1ck · June 15, 2020, 2:02am

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.

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

Many thanks.

-Eric

Er1ck · June 15, 2020, 2:13am

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

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.

jalemkul · June 15, 2020, 5:30pm

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.

Er1ck · June 15, 2020, 6:39pm

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

jalemkul · June 15, 2020, 6:53pm

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.

Er1ck · June 15, 2020, 8:18pm

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;
```
}

gongrehk · January 18, 2022, 3:30am

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.

Er1ck · January 19, 2022, 12:44am

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.

gongrehk · January 19, 2022, 11:47am

Hi @Er1ck

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

ebriand · October 16, 2023, 11:54am

For future visitors to this topic, this bug has been fixed now (see here), and a correct version of gmx traj will be in 2023.3 and later versions.

Topic		Replies	Views
About the calculation result of translational energy and rotational energy User discussions analysis-tools	0	414	October 15, 2021
Calculation of kinetic energy from velocity and mass information User discussions	2	918	February 1, 2021
Compute the kinetic energy of a specific group in Gromacs User discussions	6	460	September 6, 2023
Rotational Kinetic Energy after Electric field User discussions analysis-tools	0	138	February 20, 2024
Why it cannot calculate the "total energy" when runs "gmx mdrun -rerun"? User discussions	2	864	June 23, 2020

Negative rotational kinetic energy by gmx traj

Related topics