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
Received on Thu Dec 03 2015 - 04:30:03 PST
