# [gmx-developers] details about Heat capacity caculation

David van der Spoel spoel at xray.bmc.uu.se
Tue Nov 1 21:02:34 CET 2011

```On 2011-11-01 21:00, Haiqing Zhao wrote:
>
> Hi,
>
> 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)
> pV                              0.438734    9.5e-05 0.00136276 -0.000679764  (kJ/mol)
> Enthalpy                   -17216.8        4.2    167.296    22.4922  (kJ/mol)
>
> Yes, I did not doubt Cp=Var(Enthalpy)/k T^2.
> Actually my question is only about Enthalpy. As we see in the example, the total energy= -37.0264KJ/mol, pv = 0.4387 Kj/mol, How can enthalpy=  -17216.8 KJ/mol?
>
> sorry for the time.
You did not say how many molecules. Enthalpy is not normalized in the
g_energy output.

>
>
> On Nov 1, 2011, at 2:35 PM, David van der Spoel wrote:
>
>> On 2011-11-01 20:15, Haiqing Zhao wrote:
>>>
>>> Hi, David,
>>>
>>>
>>> I checked it. Version 4.5.5 really makes a difference from 4.5.3 which I used.
>>> In 4.5.5, Cv works exactly well by  nmol*(RMSD(E))^2/RT^2. (Even though I dont know here why need to multiply the # of molecular.)
>>>
>> This is supposed to work for liquids, hence you need to normalize.
>> How many molecules do you have, 465?
>> Heat capacity of 80 resembles water. I'm pretty sure it works correctly, since we just submitted a paper about it and I have looked into this in detail. We used:
>>
>> kB T^2 c_P = Var (Enthalpy)
>>
>> which comes straight from Allen and Tildesley.
>>
>>> But still I'm sure the Cp is not from simple Enthalpy Fluctuation,  as you see from the example below--167296 J/mol, and there is several orders mistake.
>>>
>>> Thanks.
>>>
>>> e.g.
>>>
>>> Enthalpy                   -17216.8        4.2    167.296    22.4922  (kJ/mol)
>>>
>>> You may want to use the -driftcorr flag in order to correct
>>> for spurious drift in the graphs. Note that this is not
>>> a substitute for proper equilibration and sampling!
>>>
>>> Temperature dependent fluctuation properties at T = 300.076.
>>>
>>> Heat capacities obtained from fluctuations do *not* include
>>> quantum corrections. If you want to get a more accurate estimate
>>> please use the g_dos program.
>>>
>>> a block-averaging error analysis (not implemented in g_energy yet)
>>> Volume                                   = 1.8548e-05 m^3/mol
>>> Enthalpy                                 =   -37.0255 kJ/mol
>>> Coefficient of Thermal Expansion Alpha_P =  0.0129106 (1/K)
>>> Isothermal Compressibility Kappa         = 3.00555e-10 (J/m^3)
>>> Adiabatic bulk modulus                   = 3.32718e+09 (m^3/J)
>>> Heat capacity at constant pressure Cp    =    80.3937 J/mol K
>>>
>>>
>>>
>>> On Nov 1, 2011, at 2:05 PM, David van der Spoel wrote:
>>>
>>>> On 2011-11-01 16:54, Haiqing Zhao wrote:
>>>>>
>>>>> 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).
>>>>>
>>>> Which gmx version? It's more or less correct in 4.5.4 but you really need qm correctiions. Check g_dos for that
>>>>
>>>>
>>>>
>>>>
>>>>> 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
>>>>>
>>>>
>>>>
>>>> --
>>>> David van der Spoel, Ph.D., Professor of Biology
>>>> Dept. of Cell&   Molec. Biol., Uppsala University.
>>>> Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
>>>> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
>>>> --
>>>> 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.
>>>
>>
>>
>> --
>> David van der Spoel, Ph.D., Professor of Biology
>> Dept. of Cell&  Molec. Biol., Uppsala University.
>> Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
>> spoel at xray.bmc.uu.se    http://folding.bmc.uu.se
>> --
>> 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.
>

--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se

```