Difference between revisions of "2016 AMBER tutorial with Beta Trypsin"
Stonybrook (talk | contribs) (→PMEMD) |
Stonybrook (talk | contribs) (→PMEMD) |
||
| Line 72: | Line 72: | ||
restraintmask = ':1-224 & !@H=', | restraintmask = ':1-224 & !@H=', | ||
restraint_wt = 0.1, | restraint_wt = 0.1, | ||
| + | |||
| + | 05mi.in: equilibration | ||
| + | &cntrl | ||
| + | imin = 1, maxcyc = 1000, ntmin = 2, | ||
| + | ntx = 1, ntc = 1, ntf = 1, | ||
| + | ntb = 1, ntp = 0, | ||
| + | ntwx = 1000, ntwe = 0, ntpr = 1000, | ||
| + | cut = 8.0, | ||
| + | ntr = 1, | ||
| + | restraintmask = ':1-224 & !@H=', | ||
| + | restraint_wt = 0.05, | ||
==IV. Simulation Analysis== | ==IV. Simulation Analysis== | ||
Revision as of 13:34, 11 April 2016
For additional Rizzo Lab tutorials see AMBER Tutorials.
In this tutorial, we will learn how to run a molecular dynamics simulation of a protein-ligand complex. We will then post-process that simulation by calculating structural fluctuations (with RMSD) and free energies of binding (MM-GBSA).
Contents
I. Introduction
Yaping
II. Structural Preparation
Antechamber, Parmchk, tLeap
Omar, Katie
III. Simulation using pmemd
PMEMD
Agatha, Beilei
MD input file1 01mi.in:
mkdir 01mi.in
01mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,
Copy 01md.in and rename it to 02md.in:
cp 01mi.in 02md.in
02md.in: equilibration &cntrl imin = 0, ntx = 1, irest = 0, nstlim = 50000, temp0 = 298.15, tempi = 298.15, ig = 71287, ntc = 2, ntf = 1, ntt = 1, dt = 0.001, ntb = 2, ntp = 1, tautp = 0.5, taup = 0.5, ntwx = 1000, ntwe = 0, ntwr = 1000, ntpr = 1000, cut = 8.0, iwrap = 1, ntr = 1, nscm = 100, restraintmask = ':1-224 & !@H=', restraint_wt = 5.0,
Using the same way to make 03mi.in, 04mi.in, 05mi.in, 06md.in, 07md.in, 08md.in, 09md.in, 10md.in, 11md.in:
03mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 2.0,
04mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 0.1,
05mi.in: equilibration &cntrl imin = 1, maxcyc = 1000, ntmin = 2, ntx = 1, ntc = 1, ntf = 1, ntb = 1, ntp = 0, ntwx = 1000, ntwe = 0, ntpr = 1000, cut = 8.0, ntr = 1, restraintmask = ':1-224 & !@H=', restraint_wt = 0.05,
IV. Simulation Analysis
Ptraj
Lauren, Haoyue
ptraj.strip.wat.in:
trajin ../003.PMEMD/10md.trj 1 1000 1 trajin ../003.PMEMD/11md.trj 1 1000 1 trajout 1BJU.trj.strip nobox strip :WAT
ptraj.rec.rmsd.in:
trajin 1BJU.trj.strip trajout 1BJU.com.trj.stripfit reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7 rms reference out 1BJU.rmsd.CA.dat :1-223@CA
ptraj.lig.rmsd.in:
trajin 1BJU.com.trj.stripfit reference ../002.ANTE.TLEAP/1BJU.com.gas.leap.rst7 rms reference out 1BJU.lig.rmsd.dat :224@C*,N*,O*,S* nofit
ptraj.strip.rec.in:
trajin 1BJU.com.trj.stripfit 1 2000 1 trajout 1BJU.lig.trj.stripfit strip :1-223
ptraj.strip.lig.in:
trajin 1BJU.com.trj.stripfit 1 2000 1 trajout 1BJU.rec.trj.stripfit strip :224
ptraj.dist.in:
trajin 1BJU.com.trj.stripfit distance 171OD2_224N1 :171@OD2 :224@N1 out 171OD2_224N1.dat distance 171OD1_224N2 :171@OD1 :224@N2 out 171OD1_224N2.dat
MM-GBSA Energy Calculation
begin by creating a new directory
mkdir 005.MMGBSA/
as well as an input file
vim.gbrescore.in
the following parameters will provide a good estimate of the energy of binding based on solvation enthalpy and binding enthalpy
Single point GB energy calc &cntrl ntf = 1, ntb = 0, ntc = 2, idecomp= 0, igb = 5, saltcon= 0.00, gbsa = 2, surften= 1.0, offset = 0.09, extdiel= 78.5, cut = 99999.0, nsnb = 99999, imin = 5, maxcyc = 1, ncyc = 0, /
We can now submit a job to the cluster we are working on. This script, inside a file named run.sander.rescore.csh , is designed for LIRED, however it can be modified for whatever supercomputer you are working on
#! /bin/tcsh #PBS -l nodes=1:ppn=1 #PBS -l walltime=48:00:00 #PBS -o zzz.qsub.out #PBS -e zzz.qsub.err #PBS -V #PBS -N mmgbsa set workdir = /nfs/guest46/amber_tutorial/005.mmgbsa cd $workdir sander -O -i gb.rescore.in \ -o gb.rescore.out.com \ -p ../002.tleap/4TKG.com.gas.leap.prm7 \ -c ../002.tleap/4TKG.com.gas.leap.rst7 \ -y ../004.ptraj/4TKG.com.trj.stripfit \ -r restrt.com \ -ref ../002.tleap/4TKG.com.gas.leap.rst7 \ -x mdcrd.com \ -inf mdinfo.com
sander -O -i gb.rescore.in \ -o gb.rescore.out.lig \ -p ../002.tleap/4TKG.lig.gas.leap.prm7 \ -c ../002.tleap/4TKG.lig.gas.leap.rst7 \ -y ../004.ptraj/4TKG.lig.trj.stripfit \ -r restrt.lig \ -ref ../002.tleap/4TKG.lig.gas.leap.rst7 \ -x mdcrd.lig \ -inf mdinfo.lig sander -O -i gb.rescore.in \ -o gb.rescore.out.test.rec \ -p ../002.tleap/4TKG.rec.gas.leap.prm7 \ -c ../002.tleap/4TKG.rec.gas.leap.rst7 \ -y ../004.ptraj/4TKG.rec.trj.stripfit \ -r restrt.rec \ -ref ../002.tleap/4TKG.rec.gas.leap.rst7 \ -x mdcrd.rec \ -inf mdinfo.rec exit
Then submit the job as a .csh
qsub run.sander.rescore.csh
Your ligand, receptor, and receptor-ligand complex will all have .out files. The format should look
FINAL RESULTS
NSTEP ENERGY RMS GMAX NAME NUMBER
1 5.9132E+03 2.0005E+01 1.2640E+02 C 159
BOND = 661.8980 ANGLE = 1751.7992 DIHED = 2581.7692
VDWAALS = -1696.6585 EEL = -13958.9335 EGB = -3125.9524
1-4 VDW = 747.0185 1-4 EEL = 7750.8118 RESTRAINT = 0.0000
ESURF = 11201.4791
minimization completed, ENE= 0.59132314E+04 RMS= 0.200047E+02
We can now use a system of equations to elucidate various energies related to our system. VDWAALS = ΔGvdw E(ELS) = ΔGcoulombic E(GB) = ΔGpolar SASA = ESURF
ΔGnonpolar can be estimated by:
ΔGnonpolar = (Solvent Accessible Surface Area)*0.00542 + 0.92
This equation is based on empirical observations relating a molecules surface area and solvation energy
ΔGnonpolar and ΔGmmgbsa are related by anther equation:
ΔGmmgbsa = ΔGvdw + ΔGcoul + ΔGpolar + ΔGnonpolar
all the information needed is in our output files.
ΔΔGbind can now be solved via the equation:
ΔΔGbind = ΔG(mmgbsa complex) – (ΔG(mmgbsa ligand) + ΔG(mmgbsa receptor) (3)
ΔΔGbind is a function of ligand position, so it is useful to plot it. For publications, it is standard protocol to include the mean and standard deviation for your ΔΔGbind.
Lets make a script named get.mmgbsa.csh to extract this sophisticated energetic estimation from the three initial output files of the calculations obtained solve for ΔΔGbind:
bash get.mmgbsa.csh
The output file will be .dat and we can easily use the program xmgrace to plot the information visually
xmgrace MMGBSA_vs_time.dat
