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

Andreas Kring akring at fys.ku.dk
Thu Apr 17 12:20:11 CEST 2008

Continuation of thread 
(answered by David van der Spoel)

 >> Hello all,
 >> I ran a MD simulation of TIP4P water in the NVT ensemble (976 
 >> 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 
 >> 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

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

More information about the gromacs.org_gmx-users mailing list