[gmx-developers] [Fwd: coulomb SR in 3.3cvs]
Michel Cuendet
michel.cuendet at epfl.ch
Tue Aug 9 12:20:18 CEST 2005
Hi Berk,
Yeah, we sorted this out already with David. I don't know why the
discussion went off the list, and why you were not included in the cc's.
Sorry for that. The quotation below summarizes our findings.
Basically you are right, all this comes from the computation of the
kinetic energy, which is correct in 3.3b, but not in 3.2.1. The irony is
that you and me had worked out this kinetic energy problem, along with
the integrator problem, together in mai (see quotation below). And I was
not smart enough to make the connection between the two problems directly...
The consequences of all this are quite dramatic, because they affect all
previous simulations made with gmx up to 3.3b, independently of which
thermostat and barostat is used. The difference in temperature is not so
alarming, but the difference in pressure and energy is. When rerunning
the frames of a 3.2.1 simulation with 3.3b, I find:
temperature average : 303.47 K
pressure average : 1687 bar
potential energy average difference : ~3500 KJ/mol (between two
different runs from same restart)
The good news is anyway that I can happily continue to use 3.3b on my
itaniums, because it seem to be "more right" than 3.2.1, and not the
contrary.
Bye,
Michel
>-------- Original Message --------
>
>
> Subject: Re: [Fwd: Re: results]
> Date: Mon, 08 Aug 2005 21:14:44 +0200
> From: David <spoel at xray.bmc.uu.se>
> To: gmx-developers at gromacs.org
> CC: michel.cuendet at epfl.ch
>
>Hi guys,
>
>this mail that I just got from Michel Cuendet summarizes why simulations
>with 3.3b give different results from 3.2.1. I guess Michel is right in
>that "something needs to be done" about implementing further algorithms.
>As long as we use Berendsen coupling then Leap-Frog is fine of course.
>
>
>-------- Forwarded Message --------
>From: Michel Cuendet <michel.cuendet at epfl.ch>
>Reply-To: michel.cuendet at epfl.ch
>To: David <spoel at xray.bmc.uu.se>
>Subject: Re: results
>Date: Mon, 08 Aug 2005 16:26:08 +0200
>Hi David,
>
>OK, I've got it.
>
>First check the virial, it looks pretty much ok. But :
> 321 33a
>@ s2 legend "pV" 2388.021240 3662.385742
>@ s1 legend "Pressure (bar)" 44.703457 68.559402
>
>@ s1 legend "Pres-XX (bar)" 121.256340 145.069550
>@ s2 legend "Pres-XY (bar)" 5.238311 4.988041
>@ s3 legend "Pres-XZ (bar)" -53.949970 -53.746536
>@ s4 legend "Pres-YX (bar)" 5.238311 4.987853
>@ s5 legend "Pres-YY (bar)" -91.220413 -67.740036
>@ s6 legend "Pres-YZ (bar)" 42.322739 42.367237
>@ s7 legend "Pres-ZX (bar)" -53.949970 -53.747005
>@ s8 legend "Pres-ZY (bar)" 42.322739 42.366676
>@ s9 legend "Pres-ZZ (bar)" 104.074448 128.348679
>
>@ s0 legend "Temperature" 300.639160 303.291931
>@ s0 legend "Kinetic En." 216638.687500 218550.265625
>
>-> KINETIC ENERGY IS INCONSISTENT
>-> This changes the pressure, the volume, and finally the potential energy !
>
>trjconv -f traj.trr -e 0.0 -o step0.gro -ndec 6 shows that the
>velocities are exactly the same at step 0.
>
>-> THE DISCREPANCY COMES FROM THE FOLLOWING CHANGES :
>-----------------------------------------------------------------------------
>Berk Hess 7 Apr 2005 Temperature coupling The Berendsen
>and Nose-Hoover temperature coupling now use the temperature at the half
>step. Nose-Hoover now uses a reversible integrator: Holian et al. Phys
>Rev E 52(3) : 2338, 1995
>-----------------------------------------------------------------------------
>Berk Hess 7 Apr 2005 Kinetic energy The kinetic energy in
>Gromacs is now determined as the average of the kinetic energy of the
>velocities at the half steps around the full step. This is four times
>more accurate than the old method which used the full step particle
>velocities which were obtained as the averages of the half steps.
>-----------------------------------------------------------------------------
>
>The great IRONY is that this change was originated by one of my own
>comments on gmx 3.2.1 (look for "wrong Nose-Hoover integrator" in the
>mailing list, see also the mail from Berk quoted below). ALL SIMULATIONS
>PREVIOUSLY DONE IN GROMACS (up to version 3.3.a) HAVE A "TRUE"
>TEMPERATURE HIGHER THAN THE NOMINAL TEMPERATURE. As I had done all these
>integrator tests for the case NVT, I just hadn't figured out that a
>slight temperature change could have such dramatic consequences on the
>pressure and the potential energy itself !
>
>According to the changes above, the temperature estimation is now
>correct in gmx. BUT the integrator is reversible and correct only for
>the NVT case. I haven't checked in the code what Berk has implemented
>for the NPT case, but the integrator of Holian et al. certainly doesn't
>work as-is fir the NPT case. Now that we have seen the importance of
>having a correct time reversible ("at the right 1/2 time step")
>estimation of the temperature, we should definitely check how the virial
>and the pressure are estimated !!
>
>The integrator of Holian et al. is just a fix to the regular leap-frog,
>which I found easy to hack in in gromacs to make my little tests. It is
>reversible, but I'm not sure that it has the nice phase space volume
>conservation properties. I would once again recommend, as gromacs is
>approaching a grown up age, to implement state of the art fully
>conststent integrators as described in :
>
>Martyna et al. "Explicit reversible integrators for extended systems
>dynamics", Molecular Phys. 87(5) : 1117, 1996.
>Ciccotti and Kalibaeva, "Molecular Dynamics of COmplex Systems: ...",
>Lecture Notes in Physics 640 :150, 2004.
>Kalibaeva et al. "Constant pressure - Constant Temperature Molecular
>Dynamics: a correct constrained ensemble using the molecular virial",
>Mol. Phys. 101(6) : 765, 2003
>
>Ahhhh what a relief when something like this finds a nice explaination ....
>
>Cheers,
>
>Michel
>
>
>
>
>>> -------- Original Message --------
>>> Subject: temperature fix
>>> Date: Thu, 31 Mar 2005 19:05:13 +0200
>>> From: Berk Hess <hessb at mpip-mainz.mpg.de>
>>> To: Erik Lindahl <lindahl at sbc.su.se>, David van der Spoel
>>> <spoel at xray.bmc.uu.se>, michel.cuendet at epfl.ch, Michael Shirts
>>> <mrshirts at gmail.com>
>>>
>>>
>>>
>>>Hi,
>>>
>>>I have changed the kinetic energy to 0.5(Ekin(t-1/2) + Ekin(t+1/2))
>>>and changed the Berendsen and Nose-Hoover thermostat coupling
>>>temperatures to T(t-1/2).
>>>I also made Nose-Hoover leap-frog reversible.
>>>
>>>Everything seems to be running fine. Exact continuation also works.
>>>Tomorrow I go to Gerrit's promotion and next week I'll go to a conference.
>>>I'll run some simulations in the meantime.
>>>If everything is correct I will commit next Thursday.
>>>
>>>Up till now I have only looked at water simulations.
>>>The temperature difference for SPC at 300K with dt=2 fs is 2.2 K,
>>>meaning that all simulations were actually at 302.2 K.
>>>For dt=4 fs the difference is 9 K.
>>>
>>>Berk.
>>
--------------------------------------------------------------------------------------------------------------------------------------
> -------- Original Message --------
> Subject: Re: [Fwd: Re: results]
> Date: Mon, 08 Aug 2005 19:42:12 +0200
> From: Michel Cuendet <michel.cuendet at epfl.ch>
> Reply-To: michel.cuendet at epfl.ch
> To: David <spoel at xray.bmc.uu.se>
> CC: Erik Lindahl <lindahl at sbc.su.se>, mrshirts at gmail.com
> References: <1123513128.6473.6.camel at localhost.localdomain>
>
>
>
>Hi guys,
>
>First of all, sorry for the ugly tables, I hadn't realized that my
>compose window was wider than standard.
>
>I would just like to point out that there are two distinct problems :
>
>1) the integrator itself (which needs to be reversible, and "as close to
>symplectic as possible");
>
>2) the estimation of the temperature, which must be done with
> T(t) ~ v(t-1/2)^2 + v(1+1/2)^2
> rather than
> T(t) ~ [ v(t-1/2) + v(1+1/2) ]^2
>
>Of course (1) implies (2).
>
>Point (2), which we sorted out with Berk in mai, IS ALSO A PROBLEM WITH
>BERENDSEN THERMOSTAT. It could be seen as a bias in the feedback loop of
>scaling the velocities and measuring the resulting temperature. We now
>need to make sure that there is no similar "biases" in the pressure
>control loop in 3.3b.
>
>So in my opinion it's not "fine of course" ;-) with older versions, and
>maybe neither yet with 3.3b NPT.
>
>Best,
>
>Michel
>
>
--
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Michel Cuendet, Ph.D. student
Laboratory of Computational Biochemistry and Chemistry
Swiss Federal Institute of Technology in Lausanne (EPFL)
CH-1015 Lausanne
Switzerland Phone : +41 21 693 0324
lcbcpc21.epfl.ch/~mitch +41 1 633 4194
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-developers
mailing list