Thanks for your helps! The explanation is very clear.
My amber version is Amber 11. The input file I used was:
-------- input for WAT energy calculation--------------
&cntrl
imin=1,
maxcyc=0,
ncyc=0,
ntb=0,
cut=9999.0
/
-------------------------------------------------------
By playing around the water molecule to further understand the "excluded_nonbond_list", I end up to have one more question. What I did is:
1). generate prep and frcmod files for a single water molecule with H-H disconnected. The residue name is set as "WT1" to differ from TP3 model.
2). atom types for oxygen and hydrogens are OW and HW, respectively. Then, sander calculates nonbond electrostatics with the above sander controlling input file
3). atom types for oxygen and hydrogens are OW and HO, respectively, Then, sander gives zero for nonbond electrostatics.
It looks like whenever there is HW atom, the nonbond electrostatics are calculated in PME code. Does it mean the PME code ignores the angle term for HW atom type?
Attached are my prep and frcmod files in my test.
Best,
Jinfeng
ÔÚ 2014-04-14 23:37:58£¬"Jason Swails" <jason.swails.gmail.com> дµÀ£º
>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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
- application/octet-stream attachment: WT1.prep
Received on Tue Apr 15 2014 - 21:30:02 PDT
