AMBER Archive (2008)

Subject: Re: AMBER: RE: About restart amber

From: Francesco Pietra (chiendarret_at_gmail.com)
Date: Tue Nov 25 2008 - 15:47:35 CST

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