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
Received on Thu Jun 28 2012 - 12:30:03 PDT
