Re: [AMBER] Question regarding HAS_10_12 flag in PMEMD (LINMIN Failure)

From: Ross Walker <ross.rosswalker.co.uk>
Date: Tue, 11 Oct 2011 19:19:21 -0700

Hi Brian,

I don't the HAS_10_12 option has been tested in many many moons since none
of the current force fields use it and since it requires recompilation with
an ifdef it was never built and tested as part of the standard regression
tests.

> > sander:
> >
> > NSTEP ENERGY RMS GMAX NAME
> NUMBER
> > 1 -1.8508E+04 3.9593E+02 4.3374E+04 Cl-
> 1235
> >
> > BOND = 126.0429 ANGLE = 276.9695 DIHED =
> > 704.5199
> > VDWAALS = 22134.3046 EEL = -45257.0072 HBOND =
> > -19.4807
> > 1-4 VDW = 389.3294 1-4 EEL = 3137.1801 RESTRAINT =
> > 0.0000
> >
> >
> > pmemd:
> >
> > NSTEP ENERGY RMS GMAX NAME
> NUMBER
> > 1 -1.8365E+04 3.9887E+02 4.3374E+04 Cl-
> 1235
> >
> > BOND = 126.0429 ANGLE = 276.9695 DIHED =
> > 704.5199
> > VDWAALS = 22134.3046 EEL = -45257.0072 HBOND =
> > 123.1766
> > 1-4 VDW = 389.3294 1-4 EEL = 3137.1801 RESTRAINT =
> > 0.0000

Note the ONLY difference here in these two outputs is in the HBOND energy.
Everything else matches perfectly so this should not be too difficult to
track down. I would start first by confirming that the problem is still
there in serial. Can you run sander in serial and pmemd in serial for
maxcyc=1 to just get the initial energy and see what you get. This will
determine if it is just a parallel communication issue or is a bug in the
underlying mathematical logic. Of course this says nothing about which is
actually right here, pmemd or sander or indeed they could both be wrong.
However, the limin failures, which suggest in this case that the forces
don't match the energy implies that PMEMD is incorrect. But of course sander
could still be dodgy as well. I would see if you have an earlier sander
version, version 8 or 9 will be good and run a single step in that and see
if it matches sander 11. If yes then at least we will know where to look in
PMEMD. If you can post your prmtop and inpcrd files it will certainly help
in debugging things. - Send them to me directly if they are too big for the
mailing list, or you'd rather not share them publically.

Note this is essentially the code involved: (given you have -DDIRFRC_NOVEC
defined and are, I assume, running with ntb=1).

PMEMD
#ifdef HAS_10_12
      if (ic .gt. 0) then
#endif
        r6 = delr2inv * delr2inv * delr2inv
        f6 = cn2(ic) * r6
        f12 = cn1(ic) * (r6 * r6)
        df = df + (12.d0 * f12 - 6.d0 * f6) * delr2inv
#ifdef HAS_10_12
      else
        ! This code allows 10-12 terms; in many (most?) (all?) cases, the
        ! only "nominal" 10-12 terms are on waters, where the asol and bsol
        ! parameters are always zero; hence we can skip this part.
        ic = - ic
        delr2inv = delrinv * delrinv
        r10 = delr2inv * delr2inv * delr2inv * delr2inv * delr2inv
        f10 = gbl_bsol(ic) * r10
        f12 = gbl_asol(ic) * (r10 * delr2inv)
        df = df + (12.d0 * f12 - 10.d0 * f10) * delr2inv
      end if
#endif /* HAS_10_12 */

SANDER
#ifdef HAS_10_12
         cache_r2(im_new)=delr2inv
#endif

#ifdef HAS_10_12
   delr2inv = cache_r2(im_new)
   ic = -ico(iaci+iac(j))
   r10 = delr2inv*delr2inv*delr2inv*delr2inv*delr2inv
   f10 = bsol(ic)*r10
   f12 = asol(ic)*(r10*delr2inv)
   ehb = ehb + f12 - f10
   ! -- ti decomp: add hbond-terms to vdw-energy in energy decomposition
   if(decpr .and. idecomp > 0) call decpair(3,i,j,(f12 - f10)/(nstlim/ntpr))
   df = df + (12.d0*f12 - 10.d0*f10)*delr2inv
#endif

#ifdef HAS_10_12

         ! --- this code allows 10-12 terms; in many (most?) (all?) cases,
the
         ! only "nominal" 10-12 terms are on waters, where the asol and
bsol
         ! parameters are always zero; hence we can skip the L-J part;
note
         ! that we still have to compute the electrostatic interactions

         ic = -ico(iaci+iac(j))
         r10 = delr2inv*delr2inv*delr2inv*delr2inv*delr2inv
         f10 = bsol(ic)*r10
         f12 = asol(ic)*(r10*delr2inv)
         ehb = ehb + f12 - f10
         df = (twelve*f12 - ten*f10)*delr2inv
#else
         df = zero
#endif

Not much to go on really so my initial guess would be this is something to
do with the parameters not being correctly read / stored or some issue in
parallel that works fine in serial.

Beyond that it is a case of going into the code and checking the contents of
the relevant values at each stip to see where the difference comes from.

All the best
Ross

/\
\/
|\oss Walker

---------------------------------------------------------
| Assistant Research Professor |
| San Diego Supercomputer Center |
| Adjunct Assistant Professor |
| Dept. of Chemistry and Biochemistry |
| University of California San Diego |
| NVIDIA Fellow |
| http://www.rosswalker.co.uk | http://www.wmd-lab.org/ |
| Tel: +1 858 822 0854 | EMail:- ross.rosswalker.co.uk |
---------------------------------------------------------

Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.





_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Oct 11 2011 - 19:30:02 PDT
Custom Search