AMBER Archive (2008)

Subject: Re: AMBER: RE: About restart amber

From: Robert Duke (rduke_at_email.unc.edu)
Date: Tue Nov 25 2008 - 15:55:15 CST

Sorry, but you have the pdb and prmtop, and it will take a little time to
figure it out. I find it hard to believe that with a little
experimentation, it can't be done with search/find over some range of atoms
or residues...

----- Original Message -----
From: "Francesco Pietra" <chiendarret_at_gmail.com>
To: <amber_at_scripps.edu>
Sent: Tuesday, November 25, 2008 4:47 PM
Subject: Re: AMBER: RE: About restart amber

> OK, could you please give me a hint how to change for pmemd the
> following restraints (used for the pore region in sander.MPI)
>
> density equilibration with restraint
> on protein, ligand, and POPC polar head
> &cntrl
> imin=0, irest=1, ntx=5,
> nstlim=250000, dt=0.002,
> cut=10, ntb=2, ntp=1, taup=1.0,
> ntc=2, ntf=2,
> ntpr=500, ntwx=500,
> ntt=3, gamma_ln=2.0,
> tempi=300.0, temp0=300.0,
> ntr=1,
> restraintmask=":79-341 | :POP_at_O2, P1, O3, O4, O1, C15, C11, N, C12, C13,
> C14"
> restraint_wt=32,
> /
>
> where POPC = 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine
>
> Thanks
> francesco
>
>
> On Tue, Nov 25, 2008 at 10:23 PM, Robert Duke <rduke_at_email.unc.edu> wrote:
>> Okay. Regarding your future work, the group specification stuff should
>> allow you to select specific atoms, specific residues, or search/find
>> based
>> on specific characteristics including atom type (I know atom name is
>> different). I have not done this sort of stuff extensively myself, but
>> the
>> old group specification stuff in pmemd should work as well as it did in
>> sander. I thought the description of this stuff had been dropped in
>> amber
>> 10, but it is actually now described in the manual in a separate section
>> entitled "GROUP Specification". Earlier releases (9, 8, etc.) have
>> basically the same info in the appendices. Anyway, I would be surprised
>> if
>> you can't find some way to restrain the lipid (not to say that I have not
>> been surprised before...)
>> - Bob
>> ----- Original Message ----- From: "Francesco Pietra"
>> <chiendarret_at_gmail.com>
>> To: <amber_at_scripps.edu>
>> Sent: Tuesday, November 25, 2008 3:59 PM
>> Subject: Re: AMBER: RE: About restart amber
>>
>>
>>> Thanks for this very instructive piece. However, I decided to abandon
>>> GB for this large system. As Bud Dobson warned on a parallel thread,
>>> it was already well known that GB - at its present implementation -
>>> becomes much slower than explicit solvent MD. Add the deformations
>>> introduced by GB, which would require membrane surrogates, should it
>>> work.... I can not blame than myself for having assumed GB faster
>>> without checking the archives.
>>>
>>> MD in an explicit membrane is running fast with pmemd for the pore
>>> region of this protein. Once ready, I'll try the whole system with the
>>> pore region in a membrane and the rest in water. If too slow for my
>>> means, I'll leave the input for better equipped people. Probably I'll
>>> face too high obstacles at the equilibration stage (which is quite
>>> long because of the membrane) as I can't use pmemd to restrain the
>>> lipid. Only restraintmask - as far as I know - can deal with atom
>>> names (perhaps a hint for a future version of pmemd).
>>>
>>> francesco
>>>
>>>
>>>
>>> On Tue, Nov 25, 2008 at 3:27 PM, Robert Duke <rduke_at_email.unc.edu>
>>> wrote:
>>>>
>>>> Well, it blows up after a recentering operation triggered at step 1000
>>>> by
>>>> the default nscm value... The reason it blows up is that shake can't
>>>> find a
>>>> reasonable solution for all bondlengths because some atom or atoms has
>>>> moved
>>>> too much relative to other atoms. I have no idea what is specifically
>>>> happening here but I would get a high resolution trajectory (ntwx = 50
>>>> or
>>>> 100) and look at it in vmd. It could be that partial restraints and
>>>> repositioning don't work well together, putting really unnatural forces
>>>> on
>>>> things. I would try setting nscm to a larger number (>nstlim), just to
>>>> see
>>>> what is going on, and look at the trajectories (and perhaps also
>>>> something
>>>> smaller, like 100). When I use restraints, I am always restraining
>>>> whole
>>>> molecules, so I don't have direct experience with what you are doing.
>>>> Another point: if your restraints are on a structure that was not
>>>> centered
>>>> to start, I would think some really nasty things could happen in the
>>>> first
>>>> recentering. Under ewald, repositioning is purely translational (nscm
>>>> option), and it used to fix an energy equipartition problem (Tom
>>>> Cheatham's
>>>> lab has a paper on the "ewald flying ice cube"; J Comput Chem 19, pp
>>>> 726-740). I am not sure about the generalized Born equivalent here,
>>>> but
>>>> I
>>>> would guess that turning off repositioning in a restrained run for a
>>>> few
>>>> thousand steps would not be the end of the world. I would still very
>>>> much
>>>> want to look at the trajectories though...
>>>> Regards - Bob
>>>>
>>>> ----- Original Message ----- From: "Francesco Pietra"
>>>> <chiendarret_at_gmail.com>
>>>> To: <amber_at_scripps.edu>
>>>> Sent: Tuesday, November 25, 2008 3:00 AM
>>>> Subject: Re: AMBER: RE: About restart amber
>>>>
>>>>
>>>>> Answering to myself, setting additional restraints (for the pore
>>>>> region, in the distorted conformation for one of the chains as it was
>>>>> left from previous run), was met by vlimit exceeded. I suspect that
>>>>> the procedure I followed is basically wrong. I started by assuming
>>>>> that the the extracellular part was treated adequately by GB, while I
>>>>> was not interested here in the pore region, which was distorted as to
>>>>> one of the chains.
>>>>>
>>>>> I am unable to detect any major deviation in EELEC WDW etc, so that I
>>>>> am here for help. General settings for the MD simulation are in
>>>>> previous mail on this thread, while the output, comprising the new
>>>>> restraints (for pmemd), are here:
>>>>>
>>>>> Keep pore initial chain A restrained
>>>>>
>>>>> GROUP 4 HAS HARMONIC CONSTRAINTS 32.00000
>>>>> GRP 4 RES 1 TO 41
>>>>> Number of atoms in this group = 686
>>>>> ----- READING GROUP 5; TITLE:
>>>>> Keep pore final chain A restrained
>>>>>
>>>>> GROUP 5 HAS HARMONIC CONSTRAINTS 32.00000
>>>>> GRP 5 RES 381 TO 426
>>>>> Number of atoms in this group = 725
>>>>> ----- READING GROUP 6; TITLE:
>>>>> Keep pore initial chain B restrained
>>>>>
>>>>> etc for all other chain stretches in the pore.
>>>>>
>>>>>
>>>>> --------------------------------------------------------------------------------
>>>>> 3. ATOMIC COORDINATES AND VELOCITIES
>>>>>
>>>>>
>>>>> --------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> begin time read from input coords = 118.000 ps
>>>>>
>>>>> Number of triangulated 3-point waters found: 0
>>>>> | Dynamic Memory, Types Used:
>>>>> | Reals 771184
>>>>> | Integers 1594255
>>>>>
>>>>> | Running AMBER/MPI version on 8 nodes
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> --------------------------------------------------------------------------------
>>>>> 4. RESULTS
>>>>>
>>>>>
>>>>> --------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 100 TIME(PS) = 118.200 TEMP(K) = 293.29 PRESS =
>>>>> 0.0
>>>>> Etot = -14862.0902 EKtot = 15182.4903 EPtot -30044.5805
>>>>> BOND = 4208.2193 ANGLE = 10364.8917 DIHED =
>>>>> 13636.2332
>>>>> 1-4 NB = 4615.5997 1-4 EEL = 58027.8454 AALS -9952.6334
>>>>> EELEC = -84702.6203 EGB = -27632.9571 RESTRAINT =
>>>>> 1390.8410
>>>>> EAMBER (non-restraint) = -31435.4215
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 200 TIME(PS) = 118.400 TEMP(K) = 292.26 PRESS =
>>>>> 0.0
>>>>> Etot = -14549.8615 EKtot = 15128.9530 EPtot -29678.8145
>>>>> BOND = 4249.1002 ANGLE = 10497.9349 DIHED =
>>>>> 13714.8644
>>>>> 1-4 NB = 4622.1647 1-4 EEL = 57924.9207 AALS -9907.1452
>>>>> EELEC = -84454.0644 EGB = -27791.1146 RESTRAINT =
>>>>> 1464.5248
>>>>> EAMBER (non-restraint) = -31143.3393
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 300 TIME(PS) = 118.600 TEMP(K) = 296.53 PRESS =
>>>>> 0.0
>>>>> Etot = -14495.6067 EKtot = 15350.1320 EPtot -29845.7387
>>>>> BOND = 4191.7755 ANGLE = 10416.8064 DIHED =
>>>>> 13702.0116
>>>>> 1-4 NB = 4605.7994 1-4 EEL = 57951.4365 AALS -9853.8882
>>>>> EELEC = -84508.8249 EGB = -27853.5374 RESTRAINT =
>>>>> 1502.6823
>>>>> EAMBER (non-restraint) = -31348.4210
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 400 TIME(PS) = 118.800 TEMP(K) = 298.18 PRESS =
>>>>> 0.0
>>>>> Etot = -14369.6628 EKtot = 15435.1364 EPtot -29804.7992
>>>>> BOND = 4144.3264 ANGLE = 10392.7962 DIHED =
>>>>> 13603.8654
>>>>> 1-4 NB = 4617.9524 1-4 EEL = 58109.1303 AALS -9882.6587
>>>>> EELEC = -84558.4946 EGB = -27780.5494 RESTRAINT =
>>>>> 1548.8329
>>>>> EAMBER (non-restraint) = -31353.6320
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 500 TIME(PS) = 119.000 TEMP(K) = 297.56 PRESS =
>>>>> 0.0
>>>>> Etot = -14401.4198 EKtot = 15403.3889 EPtot -29804.8086
>>>>> BOND = 4170.8469 ANGLE = 10424.8530 DIHED =
>>>>> 13619.7513
>>>>> 1-4 NB = 4575.5432 1-4 EEL = 57982.6068 AALS -9920.2602
>>>>> EELEC = -84607.0488 EGB = -27662.7872 RESTRAINT =
>>>>> 1611.6864
>>>>> EAMBER (non-restraint) = -31416.4950
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 600 TIME(PS) = 119.200 TEMP(K) = 297.37 PRESS =
>>>>> 0.0
>>>>> Etot = -14324.9156 EKtot = 15393.5581 EPtot -29718.4737
>>>>> BOND = 4068.7054 ANGLE = 10487.9072 DIHED =
>>>>> 13705.9594
>>>>> 1-4 NB = 4664.9561 1-4 EEL = 57972.6307 AALS -9837.0587
>>>>> EELEC = -84701.0042 EGB = -27673.2914 RESTRAINT =
>>>>> 1592.7218
>>>>> EAMBER (non-restraint) = -31311.1955
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 700 TIME(PS) = 119.400 TEMP(K) = 298.96 PRESS =
>>>>> 0.0
>>>>> Etot = -14226.3546 EKtot = 15475.6550 EPtot -29702.0096
>>>>> BOND = 4114.8803 ANGLE = 10518.2858 DIHED =
>>>>> 13716.4825
>>>>> 1-4 NB = 4642.9989 1-4 EEL = 57873.1198 AALS -9909.4567
>>>>> EELEC = -84559.0840 EGB = -27685.0454 RESTRAINT =
>>>>> 1585.8092
>>>>> EAMBER (non-restraint) = -31287.8188
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 800 TIME(PS) = 119.600 TEMP(K) = 297.23 PRESS =
>>>>> 0.0
>>>>> Etot = -14338.7373 EKtot = 15386.1646 EPtot -29724.9020
>>>>> BOND = 4204.7480 ANGLE = 10420.4349 DIHED =
>>>>> 13658.4428
>>>>> 1-4 NB = 4633.6685 1-4 EEL = 57948.2947 AALS -9867.9684
>>>>> EELEC = -84626.2242 EGB = -27715.7199 RESTRAINT =
>>>>> 1619.4217
>>>>> EAMBER (non-restraint) = -31344.3237
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>>
>>>>> NSTEP = 900 TIME(PS) = 119.800 TEMP(K) = 299.42 PRESS =
>>>>> 0.0
>>>>> Etot = -14227.7282 EKtot = 15499.4711 EPtot -29727.1993
>>>>> BOND = 4169.8501 ANGLE = 10520.4113 DIHED =
>>>>> 13684.0647
>>>>> 1-4 NB = 4610.5941 1-4 EEL = 57995.3457 AALS -9865.9772
>>>>> EELEC = -84707.6272 EGB = -27695.0750 RESTRAINT =
>>>>> 1561.2141
>>>>> EAMBER (non-restraint) = -31288.4135
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>> | RE_POSITION Moving by 0.499075 0.005780 1.333641
>>>>>
>>>>> NSTEP = 1000 TIME(PS) = 120.000 TEMP(K) = 299.15 PRESS =
>>>>> 0.0
>>>>> Etot = -14292.8263 EKtot = 15485.3583 EPtot -29778.1846
>>>>> BOND = 4142.2201 ANGLE = 10420.1202 DIHED =
>>>>> 13702.3912
>>>>> 1-4 NB = 4625.8817 1-4 EEL = 58004.4484 AALS -9979.7547
>>>>> EELEC = -84469.6363 EGB = -27779.4524 RESTRAINT =
>>>>> 1555.5971
>>>>> EAMBER (non-restraint) = -31333.7817
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------------
>>>>>
>>>>> vlimit exceeded for step 1025; vmax = 29.6260
>>>>>
>>>>> Coordinate resetting cannot be accomplished,
>>>>> deviation is too large
>>>>> iter_cnt, my_bond_idx, i and j are : 4 984 578 576
>>>>>
>>>>> Thanks for any advice
>>>>> francesco
>>>>>
>>>>> On Mon, Nov 24, 2008 at 10:22 PM, Francesco Pietra
>>>>> <chiendarret_at_gmail.com> wrote:
>>>>>>
>>>>>> Hi Ross:
>>>>>> As the pore region of the transmembrane protein I am dealing with
>>>>>> gets
>>>>>> deformed on long MD simulation under GB conditions (while the
>>>>>> extracellular portion behaves normally), I would like to continue
>>>>>> with
>>>>>> GB by restraining the pore region (which was not comprised on last
>>>>>> run). In fact, I am not interested in the pore region in these
>>>>>> simulations. Is it correct to use as reference (-ref ...rst) the
>>>>>> restart file of previous run while adding now the new restraints?
>>>>>>
>>>>>> From a trial with nstlim=1 all restraints (given as GROUP for pmemd)
>>>>>> are read, as from the out file. Is anything subtly wrong that can't
>>>>>> be
>>>>>> detected this way?
>>>>>>
>>>>>> Thanks
>>>>>> francesco
>>>>>>
>>>>>> On Sat, Nov 22, 2008 at 6:18 PM, Ross Walker <ross_at_rosswalker.co.uk>
>>>>>> wrote:
>>>>>>>
>>>>>>> Hi Francesco,
>>>>>>>
>>>>>>>> When restraints are imposed, which reference file should be used in
>>>>>>>> running sander or pmemd again? I guess the "-ref ...rst" pertaining
>>>>>>>> to
>>>>>>>> the crashed job. Correct? Does that involve any further
>>>>>>>> complication?
>>>>>>>
>>>>>>> If you want an actual restart then you should restrain to the same
>>>>>>> restart file as the original run. Note, however, that when running
>>>>>>> NPT
>>>>>>> I
>>>>>>> think this leads to instabilities since there is some kind of
>>>>>>> adjustment of
>>>>>>> the box coordinates done that then means the restraints are too
>>>>>>> large.
>>>>>>> I
>>>>>>> have never actually nailed this problem down but essentially it
>>>>>>> means
>>>>>>> for
>>>>>>> constant pressure that you can only ever restrain to the same
>>>>>>> restart
>>>>>>> file
>>>>>>> as you are using for the input coordinates (-i). I think the problem
>>>>>>> originates in the fact that when you restart the box coordinates in
>>>>>>> your
>>>>>>> inpcrd file do not match the box coordinates in your -ref file.
>>>>>>>
>>>>>>> At least this problem existed back in Amber 7 - maybe it got fixed,
>>>>>>> I
>>>>>>> haven't tried it in a while. I avoid running constant pressure runs
>>>>>>> with
>>>>>>> restraints anyway because I figure the restraints are effectively
>>>>>>> giving you
>>>>>>> artificial behavior anyway. E.g. suppose something wants to elongate
>>>>>>> but the
>>>>>>> restraints stop it - really the box should change size to
>>>>>>> accommodate
>>>>>>> this
>>>>>>> elongation.
>>>>>>>
>>>>>>> That is probably complicating things though. For NVE and NVT you
>>>>>>> should
>>>>>>> just be able to use the same -ref as you did for the previous run.
>>>>>>>
>>>>>>>> In view of non minor work in restarting from a crashed job while
>>>>>>>> taking everything into account (not to mention the unluky situation
>>>>>>>> that the machine was just writing the rst file), is any major
>>>>>>>> drawback
>>>>>>>> in making the various steps of MD short? If problems exist, how to
>>>>>>>> determine how much short?
>>>>>>>
>>>>>>> You should probably read the following 'very interesting' paper:
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> http://pubs.acs.org/cgi-bin/abstract.cgi/jctcce/2008/4/i10/abs/ct8002173.html
>>>>>>>
>>>>>>> This looks at problems with repeating random number sequences and
>>>>>>> Langevin / Andersen dynamics, the equivalent behavior that occurs
>>>>>>> when
>>>>>>> you
>>>>>>> continually restart a simulation so effectively running very short
>>>>>>> MD
>>>>>>> runs.
>>>>>>> The key problem, with continuous short restarts, is that AMBER (and
>>>>>>> CHARMM
>>>>>>> and most other codes I think) do not strictly do proper restarts.
>>>>>>> That
>>>>>>> is
>>>>>>> they do not preserve the state of the random number stream. Hence if
>>>>>>> you
>>>>>>> just keep blindly restarting you use the exact same set of random
>>>>>>> numbers
>>>>>>> each time. For NVE this is not a problem but for Langevin dynamics
>>>>>>> it
>>>>>>> can
>>>>>>> become acute because you introduce correlation (in the reuse of the
>>>>>>> random
>>>>>>> number stream) in the 'random' forces applied and this can lead to
>>>>>>> all
>>>>>>> kinds
>>>>>>> of weird behavior. For reasonably long runs between restarts the
>>>>>>> correlation
>>>>>>> is 'weak' enough that it does not adversely affect the dynamics but
>>>>>>> for
>>>>>>> short gaps between restarts you can get strange behavior.
>>>>>>>
>>>>>>> The main point to take home here is that you should probably
>>>>>>> 'ALWAYS'
>>>>>>> be
>>>>>>> changing the random number seed (ig=) whenever you restart a
>>>>>>> simulation - in
>>>>>>> fact you should probably never run any two MD runs with the same
>>>>>>> seed
>>>>>>> (except for testing / debugging). A simple solution to this (in
>>>>>>> AMBER
>>>>>>> 10) is
>>>>>>> to set ig=-1 and then it will use the wallclock time in microseconds
>>>>>>> to seed
>>>>>>> the random number generator.
>>>>>>>
>>>>>>> We should probably be more explicit about this in the manual /
>>>>>>> tutorials
>>>>>>> etc but the reality is that we are only now beginning to realize
>>>>>>> that
>>>>>>> such
>>>>>>> effects exist.
>>>>>>>
>>>>>>> Good luck,
>>>>>>> Ross
>>>>>>>
>>>>>>>
>>>>>>> /\
>>>>>>> \/
>>>>>>> |\oss Walker
>>>>>>>
>>>>>>> | Assistant Research Professor |
>>>>>>> | San Diego Supercomputer Center |
>>>>>>> | Tel: +1 858 822 0854 | EMail:- ross_at_rosswalker.co.uk |
>>>>>>> | http://www.rosswalker.co.uk | PGP Key available on request |
>>>>>>>
>>>>>>> Note: Electronic Mail is not secure, has no guarantee of delivery,
>>>>>>> may
>>>>>>> not be read every day, and should not be used for urgent or
>>>>>>> sensitive
>>>>>>> issues.
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> -----------------------------------------------------------------------
>>>>>>> The AMBER Mail Reflector
>>>>>>> To post, send mail to amber_at_scripps.edu
>>>>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the
>>>>>>> email)
>>>>>>> to majordomo_at_scripps.edu
>>>>>>>
>>>>>>
>>>>> -----------------------------------------------------------------------
>>>>> The AMBER Mail Reflector
>>>>> To post, send mail to amber_at_scripps.edu
>>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>>>> to majordomo_at_scripps.edu
>>>>>
>>>>
>>>> -----------------------------------------------------------------------
>>>> The AMBER Mail Reflector
>>>> To post, send mail to amber_at_scripps.edu
>>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>>> to majordomo_at_scripps.edu
>>>>
>>> -----------------------------------------------------------------------
>>> The AMBER Mail Reflector
>>> To post, send mail to amber_at_scripps.edu
>>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>>> to majordomo_at_scripps.edu
>>>
>>
>> -----------------------------------------------------------------------
>> The AMBER Mail Reflector
>> To post, send mail to amber_at_scripps.edu
>> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
>> to majordomo_at_scripps.edu
>>
> -----------------------------------------------------------------------
> The AMBER Mail Reflector
> To post, send mail to amber_at_scripps.edu
> To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
> to majordomo_at_scripps.edu
>

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
      to majordomo_at_scripps.edu