[gmx-developers] Scaling velocities in Velocity Verlet

[Shirts, Michael (mrs5pt)]

Hi, David-

If you can commit your changes, then I can check out qhop and see what's
happening.

Best,
~~~~~~~~~~~~
Michael Shirts
Assistant Professor
Department of Chemical Engineering
University of Virginia
michael.shirts at virginia.edu
(434)-243-1821


> From: David van der Spoel <spoel at xray.bmc.uu.se>
> Reply-To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Date: Sun, 4 Dec 2011 09:28:03 +0100
> To: Discussion list for GROMACS development <gmx-developers at gromacs.org>
> Subject: Re: [gmx-developers] Scaling velocities in Velocity Verlet
> 
> On 12/4/11 4:35 AM, Shirts, Michael (mrs5pt) wrote:
>> Hmm.  Which version/branch is this?
> This is the qhop branch, based on the 4.5 branch.
> 
> 
>> 
>> Best,
>> ~~~~~~~~~~~~
>> Michael Shirts
>> Assistant Professor
>> Department of Chemical Engineering
>> University of Virginia
>> michael.shirts at virginia.edu
>> (434)-243-1821
>> 
>> 
>>> From: David van der Spoel<spoel at xray.bmc.uu.se>
>>> Reply-To: Discussion list for GROMACS
>>> development<gmx-developers at gromacs.org>
>>> Date: Sat, 3 Dec 2011 22:10:57 +0100
>>> To: Discussion list for GROMACS development<gmx-developers at gromacs.org>
>>> Subject: Re: [gmx-developers] Scaling velocities in Velocity Verlet
>>> 
>>> On 2011-12-03 21:37, Shirts, Michael (mrs5pt) wrote:
>>>> Hi, all-
>>>> 
>>>>> The problem is that the change in Ekin comes one step too late. I tried
>>>>> to deduce the algorithm from the code (4.5 code base) and deduced the
>>>>> following MD algorithm:
>>>> 
>>>> It's rather difficult to follow currently because of the need to bookkeep
>>>> both md and md-vv + plus correct pressure control w/constraints. I have
>>>> high
>>>> hopes for improving this for 5.0.
>>>> 
>>>>> v(t) = v(t-dt/2) + a(t) dt/2
>>>>> 
>>>>> v'(t) = lambda v(t)<======= my addition
>>>>> 
>>>>> v(t+dt/2) = v'(t) + a(t) dt/2
>>>>> 
>>>>> r(t+dt) = r(t) + v(t+dt/2) dt
>>>>> 
>>>>> However, the Ekin that is printed is the uncorrected one even though I
>>>>> scale both ekind->ekin and ekind->tcstat[x].ekinf .
>>>>> 
>>>>> Is the algorithm correct as I write it above? Is the Ekin(t) already
>>>>> stored at the beginning of the time step?
>>>> 
>>>> If you add your modification as an option to update_tcouple, and operate on
>>>> ekind directly (not the velocities), in the same way that berendsen_tcouple
>>>> does, then the bookkeeping should take care of the velocity rescaling
>>>> appropriately.
>>> Close, but no cigar.
>>> Here's the output from g_energy:
>>> 
>>> @ s0 legend "Potential"
>>> @ s1 legend "Kinetic En."
>>> @ s2 legend "Total Energy"
>>>       0.029600  -35405.431557  6460.241804  -28945.189753
>>>       0.029800  -35407.596247  6462.418642  -28945.177606
>>>       0.030000  -35327.293038  6466.347370  -28860.945667
>>>       0.030200  -35330.858001  6387.271571  -28943.586430
>>>       0.030400  -35334.701690  6391.114750  -28943.586940
>>> 
>>> The jump in Potential is at t = 0.03, in Kinetic it is at t = 0.302. As
>>> you see the Total energy has a peak at t = 0.03, but it does come back
>>> almost to the right level at the next step. This all makes me think I
>>> should be scaling the previous velocities, hence my questions about the
>>> algorithm.
>>> 
>>> 
>>> 
>>>> 
>>>> If this _doesn't_ work, then there's some tweaks that got lost in the
>>>> different versions that are floating around, and I can either put the
>>>> changes into 4.6 alone, or just add them in with the 4.6 merge.
>>>> 
>>>> Best,
>>>> ~~~~~~~~~~~~
>>>> Michael Shirts
>>>> Assistant Professor
>>>> Department of Chemical Engineering
>>>> University of Virginia
>>>> michael.shirts at virginia.edu
>>>> (434)-243-1821
>>>> 
>>>> 
>>>> 
>>>>> v(t) = v(t-dt/2) + a(t) dt/2
>>>>> 
>>>>> v'(t) = lambda v(t)<======= my addition
>>>>> 
>>>>> v(t+dt/2) = v'(t) + a(t) dt/2
>>>>> 
>>>>> r(t+dt) = r(t) + v(t+dt/2) dt
>>>>> 
>>>>> However, the Ekin that is printed is the uncorrected one even though I
>>>>> scale both ekind->ekin and ekind->tcstat[x].ekinf .
>>>>> 
>>>>> Is the algorithm correct as I write it above? Is the Ekin(t) already
>>>>> stored at the beginning of the time step?
>>>>> 
>>>>> 
>>>>> --
>>>>> 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    http://folding.bmc.uu.se
>>>>> --
>>>>> gmx-developers mailing list
>>>>> gmx-developers at gromacs.org
>>>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>>>> Please don't post (un)subscribe requests to the list. Use the
>>>>> www interface or send it to gmx-developers-request at gromacs.org.
>>>> 
>>> 
>>> 
>>> --
>>> 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    http://folding.bmc.uu.se
>>> --
>>> gmx-developers mailing list
>>> gmx-developers at gromacs.org
>>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>>> Please don't post (un)subscribe requests to the list. Use the
>>> www interface or send it to gmx-developers-request at gromacs.org.
>> 
> 
> 
> -- 
> David.
> ________________________________________________________________________
> David van der Spoel, PhD, Professor of Biology
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,   75124 Uppsala, Sweden
> phone: 46 18 471 4205  fax: 46 18 511 755
> spoel at xray.bmc.uu.se spoel at gromacs.org   http://folding.bmc.uu.se
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> -- 
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.