QTCP-Lig
QTCP-Lig is a new method to estimate ligand-binding free energies
from QM/MM structures.
Currently, the best seems to run the standard QTCP scripts only for
the QM step and normal Amber FEP calculations for the other steps.
At the start, you should have two QM/MM structures of the two
ligands in their proteins/hosts. One will be considered the big
ligand (Bbb), the other the smaller ligand (Sss).
You need to run QM/MM-PBSA for both ligands.
You should run two sets of separate calculations for each pair of
ligands: one with the protein/host and one without.
Each QM/MM calculation of the smaller ligand must be started from
the QM/MM structure of the bigger ligand, so that they have similar
conformation and that possible water molecules in the QM system have
the same residue number.
- Ensure that you have run qmmmpbsa-ptch for the BBB ligandsand
SSS , and that the results are in the Ptch subdirectories.
- Make a directory where you want to have all the calculations.
It may be anywhere in the file system.
- Construct a leap.in file that works for the BBB and SSS
ligands (i.e. it should contain *.in and *.dat files for the
residue to be deprotonated in the protonated, deprotonate, and
dummy-atom state).
Sample
of leap commants (for more information see the equilibration
page).
The file should be placed in the directory in point 2 and be
called leap.in.
Do not forget to include possible extra bond within the
protein.
The name of the pdb file should be pdb3-cq
It should include a solvation command, like
solvateOct x TIP3PBOX 3
- Construct a file qtcp-input (in the directory in point 2) with
the following lines:
ligbigdir=direcory of the QM/MM calculations of big ligand (BBB;
with a Ptch subdirectory).
ligsldir=direcory of the QM/MM calculations of the smaller
ligand (SSS; with a Ptch subdirectory).
comres=number of
the most central residue.
solvacc=path and name of the solvacc file (if you already have a
solvacc file for thsi protein)
nmiss= the difference in the number of atoms in the two ligands
miss = a list of the atoms deleted from the big ligand in the
small ligand (e.g. miss="191 192 193")
mask should either be ' restraintmask="!:WAT & !@H='
or, if you have water molecules in the QM system, include the
residue number of them, e.g. mask=' restraintmask="!:WAT
& !@H= | :429,447,588,814"'
If you do not want to add TER within the protein, you can also
include
terval=3
- Run the script
qtcp-lig-1 >qtcp-lig-1.out
It will set up and run the first MD simulations for the BBB
state.
It will also set up the prmtop files for the SSS state.
It reads the qtcp-input file, so the job must be started from
the directory in point 2.
The script takes ~1 minute and should give no errors.
After running the script, one job should be running (qbbb1).
Check thoroughly the tleap outputs (twice)
tail -n 1 qtcp-lig-1.out
grep Error qtcp-lig-1.out
grep Added qtcp-lig-1.out
| grep -v Merz
grep usage qtcp-lig-1.out
grep 'The charge
on' qtcp-lig-1.out
grep 'before and'
qtcp-lig-1.out
grep 'More than one residue' qtcp-lig-1.out
"qtcp-lig-1 ended successfully" should be the last line
Twice: Exiting LEaP: Errors = 0
Twice Added: First many residues, second no residue
No usage (problem with added bonds)
Four lines, First line, fewer than the others. Two last, the
same number
Four lines, all with the same net charge.
Normally, the should be no 'More than one residue'.
- After starting this simulation, start also the FEP
calculations:
mkdir Fep
cd Fep
cp ../Bbb/pdb3-cq .
cp ../Bbb/leap.in .
copy the coordinates of the second ligand from ../Sss/pdb3-cq
and insert after the first ligand in the (just copied) pdb3-cq
file. Insert also also TER between the two ligands.
Change the leap.in file so that it reads in the pdb3-cq file and
that the prmtop and crd files are either prot.prm and prot.rst
(with protein) or wat.prm and wat.rst (without protein).
Run
tleap -f leap.in
and check that it was successful.
Set up a fep.in file with:
1.0 # Length of equilibration
5.0 # Length of production simulation
13 # Number of lambda values
GU2 # Residue name of big ligand
GU1 # Residue name of small ligand
# Soft-core maske of
big ligand (blank line = whole ligand)
# Soft-core maske of small
ligand (blank line = whole ligand)
0 # Number of first non-crystal
water molecule
a # computer (a=Aurora, k=kebnekaise)
n # if independent
simulations are run
# max time for prot and wat
simulations
w # w for wat simulation and
p for protein simulation
:1,2,13,17,220,378 # residues in the QM
system
Run
mkfep-lig <fep.in
Start all fep calculations with
. sub.bash
- When the MD simulation is finished (it typically take some
hours; check the file Bbb/sander.out3), run the script
qtcp-lig-2
>qtcp-lig-2.out
It will set up all the other states and run the long MD
simulations.
It reads the qtcp-input file, so the job must be started from
the directory in point 2.
mv: cannot stat `mdrest3-old': No such file or directory is a
innocent warning.
The script takes < 1 minute and should give no errors.
After running the script, 2 should be running (qbbb2 and qsss2).
tail -n 1 qtcp-lig-2.out
grep Error qtcp-lig-2.out
grep 'Final RMS'
qtcp-lig-2.out
grep 'Maximum junction'
qtcp-lig-2.out
grep 'The coordinates of'
qtcp-lig-2.out
"qtcp-lig-2 ended
successfully" should be the last line
No Error
Final RMS should be quite small (~0.1 Å)
No junction atoms should move more than 0.5 Å
The corrdinates of only the QM atoms should move (for the two
ligands).
- When both MD simulations are finished (they typically take ~48
h; check ls -l */mdcrd5-unwrap; two files of same size are
expected), run the script
qtcp-lig-3
>qtcp-lig-3.out
It will set up and run all the two MM->QM QTCP calculations.
It reads the qtcp-input file, so the job must be started from
that directory.
It takes a few minutes (but ~1 hour with the solvacc
calculation).
After running the script, 2 jobs should be running: qbbb3 and
qsss3.
Some versions of ptraj give some silly output for each of the 15
calculations:
Residue labels:
Checking coordinates:
mdcrd5-unwrap
cp: cannot stat `mos': No
such file or directory
(or alpha/beta) twice is also OK.
tail -n 1 qtcp-lig-3.out
grep Error qtcp-lig-3.out
grep octahedron qtcp-lig-3.out | wc
grep 'Max dist'
qtcp-lig-3.out
qtcp-lig-3 ended successfully should be the last line
No Error
0 0 0
The fourth one
should give 2 similar values that are quite low <70 Å.
- Collect the results and tidy up with the script
qtcp-lig-4
It will give you a file result with all results in tab-limited
format
It can also remove unnecessary files (including mdcrd5-unwrap)
and zip al files with qtcp-lig-4 -dozip