[gmx-users] Velocities with constraints

[Mark Abraham]

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