AMBER: Still puzzled about this virial computation

From: David Cerutti <>
Date: Thu, 1 Feb 2007 00:00:10 -0800 (PST)


    Thanks to Rob and Ross for their helpful inputs. I'm still puzzled a
bit about some of the computations, though.

I am using the following in my run file, but it appears to be ineffective
at removing the various long-range corrections (for which I do not yet
have an equivalent in my code) from SANDER:


1.) I can't seem to get my (real-space, no 1:4 interactions, no long-range
approximations, strict cutoff, periodic) van-der Waals energies to agree
with SANDER even when I use vdwmeth=0. This seems to be the case no
matter how many particles I use--2, 216, 1728. It's very bad with only
two particles--a difference of around 30%. I get better and better
agreement as I use a larger and larger cutoff with a bigger and bigger
system--presumably because I'm minimizing the effect of any long-range
correction, eventually obtaining that 1% discrepancy in the case of a
1728-particle system with an 18.0A cutoff.

2.) If I work hard to get the periodicity out of the system (by leaving
the configuration as is and making the box much larger so I can use a
200.0A cutoff, which is essentially doing an all:all calculation of a
non-periodic system), I get very good agreement (only 0.03% discrepancy,
which can be attributed to slight imprecision in my energy computer).

3.) In fact... I find that setting vdwmeth=0 has absolutely no effect on
the computed energy, even with a small 9.0A cutoff--it's the same as if I
do not specify a vdwmeth parameter, or explicitly set it to 1. Could this
be a bug?

4.) It seems that I get a consistent relationship between my computed
virial and the SANDER value, particularly in the case of that
"non-periodic" system--my computed virial is two times the SANDER value to
six significant digits. The system I'm working with is a collection of
monatomic particles, so I know there's no constraints contributing to the
virial. I'm not sure how I can get energies that essentially match and
yet have a virial that has some factor gone awry--which is why I really
need someone to explain the formula that AMBER uses when it computes the

5.) Like my ability to get the van-der Waals terms to agree under certain
conditions, I also get good agreement with the electrostatics under some
conditions--two water molecules in a big box or 1728 waters in a huge
box--essentially non-periodic systems. If I try to run a system where
periodicity would matter (one with a net charge or a neutral system
without a vast gulf separating the various pieces), I run into trouble.

6.) I don't know about constraints contributing to the virial, so to test
whether my code agrees with AMBER when electrostatics get into the
picture, I used a gas of "sodium" particles (I set the charges to 0.1e and
turned off the Lennard-Jones interactions to focus on the electrostatics
this time). When doing so, I finally get a hint as to what the formula
for computing that virial is: the virial comes out as negative 0.5 times
the energy; I would have expected it to be negative 1.0 times the energy,
so that's some explanation for where the factor of two in my computation
of the virial for a gas of Lennard-Jones particles came from. This
relationship between the SANDER-computed energy and virial remains if I
change the cutoff (which, as I would expect in a system without Ewald
terms, has a significant effect on the total energy computed as well as
the virial--but they vary in that 2:1 proportion).

7.) Why does the total energy seem to be greatly affected by some sort of
long-range electrostatic correction, even when I've set eedmeth=4 and
use_pme=0 (see point 5)?

There seem to be conundrums around every corner. Can anyone help
straighten this out?


The AMBER Mail Reflector
To post, send mail to
To unsubscribe, send "unsubscribe amber" to
Received on Sun Feb 04 2007 - 06:07:24 PST
Custom Search