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
Received on Tue Jun 15 2010 - 05:30:04 PDT
