On Mon, 2014-04-14 at 22:59 +0800, 黄 wrote:
> Dear Amber Users:
> I know a single water molecule shouldn't have the nonbond
> interactions. I am surprised to see sander will calculates the
> electrostatic for a single TIP3P water molecule simply if the bond
> between hydrogen atoms is removed. What I did is as follows:
> --------- In Leap ----------------------------------------------------
> source leaprc.ff99SB
> w=copy WAT
> deletebond w.1.H1 w.1.H2
> saveamberparm w w.top w.crd
> ---------------------------------------------------------------------
>
>
> The sander output is as follows:
> ---------------------------------------------------------------------------------------------
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 1 3.8140E+01 1.1876E+01 1.9923E+01 H2 3
>
>
> BOND = 0.0000 ANGLE = 0.0000 DIHED = 0.0000
> VDWAALS = 0.0000 EEL = 38.1400 HBOND = 0.0000
> 1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT = 0.0000
> ----------------------------------------------------------------------------------------------------
> The electrostatic energy for single water molecule with H-H bond removed is 38kcal/mol !
>
>
> Am I doing something wrong?
Yes, you deleted the bond between the two hydrogen atoms. Don't do
that. Use a flexible water model instead of a rigid water model if you
don't want the bond between the water hydrogen atoms and want the angle
modulated by an angle parameter instead.
I'll give a detailed explanation for the observations here in the
postscript in case people are interested.
HTH,
Jason
The reason the energy is zero for the unaltered water is because 1-2 and
1-3 interactions are not computed (the interactions between those atom
pairs are accounted for _entirely_ by the valence bond and angle terms).
For rigid water molecules, the connectivity of the atoms are assigned
via 3 bond distances (O-H1, O-H2, and H1-H2). Another way of specifying
a water geometry is by specifying the O-H1 and O-H2 distances as well as
the H1-O-H2 angle. For a rigid water model, this does not matter. But
a key here is that the rigid water model does NOT contain an angle term
between the H1-O-H2 atoms of a water model.
Now into the guts of the energy calculation: pmemd/sander construct what
is called an "exclusion" list (effectively a linked list) that contains
a list of every atom whose nonbonded interactions should be ignored
because they are connected by 1 or 2 bonds (1-4 pairs are also included
in this list because they are added afterwards after being scaled down).
The tleap program actually generates the exclusion list (see
NUMBER_EXCLUDED_ATOMS and EXCLUDED_ATOMS_LIST in the topology file) and
does so by constructing the list of all atoms separated by 1, 2, or 3
bonds. A quirky implementation detail is that the GB electrostatic code
uses the exclusion list defined in the topology file. The PME
electrostatic code ignores the prmtop exclusion list and constructs its
own. That is why Brian's suggestion of using igb=6 (vacuum
electrostatics via the GB code path) gives an energy of 0, which you
would expect. My guess is you were using ewmeth=4 in the &ewald
namelist to do vacuum electrostatics that way (you need to provide your
input file for us to help).
The PME code, however, apparently constructs the exclusion list by
looking at all atom pairs in the bond, angle, and dihedral lists (rather
than using just the bonds like tleap appears to do). As a result, since
there is no H1-O-H2 angle term defined in the topology file, H1 and H2
are not excluded from one another. The resulting energy (charge of
0.417 separated by a distance of 1.514 Angstroms) is:
0.417 * 18.2223 * 0.417 * 18.2223 / 1.514 = 38.14
The 18.2223 factor is needed to convert to the Amber-consistent set of
charges (specifically, it is the square root of the electrostatic
constant giving the final units of kcal/mol)
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon Apr 14 2014 - 09:00:02 PDT
