AMBER Archive (2009)

Subject: [AMBER] varying positional restraints with PMEMD

From: Julian Garrec (julian.garrec_at_epfl.ch)
Date: Thu Nov 05 2009 - 18:43:35 CST

Dear AMBER users,

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
TYPE='END',
/
DISANG=restraints.inp
/
constrain the protein
50.0
RES 1 206
END
END

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.

- 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.

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 :)

Any suggestion concerning the input file for PMEMD? Or could anyone tell
me which routine I need to modify to do what I want ?

Best regards,

Julian

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber