Hybrid QM/MM calculations

Now that we have seen how to carry out quantum and classical calculations within ChemShell, we can look at how to combine them to perform hybrid QM/MM calculations.

NB: as is conventional in the scientific literature we use the term QM/MM loosely to include any kind of hybrid quantum and classical calculation, regardless of whether it specifically uses a molecular mechanics forcefield for the classical part or not.

QM/MM calculations can be characterised by their level of embedding: mechanical, electrostatic or polarised. In this tutorial, we will look at the first two types which are the most commonly used. We will meet polarised embedding in later tutorials where shell model potentials are used.

The examples in this tutorial are in samples/qmmm_gulp/.

Mechanical embedding

In a mechanical embedding calculation, the QM part of the calculation is not directly affected by the MM part, and is performed as if the system is in the gas phase. The coupling of the QM and MM parts of the system is achieved purely at the MM level through assigning partial charges to the QM centres which interact electrostically with the MM centres. These charges may be already known from the forcefield definition, or may have to be derived.

In our example mech_butanol.nwchem.py and equivalents, we re-use the butanol molecule and forcefield we saw defined in the previous tutorial. However, in this case we define both a QM and MM calculation:

# Define QM and MM parts of the calculation
my_qm = NWChem(method="dft", functional="blyp", basis="3-21g")
my_mm = GULP(ff=ff, molecule=True)

As you can see, this is no different from the individual QM and MM calculations we have run before. In order to combine the two in a hybrid QM/MM calculation, we next create a QMMM() object which defines the coupling between the QM and MM parts:

# Define the combined QM/MM calculation
my_qmmm = QMMM(frag=my_mol, qm_region=range(5),
               qm=my_qm, mm=my_mm, embedding="mechanical")

Here we specify the QM and MM parts by passing the corresponding objects to the qm= and mm= options. The QM region must also be defined as a python list. Here we have used python’s range command to create a list of atoms 0-4. We also select the type of embedding to be used, in this case mechanical.

Finally, we request an energy minimisation with DL-FIND. Here we pass the QMMM object my_qmmm to the Opt() definition, rather than any of the individual QM and MM parts:

# Request an energy minimisation
my_opt = Opt(theory=my_qmmm,
            algorithm="lbfgs", tolerance=0.0003,
            maxcycle=300, maxene=400, trustradius="energy")
my_opt.run()

Notice when DL-FIND runs that each step now consists of a QM and MM calculation, and the individual and combined QM/MM energies are reported along the way, e.g.:

**********************************************************************
QM energy          =         -115.026242573788 a.u.
MM energy          =           -0.002321958949 a.u.
QM/MM total energy =         -115.028564532737 a.u.
**********************************************************************

You may find at times that either the QM or MM energy is rising during the optimisation, as it is the combined QM/MM total energy that is minimised.

Electrostatic embedding

The next level of embedding is known as electrostatic embedding, where the MM partial charges enter the QM Hamiltonian and so directly polarise the QM subsystem. This relies on the QM program used supporting calculations with additional point charges, but support for this is now almost universal.

One issue with electrostatic embedding is that the MM centres can overpolarise the QM subsystem through their proximity to link atoms that are added into the QM calculation wherever a covalent bond is broken. ChemShell handles this situation automatically through a “charge shift” scheme to reduce the polarisation effect - for details see the ChemShell publications.

In practical terms, you do not need to worry about this complexities when setting up the calculation. All that needs to be changed is the QMMM() definition as in our example elec_butanol.nwchem.py:

# Define the combined QM/MM calculation
my_qmmm = QMMM(frag=my_mol, qm_region=range(5),
               qm=my_qm, mm=my_mm, embedding="electrostatic", coupling="covalent")

Here embedding="electrostatic" specifies that electrostatic embedding should be used, and coupling="covalent" specifies that the system is covalently bonded and so the charge shift scheme should be implemented wherever a bond is broken. We will look at ionic systems, which are treated differently, in later tutorials.

Please note that the energies from calculations run at different levels of embedding cannot be directly compared, so you should not expect for example that an electrostatic embedding energy will necessarily be lower than a mechanical embedding energy calculated on the same system.