[gmx-users] Problems with calculating Cv and Cp

[David van der Spoel]

Lum Nforbi wrote:
>     Hi Andrew Paluch,
> 
> 
>       You had asked me to supply the formulas that I attempted using to 
> calculate Cv. The formula from J. M. Haile is labelled isometric heat 
> capacity and given by Cv = Nk/(N-NT*(3(N/2)-1)<Ek*^-1>) and the one in 
> Allen and Tildesley is quite complicated to write out (has some special 
> characters in it).
>        I also found in one of the volumes of gmx-users where someone had 
> used the formula
> Cv = 1/(k_B*T*T)*Var(E) and instead of having the expected 75 J/Kmol, he 
> got 86 J/Kmol. He also mentioned that Var(E) can be obtained from RMSD 
> numbers. I know that RMSD numbers are displayed from g_energy, but how 
> do you get Var(E) from these numbers?
The RMSD is the square root of the variance.

The RMSD depends strongly on the temperature coupling. Berendsen 
coupling damps the fluctuations and hence gives too low Cp/Cv.
> 
> Thanks,
> Lum
>      
> 
> 
>     Date: Wed, 9 Dec 2009 20:45:57 -0500
>     From: Andrew Paluch <apaluch at nd.edu <mailto:apaluch at nd.edu>>
>     Subject: Re: [gmx-users] Problems with calculating Cv and Cp
>     To: Discussion list for GROMACS users <gmx-users at gromacs.org
>     <mailto:gmx-users at gromacs.org>>
>     Message-ID: <8BE66D7C-EA0F-4B04-8DDE-374529EC211F at nd.edu
>     <mailto:8BE66D7C-EA0F-4B04-8DDE-374529EC211F at nd.edu>>
>     Content-Type: text/plain; charset="us-ascii"
> 
>     Lum,
> 
>     You'll have better luck if you perform a few simulations at different
>     temperatures at the same pressure and mole number, and then
>     numerically differentiate the resulting enthalpy.  I believe that the
>     default heat capacity that Gromacs will print out is only valid for
>     NVE simulations, and I am not sure of the formulas that you are
>     referring to from Allen and Tildesley & J.M. Haile.  Also, don't
>     expect your calculated heat capacity with TIP3P to agree the
>     experimental value. You should under-predict the experimental value.
> 
>     Lastly, depending on what you are interested in, you may have better
>     luck with a water model other than TIP3P.
> 
>     Hope this helps,
> 
>     Andrew
>     _______________________________________________
>     _______________________________________________
>     Andrew Paluch
>     Department of Chemical and Biomolecular Engineering
>     University of Notre Dame du Lac
>     apaluch at nd.edu <mailto:apaluch at nd.edu>
>     _______________________________________________
>     _______________________________________________
> 
>     On Dec 9, 2009, at 8:18 PM, Lum Nforbi wrote:
> 
>      > Dear all,
>      >
>      >       I have run an 8 ns NPT simulation of 2000 molecules of  TIP3P
>      > water and I have a very low Cv value of 12.4748 J/mol K (factor =
>      > 0.000164481). The result is below. I ignored this value and have
>      > tried using formulas for Cv that I found in the two books: Allen
>      > and Tildesley & J. M. Haile but I can't come out with the right
>      > answer.
>      >      Has anyone ever calculated Cv or Cp for water manually from
>      > scratch and gotten the right answer? If so, please, could you give
>      > me the details of what you did?
>      >
>      > Statistics over 4000001 steps [ 0.0000 thru 8000.0005 ps ], 10 data
>      > sets
>      > All averages are exact over 4000001 steps
>      >
>      > Energy                      Average       RMSD     Fluct.
>      > Drift  Tot-Drift
>      > ------------------------------
>      > -------------------------------------------------
>      > Potential                  -79498.8    285.331    285.331
>      > 0.000110205   0.881641
>      > Kinetic En.                 14960.5    191.868    191.854
>      > -0.0010135   -8.10799
>      > Total Energy               -64538.3    353.238    353.232
>      > -0.000903288   -7.22631
>      > Temperature                 299.962    3.84702    3.84673
>      > -2.0321e-05  -0.162568
>      > Pressure (bar)             0.945534    194.613    194.613
>      > 0.000164281    1.31425
>      > Box-X                       3.94528 0.00514348 0.00514348
>      > 0          0
>      > Box-Y                       3.94528 0.00514348 0.00514348
>      > 0          0
>      > Box-Z                       3.94528 0.00514348 0.00514348
>      > 0          0
>      > Volume                      61.4097   0.240325    0.24032
>      > 6.42593e-07 0.00514075
>      > Density (SI)                  974.3    3.80629    3.80622
>      > -1.04747e-05 -0.0837974
>      > Heat Capacity Cv:      12.4748 J/mol K (factor = 0.000164481)
>      > Isothermal Compressibility: 2.27095e-05 /bar
>      > Adiabatic bulk modulus:        44034.4  bar
>      >
>      >
>      > Thank you,
>      >
>      > Lum
>      >
>      > --
> 


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