Hi,
One of my students had this happen a few times with her REMD simulations of small peptides (<9 amino acids). By looking at each coordinate file (ie. each simulation at a different temperature) we saw that one or two hydrogens would "fly" off the molecule at the highest temperature, and only at the highest temperature. This usually occurred when we ran a simulation on a smaller system (ie. less degrees of freedom, eg. 5 amino acids) using a temperature set developed for a larger system (eg. 8 amino acids). --- I know, a new temperature set should be computed each time for each new system that has a different degrees of freedom, but this simplified the procedure for studying a large number of small peptides. --- My solution was to reduce the high temperatures by doing several simulations at various temperatures, until I had found a stable trajectory. I then recomputed the REMD temperatures range using the new highest temperature. The resulting REMD simulation was stable in that they did not have any "flying" hydrogens at any temperature. I could only rationalize this by thinking that too much energy was going into the system at the highest temperature. My simulation parameters were ntc=2, ntf=2. So, what Dr. Roitberg says concerning SHAKE makes sense to me. Also, reducing dt worked for a few systems, but that was by no means a fix for all of our systems.
This is my $0.02 worth, because I don't think I have enough for $0.03.
Cheers,
Karl
____________________________________
Karl N. Kirschner, Ph.D.
Center for Molecular Design, Co-Director
Hamilton College, Clinton NY 13323
____________________________________
----- Original Message -----
From: Adrian Roitberg <roitberg.qtp.ufl.edu>
Date: Monday, April 21, 2008 5:00 pm
Subject: Re: AMBER: hydrogens are flying in replica exchange simulations
To: amber.scripps.edu
> 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
>
-----------------------------------------------------------------------
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:35 PDT
