.. _cutting_clusters:

****************
Cutting clusters
****************

The ChemShell QM/MM procedures are most commonly used with finite clusters,
especially for highly polar materials where cutoffs cannot be used.

Many of the systems to which QM/MM models are applied are condensed phase 
in nature, and the simulation of such materials is often made tractable by 
the use of periodic boundary conditions. However, periodic boundary 
conditions are problematic for QM/MM simulations as many QM codes  
do not contain the ability to compute the interaction 
of the QM system with an infinite charge distribution. We have implemented 
a charge fitting approach which enables a finite cluster to model accurately 
a single QM defect in an infinite, perfect classical crystalline system. 
In many cases, this is actually more pertinent to the real system under 
study than a model system with a periodically repeating defect.

The approach adopted is to start from the perfect, periodic system modelled 
using the classical force field approach (Figure a), and from it cut a 
cluster into which the QM model is to be embedded (b). The electrostatic 
potential in the central active region is computed (using the classical 
charge model) for both the perfect, infinite system and for the finite 
cluster (ie (a) and (b)). The difference in these potential datasets 
represents the error introduced by truncating the cluster. 

.. image:: images/cluster_cutting.png

It is often possible to reduce the magnitude of this error by adjusting the 
charges on the atoms on the surface of the finite cluster. This is 
particularly true of covalent solids where an adjustment can be applied 
for each bond cleaved by the cluster generation.

The remaining error in the electrostatic potential can be corrected by 
defining a set of additional point charges, placed around the outside of 
the cluster. The magnitudes of these charges are found by least-squares 
fitting to the error potential so that they compensate for the missing 
long-range electrostatic terms.

Finally the QM region can be defined (c), and with it a boundary layer 
if required.


Setting up clusters with construct_cluster
==========================================

.. py:function:: construct_cluster(...)

The :py:func:`construct_cluster` method is used to generate both spherical clusters
from periodic bulk systems and hemispherical surfaces from periodic slab
systems. The nature of the resulting cluster is determined automatically
by the periodicity of the input crystal.

:py:func:`construct_cluster` can be called as::

  my_crystal = Fragment(coords=...)

  my_cluster = my_crystal.construct_cluster(...)

where ``my_crystal`` is a periodic fragment and ``my_cluster`` is the fragment
containing the generated cluster.

The arguments to :py:func:`construct_cluster` are:

.. py:attribute:: crystal_type

   Allowed values:

   * ``'ionic'``: (default) Ionic type of bonding in the crystal

   * ``'covalent'``: covalent

   * ``'molecular'``: molecular

.. py:attribute:: origin_atom

   Atom to define the centre of the cluster (centre of top face for surfaces).

.. py:attribute:: origin_point

   Coordinates to define the centre of the cluster (centre of top face for
   surfaces).

.. py:attribute:: radius_cluster

   (Default: ``40.0``) Radius (in bohr) of the cut cluster.

.. py:attribute:: radius_active

   (Default: ``20.0``) Radius (in bohr) of the active MM region in the cluster.

.. py:attribute:: bq_margin

   (Default: ``3.5``) Distance (in bohr) between the cluster boundary and the outer point charges.

.. py:attribute:: bq_density

   (Default: ``5``) Density of added outer point charges. The number of added 
   point charges is :py:attr:`(4*bq_density**2)+2`

.. py:attribute:: bq_symbol

   (Default: ``'F'``) Arbitrary element symbol to assign to BQs for reference
   in MM codes.

.. py:attribute:: adjust_charge

   Specify a model for adjusting MM charges on the cluster boundary

   Allowed values:

   * ``''``: (default)

   * ``'coordination_scaled'``: 

   * ``'valence_model'``: 

.. py:attribute:: valence_matrix

   Charge adjustment matrix for use with the ``adjust_charge='valence_model'`` scheme.


.. TODO add further text from Tcl-ChemShell documentation when module is finished

.. figure:: images/radii_cluster_cutting.png

   Region 1 - QM;
   Region 2 - QM/MM interface;
   Region 3 - Active MM;
   Region 4 - Inactive MM;
   Region 5 - Background point charges
   
Setting up QM/MM based on cut cluster
=====================================

.. py:function:: partition(...)

:py:func:`construct_cluster` can be called as::

  my_cluster = my_crystal.construct_cluster(...)

  my_qmmm_cluster = my_cluster.partition(...)
  
.. py:attribute:: label_qm

   (Default: ``'1'``) Label attached to forcefield type name for the QM region (e.g., "O1").

.. py:attribute:: label_interface

   (Default: ``'2'``) Label to attach to forcefield type name for the QM/MM interface region (e.g., "O" to "O2").

.. py:attribute:: label_active

   (Default: ``'3'``) Label to attach to forcefield type name for the active MM region (e.g., "O" to "O3").

.. py:attribute:: label_inactive

   (Default: ``'4'``) Label to attach to forcefield type name for the active MM region (e.g., "O" to "O4").

.. py:attribute:: label_bq

   (Default: ``'5'``) Label to attach to forcefield type name for the BQs (e.g., "F" to "F5").

.. py:attribute:: bq_symbol

   (Default: ``'F'``) Arbitrary element symbol to assign to BQs for reference in MM codes.

.. py:attribute:: cutoff_boundary

   (Default: ``10.0``) Atoms within this distance (in bohr) from label_qm are assigned to label_interface.

.. py:attribute:: interface_exclude

   (Default: ``[]``) List of atom symbols to exclude from the interface region

.. py:attribute:: interface_rigid

   (Default: ``[]``) List of atom symbols to be rendered rigid in the interface region. For use with the ``qmmm_interface='explicit'`` only.

.. py:attribute:: origin

   (Default: ``[0.0,0.0,0.0]``) Coordinates of the cluster origin (output from construct_cluster) or atom index (integer) of origin

.. py:attribute:: qmmm_interface

   There are two definitions for the QM/MM interface region (region label_interface) that can be selected

   Allowed values:

   * ``'implicit'``: In the implicit interface, interface atoms are represented in the QM calculation as QM atoms with ECPs attached 

   * ``'explicit'``: In the explicit interface, the interface atoms are represented as BQs with ECPs attached in the QM calculation and full core-shell MM atoms in the MM calculation (if the shell model is being used)

.. py:attribute:: radius_active

   (Default: ``20.0``) 	All atoms within this radius (in bohr) from the origin are assigned to region 3 (active MM).