[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 


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

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

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