Intrinsic point defects are indispensable in crystals at a finite temperature and can have important influences on the fundamental properties of crystalline materials^{[1, 2]}. Extrinsic elements are usually introduced into the crystalline lattice as dopants intentionally or as impurities unintentionally, which provides an extrinsic freedom for manipulating the properties of crystalline materials. The influences of defects and dopants (impurities) in semiconductors are more significant than those in metals or structural materials, because the applications of semiconductors in functional electronic and optoelectronic devices care more about the electrical and optical properties, which can be significantly changed even when a small number of defects and dopants exist in the lattice, e.g., defects or dopants with a density of 10¹⁵–10¹⁸ cm^–3 (1 defect or dopant among 10⁴–10⁷ atoms) can change the electrical and optical properties of semiconductors dramatically. For the simple elemental semiconductors such as Si and Ge, the number of intrinsic point defects is usually small^[3], so extrinsic doping is required for producing electron or hole carriers and thus achieving n-type or p-type electrical conductivity. For binary, ternary, quaternary or even multinary compound semiconductors, the number of possible point defects can be large, so their defect/dopant physics can be more complicated and thus more flexible, e.g., their electrical conductivity may be tuned through controlling the formation of different point defects^[4]. Considering the important influences of defects (dopants) on semiconductor properties and the complexity of the defect/dopant physics in new, multinary or low-symmetry semiconductors, the studies on defects and dopants are fundamental to the development of semiconductor physics and applications.

Many experimental techniques have been developed for the characterization of defects and dopants in semiconductors, such as the photoluminescence (PL) spectrum, positron annihilation spectroscopy (PAS), electron paramagnetic resonance (EPR) and deep-level transient spectroscopy (DLTS)^[5]. These techniques usually just give the measured spectra, from which the defect and dopant properties, such as the density, energy levels and carrier capture cross sections can be obtained through fitting the spectra according to a certain physical model. However, what kind of defects and dopants determine the experimental spectra (including the fitted properties of defects and dopants) is usually speculated according to the analysis on the properties of the possible defects and dopants, which is semi-empirical and imposes a limit to the accurate and quantitative defect and dopant engineering. For instance, if the origin defects of the deep-level recombination centers in optoelectronic semiconductors are not clear, it is impossible to suppress their formation and increase the minority carrier lifetime by accurately controlling the growth conditions.

Since the 1990s, significant progresses have been made in the ab-initio (first-principles) calculation of defect and dopant properties based on the density functional theory (DFT) and the supercell model^[6-8]. The calculated formation energies, charge-state transition levels and carrier capture cross sections of defects and dopants have been widely used for revealing the microscopic origin of the experimental PL, PAS, EPR and DLTS spectra, and guiding the fabrication of high-quality materials and high-performance devices^[9-11]. However, there are still unsolved issues that may cause considerable errors in the calculated defect/dopant properties or even misunderstanding of the defect/dopant physics. (i) The approximations to exchange-correlation functionals cause the underestimation of the band gap and the incorrect location and occupation of the defect/dopant levels in the band gap^[12], which can result in errors in the calculated formation energy and transition level. Such errors are serious for the local density approximation (LDA)^[13] and generalized gradient approximation (GGA)^[14], so several correction schemes were proposed^{[7, 15-18]}. Recently, hybrid functional calculations were shown to predict more reasonable band gaps for many semiconductors, so the errors in the calculated defect properties are also corrected largely (not perfectly, because the hybrid functional may give incorrect localization of defect states and also depends on the empirical exact exchange ratio parameter)^[19-21]. (ii) The calculation of defect/dopant formation energies and densities requires the determination of the allowed chemical potential ranges for all the component and dopant elements. For ternary, quaternary and multinary compound semiconductors, a large number of secondary phases competing against the host material probably exist, which may all limit the stable range of elemental chemical potentials. Therefore, if some important secondary phases are not considered in the determination of elemental chemical potentials, the calculations of defect/dopant formation energies and densities and even the stability of the compound can be wrong. (iii) For low-symmetry and multinary compound semiconductors, the possible defect/dopant types can be many and various. For each type of defect/dopant, there may be multiple non-equivalent sites and structural configurations, which can change for different charge states. If some important defect-types/sites/configurations/charge-states are not considered, the calculated results are incomplete and inaccurate. (iv) The finite supercell size corrections, such as the corrections for the electrostatic potential alignment and image charge interaction^{[12, 22, 23]}, are still controversial. Especially, different correction schemes can give rise to different results. (v) In most of the studies, only defect/dopant formation energies and transition levels were calculated, but the exact values of equilibrium defect densities, carrier densities and Fermi level are not calculated. Without the direct calculation of the defect densities, one may have misunderstandings on the influences of a certain defect, because all the defects (dopants) are correlated by the carrier densities and Fermi level, and the defect correlation may give rise to unexpected and anti-chemical-intuition changes in the defect densities^[24].

To solve those problems in a standard and comprehensive way, we developed a software package, the Defect and Dopant ab-initio Simulation Package (DASP), which can perform automated calculations of defect and dopant properties based on the supercell model and ab-initio DFT calculations. Using the software, the defect formation energies, charge-state transition levels, carrier capture cross sections and even the PL spectrum can be calculated based on the atomic structure, total energy, electronic structure, phonon spectrum and electron–phonon coupling matrix calculated by the ab-initio DFT softwares using different approximations to exchange-correlation functionals. All the possible competing secondary compounds in the materials genome database are considered for the accurate calculation of the thermodynamic stability and elemental chemical potential ranges of compound semiconductors. Various defect types, atomic sites, structure configurations and charge states can be considered, and the corrections for the electrostatic potential alignment and image charge interaction can be included automatically. With all the defects and dopants considered, the equilibrium defect densities, carrier densities and Fermi level can be calculated for the samples grown under different chemical potential conditions and different temperature. For high-density defects and dopants, their carrier dynamics properties can also be calculated, despite the heavier computational cost. Therefore, DASP is developed to be a toolbox for the automated theoretical prediction of defect and dopant properties in semiconductors.

DASP is composed of four modules: Thermodynamic Stability Calculation (TSC), Defect Energy Calculation (DEC), Defect Density Calculation (DDC), and Carrier Dynamics Calculation (CDC), as plotted in Fig. 1.

The detailed workflow of DASP is plotted schematically in Fig. 2, in which the green part shows the DEC module, the blue part shows the TSC module, the purple part shows the DDC module and the red part shows the CDC module. The only necessary input is the crystal structure file of the semiconductor. The intermediate and final output files include:

(i) TSC: the chemical potential range of component elements that stabilizes the compound semiconductor (which is also a descriptor of the thermodynamic stability of the compound) and the allowed the highest chemical potential of the dopant elements;

(ii) DEC: the formation energies of defects and dopants in different charge states, as functions of the elemental chemical potentials and Fermi level (electronic chemical potentials), from which the charge-state transition levels can be derived;

(iii) the equilibrium-state Fermi level, densities of electron and hole carriers, and densities of the defects and dopants in different charge states, as functions of the elemental chemical potentials at growth temperature and working (measuring) temperature;

(iv) the carrier capture cross sections, the radiative and non-radiative carrier recombination rates and lifetime, and the PL spectra.

Among the four modules, the most fundamental one is DEC, whose core function is for calculating the defect/dopant formation energy using the supercell model and ab-initio DFT calculations. In the supercell model, the formation energy of a point defect α in the charge state q can be calculated as^{[25, 26]},

$ \Delta {E}_{\mathrm{f}}\left(\alpha ,q\right) = E\left(\alpha ,q\right) - E\left(\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}\right) - {\sum }_{i}{n}_{i}\left({\mu }_{i} + {E}_{i}\right) + q\left({E}_{\mathrm{F}} + {\epsilon }_{\mathrm{V}\mathrm{B}\mathrm{M}}\right) + {E}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}} , $  (1)

where and are the total energies of the supercells with and without a defect (dopant), is the elemental chemical potential referenced to the total energy of the pure solid/gas elementary phase, is the Fermi level referenced to the eigenvalue of the valence band maximum (VBM) level of the bulk supercell. is the correction that accounts for the spurious interaction caused by finite supercell size and periodic boundary conditions.

In the following, we will briefly introduce the four modules:

TSC module: As shown in Eq. (1), the defect (dopant) formation energy and thus the defect density are functions of the elemental chemical potentials, so the calculation of defect properties requires the elemental chemical potentials as inputs. The TSC module is for calculating the chemical potential ranges through considering the influences of all the competing secondary compounds that can limit the thermodynamic phase stability of the host compound semiconductor. In order to make the consideration of secondary compounds as complete as possible, the material genome database (Materials Project^[27]) is used to search for the competing secondary compounds. Through combining the formation energy information of all the secondary compounds with the calculated formation energy of the host compound (equivalent inputs are used for DFT calculations), the critical secondary compounds that limit the stable chemical potential region can be determined quickly.

DEC module: The DEC module can generate the maximumly-cubic supercell structures of the host semiconductor, based on which the defect and dopant supercell structures can be generated with all the possible defect/dopant types, atomic sites and structural configurations considered. Then ab-initio software, such as Vienna Ab-Initio Simulation Package (VASP)^[28] and Quantum ESPRESSO (QE)^[29], will be called to perform atomic structural relaxation, total energy and electronic structure calculations for the defect/dopant supercells. First the neutral charge state is calculated and then different charge states will be generated according to the defect/dopant levels in the band gap and calculated. Then the formation energies of defects and dopants in different charge states can be calculated using Eq. (1). The chemical potentials outputted by TSC are used as . The correction term that accounts for the spurious interaction caused by finite supercell size and periodic boundary conditions can also be calculated automatically based on the output of ab-initio calculations. The charge-state transition levels are then calculated to be the Fermi level at which the formation energies of a defect in two charge states q and q’ are equal, = .

DDC module: Based on the results of TSC and DEC, the DDC module solves the charge-neutrality equation self-consistently to determine the equilibrium defect/dopant densities, Fermi level, and carrier densities, as shown by the purple part of Fig. 2. For a given chemical potential condition, the formation energies of defects and dopants in different charge states are treated as input for a two-step self-consistent calculation, which solves the charge-neutrality condition equations at high growth temperature and low working (measuring) temperature to determine the Fermi level of the semiconductor samples. Then the equilibrium densities of defects and dopants and the densities of electron and hole carriers can be calculated under the given chemical potential condition. Through changing the chemical potentials, the defect/dopant and carrier densities can be calculated as functions of the chemical potential conditions, which can be used for the defect/dopant engineering based on the growth condition control.

CDC module: For the high-density defects or dopants identified by DDC module, the CDC module can calculate their excited-state carrier dynamics properties. The transition energy levels from the DEC module and the Fermi level from the DDC module will be used to screen for the transitions between the defect levels and the VBM or conduction band minimum (CBM) level that may produce the PL peaks or cause the fast non-radiative recombination of carriers. Then the phonon spectrum of the defect/dopant supercell, electron-phonon coupling matrix and transition dipole moments between the defect/dopant states and VBM or CBM states will be calculated, based on which the radiative and non-radiative transition rates (carrier capture cross sections) and the lineshape of photoluminescence spectra induced by defects and dopants can be calculated. Such calculations may be performed based on the single-phonon-mode one-dimensional configuration coordinate diagram^[30] or the all-phonon-modes three-dimensional schemes^[31]. Combining the calculated transition rates and the equilibrium defect and carrier densities, the CDC module may also calculate the Shockley–Read–Hall recombination rates and the non-equilibrium carrier lifetime.

The formation of defects or the doping of extrinsic elements in semiconductor lattice involves the exchange of atoms between the semiconductor lattice and the external environment, so the formation energies of defects and dopants in Eq. (1) depend on the elemental chemical potentials, which describe the abundance of the elements (partial pressure for the gas elements) in the environment. Although the abundance of a certain element can be controlled by changing the environment, the changes of the elemental chemical potentials are actually not unlimited, because they should satisfy a series of thermodynamic conditions to make the pure-phase crystalline semiconductor stable. Here we take the quaternary compound semiconductor Cu₂ZnSnS₄ as an example to show how to determine its stable elemental chemical potential range.

Firstly, at the equilibrium state of the compound Cu₂ZnSnS₄, the weighted sum of the chemical potentials of its component elements should be equal to the formation energy of the compound,

$ 2{\mu }_{\mathrm{C}\mathrm{u}}+{\mu }_{\mathrm{Z}\mathrm{n}}+{\mu }_{\mathrm{S}\mathrm{n}}+4{\mu }_{\mathrm{S}}=\mathrm{\Delta }{H}_{\mathrm{f}}\left({\mathrm{C}\mathrm{u}}_{2}{\mathrm{Z}\mathrm{n}\mathrm{S}\mathrm{n}\mathrm{S}}_{4}\right) , $  (2)

where (Cu₂ZnSnS₄) is the compound formation energy of Cu₂ZnSnS₄ relative to the Cu, Zn, Sn and S elemental phases. That means the chemical equilibrium state is arrived for the reaction 2Cu + Zn + Sn + 4S → Cu₂ZnSnS₄, so the product compound Cu₂ZnSnS₄ is stable.

Secondly, to make the synthesized sample be pure-phase Cu₂ZnSnS₄, the formation or coexistence of the competing secondary phases, such as the binary and ternary compounds CuS, Cu₂S, ZnS, SnS, SnS₂ and Cu₂SnS₃, and the elemental phases Cu, Zn, Sn, S, should be avoided. Therefore, the weighted sum of the chemical potentials of their component elements should be lower than their corresponding formation energies (the formation energies of elemental phases are 0), as described by the following inequations,

$ {\mu }_{\mathrm{C}\mathrm{u}}+{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left(\mathrm{C}\mathrm{u}\mathrm{S}\right), $  ()

$ 2{\mu }_{\mathrm{C}\mathrm{u}}+{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left({\mathrm{C}\mathrm{u}}_{2}\mathrm{S}\right), $  ()

$ {\mu }_{\mathrm{Z}\mathrm{n}}+{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left(\mathrm{Z}\mathrm{n}\mathrm{S}\right), $  ()

$ {\mu }_{\mathrm{S}\mathrm{n}}+{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left(\mathrm{S}\mathrm{n}\mathrm{S}\right), $  ()

$ {\mu }_{\mathrm{S}\mathrm{n}}+2{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left(\mathrm{S}\mathrm{n}{\mathrm{S}}_{2}\right), $  ()

$ 2{\mu }_{\mathrm{C}\mathrm{u}}+{\mu }_{\mathrm{S}\mathrm{n}}+3{\mu }_{\mathrm{S}} < \mathrm{\Delta }{H}_{\mathrm{f}}\left({\mathrm{C}\mathrm{u}}_{2}{\mathrm{S}\mathrm{n}\mathrm{S}}_{3}\right), $  ()

$ {\mu }_{\mathrm{C}\mathrm{u}} < 0, $  ()

$ {\mu }_{\mathrm{Z}\mathrm{n}} < 0, $  ()

$ {\mu }_{\mathrm{S}\mathrm{n}} < 0, $  ()

$ {\mu }_{\mathrm{S}} < 0. $  (3)

The allowed chemical potential range of Cu, Zn, Sn and S that stabilizes the pure-phase Cu₂ZnSnS₄ is limited by these equations and inequations. If the formation energies of the host and all the competing compounds are known, the ranges can be calculated and plotted in the chemical potential space, as shown in Ref. [4]. One example about SbSeI is also shown in Section 3.2.

The four elements Cu, Zn, Sn and S can form many binary, ternary and quaternary compounds, which can all act as the competing secondary phases of Cu₂ZnSnS₄. In principle, for the accurate calculation of the chemical potential range, one should take account of all the possible competing compounds, whose number can be large, especially for quaternary or multinary compounds. In DASP, the TSC module provides an automated solution for the quick screening of the critical competing phases and thus the accurate calculation of the chemical potential range. With the chemical formula of the compound, TSC will visit the materials genome database, such as the MP database^[27], to search for all the compounds that are composed of the component elements, and then download the total energy, formation energy and DFT calculation input files of these compounds. Meanwhile, TSC will also perform ab-initio DFT calculation for the formation energy of the host compound with the input consistent with the MP database. With the formation energies of the host compound and all the competing compounds, TSC can solve the thermodynamic constraint equations and inequations to determine the stable chemical potential region and the critical competing compounds that limit the stable region directly. For the host and the critical competing compounds, TSC can further calculate their formation energies with higher accuracy and different functionals, to ensure that the calculated range of the chemical potentials are accurate. Usually these critical competing compounds are also those determining the energy above the convex hull, which shows the thermodynamic stability of the compound with respect to the phase separation into other competing compounds, so TSC can also be used for high-throughput and accurate calculation of the energy above convex hull and thermodynamic stability of compounds.

In real semiconductor lattices, the densities of defects and dopants are usually much lower than the densities of atoms, e.g., there is only one defect or dopant among 10⁴–10⁷ atoms for defects or dopants with a density of 10¹⁵–10¹⁸ cm^–3. Therefore, the distance between two defects or dopants is usually very large. However, in the periodic supercell model, the distances between the defect and the periodic image defects in the neighboring supercell are actually not very large since the supercell has only several hundred atoms. Therefore, the supercell model causes unphysical interaction between defects or dopants, which may induce errors in the calculated defect properties.

To reduce the errors caused by the finite supercell size and achieve better supercell size convergence, DASP always tends to generate the supercells that are nearly cubic, so that the defect–defect (dopant–dopant) distance is maximized and the unphysical interaction is minimized. As illustrated in Fig. 3, for the two supercells with the same volume (same number of atoms), the smallest defect–defect distance in the nearly-cubic supercell should be larger than that in the supercell derived from direct expansion of primitive cell. In the nearly-cubic supercell, the smallest defect–defect distance should be along the (100), (010) or (001) direction, and the smallest distances along the three directions are similar. However, in the supercell derived from direct expansion of primitive cell, the smallest defect–defect distance may be along the (110) or (111) directions, and the distance can be smaller than the distance along the (100), (010) or (001) direction if the supercell basis vectors are not orthogonal and one of their angles is larger than 120°. Besides the advantage in maximizing the defect–defect distance, the nearly-cubic supercell also leads to the orthogonality of the reciprocal lattices, which may accelerate the ab-initio calculation and achieve better convergence.

The generation of the nearly-cubic supercell in DASP is through the maximumly-cubic supercell generation function, which firstly transforms the primitive cell into the nearly-orthogonal supercell and then expands it to get the nearly-cubic supercell which maximizes the cubic degree of the supercell with a given number of atoms.

To search for the nearly-orthogonal supercell, we developed a traversal searching method through trying different conversion matrices that can transform the lattice vectors of the primitive cell to the nearly-orthogonal lattice vectors. The conversion matrix should be non-singular because the lattice vectors are not coplanar. In addition, only swapping the lattice vectors will not change the shape of the cell, so there is no need to consider the exchange operation during the transformation. Therefore, the matrix satisfies the LU decomposition condition, and can be written as,

$ C=\left($$\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{21}& c_{22}& c_{23}\\ c_{31}& c_{32}& c_{33}\end{array}$$\right)=LU=\left($$\begin{array}{ccc}1& 0& 0\\ l_{21}& 1& 0\\ l_{31}& l_{32}& 1\end{array}$$\right)\left($$\begin{array}{ccc}{u}_{11}& u_{12}& u_{13}\\ 0& u_{22}& u_{23}\\ 0& 0& u_{33}\end{array}$$\right), $  (4)

where , , in the upper triangular matrix describe the expansion operations, which are the repeating numbers along the three lattice vectors of the primitive cell. , , in the lower triangular matrix and , , in the upper triangular matrix describe the linear combination operations. By traversing different conversion matrices, we can get different nearly-orthogonal cells. Taking into account both the computational cost and accuracy, we only traverse the conversion matrices that meet the following conditions,

$ {c}_{ij}=0,    \pm 1, \pm 2, \pm 3, \pm 4\left(i,j=1,      2, 3\right),{u}_{kk}=1,    2, 3, 4 \left(k=1, 2, 3\right), $  (5)

where is the element in the conversion matrix , is the diagonal element in the upper triangular matrix .

With the nearly-orthogonal supercell, we can just expand it directly to get nearly-cubic supercell. The cubic degree is described by the ratio , where V_supercell is the volume of the supercell, V_cube is the volume of the cube that can envelop the whole supercell. If the degree is 1, the shape of the supercell is ideally cubic. Considering the computational cost, the number ( ) of atoms in the supercell is limited, e.g., 200–1000 atoms. Although very large can give very cubic supercell and thus tend to increase the cubic degree to 1, the cubic degree does not necessarily increase monotonically as increases from 100 to 1000. Therefore, when selecting a supercell for defect calculation, DASP considers both the cubic degree and in order to maximize the defect–defect distance under the constraint of the limited . An empirical quantity β is defined to balance the two factors and used as the criterion for choosing supercells,

$ \beta =\frac{{V}_{\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}}}{{V}_{\mathrm{c}\mathrm{u}\mathrm{b}\mathrm{e}}}+0.001{N}_{\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}}. $  (6)

The maximumly-cubic supercell is generated by DASP through maximizing the quantity β in a given range of .

After the supercell is generated, DASP can generate structures of defects and dopants in the primitive-cell region of the supercell. The common defect types will be considered automatically. For example, the intrinsic defects of the quaternary compound Cu₂ZnSnS₄ that DASP considered include the Cu, Zn, Sn and S vacancies (V_Cu, V_Zn, V_Sn and V_S), interstitials (Cu_i, Zn_i, Sn_i, S_i) and antisites (Cu_Zn, Cu_Sn, Cu_S, Zn_Cu, Zn_Sn, Zn_S, Sn_Cu, Sn_Zn, Sn_S, S_Cu, S_Zn, S_Sn). Meanwhile, for the low-energy donor and acceptor defects, they can form donor–acceptor compensated defect complexes, which can also be considered by DASP.

For low-symmetry semiconductors, the structure may have several non-equivalent atomic sites for one element, and the same-type defects on non-equivalent sites may have different properties. When generating vacancy, interstitial and antisite defects, DASP will consider all non-equivalent sites. For interstitial defects, DASP will search for the largest void region in the structure and meanwhile consider the Coulomb repulsion to determine the possible interstitial sites of cations and anions. For low-energy interstitial defects, DASP will also generate 10–20 different configurations randomly in the primitive-cell region of the supercell (different interstitial sites are ensured to be not close to each other).

For low-energy vacancy and antisite defects, there may be other distorted or metastable structure configurations, such as the DX centers^[32], which originate from large structural distortion of the antisite donor defects but act as acceptors after distortion. These distorted or metastable defects have been shown to have important impact on the electrical properties of semiconductors^[33-35]. They can be obtained by imposing a structural perturbation to the locally relaxed structure of the defect in the charge state q, because the perturbation may overcome the structural transition barrier to arrive at a new distorted structure.

DASP adds two types of structural perturbations: (i) displacing the atoms within a sphere around the defect by less than 0.5 Å randomly, as shown in Fig. 4(b); (ii) randomly distributing all the atoms within the sphere, as shown in Fig. 4(c). The default value of the sphere radius is set to 3 Å, but will be automatically increased to ensure at least 4 atoms in the sphere. The total number of the structures generated with the perturbation will be in the range of 10–20, but some of them may become identical after structural relaxation. Such kind of structural perturbations is in fact for achieving a global structural relaxation for the defect.

Ionized defects in different charge states are fundamental to the basic understanding of defect physics, but there are still several open questions during the calculation of the properties of charged defects. The first one is how to determine the range of defect charge q for an ionized defect. In the past, the octet rule was often selected as a criterion to judge the charge range of a defect, which simply adopts the nominal charge of the single element. For instance, the sulfur vacancy should take 0, 1+, 2+ charges according to the simple octet rule, however, in non-conventional semiconductors such as MoS₂, neither 1+ nor 2+ state of the sulfur vacancy is stable while the 1– charge state exists^[36].

To solve this problem, the estimation procedure of the charge range implemented in DASP is based on the calculated eigenvalues of the neutral defects at Γ point. Therefore, in the first step of DEC module, DASP will generate the structures of point defects in their neutral state, and then call ab-initio software to perform atomic relaxations and static self-consistent calculations. Once all the calculations are done, the eigenvalues of bulk and all the neutral defects will be extracted. The hybrid functional (such as HSE) is recommended for the static self-consistent calculations to obtain the more reasonable band gap and more reasonable location of the defect levels. As shown in Fig. 5, for defect A, if an occupied and an unoccupied level are found within the band gap between the VBM and CBM levels of the bulk, these two levels will be taken as defect levels and the defect A in 1+ and 1– charges will be calculated subsequently. Following this scheme, when three occupied levels and three unoccupied levels are found in the bulk band gap, the defect C in [3+, 2+, 1+, 1–, 2–, 3–] charge states will all be calculated.

In Eq. (1), the eigenvalue of bulk VBM is used as the reference of the Fermi level, i.e., E_F = 0 means the Fermi level is located at VBM. However, the eigenvalues in the calculation of bulk and defect supercells can be compared only when the electrostatic potentials of the two periodic supercells are aligned^[12]. This can be accomplished either by aligning the electrostatic potential far from the defect (e.g., the potential written in the LOCPOT file of VASP), or by aligning the deep core level of the farthest atom from the defect. The alignment term can be written as

$ \Delta {V}_{q/b}=V\left(\alpha ,q\right){|}_{\mathrm{f}\mathrm{a}\mathrm{r}}-V\left(\mathrm{b}\mathrm{u}\mathrm{l}\mathrm{k}\right){|}_{\mathrm{f}\mathrm{a}\mathrm{r}} . $  (7)

The correction of potential alignment to the formation energy of a defect in the charge state q is , and is automatically included in the term .

The spurious Coulomb interaction between charged defect and its periodic images, as well as charged defect and the neutralizing background charge should be corrected from the calculated formation energies of ionized charged defects, which are called image charge corrections. DASP can currently perform the corrections according to two schemes. The first one is the Lany–Zunger scheme^[12], which can be represented by

$ {\Delta E}_{\rm{LZ}}=[1-{c}_{\rm{sh}}(1-1/\epsilon \left)\right]\frac{{q}^{2}\alpha }{2\epsilon L}, $  (8)

where α is the lattice-dependent Madelung constant, ε is the static dielectric constant, and L is the linear supercell dimension. Since DASP will always tend to generate a nearly-cubic supercell, the default value of is set to –1/3^{[12, 22]}. Based on these definitions, the total correction term that appears in Eq. (1) can be separated into

$ {E}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}=q\Delta {V}_{q/b}+{\Delta E}_{\rm{LZ}} .$  (9)

DASP also provides an interface with the code sxdefectalign written by Freysoldt, which implements the Freysoldt–Neugebauer–Van de Walle (FNV) scheme^[23] for charged defect correction. Different from the original paper that uses the planar average of the electrostatic potential to obtain the potential offset, we here use the atomic site potential instead, as Kumagai and Oba proposed^[37]. In this scheme, the correction can be written as

$ {E}_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}}={E}_{\mathrm{l}\mathrm{a}\mathrm{t}}-q(\Delta {V}_{q/b}-{V}_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}) , $  (10)

where is the Madelung energy and is the potential obtained from the chosen model (such as Gaussian) charge density. Note that in the FNV scheme, the potential alignment defined in Sec. 2.3.4 is already included.

The two correction schemes mentioned above may be not valid for the charged defects in low-dimensional layered semiconductors, such as those in monolayer MoS₂ or WSe₂. A series of new correction schemes have been successfully developed in the past decade for the charged defect correction in layered semiconductors^[38-43], which will be implemented in the future version of DASP.

The formation of a defect or dopant causes the change of the electronic wavefunction and the redistribution of charge density in the semiconductor lattice, which produces new electronic states (defect states) in the band gap. Although the wavefunction and charge density near the defect site can be changed dramatically, the area far from the defect should be weakly influenced. Therefore, for different defects or dopants in different charge states, the wavefunction and charge density in the region far from the defect site should be similar to those of the defect-free supercell. Since the self-consistent calculations of the wavefunction and charge density will be repeated hundreds of times, it will save a large amount of computational cost if we take advantage of the slightly changed wavefunction and charge density in the region far from the defect site to reduce the number of self-consistent calculation steps.

In DASP, the initial charge density of the neutral defect/dopant supercell is generated based on the converged wavefunction and charge density of the bulk host supercell, i.e., the charge density in the region far from the defect site is the same as that of the bulk host, while the charge density in the sphere around the defect is generated by a superposition of atomic charge densities and the charge density in the interface region is smoothed and rescaled according to the number of valence electrons of the defect supercell. Then the self-consistent calculations are performed to obtain the converged charge density and wavefunction of the neutral defect/dopant supercell. For the charged defect/dopant, its initial charge density is generated based on the converged wavefunction and charge density of the neutral defect through changing the occupation number of the defect states. With the initial charge density, the ab-initio self-consistent calculations can reach convergence with much less steps. The saving of computational cost can be more significant when large supercells are used.

Under thermodynamic equilibrium, the density of a defect α in the charge state q is determined by its formation energy according to^{[25, 44]}

$ n(\alpha ,q)={N}_{\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{s}}{g}_{\mathrm{q}}{\mathrm{e}}^{(-{\Delta E}_{\mathrm{f}}/{k}_{\mathrm{B}}T)} , $  (11)

where N_sites is the density of the possible defect sites, g_q is the charge-dependent degeneracy factor, ΔE_f is the defect formation energy at the given Fermi level and elemental chemical potentials as described by Eq. (1). All the ionized defects in the charge stateq≠ 0 produce carriers. The positively charged donor defects with q > 0 produce electron carriers, and their summed charge is 0}\left[q\cdot n\left(\alpha ,q\right)\right]$?>; while the negatively charged acceptor defects with q < 0 produce hole carriers, and their summed charge is . The final densities of electron and hole carriers are contributed by both the thermal excitation and the ionization of all these defects (dopants). The equilibrium-state Fermi level can be calculated through solving the charge neutrality equation,

$ {n}_{0}+\sum _{\alpha ;q < 0}\left[(-q)\cdot n\left(\alpha ,q\right)\right]={p}_{0}+\sum _{\alpha ;q > 0}\left[q\cdot n\left(\alpha ,q\right)\right] , $  (12)

where and 0}\left[q\cdot n\left(\alpha ,q\right)\right] $?> are the summed charges of negatively charged defects and positively charged defects, weighted by the charge q. n₀ and p₀ are free carrier densities, which can be defined as

$ {n}_{0}={\int }_{{\epsilon }_{\mathrm{C}\mathrm{B}\mathrm{M}}}^{+\infty }{g}_{\mathrm{C}}\left(E\right)f\left(E\right)\mathrm{d}E, $  (13)

$ {p}_{0}={\int }_{-\infty }^{{\epsilon }_{\mathrm{V}\mathrm{B}\mathrm{M}}}{g}_{\mathrm{V}}\left(E\right)(1-f(E\left)\right)\mathrm{d}E, $  (14)

where and are the density of states (DOS) for the conduction and valence bands, respectively, and is the Fermi-Dirac occupation function. and can be calculated from both the parabolic band approximation and the exact integration of the ab-initio calculated band structure. The default setting in DASP is from the parabolic band approximation that requires the calculation of electron and hole effective masses, but it can be changed to the exact DOS from ab-initio calculations.

The semiconductors are usually grown or synthesized at a high temperature and then go through a rapid annealing process to a lower working (measuring) temperature. Therefore, the defects are usually formed at the high temperature and then the densities of different charge states will redistribute during the rapid annealing. The DDC module is developed in accordance with such fabrication process^[44-46], as shown schematically in Fig. 2. Eq. (12) is firstly solved at a high growth temperature, and then the Fermi level and the densities of each defect in different charge state q can be obtained at the high temperature. Afterwards, Eq. (12) is solved again at the lower measuring temperature, but now the densities of defects in different charge states do not follow Eq. (11). In the second step, the density summation for each defect over all charge states is fixed at the value calculated in the first step, and the density of each charge state will undergo a redistribution according to the Fermi–Dirac occupation. Then a new Fermi level can be obtained at the measuring temperature, and the redistributed defect densities and carrier densities can be calculated.

The calculation of defect and carrier densities using DDC module in DASP is quite easy. After finishing the calculation of TSC and DEC, DDC can be calculated based on the output files of TSC and DEC. The defect (dopant) densities, carrier densities and Fermi level can be plotted or saved as functions of the elemental chemical potentials and growth/measuring temperatures.

After calculating the defect densities using DDC module, one may identify a portion of high density defects (dopants), which may be critical to the optical and electrical properties of the host semiconductor. For those important defects, we implement in the CDC module three functions for studying the excited-state carrier dynamics properties based on the phonon spectrum and electron-phonon coupling calculation: (i) photoluminescence (PL) lineshape of defects; (ii) radiative carrier capture coefficient of defects; (iii) phonon-assisted nonradiative carrier capture coefficient (cross section) of defects. In the following, we will mainly introduce how to calculate the PL lineshape, while the details of the calculation on carrier capture will not be discussed here. More details can be found in Refs. [30, 31].

The photoluminescence can be caused by the transition of carriers between the defect state and the VBM or CBM state, as shown in Fig. 6(a). Its intensity at finite temperature is mainly determined by the transition dipole moment between two states and the lineshape function, which is written by^{[47, 48]},

$$$\begin{array}{l}I\left(\hslash \omega \right)={\displaystyle \frac{e^2{n}_{\mathrm{r}}{\omega }^3}{3{ϵ}_0\pi c^3\hslash }}{|⟨{\psi }_{\mathrm{i}}|\hat{\mathit{r}}|{\psi }_{\mathrm{f}}⟩|}^2\sum _m{p}_m\sum _n{|⟨{\chi }_{\mathrm{i}\mathrm{m}}|{\chi }_{\mathrm{f}\mathrm{n}}⟩|}^2\\ \times \delta (E_{\mathrm{Z}\mathrm{P}\mathrm{L}}+\hslash {\omega }_{\mathrm{i}\mathrm{m}}-\hslash {\omega }_{\mathrm{f}\mathrm{n}}-\hslash \omega ),\end{array}$$$  (15)

where e is the elementary charge, n_r is the refractive index of the bulk material, ɛ₀ is the vacuum permittivity, ћω is the corresponding energy of the emitted photon. is the transition dipole moment between the band edge state and defect state calculated at Γ point; p_m is the Boltzmann occupation of the initial vibrational state m changing with temperature; χ_im and χ_fn are the vibrational wavefunctions of the initial and final states, and ћω_im and ћω_fn are the corresponding eigenvalues. E_ZPL is corresponding to the charge-state transition level calculated in DEC module. At T = 0 K, all the charge carriers fall into the lowest-energy vibrational mode of initial state (m = 0), since there is no thermal activation for excited-state carriers. As a result, one can simply omit the summation over m in Eq. (15) to obtain the result for T = 0 K.

In order to evaluate the overlap integral in Eq. (15), one-dimensional configuration coordinate diagram is used with an example given in Fig. 6(b), which adopts a single effective phonon mode to represent all 3N phonons in the system^[49]. This phonon mode can be depicted by the distortion of the defects between two charge states; the structural difference ΔQ in the configuration coordinate can be given by^[50]

$ \mathrm{\Delta }Q=\sqrt{\sum _{\alpha }{m}_{\alpha }{({\boldsymbol{R}}_{\rm{i}\alpha }-{\boldsymbol{R}}_{\rm{f}\alpha })}^{2}}, $  (16)

where R_iα and R_fα are the Cartesian coordinates of the atom α in initial and final supercell structures, and m_α is its mass. Using configuration coordinate diagram, the effective phonon frequencies and can be obtained by fitting the potential energy surfaces, and the integral can also be numerically calculated.

Configuration coordinate diagram is also useful for further calculating the radiative and nonradiative carrier capture coefficient of such process, which requires the exact calculation of transition dipole moment and electron-phonon coupling matrix element. In the current version, DASP only supports the one-dimensional static coupling method for calculating nonradiative carrier capture^[30]; while the original static coupling theory including all 3N phonon modes will be implemented in the future version^{[31, 51]}.

In the following, we will show three examples about the application of DASP in the calculation of the defect and dopant properties of semiconductors, including the intrinsic point defect properties of the benchmark system GaN, the intrinsic point defect properties of the low-symmetry quasi-one-dimensional SbSeI, and the PL spectrum of C-doped GaN.

GaN is a well-studied wide-band-gap semiconductor whose intrinsic defects have been studied by many groups since 1990s, so it is a good benchmark system for testing our DASP software. In Fig. 7, we show the calculated formation energies of vacancies, antisites and interstitials in GaN as functions of Fermi level, which are well consistent with the published results calculated using the hybrid functional^[52-54].

The formation energies of these point defects in neutral states are relatively high in GaN, no matter under Ga-rich or N-rich condition. Among them, the donor defect, N vacancy V_N, has the lowest formation energy, but it still cannot produce a high density of electron carriers or good n-type conductivity and shift the Fermi level to close to the CBM level, even under the Ga-rich condition. Therefore, the formation of intrinsic point defects should not produce good n-type conductivity in GaN. For V_N, the DASP calculations based on the spin-polarized and hybrid functional ab-initio calculations show a (+/3+) transition level lying 0.4 eV above the VBM, and a (0/+) level close to the CBM. The (+/3+) transition level was usually not found in the early calculations based on the LDA or GGA to the exchange correlation functional^{[8, 56]}, but was found in the recent hybrid functional calculations^[52-54].

Antimony selenoiodide (SbSeI) has a quasi-one-dimensional structure, as shown in Fig. 8(a), which is similar to that of Sb₂Se₃^[24]. Its band structure is shown in Fig. 8(b), and has an indirect band gap of 1.79 eV. Its VBM is composed of the hybridized states of Sb-5s, Se-4p and I-5p, while the CBM is mainly Sb-5p. Using the TSC module, two secondary compounds, SbI₃ and Sb₂Se₃, are found to determine the allowed range of the Sb, Se and I elemental chemical potentials, together with the elemental phases of Sb, Se and I. Figs. 8(c) and 8(d) present the 3D and projected-2D phase stability diagrams in the elemental chemical potential space based on the results of the TSC calculations. The ternary compound SbSeI is stable only in the orange region surrounded by the competing phases Sb₂Se₃, SbI₃, Sb and Se. The four boundary points, A, B, C and D, are selected as the chemical potential conditions for the following DEC and DDC calculations.

Fig. 9 shows the calculated formation energies of defects in SbSeI as functions of Fermi level for the four chemical potential conditions, A, B, C and D. The densities of these defects and carriers are also calculated using the DDC module and the results for the high-density defects are plotted in Fig. 10(a) as functions of the elemental chemical potential. Meanwhile, the change of Fermi level and carrier densities with the elemental chemical potential conditions is also plotted in Fig. 10(a).

It is obvious in Fig. 9 that there are many intrinsic defects with formation energies lower than 2 eV, so the formation of defects should be energetically easy in the quasi-one-dimensional lattice of SbSeI. Among them, the antisite defects between two anions, Se_I and I_Se, should be the dominant defects with the lowest energy and highest density. Se_I is an acceptor with a (0/–) level at 0.34 eV above VBM, which is not deep. I_Se is a donor with a (0/+) level at 0.17 eV below CBM, which is relatively shallow. Figs. 10(b) and 10(c) plot the squared wavefunction of the Se_I and I_Se defect states, which show that the defect states are actually localized, although their corresponding transition levels are shallow, indicating the supercell size is sufficient for describing these defects in SbSeI. Besides them, there are also several other defects with formation energies lower than 1 eV and deep transition levels in the band gap, e.g., Se_Sb and V_I, which are both amphoteric with deep (0/+) and (0/–) transition levels.

Although the two dominant defects, the Se_I acceptor and the I_Se donor, have quite high density (10¹⁷–10¹⁸ cm^-3) and relatively shallow transition levels (should produce high density of electron and hole carriers), the final equilibrium density of carriers is not high due to the donor–acceptor compensation. Therefore, the Fermi level of the SbSeI sample with high densities of defects is always located near the middle of the band gap, i.e., at least 0.4 eV far from the CBM or VBM level as shown in Fig. 11(a), and the electrical conductivity can only be weakly n-type or weakly p-type, regardless of the elemental chemical potential.

In contrast with the shallow Se_I and I_Se whose densities are always high, the densities of deep-level defects Se_Sb and V_I change significantly as the elemental chemical potentials change. When the chemical potential points are close to the A-D line, the density of V_I is very high; when the points are close to B-C, its density decreases to lower than 10¹⁵ cm^-3. For Se_Sb, its density is very high when the chemical potential point is near the C point, but decreases quickly as the point moves away from the C point. These results suggest that the chemical potential should be controlled at the intermediate region of the B-C path, in order to suppress the formation of deep-level defects such as Se_Sb and V_I.

As we can see, Figs. 8–10 give a good example about the application of DASP for calculating the elemental chemical potential region that stabilizes the compound semiconductors, then the formation energies and transition levels of all intrinsic point defects, and finally the densities of defects and electron/hole carriers. For new semiconductors, the similar calculations should be performed to have a complete picture about their defect and dopant properties.

The PL spectrum is a widely used optical characterization method of defects in semiconductors. Usually the origin of the PL peaks is attributed to different defects according to the energy differences between the defect level and band edge levels. However, such kind of attribution is questionable, especially when there are several defects whose energy levels are nearly degenerate. DASP can calculate the PL lineshape of different defects, which provides an extra criterion for the identification of the defect origin of PL peaks.

Here we show an example about the calculated PL lineshape of the C_N dopant in GaN. Using DEC module, we find C_N has a (0/+) transition level of 0.39 eV, and a (0/–) transition level of 1.07 eV above the VBM level, which agree with the calculations of Lyons et al.^[57]. Therefore, the photo-excited electron carrier can be captured by the (0/–) level of C_N through a radiative transition, from the initial state of neutral C_N (C_N⁰) with an electron on the CBM to the final state of negatively charged C_N (C_N^–) with the electron on the (0/–) defect level. This transition is plotted in the configuration coordinate diagram in Fig. 11(b). The direct optical transition level is calculated to be 2.11 eV. The PL lineshape of the transition calculated using the CDC module is plotted in Fig. 11(c), which is in good agreement with experimental measurement and previous calculated result^[50].

Further CDC calculation shows the radiative carrier capture coefficient C_n is 1.7 × 10^–13 cm³/s, slightly larger than the previous calculated result of 0.7 × 10^–13 cm³/s^[58]. The main reason is attributed to the larger momentum matrix element that we calculated, which is |p_if|² = 0.06 a.u. between the initial and final states.

We developed a software DASP composed of four modules, TSC, DEC, DDC and CDC, for the automated calculation of defect and dopant properties in the crystalline semiconductors (insulators). DASP just needs the input of the crystal structure file of the semiconductor, then it can visit the materials genome database and call the ab-initio DFT software such as VASP to calculate the defect and dopant properties, including, (i) the chemical potential range of component elements that stabilizes the compound semiconductor (a descriptor of the thermodynamic stability of the compound) and the highest allowed chemical potential of dopant elements; (ii) the formation energies of defects and dopants in different charge states, as functions of elemental chemical potentials and Fermi level, and their charge-state transition levels; (iii) the equilibrium densities of defects and dopants, Fermi level, densities of electron and hole carriers, as functions of elemental chemical potentials at growth and working temperature; and (iv) the carrier capture cross sections, radiative and non-radiative carrier recombination rate and lifetime, and the PL spectra. DASP is designed to be an automatic toolbox for the theoretical calculation studies on defects and dopants in semiconductors, and can be used not only for interpreting the electrical and optical characterization experiments of defects and dopants, but also for the quantitative defect and dopant engineering in functional semiconductor devices.

We thank Profs. Suhuai Wei, Xingao Gong, Aron Walsh, Linwang Wang and Yuning Wu, and Drs. Zhenkun Yuan, Jiqiang Li, Zenghua Cai and Tao Zhang for their long-term collaboration and very helpful discussion. The development of the commercial version of DASP is supported by the joint project between Hongzhiwei Technology (Shanghai) Co., Ltd. and Fudan University.