### 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)

This is sander.in0 (minimisation with restraints)
&cntrl
irest=0,ntx=1,
imin=1,maxcyc=1000,drms=0.0001,
ntc=2,ntf=2,
cut=8.0,
ntpr=100,ntwx=0,ntwv=0,ntwe=0,
ipol=0,igb=0,
ibelly=0,
ntr=1,restraint_wt=100,
&end

This is sander.in1 (NTV restrained equilibration)
&cntrl
irest=0,ntx=1,
nstlim=40000,dt=0.002,
temp0=300.0,ntt=1,tautp=1.0,
ntc=1,ntf=1,
cut=8.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
ntb=1,
ipol=0,igb=0,
ibelly=0,
ntr=1,restraint_wt=500,
&end

This is sander.in2 (NTP restrained equilibration)
&cntrl
irest=0,ntx=5,
nstlim=10000,dt=0.002,
temp0=300.0,ntt=1,tautp=1.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ibelly=0,
ntr=1,restraint_wt=100,
&end

This is sander.in3 (NTP box equilibration)
&cntrl
irest=1,ntx=5,
nstlim=25000,dt=0.002,
temp0=300.0,ntt=1,tautp=1.0,ntt=1,tautp=1.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=1000,ntwx=0,ntwv=0,ntwe=0,
ntb=2,ntp=1,pres0=1.0,taup=1.0,
ipol=0,igb=0,
ibelly=0,
ntr=0
&end

This is sander.in4 (NTV equilibration)
&cntrl
irest=1,ntx=5,
nstlim=100000,dt=0.002,
temp0=300.0,ntt=1,tautp=1.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=500,ntwx=0,ntwv=0,ntwe=0,
ntb=1,
ipol=0,igb=0,
ibelly=0,
ntr=0
&end

This is sander.in5 (NTV production with sampling)
&cntrl
irest=1,ntx=5,
nstlim=5000000,dt=0.002,
temp0=300.0,ntt=1,tautp=1.0,
ntc=2,ntf=2,
cut=8.0,
ntpr=5000,ntwx=5000,ntwv=0,ntwe=5000,iwrap=1,
ntb=1,
ipol=0,igb=0,
ibelly=0,
ntr=0
&end