[gmx-developers] Scaling velocities in Velocity Verlet

[David van der Spoel]

On 2011-12-04 14:42, Shirts, Michael (mrs5pt) wrote:
> Hi, David-
>
> If you can commit your changes, then I can check out qhop and see what's
> happening.
>
Done.
I'll send you an input file off-list.

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


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