Re: AMBER: hydrogens are flying in replica exchange simulations

From: Adrian Roitberg <roitberg.qtp.ufl.edu>
Date: Mon, 21 Apr 2008 12:03:06 -0400

Ah,
but SHAKE did fail !

 From the previous emails

"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,"


So, I stand by my diagnostic ;-)

a.


Carlos Simmerling 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
>>
>
>
>

-- 
                            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
Received on Wed Apr 23 2008 - 06:07:30 PDT
Custom Search