Re: AMBER: hydrogens are flying in replica exchange simulations

From: Ashish Sangwai <ashishsangwai.gmail.com>
Date: Mon, 21 Apr 2008 19:37:51 -0500

Hello,

I agree with Dr. Simmerling that this is something deeper. Shake did not
fail in the first appearance when bond lengths were out of control. It
failed after a while.

With my experience of AMBER, shake failures are reported in normal MD on the
first failed MD step. Which is not a case in REMD, it went on for a while
until hydrogens were too far from each other.

I am using AMBER 9 compiled on AMD opteron machines with ifort from
intel/9.1 and mvapich2-intel-tcp

I am sending my files input (only the higher temperature replicas) as
attachment to this mail.

1a.prmtop = topology
restrt0 = output of minimization
inp.eq = sample input used for equilibration
inp.440 = sample REM input
groufile = indicates the temperatures this problem occurs at

Thank you very much.

- Ashish Sangwai




On Mon, Apr 21, 2008 at 10:10 AM, Carlos Simmerling <
carlos.simmerling.gmail.com> wrote:

> I agree with Adrian except that if the bond lengths get very long at high
> T due to
> not having bond forces applied, then shake should fail and it appears
> that it
> didn't in this case.
>
> we often use ntf=1 and ntc=2 with the rationale that the forces help to
> make shake convergence easier. with ntc=ntf=2, this should still work
> fine and it puzzles me that it works in MD and not REMD. I have a hard
> time imagining why this would be so since REMD and shake are pretty much
> orthogonal.
>
> maybe you could send me your inputs and I could test on our machines.
> make sure to let me know exactly which amber version you used.
> carlos
>
>
> On Mon, Apr 21, 2008 at 10:56 AM, Adrian Roitberg <roitberg.qtp.ufl.edu>
> wrote:
>
> > Dear all,
> >
> > I want to make a couple of comments just to clarify the behavior Ashish
> > was observing.
> >
> > The use of shake for bonds with H atoms (ntc=2) basically works like
> > this: the system takes one regular md step using the computed forces
> > and then 'corrects' the position of the H atoms by using Shake, which as
> > implemented (except for water) is an iterative solver.
> >
> > Since the positions of the H atoms is then more or less independent of
> > the forces acting on them, one could imagine NOT computing those forces and
> > then letting SHAKE do the work. If the H atoms are not going too far, this
> > saves time on computing forces. However, if the new positions of the H atoms
> > is such that SHAKE cannot correct for them, or if it takes too many
> > iterations to correct, then the MD fails.
> >
> > This is more likely to happen in high density fluids, high temperatures,
> > etc.
> >
> > So, Ashish was running with ntf=2, which means, according to the manual:
> >
> > "ntf
> > Force evaluation. Note: If SHAKE is used (see NTC), it is not necessary
> > to calculate
> > forces for the constrained bonds.
> > = 1 complete interaction is calculated (default)
> > = 2 bond interactions involving H-atoms omitted (use with NTC=2)"
> >
> > that the forces on the H atoms are not computed. AT regular temperatures
> > this would not be an issue. However, at the high temperatures REMD is run,
> > this might be a problem.
> >
> > So, changing ntf to 1 made the H atoms behave 'better' by computing the
> > forces on them, and then SHAKE was able to work.
> >
> > Just my 3 cents' worth of opinion.
> >
> > a.
> >
> >
> >
> >
> >
> >
> > Ashish Sangwai wrote:
> >
> > > Hello,
> > >
> > > I checked the $AMBERHOME/test/rem_gb_4rep/ in parallel. It works
> > > fine.
> > >
> > > I also checked that individual MD runs (rem=0) are stable even at
> > > temperature of 600K.
> > >
> > > Then, I observed that 'ntf' in your test run files are 1 (total
> > > potential
> > > calculation) even when ntc=2 (shake on).
> > >
> > > I tried the same in my mdin files for REM and it worked. Now the
> > > hydrogens
> > > are not dissociating.
> > >
> > > My new mdin file looks like
> > >
> > > &cntrl
> > > imin=0, ntx=5, irest=1, ntxo=1, ntpr=500, ntwr=500, ntwx=200,
> > > ntwe=500, ntf=1, ntb=0, igb=5, nstlim=10000,
> > > temp0 = 300, ntt=2, tautp=1.0,
> > > ntp=0, cut=999.0, saltcon=0.001,
> > > taup=2.0, ntc=2, ntrx=1,
> > > dt=0.001, ig=323657, numexchg=500,
> > > &end
> > > END
> > >
> > >
> > > The problem is solved for now but I am not sure whether I should trust
> > > the
> > > REM simulation.
> > >
> > > Thank you very much for your help.
> > >
> > > - Ashish Sangwai
> > > 341 Lindy Boggs Center
> > > Chemical Engineering
> > > Tulane University
> > > New Orleans LA 70118
> > >
> > >
> > >
> > >
> > > On Sat, Apr 19, 2008 at 8:49 AM, Carlos Simmerling <
> > > carlos.simmerling.gmail.com> wrote:
> > >
> > > try running a normal MD simulation at 400K or higher.
> > > > I have never seen something like this. it also seems strange
> > > > that your bonds are so long but the bond energy is not high-
> > > > is that energy output you gave for one of the structures with
> > > > dissociated hydrogen?
> > > > also check the test cases like Ross suggested.
> > > >
> > > >
> > > > On Fri, Apr 18, 2008 at 8:18 PM, Ashish Sangwai <
> > > > ashishsangwai.gmail.com>
> > > > wrote:
> > > >
> > > > This does not occur for temperatures below 380K. After that for
> > > > > 400K,
> > > > > 420K and 440K this was seem to be happening.
> > > > >
> > > > > I tried running with rem = 0 and that run is just fine. It is
> > > > > happening
> > > > > only when rem=1.
> > > > >
> > > > > The averages in out files were in following way...
> > > > >
> > > > > NSTEP = 10000 TIME(PS) = 60.000 TEMP(K) = 495.35
> > > > > PRESS
> > > > > = 0.0
> > > > > Etot = 56.1676 EKtot = 58.5688 EPtot =
> > > > > -2.4012
> > > > > BOND = 12.3056 ANGLE = 27.2254 DIHED =
> > > > > 33.3567
> > > > > 1-4 NB = 7.6295 1-4 EEL = 310.9982 VDWAALS =
> > > > > -4.9558
> > > > > EELEC = -158.7368 EGB = -230.2240 RESTRAINT =
> > > > > 0.0000
> > > > >
> > > > >
> > > > > ------------------------------------------------------------------------------
> > > > >
> > > > > Also, the simulation fails after certain point when exchange
> > > > > occurs with
> > > > > such an unbound state with replica at 300K and AMBER stops with
> > > > > message
> > > > >
> > > > > Coordinate resetting (SHAKE) cannot be accomplished,
> > > > >
> > > > > in between atom number 35 and 37 (alanine carbon and hydrogen)
> > > > >
> > > > > Also, I am doing these termini because I want to simulate a zero
> > > > > charge
> > > > > state. And this could be a small test simulation to check how that
> > > > > goes.
> > > > >
> > > > > Thank you very much,
> > > > >
> > > > > Ashish Sangwai
> > > > >
> > > > > p.s. - nice talk at ACS New Orleans in replica exchange symposium
> > > > >
> > > > > On Fri, Apr 18, 2008 at 7:02 PM, Carlos Simmerling <
> > > > > carlos.simmerling.gmail.com> wrote:
> > > > >
> > > > > that is very surprising- you might want to run the same inputs but
> > > > > > set
> > > > > > rem=0
> > > > > > and remove numexchg (just normal MD) and see if you have the
> > > > > > same
> > > > > > problem. I can't imagine why REMD would allow the hydrogens to
> > > > > > move so
> > > > > > far.
> > > > > > what are your energies like in the out files?
> > > > > >
> > > > > > as a note, your termini are unusual... are you sure that's what
> > > > > > you
> > > > > > want?
> > > > > >
> > > > > >
> > > > > > On Fri, Apr 18, 2008 at 7:43 PM, Ashish Sangwai <
> > > > > > ashishsangwai.gmail.com> wrote:
> > > > > >
> > > > > > Hello,
> > > > > > >
> > > > > > > I am trying to carry out replica exchange simulation on
> > > > > > > polyalanine
> > > > > > > (5 mer) system in Generalized Born solvent.
> > > > > > >
> > > > > > > For replicas above 440 K, the hydrogens in the system are
> > > > > > > having
> > > > > > > unbound co-ordinates.
> > > > > > >
> > > > > > > My input file looks like :
> > > > > > >
> > > > > > > &cntrl
> > > > > > > imin=0, ntx=5, irest=1, ntxo=1, ntpr=500, ntwr=500,
> > > > > > > ntwx=200,
> > > > > > > ntwe=500, ntf=2, ntb=0, igb=5, nstlim=10000,
> > > > > > > temp0 = 300, tempi=300, ntt=2, tautp=1.0,
> > > > > > > ntp=0, cut=999.0, saltcon=0.001,
> > > > > > > taup=2.0, ntc=2, ntrx=1,
> > > > > > > dt=0.001, ig=323657, numexchg=500,
> > > > > > > &end
> > > > > > > END
> > > > > > >
> > > > > > > System was prepared in xleap with following commands :
> > > > > > >
> > > > > > > source leaprc.ff03
> > > > > > > set default PBradii mbondi2
> > > > > > > 1a = sequence {ALA ALA ALA ALA ALA}
> > > > > > > saveamberparm 1a 1a.prmtop 1a.inpcrd
> > > > > > >
> > > > > > > In the most recent PDB file of REM simulation, hydrogen
> > > > > > > co-ordinates are unbound.....
> > > > > > >
> > > > > > > This PDB was generated after 38 exchange attempts for a
> > > > > > > replica at
> > > > > > > 440 K. Rest of the backbone is still bound .
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > >
> > > > > > > REMARK ALA
> > > > > > > ATOM 1 N ALA 1 4.380 1.257 0.972
> > > > > > > ATOM 2 H ALA 1 117.737 228.487 148.973
> > > > > > > ATOM 3 CA ALA 1 5.181 0.081 0.806
> > > > > > > ATOM 4 HA ALA 1 -8.373 -76.768 181.924
> > > > > > > ATOM 5 CB ALA 1 4.404 -0.715 -0.190
> > > > > > > ATOM 6 1HB ALA 1 -94.721 143.367-170.150
> > > > > > > ATOM 7 2HB ALA 1 258.561-362.545-208.892
> > > > > > > ATOM 8 3HB ALA 1 -121.572 -49.341 56.075
> > > > > > > ATOM 9 C ALA 1 6.711 0.289 0.471
> > > > > > > ATOM 10 O ALA 1 7.423 -0.676 0.212
> > > > > > > ATOM 11 N ALA 2 7.223 1.533 0.473
> > > > > > > ATOM 12 H ALA 2 -125.699 177.403 48.942
> > > > > > > ATOM 13 CA ALA 2 8.634 1.801 0.462
> > > > > > > ATOM 14 HA ALA 2 169.270-226.264-215.241
> > > > > > > ATOM 15 CB ALA 2 9.010 3.115 -0.180
> > > > > > > ATOM 16 1HB ALA 2 -311.134 87.755-407.032
> > > > > > > ATOM 17 2HB ALA 2 -64.133 420.606 351.098
> > > > > > > ATOM 18 3HB ALA 2 355.577 -48.094-155.287
> > > > > > > ATOM 19 C ALA 2 9.323 1.628 1.800
> > > > > > > ATOM 20 O ALA 2 8.714 1.354 2.846
> > > > > > > ATOM 21 N ALA 3 10.639 1.796 1.792
> > > > > > > ATOM 22 H ALA 3 55.123 15.439-100.645
> > > > > > > ATOM 23 CA ALA 3 11.555 2.060 2.927
> > > > > > > ATOM 24 HA ALA 3 -216.220 363.839 259.284
> > > > > > > ATOM 25 CB ALA 3 11.623 0.734 3.689
> > > > > > > ATOM 26 1HB ALA 3 155.952-122.441 -95.372
> > > > > > > ATOM 27 2HB ALA 3 79.353 44.635 188.102
> > > > > > > ATOM 28 3HB ALA 3 -407.256-128.682 -6.444
> > > > > > > ATOM 29 C ALA 3 12.963 2.618 2.563
> > > > > > > ATOM 30 O ALA 3 13.929 2.618 3.318
> > > > > > > ATOM 31 N ALA 4 13.114 3.032 1.325
> > > > > > > ATOM 32 H ALA 4 12.399 2.677 0.706
> > > > > > > ATOM 33 CA ALA 4 14.191 3.898 0.745
> > > > > > > ATOM 34 HA ALA 4 15.133 3.670 1.244
> > > > > > > ATOM 35 CB ALA 4 14.334 3.685 -0.741
> > > > > > > ATOM 36 1HB ALA 4 80.378-330.692 -73.240
> > > > > > > ATOM 37 2HB ALA 4 13.457 4.049 -1.277
> > > > > > > ATOM 38 3HB ALA 4 15.200 4.236 -1.108
> > > > > > > ATOM 39 C ALA 4 13.926 5.388 1.011
> > > > > > > ATOM 40 O ALA 4 14.864 6.200 1.003
> > > > > > > ATOM 41 N ALA 5 12.739 5.814 1.451
> > > > > > > ATOM 42 H ALA 5 12.045 5.089 1.560
> > > > > > > ATOM 43 CA ALA 5 12.360 7.109 2.018
> > > > > > > ATOM 44 HA ALA 5 13.242 7.714 2.229
> > > > > > > ATOM 45 CB ALA 5 11.627 7.902 0.925
> > > > > > > ATOM 46 1HB ALA 5 10.831 7.360 0.413
> > > > > > > ATOM 47 2HB ALA 5 11.254 8.883 1.220
> > > > > > > ATOM 48 3HB ALA 5 12.343 8.027 0.113
> > > > > > > ATOM 49 C ALA 5 11.526 6.921 3.332
> > > > > > > ATOM 50 O ALA 5 10.378 7.256 3.490
> > > > > > > TER
> > > > > > > END
> > > > > > >
> > > > > > >
> > > > > > > Groupfile :
> > > > > > >
> > > > > > > -O -i inp.300 -p 1a.prmtop -c restrt.300 -r r1.300 -o out.300
> > > > > > > -rem 1
> > > > > > > -x mdcrd.300
> > > > > > > -O -i inp.320 -p 1a.prmtop -c restrt.320 -r r1.320 -o out.320
> > > > > > > -rem 1
> > > > > > > -x mdcrd.320
> > > > > > > -O -i inp.340 -p 1a.prmtop -c restrt.340 -r r1.340 -o out.340
> > > > > > > -rem 1
> > > > > > > -x mdcrd.340
> > > > > > > -O -i inp.360 -p 1a.prmtop -c restrt.360 -r r1.360 -o out.360
> > > > > > > -rem 1
> > > > > > > -x mdcrd.360
> > > > > > > -O -i inp.380 -p 1a.prmtop -c restrt.380 -r r1.380 -o out.380
> > > > > > > -rem 1
> > > > > > > -x mdcrd.380
> > > > > > > -O -i inp.400 -p 1a.prmtop -c restrt.400 -r r1.400 -o out.400
> > > > > > > -rem 1
> > > > > > > -x mdcrd.400
> > > > > > > -O -i inp.420 -p 1a.prmtop -c restrt.420 -r r1.420 -o out.420
> > > > > > > -rem 1
> > > > > > > -x mdcrd.420
> > > > > > > -O -i inp.440 -p 1a.prmtop -c restrt.440 -r r1.440 -o out.440
> > > > > > > -rem 1
> > > > > > > -x mdcrd.440
> > > > > > >
> > > > > > >
> > > > > >
> > > > > > --
> > > > > >
> > > > > > ===================================================================
> > > > > > Carlos L. Simmerling, Ph.D.
> > > > > > Associate Professor Phone: (631) 632-1336
> > > > > > Center for Structural Biology Fax: (631) 632-1555
> > > > > > CMM Bldg, Room G80
> > > > > > Stony Brook University E-mail: carlos.simmerling.gmail.com
> > > > > > Stony Brook, NY 11794-5115 Web: http://comp.chem.sunysb.edu
> > > > > >
> > > > > > ===================================================================
> > > > > >
> > > > >
> > > > >
> > > > >
> > > > > --
> > > > > Ashish V. Sangwai
> > > > > 6214 York Street
> > > > > New Orleans
> > > > > LA 70125
> > > > >
> > > >
> > > >
> > > >
> > > > --
> > > > ===================================================================
> > > > Carlos L. Simmerling, Ph.D.
> > > > Associate Professor Phone: (631) 632-1336
> > > > Center for Structural Biology Fax: (631) 632-1555
> > > > CMM Bldg, Room G80
> > > > Stony Brook University E-mail: carlos.simmerling.gmail.com
> > > > Stony Brook, NY 11794-5115 Web: http://comp.chem.sunysb.edu
> > > > ===================================================================
> > > >
> > > >
> > >
> > >
> > >
> > --
> > Dr. Adrian E. Roitberg
> > Associate Professor
> > Quantum Theory Project and Department of Chemistry
> >
> > University of Florida PHONE 352 392-6972
> > P.O. Box 118435 FAX 352 392-8722
> > Gainesville, FL 32611-8435 Email adrian.qtp.ufl.edu
> >
> > ============================================================================
> >
> > To announce that there must be no criticism of the president,
> > or that we are to stand by the president right or wrong,
> > is not only unpatriotic and servile, but is morally treasonable
> > to the American public."
> > -- Theodore Roosevelt
> > -----------------------------------------------------------------------
> > The AMBER Mail Reflector
> > To post, send mail to amber.scripps.edu
> > To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
> >
>
>
>
> --
> ===================================================================
> Carlos L. Simmerling, Ph.D.
> Associate Professor Phone: (631) 632-1336
> Center for Structural Biology Fax: (631) 632-1555
> CMM Bldg, Room G80
> Stony Brook University E-mail: carlos.simmerling.gmail.com
> Stony Brook, NY 11794-5115 Web: http://comp.chem.sunysb.edu
> ===================================================================
>



-- 
Ashish V. Sangwai
6214 York Street
New Orleans
LA 70125






-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu

Received on Wed Apr 23 2008 - 06:07:40 PDT
Custom Search