QM/MM-PBSA

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:

  1. System 1 = the QM system = the active site
  2. System 2 = the MM system = the surrounding protein (possibly with a few importatnt crystallographic water molecules
  3. System 3 = water = the surrounding solvent
The method is based on the MM/PBSA method by Kollman et al.:
[2]: S. Miyamoto & P. A. Kollman, Proteins 16 (1993) 226-245
[3]: P. A. Kollman et al. Acc. Chem. Res. 33 (2000) 889-897.
In this approach, the free energy is calculated from the equation:
dG = E_elstat + E_vdW + E_int + dG_solv + dG_np - TdS, where
In the QM/MM-PBSA approach, we replace E_elstat + E_vdW + E_int by a QM/MM energy (and the entropy term outside the QM system is normally ignored). The QM/MM energy is typically
E_QM/MM = E_QM1 + E_MM12 - E_MM1
All calculations can either be performed on a single QM/MM minimised structure or on a set of structures obtained from a MD simulation. Typically, all calculations are performed with the same MM system but different QM systems (reactant and product states). Thereby, the energies become more stable.

There are several ways to do this partitioning. Two of these are described below (although the second is most developed):

  1. Electrostatics treated by QM; H junctions in MM1

  2. Electrostatics treated by MM; C junctions in MM1


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

  1. Vaccuum energies and QM ESP charges are calculated from a wavefunction optimised in vacuum.
  2. Vaccuum energies and QM ESP charges are calculated from a wavefunction polarised by the surrounding protein and solvent.
The former approach is faster and typically gives more stable energies, whereas the latter is at least theoretically more accurate.

Dependence of the results
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



How to calculate QM/MM-PBSA free energies
Simplified approach for version 2, by Ulf

  1. Run standard ComQum calculation.
    Follow the instructions in the ComQum page.

  2. If you only want QM/MM-PBSA energies on the minimised QM/MM structures, run the scripts
    qmmmpbsa-vac
    qmmmpbsa-ptch
    for each QM/MM state (in the QM/MM directory).
    It runs the vacuum and polarised calculations with version 2.
    The files and the results are put into the two subdirectories Vac and Ptch.
    No input is needed.
    It also gives updated charges and the vacuum strain energy.
    These scripts should always be run for any QM/MM calculation.

  3. If you have run protein free and want to have more stable energies (putting the QM systems into the MM systems on the other states and vice versa), run also either
    qmmmpbsa-vdir
    in the Vac subdirectory
    or
    qmmmpbsa-pdir
    in the Ptch subdirectory.
    Both programs can take some input data in the optional file qmmmpbsa-input:
    qmdir=../             # Directory for the QM Vac calculations
    mddir=../Md           # Directory for the MD simulations
    miss=0                # Number of missing atom
    The miss parameter allows there to be one more or less atom in one of the calculations (pKa estimation) and miss should be the number of this atom in the calculation with one more atoms.
    The results is put in the current directory, primarily in the file results.
    For the pdir calculation, you may also specify the directory of the pointcharges file to use:
    pcdir=$mddir/../      # Directory with the MM pointcharge file

    As usual, the results are in the file results in the current directory.
  4. 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).

  5. When finished, run either
    qmmdpbsa-vac
    or
    qmmdpbsa-ptch
    Again, you can specify non-default parameters in the file qmmmpbsa.input:
    qmdir=../../Vac     # Directory for the corresponding vacuum calculation
    mddir=..            # Directory for the MD simulations
    mdcrd=mdcrd5_wrap   # Name of the mdcrd  file to use (in .)

    prmtop=prmtop       # Name of the prmtop file to use (in .)
    first=1             # First snapshot to use
    last=' '            # Last  snapshot to use (default = all = <CR>)
    skip=0              # Number of snapshots to skip between each taken

  6. Finally, you can also use different MM and QM systems (to get more stable results.
    Unfortunately, the system is often translated and rotated by Amber. Therefore, it is not certain that the QM and MM systems are compatible. The best thing to do is to calculate the transformation matrix between the QM/MM coordinates and the final MD pdb5_wrap file and use the transformation file as an input for the qmmdpbsa scripts:

    changepdb <<EOF
    pdb3_cq

    rf              
    pdb5_wrap
    f
    b
    sander.in3
    n
    y

    y
    m
    c
    y
    transf5
    q
    EOF

  7. After this, you can run either
    qmmdpbsa-vdir
    or
    qmmdpbsa-pdir.
    Again, you can specify non-default parameters in the file qmmmpbsa.input:

    qmdir=../../Vac     # Directory for the corresponding vacuum calculation
    mddir=../           # Directory for the MD simulations
    mdcrd=mdcrd5_wrap   # Name of the mdcrd  file to use (in $mddir)
    prmtop=prmtop       # Name of the prmtop file to use (in $mddir)
    transf=tranfs5      # Name of the transformation file
    first=1             # First snapshot to use
    last=' '            # Last  snapshot to use (default = all = <CR>)
    skip=0              # Number of snapshots to skip between each taken
    miss=0              # Number of missing atom
    Please note that in this case, miss should be positive if the extra atom is in the QM system and negative if it is in the MM system.

    As usual, the results are in the file results in the current directory. However, to take advantage of the cancellation between similar terms, you need to run the program
    qmmmpbsa-2statistics
    reading in both the results file of the calculation with the QM and MM systems from the same calculation and the other QM system in the same MM system (e.g. Oo in Oo and Rw in Oo). If you are doing a redox potential or pKa calculation, you can specify this also (e = redox potential, p = pKa calculation, and c = coupled redox potential and pKa calculation). Then, it is important that the calculations with an extra electron or proton are given as the second results file.

    Quite often, this will lead to large vdWaals energy terms, which will give energetic outliers. These are automatically detected and removed by the program (6 sigma).



Old procedure (according to Markus)

  1. 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).

  2. 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).
  3. 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.pdb
    The 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.pdb
    Copy 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
  4. For version1 generate the prmtop2 file with zero charges (this should be the case already)

    changeparm

    prmtop3
    1
    w
    prmtop2
    q
  5. 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

  6. 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"

  7. 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.
  8. 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 fromprmtop3.

  9.  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 ###################"

  10. qsub run_QM_snapshots.ser

      1. TdS

    Run turbomole program aoforce on the QM system in vacuum

    there after run turbomole program freeh to get the entropy.

  11. 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.


Specific instructions for Version 2 directly on QM/MM structures

UR 5/12-08
  1. Run a standard ComQum calculation.
    Follow the instructions in the ComQum page.

  2. Construct a pdb file from the final structure:
    changeparm <<EOF
    prmtop3
    p
    pdb3
    m
    prmcrd3
    q
    EOF

  3. If you have run a protein fixed, you need to calculate charges for the QM system:

    1. mkdir Esp
    2. cp control coord basis energy mos alpha beta auxbasis comqum.dat pdb3 prmtop3 prmcrd3 Esp
    3. cd Esp
    4. Insert into control:
      $mvd
      $moments
      $pop
      $esp_fit
      shell 75000 -0.57
      shell 75000 -0.38
      shell 75000 -0.19
      shell 75000  0.00
      shell 75000  0.19
      shell 75000  0.38
      shell 75000  0.57
      shell 75000  2.46
      shell 75000  4.35
      shell 75000 11.91
      $vdw_radii
       h  3.85504
       c  4.81880
       n  4.81880
       o  4.49755
       s  5.62193
       mn 6.42507
       fe 6.42507
       co 6.42507
       ni 6.42507
       cu 6.42507
       zn 6.42507
    5. ridft -proper >logm

    6. changepdb <<EOF
      pdb3
      cc

      qcc

      t
      logm
      comqum.dat
      w
      pdb3
      q
      EOF
    7. Check that the charge of the pdb3 file is the same before and after this step
      changepdb <<EOF
      pdb3
      cc

      q
      EOF
    8. Insert these new charges into the prmtop file
      changeparm <<EOF
      prmtop3
      m
      pdb3
      w
      prmtop3
      q
      EOF
  4. Remove water molecules from the prmtop file
    changeparm<<EOF
    prmtop3
    t
    0
    y
    n
    prmtop3.nowat
    y
    m
    prmcrd3
    prmcrd3.nowat

    q
    EOF

    If you want to keep some water molecules, you have to use a belly file instead with the atoms to keep in the belly.

  5. Generate prmtop and prmcrd files for system 1 and C-junction atoms:

    changeparm <<EOF
    prmtop3
    t
    c
    n
    n
    prmtop1.cjunct
    y
    m
    prmcrd3
    prmcrd1.cjunct
    q
    EOF

  6. 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

  7. 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.

  8. 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

  9. 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

  10. Optionally, you may run a molsurf calculation to get the surface area of the protein for the non-polar energy:
    changeparm <<EOF
    prmtop3.nowat
    pqr
    m
    prmcrd3.nowat
    b
    1.4
    prmcrd3.nowat.pqr
    q
    EOF

    molsurf prmcrd3.nowat.pqr 0

    changepdb can also calculate a similar area:
    changepdb <<EOF
    prmcrd3.nowat.pqr
    so
    a
    b
    1.4
    0.2
    q
    EOF

    The non-polar energy is now the surface area * 0.0227 -3.85 kJ/mol.

  11. 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

------------------------------------------------------------------------------------------------------------------------------------------------

removewaters.belly

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

------------------------------------------------------------------

sander.in1

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

------------------------------------------------------------------

prepareversion1.scr

#!/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 .

prepareversion2_MD.scr

#!/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

------------------------------------------------------------------

ptraj_snapshots_pdb_noH2O.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

------------------------------------------------------------------

ptraj_snapshots_crd_noH2O.in

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

run_QM_snapshots.ser

#!/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}

run_tm_with_QMPC.scr

#!/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

run_tm_with_QMVP.scr

#!/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"


How to run QM/MM-PBSA without solvent accessible charges

UR 10/6-10
  1. Start from the QM/MM calculation
  2. Copy the files to a new directory

    mkdir Sae
    cd Sae
    cp ../* .

  3. Remove the solvent accessible charges and write out new prmtop3 and pointcharges files:

    changeparm <<EOF
    prmtop3
    ns
    solvacc
    w
    prmtop3
    pc
    m
    prmcrd3
    z1
    pointcharges
    EOF

  4. Optimise the wavefunction with the new pointcharges and then run qmmmpbsa

    ridft <logd
    qmmmpbsa-vac  >logv
    qmmmpbsa-ptch >logp