Hello,
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:
&ewald
use_pme=0,
eedmeth=4,
vdwmeth=1,
&end
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
virial.
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?
Thanks!
Dave
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Sun Feb 04 2007 - 06:07:24 PST
