[gmx-users] Re: question about [ pairs ] and interaction 1-4
Justin A. Lemkul
jalemkul at vt.edu
Wed Jun 16 02:41:32 CEST 2010
For what it's worth, I've had this experience with some molecules before with
phosphatidic acid derivatives. It seems to only happen during in vacuo EM.
When in the presence of solvent or bound to a protein it is less of an issue
since the surrounding tends to out-compete intramolecular interactions. The
biggest problem is in bilayers, where I've gotten around exploding molecules by
setting [exclusions] between the H and all phosphate O atoms. After solvation,
the problems tend to go away, since the solvent can satisfy the strong
electrostatic interactions in these parts of the molecule.
I don't know what your goals are, but perhaps this shared experience is useful.
-Justin
Alan wrote:
> Thanks Mark, here is the answer I got from Amber list and explains what
> I saw.
>
> Cheers,
>
> Alan
>
> ----------------------------
> This is an issue with the force field and has been noted several times
> before. I would take a look at the mailing list archive
> (http://archive.ambermd.org) - search for phosphate hydrogen and you
> will see lots of relevant answers. For example:
>
> http://archive.ambermd.org/200408/0178.html
>
> The issue is that there is no VDW radii on OH hydrogens. One possible
> fix for this is to add a small VDW radii to the H atom. I would also
> suggest using shake (even during the minimization) which will help since
> it will keep the O-H bond length fixed. Essentially none of the force
> fields were ever really designed to be used without shake so you can get
> weird things happen if you are not using shake. Phosphate groups are an
> extreme case however.
>
> Essentially, with water and hydroxyl hydrogens it is assumed the H atom
> is effectively within the VDW sphere of the oxygen atom, hence why it
> has zero VDW radii. Phosphates have a huge electrostatic interaction
> which can effectively override this. Note though that having shake stops
> the H bond increasing in length so it effectively remains within the
> oxygen VDW sphere. This may be enough to prevent the catastrophe you are
> seeing. Alternatively consider creating a new atom type for hydroxyl
> group bonded to oxygen bonded to phosphorus adding the small VDW radii
> for the H atom. Note gaff has a hp atom type:
>
> hp 0.6000 0.0157 same to hs (be careful !)
>
> Which is listed as H bonded to phosphate. Note this is different from
> hydroxyl:
>
> ho 0.0000 0.0000 OPLS Jorgensen, JACS,110,(1988),1657
>
> in that it has a VDW radii.
>
> Good luck,
>
> All the best
> Ross
>
>
>
> > -----Original Message-----
> > From: Alan [mailto:alanwilter at gmail.com <mailto:alanwilter at gmail.com>]
> > Sent: Tuesday, June 15, 2010 4:17 AM
> > To: AMBER Mailing List
> > Subject: [AMBER] to understand 1-4 interactions in Amber FF
> >
> > Hi there,
> >
> > I have this molecule ATP, deprotonated, with 46 atoms and -1 net charge
> > which topology was generated via Antechamber.
> >
> > I tested with sander, namd and gromacs (converting it via ACPYPE). The
> > latter is not discussed here although the general result is similar to
> > the
> > former 2.
> >
> > Essentially I am doing a energy minimisation (EM) following these
> > instructions (set only to emphasise the issue, not really for
> > production):
> >
> > #AMBER 11
> >
> > cat << EOF >| mdin
> > Minimization
> > &cntrl
> > imin=1, maxcyc=2000,
> > ntmin=2,
> > ntb=0,
> > igb=0,
> > cut=999,
> > /
> >
> > sander -O -i mdin -o mdout -p HATP_AC.prmtop -c HATP_AC.inpcrd
> >
> > ambpdb -p HATP_AC.prmtop < restrt > HATP_amber.pdb
> >
> > mdout OUTPUT:
> > NSTEP ENERGY RMS GMAX NAME
> > NUMBER
> > 1950 -2.2514E+08 5.9075E+13 3.4193E+14 O01
> > 4
> >
> > BOND = 35.0839 ANGLE = 124.4508 DIHED =
> > 28.9012
> > VDWAALS = -8.5795 EEL = 933.8418 HBOND =
> > 0.0000
> > 1-4 VDW = 15.4554 1-4 EEL = ************* RESTRAINT =
> > 0.0000
> >
> >
> > NSTEP ENERGY RMS GMAX NAME
> > NUMBER
> > 2000 -1.3381E+07 2.0866E+11 1.2077E+12 HHO
> > 1
> >
> > BOND = 35.0845 ANGLE = 124.4511 DIHED =
> > 28.9012
> > VDWAALS = -8.5795 EEL = 933.8418 HBOND =
> > 0.0000
> > 1-4 VDW = 15.4554 1-4 EEL = ************* RESTRAINT =
> > 0.0000
> >
> >
> > # NAMD 2.7b2
> > cat << EOF >| HATP_namd.conf
> > outputEnergies 50 # Energy output frequency
> > DCDfreq 2 # Trajectory file frequency
> > timestep 2 # in unit of fs
> > temperature 300 # Initial temp for velocity assignment
> > cutoff 999
> > switching off # Turn off the switching functions
> > PME off # Use PME for electrostatic calculation
> > amber on # Specify this is AMBER force field
> > parmfile HATP_AC.prmtop # Input PARM file
> > ambercoor HATP_AC.inpcrd # Input coordinate file
> > outputname HATP_namd # Prefix of output files
> > exclude scaled1-4
> > 1-4scaling 0.833333 # =1/1.2, default is 1.0
> > minimize 2000
> > EOF
> >
> > namd2 +p2 HATP_namd.conf >| namd.log
> >
> > LINE MINIMIZER BRACKET: DX 1.29731e-05 1.06277e-05 DU -1.70323e-05
> > 1.14306e-05 DUDX -2.62589 5.05233e-05 2.15101
> > LINE MINIMIZER REDUCING GRADIENT FROM 44.7384 TO 0.0447384
> > TIMING: 2000 CPU: 0.76958, 0.000377574/step Wall: 0.423322,
> > 0.00021419/step, 0 hours remaining, 98.023438 MB of memory in use.
> > ETITLE: TS BOND ANGLE DIHED
> > IMPRP
> > ELECT VDW BOUNDARY MISC
> > KINETIC TOTAL TEMP POTENTIAL
> > TOTAL3
> > TEMPAVG
> > ENERGY: 2000 31.2774 272.2172 33.5413
> > 0.0000
> > -2823.3313 23.6328 0.0000 0.0000
> > 0.0000 -2462.6627 0.0000 -2462.6627 -2462.6627
> > 0.0000
> >
> > vmd -parm7 HATP_AC.prmtop -rst7 HATP_AC.inpcrd HATP_namd.dcd
> >
> > I tried but I wish I could set the parameters in namd and sander so the
> > energy results could be closer. Anyway what's happening is that atoms
> > O01
> > and HHO in the phosphate have opposite charges and they are 1-4
> > neighbours
> > so after EM the Hydrogen collapse over the Oxygen (see the 1-4 EEL
> > energies
> > from sander mdout).
> >
> > With NAMD and VMD and I can reproduce and visualise this same issue
> > (happens
> > with Gromacs as well).
> >
> > I understand that in Amber 1-4 interactions are there although scaled
> > down
> > by a factor. What I just want to understand is if such a behaviour
> > would be
> > expected (although being more like an artefact) considering the
> > peculiarity
> > of the ATP molecule. And why vdW repulsion is not (apparently) acting
> > here?
> >
> > I attached HATP_amber.pdb, generated from the last frame of sander EM.
> >
> > Thanks in advance,
> >
> > Alan
> --
> Alan Wilter S. da Silva, D.Sc. - CCPN Research Associate
> Department of Biochemistry, University of Cambridge.
> 80 Tennis Court Road, Cambridge CB2 1GA, UK.
> >>http://www.bio.cam.ac.uk/~awd28<<
>
--
========================================
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
========================================
More information about the gromacs.org_gmx-users
mailing list