Re: [gmx-developers] details about Heat capacity caculation

Wed Nov 2 10:41:00 CET 2011

The best way to calculate heat capacity is to differentiate energy with respect to temperature for CV ( extract energy and corresponding temperatures from g_energy) and perform numerical differentiation to obtain Cv

Jestin B. Mandumpal Ph.D (WABRI), MRSC

Computational &amp; Physical Chemist

From: Haiqing Zhao &lt;haizhao at mtu.edu&gt;
Sent: Tue, 01 Nov 2011 21:42:38 
To: Discussion list for GROMACS development &lt;gmx-developers at gromacs.org&gt;
Subject: [gmx-developers] details about Heat capacity caculation

Dear Gmxers,

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

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

Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Average &nbsp; Err.Est. &nbsp; &nbsp; &nbsp; RMSD &nbsp;Tot-Drift

-------------------------------------------------------------------------------

Potential &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;-44.6187 &nbsp; &nbsp; &nbsp;0.007 &nbsp; 0.293318 &nbsp;0.0281861 &nbsp;(kJ/mol)

Kinetic En. &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 7.59226 &nbsp; &nbsp; 0.0032 &nbsp; &nbsp;0.19886 &nbsp;0.0201858 &nbsp;(kJ/mol)

Total Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -37.0264 &nbsp; &nbsp; 0.0091 &nbsp; 0.359776 &nbsp;0.0483719 &nbsp;(kJ/mol)

Temperature &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 300.076 &nbsp; &nbsp; &nbsp; 0.13 &nbsp; &nbsp;7.85971 &nbsp; 0.797822 &nbsp;(K)

Pressure &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 0.663728 &nbsp; &nbsp; &nbsp; 0.16 &nbsp; &nbsp;447.392 &nbsp;-0.388627 &nbsp;(bar)

Volume &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;14.3219 &nbsp; &nbsp; 0.0093 &nbsp; 0.133543 -0.0667073 &nbsp;(nm^3)

pV &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 0.438734 &nbsp; &nbsp;9.5e-05 0.00136276 -0.000679764 &nbsp;(kJ/mol)

Enthalpy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -17216.8 &nbsp; &nbsp; &nbsp; &nbsp;4.2 &nbsp; &nbsp;167.296 &nbsp; &nbsp;22.4922 &nbsp;(kJ/mol)

Temperature dependent fluctuation properties at T = 300.076. #constr/mol = 3

Isothermal Compressibility: 3.00555e-05 /bar

Adiabatic bulk modulus: &nbsp; &nbsp; &nbsp; &nbsp;33271.8 &nbsp;bar

Heat capacity at constant pressure Cp: &nbsp; &nbsp;67.9219 J/mol K

Thermal expansion coefficient alphaP: 9.66233e-05 1/K

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

Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Average &nbsp; Err.Est. &nbsp; &nbsp; &nbsp; RMSD &nbsp;Tot-Drift

-------------------------------------------------------------------------------

Total Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -33.2264 &nbsp; &nbsp; 0.0003 &nbsp;0.0900307 -0.000219237 &nbsp;(kJ/mol)

Temperature &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 300.008 &nbsp; &nbsp; 0.0063 &nbsp; &nbsp;1.95324 -0.0148441 &nbsp;(K)

Temperature dependent fluctuation properties at T = 300.008. #constr/mol = 3

Heat capacity at constant volume Cv: &nbsp; &nbsp;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

-- 

gmx-developers mailing list

gmx-developers at gromacs.org

http://lists.gromacs.org/mailman/listinfo/gmx-developers

Please don't post (un)subscribe requests to the list. Use the 

www interface or send it to gmx-developers-request at gromacs.org.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20111102/71bf5ac7/attachment.html>