A few comments embedded between ==> and <==.
Regards - Bob
----- Original Message -----
From: "David Cerutti" <dcerutti.mccammon.ucsd.edu>
To: <amber.scripps.edu>
Sent: Thursday, February 01, 2007 3:00 AM
Subject: AMBER: Still puzzled about this virial computation
> 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
>
==>
There are no long range corrections here. See below.
<==
> 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.
==>
Are you sure the vdw calc is correct? You don't say what your system is
here, but only "two particles" I presume means a couple of atoms or ions.
If you really want to do this stuff you are going to have to take the sander
code and dump info with writes from a single processor run of your model
system. For some pairs, there is no vdw (tip3p water H being a good
example). There are also 1-2, 1-3, 1-4 nonbonded calc exclusions.
<==
>
> 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?
>
==>
vdwmeth does a correction for pme. If pme is turned off, there is no pme.
So the long range vdw correction is not the source of your trouble.
<==
> 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.
>
==>
The pressure calc involves all kinds of jinxing around, especially for
systems with multi-atom molecules. I get the drift somewhere in here that
you are playing with a single atom molecules system. Good place to start.
Next start in runmd.f in sander and track backwards from the pressure calc,
seeing all that is involved. I think that having 1 atom/molecule will save
you grief on COM computation, but little else. This stuff is complicated
and it would require a lot of anyone's time to decipher it all.
<==
> 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.
>
==>
I really wonder if you are doing PBC and minimum image convention correctly.
<==
> 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).
>
==>
Probably more productive to understand the virial/pressure calcs looking at
the stuff I sent you from Bill Smith and then go through the code than to
look at the differences and guess...
<==
> 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)?
==>
When you have actual molecules rather than atoms, then there are all kinds
of complications associated with stuff that is excluded from the nonbonded
calc - these are basically 1-2, 1-3, and 1-4 interactions in amber ff's. I
really can't tell what is going on from what you describe though. I think
you are going to have to take a deep breath, start with the simplist system
possible, track through sander code to understand things, all that great
painful stuff that makes doing any of this consume time. Everybody seems to
think molecular dynamics is dead simple; heck, it's classical, and how hard
can handling a few bond, angles, dihedrals, and coulomb point charge and vdw
interactions be? Well, my favorite expression about all of this stuff is
"the devil is in the details"...
Good luck.
<==
>
> 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
>
-----------------------------------------------------------------------
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:33 PST
