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