[gmx-users] Velocities with constraints

[Miguel Caro]

Thanks a lot! I will post the solution if/when I manage it, in case it's 
useful for others.

Miguel

On 2015-11-12 16:19, Mark Abraham wrote:
> Hi,
>
> On Thu, Nov 12, 2015 at 2:49 PM Miguel Caro <miguel.caro at aalto.fi> wrote:
>
>> Dear Mark,
>>
>> I have rechecked that my equations should actually be correct: I am
>> calculating the angular velocity in the reference frame where the center
>> of mass is at rest, therefore for a rigid body the rotational axis in
>> this frame should go through the center of mass.
>>
>> I have done further analysis of the results and have found very
>> systematic errors in the angle between the velocity and the position
>> vector of the atoms belonging to the same molecule. This angle should
>> always be 90 degrees for rotational motion, however Gromacs is giving
>> slightly less (typically about 89.5 degrees). I have checked the effect
>> of single/double precision to be negligible here. My interpretation of
>> this result is that Gromacs is *estimating* what the velocity vector
>> will look like some time after the instant in which the positions are
>> calculated. If you try to picture this, I think Gromacs is printing r(t)
>> and v(t+dt) on the same line for each atom.
>>
>> Is this what is going on? If so, what is the value of dt and how can I
>> prevent this behavior? If dt is the same as the integration time I could
>> then load the positions at time = t together with the velocities
>> *printed* at t-dt. Does this make sense?
>>
> The usual integrator in GROMACS is leapfrog, which by nature computes r and
> v at different times - they are out of phase by half of the step size
> (though I forget which of them leads in GROMACS). This is reflected in the
> r and v output, because r and v are never known at the same time
> (KE-related quantities are handled specially). Estimating temperature in MD
> is a very tricky business (Eastwood, 2009)
>
> Ideally I would like to induce the behavior that r(t) and v(t) are
>> printed on the same line.
>>
> There is a velocity-Verlet integrator available, which does compute r and v
> at the same time. I suspect it writes similarly staggered velocities, but
> the integrator implementation is diabolically complex and I've never needed
> to be sure... That's something you could try out.
>
> Mark
>
> Regards,
>> Miguel
>>
>> On 2015-11-12 00:05, Mark Abraham wrote:
>>> Hi,
>>>
>>> Sounds right. Otherwise, the algorithm involved is SETTLE or RATTLE,
>> which
>>> has an analytic solution, not LINCS.
>>>
>>> Mark
>>>
>>> On Wed, Nov 11, 2015 at 10:31 PM Caro Miguel <miguel.caro at aalto.fi>
>> wrote:
>>>> I'm realizing now that the equations I'm using assume rotation about an
>>>> axis that goes through the center of mass, which needs not be the case
>>>> during the dynamics. This is probably the origin of the problem.
>>>>
>>>> Regards,
>>>> Miguel
>>>>
>>>> --
>>>> Dr. Miguel Caro
>>>> Postdoctoral researcher
>>>> Department of Electrical Engineering and Automation, and
>>>> COMP Centre of Excellence in Computational Nanoscience
>>>> Aalto University, Finland
>>>> Email: mcaroba at gmail.com
>>>> Work: miguel.caro at aalto.fi
>>>> Website: http://mcaroba.dyndns.org
>>>> ________________________________________
>>>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [
>>>> gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Miguel
>>>> Caro [miguel.caro at aalto.fi]
>>>> Sent: 11 November 2015 20:37
>>>> To: gmx-users at gromacs.org
>>>> Subject: [gmx-users] Velocities with constraints
>>>>
>>>> Dear Gromacs users,
>>>>
>>>> I am using the rigid water model TIP3P and running some dynamics. I need
>>>> very accurate velocities from the output trajectory. When I do frequency
>>>> decomposition of these velocities, I see that there is a non-negligible
>>>> vibrational component to them, that is, the atomic velocities of the
>>>> molecules are not following exactly the constraints given for a rigid
>>>> body. What I mean by this is that the velocities are not fully
>>>> consistent with what should be expected from a rigid body rotating about
>>>> an instantaneous rotation axis (plus a center of mass translation).
>>>>
>>>> I have trying debugging my code but haven't been able to find any
>>>> mistakes, hence my question: is this behavior to be expected from
>>>> Gromacs and can it be influenced for instance by tuning the LINCS
>>>> parameters?
>>>>
>>>> Many thanks for your help,
>>>> Miguel
>>>> --
>>>> *Dr. Miguel Caro*
>>>> /Postdoctoral researcher/
>>>> Department of Electrical Engineering and Automation,
>>>> and COMP Centre of Excellence in Computational Nanoscience
>>>> Aalto University, Finland
>>>> Personal email: *mcaroba at gmail.com*
>>>> Work: *miguel.caro at aalto.fi*
>>>> Website: http://mcaroba.dyndns.org
>>>> --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>> --
>>>> Gromacs Users mailing list
>>>>
>>>> * Please search the archive at
>>>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
>>>> send a mail to gmx-users-request at gromacs.org.
>>>>
>> --
>> *Dr. Miguel Caro*
>> /Postdoctoral researcher/
>> Department of Electrical Engineering and Automation,
>> and COMP Centre of Excellence in Computational Nanoscience
>> Aalto University, Finland
>> Personal email: *mcaroba at gmail.com*
>> Work: *miguel.caro at aalto.fi*
>> Website: http://mcaroba.dyndns.org
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
>> send a mail to gmx-users-request at gromacs.org.
>>

-- 
*Dr. Miguel Caro*
/Postdoctoral researcher/
Department of Electrical Engineering and Automation,
and COMP Centre of Excellence in Computational Nanoscience
Aalto University, Finland
Personal email: *mcaroba at gmail.com*
Work: *miguel.caro at aalto.fi*
Website: http://mcaroba.dyndns.org