Re: [AMBER] help for TMD

From: Jason Swails <jason.swails.gmail.com>
Date: Wed, 12 Sep 2012 08:25:56 -0400

On Tue, Sep 11, 2012 at 11:41 PM, Yanyan ZHU <zhuyanzily.gmail.com> wrote:

> I have received your explanation and advice for my puzzles ,as following:
> You should move carefully step by step. It looks like you are trying to
> equilibrate the system and run targetted MD at the same time. Try the
> following:
>
> 1. Make sure you can run ordinary dynamics. Turn off targeted dynamics,
> set ntpr=1 and nstlim=50. Get your system equilibrated to 300K. (If you
> have
> problems at this stage, then you have a much simpler debugging strategy.)
>
> 2. You input file sets ntr=1, with a very large restraint weight. For that
> part of the protein, nothing will be able to move. This will probably need
> to
> be changed.
>
> 3. Try very simple targeted MD runs once you have an equilibrated protein.
> To get started, still set ntpr=1 and nstlim to some smallish value (say
> 100).
> Use the simplest restraints you can think of and study the output carefully
> to
> try to understand exactly what is going on.
>
> 4. Your output says "old style PARM file read". This is somewhat
>
>
> I have tried to solve my problems as you suggested. Firstly, I minimize
> my system and my input file is :
>
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> initial minimisation of the protein
> &cntrl
> imin=1,
> maxcyc=50000,
> ncyc=25000,
>

This is too many minimization steps if you are just trying to get a good
starting structure for dynamics. The point of minimization is to relieve
bad contacts and start from a reasonable structure -- this can be
accomplished easily within 1000 steps. I usually use maxcyc=1000 and
ncyc=100 (and even that is probably more than enough).


> cut=12.0,
> ntb=0,
>

This does not look right. Your later input files use igb=1 and cut=999.0.
 Use those values in your minimization as well.

ntpr=20,
> ntr=1,
> restraintmask=':29-342.CA,N,C',
> restraint_wt=200,
> /
>
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>
> Secondly, I heat my system through six steps separately with a very short
> time and my input file is :
>
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> Stage 1 heating of T 0 to 50K
> &cntrl
> imin=0, irest=0, ntx=1,
> nstlim=50, dt=0.0005,
>

This heating stage takes place over 25 femtoseconds. This is way too short
for most of your heating steps. When I heat, I typically start out at 10K
to 50K (depending on how comfortable I am with my starting structure), and
ramp it up slowly to my target temperature using nmropt=1.


> ntc=2, ntf=2,
> ntt=1, tautp=1.0,
> tempi=0.0, temp0=50.0,
> ntpr=1, ntwx=1,
> ntb=0, igb=1,
> ntr=1,
> restraintmask=':29-342.CA,N,C',
> restraint_wt=200,
>

These restraints are too strong in my opinion (they're almost as stiff as a
hydrogen-heavy atom bond!) Consider using something 2 orders of magnitude
lower. I usually use restraint weights of 2 to 5 kcal/mol, which is easily
sufficient to keep the restraintmask in place.


> After doing that ,I run equilibration and my input file is as following:
>
> +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> Stage 2 equilibration
> &cntrl
> imin=0, irest=1, ntx=5,
> nstlim=10, dt=0.001,
>

This equilibration is 10 fs. So far, you've performed 6 (?) heating stages
of 25 fs and this equilibration stage of 10 fs, meaning you have done 0.16
ps of simulation so far. Your system is likely not heated or equilibrated
(even for examples you have to do more simulation than this).


> ntc=2, ntf=2,
> ntt=1, tautp=0.5,
> tempi=300.0, temp0=300.0,
> ntpr=1, ntwx=1,
> ntb=0,igb=1,ntr=1,
> restraintmask=':29-342.CA,N,C',
> restraint_wt=10,
> cut=999.,rgbmax=999.
> /
>
> +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> The corresponding output file is :
>
> +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> Amber 9 SANDER 2006
> [snip]
> begin time read from input coords = 0.150 ps
> Number of triangulated 3-point waters found: 0
>
> --------------------------------------------------------------------------------
> 4. RESULTS
>
> --------------------------------------------------------------------------------
>
> NSTEP = 1 TIME(PS) = 0.151 TEMP(K) = 28.05 PRESS =
> 0.0
> Etot = 13530.8188 EKtot = 311.1434 EPtot =
> 13219.6754
> BOND = 18079.7784 ANGLE = 2532.7846 DIHED =
> 2854.3092
> 1-4 NB = 1102.3720 1-4 EEL = 12507.0613 VDWAALS =
> -1803.3055
> EELEC = -19456.1610 EGB = -2597.1637 RESTRAINT =
> 0.0000
>
> ------------------------------------------------------------------------------
>

Look at your temperature on this step -- it is 28.05 K (and heats up to 300
K by the end of 10 ps).


> NSTEP = 10 TIME(PS) = 0.160 TEMP(K) = 300.04 PRESS =
> 0.0
> Etot = 13555.1046 EKtot = 3327.6258 EPtot =
> 10227.4788
> BOND = 14396.7904 ANGLE = 3160.5024 DIHED =
> 2906.7961
> 1-4 NB = 1133.7874 1-4 EEL = 12475.7685 VDWAALS =
> -1893.4800
> EELEC = -19404.7164 EGB = -2624.8357 RESTRAINT =
> 76.8662
> EAMBER (non-restraint) = 10150.6126
>
> ------------------------------------------------------------------------------
> [snip]
>


> When I set nstlim=90 ,my output file is:
>
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> [snip]
>
> ------------------------------------------------------------------------------
> vlimit exceeded for step 89; vmax = 23.6266
> NSTEP = 90 TIME(PS) = 0.240 TEMP(K) = 650.95 PRESS =
> 0.0
> Etot = 13425.9595 EKtot = 7219.3867 EPtot =
> 6206.5728
> BOND = 3056.6908 ANGLE = 3907.4349 DIHED =
> 3227.9425
> 1-4 NB = 1114.4126 1-4 EEL = 12253.3039 VDWAALS =
> -1723.3208
> EELEC = -18838.6379 EGB = -2808.7714 RESTRAINT =
> 6017.5183
> EAMBER (non-restraint) = 189.0545
>
> ------------------------------------------------------------------------------
>

Your temperature has skyrocketed. However, given the extremely short
timescale of your heating and equilibration phases, it does not surprise me
that your system has become unstable. My suggestions are to use much
longer time scales for heating (I would use at least 10s of picoseconds, if
not longer).

Also, I would suggest against using the Berendsen thermostat (there are
many known problems with it). Consider Langevin dynamics as an
alternative. (ntt=3)


> We can see that the results have problems from NSEEP=87. I want to know
> the reason for this problem.
> I also want to konw after finish the equilibration ,can I begin run TMD .
>

Like Dave said before, I would suggest making sure you can run standard
dynamics successfully before attempting something more advanced like TMD.


> Besides that,I also have wonder aboutthe parallel command. I use the
> following commmand ,however ,it can't run successfully: mpirun -np 4
> $AMBERHOME/exe/sander.MPI ....................
>

We need more information to even try to guess what the problem is. These
issues are likely subtle errors in how the user is using the MPI installed
on the system (either Amber wasn't compiled correctly or you're using the
wrong MPI, etc. -- there is a long list of possible issues). I also
suggest trying to tackle one problem at a time ;).

HTH,
Jason

-- 
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Sep 12 2012 - 05:30:04 PDT
Custom Search