Big-QM

The Big-QM approach is a method to post-process QM/MM calculations to get more stable energies.
It involves calculations with large-QM systems (600-1000 atoms).

The atoms are selected by the following rules:
The approach is described in


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

How to set up big-QM calculations with an existing syst1 file
This is done for each calculation of a project
  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)


Calculate mechanical-embedding correction for a vacuum calculation
The only difference from electrostatic embedding is that the charges of the QM system should not be zeroed.
Run points up to 10 in the previous scheme.
  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



Run the MM and disp calculations by hand

comqumtoturbo
cp ../prm*3 .
cqtrunc
cqzero1
cqzero3
coord2dftd
dftd3 coord-dftd3 -func tpss -bj -abc >dftd3.out
calcforce <<EOF >calcforce.out1
prmtop1
m
prmcrd1
1000
mden1
force1
EOF
calcforce <<EOF >calcforce.out2
prmtop2
m
prmcrd3
1000
mden2
force2
EOF

grep Edisp dftd3.out; grep 'Total energy' calcforce.out*