AMBER Archive (2008)

Subject: Re: AMBER: hydrogens are flying in replica exchange simulations

From: Karl Kirschner (kkirschn) (kkirschn_at_hamilton.edu)
Date: Mon Apr 21 2008 - 10:47:37 CDT

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_at_qtp.ufl.edu>
Date: Monday, April 21, 2008 5:00 pm
Subject: Re: AMBER: hydrogens are flying in replica exchange simulations
To: amber_at_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_at_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_at_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_at_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_at_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_at_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_at_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_at_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_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu
>

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