[gmx-users] Bug in g_energy for calculating heat capacity?
David van der Spoel
spoel at xray.bmc.uu.se
Wed Jan 21 09:00:47 CET 2009
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
>> > 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.
>> > ensemble, the heat capacity should be related with the fluctuation of
>> > energies.
>> > Another problem: even with the energy fluctuation algorithm for
>> > 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
> 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
> 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
> It is understood that in MD you typically conserve total linear momentum
> 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
> 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
> Bjørn Steen Sæthre
> PhD student
> Institute of Physics and Technology
> University of Bergen
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
More information about the gromacs.org_gmx-users