{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Defect calculations: ASE Tools\n\nThis section gives an (incomplete) overview of features in ASE that\nhelp in the preparation and analysis of supercell calculations as most\ncommonly employed in the computation of defect properties.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Supercell creation\n\n### Background\n\nDefect properties are most commonly investigated in the so-called\ndilute limit, i.e. under conditions, in which defect-defect\ninteractions are negligible. While alternative approaches in\nparticular embedding techniques exist, the most common approach is to\nuse supercells. To this end, one creates a supercell by a *suitable*\n(see below) repetition of the primitive unit cell, after which a\ndefect, e.g., a vacancy or an impurity atom, is inserted.\n\nThe calculation thus corresponds to a periodic arrangement of\ndefects. Accordingly, care must be taken to keep the interactions\nbetween defects as small as possible, which generally calls for large\nsupercells. Thus the typical goal for generating the simulation supercell\nfor defect calculations is to maximize the defect-defect separation in\n*all* directions, for a reasonable number of atoms (and thus computational\ncost). In principle, we can do a good job of this by using a supercell\nwith a suitable shape.\nTo illustrate this for different lattices, we build and plot\nthree 2D lattices with identical\nunit cell area but\ndifferent lattice symmetry. We build and visualize (3,3,1) supercells\nby repeating the unit cell.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n\nfrom ase.build import bulk\nfrom ase.geometry import get_distances\nfrom ase.visualize.plot import plot_atoms\n\n\ndef add_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0, 0)):\n # function for adding arrows to neighboring atoms\n vectors, distances = get_distances(\n conf.positions[centralidx], p2=conf.positions\n )\n neighboridx = [\n i for i, j in enumerate(distances[0]) if (neighborcutoff >= j)\n ]\n\n # we normalize the size of the arrows for illustration purposes\n norm = neighborcutoff / 5\n for idx in neighboridx:\n ax.arrow(\n conf.positions[centralidx][0] + arrow_offset[0],\n conf.positions[centralidx][1] + arrow_offset[1],\n vectors[0][idx][0] / norm,\n vectors[0][idx][1] / norm,\n width=0.1,\n color='k',\n )\n\n\nlattice_constant = 8\ncentralidx = 4\nneighborcutoff = 8\n\n\n# square lattice\n# build structure\nconf = bulk('Po', a=lattice_constant)\n# moving the atom in the middle of the cell\npositions = conf.get_positions()\npositions += conf.cell[0][0] / 2\nconf.set_positions(positions)\nconf = conf.repeat((3, 3, 1))\n\n# plot cs structure\nfig, ax = plt.subplots()\nadd_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0.5, 0.5))\nplot_atoms(conf, ax, offset=(0, 0)) # , rotation=('-80x,0y,0z'))\nax.set_axis_off()\nax.text(\n 0.1,\n -0.1,\n 'square lattice: r$_1$=a, Z$_1$=4',\n transform=ax.transAxes,\n fontsize=16,\n)\nplt.show()\n\n\n# rectangular lattice\n# build structure\nconf = bulk(\n 'C',\n crystalstructure='orthorhombic',\n a=lattice_constant,\n b=lattice_constant / 2,\n c=lattice_constant / 2,\n)\n\n# moving the atom in the middle of the cell\npositions = conf.get_positions()\npositions += [conf.cell[0][0] / 2, conf.cell[0][0] / 4, conf.cell[0][0] / 4]\nconf.set_positions(positions)\nconf = conf.repeat((3, 3, 1))\n\n# plot orc structure\nfig, ax = plt.subplots()\nadd_decoration(ax, conf, centralidx, neighborcutoff, arrow_offset=(0.5, 0.5))\nplot_atoms(conf, ax, offset=(0, 0))\nax.set_axis_off()\nax.text(\n 0,\n -0.2,\n 'rectangular lattice with a 2:1 aspect ratio:\\n r$_1$=a/2, Z$_1$=2',\n transform=ax.transAxes,\n fontsize=16,\n)\nplt.show()\n\n# hexagonal lattice\n# build structure\nconf = bulk('Be', a=lattice_constant)\n# here, we slice the cell to have one a 2D layer of atoms\nconfmask = [i.index for i in conf if i.position[2] < 1]\nconf = conf[confmask]\nconf = conf.repeat((3, 3, 1))\npositions = conf.get_positions()\npositions -= conf.positions[0]\nconf.set_positions(positions)\n\n# plot hpc structure\nfig, ax = plt.subplots()\nadd_decoration(\n ax,\n conf,\n centralidx,\n neighborcutoff,\n arrow_offset=(-conf.cell[1][0] + 1, 1.5),\n)\nplot_atoms(conf, ax, offset=(0, 0), rotation=('0x,0y,0z'))\nax.set_axis_off()\nax.text(\n 0.1,\n -0.1,\n 'hexagonal lattice: r$_1$=1.075a, Z$_1$=6',\n transform=ax.transAxes,\n fontsize=16,\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the case of the square lattice, each defect has $Z_1=4$\nnearest neighbors at a distance of $r_1=a_0$, where\n$a_0=\\sqrt{A}$ with $A$ being the unit cell area. By\ncomparison in a rectangular lattice with an aspect ratio of 2:1, the\ndefects are much closer to each other with $r_1 = a_0/\\sqrt{2}$ and\n$Z_1=2$, where again $a_0$ = $\\sqrt{A}$\n(the 'effective cubic length').\nThe largest defect-defect distance (at constant unit\ncell area) is obtained for the hexagonal lattice, which also\ncorreponds to the most closely packed 2D arrangement. Here, one\nobtains $r_1=\\sqrt{2}/\\sqrt[4]{3}=1.075 a_0$ and\n$Z_1=6$. For defect calculations, supercells corresponding to\nhexagonal or square lattices have thus clear advantages. This argument\ncan be extended to 3D: Square lattices in 2D correspond to cubic\nlattices (supercells) in 3D with $r_1=a_0$ and\n$Z_1=6$. The 3D analogue of the hexagonal 2D lattice are\nhexagonal and cubic close packed structures (i.e. FCC, HCP), both of which\nyield $r_1 = \\sqrt[6]{2} a_0 \\approx 1.1225 a_0$ and $Z_1=12$.\n\nIt is straightforward to construct cubic or face-centered cubic (fcc,\ncubic closed packed) supercells for cubic materials (including e.g,\ndiamond and zincblende) by using simple repetitions of the\nconventional or primitive unit cells. For countless materials of lower\nsymmetry the choice of a supercell is, however not necessarily so\nsimple. The algorithm below represents a general solution to this\nissue.\n\nIn the case of semiconductors and insulators with small dielectric\nconstants, defect-defect interactions are particularly pronounced due\nto the weak screening of long-ranged electrostatic interactions. While\nvarious correction schemes have been proposed, the most reliable\napproach is still finite-size extrapolation using supercells of\ndifferent size. In this case care must be taken to use a sequence of\nself-similar supercells in order for the extrapolation to be\nmeaningful. To motivate this statement consider that the leading\n(monopole-monopole) term $E_{mp}$, which scales with $1/r$\nand is proportional to the (ionic) dielectric constant\n$\\epsilon_0$. The $E_{mp}$ term is geometry dependent and\nin the case of simple lattices the dependence is easily expressed by\nthe Madelung constant. The geometry dependence implies that different\n(super)cell shapes fall on different lines when plotting e.g., the\nformation energy as a function of $N^{-1/3}$ (equivalent to an\neffective inverse cell size, $L^{-1} \\propto N^{-1/3}$. For\nextrapolation one should therefore only use geometrically equivalent\ncells or at least cells that are as self-similar to each other as\npossibly, see Fig. 10 in [Erhart]_ for a very clear example. In this\ncontext there is therefore also a particular need for supercells\nof a particular shape.\n\n\n### Algorithm for finding optimal supercell shapes\n\nThe above considerations illustrate the need for a more systematic\napproach to supercell construction. A simple scheme to construct\n\"optimal\" supercells is described in [Erhart]_. Optimality here\nimplies that one identifies the supercell that for a given size\n(number of atoms) most closely approximates the desired shape, most\ncommonly a simple cubic or fcc metric (see above). This approach\nensures that the defect separation is large and that the electrostatic\ninteractions exhibit a systematic scaling.\n\nThe ideal cubic cell metric for a given volume $\\Omega$ is simply\ngiven by $\\Omega^{1/3} \\mathbf{I}$, which in general does not\nsatisfy the crystallographic boundary conditions. The $l_2$-norm\nprovides a convenient measure of the deviation of any other cell\nmetric from a cubic shape. The optimality measure can thus be defined\nas\n\n\\begin{align}\\Delta_\\text{sc}(\\mathbf{h}) = ||\\mathbf{h} - \\Omega^{1/3} \\mathbf{1}||_2,\\end{align}\n\nAny cell metric that is compatible with the crystal symmetry can be\nwritten in the form\n\n\\begin{align}\\mathbf{h} = \\mathbf{P} \\mathbf{h}_p\\end{align}\n\nwhere $\\mathbf{P} \\in \\mathbb{Z}^{3\\times3}$ and\n$\\mathbf{h}_p$ is the primitive cell metric. This approach can\nbe readily generalized to arbitrary target cell metrics. In order to\nobtain a measure that is size-independent it is furthermore convenient\nto introduce a normalization, which leads to the expression\nimplemented here, namely\n\n\\begin{align}\\bar{\\Delta}(\\mathbf{Ph}_p) =\n ||Q\\mathbf{Ph}_p - \\mathbf{h}_\\text{target}||_2,\\end{align}\n\nwhere `Q = \\left(\\det\\mathbf{h}_\\text{target} \\big/\n\\det\\mathbf{h}_p\\right)^{1/3}` is a normalization factor. The\nmatrix $\\mathbf{P}_\\text{opt}$ that yields the optimal cell\nshape for a given cell size can then be obtained by\n\n\\begin{align}\\mathbf{P}_\\text{opt} =\n \\underset{\\mathbf{P}}{\\operatorname{argmin}}\n \\left\\{ \\bar\\Delta\\left(\\mathbf{Ph}_p\\right) | \\det\\mathbf{P}\n = N_{uc}\\right\\},\\end{align}\n\nwhere $N_{uc}$ defines the size of the supercell in terms of the\nnumber of primitive unit cells.\n\n\n### Implementation of algorithm\n\nFor illustration consider the following example. First we set up a\nprimitive face-centered cubic (fcc) unit cell and visualize it.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n\nfrom ase.build import bulk\nfrom ase.visualize.plot import plot_atoms\n\nconf = bulk('Au')\n\n\nfig, ax = plt.subplots()\nplot_atoms(conf, ax)\nax.set_axis_off()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then, we call\n:func:`~ase.build.find_optimal_cell_shape` to obtain a\n$\\mathbf{P}$ matrix that will enable us to generate a supercell\nwith 32 atoms that is as close as possible to a simple cubic shape:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n\nfrom ase.build import find_optimal_cell_shape\nfrom ase.build.supercells import eval_length_deviation\n\nP1 = find_optimal_cell_shape(conf.cell, 32, 'sc')\nprint(P1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More nicely rendered, this yields\n\n\\begin{align}\\mathbf{P}_1 =\n \\left(\\begin{array}{rrr} -2 & 2 & 2 \\\\ 2 & -2 & 2 \\\\\n 2 & 2 & -2 \\end{array}\\right) \\quad\n \\mathbf{h}_1 = \\left(\\begin{array}{ccc} 2 a_0 & 0 & 0 \\\\\n 0 & 2 a_0 & 0 \\\\ 0 & 0 & 2 a_0 \\end{array}\\right),\\end{align}\n\nwhere $a_0$ =4.05 \u00c5 is the lattice constant. This is indeed the\nexpected outcome as it corresponds to a $2\\times2\\times2$\nrepetition of the *conventional* (4-atom) unit cell. On the other hand\nrepeating this exercise with:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "P2 = find_optimal_cell_shape(conf.cell, 495, 'sc')\nprint(P2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "yields a less obvious result, namely\n\n\\begin{align}\\mathbf{P}_2 =\n \\left(\\begin{array}{rrr} -6 & 5 & 5 \\\\ 5 & -6 & 5 \\\\\n 5 & 5 & -5 \\end{array}\\right)\n \\quad\n \\mathbf{h}_2 = a_0\n \\left(\\begin{array}{ccc} 5 & 0 & 0 \\\\ 0.5 & 5 & 0.5 \\\\\n 0.5 & 0.5 & 5 \\end{array}\\right),\\end{align}\n\nwhich indeed corresponds to a reasonably cubic cell shape. One can\nalso obtain the optimality measure $\\bar{\\Delta}$ by executing:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "dev1 = eval_length_deviation(np.dot(P1, conf.cell))\ndev2 = eval_length_deviation(np.dot(P2, conf.cell))\nprint(f'The length deviation for P_1 is {dev1}')\nprint(f'The length deviation for P_2 is {dev2}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which yields $\\bar{\\Delta}(\\mathbf{P}_1)=0$ and\n$\\bar{\\Delta}(\\mathbf{P}_2)=0.0192$.\n\nSince this procedure requires only knowledge of the cell metric (and\nnot the atomic positions) for standard metrics, e.g., fcc, bcc, and\nsimple cubic one can generate series of shapes that are usable for\n*all* structures with the respective metric. For example the\n$\\mathbf{P}_\\text{opt}$ matrices that optimize the shape of a\nsupercell build using a primitive FCC cell are directly applicable to\ndiamond and zincblende lattices.\n\nFor illustration, the $\\bar{\\Delta}$ values for supercells of SC, FCC\nand BCC lattices with SC/FCC target shapes are shown as a function of\nthe number of unit cells $N_{uc}\\leq2000$ in the panel below (taken\nfrom :mr:`3404`). The algorithm is, however, most useful for\nnon-cubic cell shapes, for which finding several reasonably sized cell\nshapes is more challenging, as illustrated for a hexagonal material\n(LaBr\\ :sub:`3`) in [Erhart]_.\n\n\n\n\n

Note

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

\n\n### Generation of supercell\n\nOnce the transformation matrix $\\mathbf{P}$ it is\nstraightforward to generate the actual supercell using e.g., the\n:func:`~ase.build.cut` function. A convenient interface is provided by\nthe :func:`~ase.build.make_supercell` function, which is invoked as\nfollows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from ase.build import make_supercell\n\nconf = bulk('Au')\nP = find_optimal_cell_shape(conf.cell, 495, 'sc')\nsupercell = make_supercell(conf, P)\n\nfig, ax = plt.subplots()\nplot_atoms(supercell, ax)\nax.set_axis_off()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. [Erhart] P. Erhart, B. Sadigh, A. Schleife, and D. \u00c5berg.\n First-principles study of codoping in lanthanum bromide,\n Phys. Rev. B, Vol **91**, 165206 (2012),\n :doi:`10.1103/PhysRevB.91.165206`; Appendix C\n\n.. [Kavanagh] S. R. Kavanagh et al.\n doped: Python toolkit for robust and repeatable charged defect\n supercell calculations\n J. Open Source Softw, 9(**96**), 6433 (2024),\n :doi:`10.21105/joss.06433`\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.7" } }, "nbformat": 4, "nbformat_minor": 0 }