[gmx-users] Bug in g_energy for calculating heat capacity?

[David van der Spoel]

Bjorn.Sathre at ift.uib.no wrote:
> On Sunday 21. December 2008 10.28.19 David van der Spoel wrote:
>> Zhang Zhigang wrote:
>> > Hi, all,
>> >      Actually I think others have also noticed this "bug": g_energy
>> > often generate incorrect heat capacity.  When I checked the source
> code
>> > of gmx_energy.c, I found the calculation of heat capacity is only
>> > related with the fluctuation of temperature, which is actually the
>> > algorithm for NVE ensemble only (refer to Allen & Tildesley, Computer
>> > simulation of liquids, pp. 53, eqn. 2.82). For other ensembles, e.g.
> NVT
>> > ensemble, the heat capacity should be related with the fluctuation of
>> > energies.
>> >      Another problem: even with the energy fluctuation algorithm for
> the
>> > heat capacity, I found it is not easy to get an convergent value. The
>> > most important factor is that the fluctuations are sensitively related
>> > with the temperature (pressure) coupling time constants for NVT (NPT)
>> > ensembles.
>>
>> With a proper ensemble, e.g. NoseHoover, the time structure of the
>> fluctuations is influenced by the T-coupling, but not the amplitude,
>> which is what matters in the equations for heat capacity, isn't it?
> 
> Is this conclusion really valid?
> According to one of the Inventors of the extended Hamiltonian formalism:
> Nose (MolPhys 52(2) p. 255-268 (1984))  on p. 263 we find the following:
> 
> "Figure 1 shows the evolution of the instantaneous temperature through the
> first 3 runs. In the constant temperature method, the temperature reaches
> T_eq quickly and fluctuates around this value. The approximate formula for
> the fluctuation of temperature in the microcanonical ensemble[13,14],
> 
> ( He then gives the formula (3.1) which is identical to the Cv formula for
> the
> microcanonical ensemble used by GROMACS4.02 and then continues...)
> 
> .. shows that the temperature fluctuations in the canonical ensemble are
> larger than those in the microcanonical ensemble. We can readily recognize
> this in figure 1."
> 
> Looking at his figure 1 I am strongly inclined to agree. It is obvious he
> refers to the fluctuations being larger in magnitude.
> 
> Also based on my own rework of an old manuscript of ours dealing with Cv
> of a very small system of long-chained alkanes, I find that there are 
> obvious
> discrepancies, even when including
> the quantum corrections of eq (1.4) in the GROMACS 4.0 Manual. (This  is
> an odd place to put the corrected formula for Cv(QM) by the way, it 
> should be
> along with the definition of temperature and formula for counting
> classical DOFs in chapter 3)
> 
> I suggest taking Mr Nose on his word and instead using his formula (2.20)
> in the above cited article, suitable for the canonical (NVT)-ensemble:
> 
> cv(classical)=(k_B f)/2N+ Var(Potential)/(N k_B T_av^2)    (atomic units)

This can easily be verified using xmgrace. Have you tried?

> 
> where  k_B is the Bolzmann constant, f is the total number of degrees of
> freedom in the system(excluding the unphysical additional degree of
> freedom in the extended formalism) and N is the total particle number in 
> the
> system.
> It is understood that in MD you typically conserve total linear momentum 
> and
> thus have left f'=f-3.
> 
> He also has an alternative formulation for f large(his 2.24),
> 
> cv(classical)= (k_B f^2)/N* Var(s)/<s>
> 
> where s is the new dynamical variable introduced in the extended
> Hamiltonian formalism.
>  
> I do not believe, however, that the difference will be significant, except
> for quite small systems (which I am afraid I have touched upon) but it 
> would
> be nice to have GROMACS use the formally correct solution which, unless
> Nose's old paper is fundamentally flawed, I do not believe it does at 
> present.
> 
> Regards
> Bjørn Steen Sæthre
> PhD student
> Institute of Physics and Technology
> University of Bergen
> Norway
> 


-- 
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