[gmx-users] Re: question about [ pairs ] and interaction 1-4
Alan
alanwilter at gmail.com
Tue Jun 15 15:18:52 CEST 2010
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]
> 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<<
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100615/40a1bffe/attachment.html>
More information about the gromacs.org_gmx-users
mailing list