2 Modules

2.1 PREPARE Module

There are two main input files in DASP: POSCAR and dasp.in , where dasp.in is the main parameter control file.

Before calculation, PREPARE module would check if the input parameters in the dasp.in are reasonable first. If it is, the next calculation begins and marks out that the preparation work has been done in the 1prepare.out .

If it isn’t, the information of mistakes and warnings will be written in the 1prepare.out and program terminated. Users need to revise the input parameters based on the warnings and restart program.

Since DASP calls VASP to carry out the structure and electronic structure calculation, thus DASP will automatically generate the input files required for VASP calculation according to the parameters in dasp.in : POSCAR (crystal structure), INCAR (calculation parameters), KPOINTS (k-points setting), POTCAR (pseudo-potential file), and submission script.

• POSCAR

  According to the POSCAR that user input, DASP will generate a supercell that is nearly cubic within a given atoms number through the self-developed algorithm “nearly cubic supercell”. Then, fix the lattice constant and optimize the positions of all atoms in the supercell to obtain POSCAR_final file.

• KPOINTS

  For the calculations of supercells with defects, the single k point is used in DASP, namely only Gamma point is included in the KPOINTS file.

• INCAR

  DASP will generate two INCAR files, one is INCAR-relax for the structural optimization, ant the other is INCAR-static used in static calculation. The default generated is commonly used parameters. They can be modified after the PREPARE module runs if users need to modify them. According to the parameters in the dasp.in , DASP can adopt three different calculation levels (see input parameter level=1/2/3 for details). For level=2 or 3, the hybrid functional will be used. DASP can automatically determine the proportion of exchanged parts in hybrid functional by the experimental band gap that user set in dasp.in , and the matched proportion is written into INCAR.

• POTCAR

  According to the path of pseudo-potential file provided by the user, DASP will automatically generate the POTCAR required for calculation.

• Submission script

  DASP will automatically generate a submission script file based on the system name, queue name, nodes, cores and the path of VASP written by the user in dasp.in .

As mentioned above, PREPARE module will automatically call VASP to optimize the structure and determine the proportion of the exchange part in the hybrid functional. If user specifies the Lany-Zunger correction scheme in dasp.in (see correction in the introduction of input parameter for details), this scheme needs to know the Madelung constant corresponding to the supercell. According to the generated supercell, the PREPARE module will call VASP to automatically calculate the Madelung constant and write it into dasp.in . The Madelung constant here is only related to the lattice of the supercell.

All details of the above generated files and calculation process will be written into 1prepare.out . The information on status, results and errors can be queried in the file. After this module runs successfully, it will be noted the complete information at the end of the 1prepare.out , which will be detected when the subsequent TSC module starts.

2.2 TSC Module

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 which the chemical potential will be used as the input for the subsequent calculation on DEC and DDC modules.

Whether a compound semiconductor is stable depends on its formation energy. If it has more advantages in formation energy than its competitive pure phase and hetero-phase compounds, the compound can be thermodynamically stable. The judgment process is as follows.

• Judgement of thermodynamic stability of semiconductors

  A semiconductor if it is thermodynamically stable (does not decompose into simple substances or other hetero-phase, and can synthesize pure phase sample) needs to match the following three conditions:

（1）The formation of the target compound reaches thermodynamic equilibrium:

Under the equilibrium condition, the formation and decomposition of the target compound are in dynamic equilibrium. For a compound \(A_kB_lC_mD_n\) , the weighted sum of the chemical potentials \(\mu\) of its component elements should be equal to the formation energy \(E^f\) of the compound,

(1)\[k\mu_A + l\mu_B + m\mu_C + n\mu_D = E^f(A_kB_lC_mD_n)\]

（2）The formation of the competing secondary compounds concerning the host compound cannot be carried out:

For any hetero-phase compound \(A_{k'}B_{l'}D_{n'}\) , the chemical potential of each element and the formation energy of hetero-phase compound need to meet the following inequality:

(2)\[k'\mu_A + l'\mu_B + n'\mu_D < E^f(A_{k'}B_{l'}D_{n'})\]

（3）The simple substances of the constituent elements of the host compound will not form:

In order to avoid the formation of a simple substances, the chemical potential of each element satisfies the following inequality:

(3)\[\mu_A< 0\ ,\ \mu_B< 0\ ,\ \mu_C< 0\ ,\ \mu_D < 0\]

After considering all competing secondary compounds, the host compound can be thermodynamically stable if the chemical potential of each element can exist in a range that meets all the above constraints. Otherwise, the compound semiconductor is unstable and will be decomposed into simple substances or other hetero-phases, and the synthesized sample cannot guarantee a pure single-phase.

The result of thermodynamic stability given by the above conditions and processes actually has an equivalent effect to the energy above hull of compound in judging thermodynamic stability. If the energy above hull is positive, the compound is unstable, and there is no range of element chemical potential, which meets all the above conditions. If the energy above hull is negative, the compound is stable, and there is a range of element chemical potential, which meets all the above conditions.

• Two steps of thermodynamic stability calculation

  According to the above discussions, the formation energy of all competing secondary compounds is required to calculate the thermodynamic stability and the stable range of element chemical potential of compound semiconductors. For binary, ternary, quaternary, quinary and even more compounds, there are many possible competing secondary compounds, which all need to be taken into account in the calculation. If not fully considered, some unstable compounds will be predicted to be stable. Therefore, full consideration of all possible secondary compounds is critical to the accurate calculation of thermodynamic stability and elemental chemical potential range. In order to consider all secondary compounds as fully as possible, DASP will visit the Materials Project (MP) database to search for all the compounds that are composed of the component elements, and quickly determines the critical competing compounds by their formation energy. Then, for the host and critical competing compounds, TSC can calculate their formation energies with higher accuracy to ensure that the calculated range of the chemical potentials is accurate. It is divided into two steps:

• First step

  DASP will visit the materials genome database, such as Materials Project (MP) database, to search for all the compounds that are composed of the component elements, then obtained the total energy and structure of these compounds, and generate VASP calculation input files in which the parameters consistent with the MP database: INCAR , KPOINTS , POTCAR , POSCAR (copied from the files provided by user). Meanwhile, TSC will also perform a VASP calculation for the total energy and formation energy of the host compound (the calculated energy can be compared to that in MP database directly). The process of calculation is the same as that of the MP database, which is divided into twice structural optimization and one static calculation, which are relaxation1 , relaxation2 , and static under “TSC/ host homonymous directory”.

  With the formation energies of the host compound and all the competing compounds, TSC can solve the thermodynamic constraint equations and inequations to predict whether the host compound is stable, and determine the critical competing compounds that limit the stable chemical potential region. The structure and total energy of the primitive cell of the host compound are only needed to calculate at this stage, and GGA-PBE exchange-correlation is used leading to a small amount of calculation. A large number of data of secondary compounds are directly from the MP database without calculation. Therefore, all the competing secondary compounds can be considered quickly and fully to determine the critical competing compounds.

• Second step

  For the host and critical competing compounds, the unified VASP parameters and input files generated by the PREPARE module are adopted, INCAR and POTCAR , and KPOINTS and POSCAR (for the host compound, copied from the file provided by the user, for the critical competing compounds, download from the MP database) automatically generated by the TSC module, the energy and formation energy of the host and critical competing compounds are recalculated. In order to make the calculation quickly, only static calculations are made for the host and critical competing compounds, and the directory is located in “TSC/ host or secondary compounds ” static_recalc . Then according to the formation energies of the host and critical competing compounds, TSC resolves the thermodynamic constraint equations and inequations to calculate the stable chemical potential region of each element as the input of the subsequent DEC and DDC modules.

Based on the formation energies of the host and critical competing compounds, the TSC module also automatically calculates and outputs the energy above hull (ev/atom) and the most possible decomposition path of the host compound. The energy above hull also shows the thermodynamic stability of the compound with respect to the phase separation into other competing compounds. The more negative the value, the more stable the compound is, and the more positive, the more unstable it is conversely. So TSC can also be used for high-throughput and accurate calculation of the thermodynamic stability of new compounds.

For the ternary and quaternary compound semiconductors, TSC module can draw two-dimensional and three-dimensional stable region phase diagrams in the chemical potential space. Please refer to the Examples for details.

All details of TSC module operation, results and errors are output in 2tsc.out , the status of this module can be queried in this file. After this module runs successfully, it will be noted the complete information at the end of the 2tsc.out , which will be detected when the subsequent DEC module starts.

2.3 DEC Module

According to the supercell model, the formation energy and transition energy level (ionization energy level) of defects and dopants can be calculated by calling ab-initio software in the DEC module. This is also the main result since the 1990s given in most first principle research about defects and dopants.

The DEC module will generate a series of configurations with defects and dopants based on the parameters set by the user in dasp.in and the supercell produced by PREPARE module. Then, based on the VASP input file generated by the PREPARE module and the chemical potential calculated by the TSC module, the DEC module will call ab-initio software VASP to calculate the structure and electronic structure of defects and dopants. According to the result of calculations, the formation energy and transition energy level of defects and dopants will be calculated, and output the figure.

Before generating the structure, the DEC module will detect the parameters in dasp.in whether reasonable. If it is, start the next calculation, and output the running status information in real-time in 3dec.out . If it isn’t, relevant error and warning information will be output to 3dec.out file, and the program terminates. Users need to modify relevant parameters according to the error information and restart.

The workflow of DEC module is divided into the following six steps: generating neutral defects, automatically submitting tasks to calculate neutral defects, generating charged defects, submitting tasks to calculate charged defects, calculating the defect formation energy with different charge states under different chemical potentials and Fermi levels, and outputting figures.

• generating neutral defects

  The neutral defects include vacancy, antisite, and interstitial. The supercell structure POSCAR_final is generated by using PREPARE module, the vacancy and antisite on inequivalent sites are produced based on the crystal symmetry, meanwhile, the interstitial is obtained by randomly scattering at the position far away from atoms. After the defect configuration is generated, the DEC module will put the VASP input files, such as INCAR, KPOINTS, POTCAR, and submission script, generated by the PREPARE module into the directory of each defect. For the dopants, there are only two configurations: antisite (substitution) and interstitial.

• generating charged defects

  The DEC module will create the corresponding calculation directory of charged defects according to the results of neutral defects (the occupied state of the eigenvalue of neutral defects). If the calculation of neutral defects fails or does not converge, there will be no charged defects generated.

• automatically submitting jobs

  DEC module will automatically call VASP to carry out structural optimization and static calculation for all neutral and charged defects and detect whether the calculation is successful and converged. The command dec-state used under the dec directory can view the status of all jobs to be calculated at any time, including converged, does not converge, error, running, queuing, not submitted, etc. After the DEC module runs, the user can enter the corresponding directory to modify the INCAR for tasks that do not converge and have errors, and write the path of this directory into redo.in under the dec directory and re-execute the DEC module.

• formation energy calculation

  The formation energy of a point defect in the charge state q can be calculated as,

  (4)\[ΔE_{f} = E_{tot(defect)} - E_{tot(host)} - \Sigma_{i}n_{i}(\mu_{i} + E_{i}) + q(E_{F} + E_{VBM}) + E_{corr}\]

  where \(ΔE_{f}\) is the formation energy, \(E_{tot(defect)}\) and \(E_{tot(host)}\) are the total energies of the supercells with and without a defect (dopant). \(n_{i}\) is the number of atoms of element \(i\) remove from ( \(n_{i}\) ＜0) or added to ( \(n_{i}\) >0) the supercell to form the defect, and \(\mu_{i}\) is the elemental chemical potential referenced to the total energy \(E_{i}\) of the pure solid/gas elementary phase. \(q\) is the number of electrons, transferred from the supercell to the reservoirs in forming the defect cell. \(E_{F}\) is the Fermi level referenced to the eigenvalue of the valence band maximum (VBM) level of the bulk supercell. \(E_{corr}\) is the correction that accounts for the spurious interaction caused by finite supercell size and periodic boundary conditions.

  The DEC module will read the output of the first principle calculation, calculate the formation energy of each charged state of each defect according to the formula, and automatically calculate the correction value according to the correction method set by the user. If the defect calculation fails or does not converge, the formation energy will not be calculated.

• outputting figures

  The DEC module can output the figure of the defect formation energy changing with the Fermi level according to the result of the formation energy, including the data in format dat: p1.dat, p2.dat, ... , and pictures with png format: p1.png, p2.png, ... (the integer represents the number of points in the chemical potential space given by the TSC module), as well as the data and picture of the transition energy level: tl.dat and tl.png . Please draw the picture according to the instructions in the document.

• distorted defect structure calculation

  The DEC module can automatically generate the distorted defect structure based on the results of the original structure. It is worth noting that specifying the defect that will produce the distorted structure must be after the calculation of the original defect structure is completed, and running the DEC module again. For the defects that have not completed the initial structure calculation, the distorted structure will not be generated.

All details of the above generated files and calculation process will be written into 3dec.out . The information on status, results and errors can be queried in the file. After this module runs successfully, it will be noted the complete information at the end of the 3dec.out , which will be detected when the subsequent DDC module starts.

2.4 DDC Module

Defect concentration is an important parameter for the performance optimization and simulation design of semiconductor materials and devices. Experimentally, the growth conditions of the material can be controlled to inhibit the formation of harmful defects or promote favorable defects forms, to achieve the purpose of optimizing the performance of devices. In terms of calculation, a further prediction of the relationship of the defect concentration verse the growth conditions, such as chemical potential, on the basis of the defect formation energy and transition energy level results, can provide a more clear and quantitative reference for the regulation of semiconductor performance.

The work of DDC module is to read the results of formation energy and transition energy level under various chemical potentials calculated by TSC and DEC modules before, and self consistently solve the Fermi level, defect concentration, and carrier concentration of the system at a certain growth and working temperature according to the charge neutrality equation.

For a defect α in the charge state q, the concentration under equilibrium conditions can be represented as,

(5)\[N(α,q)=N_{sites}g_qexp[-ΔE_f/k_BT]\]

Where \(N_{sites}\) is the number of sites that defects can be incorporated per volume, \(g_q\) is the degeneracy factor that equals to the number of possible electron configurations for different charge states, \(ΔE_f\) is the defect formation energy. All the ionized defects in the charge state q≠0 produce carriers. The positively charged donor defects with q>0 produce electron carriers, and their summed charge is \(\sum_{\alpha,q>0} [q*n(\alpha,q)]\) ; while the negatively charged acceptor defects with q<0 produce hole carriers, and their summed charge is \(\sum_{\alpha,q<0} [(-q)*n(\alpha,q)]\) . 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,

(6)\[n_0+\sum_{\alpha,q<0} [(-q)*n(\alpha,q)]=p_0+\sum_{\alpha,q>0} [q*n(\alpha,q)]\]

where \(\sum_{\alpha,q<0} [(-q)*n(\alpha,q)]\) and \(\sum_{\alpha,q>0} [q*n(\alpha,q)]\) are the summed charges of negatively charged defects and positively charged defects, weighted by the charge q. \(n_0\) and \(p_0\) are free carrier densities, which can be defined as,

(7)\[n_0 = \int_{\varepsilon_c}^{+\infty}g(E)f(E)dE\]

(8)\[p_0 = \int_{-\infty}^{\varepsilon_v}g(E)(1-f(E))dE\]

where \(g(E)\) is the density of states (DOS), and \(f(E)\) is the Fermi-Dirac distribution function.

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 it can be assumed that the atomic structure will not change and new defects will not form in the subsequent rapid annealing process. But the densities of different charge states will redistribute during the rapid annealing. The DDC module is developed in accordance with such fabrication process, i.e. self-consistently solving the charge neutrality equation at a high growth temperature firstly, and then the Fermi level and the densities of each defect in different charge states can be obtained at the high temperature. Afterward, self-consistently solves the charge neutrality equation again at the lower working temperature, but now the formula for the concentrations calculation of defects in different charge states should be modified. In the second self-consistently solves step, the density summation for each defect over all charge states is fixed, and the density of each charge state will undergo a redistribution according to the Fermi-Dirac distribution. Then a new Fermi level can be obtained at the working temperature, and the redistributed defect densities and carrier densities can be calculated.

DDC module can calculate the Fermi level, carrier densities, and concentrations of differently charged defects under various chemical potentials, growth and working temperature based on the results and output files of TSC and DEC. It mainly includes the following steps:

• summary of formation energies

  Based on the defects calculation in DEC, DDC will automatically search the data, such as formation energy and transition energy level, that output under each defect directory, then summarize and output them in the file DefectParams.txt .

• self-consistent calculation of Fermi level

  DDC self-consistently solves the charge-neutrality equation to determine the Fermi level at growth and working temperature respectively, and the results will be written into 4ddc.out and Fermi.dat .

• calculation of defect and carrier concentration

  Based on the Fermi level, the corresponding carrier densities and differently charged defect concentration can be calculated, and output the data file Carrier.dat and Defect_charge.dat , as well as the png image density.png .

All details of the above generated files and calculation process will be written into 4ddc.out file. Relevant status, results, and error information can be queried in this file.

2.5 CDC Module

The CDC can calculate the following properties: (1) radiative carrier capture coefficient of defects, (2) lineshape of photoluminescence spectra, and (3) phonon-assisted nonradiative carrier capture coefficient (cross section).

1. radiative carrier capture coefficient

Assuming that defect \(A\) has different charge states q=0 and q=+1, and the \((0/+)\) transition energy level located in the band gap, so the defect \(A^+\) with \(+1\) charge state can capture electrons from conduction band minimum (CBM) and convert to \(A^0\) . If this process is the radiative transition (releasing photons), its capture coefficient can be expressed as:

(9)\[C_n = fη_{spin}V_{supercell}\frac{n_{r}e^{2}}{3m^{2}ε_{0}πc^{3}ħ^{2}}|\langle\psi_i|p|\psi_f\rangle|^{2}E_{opt}\]

where \(f\) is the Sommerfeld factor used to express the Coulomb interaction between the charged defect before capture and the trapped carrier. As it is Coulomb attraction for \(A^+\) to capture an electron, the order of magnitude of \(f\) is about 5 to 10. Conversely, when \(A^0\) captures a hole becomes \(A^+\) , \(f\) will be 1 due to the neutral defect neither attracting nor repelling positive and negative charges. \(η_{spin}\) is the selection rule of spin, \(V_{supercell}\) is the volume of the supercell, \(n_{r}\) represents the refractive index, \(m\) is the mass of free electrons, and \(\langle\psi_i|p|\psi_f\rangle\) is the momentum matrix element, which needs to modify the VASP source program and recompile to output. \(E_{opt}\) is the optical transition level of defects, which can be obtained by subtracting the lattice relaxation energy from the thermodynamic transition level of defects.

CDC module can call VASP automatically in batches to do first principle calculation based on the directory and output of DEC module (such as structure, transition energy level, etc.) to obtain the radiative capture coefficient in \(cm^{3}s^{-1}\) . Based on radiative capture coefficient \(C_n\) , it is necessary to multiply the coefficient with the real defect concentration \(N_D\) (refer to the output of DDC module) to obtain the capture rate \(r\) , and then take the reciprocal to obtain the lifetime of radiative transition:

(10)\[τ = \frac{1}{r} = \frac{1}{C_{n}N_{D}}\]

2. lineshape of photoluminescence spectra

Steady-state photoluminescence spectroscopy is one of the important experimental measurements to characterize defects. Previous theoretical studies often only compared the calculated transition energy level with the PL peak to speculate the possible defect configuration, while the CDC module of DASP can simulate the shape of PL spectrum, so as to provide physical quantities such as the peak position, FWHM, zero phonon line, and Huang-Rhys factor of PL spectrum for comparison with experiments. The simulation of PL spectrum is mainly dependent on the calculation of “spectral function”, that is, the intensity is written as a function of the radiant photon energy \(ħω\) :

(11)\[I(ħω) = Nω^{3}\sum_{m}p_{m}\sum_{n}|\langleχ_{im}|χ_{fn}\rangle|^{2}δ(E_{ZPL}+E_{im}-E_{fn}-ħω)\]

Where N is the normalization factor, \(χ_{im}\) and \(χ_{fn}\) are the vibrational wavefunctions of the initial and final states, and \(E_{im}\) and \(E_{fn}\) are the corresponding eigenvalues. \(E_{ZPL}\) is the zero phonon line of the transition process. In order to evaluate the accurate vibrational wavefunction and the eigenvalues, one-dimensional configuration coordinate diagram is used in CDC to plot the potential energy surfaces, and the steady-state Schrodinger equation can be solved to obtain the numerical solutions of the wavefunction and eigenvalues.

CDC module can call VASP automatically in batches to do first principle calculation based on the directory and output of DEC module (such as structure, transition energy level, etc.) to calculate the configuration diagram, wavefunction, eigenvalues, and numerically solve the above overlapping integral, and finally output the intensity with the radiation photon energy \(ħω\) , that is PL spectrum.