[gmx-developers] details about Heat capacity caculation

[David van der Spoel]

On 2011-11-01 20:15, Haiqing Zhao wrote:
>
> Hi, David,
>
> thanks for reply.
>
> 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.
>
> Please help me check it.
> 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.
>
> WARNING: Please verify that your simulations are converged and perform
> 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