Atomistic Corrective Scheme for Supercell Density Functional Theory Calculations of Charged Defects

A new method to correct formation energies of charged defects obtained by supercell density-functional calculations is presented and applied to bulk, surface, and low-dimensional systems. The method relies on atomistic models and a polarizable force field to describe a material system and its dielectric properties. The polarizable force field is based on a minimal set of fitting parameters, it accounts for the dielectric screening arising from ions and electrons separately, and it can be easily implemented in any software for atomistic molecular dynamics simulations. This work illustrates both technical aspects and applications of the new corrective scheme. The method is tested on systems in vacuo to validate the energy scheme. It is applied to charged defects in the bulk and at the surface of realistic materials to achieve comparison with published results obtained by using available corrective schemes based on continuum electrostatics treatments. Moreover, to demonstrate its generality, the method is used to correct the formation energy obtained by DFT of a singly negatively charged S vacancy in monolayer, bilayer, trilayer and bulk MoS2.

models for a defected dielectric system, and either analytical 3 or numerical 6 solutions of the Poisson equation. Although all valid, these schemes are limited either in scope to a selected class of defects and systems 1,3,9 or (in our view) by the stringent requirement of using a structureless jellium to represent the defected material system. Here, we present a new, alternative, and general method to correct formation energies of charge defects obtained by DFT. At variance with previous approaches, our method is based on the use of atomistic model structures of a defected material, a simple polarizable force field, and a self-consistent treatment of the dielectric screening arising separately from ions and electrons. It is to be noted that since our method is based on atomistic models and force fields, it can be easily implemented in any software for classical molecular dynamics simulations, and therefore it is readily accessible. In the following, we discuss technical aspects and applications of our method to calculate ΔE L .

Results and Discussion
Method description. Our method relies on the use of periodic atomistic model structures to represent a dielectric material. Each atomic site encompasses an ionic charge, Q, connected to an equilibrium position by a harmonic spring, K, and an electronic charge, q, treated either as a spherical shell or a Gaussian charge distribution of width σ. Ionic and electronic charges belonging to the same atomic site interact only via a harmonic spring, k, whereas inter-site electrostatic interactions are calculated by using the Ewald method 13 . In mathematical terms, the potential energy function, U, of a periodic supercell encompassing N ions and n spherical shells takes the following form:  where, in the equations above, V is the volume of the supercell, n and R n are indexes and corresponding lattice vector of the direct space of supercells, G is a vector of the reciprocal space,  q indicates either ionic (Q) or electronic (q) charges,  r refers to the vector position of either an ion (r) or an electronic charge (s), σ ew is the parameter controlling the convergence of the Ewald sums, Δr refers to the displacement of an ion from its equilibrium position, and Δs indicates the distance between ion and the respective electronic charge (see Fig. 1). In Eq. (4), S(G) is the structure factor of the ionic and electronic charges, the first sum excludes the interactions between ions and electronic charges belonging to same sites, and the last sum is used to eliminate these same interactions from the sum in the reciprocal space. In case of electronic charges treated as Gaussian charge distributions, Eq. (2) includes the following additional term: Figure 1. Left panel, ball and stick image of an atomistic model structure. In our scheme, each atomic site is described as shown in the 2D schematics enclosed in the dotted ellipse, i.e. an ionic charge, Q (red disc), and an electronic charge, q (wide blue disc). The ionic charge is connected to an equilibrium position via a harmonic spring K, and the electronic charge is connected to the ion via a harmonic spring k. The electronic charge can be treated either as a spherical shell or a Gaussian distribution of width σ. These parameters are calibrated to reproduce high-and low-frequency dielectric constants of the material system.  where the σ i 's correspond to the widths of the Gaussian distributions assigned to the electronic charge.
To determine equilibrium displacements of both ions and electronic charges (due to the presence of either a static electric field or a charge defect), we assign fictitious masses to both ionic and electronic charges and we use conventional damped molecular dynamics techniques 13 to determine the ground state energy and dielectric displacements of the supercell. Periodic model structures incorporating a charged defect are neutralized by using a uniform charge background, and to estimate the correction energy ΔE L = E ∞ − E L to be used in Eq. (1), the polarizable energy scheme in Eqs (2) and (5) is used to calculate E L and, by extrapolation, E ∞ . Unless otherwise specified, = L V 3 . The parameters {K, Q} and {k, q, σ} for each atomic site require to be calibrated to reproduce the high-and low-frequency dielectric constants of the material system. In case of homogeneous (periodic and simple aperiodic) systems, all cationic/anionic sites are equivalent, and thus our scheme requires to determine the values of only a few parameters. In this case, the fitting procedure involves three simple steps. First, selecting the value for the ionic and electronic charges, and optionally of the width of the Gaussian distributions; this task is easily accomplished by relying on chemical intuition and formal oxidation numbers. Second, using the Clausius-Mossotti relationships to obtain tentative values for the harmonic spring constants k and K. Third, carrying out a few calculations with the polarizable energy scheme to refine these tentative values and obtain the optimal spring constants yielding the desired values for both the high-and low-frequency dielectric constants of the material system of interest. In the general case of non-equivalent sites, such as for complex surfaces or interfaces, or multicomponent systems, the fitting procedure may, in principle, become involved. In this case, however, DFT calculations of dielectric-constant spatial profiles may facilitate the task, as it has been recently shown for a corrective scheme based on solving the Poisson equation for continuum models of inhomogeneous materials 6 .
Applications to model systems. We first apply our method to the following charge distributions in cubic supercells of vacuum: (i) a unit point charge, (ii) a unit charge distributed according to a Gaussian distribution of width 4 a 0 , and (iii) two Gaussian charge distributions of width 1 a 0 placed on the diagonal of the supercell at a distance of a 10 3 0 . In all cases, we use our scheme to calculate the energy of the array of charges interacting with a compensating uniform charge background.
The energy of a cubic array of point charges scales linearly with L −1 , according to the scaling function shown as inset in Fig. 2, with a Madelung constant derived by fitting our data in excellent agreement with published  10 . Our scheme yields also the correct behavior for the energy of an array of Gaussian charges (Fig. 2). In this case, in addition to terms scaling as L −1 , the energy depends also on terms scaling as L −3 , accounting for both the size and quadrupole moment of the charge distribution in a supercell ( Fig. 2) 1,10 . These results demonstrate that our method encompasses correction schemes relying on analytical formulas, applicable to charged defects of finite size in homogeneous bulk systems 1, 3 .
To illustrate key technical aspects of our method, we consider a charged defect in a fictitious material with a static dielectric permittivity of 4. The material is described by using cubic lattices with a spacing of 2 Å. Each site of the lattice encompasses an immobile ion carrying a positive unit charge, and an electronic spherical shell carrying a negative unit charge, connected to the ion via a harmonic spring k. To determine the value of k yielding the desired permittivity, we use the Clausius-Mossotti relation to estimate a tentative value (in this case k = 15.1 eV/ Å 2 ), and then we perform a few trial-and-error calculations to determine the optimal value k = 11.0 eV/Å 2 . We use this model dielectric material to calculate values of E ∞ − E L with supercells hosting a unit charge defect in the center of a cubic interstice of the lattice. Figure 3 shows that, in the case of a point charge defect, the energy scales linearly with L −1 , with a prefactor proportional to the Madelung constant of a cubic lattice and inversely proportional to the permittivity of the material. In the case of a defect charge distributed according to a Gaussian distribution of width 2 Å, our calculations show that in supercells of small to moderate size, the charge is only partially screened by the dielectric material, and thus that the linear scaling, or else the continuum limit, is approached only at large L (Fig. 3), when dielectric displacements are not constrained by periodic boundary conditions. These results demonstrate that, thanks to the atomistic representation and the self-consistent treatment of dielectric screening, our method accounts, at variance with previous schemes, for size effects associated with the screening of charged defects of finite size.
Applications to realistic systems. To show applications to realistic materials, we consider a Cl − vacancy defect in the bulk and at the (100) surface of NaCl 6, 14 , and a doubly negatively charged B N antisite defect in the bulk of the hexagonal form of BN 9 . Here we use our method to calculate the energy correction, E ∞ − E L , and achieve comparison with recent computational studies of these defects 6,9,14 . To represent NaCl, we use rock-salt-type lattices, whereas for h-BN, we use layered atomistic structures with the AA′ stacking pattern. In both cases, we use experimental lattice constants. With the parameters reported in Table 1, our energy scheme yields, in agreement with experimental and DFT results, a dielectric constant of 5.9 ( = .   Table 1. Anionic and cationic parameters used in our atomistic polarizable energy scheme to reproduce highand low-frequency dielectric constants of NaCl, h-BN, and monolayer MoS 2 . Harmonic spring constants, K and k, are expressed in eV/Å 2 , the Gaussian width, σ, in Å, and atomic and electronic charges, Q and q, in electron charge unit.
method yields an energy correction of about 0.8 eV (Fig. 4), in agreement with the result obtained by extrapolating a few energy values calculated by DFT using supercells with L up to a value of 5l 0 14 . Also, our method yields an energy correction of 0.7 eV for a Cl − vacancy defect at the (001) surface of rigid NaCl crystals (with = × × L l 2 2 4 0 3 ). This result is in agreement with the energy correction of about 0.6 eV obtained by solving the Poisson equation for continuum models of the defected surface 6 . Figure 4 shows the energy values obtained for a Cl − vacancy defect in both rigid and fully polarizable crystals, i.e. in which both ions and electronic shells participate in screening the charged defect. In both cases, and for both bulk and surface, at large L the energy scales linearly with L −1 , with a scaling function proportional to the Madelung constant of the lattice of supercells and inversely proportional to the static permittivity of the periodic system (Fig. 4). By fitting the energy values of the defected surface, we find, as expected, a Madelung constant equal to 2.275, and a dielectric constant equal to + ( 1)/2 0  , corresponding to the Madelung constant and permittivity of an array of tetragonal supercells, with c = 2a and encompassing equal parts of dielectric ( = . 5 9 0  ) and vacuum (Fig. 4).
A good agreement with a recent DFT study 9 is also obtained in the case of the − B N 2 defect in h-BN. In fact, our method gives E ∞ − E L = 0.85 eV to correct a formation energy obtained by DFT using a supercell of 8 × 8 × 3 unit cells, in close agreement with the energy value of about 0.75 eV obtained by Kumagai and Oba 9 . Also, we remark that, at variance with previous cases, the linear scaling function interpolating the energy values at large L (Fig. 4) is not related trivially to the symmetry and dielectric constant of the material. In case of anisotropic materials such as h-BN, the numerical approach is the only viable strategy to estimate E ∞ − E L 9 . Overall, these results demonstrate that our method is sound and applicable to a variety of defected systems, including aperiodic and anisotropic materials.
Application to a charged defect in a low-dimensional material. To demonstrate the generality of our method, we consider also the case of a charged defect in low-dimensional materials. In particular, we calculate the energy correction resulting from periodic DFT calculations of a singly negatively charged S vacancy in monolayer, bilayer, and tri-layer (and for completeness also in bulk) MoS 2 15, 16 . Our periodic DFT calculations 17 are carried out using normconserving pseudopotentials 18 for both Mo and S, a generalized gradient approximation for the exchange and correlation energy 19 , a plane-wave energy-cutoff of 80 Ry, and a semi-empirical corrective scheme to account for London dispersion interactions 20 . To describe the MoS 2 systems, we use orthorhombic supercells and atomic positions in accord to the 2H crystalline phase of MoS 2 with a AB layer stacking pattern. Single layers include 30 MoS 2 formula units, and 2D films are separated by vacuum regions of about 12 Å. We use DFT to calculate structural, dielectric, and electronic properties, as well as the formation energy of the neutral and singly negatively charged S vacancy defect. Lattice parameters of monolayer and bulk MoS 2 are found ~3% larger than experimental values, and in the case of bilayer and trilayer MoS 2 , the intra-layer spacing between S planes is found to be 3.17 Å, and the inter-layer separation between Mo planes is found to be 6.21 Å. Overall, structural, dielectric, and electronic properties computed in this study are in agreement with recent DFT study of the MoS 2 systems 16 . To calculate formation energies, we use half the energy of molecular S 2 and a Fermi level equal to the valence-band edge as chemical potentials for S and electrons, respectively. Relevant DFT results are reported in Table 2.
With the parameters in Table 1, our polarizable energy scheme gives, in agreement with our DFT calculations (Table 2), an in-plane and out-of-plane dielectric constant of 14.1 and 2.0, respectively, for monolayer MoS 2 , and values of 14.5 and 2.2 for bilayer MoS 2 , 14.7 and 2.5 for the tri-layer film, and 15.0 and 6.6 for the bulk phase. This polarizable energy scheme and atomistic model structures of the MoS 2 systems are thus used to estimate the energy corrections for the formation energies of the charged S vacancy defect obtained by DFT (Table 2). In these calculations, the defect is described by using a negative unit charge in the center of a triangular interstice formed by nearest neighboring Mo ions. Point-or Gaussian-like (with a width of 1 Å) distributions for the unit charge yield results differing by only 0.02 eV. Figure 5 shows the results obtained in the case of monolayer MoS 2 ; similar energy vs. 1/L curves are obtained also in the case of bilayer and trilayer MoS 2 . The correction energies obtained by extrapolation of the energy values to infinite volumes are shown in Table 2.
It has to be noted that due to the low dimensionality of the material, the defect energy in Fig. 5 scales non-trivially with L −1 . In the present case, we find that extrapolation of the energy values to large volumes is aided  Table 2. Results obtained from DFT calculations for the dielectric constants (perpendicular, ⊥  , and parallel,  , to the layers), work function (Φ), and formation energies of the neutral and singly negatively S vacancy in monolayer (1L), bilayer (2L), tri-layer (3L), and bulk MoS 2 . The last column reports the energy corrections for the formation energies of the charged defects obtained by using the method presented in this work. Energy values are in eV. , where 0  is the component of the dielectric tensor parallel to the monolayer. Discs show energy values obtained by using our energy scheme, and the dotted line shows the linear interpolation of energy values at large L obtained by using supercells with l z = l x . Energy values are referred to that one obtained with the same supercell used to calculate the defect formation energy by DFT. The inset shows a schematics of the defected 2D material. by introducing a rescaled linear dimension L * , related to L, thickness and in-plane dielectric constant of the 2D MoS 2 systems as described in the caption of Fig. 5. In agreement with recent DFT studies 8,11,12,16 , these results demonstrate that the task of calculating the formation energy of a charged defect in a low-dimensional material cannot be accomplished by using only DFT. A numerical scheme needs to be used to estimate ΔE L = E ∞ − E L , which in the case of a singly negatively charged S vacancy in monolayer MoS 2 amounts to about 10% of the energy value obtained by DFT (Table 2).

Conclusion
We have introduced and applied a novel general method to correct the energy of charged defects obtained by DFT. At variance with previous methods, our approach is based on a polarizable energy scheme, atomistic representations, and a self-consistent treatment of the dielectric screening. The force field is based on a minimal set of fitting parameters, which can be easily calibrated by relying on either experimental data or, in case of more complex cases such as aperiodic or multicomponent systems, polarizability profiles derived from DFT calculations 6 . Being based on atomistic models, force fields, and Ewald sums, our method can be easily implemented in any available software for molecular dynamics simulations of materials or biological systems. Therefore, besides being new and general, our method is readily accessible to the community interested in DFT calculations and charged defects.