[AMBER-Developers] Possible bug with ntt=3 and ntr=1

From: Ross Walker <ross.rosswalker.co.uk>
Date: Tue, 7 Jul 2009 07:15:00 +0100

Hi All,

I am trying to track down what I think is a bug with GB runs where ntt=3 and
harmonic restraints (ntr=1) are in use. I am currently seeing weird behavior
with simulations of the nucleosome in GB where I am restraining just the
DNA. With gamma_ln >= 1.0 the system is stable but the dynamics is plain
wrong. As I reduce gamma_ln the system becomes more and more unstable. This
is strange since it would seem to me that as gamma_ln tends towards zero is
should approximate NVE which in this case is completely stable.

My input file looks like this:

&cntrl
  imin=0,irest=1,ntx=5,
  nstlim=250000,dt=0.002,
  ntf=2,ntc=2,ntb=0,
  igb=5,cut=10.0,rgbmax=10.0,
  ntpr=100,ntwx=100,ntwr=10000,
  ntr=1,saltcon=0.1,
  tempi=310.0,temp0=310.0,
  ntt=3,gamma_ln=0.001,
  ig=1037891,
 /
RESTRAIN DNA
0.5
RES 1 294
END
END

In this case with gamma_ln=0.001 the system crashes after 3100 steps with a
temperature of 4800K!. This is true both of sander and pmemd.

Interestingly if I turn off restraints all is fine. Similarly if I set ntt=0
(NVE) all is fine. Example outputs are attached.

       mdout.orig = pmemd with ntt=3, ntr=1
mdout.orig.sander = sander with ntt=3, ntr=1
  mdout.orig.ntr0 = pmemd with ntt=3, ntr=0
   mdout.orig.nve = pmemd with ntt=0, ntr1

Interestingly everything looks good until NSCM is triggered at step 1000.

NSTEP = 900 TIME(PS) = 101.800 TEMP(K) = 308.05 PRESS =
0.0
 Etot = -66743.9688 EKtot = 19528.7840 EPtot =
-86272.7528
 BOND = 5722.9654 ANGLE = 13508.0515 DIHED =
17017.8931
 1-4 NB = 5619.2367 1-4 EEL = 1382.5208 VDWAALS =
-13061.4667
 EELEC = -73648.3380 EGB = -43368.0253 RESTRAINT =
554.4096
 EAMBER (non-restraint) = -86827.1624
 
---------------------------------------------------------------------------- 
--
| RE_POSITION Moving by  -0.484242 -0.548000 -0.613646
 NSTEP =     1000   TIME(PS) =     102.000  TEMP(K) =   306.89  PRESS =
0.0
 Etot   =    -66749.0252  EKtot   =     19454.8086  EPtot      =
-86203.8338
 BOND   =      5765.1339  ANGLE   =     13518.0488  DIHED      =
17010.2508
 1-4 NB =      5557.7458  1-4 EEL =      1277.5502  VDWAALS    =
-13007.4986
 EELEC  =    -74808.4224  EGB     =    -42038.1373  RESTRAINT  =
521.4950
 EAMBER (non-restraint)  =    -86725.3288
 
----------------------------------------------------------------------------
--
 NSTEP =     1100   TIME(PS) =     102.200  TEMP(K) =   365.32  PRESS =
0.0
 Etot   =    -62356.5394  EKtot   =     23159.5258  EPtot      =
-85516.0652
 BOND   =      5762.5203  ANGLE   =     13597.3821  DIHED      =
17100.7761
 1-4 NB =      5590.5033  1-4 EEL =      1324.8176  VDWAALS    =
-12951.1683
 EELEC  =    -77138.7949  EGB     =    -39800.5919  RESTRAINT  =
998.4905
 EAMBER (non-restraint)  =    -86514.5557
 
----------------------------------------------------------------------------
--
The restraint energy then takes a jump and the system never recovers and the
temperature just keeps rising.
Thus I think the problem is related to NSCM recentering the system and not
allowing for restraints properly.
Running the simulation with nscm explicitly set to 0 also appears to fix the
problem.
mdout.orig.nscm0.bz2
Inspection of the re-centering routine (in both PMEMD and sander) shows that
there is a call to re_position:
which should execute the following code:
  if (master) &
    write(mdout, '(a, 3f10.6)')"| RE_POSITION Moving by ", xd, yd, zd
  vec(1) = vec(1) + xd
  vec(2) = vec(2) + yd
  vec(3) = vec(3) + zd
  do i = 1, n
    x(1,i) = xd + x(1,i)
    x(2,i) = yd + x(2,i)
    x(3,i) = zd + x(3,i)
  end do
  if (ntr .gt. 0) then
    do i = 1, n
      xr(1,i) = xd + xr(1,i)
      xr(2,i) = yd + xr(2,i)
      xr(3,i) = zd + xr(3,i)
    end do
  end if
However, in both sander and PMEMD ntr has been hard coded to be zero inside
this routine:
          call re_position(atm_cnt, 0, crd, crd, vcm(1), vcm(2), vcm(3), &
                           sys_x, sys_y, sys_z, sys_range, 0)
subroutine re_position(n, ntr, x, xr, xc, yc, zc, x0, y0, z0, vec, verbose)
Thus I tried changing this to execute the above code when ntr=1 and ntt=3
which it would seem should fix the problem but alas it does not:
mdout.mod.bz2
Thus it would seem that ntr=1 with ntt=3 and nscm = 1000 (default) is
broken. However, the fact that the re_positioning code exists for ntr=1 and
ntt=3 but has been hard wired to not be used means that someone at least was
aware of the situation.
Comments / suggestions for fixes anyone?
All the best
Ross
/\
\/
|\oss Walker
| Assistant Research Professor |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross.rosswalker.co.uk | 
| http://www.rosswalker.co.uk | PGP Key available on request |
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.
Received on Tue Jul 07 2009 - 01:07:23 PDT
Custom Search