
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples_generated/tutorials/defects.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_generated_tutorials_defects.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_generated_tutorials_defects.py:

.. _defects:

==============================
Defect calculations: ASE Tools
==============================

This section gives an (incomplete) overview of features in ASE that
help in the preparation and analysis of supercell calculations as most
commonly employed in the computation of defect properties.

.. contents::

.. GENERATED FROM PYTHON SOURCE LINES 14-41

Supercell creation
==================

Background
----------

Defect properties are most commonly investigated in the so-called
dilute limit, i.e. under conditions, in which defect-defect
interactions are negligible. While alternative approaches in
particular embedding techniques exist, the most common approach is to
use supercells. To this end, one creates a supercell by a *suitable*
(see below) repetition of the primitive unit cell, after which a
defect, e.g., a vacancy or an impurity atom, is inserted.

The calculation thus corresponds to a periodic arrangement of
defects. Accordingly, care must be taken to keep the interactions
between defects as small as possible, which generally calls for large
supercells. Thus the typical goal for generating the simulation supercell
for defect calculations is to maximize the defect-defect separation in
*all* directions, for a reasonable number of atoms (and thus computational
cost). In principle, we can do a good job of this by using a supercell
with a suitable shape.
To illustrate this for different lattices, we build and plot
three 2D lattices with identical
unit cell area but
different lattice symmetry. We build and visualize (3,3,1) supercells
by repeating the unit cell.

.. GENERATED FROM PYTHON SOURCE LINES 41-162

.. code-block:: Python


    import matplotlib.pyplot as plt

    from ase.build import bulk
    from ase.geometry import get_distances
    from ase.visualize.plot import plot_atoms


    def add_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0, 0)):
        # function for adding arrows to neighboring atoms
        vectors, distances = get_distances(
            conf.positions[centralidx], p2=conf.positions
        )
        neighboridx = [
            i for i, j in enumerate(distances[0]) if (neighborcutoff >= j)
        ]

        # we normalize the size of the arrows for illustration purposes
        norm = neighborcutoff / 5
        for idx in neighboridx:
            ax.arrow(
                conf.positions[centralidx][0] + arrow_offset[0],
                conf.positions[centralidx][1] + arrow_offset[1],
                vectors[0][idx][0] / norm,
                vectors[0][idx][1] / norm,
                width=0.1,
                color='k',
            )


    lattice_constant = 8
    centralidx = 4
    neighborcutoff = 8


    # square lattice
    # build structure
    conf = bulk('Po', a=lattice_constant)
    # moving the atom in the middle of the cell
    positions = conf.get_positions()
    positions += conf.cell[0][0] / 2
    conf.set_positions(positions)
    conf = conf.repeat((3, 3, 1))

    # plot cs structure
    fig, ax = plt.subplots()
    add_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0.5, 0.5))
    plot_atoms(conf, ax, offset=(0, 0))  # , rotation=('-80x,0y,0z'))
    ax.set_axis_off()
    ax.text(
        0.1,
        -0.1,
        'square lattice: r$_1$=a, Z$_1$=4',
        transform=ax.transAxes,
        fontsize=16,
    )
    plt.show()


    # rectangular lattice
    # build structure
    conf = bulk(
        'C',
        crystalstructure='orthorhombic',
        a=lattice_constant,
        b=lattice_constant / 2,
        c=lattice_constant / 2,
    )

    # moving the atom in the middle of the cell
    positions = conf.get_positions()
    positions += [conf.cell[0][0] / 2, conf.cell[0][0] / 4, conf.cell[0][0] / 4]
    conf.set_positions(positions)
    conf = conf.repeat((3, 3, 1))

    # plot orc structure
    fig, ax = plt.subplots()
    add_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0.5, 0.5))
    plot_atoms(conf, ax, offset=(0, 0))
    ax.set_axis_off()
    ax.text(
        0,
        -0.2,
        'rectangular lattice with a 2:1 aspect ratio:\n r$_1$=a/2, Z$_1$=2',
        transform=ax.transAxes,
        fontsize=16,
    )
    plt.show()

    # hexagonal lattice
    # build structure
    conf = bulk('Be', a=lattice_constant)
    # here, we slice the cell to have one a 2D layer of atoms
    confmask = [i.index for i in conf if i.position[2] < 1]
    conf = conf[confmask]
    conf = conf.repeat((3, 3, 1))
    positions = conf.get_positions()
    positions -= conf.positions[0]
    conf.set_positions(positions)

    # plot hpc structure
    fig, ax = plt.subplots()
    add_decoration(
        ax,
        conf,
        centralidx,
        neighborcutoff,
        arrow_offset=(-conf.cell[1][0] + 1, 1.5),
    )
    plot_atoms(conf, ax, offset=(0, 0), rotation=('0x,0y,0z'))
    ax.set_axis_off()
    ax.text(
        0.1,
        -0.1,
        'hexagonal lattice: r$_1$=1.075a, Z$_1$=6',
        transform=ax.transAxes,
        fontsize=16,
    )
    plt.show()





.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_defects_001.png
         :alt: defects
         :srcset: /examples_generated/tutorials/images/sphx_glr_defects_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_defects_002.png
         :alt: defects
         :srcset: /examples_generated/tutorials/images/sphx_glr_defects_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples_generated/tutorials/images/sphx_glr_defects_003.png
         :alt: defects
         :srcset: /examples_generated/tutorials/images/sphx_glr_defects_003.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 163-270

In the case of the square lattice, each defect has :math:`Z_1=4`
nearest neighbors at a distance of :math:`r_1=a_0`, where
:math:`a_0=\sqrt{A}` with :math:`A` being the unit cell area. By
comparison in a rectangular lattice with an aspect ratio of 2:1, the
defects are much closer to each other with :math:`r_1 = a_0/\sqrt{2}` and
:math:`Z_1=2`, where again :math:`a_0` = :math:`\sqrt{A}`
(the 'effective cubic length').
The largest defect-defect distance (at constant unit
cell area) is obtained for the hexagonal lattice, which also
correponds to the most closely packed 2D arrangement. Here, one
obtains :math:`r_1=\sqrt{2}/\sqrt[4]{3}=1.075 a_0` and
:math:`Z_1=6`. For defect calculations, supercells corresponding to
hexagonal or square lattices have thus clear advantages. This argument
can be extended to 3D: Square lattices in 2D correspond to cubic
lattices (supercells) in 3D with :math:`r_1=a_0` and
:math:`Z_1=6`. The 3D analogue of the hexagonal 2D lattice are
hexagonal and cubic close packed structures (i.e. FCC, HCP), both of which
yield :math:`r_1 = \sqrt[6]{2} a_0 \approx 1.1225 a_0` and :math:`Z_1=12`.

It is straightforward to construct cubic or face-centered cubic (fcc,
cubic closed packed) supercells for cubic materials (including e.g,
diamond and zincblende) by using simple repetitions of the
conventional or primitive unit cells. For countless materials of lower
symmetry the choice of a supercell is, however not necessarily so
simple. The algorithm below represents a general solution to this
issue.

In the case of semiconductors and insulators with small dielectric
constants, defect-defect interactions are particularly pronounced due
to the weak screening of long-ranged electrostatic interactions. While
various correction schemes have been proposed, the most reliable
approach is still finite-size extrapolation using supercells of
different size. In this case care must be taken to use a sequence of
self-similar supercells in order for the extrapolation to be
meaningful. To motivate this statement consider that the leading
(monopole-monopole) term :math:`E_{mp}`, which scales with :math:`1/r`
and is proportional to the (ionic) dielectric constant
:math:`\epsilon_0`. The :math:`E_{mp}` term is geometry dependent and
in the case of simple lattices the dependence is easily expressed by
the Madelung constant. The geometry dependence implies that different
(super)cell shapes fall on different lines when plotting e.g., the
formation energy as a function of :math:`N^{-1/3}` (equivalent to an
effective inverse cell size, :math:`L^{-1} \propto N^{-1/3}`. For
extrapolation one should therefore only use geometrically equivalent
cells or at least cells that are as self-similar to each other as
possibly, see Fig. 10 in [Erhart]_ for a very clear example. In this
context there is therefore also a particular need for supercells
of a particular shape.


Algorithm for finding optimal supercell shapes
----------------------------------------------

The above considerations illustrate the need for a more systematic
approach to supercell construction. A simple scheme to construct
"optimal" supercells is described in [Erhart]_. Optimality here
implies that one identifies the supercell that for a given size
(number of atoms) most closely approximates the desired shape, most
commonly a simple cubic or fcc metric (see above). This approach
ensures that the defect separation is large and that the electrostatic
interactions exhibit a systematic scaling.

The ideal cubic cell metric for a given volume :math:`\Omega` is simply
given by :math:`\Omega^{1/3} \mathbf{I}`, which in general does not
satisfy the crystallographic boundary conditions. The :math:`l_2`-norm
provides a convenient measure of the deviation of any other cell
metric from a cubic shape. The optimality measure can thus be defined
as

.. math::
  \Delta_\text{sc}(\mathbf{h}) = ||\mathbf{h} - \Omega^{1/3} \mathbf{1}||_2,

Any cell metric that is compatible with the crystal symmetry can be
written in the form

.. math:: \mathbf{h} = \mathbf{P} \mathbf{h}_p

where :math:`\mathbf{P} \in \mathbb{Z}^{3\times3}` and
:math:`\mathbf{h}_p` is the primitive cell metric.  This approach can
be readily generalized to arbitrary target cell metrics. In order to
obtain a measure that is size-independent it is furthermore convenient
to introduce a normalization, which leads to the expression
implemented here, namely

.. math::
  \bar{\Delta}(\mathbf{Ph}_p) =
  ||Q\mathbf{Ph}_p - \mathbf{h}_\text{target}||_2,

where `Q = \left(\det\mathbf{h}_\text{target} \big/
\det\mathbf{h}_p\right)^{1/3}` is a normalization factor.  The
matrix :math:`\mathbf{P}_\text{opt}` that yields the optimal cell
shape for a given cell size can then be obtained by

.. math:: \mathbf{P}_\text{opt} =
  \underset{\mathbf{P}}{\operatorname{argmin}}
  \left\{ \bar\Delta\left(\mathbf{Ph}_p\right) | \det\mathbf{P}
  = N_{uc}\right\},

where :math:`N_{uc}` defines the size of the supercell in terms of the
number of primitive unit cells.


Implementation of algorithm
---------------------------

For illustration consider the following example. First we set up a
primitive face-centered cubic (fcc) unit cell and visualize it.

.. GENERATED FROM PYTHON SOURCE LINES 271-284

.. code-block:: Python


    import matplotlib.pyplot as plt

    from ase.build import bulk
    from ase.visualize.plot import plot_atoms

    conf = bulk('Au')


    fig, ax = plt.subplots()
    plot_atoms(conf, ax)
    ax.set_axis_off()




.. image-sg:: /examples_generated/tutorials/images/sphx_glr_defects_004.png
   :alt: defects
   :srcset: /examples_generated/tutorials/images/sphx_glr_defects_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 285-289

Then, we call
:func:`~ase.build.find_optimal_cell_shape` to obtain a
:math:`\mathbf{P}` matrix that will enable us to generate a supercell
with 32 atoms that is as close as possible to a simple cubic shape:

.. GENERATED FROM PYTHON SOURCE LINES 289-297

.. code-block:: Python

    import numpy as np

    from ase.build import find_optimal_cell_shape
    from ase.build.supercells import eval_length_deviation

    P1 = find_optimal_cell_shape(conf.cell, 32, 'sc')
    print(P1)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [[-2  2  2]
     [ 2 -2  2]
     [ 2  2 -2]]




.. GENERATED FROM PYTHON SOURCE LINES 298-310

More nicely rendered, this yields

.. math:: \mathbf{P}_1 =
  \left(\begin{array}{rrr} -2 & 2 & 2 \\ 2 & -2 & 2 \\
  2 & 2 & -2 \end{array}\right) \quad
  \mathbf{h}_1 = \left(\begin{array}{ccc} 2 a_0 & 0 & 0 \\
  0 & 2 a_0 & 0 \\ 0 & 0 & 2 a_0 \end{array}\right),

where :math:`a_0` =4.05 Å is the lattice constant. This is indeed the
expected outcome as it corresponds to a :math:`2\times2\times2`
repetition of the *conventional* (4-atom) unit cell. On the other hand
repeating this exercise with:

.. GENERATED FROM PYTHON SOURCE LINES 311-315

.. code-block:: Python


    P2 = find_optimal_cell_shape(conf.cell, 495, 'sc')
    print(P2)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    [[-5  5  5]
     [ 5 -4  5]
     [ 5  5 -4]]




.. GENERATED FROM PYTHON SOURCE LINES 316-328

yields a less obvious result, namely

.. math:: \mathbf{P}_2 =
  \left(\begin{array}{rrr} -6 & 5 & 5 \\ 5 & -6 & 5 \\
  5 & 5 & -5 \end{array}\right)
  \quad
  \mathbf{h}_2 = a_0
  \left(\begin{array}{ccc} 5 & 0 & 0 \\ 0.5 & 5 & 0.5 \\
   0.5 & 0.5 & 5 \end{array}\right),

which indeed corresponds to a reasonably cubic cell shape. One can
also obtain the optimality measure :math:`\bar{\Delta}` by executing:

.. GENERATED FROM PYTHON SOURCE LINES 328-333

.. code-block:: Python

    dev1 = eval_length_deviation(np.dot(P1, conf.cell))
    dev2 = eval_length_deviation(np.dot(P2, conf.cell))
    print(f'The length deviation for P_1 is {dev1}')
    print(f'The length deviation for P_2 is {dev2}')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    The length deviation for P_1 is 3.8459253727671276e-16
    The length deviation for P_2 is 0.019161199641064268




.. GENERATED FROM PYTHON SOURCE LINES 334-382

which yields :math:`\bar{\Delta}(\mathbf{P}_1)=0` and
:math:`\bar{\Delta}(\mathbf{P}_2)=0.0192`.

Since this procedure requires only knowledge of the cell metric (and
not the atomic positions) for standard metrics, e.g., fcc, bcc, and
simple cubic one can generate series of shapes that are usable for
*all* structures with the respective metric. For example the
:math:`\mathbf{P}_\text{opt}` matrices that optimize the shape of a
supercell build using a primitive FCC cell are directly applicable to
diamond and zincblende lattices.

For illustration, the :math:`\bar{\Delta}` values for supercells of SC, FCC
and BCC lattices with SC/FCC target shapes are shown as a function of
the number of unit cells :math:`N_{uc}\leq2000` in the panel below (taken
from :mr:`3404`). The algorithm is, however, most useful for
non-cubic cell shapes, for which finding several reasonably sized cell
shapes is more challenging, as illustrated for a hexagonal material
(LaBr\ :sub:`3`) in [Erhart]_.

.. image:: https://gitlab.com/-/project/470007/uploads/5c52f1b09cfd8f82c3b8453f45762d4f/image.png


.. note::
   For unit cells with more complex space groups,
   this approach can be cumbersome due
   to the implementation which loops over many possible
   transformation matrices. The
   `find_optimal_cell_shape <https://doped.readthedocs.io/en/latest/doped.utils.html
   doped.utils.supercells.find_optimal_cell_shape>`_
   function in `doped <https://doped.readthedocs.io>`_ implements
   the same algorithm with
   some efficiency improvements (~100x compute time speedup),
   and offers an efficient
   `algorithm <https://doped.readthedocs.io/en/latest/doped.utils.html
   doped.utils.supercells.find_ideal_supercell>`_
   for *directly* optimising the periodic defect-defect distance
   (~10-50% improvements);
   see [Kavanagh]_ or the ``doped`` `tutorials
   <https://doped.readthedocs.io/en/latest/generation_tutorial.html>`_.

Generation of supercell
-----------------------

Once the transformation matrix :math:`\mathbf{P}` it is
straightforward to generate the actual supercell using e.g., the
:func:`~ase.build.cut` function. A convenient interface is provided by
the :func:`~ase.build.make_supercell` function, which is invoked as
follows:

.. GENERATED FROM PYTHON SOURCE LINES 382-392

.. code-block:: Python


    from ase.build import make_supercell

    conf = bulk('Au')
    P = find_optimal_cell_shape(conf.cell, 495, 'sc')
    supercell = make_supercell(conf, P)

    fig, ax = plt.subplots()
    plot_atoms(supercell, ax)
    ax.set_axis_off()



.. image-sg:: /examples_generated/tutorials/images/sphx_glr_defects_005.png
   :alt: defects
   :srcset: /examples_generated/tutorials/images/sphx_glr_defects_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 393-403

.. [Erhart] P. Erhart, B. Sadigh, A. Schleife, and D. Åberg.
  First-principles study of codoping in lanthanum bromide,
  Phys. Rev. B, Vol **91**, 165206 (2012),
  :doi:`10.1103/PhysRevB.91.165206`; Appendix C

.. [Kavanagh] S. R. Kavanagh et al.
  doped: Python toolkit for robust and repeatable charged defect
  supercell calculations
  J. Open Source Softw, 9(**96**), 6433 (2024),
  :doi:`10.21105/joss.06433`


.. _sphx_glr_download_examples_generated_tutorials_defects.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: defects.ipynb <defects.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: defects.py <defects.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: defects.zip <defects.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
