.. _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 `_). 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 `_. 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``.