Big-QM

• All atoms within 4.5-6 Å of a minimal QM system
• All buried charges in the protein
• Move all junctions two residues away from the minimal QM system
  

• L. Hu, P. Söderhjelm, U. Ryde "Accurate reaction energies in proteins obtained by combining QM/MM and large QM calculations" J. Chem. Theory Comput. 2013, 9, 640-649; DOI 10.1021/ct3005003.
• Sumner, S.; Söderhjelm, P.; Ryde, U., "Effect of geometry optimisations on QM-cluster and QM/MM studies of reaction energies in proteins", J. Chem. Theory Comput., 2013, 9, 4205-4214; DOI: 10.1021/ct400339c.

How to set up the syst1 file for Big-QM calculations

This typically done only once for each project/protein.
It is essential to use exactly the same atoms in the big-QM system for all calculations of the same project.

1. Start from the pdb3 file from a QM/MM calculation with a minimal (or reasonable) QM system, together with the syst1 file for that calculation (called s1 below).
   
   
2. Run changepdb, command bc to find the buried charges of the protein. Use the actual charges.
   For each residue in the list of buried charges, add a note if it forms an ionic pair with another residue in the protein (also given in the output).
   Check all residue that do not form such ionic pairs or are not part of the QM system or binds to a metal site with a pdb viewer.
   Decide whether they should be charged or not (in particular, check for water-mediated H bonds)
   and whether they are buried or not.
   All this should in principle be done already when the MD simulation is set up (to get proper protonation state of all residues).
   
   
3. Run changepdb on the same protein, command bq.
   This is best done by setting up an input text file:
   changepdb <bigqm.in >bigqm.out
   
   Enter the protein file
   bq
   s1
   radius (4.5-6)
   3
   enter a list of the buried residues (residue numbers, one on each line, finished by a blank line)
   optionally add additional atoms (e.g. bridging water molecules, atom numbers one on each line)
   Hit return (possibly a few times) until finished
   syst1-bigqm
   q
   
   It is quite hard to find the charge of the QM system
   This gives the number of positive and negative amino acids
   
   grep ASP bq.pdb | grep CG; grep ASP bq.pdb | grep CG | wc
   grep GLU bq.pdb | grep CD; grep GLU bq.pdb | grep CD | wc
   grep ARG bq.pdb | grep CZ; grep ARG bq.pdb | grep CZ | wc
   grep LYS bq.pdb | grep NZ; grep LYS bq.pdb | grep NZ | wc
   grep HIP bq.pdb | grep NE; grep HIP bq.pdb | grep NE | wc
   grep OXT bq.pdb; grep OXT bq.pdb|wc
   grep H3 bq.pdb; grep H3 bq.pdb|wc
   
   To these numbers, you should add the charge of the minimal QM system, as well as other metals or non-standard residues.
   
   
   
4. Set up the big-QM calculations by pdbtocomqum
   
   pdbtocomqum <<EOF
   logfil
   comqum.pdb
   
   
   syst1-bigqm
   0
   0
   n
   n
   1000
   a
   -1
   n
   n
   n
   EOF
   
   The charge of the QM system is quite hard to estimate automatically.
   Normally, the column best gives the best estimate. Insert it into a spread sheet and check each charge thoroughly before summing all charges.
   However, check it by calculating the charge by hand also.
   
   
5. Run comqumtoturbo
6. Run rasp to check the QM system
7. Run comqumtoamber
   
8. Copy the prmtop3 and prmcrd3 files from the previous QM/MM calculation (or set up new ones as described in the comqum page).
9. Run cqprep to set up the QM/MM calculation. No cap is needed and sander.out3 can be ignored.
   Set up the dft calculation as usual, but skip ff and add also marij.
10. Change the dispersion line in the control file to
    $disp3 -bj -abc
    
11. Run ridft
12. Read the standard MM correction from the calcforce.out1 and calcforce.out2 files.
    grep 'Total energy               =' calcforce.out?
    The final big-QM energy is then:
    Etot = E(QM+D3bj123) + E(MM12 from calcforce.out2) - E(MM1 from calcforce.out1)

1. Run pdbtocomqum
   
   pdbtocomqum <<EOF
   logfil
   comqum.pdb
   
   
   syst1-bigqm
   0
   0
   n
   n
   1000
   a
   -1
   n
   n
   n
   EOF
   
   The charge of the QM system is quite hard to estimate automatically.
   Normally, the column best gives the best estimate. Insert it into a spread sheet and check each charge thoroughly before summing all charges.
   However, check it by calculating the charge by hand also.
   
   
2. Run comqumtoturbo
3. Possibly, run rasp to check the QM system
4. Run comqumtoamber
   
5. Copy the prmtop3 and prmcrd3 files from the QM/MM calculation (or set up new ones as described in the comqum page).
6. Run cqprep to set up the QM/MM calculation.
   No cap is needed and sander.out3 can be ignored.
   Set up the dft calculation (define) as usual, but skip ff and add also marij.
7. Change the dispersion line in the control file to
   $disp3 -bj -abc
   
8. Run ridft
9. If the calcforce.out1/2 files do not exist, run calcforce:
   calcforce <<EOF >calcforce.out1
   prmtop1
   m
   prmcrd1
   1000
   mden1
   force1
   EOF
   less calcforce.out1
   calcforce <<EOF >calcforce.out2
   prmtop2
   m
   prmcrd3
   1000
   mden2
   force2
   EOF
   less calcforce.out2
   
10. Read the standard MM correction from the calcforce.out1 and calcforce.out2 files.
    grep 'Total energy               =' calcforce.out?
    The final big-QM energy is then:
    Etot = E(QM+D3bj123) + E(MM12 from calcforce.out2) - E(MM1 from calcforce.out1)

1. Delete the line in the control file:
   $point_charges  file=pointcharges
   
2. rm pointcharges
3. changeparm <<EOF
   prmtop3
   tr
   comqum.dat
   n
   y
   y
   prmtop1-me
   y
   m
   prmcrd3
   prmcrd1-me
   q
   EOF
4. calcforce <<EOF >calcforce.out1
   prmtop1
   m
   prmcrd1
   1000
   mden1
   force1
   EOF
   less calcforce.out1
   calcforce <<EOF >calcforce.out2
   prmtop2
   m
   prmcrd3
   1000
   mden2
   force2
   EOF
   less calcforce.out2