[gmx-developers] Re: Integrators

[Berk Hess]

Hi,

Constraining velocities with any algorithm (shake, (p-)lincs or settle) 
is easy now
there is just one constrain call that does everything.
So it should be easy to implement the iterative algorithm.
We can add an mdp option to turn off iteration if the user does not want 
this.
In parallel the performance hit might not be too bad, since for P-Lincs only
one communication step is required.
The question is how accurate lincs or shake should be to really see an
accuracy improvement due to the iteration.

Berk

Michael Shirts wrote:
> Well, unfortunately, it looks like the case with constraints is more
> complicated than I had initially thought.  The issue is that there's
> coupling between the Trotter steps and the constraint steps, which
> means there needs to be iteration over the constraints.  Fortunately,
> this is discussed in Appendix D of the 1996
> Martyna-Tuckerman-Tobias-Klein paper.
>
> The procedure presented there goes like this (using a = F/m, and h = \delta t).
>
> 1)  Trotter half-step [using P(t-1) and K(t-1)], using
>     the Lagrangian multipliers from RATTLE of the previous step.
> 2)  r(t+h)  = r(t) + h*v(t) + h^2/2*a(t)
> 3)  v(t+h/2) = v(t) + h/2*a(t)
> 4)  Constrain r(t+h), and recalculate the Lagrangian multipliers
> 5)  Iterate 1-4
> 6)  Calculate a(t+h)
> 7)  v(t+h)  = v(t+h/2) + h/2*a(t)
> 8)  Second Trotter half-step, using the P(t) with the forces from the Lagrangian
>     multiplier obtained after step 5.
> 9)  Constrain v(t+h) using RATTLE, to determine the full step velocities.
> 10) Repeat steps 7-9 until the constraints are satisfied, and the box
>     speed is consistent.
>
> Note that the even though the constraints can be satisfied with one
> step, the Trotter step is no longer correct in this case, because it
> requires the velocities -- which were not correctly constrainted when
> we started.  So we must iterate.
>
> Translating this into the current GROMACS bookkeeping, we get:
>
> 1)  Calculate a(t)
> 2)  v(t)   = v(t-h/2) + h/2*a(t)
> 3)  A Trotter half-step, using the P with the forces from the Lagrangian
>     multiplier obtained after step 5.
> 4)  Constrain v(t+h) using RATTLE, to determine v(t).
> 5)  Repeat steps 2-4 until the constraints are satisfied, and the box
> speed is consistent.
> 6)  Output r(t),v(t),a(t)
> 7)  Another Trotter half-step, using the Lagrangian multipliers from
> RATTLE of the previous step.
> 8)  r(t+h)  = r(t) + h*v(t) + h^2/2*a(t)
> 9)  v(t+h/2) = v(t) + h/2*a(t)
> 10) Constrain r(t+h), and recalculate the Lagrangian multipliers.
> 11) Iterate 7-10.
> 12) Output P(t), K(t), and U(t).
>
> This is obviously a bit more complicated.  We must do several
> iterations over the constraints, and several summations of the kinetic
> energy each time we do the constraints, leading to a lot of communication.
>
> I also still need to translate this into the constraint formalism used
> in Gromacs currently, where the Lagrangian multipliers are not
> specifically stored.
>
> One thing that might not totally kill the scaling is that it -looks-
> like we still don't have to pressure / temperature couple after every
> step.  Perfectly legitimate dynamics can be constructed by only
> pressure coupling every N steps, though the bounds to the
> approximation needs to be tested pretty strictly -- it could be that
> it starts having issues if it's more than 2-3 steps at a time between
> pressure "bumps".  Kinetic energy and multiple constraint iterations
> don't need to be done in between.
>
> Any thoughts as to what the right order of putting things might be?  I
> think it's worth getting this in, even if it's not going to be the
> default for large scale simulations.   But it will take a bit of time.
>  I can probably go ahead and get everything except the iteration over
> the pressure coupling pretty easily -- it's pretty much in, I just
> need to tie together the loose ends I left by trying to get
> constraints working and debug a bit.
>
> Anything beyond that will take a bit longer. . .
>
> Cheers,
> Michael
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.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.
>