Calc QM/MM-2QM

This is an approach to calculate redox potentials and reorganisation energies by QM/MM-2QM calculations, based on MD simulations.

When using this implementation, cite this publication:
Hu, L.; Farrokhnia, M.; Heimdal, J.; Shleev, S.; Rulíšek, L.; Ryde, U. "Reorganisation energy for internal electron transfer in multicopper oxidases", J. Phys. Chem. B, 2011, 115, 13111-13126; DOI 10.1021/jp205897z.

1. Perform a PON parametrisation of the two electron-transfer sites of interest in both redox states (red and ox).
   
   
2. Run QM/MM-2QM optimisation of the two electron-transfer sites in both combinations of redox states (ox/red and red/ox; mainly to set up the needed QM files and get proper charges of the QM sites).
   
   
3. Run a MD simulation of the protein with the two electron-transfer sites in both combinations of redox states (ox/red and red/ox). Use periodic boundary conditions.
   Run ~10 ns and save snapshots every ~10 ps (i.e. in total ~1000 snapshots).
   Sample sander.in files are given below.
   
   
4. Set up the calcqmmm-2qm calculations.
   In the following we will have two sets of redox states, which we call 0 and 1.
   Set the directories of the two MD simulations and the two QM/MM calculations: $md0, $md1, $qm0, $qm1.
   Insert these four commands into a file call dirs (so that you remember what files were used):
   md0=/lunarc/nobackup/users/geng/Azar/Md/Ro
   md1=/lunarc/nobackup/users/geng/Azar/Md/Or
   qm0=/lunarc/nobackup/users/geng/Azar/Qmmm/Ro
   qm1=/lunarc/nobackup/users/geng/Azar/Qmmm/Or
   
   Check that the directories are correc:
   ls $md0 $md1 $qm0 $qm1
   
   
5. Then do all the following calculations in the new directory where the calculations will be done:
   \cp $qm0/comqum.dat? .
   \cp $qm0/basis? .
   \cp $qm0/auxbasis? .
   \cp $qm0/coord? .
   \cp $qm0/control1 control01
   \cp $qm0/control2 control02
   \cp $qm1/control1 control11
   \cp $qm1/control2 control12
   \cp $qm0/pointcharges1 pointcharges01
   \cp $qm0/pointcharges2 pointcharges02
   \cp $qm1/pointcharges1 pointcharges11
   \cp $qm1/pointcharges2 pointcharges12
   \cp $qm0/mos1 mos01
   \cp $qm0/mos2 mos02
   \cp $qm1/mos1 mos11
   \cp $qm1/mos2 mos12
   \cp $qm0/alpha1 alpha01
   \cp $qm0/alpha2 alpha02
   \cp $qm1/alpha1 alpha11
   \cp $qm1/alpha2 alpha12
   \cp $qm0/beta1 beta01
   \cp $qm0/beta2 beta02
   \cp $qm1/beta1 beta11
   \cp $qm1/beta2 beta12
   \cp $md0/prmtop prmtop02
   \cp $md1/prmtop prmtop12
   
   Some warnings for missing mos, alpha, or beta files are expected.
   
   \cp $md0/mdcrd .
   
   cat <<EOF >sander_fep1
   This is sander_fep1
    &cntrl
     irest=0,ntx=1,
     imin=1,ntmin=2,maxcyc=0,
     temp0=300.0,ntt=1,tautp=1.0,
     ntc=2,ntf=2,
     cut=1000.0,
     ntpr=1,ntwx=0,ntwv=0,ntwe=0,
     ntb=0,ipol=0,igb=0,
     ntr=0,ibelly=0
    &end
   EOF
   
   cat <<EOF >sander_fep2
   This is sander_fep2
    &cntrl
     irest=0,ntx=1,
     imin=1,ntmin=2,maxcyc=0,
     temp0=300.0,ntt=1,tautp=1.0,
     ntc=2,ntf=2,
     cut=8.0,
     ntpr=1,ntwx=0,ntwv=0,ntwe=0,
     ntb=1,ipol=0,igb=0,
     ntr=0,ibelly=0
    &end
   EOF
   # Note that this file is with PBC and PME; if possible, it is probably better to let sander_fep2 = sander_fep1
   
   sed -i "s/energy1/energy01/" control01
   sed -i "s/pointcharges1/pointcharges01/" control01
   sed -i "s/mos1/mos01/" control01
   sed -i "s/alpha1/alpha01/" control01
   sed -i "s/beta1/beta01/" control01
   
   sed -i "s/energy2/energy02/" control02
   sed -i "s/pointcharges2/pointcharges02/" control02
   sed -i "s/mos2/mos02/" control02
   sed -i "s/alpha2/alpha02/" control02
   sed -i "s/beta2/beta02/" control02
   
   sed -i "s/energy1/energy11/" control11
   sed -i "s/pointcharges1/pointcharges11/" control11
   sed -i "s/mos1/mos11/" control11
   sed -i "s/alpha1/alpha11/" control11
   sed -i "s/beta1/beta11/" control11
   
   sed -i "s/energy2/energy12/" control12
   sed -i "s/pointcharges2/pointcharges12/" control12
   sed -i "s/mos2/mos12/" control12
   sed -i "s/alpha2/alpha12/" control12
   sed -i "s/beta2/beta12/" control12
   
   # Zero the charges of the QM systems in the two full prmtop files
   ln -fs comqum.dat1 comqum.dat
   changeparm <<EOF
   prmtop02
   1
   w
   prmtop02
   q
   EOF
   changeparm <<EOF
   prmtop12
   1
   w
   prmtop12
   q
   EOF
   ln -fs comqum.dat2 comqum.dat
   changeparm <<EOF
   prmtop02
   1
   w
   prmtop02
   q
   EOF
   changeparm <<EOF
   prmtop12
   1
   w
   prmtop12
   q
   EOF
   
   # Construct the truncated prmtop files
   ln -fs comqum.dat1 comqum.dat
   changeparm <<EOF
   prmtop02
   tr
   comqum.dat
   n
   n
   prmtop011
   n
   q
   EOF
   changeparm <<EOF
   prmtop12
   tr
   comqum.dat
   n
   n
   prmtop111
   n
   q
   EOF
   ln -fs comqum.dat2 comqum.dat
   changeparm <<EOF
   prmtop02
   tr
   comqum.dat
   n
   n
   prmtop012
   n
   q
   EOF
   changeparm <<EOF
   prmtop12
   tr
   comqum.dat
   n
   n
   prmtop112
   n
   q
   EOF
   \rm comqum.dat
   
   # Construct the input file
   cat <<EOF >calcqmm.in
   Title
   mdcrd
   comp
   1
   
   0
   o     (octahedral or rectangular box?)
   EOF
   
   
6. Run the program calcqmmm-2qm
   If run in a batch queue, it is normally wise to copy the mdcrd to the Calcqmmm directory and run everything on the local disk.
   calcqmmm-2qm <calcqmm.in >$SLURM_SUBMIT_DIR/calcqmm.out
   
   Also include the comp file in the slurm_save_files file.
   echo comp >slurm_save_files
   
   
7. When running the opposite redox state, all files for calcqmmm-2qm could be copied, except the mdcrd files, which should be replaced by those of the $md1 stat.
   
   
8. The final result are the averages over the columns in the comp file, especially the first column, total.
   You need to run calcqmmm-2qm for both the red/ox and ox/red MD simulations.
   E0 = ( Av(ro) + Av(or) )/2
   ReorgE = ( Av(ro) - Av(or) )/2
   
   For electrostatic embedding (the pointcharges files are not empty, as it is set up in these instructions), all bond, angle and dihedral energies should be small (< 1kJ/mol), the van der Waals energy is sizeable (-3 to 24 kJ/mol) and the electrostatic energies are exactly = 0. However, if you use periodic boundary conditions and Ewald summation, the electrostatic term can be 1-2 kJ/mol.
   
   For mechanical embedding (the pointcharges files are empty), all bond, angle and dihedral energies should be very small (+- 0.001kJ/mol), the van der Waals energies are rather big (-9 to 1 kJ/mol) and the electrostatic energies are very large = (172-303 kJ/mol).
   
   If the bonded terms are not < 1 kJ/mol, you have probably differences between the two systems outside the QM systems.
   
   If the bonded terms are exactly = 0, you probably have no difference in the force field of the two QM sites, besides the charges.
   
   If elstat term is not zero, build pdb files from the prmtop02 and prmtop12 files and compare the charges.
   
   The result in the comp file can be reproduced by those in the calcqmm.out file:
   It is always the difference between the 1 and 0 calculations.
   Tot is the total QM/MM energy, i.e. QM1+QM2+MM2-MM11-MM12
   MM11=Bond11+Angle11+Dihed11+vdWaals11+Elstat11
   MMtot=MM2-MM11-MM12 (and similar for Bond, Angle, Dihed, vdWaals, and elstat)
   So the Bon, Ang, Dih, vdW, and Elstat energies are got in that way.
   

Sample sander files (essentially taken from QTCP)