Re: [AMBER] NaN with GTX580 pmemd.cuda simulation with and without iwrap = 1

From: Ross Walker <ross.rosswalker.co.uk>
Date: Thu, 28 Jun 2012 14:05:55 -0700

Hi Fabricio,

I suggest 'NOT' using parmed in this situation - if you start hacking the
prmtop file you will most likely forget what you did, or not interpret it
properly and when it comes time to publish and someone asks for the force
field parameters file and setup files to reproduce your work you won't have
them. Instead I would go back and modify the original parameterization files
you used in leap. This way it is clearly documented what you did and what
you changed. The hacking approach is useful to quickly test a hypothesis but
is no excuse for not being diligent in the approach in the first place.

The VDW Radii do not come from the pdb file. They are VDW parameters
specified for the atom type in the parmXX.dat file or frcmod file. Here's an
example:

http://ambermd.org/tutorials/advanced/tutorial1_adv/

plc.frcmod:

# modifications to force field for poplar plastocyanin

MASS
SM 32.06
CU 65.36

BOND
NB-CU 70.000 2.05000 #kludge by JRS
...

ANGLE
CU-NB-CV 50.000 126.700 #JRS estimate
...

DIHE
X -NB-CU-X 1 0.000 180.000 3.000
...

NONBON
CU 2.20 0.200
SM 2.00 0.200

The values in the nonbond section correspond to the following:

                    LTYNB , R , EDEP

         LTYNB Atom symbol.

         R The van der Waals radius of the atoms having the symbol
                     "LTYNB" (Angstoms)

         EDEP The 6-12 potential well depth. (kcal/mol)

If you then look at $AMBERHOME/dat/leap/parm/parm99.dat you will see:

  H 0.6000 0.0157 !Ferguson base pair geom.
  HO 0.0000 0.0000 OPLS Jorgensen,
JACS,110,(1988),1657
  HS 0.6000 0.0157 W. Cornell CH3SH --> CH3OH FEP

And there is your problem. The HO has no VDW interaction and thus can
collapse on other atoms. It doesn't if it is bonded or angled with them
since nonbond interactions are excluded for 1-2 and 1-3 terms but if it is a
1-4 term or a 1-5 term or more than it has a charge interaction which could
be strong if they are different signs and pull the HO atom into it.

One approach is what they do for hydroxyl's in Charmm All22:

!V(Lennard-Jones) = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
...
HT 0.000000 -0.046000 0.224500 ! ALLOW WAT

Or you can look at what they do for GLYCAM_06H which is to define a Ho atom
type and give it a small VDW term:

  Ho 0.2000 0.0300 M.B. Tessier 2011, use by
permission only

This is how they stop the hydroxyl hydrogens collapsing onto hydroxyl
oxygens in sugar rings.

These are both hacks I believe, I don't think they actually fitted the
parameters here but you might want to ask the authors of those force fields
for details behind it. Then for your system you'll take your frcmod file
that you build during the parameterization and also the mol2 file (or lib
file) edit the atom types for the hydroxyl's to be something unused and then
provide the parameters (mass, bond, angle, dihedral) for this in the frcmod
file and finally at the end include the small VDW parameters that you want
for that atom.

Good luck,

All the best
Ross


> -----Original Message-----
> From: Fabrício Bracht [mailto:bracht.iq.ufrj.br]
> Sent: Thursday, June 28, 2012 1:29 PM
> To: AMBER Mailing List
> Subject: Re: [AMBER] NaN with GTX580 pmemd.cuda simulation with and
> without iwrap = 1
>
> Hello Ross. Thank you for the reply. I´ll restart the parameterization
> procedure. You suggest I create a new H type. I understand that
> setting a new type for this hydrogen would require me to include
> vdwradii for this atom. The problem is that I am not sure how I would
> do this. I would just have to rename the pdb source file? But then,
> how would I set the vdwradii? Can´t I just set the vdwradii for this
> particular atom in tleap, after having performed the parameterization?
> Wouldn´t the lmodprmtop program be of use in this case?.
> Thank you in advance
> Fabrício Bracht
>
> 2012/6/28 Ross Walker <ross.rosswalker.co.uk>:
> > Hi Fabricio,
> >
> > So you are running into one of the problems with the AMBER force
> field that
> > people have ignored for years and it was never too much of a problem
> since
> > people didn't typically run for long enough to see it. That is that
> hydroxyl
> > hydrogens have a zero VDW radii. In TIP3P water this is not a problem
> > because it is rigid and the H sits inside the VDW radii of the
> oxygen.
> > However, when you start moving to highly charged systems, a classic
> example
> > being phosphates you can end up with the H collapsing into one of the
> other
> > negatively charged atoms for which there is not a bond or angle term
> to
> > prevent it and as soon as that distance becomes very small so you
> divide by
> > zero and get NAN. The CPU is a little more tolerant of division by
> VERY
> > small numbers than the GPU and that's why it manifests itself more as
> a blow
> > up (trapped by vlimit) on the CPU and you get NAN on the GPU. Either
> way the
> > simulation is bust because this event is non physical. 1fs can help,
> since
> > the thing won't shoot off as much but it still won't really fix the
> problem.
> >
> > I would take a VERY careful look at the parameterization - that is
> the
> > underlying problem. To begin with you probably should not use type HO
> for
> > the hydrogen and instead create a new type. Then make sure this new H
> type
> > has at least a small VDW radii, a very small radii will help (look in
> the
> > frcmod file) but ideally one should probably attempt to refit the
> radii
> > properly for such cases but this is tricky. In your case I would just
> put a
> > small radii on the hydrogen (take a look at what the charmm people do
> with
> > TIP3P - which in charmm should really strictly be called TIP3P-CHARMM
> since
> > it isn't the original TIP3P parameterization due to the addition of
> small
> > radii on the hydroxyl hydrogens.)
> >
> > What this does highlight though is that it is important that people
> still
> > test things on the CPU and not just rely on the GPU code for
> everything. If
> > you see a problem with a GPU run then switch to CPU and analyse
> things
> > carefully. Don't expect the CPU code to also explode since it is more
> > tolerant of marginal stability cases but at least be aware that you
> might be
> > simulating something that is borderline stable at best.
> >
> > All the best
> > Ross
> >
> >> -----Original Message-----
> >> From: Fabrício Bracht [mailto:bracht.iq.ufrj.br]
> >> Sent: Thursday, June 28, 2012 12:09 PM
> >> To: AMBER Mailing List
> >> Subject: Re: [AMBER] NaN with GTX580 pmemd.cuda simulation with and
> >> without iwrap = 1
> >>
> >> Hello. Thank you for the reply. In fact, a time step of 1fs seems to
> >> work around the problem. At the same time as i was writing the last
> >> email, I had submitted a job with dt = 0.001 and it finished without
> >> any problems. I haven´t tested yet with pmemd.cuda, but I´ll do that
> >> later on.
> >> The two atoms involved in this case are O and H of the hydroxyl
> >> residue linked to the zinc atom of the active site. Where do I find
> >> the Lennard-Jones rmin values for them? The hole active site was
> >> parameterized following the steps on the MTK++ manual.
> >> I´ll follow your instructions as to getting a restart file from a
> >> previous step and test it out. I´ll inform you as soon as I have the
> >> results.
> >> Thank you
> >> Fabrício Bracht
> >>
> >> 2012/6/28 David A Case <case.biomaps.rutgers.edu>:
> >> > On Thu, Jun 28, 2012, Fabrício Bracht wrote:
> >> >
> >> >> Hi. I did what you asked. I ran the simulation again, in order to
> >> make
> >> >> sure that the error happened again. This time, it dies a lot
> earlier
> >> >> with the same warnings as last time (vlimit exceed....etc). I
> >> >> continued the simulation using as input the last restrt file
> written
> >> >> using the serial version of pmemd. The simulation crashed after a
> >> few
> >> >> hours , the last output lines are:
> >> >>
> >> >>
> >> >>
> >> >> vlimit exceeded for step 112126; vmax =    91.1788
> >> >> vlimit exceeded for step 112176; vmax =    45.6806
> >> >>
> >> >>      Coordinate resetting cannot be accomplished,
> >> >>      deviation is too large
> >> >>      iter_cnt, my_bond_idx, i and j are
> >> :       2    2690    5315    5316
> >> >>
> >> >
> >> > Thanks for the info.  It looks to me like the problem has nothing
> to
> >> do (per
> >> > se) with GPU's, since you can reproduce it on a serial CPU using
> >> pmemd.
> >> > Unfortuantely, problems like this (generally with the force field)
> >> can be hard
> >> > to track down.  Your mdin file looks fine to me, although it is
> safer
> >> to use
> >> > dt=0.001 rather than 0.002, especially for constant pressure
> >> simulations, and
> >> > if you system is not *quite* equilibrated.
> >> >
> >> > The information above indicates that the problem occurs around
> atoms
> >> 5315
> >> > and 5316.  You might see what those are, as this might (or might
> not)
> >> give
> >> > a clue to the origin of the problem.  Look hard (and/or let us
> know
> >> about)
> >> > any non-standard residues or molecules that are in your
> >> simulation.  Do
> >> > you have any Lennard-Jones rmin values that are either zero or
> small
> >> (i.e.
> >> > less than typical values for that type of atom)?  You may have to
> get
> >> a
> >> > restart file someplace a little before step 112126; (don't use
> ig=-1
> >> to do
> >> > this: copy the actual seed that is printed in your output file, so
> >> you will
> >> > get an identical trajectory.)  Then you can set ntwx=1, and
> examine
> >> the
> >> > coordinates at every step as the system crashes....one hopes that
> >> this will
> >> > identify the source of the problem.
> >> >
> >> > ...good luck...dac
> >> >
> >> >
> >> > _______________________________________________
> >> > 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
> >
> >
> > _______________________________________________
> > 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


_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jun 28 2012 - 14:30:02 PDT
Custom Search