.. _setup_cluster:

****************************
Setting up the QM/MM Cluster
****************************

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

Preparation of a surface cluster
================================

This section of the tutorial will outline how to set up the QM/MM cluster 
using the input script ``construct_cluster.py``. 

The periodic crystalline structure of interest is first read in  
to ChemShell, which is stored as a 2D slab in the file 
``mgo_shells_periodic.pun`` in ChemShell's punch file format. 
This is imported by ``construct_cluster.py`` as 
follows:

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/construct_cluster.py
    :lines: 6-7

where ``mgo`` is a Fragment containing the periodic input structure in 
``mgo_shells_periodic.pun``. The connectivity (denoted as ``connect_mode``) 
should be set to ionic, as MgO is an ionic system.

.. note:: Periodic inputs can also be read in from CIF files using the ChemShell/ASE interface, with an example given in ``mgo_cluster_cif.py``.

The cluster is then built using the ``construct_cluster()`` method:

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/construct_cluster.py
    :lines: 9-14

When this command is run ChemShell will first cut the cluster from the crystal,
which is controlled by the following options:

* ``origin`` dictates where the cluster is centred on, normally one of the 
  atoms at the surface (note that in Py-ChemShell atoms are counted from zero).
  The atoms in the resulting cluster will be ordered by increasing distance 
  from the origin.
* ``radius_cluster`` determines the overall size of the cluster in atomic units.

Next, the ionic charges on the outermost layer of the generated cluster are 
adjusted to account for the loss of their near neighbours when they were cut 
from the slab. As MgO is a straightforward two species system, the charges can 
simply be scaled by their coordination number, which is selected by the 
``adjust_charge`` option.

Finally, additional point charges are fitted around the outside of the cluster, 
in order to mimic as closely as possible the periodic electrostatic 
interactions that are lost when the cluster is cut out of the slab.
A number of options control this fit:

* ``radius_active`` sets the size of the electrostatic sampling region from 
  which potential and electric field information is taken for the fit. This 
  should always be set to be at least as large as the largest geometry 
  optimisation region you intend you use.
* ``bq_density`` relates to how many point charge are placed in the shell 
  of charges immediately surrounding the cluster. 
  This should be chosen to ensure a good fit, 
  but not so high that the point charges are overfitted (see below).
* ``bq_margin`` sets the distance of the shell of fitted point charges
  fitted point charges around the edge of the cluster. It is OK for the 
  point charges to be quite close as the outer edge of ions in the cluster 
  is not part of the active (potential sampling) region.
* For the surface case, ``bq_layer`` adds additional point charge layers above 
  the surface level, to ensure that the electrostatic environment is well 
  represented at the level at which reactions will take place.

There are no hard and fast rules for determining these settings, but the fit 
procedure will output a number of measures that indicate how well the 
resulting potential and fields match the original periodic structure. You 
should examine the 
fitted charges in the output (under the header "BQ charges") to check the fit 
has produced sensible results. 
The first six are located at a significant distance from the cluster in order 
to account for long range dipoles, and so can be large. 
The remainder should, in general, be reasonably small if the fit is of good 
quality.

A good strategy is to try the default values first, and then to check how 
sensitive the results are to changes in the input parameters.

Once the cluster has been constructed, the ``construct_cluster.py`` script 
will save it to a new punch file, ``mgo_shells_cluster.pun``, for use in the 
next stage of the process:

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/construct_cluster.py
    :lines: 16-17

Note that by default the point charges are denoted (arbitrarily) as **F**, so 
they can be understood by external MM software in subsequent calculations. For 
validation, compare against the pre-calculated example 
``mgo_shells_cluster.pun.REF``.


Partitioning the cluster
========================

Next, the cluster needs to be partitioned into regions. The ionic embedding 
model consists of five numbered regions corresponding to the QM region, QM/MM 
boundary, active MM, inactive MM, and fitted point charges. 

To proceed, it is 
required that the user identifies the **QM** region they wish to use. This can 
either be done manually by relabelling the relevant atoms in the punch file 
with a suffix of **"1"**, e.g. ``Mg`` becomes ``Mg1``, or by providing a list 
of atoms to the partitioning command using the option ``qm_region``.
In our example we take the latter approach, taking the 14 atoms closest to the 
origin as our QM region:

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/partition_cluster.py
   :lines: 13-15

The embedded cluster is then split into regions using the ``partition()`` 
method:

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/partition_cluster.py
   :lines: 17-20

Beyond specifying the QM region, the following options should be set to control 
how the other regions are selected:

* ``cutoff_boundary`` determines which ions are contained within the QM/MM 
  boundary region (region 2). Here we are using the so-called "explicit" 
  interface (specified by ``qmmm_interface='explicit'``), in which the 
  interface ions will be explicitly represented in the MM calculation 
  by charged particles. The interface ions will enter the QM calculation 
  as point charges with effective core potentials attached. The alternative 
  approach, known as an "implicit" interface, is not supported by NWChem and 
  so we do not consider it in this tutorial.
* ``interface_exclude`` is used to exclude certain ion types from the boundary 
  region. In our example we exclude oxygen ions, which then automatically go 
  into the active MM region (region 3). We do this because effective core 
  potentials for the oxygen anions are not readily available, and it is 
  sufficient to have only cations in the boundary region.
* ``radius_active`` defines the region that you want to optimise, i.e. the 
  active MM region (region 3). This can be the same size or smaller than the 
  sampling region used to fit charges around the cluster, but no bigger.

Finally the script will write a new punch file where each atom is labelled 
with a region number: 

.. literalinclude:: ../../samples/mgo_surface/1_cut_cluster/partition_cluster.py
   :lines: 22-23

As well as relabelling, the partitioning routine also deletes the superfluous
classical charges on the QM region atoms, as well as any corresponding shells,
and instead provides an overall charge to be used in QM calculations
(in this case -8).

.. note:: As MM partial charges are floating point numbers the calculated 
   QM charge will also be a floating point number, and should be rounded to 
   an integer before use in a QM calculation.