Dear Ross,
Your answer is perfect fine and completely makes sense for me. I guess I
never saw this problem before is because a user spotted it here by using a
very "loose" EM protocol (vacuum and no constraints).
I agree that shake, lincs (gromacs) would solve this issue as I always did
for myself.
Thank you very much for you detailed explanation.
Best,
Alan
On Tue, Jun 15, 2010 at 13:24, Ross Walker <ross.rosswalker.co.uk> wrote:
> Hi 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.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<<
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
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<<
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Jun 15 2010 - 06:30:04 PDT
