Erratum : Charged Point Defects in the Flatland : Accurate Formation Energy Calculations in Two-Dimensional Materials

DOI: 10.1103/PhysRevX.8.039902 The correction method presented in the original paper involves calculating the electrostatic energy difference between model systems where charge is placed in (i) a small supercell, of the same size as used in the density functional theory (DFT) calculations, and subject to periodic boundary conditions Eperiodic, and (ii) isolated charge Eisolated. We proposed that the latter can be efficiently obtained by extrapolation from the former. While conceptually true, the extrapolations presented in the original paper were not evaluated from sufficiently large supercell sizes, which led to underestimation of the correction. In this Erratum, we briefly analyze the origin of the error, show how to circumvent it, and provide corrected energy values. The reason for this oversight was that, unlike in bulk systems with a homogeneous dielectric constant, the electrostatic energy for 2D systems shows a sudden change in the asymptotic behavior at large supercell sizes, evidencing the existence of two characteristic length scales. This is illustrated in Fig. 1, which shows the extrapolation of Eperiodic as a function of the inverse of the supercell size scaling parameter α in the case of the MnMoðþ1Þ defect in MoS2. Here, α 1⁄4 1 corresponds to a 4 × 4 supercell and c 1⁄4 12 Å, and all cell dimensions increase uniformly upon increasing α (Lx;y;z → αLx;y;z). In the original paper, the extrapolation was performed using a second-order polynomial with respect to 1=α, fitted to energies from the five smallest supercell sizes, and it was found to result in a good fit (cf. the dashed-dotted red line and solid circles in Fig. 1). For 1=α ≲ 0.25, corresponding to supercells larger than 16 × 16, there is a sudden upturn. This issue was first pointed out by Noh et al. [1] (and later on reported at least in Refs. [2,3]), who suggested using a fifth-order polynomial for the fitting. From the data shown in Fig. 1 and the higher-order polynomial fits, it is still not clear if they reproduce the correct extrapolated value. On the other hand, an efficient method for calculating Eisolated was recently presented in Ref. [4] that

The correction method presented in the original paper involves calculating the electrostatic energy difference between model systems where charge is placed in (i) a small supercell, of the same size as used in the density functional theory (DFT) calculations, and subject to periodic boundary conditions E periodic , and (ii) isolated charge E isolated .We proposed that the latter can be efficiently obtained by extrapolation from the former.While conceptually true, the extrapolations presented in the original paper were not evaluated from sufficiently large supercell sizes, which led to underestimation of the correction.In this Erratum, we briefly analyze the origin of the error, show how to circumvent it, and provide corrected energy values.
The reason for this oversight was that, unlike in bulk systems with a homogeneous dielectric constant, the electrostatic energy for 2D systems shows a sudden change in the asymptotic behavior at large supercell sizes, evidencing the existence of two characteristic length scales.This is illustrated in Fig. 1, which shows the extrapolation of E periodic as a function of the inverse of the supercell size scaling parameter α in the case of the Mn Mo ðþ1Þ defect in MoS 2 .Here, α ¼ 1 corresponds to a 4 × 4 supercell and c ¼ 12 Å, and all cell dimensions increase uniformly upon increasing α (L x;y;z → αL x;y;z ).In the original paper, the extrapolation was performed using a second-order polynomial with respect to 1=α, fitted to energies from the five smallest supercell sizes, and it was found to result in a good fit (cf. the dashed-dotted red line and solid circles in Fig. 1).For 1=α ≲ 0.25, corresponding to supercells larger than 16 × 16, there is a sudden upturn.This issue was first pointed out by Noh et al. [1] (and later on reported at least in Refs.[2,3]), who suggested using a fifth-order polynomial for the fitting.
From the data shown in Fig. 1 and the higher-order polynomial fits, it is still not clear if they reproduce the correct extrapolated value.On the other hand, an efficient method for calculating E isolated was recently presented in Ref. [4]  eliminates the need for extrapolation.It relies on Bessel expansion in lateral directions and a truncated Green's function approach in the perpendicular direction.It can be readily generalized to an anisotropic dielectric tensor profile and to a charge density comprised of several Gaussians with different widths σ in the lateral and perpendicular directions.The value obtained with this approach is found to reside 30 meV higher than the fifth-order polynomial fit (black square in Fig. 1).While direct evaluation of E isolated should usually be preferred, if an extrapolation procedure is nevertheless needed, one can make use of the fact that the energies at x ¼ 1=α → 0 simply correspond to the Madelung energy E M of a point charge in a homogeneous background charge and no screening (ε ¼ 1).The lack of screening is clear for interactions between charges in different sheets, going largely through vacuum.While the interactions between charges in the same sheet are partially screened, at increasing distance, the screening becomes weaker as a larger fraction of the field lines pass through the vacuum, and the effective interaction approaches that of an unscreened one [5].This change in the screening is the cause for the two length scales.One could then fit only the energies in the weak screening regime with a low-order polynomial, if such large supercells are computationally tractable [3].Alternatively, we found that the energies can be well fitted in both regimes when the second-order polynomial is appended with an exponential term: where c i are the fitting parameters, and d ¼ ðc 1 − ∂E M =∂xÞ=c 3 guarantees that the energy gradient is similar to the Madelung energy at x ¼ 1=α → 0, i.e., ∂E=∂x ¼ ∂E M =∂x.As seen in Fig. 1, this works better than the fifth-order polynomial with the given number of data points reaching supercell sizes of 40 × 40.Naturally, as more data points are calculated, the extrapolations from different functions should converge.Similarly, good agreement with the E isolated and the extrapolations using Eq.(1) were obtained in all cases, giving us confidence that we indeed have reached the correct value for E isolated .Finally, we note that the Madelung energy scaled with the average dielectric constant yields fairly good corrected values.The magnitude of errors induced by this are (i) 0.185 eV in the case of C B (and adatoms C ad , B ad , N ad at þ1 and −1 charge states via the special vacuum construction) in h-BN; (ii) 0.177 eV in the case of C N in h-BN; (iii) 0.101 eV in the case of Mn Mo in MoS 2 ; and (iv) 0.099 eV in the case of F S in MoS 2 .Correspondingly, a rigid shift should be applied to all corrected values in Figs.2-4, both formation energies and charge transition levels.Here, it is sufficient to present revised versions of Tables I and II.While the errors are non-negligible, they are relatively small when compared to the calculated band gaps of these materials: 4.67 eV for h-BN and 1.69 eV for MoS 2 .As a result of this shift, extrapolations from DFT values using a second-order polynomial will not match the corrected values, due to the same problems as discussed above.
FIG.1.Energy of the periodic system E periodic , as a function of the inverse supercell scaling parameter α.Only the first five calculated data points (marked by filled circles) were used for fitting of the second-order polynomial.
Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License.Further distribution of this work must maintain attribution to the author(s) and the published articles title, journal citation, and DOI.

TABLE I .
Special vacuums L s z (in Å) for different lateral supercell sizes (in multiples of the primitive cell) of monolayer h-BN and MoS 2 .

TABLE II .
Formation energies and charge transition levels for all considered defects in 2D h-BN and MoS 2 (in eV).For formation energies of charged defects, the Fermi level is set to the valence band maximum.Charge transition levels are referenced to the valence band maximum.Substitutional defect formation energies are taken from the extrapolation, whereas adatom energies are taken from the special vacuum calculation with a 6 × 6 supercell.Three dots denote cases where the relevant charge state is unstable.