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