[gmx-users] Time reversibility and settle

[Francesco Mercuri]

> > Hello and thank you for your reply.
> > However, my question about the time reversibility of a leap-frog integrator
> > with settle (e.g. for TIP3P water molecules) was concerned with the
> > actual implementation in Gromacs, rather than the analytical point
> > of view. In other words, is this implementation of leap-frog +
> > settle still time-reversible in the limit of infinite numerical
> > accuracy?
> > This question arises since I tried different numerical
> > accuracies (single, double and quadruple precision numbers and
> > operations; for the latter I had to rewrite large parts of the
> > code...) with, in all cases, similar (and pretty large)
> > deviations from the "forward" trajectory when reversing the time.
> > On the other hand, just bypassing the "settle" subroutine
> > (e.g. by commenting all the code related to the original
> > implementation of settle by Miyamoto and Kollman), it works as
> > expected: the algorithm is almost perfectly time-reversible
> > and the only "noise" is due to numerical inaccuracies, with
> > deviations comparable, in the three different cases (single,
> > double, long double) to the accuracy of numbers.
> > Thus, I'm trying to understand whether just resetting the position
> > of atoms, as done in the current implementation of "settle",
> > still leads to a time-reversible MD (in the limit of infinite
> > machine accuracy).
> > Any suggestion about that?
> 
> In the limit is should be time-reversible, unless there is a bug in
> the settle algorithm, which I consider very unlikely.
> Have you tried to use shake or lincs (with high precision settings)?
> 
> There could be some other issues, but when the unconstrained
> dynamics is reversible you have probably taken care of these:
> You should not use temperature and pressure coupling (the new
> version will have reversible Nose-Hoover coupling).
> When reversing you should take care to use the proper velocities,
> as x(t), v(t-t/2dt) is stored you need to take the velocities from
> the next step when reversing.
> 
> Berk.

Hi again and thanks for your suggestions.
Still in the search of better time-reversibility properties, I'd like
to ask a few more questions.
If I understood correctly, the leap-frog + settle algorithm is implemented
in Gromacs according to the following steps:

v'(t+dt/2) = v(t-dt/2) + dt*f(t)/m   ("uncorrected" velocities)
r'(t+dt) = r(t) + dt*v'(t+dt/2)      (unconstrained move)
r(t+dt) = settle[r'(t+dt)]           (reset of atomic position calling "csettle")
v(t+dt/2) = [r(t+dt)-r(t)]/dt        (corrected velocities)

is that correct?
In this case, is not very clear to me how is it possible to "reverse" it,
e.g. restarting from a [r(t),v(t+dt/2)] configuration and inverting the
sign of dt, since the settle procedure modifies the atomic position in a
way that is, apparently, not "reversible".
Is any modification of the velocities needed?

Another question concerns the preparation of restart files containing
configurations like the one cited above for time-reversing purposes,
i.e. [r(t),v(t+dt/2)] (instead of the default [r(t),v(t-dt/2)] ).
I tried the quick and dirty way by just running one more MD step,
converting the .trr file in ascii format, manually editing the
ascii file in order to get the [r(t),v(t+dt/2)] configuration,
pasting into a .gro file and converting back the .gro file into
a .trr file for restarting. However, even if I modified the routine
pr_rvecs in txtdump.c in order to have more decimal digits in the
ascii files, the accuracy seems to be limited to the order of
1.e-8, whereas I expected a better accuracy for double precision
numbers. Is there any better way to do that? Which is the accuracy
with which numbers are stored in the .trr file (with double
precision compilation)?
Thank you again and best regards,

             Francesco

-- 
___________________________________________________________
Sign-up for Ads Free at Mail.com
http://promo.mail.com/adsfreejump.htm