Dear Fabian,
setup protocols vary from person to person. In some cases that person
may even have a reason why they do things in a specific way :-) What
matters at the end is that you can do a stable MD and do not create
"funny" conformations during setup.
What I would do in your case is to carry out maybe 100 or 200 steps of
minimisation but that may not even be necessary. Heating I would do
over maybe just 5 or 10 thousand steps (and set istep2 less than
nstlim by a few thousand steps because, from experience, the final
temperature will not be reached. At least that's the case with
Berendsen).
What I would recommend is to run MD at constant pressure even if you do
not plan to do production at const p. The reason for this is to make
sure that your final density is sensible (when you use the TIP3P box
in tleap your density will probably be between 0.6 and 0.7 for a
system like yours). Say, 10-20000, should be fine for this but check
convergence of the density.
If you used restraints you would probably want to also discard the
first "few" steps after having released those. But the decision what is
production run and what is not should really be decided dynamically,
i.e. depending on the convergence behaviour of your data of interested.
Hope that helps,
Hannes.
On Thu, 3 Dec 2015 14:15:17 +0200
Fabian gmail <fabian.glaser.gmail.com> wrote:
> Dear Hannes,
>
> thanks a lot, to be honest it will greatly help you if you go over
> the input files I prepared for the heating, equil and and prod for
> one lambda, and let me now what to change…. I just was trying to
> following the protocol on Kaus et al (dx.doi.org/10.1021/ct400340s |
> J. Chem. Theory Comput. 2013, 9, 4131−4139) since they do exactly
> what I want…. but I don't’ have experience at all with TI…. so any
> suggestion will be more than welcomed.
>
> If you have the time, I would be very grateful if you could take a
> quick look and tell me what to correct or what is the protocol your
> recommend in this case. I already reduced the minimisation steps to
> 500.
>
> So regarding the equilibration time I should increment the time by
> how much ? I start with 500 ps for lambda 0, then 1 ps lambda 0.1
> then 1.5 0.2 etc.? Is that a good strategy?
>
> I paste here the four input files I prepared following Kaus paper
> (again my goal is to try and calculate solvation energy of a drug). I
> would be very grateful for any correction….This inputs are for
> clambda = 0, but I thought to use them for every lambda, now that you
> mention the void issue, I understand I need to increase the
> equilibration time.
>
> Here they are, in the order I thought to use it for every lambda,
> obviously changing the clambda value in each case see PBS file below.
>
> I will highly appreciate also if you could recommend how to make the
> analysis later, I understand they are several tools….
>
> Thanks!!
>
> Fabian
>
>
> ===========
> minimisation 500 steps
> &cntrl
> imin = 1,
> ntmin = 2,
> maxcyc = 500,
> ntpr = 20,
> ntwe = 20,
> ntb = 1,
> cut = 8.0
>
> clambda = 0.0,
> scalpha = 0.5,
> scbeta = 12.0,
> logdvdl = 0,
>
> icfe = 1,
> ifsc = 1,
> timask1='',
> timask2=':1',
> scmask1='',
> scmask2=':1',
> /
> &ewald
> /
>
>
> ============
> heating constant vol 500 ps
> &cntrl
> imin = 0,
> nstlim = 250000,
> irest = 0,
> ntx = 1,
> dt = 0.002,
> nmropt = 1,
>
> ntb=2,
> ntp=1,
> ntt = 3,
> gamma_ln = 2,
> pres0 = 1.0,
>
> temp0 = 300.0,
> tempi = 0.0,
> tautp = 1.0,
> cut = 8.0,
>
> ntc = 2,
> ntf = 1,
> ntxo=2,
> ioutfm = 1,
> iwrap = 1,
>
> ntwe = 1000,
> ntwx = 1000,
> ntpr = 1000,
> ntwr = 5000,
>
> clambda = 0.0,
> scalpha = 0.5,
> scbeta = 12.0,
> logdvdl = 0,
> icfe = 1,
> ifsc = 1,
> timask1='',
> timask2=':1',
> scmask1='',
> scmask2=':1',
>
>
> /
> &ewald
> /
>
> /
> &wt TYPE='TEMP0', istep1=0, istep2=25000,
> value1=0.1, value2=300.0,
> /
> &wt TYPE='END'
> /
>
> ============
> heating constant vol 500 ps
> &cntrl
> imin = 0,
> nstlim = 250000,
> irest = 0,
> ntx = 1,
> dt = 0.002,
> nmropt = 1,
>
> ntb=2,
> ntp=1,
> ntt = 3,
> gamma_ln = 2,
> pres0 = 1.0,
>
> temp0 = 300.0,
> tempi = 0.0,
> tautp = 1.0,
> cut = 8.0,
>
> ntc = 2,
> ntf = 1,
> ntxo=2,
> ioutfm = 1,
> iwrap = 1,
>
> ntwe = 1000,
> ntwx = 1000,
> ntpr = 1000,
> ntwr = 5000,
>
> clambda = 0.0,
> scalpha = 0.5,
> scbeta = 12.0,
> logdvdl = 0,
> icfe = 1,
> ifsc = 1,
> timask1='',
> timask2=':1',
> scmask1='',
> scmask2=':1',
>
>
> /
> &ewald
> /
>
> /
> &wt TYPE='TEMP0', istep1=0, istep2=25000,
> value1=0.1, value2=300.0,
> /
> &wt TYPE='END'
> /
>
> ============
> Constant Volume production 5ns
> &cntrl
> imin = 0,
> nstlim = 2500000,
> irest = 1,
> ntx = 5,
> dt = 0.002,
>
> ntt = 3,
> temp0 = 300.0,
> gamma_ln = 2.0,
> ig = -1,
> ntc = 2,
> ntf = 1,
>
> ntb=1,
> ntp=0,
> ntxo=2,
> ioutfm = 1,
> iwrap = 1,
>
> ntwe = 1000,
> ntwx = 10000,
> ntpr = 10000,
> ntwr = 20000,
>
> icfe = 1,
> ifsc = 1,
> clambda = 0.0,
> scalpha = 0.5,
> scbeta = 12.0,
> logdvdl = 0,
>
> ifmbar = 1,
> bar_intervall = 100,
> bar_l_min = 0.0,
> bar_l_max = 1.0,
> bar_l_incr = 0.1,
>
> timask1='',
> timask2=':1',
> scmask1='',
> scmask2=':1',
> /
> &ewald
> /
> ================
>
> #!/bin/bash
>
> #PBS -l nodes=1:ppn=8
> #PBS -N lam0
>
> cd $PBS_O_WORKDIR
>
> pmemd=/opt/amber14/bin/pmemd.MPI
> mpirun="/usr/local/bin/mpirun -np 8"
>
> prmtop="../../prmtop/SKF_avo_solv.prmtop"
> inpcrd="../../prmtop/SKF_avo_solv.inpcrd"
>
>
> echo "Minimising..."
> $mpirun $pmemd -i min0.in -p $prmtop -c $inpcrd -O -o min0.out -e
> min0.en -inf min0.info -r min0.rst -l min0.log
>
> echo "Heating constant pressure 0.5 ns ..."
> $mpirun $pmemd -i heat0.in -p $prmtop -c min0.rst -O -o heat0.out -e
> heat0.en -inf heat0.info -r heat0.rst -x heat0.nc -l heat0.log
>
> echo "Equil constant volume 0.5 ns..."
> $mpirun $pmemd -i equil0.in -p $prmtop -c heat0.rst -O -o equil0.out
> -e equil0.en -inf equil0.info -r equil0.rst -x equil0.nc -l equil0.log
>
> echo "Production constant volume 5 ns..."
> $mpirun $pmemd -i prod0.in -p $prmtop -c equil0.rst -O -o prod0.out
> -e prod0.en -inf prod0.info -r prod0.rst -x prod0.nc -l prod0.log
>
>
> =====================
>
>
> Dr. Fabian Glaser
> Head of the Structural Bioinformatics section
>
> Bioinformatics Knowledge Unit - BKU
> The Lorry I. Lokey Interdisciplinary Center for Life Sciences and
> Engineering Technion - Israel Institute of Technology, Haifa 32000,
> ISRAEL
>
> fglaser at technion dot ac dot il
> Tel: +972 4 8293701
> http://bku.technion.ac.il
>
>
> > On 3 Dec 2015, at 1:46 PM, Hannes Loeffler
> > <Hannes.Loeffler.stfc.ac.uk> wrote:
> >
> > On Wed, 2 Dec 2015 16:49:20 +0200
> > Fabian gmail <fabian.glaser.gmail.com
> > <mailto:fabian.glaser.gmail.com>> wrote:
> >
> >> Dear Jason,
> >>
> >> Thanks a lot for the detailed explanation, so in the paper I am
> >> following, they minimised once, and then use the same minimised
> >> molecule (in this case a drug Improving the Efficiency of Free
> >> Energy Calculations in the Amber Molecular Dynamics Package Joseph
> >> W. Kaus,*,† Levi T. Pierce,†,‡ Ross C. Walker,†,‡ and J. Andrew
> >> cCammon†,§,∥,⊥) and then perform IT calculation.
> >>
> >> So the protocol I will be using will something like this:
> >>
> >> 0) minimize 20,000 steps ONCE for lambda = 0. 5 and then for each
> >> lambda for lambda 0, 0,1, 0,2, etc. … 1 do the following:
> >
> > That's just overkill. All you want to do with the minimisation
> > step is to eliminate close contacts which could lead to high forces
> > in MD and thus potential instability. A few hundred steps should be
> > enough.
> >
> > You could do the same protocol for every lambda as it costs very
> > little. Keep in mind that you are extrapolating away from the
> > references state and thus "equilibration" of later lambdas
> > potentially needs more time. Not sure what the convergence
> > behaviour for decoupling is but at the end-point your molecule will
> > have created a void and will be filled with water.
> >
> > You could take over the data from the previous neighbour lambda if
> > you wanted to but in practice this means that you would have to run
> > all lambdas one after another.
> >
> >
> >> 1) heating to 300 K during 500 ps at constant pressure
> >> 2) equilibrate at constant volume during 500 ps
> >> 3) production for 5 ns at constant volume x 3
> >>
> >> This protocol is different than in the A9 tutorial, but sounds more
> >> right for my purpose, does it sounds a right protocol? This is to
> >> calculate the DG of solvation (disspereance of a drug in water).
> >
> > With constant volume you will get the Helmholtz free energy, not
> > Gibbs.
> >
> >
> > Cheers,
> > Hannes.
> >
> >
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org <mailto:AMBER.ambermd.org>
> > http://lists.ambermd.org/mailman/listinfo/amber
> > <http://lists.ambermd.org/mailman/listinfo/amber>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Dec 04 2015 - 02:00:03 PST
