AMBER Archive (2008)

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

From: Adrian Roitberg (roitberg_at_qtp.ufl.edu)
Date: Mon Apr 21 2008 - 09:56:33 CDT

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