.. _qm_mm_spoint:

*******************************
Single-Point Energy Calculation
*******************************

.. note:: The scripts for this step of the tutorial can be found in ``mgo_surface/2_mgo_sp``.

Now that the embedded cluster model has been set up, it is ready for QM/MM 
simulations. In this tutorial, we demonstrate how 
a single-point QM/MM energy evaluation can be performed on the model structure.

`1. Quantum Mechanics (QM)`_
============================

.. note:: The QM calculations in this example use NWChem.

As with all QM packages, an NWChem calculation requires a basis set to be 
specified. Basis sets for the MgO system can be 
obtained from the `Basis Set Exchange <http://www.basissetexchange.org>`_. 

For Mg, we will use the Stuttgart relativistic, large core basis which 
is designed to be used with an effective core potential (ECP). 
The ECP is a representative potential which replaces some core electrons,
thereby meaning a smaller number of basis functions can be used to represent 
the system.

.. note:: Technical note: the Stuttgart basis was modified to remove the most
  diffuse functions in order to help prevent artificial spreading of charge density 
  outside the QM region.

For O, we will use a standard Ahlrichs pVDZ double zeta basis with polarisation 
functions. We will also use this basis for the adorbate in later stages 
of the tutorial.

The basis set an ECP specifications are provided to ChemShell as separate 
**.basis** and **.ecp** files in NWChem format.

Below are the ``mgo_nwchem.basis`` and ``mgo_nwchem.ecp`` files used:

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_nwchem.basis
.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_nwchem.ecp

The ECP fulfills two purposes in the ionic embedded cluster model:

 * In the inner QM region (Mg1), it is used in a conventional way
   to replace core electrons 
   (specified by ``nelec 10``, calculated as 12 Mg electrons minus 
   2 from the valence orbital) 
 * In the QM/MM boundary region (bq_Mg2_e), 
   the ECP acts as a boundary potential to prevent charge density diffusing 
   outside the QM region. In the boundary region, no full QM atoms are present 
   so no core electrons are replaced 
   (as specified by ``nelec 0``).

With the exception of the number of electrons replaced, the ECP definitions
for region 1 and region 2 are identical.

`2. Molecular Mechanics (MM)`_
==============================

.. note:: The MM calculations in this example use GULP. We use the term MM 
   here loosely, as is usual in the QM/MM literature. The interatomic 
   potentials used in the ionic embedding model are quite different in form
   to molecular mechanics forcefields used for covalently-bonded systems,
   but they fulfill the same purpose.

For the classical environment, we need to specify a forcefield in GULP format.
This can be sourced in one of two ways: either using a pre-existing forcefield
from the literature or generated manually using the approach described in the
:doc:`ffield` tutorial. The force field used in this tutorial, 
``mgo_2body.ff``, is shown below:

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_2body.ff

If you are familiar with GULP calculations, you will recognise standard 
forcefield terms such as Buckingham ``buck``, Lennard-Jones ``lennard``, 
and Coulomb ``coulomb`` interatomic potentials. There are also Spring ``spring`` 
terms to control the interaction between classical cores and shells. 
More information on setting these terms can be found in the 
`GULP manual <https://gulp.curtin.edu.au/gulp/help/manuals.cfm>`_.

The forcefield description also depends on which region the atoms are in.
The parameters for each term do not vary by region, but certain terms 
are not necessary. For example, there are no O ions in region 2, so there 
are no region 2-2 terms (as there are no Mg-Mg terms in the forcefield).
All terms involving O ions are between shells, except for region 1 O ions
(where shells have been deleted). Similarly, spring terms are only necessary 
for region 3. Finally, interactions within regions 4 and 5 (including shell
interactions) are excluded as these regions are fixed.


`3. Calculating the Single-Point Energy`_
=========================================

The input structure used in this calculation is the ``mgo_shells_region.pun`` 
created previously. The coordination table used in the 
cluster creation step is no longer required, and so it can (optionally) 
be deleted by setting ``connect_mode`` to **None**.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 7-8

Next, we read in the QM region and charge as defined in the saved 
fragment.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 10-15

The ``getRegion()`` method is used to obtain a list of atoms with the 
specified region label (here region 1, the QM region).

For the QM charge, we round the floating point value in the saved 
fragment to the nearest integer for use in the QM calculation.

The settings for the QM/MM calculation can then be defined.
First, we set the options for the QM (NWChem) calculation.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 18-23

Here, we have specified a density functional theory calculation using 
the well-known B3LYP functional, with the basis set and ECP files we
discussed above. The saved QM charge is passed in. As this is a 
singlet calculation, the spin multiplicitly can be left as the default (1).

Special consideration should be given to the ``path`` option, which 
we have set here to call ``mpirun`` in order to run NWChem in parallel with 
4 cores. If you do not have an MPI environment set up on your  machine, 
you should remove the 
``path`` option so that NWChem can be run directly. The calculation 
will of course take longer as a result. On the other hand, if you have 
more than 4 cores available, you could experiment with increasing the 
number of cores passed to ``mpirun``.

.. note:: These instructions assume you are running ChemShell and NWChem 
   as separate executables. If you are running ChemShell with NWChem directly-linked, 
   however (as will be the case on most HPC platforms), the path 
   option should be removed as the NWChem calculation will use the same 
   MPI environment as ChemShell.

Next, the options for the GULP part of the calculation are given.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 25-26

These are quite straightforward, specifying the file containing the 
forcefield definition, and the optimisation algorithm for shell 
relaxation.

Finally, the options for the QM/MM driver are given.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 28-33

Here, the QM and MM theory options specified above are passed in, together 
with the QM region definition obtained above, and the QM/MM coupling model 
used (in the case of MgO, the ionic embedding model).
Finally, we set the maximum number of shell relaxation cycles to the 
relatively high value of 20, to ensure that the shells are fully converged 
for the final energy evaluation.

The QM/MM calculation is now ready to run. We set up a single-point 
calculation by passing the QM/MM options to the  ``SP()`` class, and then call
the ``run`` method to run the calculation.

.. literalinclude:: ../../samples/mgo_surface/2_mgo_sp/mgo_sp_energy.py
   :lines: 35-37

The remainder of the input file is a validation step to check that the 
calculation you obtain is the same as our pre-calculated value.
This is purely for testing purposes and would not normally form part 
of a production calculation.

.. note:: If your calculated single-point energy differs significantly from 
   the pre-calculated value and so fails the test, it may indicate 
   that there is a problem with your ChemShell, NWChem or GULP installations.
