Dear AMBER users,
I was hoping for another set of eyes to gaze upon my input files. My simulation works and I just want to check whether I am missing something and/or over/under-complicating things as I am the only AMBER user in my department and thus no-one can point out any flaws/mistakes in what I have done.
I am using AMBER to minimise and equilibrate a homology model of a Cytochrome P450 enzyme in explicit solvent. I used parameters in Shahrokh et al. 2012 for the heme and Cys-Fe bond. The overall aim is to carry out docking studies to identify residues important for docking and/or catalysis.
My leap file input is:
source leaprc.ff14SB
source leaprc.gaff
loadamberparams frcmod.ionsjc_tip3p
loadamberparams frcmod.hem
HEM = loadmol2 hem.mol2
CYP = loadmol2 cyp.mol2
2j2 = loadpdb CYP2J2.pdb
set 2j2.460.28 element "Fe"
bond 2j2.460.28 2j2.405.8 #BOND CYS-S to HEM-FE
set 2j2.405 connect0 2j2.405.1
set 2j2.405 connect1 2j2.405.9
bond 2j2.404.9 2j2.405.1 #BOND CYS to the protein (new CYP residue)
bond 2j2.405.9 2j2.406.1 #As above
savemol2 2j2 CYP2J2_FE.mol2 0
saveamberparm 2j2 CYP2J2_FE_gas.prmtop CYP2J2_FE_gas.rst7
addions 2j2 Na+ 0
solvateoct 2j2 TIP3PBOX 10
saveamberparm 2j2 CYP2J2_FE_box.prmtop CYP2J2_FE_box.inpcrd
The system then has 460 residues and 13637 waters. I then minimise while restraining various parts in successive minimisations:
Initial Minimisation with solute fixed
&cntrl
imin=1, irest=0, ntx=1, !minimise, no restart, no vels
ntmin=1, maxcyc=2000, ncyc=500, !steepest followed by conjugated gradient minimisation
ntb=1, !PBCs
ntr=1, restraintmask=':1-460', restraint_wt=100.0, !apply restraint mask
cut=8.0, !set nonbond cutoff at 8.0A
ntxo=1 !set output to ASCII
/
CYP2J2: Minimisation 2 with heavy atoms fixed
&cntrl
imin=1, irest=0, ntx=1, !minimise, no restart, no vels
ntmin=1, maxcyc=2000, ncyc=500,!steepest followed by conjugated gradient minimisation
ntb=1, !PBCs
ntr=1, restraintmask=':1-460 & !.H=', restraint_wt=10.0 !restrain heavy atoms
cut=8.0, !set nonbond cutoff at 8.0A
ntxo=1 !set output to ASCII
/
CYP2J2: Minimisation 3 with backbone fixed
&cntrl
imin=1, irest=0, ntx=1, !minimise, no restart, no vels
ntmin=1, maxcyc=2000, ncyc=500,!steepest followed by conjugated gradient minimisation
ntb=1, !PBCs
ntr=1, restraintmask=':1-459 & .CA,C,O,N,H', restraint_wt=10.0 !restrain backbone
cut=8.0, !set nonbond cutoff at 8.0A
ntxo=1 !set output to ASCII
/
CYP2J2: Minimisation with no restraints
&cntrl
imin=1, irest=0, ntx=1, !minimise, no restart, no vels
ntmin=1, maxcyc=10000, ncyc=3000,!steepest followed by conjugated gradient minimisation
ntb=1, !PBCs
ntr=0, !no restraints
cut=8.0, !set nonbond cutoff at 8.0A
ntxo=1 !set output to ASCII
/
I then heat:
CYP2J2: Heat while restraining solute
&cntrl
imin=0, irest=0, ntx=1, !No minimise, no restart, no vels
nstlim=100000, dt=0.002, !100,000 steps of 0.002ps
ntf=2, ntc=2, !Set SHAKE to constrain H
tempi=0.0, temp0=300.0, !Set initial and final temps
ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
ntb=1, ntp=0, !PBCs, no pressure scaling,
ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random seed
nmropt=1, !Use varying conditions
cut=8.0, !Set nonbond cutoff to 8.0A
ntr=1, restraintmask=':1-460', restraint_wt=10.0, !Apply restraint to solute
ntxo=2 !Output in netCDF format
/
&wt type='TEMP0', istep1=0, istep2=90000, value1=0.0, value2=300.0 /
&wt type='TEMP0', istep1=90001, istep2=100000, value1=300.0, value2=300.0 /
&wt type='END' /
and equilibrate, again by restraining various bits in successive runs:
CYP2J2: Equilibration phase
&cntrl
imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
nstlim=10000, dt=0.002, !10,000 steps of 0.002ps
ntf=2, ntc=2, !Set SHAKE to constrain H
temp0=300.0, !Set temp
ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
ntwr=10000, !Write restart every 10000 steps
ioutfm=1, !Output mdcrd and mdvel in netCDF format
ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random seed
cut=8.0, !Set nonbond cutoff t 8.0A
ntr=1, restraintmask=':1-460', restraint_wt=10.0, !Apply restraint
ntxo=2 !Output in netCDF format
/
CYP2J2: Equilibration phase
&cntrl
imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
nstlim=10000, dt=0.002, !10,000 steps of 0.002ps
ntf=2, ntc=2, !Set SHAKE to constrain H
temp0=300.0, !Set temp
ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
ntwr=10000, !Write restart every 10000 steps
ioutfm=1, !Output mdcrd and mdvel in netCDF format
ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random seed
cut=8.0, !Set nonbond cutoff t 8.0A
ntr=1, restraintmask=':1-459 & .CA,C,O,N,H', restraint_wt=10.0, !Apply restraint
ntxo=2 !Output in netCDF format
/
CYP2J2: Equilibration phase
&cntrl
imin=0, irest=1, ntx=5, !No minimise, no restart, no vels
nstlim=50000, dt=0.002, !50,000 steps of 0.002ps
ntf=2, ntc=2, !Set SHAKE to constrain H
temp0=300.0, !Set temp
ntpr=500, ntwx=500, !Write energy and coordinates every 500 steps
ntwr=10000, !Write restart every 10000 steps
ioutfm=1, !Output mdcrd and mdvel in netCDF format
ntb=2, ntp=1, pres0=1.0, taup=2.0, !PBCs auto ntb, pressure scaling,
ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random seed
cut=8.0, !Set nonbond cutoff t 8.0A
ntxo=2 !Output in netCDF format
/
And finally the production run:
CYP2J2: Production Phase
&cntrl
imin=0, irest=1, ntx=5, !No minimise, restart, vels
nstlim=10000000, dt=0.002, !10,000,000 steps of 0.002ps
ntf=2, ntc=2, !Set SHAKE to constrain H
temp0=300.0, !Set temp
ntpr=1000, ntwx=1000, !Write energy and coordinates every 1000 steps
ntwr=100000, !Write restart every 100000 steps
ioutfm=1, !Output mdcrd and mdvel in netCDF format
ntb=2, ntp=1, pres0=1.0, taup=1.0, !PBCs auto ntb, pressure scaling
ntt=3, gamma_ln=2.0, ig=-1, !Langevin dynamics, collision freq, random seed
cut=8.0, !Set nonbond cutoff to 8.0A
ntxo=2 !Output in netCDF format
/
My submission script then looks something like this:
cd CYP2J2_Sim1
RunName=CYP2J2_FE_box
cd 001.leap
cp -t ../002.min $RunName.prmtop $RunName.rst7
cd ../002.min
$AMBERHOME/bin/pmemd.cuda -O -i min1.in -o min1.out -p $RunName.prmtop -c $RunName.rst7 -r min1.rst7 -ref $RunName.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min2.in -o min2.out -p $RunName.prmtop -c min1.rst7 -r min2.rst7 -ref min1.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min3.in -o min3.out -p $RunName.prmtop -c min2.rst7 -r min3.rst7 -ref min2.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i min4.in -o min4.out -p $RunName.prmtop -c min3.rst7 -r min4.rst7 -inf ../mdinfo
cp -t ../003.heat $RunName.prmtop min4.rst7
cd ../003.heat
$AMBERHOME/bin/pmemd.cuda -O -i heat.in -o heat_res.out -p $RunName.prmtop -c min4.rst7 -r heat_res.rst7 -x heat_res.mdcrd -ref min4.rst7 -inf ../mdinfo
cp -t ../004.equil $RunName.prmtop heat_res.rst7
cd ../004.equil
$AMBERHOME/bin/pmemd.cuda -O -i equil1.in -o equil1.out -p $RunName.prmtop -c heat_res.rst7 -r equil1.rst7 -x equil1.mdcrd -ref heat_res.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i equil2.in -o equil2.out -p $RunName.prmtop -c equil1.rst7 -r equil2.rst7 -x equil2.mdcrd -ref equil1.rst7 -inf ../mdinfo
$AMBERHOME/bin/pmemd.cuda -O -i equil3.in -o equil3.out -p $RunName.prmtop -c equil2.rst7 -r equil3.rst7 -x equil3.mdcrd -inf ../mdinfo
cp -t ../005.md $RunName.prmtop equil3.rst7
cd ../005.md
$AMBERHOME/bin/pmemd.cuda -O -i md.in -o md.out -p $RunName.prmtop -c equil3.rst7 -r md.rst7 -x md.mdcrd -inf ../mdinfo
I just want to point out again that it runs without any discernible errors, and the RMSD plots and energies look reasonable. From my limited experience in simulations though, I found that this does not always imply a healthy simulation. I will appreciate any input and criticisms.
Many thanks,
Kavin
-----------------------------
Kavin Abelak
PhD Student
Comparative Biomedical Sciences
Royal Veterinary College
4 Royal College St
London NW1 0TU
[
http://www.rvc.ac.uk/media/default/images/social-media/change.gif]
[RVC Logo - link to RVC Website]<
http://www.rvc.ac.uk> [Twitter icon - link to RVC (Official) Twitter] <
http://twitter.com/RoyalVetCollege> [Facebook icon - link to RVC (Official) Facebook] <
http://www.facebook.com/theRVC> [YouTube icon - link to RVC YouTube] <
http://www.youtube.com/user/RoyalVetsLondon?feature=mhee> [Pinterest icon - link to RVC Pinterest] <
http://pinterest.com/royalvetcollege/> [Instagram icon - link to RVC Instagram] <
http://instagram.com/royalvetcollege>
This message, together with any attachments, is intended for the stated addressee(s) only and may contain privileged or confidential information. Any views or opinions presented are solely those of the author and do not necessarily represent those of the Royal Veterinary College (RVC). If you are not the intended recipient, please notify the sender and be advised that you have received this message in error and that any use, dissemination, forwarding, printing, or copying is strictly prohibited. Unless stated expressly in this email, this email does not create, form part of, or vary any contractual or unilateral obligation. Email communication cannot be guaranteed to be secure or error free as information could be intercepted, corrupted, amended, lost, destroyed, incomplete or contain viruses. Therefore, we do not accept liability for any such matters or their consequences. Communication with us by email will be taken as acceptance of the risks inherent in doing so.
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Oct 09 2015 - 04:00:03 PDT
