Thanks very much for the reply, and the detailed explanation.
When my IMIN=5 values did not match the original energies, I had at
first feared that I had corrupted trajectory files, bad input files,
version-dependent differences, etc..
Going forward, I'll just eyeball the IMIN=5 scan and the original
numbers for rough similarity only.
At least, you have now be documented this on the reflector, for other
newcomers.
For AMBER 10, I doubt it is worth the hassle of keeping an extra set of
coordinate around.
As I mentioned in another thread, it would be nice if AMBER 10 allowed
IMIN=5 scanning of netcdf trajectory files.
Thanks
Chris
David A. Case wrote:
> 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.scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
>
>
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Wed Apr 11 2007 - 06:07:24 PDT
