I thought nscm >0 did not work with ntr. if you change cm it doesn't
do the same to refc. can't check the code now but that seems likely.
On 7/7/09, Ross Walker <ross.rosswalker.co.uk> wrote:
> 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.
>
>
>
--
===================================================================
Carlos L. Simmerling, Ph.D.
Professor, Department of Chemistry
CMM Bldg, Room G80 Phone: (631) 632-1336 Fax: 632-1555
Stony Brook University E-mail: carlos.simmerling.gmail.com
Stony Brook, NY 11794-5115 Web: http://comp.chem.sunysb.edu
===================================================================
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Tue Jul 07 2009 - 10:07:53 PDT
