bjorn.sat...@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.
sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se
_______________________________________________
gmx-users mailing list gmx-users@gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
