AMBER Archive (2007)

Subject: Re: AMBER: SANDER PMEMD ntpr ntwx discrepancy (?bug? report)

From: David A. Case (case_at_scripps.edu)
Date: Fri Apr 06 2007 - 17:37:46 CDT


On Fri, Apr 06, 2007, Chris Moth wrote:

> When I set ntpr=1000 and ntwx=1000 in my sander MD input, I would expect
> to get energy readings, every ps, that corresponded to the coordinates
> that are written, every ps, to my coordinate output file.
>
> This is not the case. I believe there is an off-by-one error in sander
> (and pmemd).
>
> When I post-process the trajectory (source ./run_analyze.sh), I get
> energies that correspond to NSTEP=1001, 2001, 3001, etc. and NOT the
> energies at NSTEP=1000,2000,3000...
>

See the comments at line 670 of runmd.f. (I thought this was in the manual
someplace, but maybe not....)

Basically, this all arises from the leapfrog scheme that Amber uses to
integrate Newton's equations. In order to estimate the kinetic energy
at time t (and hence, print out energy values), one needs to know the
velocities at time t + (1/2)dt, where "dt" is the time step. In order to know
those velocities, one has to compute the coordinates at time t + dt. So, by
the time the energies are printed, the coordinates have been advanced by a
full time step.

There are lots of ways of getting around this problem, but they all involve
breaking backwards compatibility with earlier versions of the code. So it was
never done.

>
> While post-processing a trajectory with sander using imin=5, I found
> that a simple restraint energy differed significantly between my
> original ntpr=1000 energy outputs, and my new post-processing calculations.

As you note, one "work-around" is to ignore the values that are printed
in the mdout file, and use imin=5 to re-process the snapshots in the trj
file.

Probably better would be to modify the codes to save the coordinates
from one earlier step, and send those to the trajectory file, so that they
would then be in-sync with the printed energies. We'll think about making
this sort of change for Amber 10 (comments are welcome here). We might also
actually move to something like velocity Verlet, where the velocities and
coordinates are kept synchronized, but that's an even bigger departure from
past practice.

...regards...dac

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu