[gmx-users] Re: Wrong heat capacity (water)

[Andreas Kring]

Continuation of thread 
http://www.gromacs.org/index.php?option=com_wrapper&Itemid=84
(answered by David van der Spoel)

 >> Hello all,
 >>
 >> I ran a MD simulation of TIP4P water in the NVT ensemble (976 
molecules,
 >> density=0.98 g/cm3, T=300K (Nose-Hoover), SETTLE and PME enabled). This
 >> works out fine, i.e. I get a correct radial distribution function
 >> (g_OO), but the heat capacity (C_V) is wrong. I calculated it using the
 >> command:
 >>
 >> $  g_energy -f ener.edr -ndf 6 -nmol 976
 >>
 >> and get the value below (which sould have been approximately 75 
J/(mol*K)).
 >>
 >> Is the command above incorrect?'

 >Not necessarily. It would be good if you check the following things:
 >
 >- is the simulation in equilibrium?
 >- how large are the fluctuations for the equilibrium part?
 >- is the density correct?
 >
 >You can compute the cV from the fluctuations manually as well. Does that
 >match the value computed by the program?
 >

Thanks for you answer. Yes, you were right. The simulation was not in 
equilibrium. I've run a new simulation which I have made sure is in 
equilibrium. I have also checked that the density is correct, as you
suggest above, using the g_density command.

The numbers for the new simulation can be seen below (same simulation 
conditions as stated in the top of this mail), and the answer is still 
not right? With the command

$ g_energy -f ener.edr -ndf 6 -nmol 976

I get Cv = 305.413 J/(mol*K), while a manual calculation gives

Cv = 1/(k_B*T*T)*Var(E) = 85.6 J/(mol*K), where Var(E) was found using 
the RMSD numbers remembering that the unit of these are kJ/mole. The 
table value is Cv equal to approximately 75 J/(mole*K).

Any ideas?

Kind regards
Andreas



Energy                      Average       RMSD     Fluct.      Drift 
Tot-Drift
-------------------------------------------------------------------------------
Potential                  -41.3012   0.219318          0    8.76027 
87.6114
Kinetic En.                 7.48311   0.135439          0    4.46649 
44.6693
Total Energy               -33.8181   0.256267          0    13.2266 
  132.28
Temperature                 300.002    5.42982    5.40392   0.183467 
1.83485
Pressure (bar)             -269.723    384.283    384.133    3.71038 
37.1075
Heat Capacity Cv:      305.413 J/mol K (factor = 0.639443)