Entropy calculation with gmx anaeig: inconsistent results with QH analysis
Summary
It seems that there is an issue with the entropy calculation by gmx anaeig, where the results calculated using Schlitter's formula and QH analysis differ by even one order of magnitude and are thus inconsistent. The issue is with the QH calculation.
GROMACS version
Verified to exist in 2022.5 and 2024.1
Steps to reproduce
This issue appears to happen when using gmx covar + gmx anaeig to perform entropy calculations on trajectories. Specifically, gmx anaeig -entropy outputs two entropy values, one calculated using Quasi-Harmonic analysis and one using Schlitter's formula, which appear to be inconsistent, and differ by up to an order of magnitude. I was able to reproduce this with different molecular systems, e.g. a simple protein simulation in water in the NPT ensemble, as well as simulations of nanoparticles such and dendrimers (different systems, different force fields, and no obvious issues or instabilities in the simulations themselves). To reproduce the bug, it seems to be sufficient to calculate eigenvalues and eigenvectors of a given system using gmx covar and then performing entropy analysis using gmx anaeig with the -entropy flag.
What is the current bug behavior?
The entropy result by gmx anaeig using the two calculation methods is inconsistent. In one example protein simulation in water (CHARMM36m FF, NPT ensemble), the result of gmx anaeig -entropy is:
Schlitter: 3228.01 J/molK
QH: 34517.5 J/molK
Again, I was able to confirm this behaviour on very different simulations, so it should not depend on the specific simulation being analyzed.
What did you expect the correct behavior to be?
Entropy values should be consistent, or at the very least differ by a small value (as reported e.g. in https://doi.org/10.1063/1.1401821 )
Possible fixes
I suspect that this issue might be connected with the one that was briefly reported in issue #3481 (closed). By looking at the gromacs source code, the calculation of the entropy starting from the eigenvalues is performed in /src/gromacs/gmxana/thermochemistry.cpp. It seems to me that the calculation of frequencies starting from the eigenvalues performed in the function eigval_to_frequency is not correct. In recent versions (i think 2019 onwards), the conversion is performed by taking the square root of the product of the eigenvalues times a factor (called "factor_gmx_to_omega2") equal to 1e21/(NA*AMU), where NA is Avogadro's number and AMU is the conversion between a.m.u. and kg. In the current versions, the calculation of the frequency is thus performed as:
std::sqrt(std::max(0.0, eigval) * factor_gmx_to_omega2)
I think the issue might be that the correct way of calculating omega is by first calculating the factor as:
factor_gmx_to_omega2 = kT1e18 / AMU
where k is the Boltzmann constant, T the temperature, 1e18 is the reciprocal of the conversion factor between nm^2 (unit of the non-mass-weighted eigenvalues as produced by gmx covar) and m^2 (SI). Omega should be then given by using this corrected factor in:
w = sqrt(1/lambda * factor_gmx_to_omega2)
where lambda is the eigenvalue in nm^2. This should be the correct implementation of the equation relating eigenvalues to frequencies as reported e.g. in https://doi.org/10.1063/1.1401821, and indeed gives realistic values.
I think the issue might not be present in versions before 2018, since there the calculation is performed in /src/gromacs/gmxana/gmx_anaeig.cpp and the corresponding function, calc_entropy_qh, indeed calculates the frequencies as follows:
lambda = eigval[i]*AMU;
w = std::sqrt(BOLTZMANN*temp/lambda)/NANO;
which should be the same as the suggested solution above.