[gmx-developers] do_update_md

[Zachmann, Martin]

I am referring to the original Berendsen paper (The Journal of Chemical Physics 81, 3684 (1984); doi: 10.1063/1.448118).


     "In any algorithm, you should scale the velocity at the exact same time you computed the kinetic energy for. That is what the code since 4.6 does, right?"

Yes, since 4.5. scaling is done like this.


What i meant was that in this paper scaling is NOT in this way, but rather v(t+dt/2) is scaled with lambda(t-dt/2), noting:

"Note. Although lambda is based on the temperature at t - dt/2, its

value can be used to scale the velocity at t + dt/2 because of the slow variation of lambda"

Thanks for your help, martin

________________________________
From: gromacs.org_gmx-developers-bounces at maillist.sys.kth.se <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on behalf of Berk Hess <hess at kth.se>
Sent: Friday, June 9, 2017 3:38:30 PM
To: gmx-developers at gromacs.org
Subject: Re: [gmx-developers] do_update_md

Which paper are you referring to? And what do you mean with current v?
In any algorithm, you should scale the velocity at the exact same time you computed the kinetic energy for. That is what the code since 4.6 does, right?
But note that as tau is usually much larger than dt, the difference will be extremely small. But the v-rescale thermostat is implemented such that one can use tau=dt, or even tau=0, and get a fully correct ensemble.

Cheers,

Berk

On 2017-06-09 15:10, Zachmann, Martin wrote:

Thank you so much for your fast answers.

So you are saying the difference from the code to the original paper (as is now also the case in the current master-branch
in updateMdLeapfrogSimple()) is on purpose?
It says in the paper that scaling the current v with a lambda calculated from minus half-step v is incorrect, yet
acceptable since lambda changes slowly.
regards Martin

________________________________
From: gromacs.org_gmx-developers-bounces at maillist.sys.kth.se<mailto:gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se><mailto:gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on behalf of Berk Hess <hess at kth.se><mailto:hess at kth.se>
Sent: Friday, June 9, 2017 2:51:46 PM
To: gmx-developers at gromacs.org<mailto:gmx-developers at gromacs.org>
Subject: Re: [gmx-developers] do_update_md

On 2017-06-09 14:29, Mark Abraham wrote:
Hi,

That looks like a bug that was silently fixed during the adoption of the velocity-rescale thermostat of Bussi (see commit 333641e040ae824badf9d9531290a46ad581ed53). It may be that the issue wasn't noticed because people often used parameters that produce lg == 1...
Is debatable if one should call that a bug. The original Berendsen thermostat that this code was written for did not produce the correct ensemble anyhow, so it was likely purposely coded this way to give better/more scaling with consistently too high temperature, as would happen with the group cut-off scheme and no pair list buffering. With the V-rescale thermostat you should scale the velocity you computed the temperature for, which are the minus half-step velocities, as is now correctly done.
In practice it will be difficult to detect a difference between these two codes.

I would strongly advise avoiding doing development work on unsupported versions of the software more than five years old... ;-) Get the git master branch and stay up to date with the hundreds of bugs that have been fixed since :-)
His mail didn't imply he is using an old version for development, it could just be for trying to understand the code or porting changes from an old version to master, which is good!

Cheers,

Berk

Mark

On Fri, Jun 9, 2017 at 2:01 PM Zachmann, Martin <m.zachmann at tum.de<mailto:m.zachmann at tum.de>> wrote:

Hello community!

i hope someone can help me understand some code changes in the do_update_md() function in update.cpp.

In the case where it advances r and v in the leap-frog scheme together with the berendsen thermostat

(last two else if { } and else { }) in version 4.0.7. it looks like this

vv = lg*(vn + f[n][d]*w_dt);

where lg is the lambda-scaling factor from the thermostat and w_dt is dt/m.

Instead in version 4.5. and higher:

 vv = lg*vn + f[n][d]*w_dt;


The original paper (The Journal of Chemical Physics 81, 3684 (1984); doi: 10.1063/1.448118) seems to scale

both old velocity and force. Btw in do_update_visc() scaling is done to both.

What is the purpose of this?

any help is much appreciated,

regards Martin


--
Gromacs Developers mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers or send a mail to gmx-developers-request at gromacs.org<mailto:gmx-developers-request at gromacs.org>.






-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170609/501317eb/attachment-0001.html>