QM/MM PBSA is a method to calculate the total free energy of
the protein, using only simulations of the end states of the
reactions. It can be seen as a postprocessing of QM/MM minimisation, in
order to get more stable and reliable energies.
The method is originally described in
[1]: M.Kaukonen, P. Söderhjelm, J. Heimdal &
U. Ryde, "A QM/MM-PBSA method for estimates of free energies in
proteins", J. Phys. Chem. B, 2008, 112,
12537-12548;
DOI: 10.1021/jp802648k;
free
reprints.
Please, note that similar approaches have been suggested earliear, but
primarily for ligand-binding studies, whereas our approach is directed
towards studies of enzyme active sites, for which the QM system is
covalently bound to the MM system.
In our
implementation, the QM part is calculated by Turbomole, and MM and
solvation part by Amber (mmpbsa and pbsa).
In the following, we will discuss three systems:
There are several ways to do this partitioning. Two of these are described below (although the second is most developed):
Version 1: 1-2
electrostatics by QM, H-junctions in MM1
In this variant, the QM system is polarised by point charges from
the protein only, the 1-2 electrostatics is calculated in the QM
calculation, and there are H junctions atoms in the MM1
calculation.
The free energy expression is:
Atot=
EQM+ptch + EMM12 - EMM1
+ dGsol+ dGnp
- TdS
11 all 11 noel 11 noel 13
el 13 other vibrational entropy
12 el 12
noel
23
el 23 other
22
all
33
all
system1=QM, system2=explicit MM
(typically protein), system3=implicit
solvent
all=all interactions, el=electrosatics,
noel=bonded+vdw, other=hydrophobic interactions
The
first and last terms (EQM+ptch and - TdS) are
calculated by Turbomole, the others by Amber (sander and mm_pbsa.pl).
Variant
2: 1-2
electrostatics by MM, C-junctions in MM
In this variant, the 1-2 electrostatics are calculated by the MM
program and there are no H junction atoms in the MM1
calculation.
Atot= EQM,1SCF + EMM12
-
EMM1 + dGsol+ dGnp
- TdS
11 all
11
all 11 all 13 el 13
other vibrational entropy
12
all
23
el 23 other
22
all
33
all
system1=QM, system2=explicit MM (typically
protein), system3=implicit solvent
all=all
interactions, el=electrosatics, noel=bonded+vdw,
other=hydrophobic interactions
The
first and last terms (EQM+ptch and - TdS) are
calculated by Turbomole, the others by Amber (sander and mm_pbsa.pl).
The dGsol, dGnp, and TdS terms are identical between the two versions.
There are two additional variants of this approach, depending on how
the QM energy and the charges of the QM system are calculated
| Term |
Calculation |
System |
junctions |
Results
depend on |
|||
| 1
coord |
1
charges |
2
coord |
2
charges |
||||
| Variant 2 in vacuum | |||||||
| E_QM |
QM in vacuum |
1 |
H |
x |
|||
| E_MM1 |
MM1 with charges | 1 |
H |
x |
x |
||
| E_MM12 |
MM12 with charges | 12 |
C |
x |
x |
x |
x |
| dG_sol |
PB or GB |
12 |
C |
x |
x |
x |
x |
| dG_np |
Surface area |
12 |
C |
x |
x |
||
| Variant 2 polarised |
|||||||
| E_QM | QM, polarised with point charges |
1+2 |
H |
x |
x |
x |
|
| E_MM1 | MM1 with charges |
1 |
H |
x |
x |
||
| E_MM12 | MM12 with charges | 12 |
C |
x |
x |
x |
x |
| dG_sol | PB or GB | 12 |
C |
x |
x |
x |
x |
| dG_np | Surface area | 12 |
C |
x |
x |
||
Run standard ComQum calculation.
Follow the instructions
in the ComQum page.
If you intend to use MD snapshots rather than a single QM/MM structure, follow the instructions for QTCP calculations until the point where you have run MD trajectories for all states of interest (typically only the reactants and products, because no intermediate states are needed).
Run standard ComQum calculation.
Follow the instructions
in the ComQum page.
Ensure that you have updated ESP charges for the QM system (especially
if you run protein fixed).
Prepare input files, here version 1
and
all waters stripped, using comqum relaxed structure.
If
you have done only comqum fixed, make sure that you have correct
charges in your prmtop3 and prmtop2 files for the
QM system (and other metal centers as well).
A script to perform this item 3 is:
prepareVersion1.scr (with no water,
version 1)
prepareVersion2_MD.scr(with
MD snapshots, version 2)
Generate a 'basis' pdb file with the correct QM structure and the
correct number of waters in it (depending on whether it is 1) static/MD
with no crystal water 2) static/MD with crystal water )
Change
the indexes of waters in comqum*.dat
if needed
cd COMQUM_directory
Strip waters by changeparm
prmtop3
Truncate
b
/away/bio/Mka/removewaters.belly
n
n
nowat.prmtop
y
p
prmtop3.prmcrd3.pdb
nowat.prmcrd
ambpdb -pqr -p nowat.prmtop < nowat.prmcrd
> nowat.pdbThe file '
removewaters.belly' is a standard
sander input file which contains the atoms to be kept as free atoms in
the belly group
For instance the following keeps atoms 1-9052 fixed in position:
.....
&end
Belly
ATOM 1 9052
END
END
In case you want to keep crystal waters just change the last atom index
accordingly in the belly group
For instance
.....
&end
Belly
ATOM 1 9940
END
END
and
ambpdb -pqr -p crystwat.prmtop < crystwat.prmcrd >
crystwat.pdbCopy files to the directory you are about to run QM/MM-PBSA
mkdir Nowat_version 1
Turbomole files:
cp alpha beta mos
coord control auxbasis basis pointcharges Nowat_version 1;
cp nowat.pdb Nowat_version1/charges.pdb;
cp comqum.dat Nowat_version1
(charges of file 'charges.pdb' are used to generate charges in the 'pointcharges' file of turbomole,
coordinates for pointcharges are taken from ${MMDIR}/${recogpdb} of the jobscript run_QMPC_snapshots.ser below )
comqum.dat is needed when generating varying pointcharges (to avoid including QM in pointcharges)
Amber files:
cp nowat.prmtop Nowat_version1/prmtop3; cp nowat.prmcrd Nowat_version1/prmcrd3;
cp prmtop1 prmcrd1 Nowat_version1
cp /away/bio/Mka/sander.in1 Nowat_version1
(varying coordinates are taken from ${MMDIR}/${recogcrd} of the jobscript run_QMPC_snapshots.ser below )
Amber-PBSA files:
cp /away/bio/Mka/mm_pbsa.in Nowat_version1
(varying coordinates are taken from ${MMDIR}/${recogcrd} of the jobscript run_QMPC_snapshots.ser below )
cd Nowat_version1
For version1 generate the prmtop2 file with zero charges (this should be the case already)
changeparm
prmtop3
1
w
prmtop2
q
For version2: modify the link atoms H->C
changeparm
prmtop3
tr
comqum.dat
n
n
prmtop1
y
mdrest
prmcrd3
prmcrd1
q
# for version2 no charges are zeroed in prmtop2 so prmtop2==prmtop3
rm -f prmtop2; cp prmtop3 prmtop2
For method 1/2 run the script: run_QM_snapshots.ser
You must modify the script based on whether you use method 1 or method 2, by commenting unnecessary lines after line
"comment THIS LINE in the script, the following exit line, and comment the not-needed method below"
If you run a md snapshots run go to QTCP directory where you have trajectories and run
ptraj prmtop /away/bio/Mka/ptraj_snapshots_pdb_noH2O.in
ptraj prmtop /away/bio/Mka/ptraj_snapshots_crd_noH2O.in
To generate coordinate snapshots in pdb and amber format. Copy these files to the directory you submit the job
Change in the file run_QM_snapshots.ser
tags to
recogpdb=mdcrd5_wrap_noH2O.pdb*
recogcrd=mdcrd5_wrap_noH2O.crd*
In the script (run_QMPC_snapshots.ser or run_QMVP_snapshots.ser) the MMDIR is the directory where the QTCP MD run was carried out and you issued the ptraj commands above.
Charges for turbomole are taken from file charges.pdb,
positions from (varying) pdb files
Charges for sander are taken from prmtop1 (EMM1)
andprmtop2 (EMM12 ), Charges for PBSA are taken from
prmtop3.
Before submitting the job check that you have correct info in the job script "run_QM_snapshots.ser" between lines
"#################CHANGE BEGIN###################"
and
"#################CHANGE END ###################"
qsub run_QM_snapshots.ser
Run turbomole program aoforce on the QM system in vacuum
there after run turbomole program freeh to get the entropy.
To analyse result:
avesum_and_std_x-y.pyt indir1
indir2 skip TAG1 TAG2
Note that Amber runs are without explicit water (except crystal
water) and they are carried out with no
peridiocity.
Run a standard ComQum calculation.
Follow the instructions
in the ComQum page.
If you have run a protein fixed, you need to calculate charges for the QM system:
changeparm <<EOF
prmtop3
t
c
n
n
prmtop1.cjunct
y
m
prmcrd3
prmcrd1.cjunct
q
EOF
Run the QM calculation:
Change in the control file the following three keywords:
#$esp_fit
(commented out)
$scfiterlimit 1 (instead of 300)
gridsize 3 (instead of
m3; under $dft)
Then run
ridft>logd&
The energy is only written in the log file, not in the energy file:
grep 'total energy' logd
This
calculation is not necessary, because the same energy terms are also
obtained in the PB or GB calculations below. Skip to next point.
Run the first (MM12) calculation:
cat >sander.in <<EOF
sander file for QM/MM-PBSA, UR 5/12-08
&cntrl
irest=0,ntx=1,
nstlim=1,dt=0.0,
imin=0,maxcyc=1,drms=0.001,
temp0=300.0,ntt=1,tautp=0.2,
ntc=1,ntf=1,
nsnb=25,cut=1000.0,dielc=1.0,
ntpr=1,ntwx=0,ntwv=0,ntwe=1,
ntb=0,ntp=0,taup=0.2,
scnb=2.0,scee=1.2,
ibelly=0
&end
EOF
sander -O -i sander.in -o sander.out2 -p prmtop3.nowat -c prmcrd3.nowat
The energy is the EPtot in sander.out2 (first three)
grep EPtot sander.out2
The same energy can be obtained with changeparm, command e.
Calculate the second MM1 energy:
sander -O -i sander.in -o sander.out1 -p prmtop1.cjunct -c
prmcrd1.cjunct
The energy is the EPtot in sander.out1 (first three)
grep EPtot sander.out1
Calculate the polar solvation
energy with PB
cat >sander.in3 <<EOF
pbsa file for
QM/MM-PBSA PB calculation; from MM/PBSA, UR 5/12-08
&cntrl
ntf
=
1, ntb = 0,
igb
=
10, dielc = 1.0,
cut
=
999.0, nsnb = 99999,
scnb =
2.0, scee = 1.2,
imin =
1, maxcyc =
0, ntmin = 2,
&end
&pb
epsin =
1.0, epsout = 80.0,
istrng =
0, radiopt = 1,
sprob =
1.6, space = 0.5,
maxitn = 500
npbverb= 1,
cutres =
12,fillratio=3
&end
EOF
sander -O -i sander.in3 -o sander.out3 -p prmtop3.nowat -c prmcrd3.nowat
grep EPB sander.out3
Alternatively, you may run a GB calculation instead (it is faster and
often more stable):
cat >sander.in4 <<EOF
pbsa file for
QM/MM-PBSA PB calculation; from MM/PBSA, UR 5/12-08
&cntrl
ntf
=
1, ntb =
0, dielc = 1.0,
idecomp= 0,
igb
= 5, saltcon= 0.00,
offset = 0.09,
intdiel=
1.0, extdiel= 80.0,
gbsa =
0, surften= 1.0,
cut
=
999.0, nsnb = 99999,
scnb =
2.0, scee = 1.2
imin =
1, maxcyc =
1, ncyc = 0,
&end
EOF
sander -O -i sander.in4 -o sander.out4 -p prmtop3.nowat -c prmcrd3.nowat
grep EGB sander.out4
The entropy of the QM system is obtained by the Turbomole
program aoforce.
For this, you first need to optimise the QM system in vacuum.
The programs describe or freeh is the used to get the entropy
------------------------------------------------------------------------------------------------------------------------------------------------
Just a dummy sander input file containing the atoms to be kept in the belly group (other atoms will be deleted)
&cntrl
irest=0,ntx=1,
nstlim=1,dt=0.0,
imin=0,maxcyc=1,drms=0.001,
temp0=300.0,ntt=1,tautp=0.2,
ntc=1,ntf=1,
nsnb=25,cut=1000.0,dielc=1.0,
ntpr=1,ntwx=0,ntwv=0,ntwe=1,
ntb=0,ntp=0,taup=0.2,
scnb=2.0,scee=1.2,
iwforc=1,ibelly=1
&end
Belly
ATOM 1 9052
END
END
------------------------------------------------------------------
Title
&cntrl
irest=0,ntx=1,
nstlim=1,dt=0.0,
imin=0,maxcyc=1,drms=0.001,
temp0=300.0,ntt=1,tautp=0.2,
ntc=1,ntf=1,
nsnb=25,cut=1000.0,dielc=1.0,
ntpr=1,ntwx=0,ntwv=0,ntwe=1,
ntb=0,ntp=0,taup=0.2,
scnb=2.0,scee=1.2,
ibelly=0
&end
------------------------------------------------------------------
#!/bin/bash
# prepare files for a no-water QM/MM-PBSA run
# Method 1 (PC)
gunzip -f *
# truncate the prmtop file
# write it and also the truncated coordinates in amber format.
changeparm <<LOPPU1
prmtop3
Truncate
b
/away/bio/Mka/removewaters.belly
n
n
nowat.prmtop
y
p
prmtop3.prmcrd3.pdb
nowat.prmcrd
q
LOPPU1
#output the trncated system in pdb format
ambpdb -pqr -p nowat.prmtop < nowat.prmcrd > nowat.pdb
mkdir Nowat_version1
#copy Turbomole files:
cp alpha beta mos coord control auxbasis basis pointcharges \
Nowat_version1
cp nowat.pdb Nowat_version1/charges.pdb
cp comqum.dat Nowat_version1
#(charges of file 'charges.pdb' are used to generate charges
# in the 'pointcharges' file of turbomole,
# coordinates for pointcharges are taken from
# ${MMDIR}/${recogpdb} of the jobscript run_QMPC_snapshots.ser below )
# comqum.dat is needed when generating varying pointcharges
# (to avoid including QM in pointcharges)
#COPY THE Amber files:
cp nowat.prmtop Nowat_version1/prmtop3
cp nowat.prmcrd Nowat_version1/prmcrd3
cp prmtop1 prmcrd1 Nowat_version1
cp /away/bio/Mka/sander.in1 Nowat_version1
#(varying coordinates are taken from
# ${MMDIR}/${recogcrd} of the jobscript run_QMPC_snapshots.ser below )
#Amber-PBSA files:
cp /away/bio/Mka/mm_pbsa.in Nowat_version1
#(varying coordinates are taken from ${MMDIR}/${recogcrd} of the
# jobscript run_QM_snapshots.ser below )
cd Nowat_version1
#For version1_PC generate the prmtop2 file with zero charges:
changeparm <<LOPPU2
prmtop3
1
w
prmtop2
q
LOPPU2
cp /away/bio/Mka/run_QM_snapshots.ser .
#!/bin/bash
# prepare files for a no-water QM/MM-PBSA run
# Method 2 (VP)
# You should be one directory above QTCP one
rm -rf MD_version2
mkdir MD_version2
cd QTCP
gunzip auxbasis.gz basis.gz c0.gz comqum.dat.gz control.gz mdrest4_wrap.1.gz mos
.gz alpha.gz beta.gz prmtop.gz
cp auxbasis basis c0 comqum.dat control mdrest4_wrap.1 mos alpha beta prmtop ../
MD_version2
cd ../MD_version2
mv c0 coord
# truncate the prmtop file
# write it and also the truncated coordinates in amber format.
changeparm <<LOPPU1
prmtop
Truncate
b
/away/bio/Mka/removewaters.belly
n
n
prmtop3
y
m
mdrest4_wrap.1
prmcrd3
q
LOPPU1
#output the trncated system in pdb format
ambpdb -pqr -p prmtop3 < prmcrd3 > prmtop3.prmcrd3.pdb
cp prmtop3.prmcrd3.pdb charges.pdb
#(charges of file 'charges.pdb' are used to generate charges
# in the 'pointcharges' file of turbomole,
# coordinates for pointcharges are taken from
# ${MMDIR}/${recogpdb} of the jobscript run_QMVP_snapshots.ser below )
# comqum.dat is needed when generating varying pointcharges
# (to avoid including QM in pointcharges)
#COPY THE Amber files:
cp ../Nowat_version2/prmtop1 .
cp ../Nowat_version2/prmcrd1 .
cp /away/bio/Mka/sander.in1 .
#(varying coordinates are taken from
# ${MMDIR}/${recogcrd} of the jobscript run_QM_snapshots.ser below )
#Amber-PBSA files:
cp /away/bio/Mka/mm_pbsa.in .
#(varying coordinates are taken from ${MMDIR}/${recogcrd} of the
# jobscript run_QM_snapshots.ser below )
#for version2 (VP) modify the link atoms H->C
changeparm <<EOF2
prmtop3
tr
comqum.dat
n
n
prmtop1
y
mdrest
prmcrd3
prmcrd1
q
EOF2
# for version2 no charges are zeroed in prmtop2 so prmtop2==prmtop3
rm -f prmtop2; cp prmtop3 prmtop2
cp /away/bio/Mka/run_QM_snapshots.ser .
#remove parallel stuff from the control file
cp control control.orig
rm -f control.tmp
grep -v '$parallel_platform' control| grep -v '$numprocs' > control.tmp
mv control.tmp control
mkdir ../MD_version2_envOTHER
cp * ../MD_version2_envOTHER
#generate trajectory files
rm ptraj.stripH2O.crd.in
echo "trajin ../QTCP/mdcrd4_wrap.gz 9 200 10" > ptraj.stripH2O.crd.in
echo "trajout noH2Ocrd restart" >> ptraj.stripH2O.crd.in
echo "strip :603-99999" >> ptraj.stripH2O.crd.in
echo "go" >> ptraj.stripH2O.crd.in
gunzip ../QTCP/prmtop.gz
ptraj ../QTCP/prmtop ptraj.stripH2O.crd.in
rm ptraj.stripH2O.pdb.in
echo "trajin ../QTCP/mdcrd4_wrap.gz 9 200 10" > ptraj.stripH2O.pdb.in
echo "trajout noH2Opdb pdb" >> ptraj.stripH2O.pdb.in
echo "strip :603-99999" >> ptraj.stripH2O.pdb.in
echo "go" >> ptraj.stripH2O.pdb.in
gunzip ../QTCP/prmtop.gz
ptraj ../QTCP/prmtop ptraj.stripH2O.pdb.in
------------------------------------------------------------------
trajin mdcrd5.gz 0 200 10
trajout mdcrd5_wrap_noH2O.pdb pdb
image center :WAT byres familiar com :1-1005
strip :1007-1185,:1187-1365,1367-99999
go
------------------------------------------------------------------
trajin mdcrd5.gz 0 200 10
trajout mdcrd5_wrap_noH2O.crd restart
image center :WAT byres familiar com :1-1005
strip :1007-1185,:1187-1365,1367-99999
go
--------------------------------------------------------
RUN SCRIPTS FOR METHOD 1/2
#!/bin/bash
#PBS -l nodes=1
#PBS -l walltime=84:00:00
##PBS -j e
#PBS -N QMMMPBSA_N_HIP_nowat
#PBS -v TURBODIR,AMBERHOME,PATH
. use_modules
module load intel/8.1
#-------- Amber------------------
export AMBERHOME=/sw/pkg/bio/Amber9
export PATH=$PATH:$AMBERHOME/exe
# ------- Ulf's programs ---------
export PATH=/sw/pkg/bio/Bin:$PATH
#-------- Turbomole --------------
export TURBODIR=/sw/turbomole/Turbo5.9
export PATH=$TURBODIR/scripts:$PATH
export PATH=$TURBODIR/bin/`sysname`:$PATH
cd $PBS_O_WORKDIR
#################CHANGE BEGIN###################
# 1) change pattern to recognise crd and pdb snapshot (static) files:
#recogpdb=strip.pdb
#recogcrd=strip.crd
recogpdb=noH2Opdb
recogcrd=noH2Ocrd
#
# 2) MMDIR is the directory from which the MM2 coordinates are taken
# for static ones set MMDIR=$PBS_O_WORKDIR
#MMDIR=${WRKDIR}Hydrogenase_QM-PBSA_QMPC/Comqum_unoohidcys588_PBSAfree_nowater
MMDIR=$PBS_O_WORKDIR
#################CHANGE END ###################
###### SOME CHECKING THET FILES ARE THERE #########
# check that most of the QM files exist
for file in basis control comqum.dat ${MMDIR}/${recogpdb}* charges.pdb ;do
# echo $file
if !(test -e ${file} ); then
echo "SOME FILES FOR QM are MISSING"
echo "MISSING FILE:"
echo ${file}
exit
fi
done
# check that most of the Amber files exist
for file in sander.in1 prmtop1 prmtop2 ${MMDIR}/prmtop2 prmcrd1 prmcrd3;do
# echo $file
if !(test -e ${file} ); then
echo "SOME FILES FOR AMBER are MISSING"
echo "MISSING FILE:"
echo ${file}
exit
fi
done
# check that most of the PBSA files exist
for file in ${MMDIR}/${recogcrd}* comqum.dat mm_pbsa.in prmtop3 prmcrd3;do
# echo $file
if !(test -e ${file} ); then
echo "SOME FILES FOR PBSA are MISSING"
echo "MISSING FILE:"
echo ${file}
exit
fi
done
##############################################################
# run turbomole energies with point charges
# assumes files:
#
# ${MMDIR}/${recogpdb} for variying coordinates and
# charges.pdb for the correct charges
#
# auxbasis basis comqum.dat control coord mos alpha beta
# (mos or [alpha beta], auxbasis with ridft
# results in tm.energies.with.varying.pointcharges
##############################################################
echo "MAIN", ${recogpdb} ${MMDIR}
echo "comment THIS LINE in the script, the following exit line, and comment the not-needed method below"
exit
#FOR METHOD 1 SELECT THIS
/away/bio/Mka/run_tm_with_QMPC.scr ${recogpdb} ${MMDIR}
#FOR METHOD 2 SELECT THIS
/away/bio/Mka/run_tm_with_QMVP.scr ${recogpdb} ${MMDIR}
##############################################################
# run Amber energies :
# sander.in1: 11NOEL, 12NOEL, 22ALL (charges in 1 are zeroed)
# sander.in1: 11NOEL
# assume files:
# ${MMDIR}/${recogcrd}
# sander.in1
# prmtop1 ${MMDIR}/prmtop2 prmcrd1 prmcrd3
#
# qm coordinates are taken from prmcrd3
# mm coordinates are taken from ${MMDIR}/${recogcrd}
#
# results in sander.energies.with.varying.pointcharges
# and sander.Eamber11.noel
##############################################################
/away/bio/Mka/run_sander_snapshots.scr ${recogcrd} ${MMDIR}
##############################################################
# run Amber/DELPHI PBSA energies
#
# assume files:
# ${MMDIR}/${recogcrd}
# mm_pbsa.in prmtop3 comqum.dat prmcrd3
#
# qm coordinates are taken from prmcrd3
# mm coordinates are taken from ${MMDIR}/${recogcrd}
#
##############################################################
/away/bio/Mka/run_pbsa_snapshots.scr ${recogcrd} ${MMDIR}
#!/bin/bash
# run turbomole with changing point charges
# the coordinates of point charges are extracted from
# pdb files by script changepointcharges.pyt
#
# It is assumed that the values of the point charges DO NOT CHANGE
# (only the coordinates change)
#
recogpdb=$1
MMDIR=$2
rm tm.energies.with.varying.pointcharges
rm dscf_problem
touch tm.energies.with.varying.pointcharges
echo "runtm:", ${recogpdb} ${MMDIR}
# is it dscf or ridft?
if !(test -e auxbasis ); then
name=dscf
# echo "DSCF CALCULATION !!!"
else
name=ridft
# echo "RIDFT CALCULATION !!!"
fi
for i in $( ls ${MMDIR}/${recogpdb}*); do
echo item: $i
rm -f ./pointcharges
/away/bio/Mka/changepointcharges.pyt charges.pdb $i pointcharges
rm ${name}.out.tmp
${name} > ${name}.out.tmp
grep 'total energy =' ${name}.out.tmp >> tm.energies.with.varying.pointcharges
if (test -e dscf_problem ); then
echo "Problem with turbomole run in directory:"
pwd
mv tm.energies.with.varying.pointcharges tm.energies.with.varying.pointcharges.problem
mv ${name}.out.tmp ${name}.out.tmp.problem
exit
fi
done
#!/bin/bash
# run turbomole FIRST with point charges
# the coordinates of point charges are extracted from
# pdb files by script changepointcharges.pyt
#
# It is assumed that the values of the point charges DO NOT CHANGE
# (only the coordinates change)
#
# Thereafter run one iteration of dscf/ridft to get energy in vacuum with
# point charge polarised wavefunctions
recogpdb=$1
MMDIR=$2
rm tm.energies.with.varying.pointcharges
touch tm.energies.with.varying.pointcharges
rm -f control.org
cp control control.org
# is it dscf or ridft?
if !(test -e auxbasis ); then
name=dscf
# echo "DSCF CALCULATION !!!"
else
name=ridft
# echo "RIDFT CALCULATION !!!"
fi
for i in $( ls ${MMDIR}/${recogpdb}*); do
echo item: $i
rm -f ./pointcharges
/away/bio/Mka/changepointcharges.pyt ./charges.pdb $i pointcharges
cp control control.pointcharges
rm -f ${name}.out.tmp.PC
${name} > ${name}.out.tmp.PC
# copy the point charge wf for the initial guess on the next round
rm -f alpha.pc beta.pc mos.pc
cp alpha alpha.pc
cp beta beta.pc
cp mos mos.pc
rm -f control.tmp
grep -v '$point_charges' control| grep -v '$scfiterlimit' | grep -v '$scfconv'|sed -e 's/m3/3/g' > control.tmp
rm -f control
echo "\$scfiterlimit 1" > control
# echo "\$scfconv -6">> control
cat control.tmp >> control
rm ${name}.out.tmp
${name} > ${name}.out.tmp
grep 'total energy =' ${name}.out.tmp >> tm.energies.with.varying.pointcharges
rm -f control
cp control.org control
# copy the point charge wf for the initial guess on the next round
rm -f alpha beta mos
cp alpha.pc alpha
cp beta.pc beta
cp mos.pc mos
done
mm_pbsa.in
#
# Input parameters for mm_pbsa.pl
# modified by markus.kaukonen@iki.fi by
# instructions in
#www.teokem.lu.se/~ulf/Methods/mm_pbsa.html#mm_pbsa.in_sample_file_delphi
#
# Holger Gohlke
# 08.01.2002
#
################################################################################
@GENERAL
#
# General parameters
# 0: means NO; >0: means YES
#
# mm_pbsa allows to calculate (absolute) free energies for one molecular
# species or a free energy difference according to:
#
# Receptor + Ligand = Complex,
# DeltaG = G(Complex) - G(Receptor) - G(Ligand).
#
# PREFIX - To the prefix, "{_com, _rec, _lig}.crd.Number" is added during
# generation of snapshots as well as during mm_pbsa calculations.
# PATH - Specifies the location where to store or get snapshots.
#
# COMPLEX - Set to 1 if free energy difference is calculated.
# RECEPTOR - Set to 1 if either (absolute) free energy or free energy
# difference are calculated.
# LIGAND - Set to 1 if free energy difference is calculated.
#
# COMPT - parmtop file for the complex (not necessary for option GC).
# RECPT - parmtop file for the receptor (not necessary for option GC).
# LIGPT - parmtop file for the ligand (not necessary for option GC).
#
# GC - Snapshots are generated from trajectories (see below).
# AS - Residues are mutated during generation of snapshots from trajectories.
# DC - Decompose the free energies into individual contributions
# (only works with MM and GB).
#
# MM - Calculation of gas phase energies using sander.
# GB - Calculation of desolvation free energies using the GB models in sander
# (see below).
# PB - Calculation of desolvation free energies using delphi (see below).
# MS - Calculation of nonpolar contributions to desolvation using molsurf
# (see below).
# If MS == 0, nonpolar contributions are calculated with the LCPO method
# in sander.
# NM - Calculation of entropies with nmode.
#
PREFIX pbsa
PATH ./
#
COMPLEX 1
RECEPTOR 0
LIGAND 0
#
COMPT prmtop3
RECPT
LIGPT
#
GC 0
AS 0
DC 0
#
MM 1
GB 1
PB 1
MS 1
#
NM 0
#
#
################################################################################
@PB
#
# PB parameters (this section is only relevant if PB = 1 above)
#
# The following parameters are passed to the PB solver.
# Additional parameters (e.g. SALT) may be added here.
# For further details see the delphi and pbsa documentation.
#
# PROC - Determines which method is used for solving the PB equation:
# If PROC = 1, the delphi program is applied. If PROC = 2,
# the pbsa program of the AMBER suite is used.
# REFE - Determines which reference state is taken for PB calc:
# If REFE = 0, reaction field energy is calculated with EXDI/INDI.
# Here, INDI must agree with DIELC from MM part.
# If REFE > 0 && INDI > 1.0, the difference of total energies for
# combinations EXDI,INDI and 1.0,INDI is calculated.
# The electrostatic contribution is NOT taken from sander here.
# INDI - Dielectric constant for the molecule.
# EXDI - Dielectric constant for the surrounding solvent.
# SCALE - Lattice spacing in no. of grids per Angstrom.
# LINIT - No. of iterations with linear PB equation.
# PRBRAD - Solvent probe radius in A (e.g. use 1.4 with the PARSE parameter se
t
# and 1.6 with the radii optimized by R. Luo)
#
# Parameters for pbsa only
#
# RADIOPT - Option to set up atomic avity radii for molecular surface calculat
ion
# and dielectric assignment. A value of 0 uses the cavity radii from the prm
top file.
# A value of 1 sets up optimized cavity radii at the pbsa initialization pha
se.
# The latter radii are optimized for model compounds of proteins only; use c
autions
# when applying these radii to nucleic acids.
#
# Parameters for delphi only
#
# FOCUS - If FOCUS > 0, subsequent (multiple) PERFIL and SCALE parameters are
# used for multiple delphi calculations using the focussing technique.
# The # of _focussing_ delphi calculations thus equals the value of FOCUS.
# PERFIL - Percentage of the lattice that the largest linear dimension of the
# molecule will fill.
# CHARGE - Name of the charge file.
# SIZE - Name of the size (radii) file.
#
# SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
# the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
#
#
PROC 2
REFE 0
INDI 1.0
EXDI 80.0
SCALE 2
LINIT 1000
PRBRAD 1.4
ISTRNG 0.0
RADIOPT 0
NPOPT 1
CAVITY_SURFTEN 0.0072
CAVITY_OFFSET 0.00
#
SURFTEN 0.0072
SURFOFF 0.00
#
################################################################################
@MM
#
# MM parameters (this section is only relevant if MM = 1 above)
#
# The following parameters are passed to sander.
# For further details see the sander documentation.
#
# DIELC - Dielectricity constant for electrostatic interactions.
# Note: This is not related to GB calculations.
#
DIELC 1.0
#
################################################################################
@GB
#
# GB parameters (this section is only relevant if GB = 1 above)
#
# The first group of the following parameters are passed to sander.
# For further details see the sander documentation.
#
# IGB - Switches between Tsui's GB (1), Onufriev's GB (2, 5).
# GBSA - Switches between LCPO (1) and ICOSA (2) method for SASA calc.
# Decomposition only works with ICOSA.
# SALTCON - Concentration (in M) of 1-1 mobile counterions in solution.
# EXTDIEL - Dielectricity constant for the solvent.
# INTDIEL - Dielectricity constant for the solute
#
# SURFTEN / SURFOFF - Values used to compute the nonpolar contribution Gnp to
# the desolvation according to Gnp = SURFTEN * SASA + SURFOFF.
#
IGB 2
GBSA 2
SALTCON 0.00
EXTDIEL 80.0
INTDIEL 1.0
#
SURFTEN 0.0072
SURFOFF 0.00
#
################################################################################
@MS
#
# Molsurf parameters (this section is only relevant if MS = 1 above)
#
# PROBE - Radius of the probe sphere used to calculate the SAS.
# Since Bondi radii are already augmented by 1.4A, PROBE should be 0.0
#
PROBE 0.0
#
################################################################################
@PROGRAMS
#
# Program executables
#
avesum_and_std_x-y.pyt
#!/usr/bin/env python
# get average values of X-Y, where X and Y are the same quantity from
# different calculations.
# Get also std of X-Y
# Print results
DATE='18.10.2006'
VERSION='v0.00'
import sys, string, math, glob
split=string.split
def header():
print " "+VERSION+" by Markus Kaukonen, "+DATE+"\n"
print " "+"Read data x from file1, data y from file 2 \n"
print " "+"calculated ave(x-y) and std(x-y)\n"
def usage():
print "Usage: "
print sys.argv[0]+" indir1 indir2 skip TAG1 TAG2"
print "data is in one column in "
print "indir1/PBSA_and_PTCH indir1/QM_and_PTCH indir1/SANDER_and_PTCH"
print "indir2/PBSA_and_PTCH indir2/QM_and_PTCH indir2/SANDER_and_PTCH"
# print "(2625.5002 H->kJ/mol, 4.184 kcal/mol -> kJ/mol)"
# print "tag: tag for lines containing the energy input"
# print "column: which column in data is the energy one, 0,1,2,..."
print "skip: skip z values from beginning, 0,1,2,..."
def ave(x):
result=0
#print len(x),len(y)
for i in range(len(x)):
result = result + x[i]
if (len(x)>0):
result = result / len(x)
return result
def dif(x,y,skip=0,name=""):
result=[]
#print len(x),len(y)
if (len(x)<>len(y)):
print "lens differ", name,len(x),len(y)
sys.exit()
for i in range(len(x)):
#print i
if (i >= skip):
#print i
result.append(y[i]-x[i])
# etot=etot+(ene1[j]-ene2[j])
# print ene1[j]-ene2[j]
return result
def std(ave,x):
result=0
if (len(x)>1):
for i in range(len(x)):
result=result+(x[i]-ave)*(x[i]-ave)
result=result/float(len(x)-1)
result=math.sqrt(result)
else:
result=0
return result
args=sys.argv
if (len(args)<4):
print "your number of arguments is", len(args)
sys.stderr.write("Error: Wrong number of arguments \n,len(args) ")
header()
usage()
sys.exit(1)
if ((len(args) == 4) or (len(args) == 5)):
tag1=" "
tag2=" "
else:
tag1=args[4]
tag2=args[5]
#infile2=open(args[2],'r')
# read data from dir1, QM
name1=glob.glob(args[1]+'/tm.energies.with.varying.pointcharges')
infile1=open(name1[0],'r')
if (len(name1)> 1):
print "too many files match", name1
sys.exit()
line1=infile1.readline()
col1=split(line1)
i=0
ene1qm=[]
while line1:
if string.find(line1,"energy")<>-1:
col1=split(line1)
# print line1
#print col1[col], col1[4],i
i=i+1
# print i, imin, imax
ene1qm.append(float(col1[4])*2625.5002)
line1=infile1.readline()
infile1.close()
#read data from dir1, MM12
name1=glob.glob(args[1]+'/sander.energies.with.varying.pointcharges')
infile1=open(name1[0],'r')
if (len(name1)> 1):
print "too many files match", name1
sys.exit()
line1=infile1.readline()
col1=split(line1)
i=0
ene1mm12=[]
while line1:
if string.find(line1,"Etot")<>-1:
col1=split(line1)
# print line1
# print col1[col], col1[2]
i=i+1
# print i, imin, imax
ene1mm12.append(float(col1[2])*4.184)
line1=infile1.readline()
infile1.close()
# read MM1 data
name1=glob.glob(args[1]+'/sander.Eamber11.noel')
infile1=open(name1[0],'r')
if (len(name1)> 1):
print "too many files match", name1
sys.exit()
line1=infile1.readline()
col1=split(line1)
ene1mm1 = (float(col1[2])*4.184)
infile1.close()
#read pbsa data from the first set
name1=glob.glob(args[1]+'/pbcal.energies.with.varying.pointcharges')
infile1=open(name1[0],'r')
if (len(name1)> 1):
print "too many files match", name1
sys.exit()
line1=infile1.readline()
col1=split(line1)
i=0
ene1pbcal=[]
while line1:
if string.find(line1,"PBCAL")<>-1:
col1=split(line1)
# print line1
# print col1[col], col1[1]
i=i+1
# print i, imin, imax
ene1pbcal.append(float(col1[1])*4.184)
line1=infile1.readline()
infile1.close()
# read data from dir2, QM
name2=glob.glob(args[2]+'/tm.energies.with.varying.pointcharges')
infile2=open(name2[0],'r')
if (len(name2)> 1):
print "too many files match", name2
sys.exit()
line2=infile2.readline()
col2=split(line2)
i=0
ene2qm=[]
while line2:
if string.find(line2,"energy")<>-1:
col2=split(line2)
# print line2
# print col2[col], col2[4]
i=i+1
# print i, imin, imax
ene2qm.append(float(col2[4])*2625.5002)
line2=infile2.readline()
infile2.close()
# read data from dir2, MM12
name2=glob.glob(args[2]+'/sander.energies.with.varying.pointcharges')
infile2=open(name2[0],'r')
if (len(name2)> 1):
print "too many files match", name2
sys.exit()
line2=infile2.readline()
col2=split(line2)
i=0
ene2mm12=[]
#print infile2
while line2:
if string.find(line2,"Etot")<>-1:
col2=split(line2)
# print line2
# print col2[col], col2[2]
i=i+1
# print i, imin, imax
ene2mm12.append(float(col2[2])*4.184)
line2=infile2.readline()
infile2.close()
# read MM1 data
name2=glob.glob(args[2]+'/sander.Eamber11.noel')
infile2=open(name2[0],'r')
if (len(name2)> 1):
print "too many files match", name2
sys.exit()
line2=infile2.readline()
col2=split(line2)
ene2mm1 = (float(col2[2])*4.184)
infile2.close()
#read pbsa data from the first set
name2=glob.glob(args[2]+'/pbcal.energies.with.varying.pointcharges')
infile2=open(name2[0],'r')
if (len(name2)> 1):
print "too many files match", name2
sys.exit()
line2=infile2.readline()
col2=split(line2)
i=0
ene2pbcal=[]
while line2:
if string.find(line2,"PBCAL")<>-1:
col2=split(line2)
# print line2
# print col2[col], col2[1]
i=i+1
# print i, imin, imax
ene2pbcal.append(float(col2[1])*4.184)
line2=infile2.readline()
infile2.close()
# calculate averages and STD:s
skip=int(args[3])
eqm_dif=dif(ene1qm,ene2qm,skip,"QM")
eqm_ave=ave(eqm_dif)
eqm_std=std(eqm_ave,eqm_dif)
emm12_dif=dif(ene1mm12,ene2mm12,skip,"MM12")
emm12_ave=ave(emm12_dif)
emm12_std=std(emm12_ave,emm12_dif)
epbcal_dif=dif(ene1pbcal,ene2pbcal,skip,"PBCAL")
epbcal_ave=ave(epbcal_dif)
epbcal_std=std(epbcal_ave,epbcal_dif)
if (len(eqm_dif) <> len(emm12_dif)):
print "lens differ QM: MM12:",len(eqm_dif),len(emm12_dif)
if (len(eqm_dif) <> len(epbcal_dif)):
print "lens differ QM: PBCAL:",len(eqm_dif),len(epbcal_dif)
etot_dif=[]
for i in range(len(eqm_dif)):
etot_dif.append(eqm_dif[i]+emm12_dif[i]-ene2mm1+ene1mm1+epbcal_dif[i])
etot_ave=ave(etot_dif)
etot_std=std(etot_ave,etot_dif)
print "data contains",len(ene1qm),len(ene1mm12),len(ene1pbcal),\
"skipped first",args[3]
#print tag2+"-"+tag1,"Init:dir1 ",args[1]," Final:dir2 ",args[2]," ,DIR2-DIR1"
#print tag2+"-"+tag1,args[2]+"-"+args[1]
print tag2+"-"+tag1," skip "+ args[2]+"-"+args[1],\
"Init:dir1 ",args[1]," Final:dir2 ",args[2]
print '%-12s%-12s%-12s%-12s%-12s%-12s%-12s%-12s'\
%("dQM","std(dQM)","dMM12-dMM1","std(dMM12)","dPBCAL",\
"std(dPBCAL)","dEtot","std(dEtot)")
#print eqm_ave, eqm_std, emm12_ave, emm12_std, \
# ene2mm1-ene1mm1, epbcal_ave, epbcal_std
print '%-12.2f%-12.2f%-12.2f%-12.2f%-12.2f%-12.2f%-12.2f%-12.2f'\
%(eqm_ave, eqm_std, emm12_ave-ene2mm1+ene1mm1, emm12_std,\
epbcal_ave, epbcal_std,etot_ave,etot_std)
print "=========================================================================
======================\n"