defects_theory.rst 9.61 KB
 Sten Haastrup committed Feb 04, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 .. _defects_theory: ======================================================================= Localised electrostatic charges in non-uniform dielectric media: Theory ======================================================================= For examples, see :ref:defects_tutorial. Introduction ============ The purpose of this section is to enable us to calculate the formation energy of charged defects. In the Zhang-Northrup formula, this is given by .. math:: E^f[X^q] = E[X^q] - E_0 - \sum_i\mu_in_i + q (\epsilon_v + \epsilon_F) In this formula, X labels the type of defect and q its charge state, i.e. the net charge contained in some volume surrounding the defect. q is defined such that q=-1 for an electron. E[X^q] is the total energy of the sample with the defect, and E_0 the energy of the pristine (bulklike) sample. These quantities are usually calculated in a supercell approach with periodic boundary conditions. That is, we create the defect we are interested in, and place it in an environment containing many repetitions of the pristine unit cell. In the infinite supercell limit, we can describe the properties of the isolated defect. When we employ periodic boundary conditions, our calculation includes spurious electrostatic interactions between localised charge distribution of the defect state, and all its periodically repeated images. These long-ranged interactions mean that the convergence with respect to supercell size is slow. To accelerate this convergence, we employ an electrostatic correction scheme, and write .. math:: E^f[X^q]_{\mathrm{corrected}} = E^f[X^q]_{\mathrm{uncorrected}} - E_{\mathrm{periodic}} + E_{\mathrm{isolated}} + q\Delta V This assumes that DFT describes the bonding and energies of the system well, but contains errors in the electrostatics. The correction thus consists in subtracting the spurious interactions between the periodic images (the term included in the DFT calculation) and adding in the energy of an isolated charge distribution in the given dielectric environment. Finally, the potentials between the charged and neutral states must be aligned to the same reference. Implementing this scheme consists of the following steps, which will be described in further detail in the sections below. First, we determine the z-dependent dielectric function of the 2D layers. Knowing this, we use a model charge distribution to emulate the behaviour of the defect charge state, and find the potential associated with this model distribution by solving the Poisson equation. We do this twice: first with periodic boundary conditions, and then with zero boundary conditions. The electrostatic energy associated with the charge distribution in a given potential can then be found from the usual formula, .. math:: U = \frac{1}{2}\int_{\Omega} \rho V. Defining the dielectric response of 2D layers ============================================= In two dimensions, the bulk dielectric response is poorly defined, and some model must be used. The approach used here, following Ref [#Komsa]_ is to assume that the dielectric function of the isolated layer is isotropic in plane, and varies only in the z direction. Additionally, the screening is assumed to follow the density distribution of the system so that we can write .. math:: \varepsilon^{i}(z) = k^i\cdot n(z) + 1, Where i varies over "in-plane" and "out-of-plane", and n is the in-plane averaged density of the system. The normalization constants, k^i, are chosen such that .. math:: \frac{1}{L} \int \mathrm{d} z\, \varepsilon^{\parallel}(z) &= \varepsilon^{\parallel}_{\mathrm{DFT}} \\ \frac{1}{L} \int \mathrm{d} z\, \left(\varepsilon^{\perp}(z)\right)^{-1} &= \left(\varepsilon^{\perp}_{\mathrm{DFT}}\right)^{-1} Calculating the energy of the periodic images ============================================= Once we have the dielectric response, we proceed by solving the poisson equation .. math:: \nabla \cdot \boldsymbol{\varepsilon}(z) \odot \nabla \phi(\mathbf r) = -\rho(\mathbf r). Here \odot is the elementwise (Hadamard) product, and \boldsymbol{\varepsilon} = (\varepsilon^{\parallel}, \varepsilon^{\parallel}, \varepsilon^{\perp}). \rho(\mathbf r) is a model charge distribution that we use to describe the charge distribution of the defect state. For convenience, a Gaussian is often chosen, so that .. math:: \rho(\mathbf r) = \frac{1}{\left(\sqrt{2\pi}\sigma\right)^3} e^{-(\mathbf r - \mathbf r_0)^2/(2\sigma^2)}. In terms of the coordinate axes, the poisson equation reduces to .. math:: \varepsilon^{\parallel}(z)\frac{\partial^2 V}{\partial x^2} + \varepsilon^{\parallel}(z)\frac{\partial^2 V}{\partial y^2} + \varepsilon^{\perp}(z)\frac{\partial^2 V}{\partial z^2} + \frac{\partial \varepsilon^{\perp}}{\partial z} \frac{\partial V}{\partial z} = -\rho(\mathbf r) Since we wish to solve this equation with periodic boundary conditions, we Fourier transform the above equation, giving .. math:: \varepsilon^{\parallel}(G_z) * \left[(G_x^2 + G_y^2)V(\mathbf G)\right] + \varepsilon^{\perp}(G_z) * \left[G_z^2 V(\mathbf G)\right] + \left[G_z \varepsilon^{\perp}(G_z)\right] * \left[ G_z V\mathbf(G)\right] = \rho(\mathbf G), where * denotes a convolution along the z axis. Writing out these convolutions, we finally arrive at the expression .. math:: \rho_{G_x,G_y,G_z} = \sum_{G_z'} \left[\varepsilon^\parallel_{G_z - G_z'}\left(G_x^2 + G_y^2\right) + \varepsilon^{\perp}_{G_z - G_z'}G_zG_z'\right] V_{G_x, G_y, G_z'}. For each value of (G_x, G_y), we can thus calculate the corresponding potential V_{G_z} through a matrix inversion, and use that to calculate the energy of the model charge distribution. We can also use the potential V_{G_z} to calculate the alignment term, \Delta V. We can Fourier transform this to get real-space potential of the model charge distribution. If we have described the electrostatics of the system well, this potential should be similar to the true potential of the defect charge distribution, up to a constant shift. Defining .. math:: \Delta V(\vec{r}) = V(\vec{r}) - [V^{X^q}_\mathrm{el}(\vec{r}) - V^{0}_\mathrm{el}(\vec{r}) ], We set Calculating the energy of the isolated system ============================================= We start as before, with the Poisson equation, but since we would like to describe the energy of the isolated defect, we do not impose periodic boundary conditions and Fourier transform. Instead, following Ref. [#Ping]_ we can exploit the in-plane symmetry of the problem and expand \phi using cylindrical Bessel functions. .. math:: \phi(\mathbf r) = \int_0^\infty \mathrm{d}k'\, 2qe^{-k'^2\sigma^2/2} \varphi_{k'}(z) J_0(\rho k')\right Inserting this into the above equation and using the orthogonality relation \int \rho\mathrm{d}\rho J_0(\rho k)J_0(\rho k') = \delta(k - k') / k we find that \varphi_k must obey the Poisson equation .. math:: -\frac{\partial}{\partial z}\left(\varepsilon^{\perp}(z) \frac{\partial \varphi_k(z)}{\partial z}\right) + k^2\varepsilon^{\parallel}(z) \varphi_k(z) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)}, where z_0 is the center of the gaussian density along the z direction. The normalization of \varphi_k defined above was chosen precisely so that the right hand side of this equation is a normalized gaussian along the z direction. We solve this equation by separating the response into two components: The bulk response, describing the screening far away from the material, and the remaining z -dependent response close to the system. We thus define \Delta \varepsilon^i(z) = \varepsilon^i(z) - \varepsilon^{i}_{\mathrm{bulk}} and the Green's function of the bulk response \hat K = (-\varepsilon^{\perp}_{\mathrm{bulk}} \frac{\partial^2}{\partial z^2} + k^2\varepsilon^{\parallel}_{\mathrm{bulk}})^{-1}. As an implementation detail, we note that for 2D materials, the bulk response is generally 1. The equation for \varphi_k can then be written as .. math:: \hat{K}^{-1} \varphi_k - \frac{\partial}{\partial z}\left(\Delta\varepsilon^{\perp} \frac{\partial \varphi_k}{\partial z}\right) + k^2\Delta\varepsilon^{\parallel}\varphi_k = \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)}. Only the first term on the left hand side is affected by the boundary conditions on \varphi_k. We can solve this by Fourier transforming along the z axis and wigner-seitz truncating the Green's function, which yields the following equation .. math:: \sum_{G_z} D_{G_zG_z'} \left(\varphi_k\right)_{G_z} = e^{-i G_z'z_0 - G_z'^2\sigma^2/2}, with the matrix D given by .. math:: \frac{1}{L}D_{G_z'G_z} = \frac{\varepsilon^{\parallel}_{\mathrm{b}}k^2 + \varepsilon^{\perp}_{\mathrm{b}}G_z^2}{1 - e^{-kL/2}\cos(G_zL/2)}\delta_{G_zG_z'} + \Delta\varepsilon^{\parallel}_{G_z - G_z'}k^2 + \varepsilon^{\perp}_{G_z - G_z'}G_zG_z'. Finding \varphi_k is thus just a simple matrix inversion. Once we have solved the poisson equation, we calculate the total energy. .. math:: U &= \frac{1}{2} \int_{\Omega} \rho(\mathbf r) \phi(\mathbf r) \\ &= q^2 \int k \mathrm{d}k e^{-k^2 \sigma^2} U_k, with .. math:: U_k \int \mathrm{d}z\, \varphi_k(z) \frac{1}{\sqrt{2\pi}\sigma}e^{-\left(z - z_0\right)^2/\left(2\sigma^2\right)} Using the solution to the Poisson equation, this reduces to .. math:: U_k = \sum_{G_z,G_z'} e^{i(G_z - G_z')z_0 - (G_z^2 + G_z'^2)\sigma^2 / 2} \left(D^{-1}\right)_{G_zG_z'}, With D defined as above. References ========== .. [#Komsa] H.-P. Komsa, T. T. Rantala and A. Pasquarello *Phys. Rev. B* **86**, 045112 (2012) .. [#Ping] R. Sundararaman and Y. Ping *J. Chem. Phys.* **146**, 104109 (2017)