[gmx-developers] details about Heat capacity caculation

[Haiqing Zhao]

Dear Gmxers,

I have one question about how Gromacs calculate the heat capacity. I dont think it is only Enthalpy fluctuation÷(kT^2) in NPT ensemble. And also I checked it in NVT, it is not simply Total energy(/internal energy) fluctuation÷(kT^2).

One data for NPT is: ( -nmol   465  -nconstr    3 )

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Potential                  -44.6187      0.007   0.293318  0.0281861  (kJ/mol)
Kinetic En.                 7.59226     0.0032    0.19886  0.0201858  (kJ/mol)
Total Energy               -37.0264     0.0091   0.359776  0.0483719  (kJ/mol)
Temperature                 300.076       0.13    7.85971   0.797822  (K)
Pressure                   0.663728       0.16    447.392  -0.388627  (bar)
Volume                      14.3219     0.0093   0.133543 -0.0667073  (nm^3)
pV                         0.438734    9.5e-05 0.00136276 -0.000679764  (kJ/mol)
Enthalpy                   -17216.8        4.2    167.296    22.4922  (kJ/mol)

Temperature dependent fluctuation properties at T = 300.076. #constr/mol = 3
Isothermal Compressibility: 3.00555e-05 /bar
Adiabatic bulk modulus:        33271.8  bar
Heat capacity at constant pressure Cp:    67.9219 J/mol K
Thermal expansion coefficient alphaP: 9.66233e-05 1/K

And data for NVT, ( -nmol   7631  -nconstr    3 )

Energy                      Average   Err.Est.       RMSD  Tot-Drift
-------------------------------------------------------------------------------
Total Energy               -33.2264     0.0003  0.0900307 -0.000219237  (kJ/mol)
Temperature                 300.008     0.0063    1.95324 -0.0148441  (K)

Temperature dependent fluctuation properties at T = 300.008. #constr/mol = 3
Heat capacity at constant volume Cv:    70.1813 J/mol K

(BTW, I found that (1)if you choose Enthalpy as output, g_energy will always give Cp, no matter what ensemble it is. (2)And only you choose tot energy and temp as output, no other terms, g_energy will give you Cv. If not, you cannot get Cv by g_energy.)

I appreciate if anybody could help me to check it.
Thanks!!

---------------------------------------------
Haiqing Zhao
PH.D. candidate
in Computational Biophysics
Michigan Technological University
http://www.phy.mtu.edu/~haizhao