ComQum-2QM is a special version of the ComQum QM/MM software, designed for the optimisation of two QM systems in the same protein, but treating them in separate QM systems (typically because otherwise a certain electronic state cannot be found).

The method is published in: Hu, L.; Farrokhnia, M.; Heimdal, J.; Shleev, S.; Rulísek, L.; Ryde, U. "Reorganisation energy for internal electron transfer in multicopper oxidases", J. Phys. Chem. B, 2011, 115, 13111-13126;DOI 10.1021/jp205897z.

The interface programs are located in:

• /temp4/bio/Bin/
  

To run ComQum, you also have to be able to run Turbomole and Amber .

The only available documentation is this file, the above reference, the ComQum page and a detailed specification of the ComQum file formats: comqum01.pdf .

1. Start with an equilibrated crystal structure (with hydrogen atoms) stored in pdb format with charges.
 
2. Change the name of metal-bound water molecules to WAM and move them before the other water molecules.
   
   
3. When selecting the quantum system, note that if there is a hydrogen bond from the back-bone NH or CO groups of a junction residue to the quantum system, this group must be included in the quantum system (because they otherwise get strange charges).
   
   
4. Ensure that you have the desired junctions in the file junctfactor
   (cp /away/bio/Data/Comqum/junctfactor .)
   
   For each junction, it should contain the following entries:
   Residue name, junction atom, junction neighbour, desired Amber atom type of junction atom;
   Next row: jtype (indicating method and basis set) and ideal bond lenght with this method
   The format is free, but fields must be delimited by spaces.
   
   $version
     2006
   
   ASP CB   CA   HC  # Asp OOC-CH3 junction, ff94
   7   1.1065d0      # Acetate, BPRi/6-31G*, Cs, C-H in plane
   
   
5. Make two directories, one for each QM system and copy the pdb file to both.
   
   
6. For each of the QM systems: Run pdbtocomqum
   
1. Select a logfile
2. Enter the name of the pdb file
3. Do not search for short contatcts
4. Enter a new title (can be done in a file, see next point).
5. Select the quantum system, including the junction atoms that will be changed to hydrogens (you first give the number of the amino acid, then the number of the atoms in this amino acid; one <CR> let you select a new amino acid; another <CR> ends the selection).
This can optionally be done by entering the name of a file with the atom numbers (see below). This file should also contain the new title as the first line.
6. Select the radius for system 2 (MM system that is optimised). Typically ~6 Å.
Note that it is absolutely essential that system 2 contains the same residues in all QM/MM structures if energies should be comparable. Therefore, you may add or remove residues from the originally suggested list.
7. Select the radius for system 3 (MM system that is fixed). Typically the rest of the protein.
8. Select amber (default)
9. Select the basis set (i.e. the type of junction, typically 6 (B3LYP) or 7 (BP/6-31G*)).
10. Enter the junction atoms (those that will become hydrogen atoms in the quantum system)
Check that they are not unknown, i.e. that the junction type is not -1 at the end (if so, look at point 6 below how to insert a new fragment into the libraries).
11. Remove the charge of neighbours.
12. Check that the net charge outside system 1 is correct (for normal side-chains, it should be 0, but for backbone fragments and strange groups, it may be <>0).
13. Check that the change in the charges is not too large. If so, remove more neighbours. For example, you may want to remove charges of hydrogens if you have already removed the charge of their parent atom.
 
7. For each of the QM systems: Run comqumtoturbo and  comqumtoamber
   
   
8. For each of the QM systems: Run define.
   First check that metal names are OK in file coord.
   In define:
   Run ff in the coordinate menu (and quit it by * and then no),
   Set the basis set
   Set the wavefunction (eht)
   Possibly change the DFT functional, include RI and turn on dispersion (dsp).
   
   Check that $grad is still present in the control file.
   
   Insert $esp_fit kollman into the control file
   
   Change all file references (i.e. file = xxx) in the control file, by adding 1 to all files in the first QM system and 2 in the second QM system, e.g.
   $coord    file=coord1
   (~12 times)
   
   
9. For only the first of the QM systems: Edit the file leap_commands to insert unusual residues (i.e. insert extra loadAmberPrep and  loadAmberParams statements).
   Copy all *.in and *.dat files to the ComQum directory from /away/bio/Amber/Data  (these files tend to change).
   Then run the tleap using these commands.
   tleap -f leap.in
   Check that no new atoms are defined.
   
 
10. For the same QM system: If you use non-standard charges for any residue (e.g. metal sites), change these charges by changeparm command m on prmtop3, using a pdb file with the correct charges (probably comqum.pdb). This is necessary only if you will run with the protein free. Note that leap will change back to standard charges, even if the charges were modified in the input pdb file.
    
    changeparm <<EOF
    prmtop3
    m
    comqum.pdb
    w
    prmtop3
    q
    EOF
    
    
11. Run cq2qmprep in the base directory and with the directories of the two QM systems as options, e.g.
    cq2qmprep T1 T23
    Follow the instructions.
    You will do the following:
    1. Construct prmtop1 and prmcrd1 for the first QM system.
    2. Zero charges in the first prmtop1
    3. Zero the charges of the first QM system in prmtop2
    4. Copy files to main directory
    5. Construct prmtop1 and prmcrd1 for the second QM system.
    6. Zero charges in the second prmtop1
    7. Zero the charges of the second QM system in prmtop2
    8. Copy files to main directory
    9. Insert cap into prmtop3 (this is done manually and you need to know the values to insert from the equilibration). It is important that you do not have any cap before running this command.
    10.Construct the sander.in2 and sander.in3 files by merging the bellies.
    11. Run a test of sander.in11 (write down the number of atoms and the energy)
    12. Run a test of sander.in12 (write down the number of atoms and the energy)
    13. Run a test of sander.in2 (write down the number of atoms and the energy)
    14. Run checkcoord to see that the correspondance list works (check that there are no warnings (Amber may change the order of the atoms; if so, you have to reorder them also in the original pdb file and rerun everyting from pdbtocomqum (point 3). There should be no warnings. If you get warnings forr HB-atoms, e.g. in CYM, run fixcoord1 (fixcoord1_amberin; fixcoord1_turboin; fixcoord1; fixcoord1_amberout) and see if the program gives errors. If so, you have to swap coordinates in the prmcrd1 file.
    15. Run a test of sander.in3 (This job needs to bee killed after ~30 seconds. Write down the number of atoms and the energy)
    
    If some of the QM systems does not contain junctions, the construction of the prmtop1 file will fail (and therefore also sander and checkcoord.
    Then, you have to run cq2qmprep by hand instead.
     
12. Zip unneccesary files (an make a backup) by
    gzip */* &
    
    
13. Start the comqum calculation with:
    comqum-2qm -backup -c 200
    or (with the RI approximation):
    comqum-2qm -ri -backup -c 200
    
    
14. When the job is finished, the final coordinates are in prmcrd3, coord1, and coord2.
    You may build pdf files (pdb11, pdb12, and pdb3) from these by also done automatically in point 16)
    cq2topdb
    
    
15. If you want run a protein free calculations (i.e. relax the surrounding protein in system 2), you set
    $protein free
    in comqum.dat
    and then run the comqum calculation again.
    
    
16. You should typically run qmmmpbsa on all ComQum results (it takes only a few hours).
    This is done by
    qmmmpbsa-vac-2qm >logv
    qmmmpbsa-ptch-2qm>logp
    Note that you need the files, so you cannot copy the calculation to the temp disk of the batch computer.
    
    

# Start with a scf step then loop through the following four steps

# Gradient step (gradient calculation)
ln -fs control1 control
grad  >logg1
ln -fs control2 control
grad  >logg2
ln -fs force11 force1
fixforce_amberin >fixf
ln -fs control1    control
ln -fs comqum.dat1 comqum.dat
cp gradient1 grad1
offsetf  >>fixf
fixforce_turboin>>fixf
fixforce>>fixf
fixforce_turboout>>fixf
mv gradient g11
mv grad1 gradient1

ln -fs force12 force1
fixforce_amberin >>fixf
ln -fs control2    control
ln -fs comqum.dat2 comqum.dat
cp gradient2 grad1
offsetf  >>fixf
fixforce_turboin>>fixf
fixforce>>fixf
fixforce_turboout>>fixf
mv gradient g12
mv grad1 gradient2

# Relaxation step
ln -fs control1 control
relax  >logr
ln -fs control2 control
relax  >logr
ln -fs control1    control
ln -fs comqum.dat1 comqum.dat
ln -fs prmcrd11 prmcrd1
ln -fs prmtop11 prmtop1
fixcoord1_amberin >fix1
fixcoord1_turboin>>fix1
fixcoord1>>fix1
fixcoord1_amberout>>fix1
ln -fs control2    control
ln -fs comqum.dat2 comqum.dat
ln -fs prmcrd12 prmcrd1
ln -fs prmtop12 prmtop1
fixcoord1_amberin >fix1
fixcoord1_turboin>>fix1
fixcoord1>>fix1
fixcoord1_amberout>>fix1
# With 2 QM, we must also run fixcoord2 after both fixcoord1, because the QM systems are part of the MM system of the other calculation
\cp prmcrd3 mdrest3
ln -fs control1    control
ln -fs comqum.dat1 comqum.dat
ln -fs pointcharges1 pointcharges
fixcoord2_amberin >fix2
fixcoord2_turboin>>fix2
fixcoord2>>fix2
fixcoord2_turboout>>fix2
ln -fs control2    control
ln -fs comqum.dat2 comqum.dat
ln -fs pointcharges2 pointcharges
fixcoord2_amberin >fix2
fixcoord2_turboin>>fix2
fixcoord2>>fix2
fixcoord2_turboout>>fix2

# Molmech step
# if  protfree then
   sander -O -i sander.in3 -o sander.out3 -r mdrest3 -c prmcrd3 -p prmtop3
   cp mdrest3 prmcrd3
   ln -fs control1    control
   ln -fs comqum.dat1 comqum.dat
   ln -fs pointcharges1 pointcharges
   fixcoord2_amberin >fix2
   fixcoord2_turboin>>fix2
   fixcoord2>>fix2
   fixcoord2_turboout>>fix2
   ln -fs control2    control
   ln -fs comqum.dat2 comqum.dat
   ln -fs pointcharges2 pointcharges
   fixcoord2_amberin >fix2
   fixcoord2_turboin>>fix2
   fixcoord2>>fix2
   fixcoord2_turboout>>fix2
#endif

# Scf step (energy calculation)
\rm forcedump.dat
sander -O -i sander.in1 -o sander.out11 -r mdrest1 -e mden11 -c prmcrd11 -p prmtop11
\mv forcedump.dat force11
\cp force11 mden11
sander -O -i sander.in1 -o sander.out12 -r mdrest1 -e mden12 -c prmcrd12 -p prmtop12
\mv forcedump.dat force12
\cp force12 mden12
sander -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
\mv forcedump.dat force2
\cp force2 mden2
ln -fs control1 control
ridft>logd1
ln -fs control2 control
ridft>logd2

# ComQum: Fix Charge; For 2QM, charges should always be udated, even with protein fixed (because they are part of MM2 for the other QM system)
ln -fs control1 control
ln -fs comqum.dat1 comqum.dat
sed -i 's/$point_charges/#point_charges/' control1
sed -i 's/#esp_fit/$esp_fit/' control1
ridft -proper>mulliken1
sed -i 's/#point_charges/$point_charges/' control1
sed -i 's/$esp_fit/#esp_fit/' control1
ln -fs control2 control
ln -fs comqum.dat2 comqum.dat
sed -i 's/$point_charges/#point_charges/' control2
sed -i 's/#esp_fit/$esp_fit/' control2
ridft -proper>mulliken2
sed -i 's/#point_charges/$point_charges/' control2
sed -i 's/$esp_fit/#esp_fit/' control2
fixcharge-2qm >>$output

ln -fs comqum.dat1 comqum.dat
fixenergy_2qm >fixe
checkconv

Sample comqum.dat file

     1      2     1.09800000000000 H1 
$correspondance_list
    15
  1976  1978  1979  1980  1981  1982  1983  4865  4866  4867  4868  4869
  4870  4871  4872
$junction_atoms
     1     2    1.38979963570127
$end

Note that $junctions is added by pdbtocomqum, whereas $junction_atoms is added by changeparm, command tr.