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