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

Atomistic Corrective Scheme for Supercell Density Functional Theory Calculations of Charged Defects www.nature.com/scientificreports OPEN Atomistic Corrective Scheme for Supercell Density Functional Theory Calculations of Charged Received: 20 March 2017 Defects Accepted: 20 April 2017 Published: xx xx xxxx 1 1,2,3 Tengfei Cao & Angelo Bongiorno 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 MoS . Charged defects pose serious challenges to density functional theory (DFT) calculations relying on the use of 1–3 1, 2 supercells . The electrostatic energy of a periodic array of charged defects diverges , and to eliminate the singu- 1, 2 larity, the conventional approach consists in neutralizing the supercell by using a uniform charge background . Periodic DFT calculations rely on this solution to study charged systems and to calculate the formation energy of 3–7 a charged defect in a dielectric system . Unfortunately, the use of supercells and a uniform charge background introduce undesired interactions (between the charged defect and its replicas, and between the charged defect 1, 2 and the uniform background) that are fairly large and turn off slowly with the dimensions of the supercell . To account for these size effects, the formation energy, E , of a defect carrying a charge q is calculated by DFT as follows: DFTDFT EE =−− () qE (0), μμ −+ qE Δ fL L L (1) Ie DFT DFT where ΔE is the corrective energy term canceling the spurious interactions, Eq () and E (0) are the DFT L L L energies of the defected and reference systems, respectively, computed by using supercells of linear dimension L, and μ and μ are the chemical potentials of the ions and electrons forming the defect. In Eq. (1), ΔE is equal to I e L 1, 3, 6, 8, 9 E − E , that is the difference in electrostatic energy between the isolated defect in the infinite material ∞ L system, E , and the periodic array of defected supercells of linear dimension L compensated by a uniform charge background, E . Here, we present a new and original method to calculate ΔE = E − E . L L ∞ L e w Th ell-known formula introduced by Makov and Payne is valid and can be used to calculate Δ E only in the case of a point-charge defect in an isotropic and homogeneous bulk dielectric material . In recent years, various corrective schemes have been put forward to go beyond the Makov and Payne scheme and calculate ΔE in a 7, 8, 10–12 3 9 6 general situation , such as for a charged defect of finite size in an anisotropic , inhomogeneous , or ape- riodic system . All these previous corrective schemes rely on continuum electrostatics, structureless jellium-like 1 2 Department of Chemistry, College of Staten Island, Staten Island, NY, 10314, USA. Ph.D. Program in Chemistry, The Graduate Center of the City University of New York, New York, NY, 10016, USA. Ph.D. Program in Physics, The Graduate Center of the City University of New York, New York, NY, 10016, USA. Correspondence and requests for materials should be addressed to A.B. (email: angelo.bongiorno@csi.cuny.edu) Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 1 www.nature.com/scientificreports/ 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. 3 6 models for a defected dielectric system, and either analytical or numerical solutions of the Poisson equation. 1, 3, 9 Although all valid, these schemes are limited either in scope to a selected class of defects and systems 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 . 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 distribu- tion 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 . In mathematical terms, the potential energy function, U, of a periodic supercell encompassing N ions and n spherical shells takes the following form: UU =+ U , SpringsEwald (2) with n N 1 1 2 2 Uk =ΔsK +Δr , Springs ∑∑ ii jj 2 2 i== 1 j 1 (3) and Nn + Nn + qq    rr  −− R 1 ij ij n   U = erfc  Ewald ∑∑ ∑ 2  2 σ rr  −− R   n i ji ≠ ij n ew   σ G ew exp−  Nn + 2π   1 ++ Sq () G −  ∑∑ 2 i V G 2πσ G≠0 ew i   qQ Δs i i  +− erf , Δs  2 σ    i i ew (4) where, in the equations above, V is the volume of the supercell, n and R 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 elec- r tronic (q) charges, refers to the vector position of either an ion (r) or an electronic charge (s), σ is the parame- ew ter 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: Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 2 www.nature.com/scientificreports/ −1 Figure 2. (a) Energy vs. L of (black discs) a cubic array of unit point charges and (red discs) unit Gaussian −1 charges of width 4 a . (b) Energy vs. L of (black discs) a point charge with q = +2 and (blue discs) two unit Gaussian charges of width 1 a separated by a distance of 103 a . For clarity, schematics of the charged periodic systems, as well as the Makov-Payne formula and the Madelung constant obtained by fitting the energy values at large L are shown as insets. n N qQ  Sr −− R  j ij n i  U =− erfc  ∑∑∑  Gaussians  2 σ sr −− R   n i j i ij n   n n qq  ssR −−  ij ij n +− erfc , ∑∑∑  ssR −−   n i ji ≠ 2(σσ + ) ij n  ij    (5) where the σ ’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 to determine the ground state energy and dielectric dis- placements 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 = E − E to be used in Eq. (1), the polar- L ∞ L izable energy scheme in Eqs (2) and (5) is used to calculate E and, by extrapolation, E . Unless otherwise speci- L ∞ fied, LV = . 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 ape- riodic) systems, all cationic/anionic sites are equivalent, and thus our scheme requires to determine the val- ues 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, car - rying 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 inter- faces, 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 . 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 , and (iii) two Gaussian charge distributions of width 1 a placed on the diagonal of the supercell at a 0 0 distance of 103 a . In all cases, we use our scheme to calculate the energy of the array of charges interacting with a compensating uniform charge background. −1 e en Th ergy of a cubic array of point charges scales linearly with L , 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 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 3 www.nature.com/scientificreports/ −1 Figure 3. Energy vs. L of a cubic array of (black discs) unit point charges and (red discs) a unit Gaussian charge of width 2 Å in a model dielectric material with a permittivity of 4. For clarity, schematics of the defected supercells and energy scaling function are shown as insets. Material Cation, (Q, K); (q, k, σ) Anion, (Q, K); (q, k, σ) NaCl Na, (1.0, ∞); (0, −, −) Cl, (0.0, 2.95); (−1.0, 4.00, −) h-BN B, (1.0, 5.50); (0, −, −) N, (0.25, 6.0); (−1.25, 6.0, 1.0) 1L-MoS Mo, (0.0, ∞); (0, −, −) S, (2.0, ∞); (−2.0, 12.4, 0.5) Table 1. Anionic and cationic parameters used in our atomistic polarizable energy scheme to reproduce high- and low-frequency dielectric constants of NaCl, h-BN, and monolayer MoS . Harmonic spring constants, K and k, are expressed in eV/Å , the Gaussian width, σ, in Å, and atomic and electronic charges, Q and q, in electron charge unit. results . Our scheme yields also the correct behavior for the energy of an array of Gaussian charges (Fig. 2). In −1 −3 this case, in addition to terms scaling as L , the energy depends also on terms scaling as L , accounting for both 1, 10 the size and quadrupole moment of the charge distribution in a supercell (Fig. 2) . These results demonstrate that our method encompasses correction schemes relying on analytical formulas, applicable to charged defects of 1, 3 finite size in homogeneous bulk systems . 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 carry- ing 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 2 Å ), and then we perform a few trial-and-error calculations to determine the optimal value k = 11.0 eV/Å . We use this model dielectric material to calculate values of E − E with supercells hosting a unit charge defect in the ∞ L center of a cubic interstice of the lattice. −1 Figure 3 shows that, in the case of a point charge defect, the energy scales linearly with L , 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 calcu- lations 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 6, 14 defect in the bulk and at the (100) surface of NaCl , and a doubly negatively charged B antisite defect in the bulk of the hexagonal form of BN . Here we use our method to calculate the energy correction, E − E , and ∞ L 6, 9, 14 achieve comparison with recent computational studies of these defects . 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 ( =. 24) for NaCl, and of 6.6 and 3.5 for the components of the dielectric tensor perpendicular and parallel to the c axis of h-BN, respec- 9, 14 − tively ( =. 43 and  =. 26) . To describe the charged defects, we remove a Cl site from the lattice of bulk ∞ ∞ 2 − models, and at the terminal layer of surface models of NaCl. In the case of the B antisite defect in h-BN, we replace a N site with an atomic site carrying an ionic charge of 0.25e and an electronic charge of −3.25e, distrib- uted according to a Gaussian distribution of width 1.5 Å; both charges are fixed at the equilibrium position. Recently, Chen et al. used a periodic DFT scheme to calculate the formation energy of a Cl vacancy defect in rigid supercells of NaCl, i.e. consisting of immobile ions . Based entirely on DFT, this study obtained a value for E − E of about 0.7 eV, with L equal to twice the lattice constant, l , of NaCl . For the same value of L, our ∞ L 0 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 4 www.nature.com/scientificreports/ −1 − Figure 4. (a) Energy vs. L for a Cl vacancy defect in (red) rigid and (black) fully polarizable lattice models of (vivid colors) bulk NaCl and (pale colors) the NaCl(001) surface. Supercells of bulk NaCl are cubic, whereas surface models consist of tetragonal supercells with c = 2a and equal parts of material and vacuum along the c axis. Energy values are refereed to that of a defected supercell with a dimension along the a axis equal to 2l , −1 −2 where l is the experimental lattice constant. (b) Energy vs. L for a antisite defect in orthorhombic model structures of h-BN. Discs show energy values calculated by using our polarizable energy scheme, and dotted lines show linear interpolations. Insets show schematics of the defected materials systems. method yields an energy correction of about 0.8 eV (Fig. 4), in agreement with the result obtained by extrapolat- ing a few energy values calculated by DFT using supercells with L up to a value of 5l . 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 ). This result is in agreement with the energy correction of about 0.6 eV obtained by solving Ll =× 224 × the Poisson equation for continuum models of the defected surface . 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 −1 bulk and surface, at large L the energy scales linearly with L , 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, corresponding to the Madelung constant and permit- tivity of an array of tetragonal supercells, with c = 2a and encompassing equal parts of dielectric ( =. 59) and vacuum (Fig. 4). 2 − A good agreement with a recent DFT study is also obtained in the case of the B defect in h-BN. In fact, our method gives E − E = 0.85 eV to correct a formation energy obtained by DFT using a supercell of 8 × 8 × 3 unit ∞ L cells, in close agreement with the energy value of about 0.75 eV obtained by Kumagai and Oba . 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 . Overall, these results ∞ L 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 mon- 15, 16 olayer, bilayer, and tri-layer (and for completeness also in bulk) MoS . 17 18 Our periodic DFT calculations are carried out using normconserving pseudopotentials for both Mo and S, a generalized gradient approximation for the exchange and correlation energy , a plane-wave energy-cutoff of 80 Ry, and a semi-empirical corrective scheme to account for London dispersion interactions . To describe the MoS systems, we use orthorhombic supercells and atomic positions in accord to the 2H crystalline phase of MoS with a AB layer stacking pattern. Single layers include 30 MoS formula units, and 2D films are separated 2 2 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 are found ~3% larger than experimental values, and in the case of bilayer and trilayer MoS , the intra-layer spacing between S planes is found to be 3.17 Å, and the inter-layer separation between Mo Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 5 www.nature.com/scientificreports/ 0 −1 MoS ε Φ E (V ) E (V ) E -E 2 ⊥ f S f S ∞ L 1L 1.9 14.1 6.1 1.65 3.29 0.30 2L 2.6 14.8 5.5 1.87 2.67 0.17 3L 2.9 14.8 5.3 1.88 2.61 0.08 bulk 6.0 14.7 — 1.92 2.47 0.10 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 . 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. −1 Figure 5. (a) Energy vs. L for a singly negatively charged S vacancy in monolayer MoS calculated by using orthorhombic supercells with edges parallel to the monolayer equal to l and , and with a l = 3/2l yx perpendicular edge equal to (red) l = 2l , (black) l = l , or (blue) l = l + d (d = 6.1 Å is an approximate value z x z x z x for the monolayer thickness). Energy values are obtained with supercells of increasing l and plotted against L = lll . (b) Same as (a), with energy values plotted versus the rescaled linear dimension LL =  , with xy z SC , where is the component of the dielectric tensor parallel to the monolayer. Discs show  =− (1  )/ dl + 1 0 SC 0 z 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 = l . Energy values are referred to that one obtained with the z x same supercell used to calculate the defect formation energy by DFT. The inset shows a schematics of the defected 2D material. 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 systems . To calculate formation energies, we use half the energy of molecular S 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 , and values of 14.5 and 2.2 for bilayer MoS , 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 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 ; similar energy vs. 1/L curves are obtained also in the case of bilayer and trilayer MoS . e Th 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 −1 non-trivially with L . In the present case, we find that extrapolation of the energy values to large volumes is aided Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 6 www.nature.com/scientificreports/ by introducing a rescaled linear dimension L , related to L, thickness and in-plane dielectric constant of the 2D 8, 11, 12, 16 MoS systems as described in the caption of Fig. 5. In agreement with recent DFT studies , 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 = E − E , L ∞ L which in the case of a singly negatively charged S vacancy in monolayer MoS 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 rep- resentations, and a self-consistent treatment of the dielectric screening. e Th 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 calcula- tions . 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. References 1. Makov, G. & Payne, M. C. Periodic boundary conditions in ab initio calculations. Phys. Rev. B 51, 4014–4022, doi:10.1103/ PhysRevB.51.4014 (1995). 2. Leslie, M. & Gillan, N. J. The energy and elastic dipole tensor of defects in ionic crystals calculated by the supercell method. Journal of Physics C: Solid State Physics 18, 973–982, doi:10.1088/0022-3719/18/5/005 (1985). 3. Freysoldt, C., Neugebauer, J. & Van de Walle, C. G. Fully Ab Initio finite-size corrections for charged-defect supercell calculations. Phys. Rev. Lett. 102, 016402, doi:10.1103/PhysRevLett.102.016402 (2009). 4. Castleton, C. W. M., Höglund, A. & Mirbt, S. Managing the supercell approximation for charged defects in semiconductors: Finite- size scaling, charge correction factors, the band-gap problem, and the ab initio dielectric constant. Phys. Rev. B 73, 035215, doi:10.1103/PhysRevB.73.035215 (2006). 5. Lany, S. & Zunger, A. Assessment of correction methods for the band-gap problem and for finite-size effects in supercell defect calculations: Case studies for zno and gaas. Phys. Rev. B 78, 235104, doi:10.1103/PhysRevB.78.235104 (2008). 6. Komsa, H.-P. & Pasquarello, A. Finite-size supercell correction for charged defects at surfaces and interfaces. Phys. Rev. Lett. 110, 095505, doi:10.1103/PhysRevLett.110.095505 (2013). 7. Komsa, H.-P., Rantala, T. T. & Pasquarello, A. Finite-size supercell correction schemes for charged defect calculations. Phys. Rev. B 86, 045112, doi:10.1103/PhysRevB.86.045112 (2012). 8. Wang, D. et al. Determination of formation and ionization energies of charged defects in two-dimensional materials. Phys. Rev. Lett. 114, 196801, doi:10.1103/PhysRevLett.114.196801 (2015). 9. Kumagai, Y. & Oba, F. Electrostatics-based finite-size corrections for first-principles point defect calculations. Phys. Rev. B 89, 195205, doi:10.1103/PhysRevB.89.195205 (2014). 10. Dabo, I., Kozinsky, B., Singh-Miller, N. E. & Marzari, N. Electrostatics in periodic boundary conditions and real-space corrections. Phys. Rev. B 77, 115139, doi:10.1103/PhysRevB.77.115139 (2008). 11. Chan, T.-L., Zhang, S. B. & Chelikowsky, J. R. Charged dopants in semiconductor nanowires under partially periodic boundary conditions. Phys. Rev. B 83, 245440, doi:10.1103/PhysRevB.83.245440 (2011). 12. Komsa, H.-P., Berseneva, N., Krasheninnikov, A. V. & Nieminen, R. M. Charged point defects in the flatland: Accurate formation energy calculations in two-dimensional materials. Phys. Rev. X 4, 031044 (2014). 13. Allen, M. & Tildesley, D. Computer Simulation of Liquids (Oxford: Clarendon Pr, 1987). 14. Chen, W., Tegenkamp, C., Pfnür, H. & Bredow, T. Color centers in nacl by hybrid functionals. Phys. Rev. B 82, 104106, doi:10.1103/ PhysRevB.82.104106 (2010). 15. Komsa, H.-P. & Krasheninnikov, A. V. Native defects in bulk and monolayer mos from first principles. Phys. Rev. B 91, 125304, doi:10.1103/PhysRevB.91.125304 (2015). 16. Noh, J.-Y., Kim, H. & Kim, Y.-S. Stability and electronic structures of native defects in single-layer mos . Phys. Rev. B 89, 205417, doi:10.1103/PhysRevB.89.205417 (2014). 17. Giannozzi, P. et al. Quantum espresso: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009). 18. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 43, 1993–2006, doi:10.1103/ PhysRevB.43.1993 (1991). 19. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868, doi:10.1103/PhysRevLett.77.3865 (1996). 20. Grimme, S. Semiempirical gga-type density functional constructed with a long-range dispersion correction. Journal of Computational Chemistry 27, 1787–1799, doi:10.1002/jcc.20495 (2006). Acknowledgements We acknowledge the support of the NSF grant CMMI-1436375 and the computational resources provided by the CUNY High Performance Computing Center (NSF Grants: CNS-0855217, CNS-0958379 and ACI-1126113). Author Contributions T.C. implemented the corrective scheme and carried out the calculations. A.B. devised the conceptual basis of the new corrective method and contributed to its technical implementation. Additional Information Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 7 www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre- ative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not per- mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 8 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Scientific Reports Springer Journals

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

Free
8 pages
Loading next page...
 
/lp/springer_journal/atomistic-corrective-scheme-for-supercell-density-functional-theory-6ENTP0Yrpf
Publisher
Nature Publishing Group UK
Copyright
Copyright © 2017 by The Author(s)
Subject
Science, Humanities and Social Sciences, multidisciplinary; Science, Humanities and Social Sciences, multidisciplinary; Science, multidisciplinary
eISSN
2045-2322
D.O.I.
10.1038/s41598-017-02986-5
Publisher site
See Article on Publisher Site

Abstract

www.nature.com/scientificreports OPEN Atomistic Corrective Scheme for Supercell Density Functional Theory Calculations of Charged Received: 20 March 2017 Defects Accepted: 20 April 2017 Published: xx xx xxxx 1 1,2,3 Tengfei Cao & Angelo Bongiorno 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 MoS . Charged defects pose serious challenges to density functional theory (DFT) calculations relying on the use of 1–3 1, 2 supercells . The electrostatic energy of a periodic array of charged defects diverges , and to eliminate the singu- 1, 2 larity, the conventional approach consists in neutralizing the supercell by using a uniform charge background . Periodic DFT calculations rely on this solution to study charged systems and to calculate the formation energy of 3–7 a charged defect in a dielectric system . Unfortunately, the use of supercells and a uniform charge background introduce undesired interactions (between the charged defect and its replicas, and between the charged defect 1, 2 and the uniform background) that are fairly large and turn off slowly with the dimensions of the supercell . To account for these size effects, the formation energy, E , of a defect carrying a charge q is calculated by DFT as follows: DFTDFT EE =−− () qE (0), μμ −+ qE Δ fL L L (1) Ie DFT DFT where ΔE is the corrective energy term canceling the spurious interactions, Eq () and E (0) are the DFT L L L energies of the defected and reference systems, respectively, computed by using supercells of linear dimension L, and μ and μ are the chemical potentials of the ions and electrons forming the defect. In Eq. (1), ΔE is equal to I e L 1, 3, 6, 8, 9 E − E , that is the difference in electrostatic energy between the isolated defect in the infinite material ∞ L system, E , and the periodic array of defected supercells of linear dimension L compensated by a uniform charge background, E . Here, we present a new and original method to calculate ΔE = E − E . L L ∞ L e w Th ell-known formula introduced by Makov and Payne is valid and can be used to calculate Δ E only in the case of a point-charge defect in an isotropic and homogeneous bulk dielectric material . In recent years, various corrective schemes have been put forward to go beyond the Makov and Payne scheme and calculate ΔE in a 7, 8, 10–12 3 9 6 general situation , such as for a charged defect of finite size in an anisotropic , inhomogeneous , or ape- riodic system . All these previous corrective schemes rely on continuum electrostatics, structureless jellium-like 1 2 Department of Chemistry, College of Staten Island, Staten Island, NY, 10314, USA. Ph.D. Program in Chemistry, The Graduate Center of the City University of New York, New York, NY, 10016, USA. Ph.D. Program in Physics, The Graduate Center of the City University of New York, New York, NY, 10016, USA. Correspondence and requests for materials should be addressed to A.B. (email: angelo.bongiorno@csi.cuny.edu) Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 1 www.nature.com/scientificreports/ 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. 3 6 models for a defected dielectric system, and either analytical or numerical solutions of the Poisson equation. 1, 3, 9 Although all valid, these schemes are limited either in scope to a selected class of defects and systems 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 . 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 distribu- tion 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 . In mathematical terms, the potential energy function, U, of a periodic supercell encompassing N ions and n spherical shells takes the following form: UU =+ U , SpringsEwald (2) with n N 1 1 2 2 Uk =ΔsK +Δr , Springs ∑∑ ii jj 2 2 i== 1 j 1 (3) and Nn + Nn + qq    rr  −− R 1 ij ij n   U = erfc  Ewald ∑∑ ∑ 2  2 σ rr  −− R   n i ji ≠ ij n ew   σ G ew exp−  Nn + 2π   1 ++ Sq () G −  ∑∑ 2 i V G 2πσ G≠0 ew i   qQ Δs i i  +− erf , Δs  2 σ    i i ew (4) where, in the equations above, V is the volume of the supercell, n and R 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 elec- r tronic (q) charges, refers to the vector position of either an ion (r) or an electronic charge (s), σ is the parame- ew ter 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: Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 2 www.nature.com/scientificreports/ −1 Figure 2. (a) Energy vs. L of (black discs) a cubic array of unit point charges and (red discs) unit Gaussian −1 charges of width 4 a . (b) Energy vs. L of (black discs) a point charge with q = +2 and (blue discs) two unit Gaussian charges of width 1 a separated by a distance of 103 a . For clarity, schematics of the charged periodic systems, as well as the Makov-Payne formula and the Madelung constant obtained by fitting the energy values at large L are shown as insets. n N qQ  Sr −− R  j ij n i  U =− erfc  ∑∑∑  Gaussians  2 σ sr −− R   n i j i ij n   n n qq  ssR −−  ij ij n +− erfc , ∑∑∑  ssR −−   n i ji ≠ 2(σσ + ) ij n  ij    (5) where the σ ’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 to determine the ground state energy and dielectric dis- placements 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 = E − E to be used in Eq. (1), the polar- L ∞ L izable energy scheme in Eqs (2) and (5) is used to calculate E and, by extrapolation, E . Unless otherwise speci- L ∞ fied, LV = . 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 ape- riodic) systems, all cationic/anionic sites are equivalent, and thus our scheme requires to determine the val- ues 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, car - rying 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 inter- faces, 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 . 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 , and (iii) two Gaussian charge distributions of width 1 a placed on the diagonal of the supercell at a 0 0 distance of 103 a . In all cases, we use our scheme to calculate the energy of the array of charges interacting with a compensating uniform charge background. −1 e en Th ergy of a cubic array of point charges scales linearly with L , 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 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 3 www.nature.com/scientificreports/ −1 Figure 3. Energy vs. L of a cubic array of (black discs) unit point charges and (red discs) a unit Gaussian charge of width 2 Å in a model dielectric material with a permittivity of 4. For clarity, schematics of the defected supercells and energy scaling function are shown as insets. Material Cation, (Q, K); (q, k, σ) Anion, (Q, K); (q, k, σ) NaCl Na, (1.0, ∞); (0, −, −) Cl, (0.0, 2.95); (−1.0, 4.00, −) h-BN B, (1.0, 5.50); (0, −, −) N, (0.25, 6.0); (−1.25, 6.0, 1.0) 1L-MoS Mo, (0.0, ∞); (0, −, −) S, (2.0, ∞); (−2.0, 12.4, 0.5) Table 1. Anionic and cationic parameters used in our atomistic polarizable energy scheme to reproduce high- and low-frequency dielectric constants of NaCl, h-BN, and monolayer MoS . Harmonic spring constants, K and k, are expressed in eV/Å , the Gaussian width, σ, in Å, and atomic and electronic charges, Q and q, in electron charge unit. results . Our scheme yields also the correct behavior for the energy of an array of Gaussian charges (Fig. 2). In −1 −3 this case, in addition to terms scaling as L , the energy depends also on terms scaling as L , accounting for both 1, 10 the size and quadrupole moment of the charge distribution in a supercell (Fig. 2) . These results demonstrate that our method encompasses correction schemes relying on analytical formulas, applicable to charged defects of 1, 3 finite size in homogeneous bulk systems . 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 carry- ing 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 2 Å ), and then we perform a few trial-and-error calculations to determine the optimal value k = 11.0 eV/Å . We use this model dielectric material to calculate values of E − E with supercells hosting a unit charge defect in the ∞ L center of a cubic interstice of the lattice. −1 Figure 3 shows that, in the case of a point charge defect, the energy scales linearly with L , 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 calcu- lations 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 6, 14 defect in the bulk and at the (100) surface of NaCl , and a doubly negatively charged B antisite defect in the bulk of the hexagonal form of BN . Here we use our method to calculate the energy correction, E − E , and ∞ L 6, 9, 14 achieve comparison with recent computational studies of these defects . 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 ( =. 24) for NaCl, and of 6.6 and 3.5 for the components of the dielectric tensor perpendicular and parallel to the c axis of h-BN, respec- 9, 14 − tively ( =. 43 and  =. 26) . To describe the charged defects, we remove a Cl site from the lattice of bulk ∞ ∞ 2 − models, and at the terminal layer of surface models of NaCl. In the case of the B antisite defect in h-BN, we replace a N site with an atomic site carrying an ionic charge of 0.25e and an electronic charge of −3.25e, distrib- uted according to a Gaussian distribution of width 1.5 Å; both charges are fixed at the equilibrium position. Recently, Chen et al. used a periodic DFT scheme to calculate the formation energy of a Cl vacancy defect in rigid supercells of NaCl, i.e. consisting of immobile ions . Based entirely on DFT, this study obtained a value for E − E of about 0.7 eV, with L equal to twice the lattice constant, l , of NaCl . For the same value of L, our ∞ L 0 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 4 www.nature.com/scientificreports/ −1 − Figure 4. (a) Energy vs. L for a Cl vacancy defect in (red) rigid and (black) fully polarizable lattice models of (vivid colors) bulk NaCl and (pale colors) the NaCl(001) surface. Supercells of bulk NaCl are cubic, whereas surface models consist of tetragonal supercells with c = 2a and equal parts of material and vacuum along the c axis. Energy values are refereed to that of a defected supercell with a dimension along the a axis equal to 2l , −1 −2 where l is the experimental lattice constant. (b) Energy vs. L for a antisite defect in orthorhombic model structures of h-BN. Discs show energy values calculated by using our polarizable energy scheme, and dotted lines show linear interpolations. Insets show schematics of the defected materials systems. method yields an energy correction of about 0.8 eV (Fig. 4), in agreement with the result obtained by extrapolat- ing a few energy values calculated by DFT using supercells with L up to a value of 5l . 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 ). This result is in agreement with the energy correction of about 0.6 eV obtained by solving Ll =× 224 × the Poisson equation for continuum models of the defected surface . 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 −1 bulk and surface, at large L the energy scales linearly with L , 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, corresponding to the Madelung constant and permit- tivity of an array of tetragonal supercells, with c = 2a and encompassing equal parts of dielectric ( =. 59) and vacuum (Fig. 4). 2 − A good agreement with a recent DFT study is also obtained in the case of the B defect in h-BN. In fact, our method gives E − E = 0.85 eV to correct a formation energy obtained by DFT using a supercell of 8 × 8 × 3 unit ∞ L cells, in close agreement with the energy value of about 0.75 eV obtained by Kumagai and Oba . 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 . Overall, these results ∞ L 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 mon- 15, 16 olayer, bilayer, and tri-layer (and for completeness also in bulk) MoS . 17 18 Our periodic DFT calculations are carried out using normconserving pseudopotentials for both Mo and S, a generalized gradient approximation for the exchange and correlation energy , a plane-wave energy-cutoff of 80 Ry, and a semi-empirical corrective scheme to account for London dispersion interactions . To describe the MoS systems, we use orthorhombic supercells and atomic positions in accord to the 2H crystalline phase of MoS with a AB layer stacking pattern. Single layers include 30 MoS formula units, and 2D films are separated 2 2 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 are found ~3% larger than experimental values, and in the case of bilayer and trilayer MoS , the intra-layer spacing between S planes is found to be 3.17 Å, and the inter-layer separation between Mo Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 5 www.nature.com/scientificreports/ 0 −1 MoS ε Φ E (V ) E (V ) E -E 2 ⊥ f S f S ∞ L 1L 1.9 14.1 6.1 1.65 3.29 0.30 2L 2.6 14.8 5.5 1.87 2.67 0.17 3L 2.9 14.8 5.3 1.88 2.61 0.08 bulk 6.0 14.7 — 1.92 2.47 0.10 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 . 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. −1 Figure 5. (a) Energy vs. L for a singly negatively charged S vacancy in monolayer MoS calculated by using orthorhombic supercells with edges parallel to the monolayer equal to l and , and with a l = 3/2l yx perpendicular edge equal to (red) l = 2l , (black) l = l , or (blue) l = l + d (d = 6.1 Å is an approximate value z x z x z x for the monolayer thickness). Energy values are obtained with supercells of increasing l and plotted against L = lll . (b) Same as (a), with energy values plotted versus the rescaled linear dimension LL =  , with xy z SC , where is the component of the dielectric tensor parallel to the monolayer. Discs show  =− (1  )/ dl + 1 0 SC 0 z 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 = l . Energy values are referred to that one obtained with the z x same supercell used to calculate the defect formation energy by DFT. The inset shows a schematics of the defected 2D material. 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 systems . To calculate formation energies, we use half the energy of molecular S 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 , and values of 14.5 and 2.2 for bilayer MoS , 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 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 ; similar energy vs. 1/L curves are obtained also in the case of bilayer and trilayer MoS . e Th 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 −1 non-trivially with L . In the present case, we find that extrapolation of the energy values to large volumes is aided Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 6 www.nature.com/scientificreports/ by introducing a rescaled linear dimension L , related to L, thickness and in-plane dielectric constant of the 2D 8, 11, 12, 16 MoS systems as described in the caption of Fig. 5. In agreement with recent DFT studies , 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 = E − E , L ∞ L which in the case of a singly negatively charged S vacancy in monolayer MoS 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 rep- resentations, and a self-consistent treatment of the dielectric screening. e Th 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 calcula- tions . 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. References 1. Makov, G. & Payne, M. C. Periodic boundary conditions in ab initio calculations. Phys. Rev. B 51, 4014–4022, doi:10.1103/ PhysRevB.51.4014 (1995). 2. Leslie, M. & Gillan, N. J. The energy and elastic dipole tensor of defects in ionic crystals calculated by the supercell method. Journal of Physics C: Solid State Physics 18, 973–982, doi:10.1088/0022-3719/18/5/005 (1985). 3. Freysoldt, C., Neugebauer, J. & Van de Walle, C. G. Fully Ab Initio finite-size corrections for charged-defect supercell calculations. Phys. Rev. Lett. 102, 016402, doi:10.1103/PhysRevLett.102.016402 (2009). 4. Castleton, C. W. M., Höglund, A. & Mirbt, S. Managing the supercell approximation for charged defects in semiconductors: Finite- size scaling, charge correction factors, the band-gap problem, and the ab initio dielectric constant. Phys. Rev. B 73, 035215, doi:10.1103/PhysRevB.73.035215 (2006). 5. Lany, S. & Zunger, A. Assessment of correction methods for the band-gap problem and for finite-size effects in supercell defect calculations: Case studies for zno and gaas. Phys. Rev. B 78, 235104, doi:10.1103/PhysRevB.78.235104 (2008). 6. Komsa, H.-P. & Pasquarello, A. Finite-size supercell correction for charged defects at surfaces and interfaces. Phys. Rev. Lett. 110, 095505, doi:10.1103/PhysRevLett.110.095505 (2013). 7. Komsa, H.-P., Rantala, T. T. & Pasquarello, A. Finite-size supercell correction schemes for charged defect calculations. Phys. Rev. B 86, 045112, doi:10.1103/PhysRevB.86.045112 (2012). 8. Wang, D. et al. Determination of formation and ionization energies of charged defects in two-dimensional materials. Phys. Rev. Lett. 114, 196801, doi:10.1103/PhysRevLett.114.196801 (2015). 9. Kumagai, Y. & Oba, F. Electrostatics-based finite-size corrections for first-principles point defect calculations. Phys. Rev. B 89, 195205, doi:10.1103/PhysRevB.89.195205 (2014). 10. Dabo, I., Kozinsky, B., Singh-Miller, N. E. & Marzari, N. Electrostatics in periodic boundary conditions and real-space corrections. Phys. Rev. B 77, 115139, doi:10.1103/PhysRevB.77.115139 (2008). 11. Chan, T.-L., Zhang, S. B. & Chelikowsky, J. R. Charged dopants in semiconductor nanowires under partially periodic boundary conditions. Phys. Rev. B 83, 245440, doi:10.1103/PhysRevB.83.245440 (2011). 12. Komsa, H.-P., Berseneva, N., Krasheninnikov, A. V. & Nieminen, R. M. Charged point defects in the flatland: Accurate formation energy calculations in two-dimensional materials. Phys. Rev. X 4, 031044 (2014). 13. Allen, M. & Tildesley, D. Computer Simulation of Liquids (Oxford: Clarendon Pr, 1987). 14. Chen, W., Tegenkamp, C., Pfnür, H. & Bredow, T. Color centers in nacl by hybrid functionals. Phys. Rev. B 82, 104106, doi:10.1103/ PhysRevB.82.104106 (2010). 15. Komsa, H.-P. & Krasheninnikov, A. V. Native defects in bulk and monolayer mos from first principles. Phys. Rev. B 91, 125304, doi:10.1103/PhysRevB.91.125304 (2015). 16. Noh, J.-Y., Kim, H. & Kim, Y.-S. Stability and electronic structures of native defects in single-layer mos . Phys. Rev. B 89, 205417, doi:10.1103/PhysRevB.89.205417 (2014). 17. Giannozzi, P. et al. Quantum espresso: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009). 18. Troullier, N. & Martins, J. L. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B 43, 1993–2006, doi:10.1103/ PhysRevB.43.1993 (1991). 19. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868, doi:10.1103/PhysRevLett.77.3865 (1996). 20. Grimme, S. Semiempirical gga-type density functional constructed with a long-range dispersion correction. Journal of Computational Chemistry 27, 1787–1799, doi:10.1002/jcc.20495 (2006). Acknowledgements We acknowledge the support of the NSF grant CMMI-1436375 and the computational resources provided by the CUNY High Performance Computing Center (NSF Grants: CNS-0855217, CNS-0958379 and ACI-1126113). Author Contributions T.C. implemented the corrective scheme and carried out the calculations. A.B. devised the conceptual basis of the new corrective method and contributed to its technical implementation. Additional Information Competing Interests: The authors declare that they have no competing interests. Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 7 www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Cre- ative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not per- mitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. © The Author(s) 2017 Scientific Repo R ts | 7: 2834 | DOI:10.1038/s41598-017-02986-5 8

Journal

Scientific ReportsSpringer Journals

Published: Jun 6, 2017

References

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off