AMBER Archive (2009)

Subject: Re: [AMBER] varying positional restraints with PMEMD

From: case (
Date: Sun Nov 08 2009 - 11:37:53 CST

On Fri, Nov 06, 2009, Julian Garrec wrote:
> I'm wondering if it is possible to make _positional_ restraints varying
> linearly with the simulation time during a single PMEMD run.
> Here is one example of input file I have:
> title
> &cntrl
> imin=0, dt=0.0015, nstlim=200000,
> ntpr=100, ntwx=1000, ntwr=100,
> ioutfm=1,
> irest=1, ntx=5,
> ntc=2, ntf=2,
> cut=10,
> ntb=2, ntp=1, pres0=1.0, taup=2.0,
> ntt=3, gamma_ln=1.0, ig=-1, temp0=300.0,
> ntr=1, nmropt=1
> /
> &wt
> /
> DISANG=restraints.inp
> /
> constrain the protein
> 50.0
> RES 1 206
> I would like to decrease slowly the force constant from 50 to 0
> kcal/mol. Their are several paper of simmerling et al. whre they use a
> similar approach to equilibrate their system with PMEMD. (see, eg, J.
> Am. Chem. Soc. (2008) 130:7184)

> I think using a series of short runs with decreasing force constant is
> not safe for two reasons:
> - first, according to my little experience in the field, using always
> the same reference file (eg, the cristallographic structure) for a set
> of md run in AMBER leads systematically to a very short but huge
> increase of the kinetic energy at the beginning of each run (the
> temperature raises to say 600-800 K and then re-stabilizes quickly to
> the target value). I guess this comes from the finite precision of the
> coordinates that are stored in the .rst file. One solution is then to
> use the coordinates of the final step of run number n as a reference for
> run number n+1. However, this would make the reference changing
> arbitrarily and without any control, especially when the force constant
> will become little.

This probably has nothing to do with the finite precision of the restart file,
but rather with the fact that you are presumably using constant pressure
simulations, where the moving protein is being scaled, but the crystal
structure reference remains constant. You probably have to choose your poison

> - second, I am using Langevin dynamics. Even if I request the seed for
> the random number generator to be based on clock time setting ig=-1, the
> state of the thermostat is not sent from one job (n) to another (n+1).
> So to me using very short Langevin runs, whose duration would be of the
> same order of magnitude than the period between two random collisions,
> does not sound very safe.

I don't understand your argument here. The idea of the "frequency of random
collisions" is just a fiction -- a random kick is given to all atoms at every
step. It's not clear that there is a problem with using short Langevin runs.

> One alternative would be to use sander to do this, and then switch to
> PMEMD to run without restraints. However, as far as I know, PME is
> implemented differently in pmemd and sander (or, at least, default
> parameters are different), so that switching from one code to another
> should be done carefully and would require additional relaxation of the
> system. This I don't want :)

This argument is not correct: pmemd and sander should produce identical
results, and there should be neither any change in default parameters nor any
need to do additional relaxation. [Nor am I sure exactly what would be
possible in terms of restraints with sander vs. pmemd -- the popular sorts of
restraints are implemented in both codes.]


AMBER mailing list