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<br><br>Jestin B. Mandumpal Ph.D (WABRI), MRSC<br>
Computational &amp; Physical Chemist<br><br><br>From: Haiqing Zhao &lt;haizhao@mtu.edu&gt;<br>Sent: Tue, 01 Nov 2011 21:42:38 <br>To: Discussion list for GROMACS development &lt;gmx-developers@gromacs.org&gt;<br>Subject: [gmx-developers] details about Heat capacity caculation<br><br>
Dear Gmxers,<br>
<br>
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).<br>
<br>
One data for NPT is: ( -nmol &nbsp; 465 &nbsp;-nconstr &nbsp; &nbsp;3 )<br>
<br>
Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Average &nbsp; Err.Est. &nbsp; &nbsp; &nbsp; RMSD &nbsp;Tot-Drift<br>
-------------------------------------------------------------------------------<br>
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)<br>
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)<br>
Total Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -37.0264 &nbsp; &nbsp; 0.0091 &nbsp; 0.359776 &nbsp;0.0483719 &nbsp;(kJ/mol)<br>
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)<br>
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)<br>
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)<br>
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)<br>
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)<br>
<br>
Temperature dependent fluctuation properties at T = 300.076. #constr/mol = 3<br>
Isothermal Compressibility: 3.00555e-05 /bar<br>
Adiabatic bulk modulus: &nbsp; &nbsp; &nbsp; &nbsp;33271.8 &nbsp;bar<br>
Heat capacity at constant pressure Cp: &nbsp; &nbsp;67.9219 J/mol K<br>
Thermal expansion coefficient alphaP: 9.66233e-05 1/K<br>
<br>
And data for NVT, ( -nmol &nbsp; 7631 &nbsp;-nconstr &nbsp; &nbsp;3 )<br>
<br>
Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;Average &nbsp; Err.Est. &nbsp; &nbsp; &nbsp; RMSD &nbsp;Tot-Drift<br>
-------------------------------------------------------------------------------<br>
Total Energy &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; -33.2264 &nbsp; &nbsp; 0.0003 &nbsp;0.0900307 -0.000219237 &nbsp;(kJ/mol)<br>
Temperature &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 300.008 &nbsp; &nbsp; 0.0063 &nbsp; &nbsp;1.95324 -0.0148441 &nbsp;(K)<br>
<br>
Temperature dependent fluctuation properties at T = 300.008. #constr/mol = 3<br>
Heat capacity at constant volume Cv: &nbsp; &nbsp;70.1813 J/mol K<br>
<br>
(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.)<br>
<br>
I appreciate if anybody could help me to check it.<br>
Thanks!!<br>
<br>
---------------------------------------------<br>
Haiqing Zhao<br>
PH.D. candidate<br>
in Computational Biophysics<br>
Michigan Technological University<br>
http://www.phy.mtu.edu/~haizhao<br>
<br>
-- <br>
gmx-developers mailing list<br>
gmx-developers@gromacs.org<br>
http://lists.gromacs.org/mailman/listinfo/gmx-developers<br>
Please don't post (un)subscribe requests to the list. Use the <br>
www interface or send it to gmx-developers-request@gromacs.org.<br>
<br><A HREF="http://sigads.rediff.com/RealMedia/ads/click_nx.ads/www.rediffmail.com/signatureline.htm@Middle?" target="_blank"><IMG SRC="http://sigads.rediff.com/RealMedia/ads/adstream_nx.ads/www.rediffmail.com/signatureline.htm@Middle"></A><br><table width="578" border="0" cellspacing="0" cellpadding="0"><tr><td><span style="font-family:Arial, Helvetica, sans-serif; font-size:12px; color:#393939;">Follow&nbsp;<span style="color:#0000CC;"><b><u><a href="http://track.rediff.com/click?url=___http://dealhojaye.rediff.com?sc_cid=rediffmailsignature___&cmp=signature&lnk=rediffmailsignature&newservice=deals" target="_blank">Rediff Deal ho jaye!</a></u></b></span>&nbsp;to get exciting offers in your city everyday.</span></td></tr></table>