[gmx-developers] Verlet Algorithm
Michael R Shirts
Michael.Shirts at Colorado.EDU
Thu Mar 2 03:03:46 CET 2017
Generally, velocity verlet is performed the following way (from wikipedia https://en.wikipedia.org/wiki/Verlet_integration)
1. Calculate v(t+dt/2) = v(t) + dt/2 a(t)
2. Calculate x(t+dt) = x(t) + dt*v(t+dt/2)
3. Calculate a(t+dt) from x(t+dt)
4. Calculate v(t+dt) = v(t+dt/2) + dt/2 a(t+dt)
After the first step, it’s only one force call per step because you can reuse the force in step 4 in the next step 1.
The logic is not always clear in the gromacs code, since the bookkeeping needs to be rather complicated to do both leapfrog and velocity.
Best,
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Michael Shirts
Associate Professor
michael.shirts at colorado.edu<mailto:michael.shirts at colorado.edu>
http://www.colorado.edu/lab/shirtsgroup/
Phone: (303) 735-7860
Office: JSCBB C123
Department of Chemical and Biological Engineering
University of Colorado Boulder
From: <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on behalf of Elena della Valle <elena.dv46 at yahoo.it>
Reply-To: "gmx-developers at gromacs.org" <gmx-developers at gromacs.org>
Date: Wednesday, March 1, 2017 at 12:19 PM
To: "gromacs.org_gmx-developers at maillist.sys.kth.se" <gromacs.org_gmx-developers at maillist.sys.kth.se>
Subject: Re: [gmx-developers] Verlet Algorithm
Hi, sorry if I bother you again with the velocity verlet algorithm, but I have another question, in the Gromacs Manual for the Velocity verlet algorithm for the position is written that r (t + deltaT) is equal to:[cid:image001.png at 01D292CF.4A8C25F0]
But if I look in the gromacs code, in the update.cpp file I see that the position for the verlet algorithm the update of the positions is due as:
[cid:image002.png at 01D292CF.4A8C25F0]
and I don't see the therm of the Force multiplied for deltaT^2/2*m, is because in the code the therm v[n] is interpreted as v(t + delta)?
thanks in advance
Best Regards
Elena della Valle
Il 06/10/2016 12:03, Berk Hess ha scritto:
do_update_vv_vel() is called twice during each step at different points. Thus the updates for F(t) and F(t+dt) are done separately.
Berk
On 2016-10-06 11:53, Elena della Valle wrote:
The question is that on the gromacs manual I read that the verlet algorithm for the velocity corresponds to [cid:image003.png at 01D292CF.4A8C25F0] , now in the gromacs code in the verlet part I don't find the corresponds for the term [1/2*(∆t/m)* F(t + ∆t)], is not implemented?
Sorry if I'm not so clear
Thanks
Elena della Valle
Il 06/10/2016 11:39, Berk Hess ha scritto:
I don't understand our question.
accel[ga] is a fixed acceleration per atom for non-equilibrium dynamics. This will be 0, unless you are doing non-equilibrium dynamics. You can ignore or remove this term, since you will not need it.
Cheers,
Berk
On 2016-10-06 11:25, Elena della Valle wrote:
Okey but in the update.c file line 300, in the following equation : v[n][d] = mv1*(mv1*v[n][d] + 0.5*(w_dt*mv2*f[n][d]))+0.5*accel[ga][d]*dt;
with the term "0.5*accel[ga][d]*dt" it corresponds to the term "[1/2*(∆t/m)* F(t + ∆t)] " in the verlet algorithm ?
Thanks in advance for helping me
Best Regards
Elena della Valle
Il 06/10/2016 09:30, Berk Hess ha scritto:
Hi,
I guess you mean ga = cACC[n].
cACC[n] contains the acceleration group index for atom n. This is for non-equilibrium dynamics where groups of atom are accelerated.
(apologies for the cryptic notation without documentation)
Cheers,
Berk
On 2016-10-06 07:44, Elena Della valle wrote:
Hi, I have just a last question what gromacs means exactly with the term: cACC[ga]?
Thanks in advance
Best Regards
Elena della Valle
Il giorno 04 ott 2016, alle ore 17:58, Berk Hess <hess at kth.se><mailto:hess at kth.se> ha scritto:
On 10/04/2016 05:36 PM, David van der Spoel wrote:
On 04/10/16 17:00, Elena Della valle wrote:
I used the gromacs 4.6.5
Best regards
Elena della Valle
Crudely put we do not support that version anymore, so if you would like to have this considered for inclusion in gromacs you would have to implement it in the master branch and upload it to gerrit. Whether it is correct you should test your self using some known cases and carefully constructed inputs.
In addition, even when you are not interested in contributing the code, all developer work with the most recent versions, so it's easier for us to help you with the master version, especially for code that has changed a lot over the past years like the Velocity Verlet integrator. (And I think anyone should prefer the latest version for development, so one doesn't end up with functionality that only works in 5 year old code.)
Cheers,
Berk
Maybe you can also answer the question when a magnetic effect would be important?
Il giorno 04 ott 2016, alle ore 16:49, David van der Spoel <spoel at xray.bmc.uu.se><mailto:spoel at xray.bmc.uu.se> ha scritto:
On 04/10/16 15:15, Elena della Valle wrote:
Hi to all,
I'm Elena della Valle, a Ph.D. student coming from la Sapienza
University of Rome.
I'm writing because I have some questions about some modifications that
I have done to the update of velocities and positions in the verlet
algorithm. My aim is to implement the magnetic field in gromacs by
introducing the therm of the larmor frequency to the velocities and
positions. By literature I read that this has been done by the update of
the verlet velocities and positions with the frequency larmor therm.
These are the equations in order to be more clear on what i wanted to do:
I did this by modifying in the update.c file in gromacs the
update_do_vv_vel and update_do_vv_pos ad follows:
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
v[n][0] = mv1*(mv1*v[n][0] +
0.5*(w_dt*mv2*f[n][0]))+0.5*accel[ga][0]*dt +
w_dt*charge[n]*campoB*v[n][1]*mv1*mv1 +
0.25*dt*invmass[n]*charge[n]*campoB*((w_dt*mv2*mv1*f[n][1]) +
accel[ga][1]*dt -2*w_dt*charge[n]*campoB*mv1*mv1*v[n][0]);
v[n][1] = mv1*(mv1*v[n][1] +
0.5*(w_dt*mv2*f[n][1]))+0.5*accel[ga][1]*dt +
w_dt*charge[n]*campoB*v[n][0]*mv1*mv1 -
0.25*dt*invmass[n]*charge[n]*campoB*((w_dt*mv2*mv1*f[n][0]) +
accel[ga][0]*dt +2*w_dt*charge[n]*campoB*mv1*mv1*v[n][1]);
v[n][2] = mv1*(mv1*v[n][2] +
0.5*(w_dt*mv2*f[n][2]))+0.5*accel[ga][2]*dt;
//printf("frame %lf \n", mv1);
// printf("frame %d: %lf\t%lf%lf\n", n, v[n][0], v[n][1],
v[n][2]);
}
in the do_update_vv_pos:
if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
{
xprime[n][0] = mr1*(mr1*x[n][0]+mr2*dt*v[n][0]) +
mr1*0.5*dt*(w_dt*mr2*f[n][0] + w_dt*charge[n]*campoB*v[n][1]*mr1);
xprime[n][1] = mr1*(mr1*x[n][1]+mr2*dt*v[n][1]) +
mr1*0.5*dt*(w_dt*mr2*f[n][1] - w_dt*charge[n]*campoB*v[n][0]*mr1);
xprime[n][2] = mr1*(mr1*x[n][2]+mr2*dt*v[n][2]) +
mr1*0.5*dt*(w_dt*mr2*f[n][2]);
}
After that in the grompp file I modified the integrator in md-vv and I
run the simulation and it seems to use the position and velocities
modifications.
I want know if you think that this way is correct, if I implemented in
the right way
And Also I would like to know if in the grompp file I have to modify
other parameters (I use a beredensed coupling)
Thanks in advance
Sorry for bothering you
Best Regards
Elena della Valle
--
Elena della Valle
Ph.D. Student in Electronic Engineering
Department of Information Engineering, Electronics and Telecommunications
Sapienza, University of Rome
via Eudossiana, 18 00184 Rome
Just a quick question: which version did you use? All development should be in the master version and this looks like uses something older.
--
David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205.
spoel at xray.bmc.uu.se<mailto:spoel at xray.bmc.uu.se> http://folding.bmc.uu.se
--
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>.
--
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>.
--
Elena della Valle
Ph.D. Student in Electronic Engineering
Department of Information Engineering, Electronics and Telecommunications
Sapienza, University of Rome
via Eudossiana, 18 00184 Rome
--
Elena della Valle
Ph.D. Student in Electronic Engineering
Department of Information Engineering, Electronics and Telecommunications
Sapienza, University of Rome
via Eudossiana, 18 00184 Rome
--
Elena della Valle
Ph.D. Student in Electronic Engineering
Department of Information Engineering, Electronics and Telecommunications
Sapienza, University of Rome
via Eudossiana, 18 00184 Rome
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170302/741aaf69/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 5608 bytes
Desc: image001.png
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170302/741aaf69/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image002.png
Type: image/png
Size: 9374 bytes
Desc: image002.png
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170302/741aaf69/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image003.png
Type: image/png
Size: 5257 bytes
Desc: image003.png
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20170302/741aaf69/attachment-0005.png>
More information about the gromacs.org_gmx-developers
mailing list