.. _li_doping:

****************************
Localising the electron hole
****************************

With the MgO cluster now constructed, the structure can now be modified
to produce the lithium doped cluster.
When MgO is doped with lithium, the imbalance of ionic charges which 
manifests from substituting Mg :sup:`2+` with Li :sup:`+` creates an 
extra unoccupied hole. Literature indicates that this hole is localized 
to solely one oxygen, as such the active site is a [Li :sup:`+` O :sup:`-` ] 
center (see `Wang-1986 <https://doi.org/10.1021/j100280a084>`_).
Therefore, to create a suitable model for the subsequent calculations, 
it is necessary to localise the electron hole on any given oxygen. 

To localise the hole, it is first necessary to generate a Li doped MgO cluster. 
The Li can then be displaced, and a geometry optimisation applied to find a
local minima corresponding to the applied displacement.
The geometry optimisation can be used in tandem with a Mulliken population analysis 
to determine the partial atomic charges on each oxygen and thus determine whether the electron hole 
is localised on one oxygen or not.

`1. Localising the electron hole`_
=========================================

.. note:: The scripts for this section of the tutorial can be found in ``li_mgo/2_li-mgo_opt/O7_localised`` .
  
To generate the Li doped MgO cluster, the first step is to read the previously constructed 
``ss-mgo_shells_regions.pun`` structure into the script ``Displace_Li.py``.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Displace_Li.py
   :lines: 5-6

To dope the cluster with lithium, the central Mg :sup:`2+` is substitued with Li :sup:`+`. 
The atom names and atomic numbers in the cluster can be redefined simply by changing the 
``names`` and ``znums`` labels.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Displace_Li.py
   :lines: 8-10

A small displacement of 0.4 bohr in the z axis is subsequently applied to the
central Li ion. This displacement is achieved using a numpy vector which takes 
the cartesian coordinates of Li, and adds the corresponding numpy vector. 

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Displace_Li.py
   :lines: 12-13

This structure is subsequently saved in Chemshell punch format for further manipulation.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Displace_Li.py
   :lines: 15-16


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

For the subsequent calculations, NWChem requires that the basis set is specified. 
The basis sets used for this tutorial are modified from the basis sets specified in the 
section :ref:`qm_mm_spoint` of the **'Modelling an MgO surface'** tutorial. 
The ``mgo_nwchem.basis`` and ``mgo_nwchem.ecp`` files were modified to include the basis set 
for Li. Similar to Mg, a Stuttgart relativisitic large core basis set and ecp are used. It is worth
noting that the most diffuse functions were removed to prevent artificial spreading of charge density 
into the active MM region.

The ``li-mgo_nwchem.basis`` and ``li-mgo_nwchem.ecp`` files used for all subsequent calculations in 
this tutorial are shown below:

**Basis sets**

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/li-mgo_nwchem.basis
   :lines: 1-43
   
**Effective core potentials**

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/li-mgo_nwchem.ecp
   :lines: 1-28


`3. Geometry optimisation`_
=========================================

The input structure used in this calculation is the ``Li_MgO-cluster.pun`` file created previously. This structure is 
initially read into the script ``Li_MgO-cluster_optimization.py``.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 6-7
   
The QM region and charge of the fragment is then read in using ``getRegion(1)`` and ``totalcharge`` respectively.
The active atoms which will be relaxed for the geometry optimisation is specified using ``getRegion``, specifying
regions 1-3 (the QM region, QM/MM interface and active MM region).

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 9-15

The theory for the QM/MM calculation can now be defined. The settings for the QM calculation use NWChem as the calculator,
the PBE0 hybrid functional, the charge of the QM region and the newly modified basis set and ecp files. The 
argument ``maxiter=500`` specifies for there to be up to 500 self-consistent field (SCF) cycles, instead of the default value
of 100 to help the calculation converge. A Mulliken analysis is then specified using the keyword ``pop_analysis='mulliken'``.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 18-20

.. caution:: The ``path`` argument should be removed if you are running ChemShell and NWChem directly-linked. This is
   the case for most of the high-performance computing (HPC) facilities, such as Hawk or ARCHER2. In which case,
   you should specify the requested resources in the bash submission script.
   
The MM settings use GULP as the calculator with the ``mgo_2body.ff`` forcefield defined in the :ref:`qm_mm_spoint` section
of **'Modelling an MgO surface'** tutorial.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 22

The QM/MM settings are then defined. This takes the ``qm_theory`` and ``mm_theory`` options defined above as arguments. The punch file of the
structure, the coupling model (ionic) and the shell relaxation cycles are then specified. The ``shell maxcycles`` is set to a small value of 5 to 
ensure timely completion of the calculation.

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 24-26
   
The geometry optimisation is carried out using the standard optimisation module in PyChemShell, DL-FIND, with the settings defined in the ``Opt`` method.
   
.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 28-34
   
The arguments specified in the ``Opt()`` method are defined below.

* ``theory`` determines the theory requested for the geometry optimisation (``qmmm`` in this case).

* ``algorithm`` specifies the optimisation algorithm, which in this case is the standard lbfgs algorithm.
 
* ``trust_radius`` specifies the trust radius approach (``'energy'`` bases this on the energy change).

* ``maxcycle`` determines the maximum number of optimisation cycles.

* ``maxene`` determines the maximum number of energy evaluations.

* ``active`` specifies the atoms in the geometry optimisation that will be relaxed. In this case, the ``active_region`` variable defined earlier is specified.

The geometry optimisation is then run using ``opt.run()`` and the optimised structure is saved as ``Li_MgO-cluster_optimized`` in both punch and xyz format for further 
manipulation. 

.. literalinclude:: ../../samples/li_mgo/2_li-mgo_opt/O7_localised/Li_MgO-cluster_optimization.py
   :lines: 36-38

`4. Characterisation of the electron hole location`_
=====================================================

The Mulliken analysis results are specified in the ``_nwchem.out`` output file under the heading **'Total Density - Mulliken Population Analysis'**. 
These results show the atom names, atomic numbers, charges and shell charges for each of the atoms in the QM region. The key attribute to note is the 
charges on the oxygens. You should find that the charge on the oxygen with atom index 7 is significantly lower than the other oxygens. This indicates 
that only oxygen 7 is electron deficient, and the electron hole is now localised.

As the partial atomic charges of the Mulliken analysis are sensitive on the choice of basis set, a more robust method of confirming the localisation
of the electron hole is to view the atomic orbital contributions to the molecular orbital which corresponds to the electron hole. If you look under the heading 
**'DFT Final Beta Molecular Orbital Analysis'** in the ``_nwchem.out`` output, you will find there are a number of headings. The **'Vector'** heading gives 
the molecular orbital numbers, **'Occ'** gives the orbital occupation and **'E'** gives the energy of the orbital.

The first 30 alpha and 29 beta orbitals are occupied with one electron. Therefore, the unoccupied beta spin LUMO (orbital 30) is where the electron hole is located. 
You should find that the atomic orbital contribution with the highest coefficient in beta spin orbital 30 is the pz function on oxygen 7. This result is consistent 
with the Mulliken analysis and indicates that the hole has been successfully localised.

.. image:: orbitals.png
   :width: 596px
   :height: 438px
   :scale: 50 %
   :alt: "orbital-occupation_diagram"
   :align: center

The beta spin 30 molecular orbital can also be plotted using `DPLOT <https://github.com/nwchemgit/nwchem/wiki/DPLOT>`_. First, the ``_nwchem.inp``, ``_nwchem.db``
and ``_nwchem.movecs`` input files are generated using a Py-ChemShell geometry optimisation as outlined earlier except the optimized structure is used as the input 
structure. These files can be found in the ``DPLOT`` subdirectory. To the ``_nwchem.inp`` file, the following is added beneath the ``end`` keyword of the ``dft`` 
section::

  dplot
    TITLE o30_beta  
    vectors _nwchem.movecs
     LimitXYZ
   -24.0 24.0 150
   -24.0 24.0 150
   -24.0 24.0 150
    spin beta
    orbitals view; 1; 30
    gaussian
    output o30_beta.cube
  end
  task dplot

The calculation can subsequently be run using NWChem. On a local machine, the command for this is :guilabel:`nwchem _nwchem.inp`. This generates a Gaussian cube file
(``o30_beta.cube``) which can be visualised using Jmol or other suitable open source molecular modelling packages. 

 .. image:: Li_MgO-cluster_beta-30.jpg
   :scale: 50 %
   :alt: "beta_orbital-30_cube-plot"
   :align: center

The electron hole can be localised on the other oxygens simply by changing the direction of the Li displacement before carrying out the geometry optimisation. The scripts for these
are contained in the other five subdirectories in the directory ``li_mgo/2_li-mgo_opt``. 

