Re: AMBER: hydrogens are flying in replica exchange simulations

From: Carlos Simmerling <carlos.simmerling.gmail.com>
Date: Tue, 22 Apr 2008 14:57:11 -0400

Ashish,
this is missing some of the input files, such as the coordinates
inp.380, inp.400 and inp.420. can you check and send
all of the files that you used? send directly to me, not to the whole list.
carlos

2008/4/21 Ashish Sangwai <ashishsangwai.gmail.com>:
> 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



-- 
===================================================================
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
===================================================================
-----------------------------------------------------------------------
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:56 PDT
Custom Search