AMBER Archive (2007)

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

From: Chris Moth (chris.moth_at_vanderbilt.edu)
Date: Fri Apr 06 2007 - 16:07:17 CDT

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). If this is not easy to correct, it might be nice to
address this in the manual.

How I found the problem:

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.

Further analysis of the precise coordinates of the restrained atoms
proved that that the post-processing calculations were correct. I.e.,
the original energy outputs from sander where not exactly reflective of
the output coordinate sets.

To duplicate the "bug" for you, I setup a minimal 2 "Xenon" (zero charge
spheres with artifical van der Waals radii and epsilon parameters) atom
simulation to probe this farther, and I have placed the files on my
website at:

http://structbio.vanderbilt.edu/~mothcw/amber_bug/

In the simulation (source ./run_md.sh) , I have set ntpr=1 and
ntwx=1000. This way, I get energy output every fs, in addition to the
outputs I would normally see at NSTEP=1000,2000,3000 etc.

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

Specifically, here are the energies from post_processing 10ps (10
frames) of trajectory coordinates. I've checked the coordinates in a
spreadsheet - and they energies are correct for them:

socrates 109> grep VDWAALS analyze.out | uniq
VDWAALS = -0.2924 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.2650 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.1436 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.2032 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.1089 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.0402 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.0626 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.4988 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.2807 EEL = 0.0000 EGB =
0.0000
VDWAALS = -0.4684 EEL = 0.0000 EGB =
0.0000

Below are the energies at each step NSTEP=1001, 2001, 3001, 4001 from
the original simulation. These are correct to within the rounding errors
caused by the truncated ASCII output in the coordinate file. Note they
MATCH the post-processing energies.... but I believe they should NOT match.

socrates 115> grep "NSTEP = .*001 .*TIME" --after=5 md.out | grep VDWAALS
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2924
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2652 (larget variance - off by .0002 rounding error)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.1436
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2033
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.1088
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.0402
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.0627
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.4988
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2807
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.4683

Below are the energies at each step NSTEP = 1000,2000,3000,4000 from
the original simulation.
These are the energies I would see if I simply had ntpr=1000.

I believe that these should match the energies of the coordinate sets
output with ntwx=1000. But, they diverge considerably more than the
energies at N001.

socrates 116> grep "NSTEP = .*000 .*TIME" --after=5 md.out | grep VDWAALS
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2922 (off by .0002)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2658 (off by .0008)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.1433 (off by .0003)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2035 (off by .0003)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.1089
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.0402
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.0625
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.4989
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.2812 (off by .0005)
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.4651 (off by .0033)

Apologies that the numbers above are small. With only 2 atoms, I can
restrain one, but I cannot remove angular momentum. So, I suspect
temperature is manifesting itself as spinning.

Perhaps I could have been smarter with this minimal simulation. In
larger systems, I am seeing energies vary by as much as hundreds of
kCals. If you would like me to change the demo simulation in some way,
I am happy to do so.

Thanks for looking into this. It is not the end-of-the-world if it
can't be easily fixed. But, folks should be aware of the discrepancy,
and its source.

Chris

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