ADAQ: Automatic workflows for magneto-optical properties of point defects in semiconductors

Automatic Defect Analysis and Qualification (ADAQ) is a collection of automatic workflows developed for high-throughput simulations of magneto-optical properties of point defect in semiconductors. These workflows handle the vast number of defects by automating the processes to relax the unit cell of the host material, construct supercells, create point defect clusters, and execute calculations in both the electronic ground and excited states. The main outputs are the magneto-optical properties which include zero-phonon lines, zero-field splitting, and hyperfine coupling parameters. In addition, the formation energies are calculated. We demonstrate the capability of ADAQ by performing a complete characterization of the silicon vacancy in silicon carbide in the polytype 4H (4H-SiC).


I. INTRODUCTION
Point defects in wide-bandgap semiconductors have spanned a wide range of applications, including but not limited to qubit realizations [1][2][3][4][5], biosensors [6][7][8], accurate chemical sensors [9], nanoscale electric field and strain sensors [10], and nano thermometers [11]. Most of these applications have been realized with the NV center in diamond [12][13][14] -a point defect cluster in diamond consisting of a carbon vacancy and a nitrogen substitution. Recently, other point defects in diamond (such as the silicon vacancy cluster and the related group 14 vacancy clusters [15,16]), as well as point defects in other materials (such as divacancy and silicon vacancy in silicon carbide (SiC) [17][18][19] and boron vacancy in boron nitride (BN) [20,21]), have attracted remarkable interest.
Other not yet discovered point defects in various semiconductor hosts may have even more attractive properties for existing as well as novel applications.
Due to the vast number of possible point defects, to discover novel potentially interesting candidates is a challenging inverse design problem [22]. The latter means selecting the desired properties and letting the structure and materials vary. A potential starting point is to use high-throughput calculations and collect the results in a database. Previous high throughput works in this direction focused on the energetics of point defects [23,24]. Furthermore, these high-throughput workflows handle only single defects, whereas defect clusters, such as pair defects, are among the most studied defects for quantum applications.
To find and identify point defects is a labor-intensive process. Point defect identification in semiconductors may be achieved by comparing calculated magneto-optical properties with experimental values [17]. The considered magneto-optical properties may include zero phonon line (ZPL), hyperfine coupling parameters, and zero field splitting (ZFS). For the ZPL, both the energy, polarization, and intensity of the line provide information about the workings of the point defect.
A point defect can exist in different configurations depending on the host material. For reliable identifications, one should look at different configurations and different charge and spin states of each defect. A suitable automatic workflow can take all this into account and fully characterize a point defect. In addition, the formation energy can be calculated to characterize the stability of different point defects. All these data can assist the researcher to better understand the point defect and evaluate if it is useful for a given application.
In this paper, we present the full characterization workflow, which is the primary component of ADAQ, that allows for high-throughput calculations of magneto-optical properties of point de-fects and their clusters of arbitrary size for finding and identifying potentially interesting systems. Furthermore, ADAQ also contains a simplified workflow for quick estimates of the ZPL and energy of point defects intended for high-throughput screening of defects to identify candidates for which the full characterization workflow can be deployed. The present work goes beyond similar prior efforts in its focus on magneto-optical properties and point defect clusters, both of which have not been considered earlier in the context of automated workflows. We demonstrate how these properties are obtained with ADAQ considering a well-known defect: the silicon vacancy in 4H-SiC [18].
The outline of the paper is the following. Section II introduces the properties calculated in ADAQ and the theory behind the calculations. Section III describes the software used to set up and execute the calculations, as well as the settings for the first principle software used. In addition, an overview of ADAQ and its default input parameters, such as the size of the supercell, are presented, as well as details for the unit cell and host workflows. The full characterization workflow is described in Section IV. Here, overviews are presented for both the ground and excited state workflows as well as details about each step in the workflows. The workflow results are stored in a database, and Section V outlines which properties are stored and how. Section VI shows the information that ADAQ produces for a point defect, illustrated by the silicon vacancy in 4H-SiC.
The strengths and limitations of the full characterization workflow are examined in Section VII and conclusions are presented in Section VIII. Appendix A shows the algorithm for defect generation which constructs point defect clusters up to an arbitrary size. Appendix B outlines the screening workflow that produces quick estimates of the ZPL and energy of the point defect. There is also a list of acronyms at the end of the paper.

II. THEORY
The subsections below present the theoretical background on first-principles calculations of the point defect properties built-in into ADAQ.

A. Photoluminescence
Point defects in semiconductors may introduce states in the band gap. If an optical transition between these states is allowed and the non-radiative decay rate is low, a ZPL will appear in the photoluminescence spectrum. One of the main tasks of ADAQ is to calculate whether a ZPL exists for a given defect and predict its energy. Figure 1 a)  parabola represents the ground (excited) state energy surface. The displacement between different ion minima in the ground and excited state is ∆Q.

ZPL
The sharp peaks in the photoluminescence spectra, like the one in Figure 1a), are called ZPL, which arise from the direct transition between the defect states in the band gap. They can be used to identify the point defect present in the material. As seen in Figure 1b), the excitation cycle starts with the absorption of a photon that excites the system from the ground state to the excited state, followed by a relaxation of the ions to the excited state equilibrium configuration. In the excited state, there are two possibilities for radiative decay back to the ground state. The state can either decay straight back to the ground state minimum with no phonon contribution, which produces the ZPL, or decay to a higher energy atomic configuration and relax to the lowest energy configuration through phonon emission. Given the energy for the ground and excited states, the ZPL energy is defined as: where E e,min (E g,min ) is the excited (ground) state in its corresponding minima. Here, the two states are assumed to have similar phonon spectra, and the zero-point energy of the phonon modes cancels out between the ground and excited states. As mentioned before, the ZPL exists as long as the non-radiative decay is small. The non-radiative effects are not taken into account in the calculations. The rate of non-radiatively relaxation increases exponentially as the energy difference between excited and ground states decreases. Hence, predicted ZPLs below 0.6 eV are considered to be uncertain.
If a ZPL exists, additional factors could be used to determine if the defect is a promising candidate for a particular application, such as a single photon emitter. Three additional properties that can be used for the characterization include the transition dipole moment (TDM) that describes the polarization and intensity of the ZPL [25], the ion relaxation between the ground and excited state which indicates the ratio between the ZPL and the phonon sideband emission, and the spontaneous macroscopic polarization which tells if the ZPL is stable against electric fields.

Intensity of ZPL
The TDM is calculated between the single particle orbitals of the transition and is defined as: Here, qr is the dipole operator between the initial ψ i and final state ψ f , which can be rewritten into reciprocal space with the momentum operator p instead and a prefactor with the eigenvalue difference between the final and initial states, f,k − i,k . The constants used in this formula are the Planck constant and the mass of the electron m. The µ is calculated for all optical transitions in the excitation cycle and from it, the optical polarization of absorption, ZPL, and emission are extracted as well as the intensity |μ| 2 . The optical lifetime τ is calculated from the intensity |μ| 2 [25] with the following equation where the 0 and c are the vacuum permittivity and the speed of light, respectively, ν is the transition frequency for the specific transition in question (absorption, ZPL, or emission) and n is the refractive index, which for 4H-SiC is 2.6473. This method has been applied to the divacancy in 4H-SiC, where the inclusion of the ion relaxation and the corresponding electronic change of the excited state is crucial to include when calculating µ for the ZPL [25].

Ion relaxation
An additional desired property for a bright single photon emitter is a large ratio between the photon count of the ZPL and the phonon emission sideband. This ratio can be determined from the Huang-Rhys factor, which calculates the coupling between the electronic and vibronic states [26].
This factor is expensive to calculate. Hence a 1D model with two measures of ion relaxation between the ground and excited state has been introduced [27] and tested [28]. These two measures are the squared displacement of the ions ∆R, and a parameter in which each displacement is scaled with the atom weight ∆Q, which is shown in Figure 1b). ∆R is calculated as: whereR ea andR ga are the ion positions in the excited and ground state, respectively. The summation runs over all ions in the supercell. ∆Q is calculated as: where m a is the atomic mass.

Spontaneous macroscopic polarization
For an ideal single photon emitter, the ZPL should be stable against electric fields, i.e., show no spectral diffusion [29]. This can be achieved if the spontaneous macroscopic polarization, that couples the defect to external electric field fluctuations, is small. The spontaneous macroscopic polarization is the difference between the ground and excited state macroscopic polarization which is calculated from the Berry phase according to the modern theory of polarization [30][31][32]. However, one may need to do additional calculations to ensure that no wrap-around error is present.

Defect ionization and affinity energy
Depending on the positions of the defect states in the band gap with respect to the conduction and valence band edges, it is possible that an optical excitation can change the charge state of the point defect. This can happen either if an electron moves from a defect state to the conduction band edge (bound-to-free transition), i.e. defect ionization energy, or if an electron moves from the valence band edge to a defect state (free-to-bound transition), i.e. defect affinity energy. If any of these transitions require lower energy than the ZPL, the point defect will change the charge state before emitting a ZPL.

B. Hyperfine coupling parameters and ZFS
For point defects with spin (unpaired electrons), additional interactions occur. The unpaired electrons, which have magnetic moments, interact with magnetic fields and other magnetic moments. If some ions in the material have spin (a non-zero nuclear g-factor, which are listed at easyspin website [33] for all isotopes), the magnetic moments of the ion and electron can couple and give rise to hyperfine interaction. The hyperfine interaction [34] is defined as: HereŜ is the electron spin operator,Î is the nuclear spin operator, and A is the hyperfine tensor.
We define A xx , A yy , and A zz as the diagonal components of the A tensor and A z , the projection on the z-axis, as the hyperfine splitting. The tensor only exists if both the spin of the defect and ions are non-zero. In practice, this tensor is approximated by calculating the Fermi-contact interaction, which depends on the spin density at the center of the nucleus, and the dipole-dipole interaction.
If the point defect has more than one unpaired electron (at least spin-1), the unpaired electrons interact with each other and separate states with different absolute spin quantum numbers without a magnetic field present, producing a ZFS. The ZFS comes from the D-tensor. The D-tensor describes the interaction for the total effective spin, which can be approximated with spin-spin dipole interaction in semiconductors of light elements [35] and is defined as: HereŜ is the electron spin operator and D is the D-tensor which is traceless and symmetric. Usually, the diagonal elements of the D-tensor are ordered in increasing order, where the z component is the largest. From the D-tensor, the ZFS is calculated as 3 2 D z for spin-1 defects. A detailed description of the methodology for the D-tensor calculation can be found in Ref. 35.

C. Formation energy and charge-state transition levels
The formation energy is the energy needed to create the defect. It can be useful when comparing different defects to see which one is the most stable. The charge-state transition levels tell us where the Fermi energy needs to be in the material in order to keep the defect in a given charge state. The formation energy for charged defects is defined as [23,36]: The formation energy ∆H D,q (E f , µ) depends on the Fermi energy E f and the chemical potential µ, both is discussed further below. E D,q and E H are the total energies for the charged defected supercell and the perfect host supercell, respectively. The variable n i keeps track of the atoms of chemical element i added (n i < 0) or removed (n i > 0), q is the charge of the defect supercell and E corr is the correction term needed to minimize finite-size effects.
From the formation energy, charge-state transition levels can be found [36]. These levels show where a transition from one charge state to another occurs and are defined as: where E f is set to zero.

Chemical potential
In Eq. (8), the chemical potentials for the elements involved are needed. The chemical potentials used in this paper are calculated with the same settings as in the workflow for all elements.
This calculation gives the total energy (calculated at zero temperature) for each element in the periodic table. The structures for the elements are taken from Ref. 37, which has a list of all the elements in their periodic ground state structure at zero temperature. Therefore, all elements are treated equally, and thus the reference states for e.g., N2 and O2 are not isolated molecules, but the 0K elemental structures (which for these are dimers in a periodic structure). Furthermore, we use the energies as calculated by VASP without molecular energy correction, which is often used in similar phase diagrams (see e.g., Ref. 38). The total energy will give an upper limit to the chemical potential. For SiC, the silicon and carbon ground state will give an upper limit to the chemical potential, denoted a rich phase (µ rich ). In addition, the following relation holds: µ SiC = µ Si + µ C .
Using the upper limit for one chemical element, the lower limit for the other chemical element can be obtained from this relation. Here, one assumes that there is only one stable phase (SiC) between the elements (Si, C) on the convex hull. If there would be other stable phases on the convex hull, like SiC 2 , one should consider these like in Ref. 24.

Potential alignment
The formation energy is plotted against the Fermi energy, and most often, the valence band maximum (VBM) is set to zero. This convention means that the endpoint of the Fermi energy is the bandgap energy, i.e., conduction band minimum (CBm). However, when comparing the host supercell to a charged defect supercell, the VBM and CBm may be shifted. This shift can be accounted for by comparing the average potential far away from the defect. In this work, the supercells are large enough that potential alignment is negligible. Thus the VBM and CBm are taken straight from the host supercell.

Charge correction
When comparing the energy of periodic cells of different charge states, one has to account for the self-interaction energy contribution of the extra charge. Different charge correction schemes have been suggested: Makov-Payne (MP) [39], Lany-Zunger (LZ) [40], and Freysoldt-Neugebauer-Van de Walle (FNV) [41]. FNV is the most accurate but non-trivial to use for defect clusters and might run into computational difficulties [42], making this correction challenging to use in highthroughput frameworks. Hence, we choose the LZ correction scheme that gives the same correction for all defects. The LZ correction E corr is defined as where f is a proportionality factor linking the third-order correction to the first-order, q is the charge, α M is the Madelung constant, = 4π 0 r where r is the dielectric constant, and L is the length of the supercell (L = V 1/3 ).

III. METHODOLOGY
The High-Throughput Toolkit (httk) [43,44] is used for automatic control of the ADAQ calculations. It is a framework that handles transferring calculations between a local computer and a supercomputer and executing the runs by running taskmanagers. The taskmanager includes checkers that monitor the runs, ensure that the runs converge as intended, and cancel any run that breaks the predefined rules. Any software can be controlled by httk, but for ADAQ, all workflows are presently implemented only for the Vienna Ab initio Simulation Package (VASP) [45,46].
VASP runs density functional theory (DFT) [47][48][49] calculations with projector augmented wave (PAW) [50,51] method. For the excited state calculation, the constrained occupation DFT method [52,53] is used. This method is also known as ∆SCF [52,54,55] when using the total energy difference between the ground and excited state. Formally, one should use the DFT formalism of Görling (generalized adiabatic connection [56]) to include stationary densities and handle the excited states [57]. However, this formalism introduces an orbital-dependent exchange-correlation functional. In practical calculations, this functional is approximated by the exchange-correlation functional approximation from the standard DFT formalism. The exchange-correlation effects are described by the semi-local functional of Perdew, Burke, and Ernzerhof (PBE) [58]. For calculating ZPL, we use the method described in detail and tested for the divacancy point defect in 4H-SiC in Ref. 17. In this previous paper, we concluded that using the PBE functional for a 576 atom supercell with 2 × 2 × 2 k-point set provides a good compromise between the accuracy and efficiency, which is suitable for high-throughput calculations. Due to the use of the PBE functional, we have found a systematic underestimation of about 0.2 eV for the calculated ZPL compared with experiment for 4H-SiC. The same offset can be seen for the neutral divacancy and charged silicon vacancy in 6H-SiC [19]. The hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE06) [59,60] corrects this error but is not suitable for high-throughput calculations at present due to its high computational cost [17]. To compare the ADAQ results with experimental data, one should add this systematic shift of 0.2 eV to the estimated ZPL. If higher accuracy is needed, HSE calculations can be run on top of the PBE results which produces accurate results as tested in Ref. 17. For even greater accuracy, fully self-constituent HSE or other higher-order methods, such as GW, are possible. Workflows using these higher-order methods are presently not part of ADAQ, however, such an extension would be relatively straightforward.
When running high-throughput calculations for a wide range of different elements, convergence settings must be chosen to match the requirements of the most demanding chemical elements. To ensure sufficient convergence for all elements, the plane wave energy cutoff is set to 600 eV and kinetic energy cutoff to 900 eV. These values are in the middle of the range recommended for the more demanding PAW pseudopotential and are slightly larger than those used in the Materials Project (520 eV) [61,62]. Unless specified otherwise, we apply the following settings: • The Fast Fourier transform (FFT) grid is set to twice the largest wave vector.
• When only Γ-point is used, Fermi smearing with a width of 1 meV, otherwise the tetrahedron method with Blöchl corrections [64] is used.
• For defect calculations, symmetry is not used.
• Non-spherical contributions of the PAW spheres are included.
• Projection is done in real space with automatic optimization.
• The ions are updated using the quasi-Newton method.
Using httk and VASP, the workflows are constructed for ADAQ. Figure 2 shows an overview of all the workflows included in ADAQ with their corresponding inputs. Computational details for the two smaller workflows, the unit cell and host workflows, are presented below. Details of the more extensive workflows, ground and excited state workflows, are presented later in this paper.
Here are the details for the two smaller workflows in ADAQ, the unit cell and host. First, the unit cell of the chosen host material is relaxed. The unit cell workflow carries out volume relaxations with PBE functional as default. To ensure accurate volume, these calculations are executed with a k-point set 10 × 10 × 10, Gaussian smearing with a width of 1 meV is used, the plane wave energy cutoff is set to 700 eV, the plane wave kinetic energy cutoff is 1400 eV, the electronic convergence criterion is 10 −6 eV, the ion convergence criterion is 5 · 10 −5 eV with respect to energy, and the projection is done in reciprocal space. When the volume of the unit cell change, the plane wave basis set might not be as accurate as the starting settings. To handle this, the unit cell workflow rerelaxes the structure until the energy difference between iterations is smaller than 5 · 10 −4 eV.
After the unit cell has been optimized, the supercell is created. The supercell size is set to be approximately 20Å in every direction as default to ensure a low defect-defect interaction. To preserve the symmetry of the crystal, the unit cell is copied until the size criterion is met. To get the energy of the host supercell, it is processed in the host workflow that is similar to the ground state workflow (Sec. IV A) but only runs the first four steps but with symmetry turned on. Then the point defect supercells, which will be denoted defect cells from now on, are generated as described in Appendix A. Now, the defect cells can either be screened -as described in Appendix B -or run through the full characterization workflow directly.

IV. FULL CHARACTERIZATION WORKFLOW
After the host material has been selected, the unit cell relaxed, and the defect cells have been generated, the full characterization workflow can start. Figure 3 shows an overview of the different steps involved to fully characterize any point defect cluster. Hereinafter, we refer to these steps as neutral, charge, and spin step, which include neutral and charged ground state calculations and alternative spin calculations, respectively. These steps are processed with the ground state workflow, see Sec. IV A. First, the neutral defect is processed. After this step has finished and there are defect states present in the band gap, the charge step follows. The charge step runs the ground state workflow again but now with charged supercells. The default settings remove and add up to two electrons. These charge runs, a maximum of four, are processed in parallel. When the charge step is finished, alternative spin states are processed with the ground state workflow.
For this step, the ground state electronic occupation is changed (for example from spin-3/2 to spin-

Start
First, the data files with structure and computational parameters are copied into the running directory, and httk selects the PAW pseudopotentials. For charge and spin calculations, the number of electrons and fixed occupation is added to the input, respectively.

Prerelax
This first step is a fast and coarse ion relaxation to ensure the ions are not too close. Here, the convergence settings are reduced: electronic convergence criterion is 10 −4 eV; ion convergence criterion is 0.3 eV with respect to forces and includes 20 ionic relaxation steps; spin polarization is turned off; the FFT grid is set to 3/2 of the largest wave vector; and only the Γ-point is used for the integration over the Brillouin zone. When the calculation is finished, the wave function is saved and used as an input for the next step. Since this step only uses the Γ point, it could be run with the gamma-only version of the VASP software. However, the wave function from this version cannot be loaded into the standard version of VASP, and our tests show that the overall time is lower if the wave function is propagated from the prerelax to relax step. This step is only carried out for the neutral defect.

Relax
Here, the charge density and wave function from the prerelax step is loaded. The workflow starts here for charge and spin steps and loads the wave function from the neutral prerelax step. In this step, the convergence parameters are increased: ion convergence criterion is 5 · 10 −4 eV with respect to energy and 30 ionic relaxation steps; spin polarization is turned on; the FFT grid is set to twice of the largest wave vector; and 2 × 2 × 2 k-point grid is used.

Final relax
In this final ion relaxation step, the electronic convergence criterion is increased to 10 −6 eV, and the ion convergence criterion is changed 5 · 10 −5 eV, still with respect to forces and up to 30 ionic relaxation steps are included. If the calculations fail to relax the ions during any of these three relaxation steps, the httk checkers try to find and correct the error to converge the calculation in this final step. If this is not possible, it aborts and proceeds to the cleanup step. When the ion relaxation is completed, the final wave function is saved and analyzed in the next step.

Analyze
First, if fractional occupation occurs, the nearest non-fractional occupation is forced, and the final relax step is repeated. Next, one needs to identify if defect states are present in the band gap. This identification uses the inverse participation ratio (IPR) [65,66], which is a measure of localization. In practice, the IPR is calculated with a python library PyVaspwfc [67] for the 30 bands closest to the Fermi level in each spin channel and averaged over all k-points. If a band has an IPR larger than the threshold (10 −4 ), it is considered a defect state. Sometimes, stray band appears in the valence or conduction band which goes above this threshold suggesting it may be a defect state. The identified defect bands should be continuous, and thus any outliers are removed, these are saved in a separate file since they can provide further information about the point defect but cannot be handled in the workflow. If no local states are found, the workflow ends.
Once the defect states have been identified, the input for the rest of the workflow can be con- where the first term is the binomial coefficient, the second term excludes the ground state, and the

Partial densities
The first post-processing step calculates and saves the partial densities for the local states as well as the top and bottom of the valence and conduction band, respectively.

Hyperfine coupling parameters
If the spin of the defect and the nuclear spin of the chemical elements are non-zero, hyperfine tensors are calculated. The nuclear spin depends on the nuclear g-factor, which in turn depends on the isotope. Only certain ions have a non-zero nuclear g-factor, for example, 13 C has a non-zero value while 12 C and 14 C do not. The g-factors are extracted from the easyspin website [33]. The algorithm calculates the hyperfine interaction for all possible paramagnetic isotopes. An intrinsic defect in SiC has only 13 C and 29 Si with non-zero values, hence only one hyperfine interaction exists. But if one would dope with B and N, which have non-zero values for both isotopes, then four hyperfine interactions would exist. One calculation of the hyperfine coupling parameters is carried out in VASP with g-factors set to unity, and the results are multiplied with the different g-factor to get all the hyperfine interactions.

Zero-field splitting
If the spin of the defect is larger than one-half, the D-tensor is calculated as described in Ref. 35.

Macroscopic polarization
The final post-processing step calculates the macroscopic polarization of the defect cell using Berry phase calculation.  iterations or if the final energy difference between this state and final relax ground state is lower than the experimental limit (0.4 eV, see Sec. III), the run is stopped, and the workflow proceeds to the cleanup step.

Excitation relax
In this step, the same settings are used as in the relax step in the ground state workflow, see Sec. IV A 3, but with excited state occupation.

Excitation final relax
In this step, the same settings are used as in the final relax step in the ground state workflow, see Sec. IV A 4, but with excited state occupation.

Macroscopic polarization
If the excitation is local, the macroscopic polarization of the excited state is calculated the same way as in Sec. IV A 9.

Emission
When the final excited state geometry is found, the occupation is set to the ground state occupation. This step is the counterpart to the absorption step; here, the excited state geometry is used with the ground state occupation.

Post-processing
Between each step in the excited state workflow, the following post-processing steps take place.
The partial charge densities are calculated the same way as in the ground state workflow. Here, the TDM is also calculated in two ways. First, between each step in the excited state workflow, for example, from absorption to excitation relax step. Second, between the excited state step and its corresponding ground state step, for example, from excitation relax step in excited state workflow to relax step in ground state workflow. To calculate the TDM, the WAVECARs are post-processed using the PyVaspwfc python library [67], which we have modified to handle two WAVECARs and calculate the TDM between the defect orbitals.

B. Hyperfine coupling parameters and ZFS
For the hyperfine interaction, the three different tables related to hyperfine coupling parameters, as described in the VASP manual [68], are extracted from the run and post-processed. Since the calculations are run with the g-factors set to unity, each atom is multiplied with the corresponding g-factor to produce all hyperfine interactions (Sec. IV A 7). For each atom, the bipolar hyperfine coupling parameter matrix is extracted, and the Fermi contact coupling parameter (A 1c and A tot ) diagonal matrix is constructed. The eigenvalues and eigenvectors of these added matrices are calculated as well as the eigenvectors' angles to the z-axis. Since symmetrical equivalent atoms have the same hyperfine coupling parameters, the multiplicity of each hyperfine value is also marked. The hyperfine coupling parameters are only saved if the hyperfine splitting (A z ) is larger than 3 MHz.
The D-tensor is extracted from the output, both the calculated and diagonalized version with corresponding eigenvectors.

C. Bands with IPR
The 30 bands closest to the highest occupied band in both spin channels are saved to the database. This includes the local bands, if any exist, and a few bands in the VB and CB, respectively. For easy identification of local bands, the IPR value (discussed in detail in Sec. IV A 5) is also stored for each band and both spin channels.

D. Formation energy and charge-state transition levels
Formation energy for all charge and spin states as well as charge-state transition levels between all different charge and spin states are calculated and stored.

VI. COMPUTATIONAL RESULTS
To demonstrate the results produced by ADAQ, the silicon vacancy (V Si ) in 4H-SiC is used as an example. SiC is a technologically mature material where it is possible to combine quantum and classical applications in the same device [69]. However, it is a material with many different polytypes and multiple unidentified defects. The silicon vacancy has two configurations in 4H-SiC, denoted h and k. Given that these configurations have been identified earlier [69], both these configurations are processed directly through the full characterization workflow. If these had been unknown defects, the best approach would be to generate an array of different point defect clusters and run these defects through the screening workflow (Appendix B), a scaled-down version of the full characterization workflow. After this step, the most likely candidates are processed through the full characterization workflow to identify the unknown defects. We present data for the silicon vacancy from both the full characterization and screening workflow, starting with the former.

A. Full characterization workflow results
First, Figure 6 shows the formation energy with charge transition levels for both configurations of the V Si , with the chemical potential for the silicon: µ rich Si = −5.42 eV. Using the charge correction formula in Eq. (10) with the values (1 + f ) = 0.65, α M = 2.8373 (Madelung constant for simple cubic [70] since the 576 supercell is close to cubic), r = 9.6 (taken from experiment [71]), and L = V 1/3 = 18.21Å for the 576 atom supercell gives a correction of 0.076 eV for q = ±1 and 0.304 eV for q = ±2.   For the negatively charged spin-3/2 states, Table I   Additional data about hyperfine coupling parameters and ZFS for the negatively charged spin-3/2 states are presented in Table II and Table III. For both configurations of the V Si , the first row is the largest hyperfine splitting related to the carbon atom above the silicon vacancy, and the second row is the three carbons below the silicon vacancy. Both the hyperfine coupling parameters and ZFS show a small but perceivable difference between the two configurations. These results and trends are comparable with HSE06 hyperfine and ZFS data [18].  As a demonstration of the screening workflow described in Appendix B, the silicon vacancy is also processed through this workflow and the results are presented in Table IV. Eigendifference stands for eigenvalue difference.

VII. DISCUSSION
In this section, we discuss the strengths and limitations of the screening and full characterization workflows, starting with the latter.

A. Full characterization workflow comments
We start with the formation energy in Figure 6. The negative charge spin-3/2 state is the only state with a ZPL for both configurations in the reported experimental range, k configuration is 1.166 eV and h configuration is 1.240 eV. Note that these ZPL calculated with the PBE functional are systematically underestimated by 0.2 eV compared with experiment but leaves the order unaffected [17]. Adding this 0.2 eV correction to the calculated results brings the ZPL almost on top of the experimental values, and the given order (k is lower than h) agrees with previous identification [18]. When comparing the corrected ZPL results from ADAQ with experiment, one gets an accurate picture of which ZPL belongs to which charge and spin state of the defect. This is what ADAQ is designed to do.
Other useful data, which can help experimentalists, are the other lines presented in Figure 7.
The tiny red line close to 1.4 eV in Figure 7 k) shows valence band edge to defect state transition, corresponding to the defect affinity energy. It gives the limit of the laser energy used to excite the electron. If an excitation energy larger than 1.6 eV (1.4+0.2 eV) is used, the defect will change to a double negative charge state as discussed in Sec. II A 5, thus removing the visible ZPL from the spectra. The 0.2 eV correction may be off for the free-to-bound and bound-to-free transitions since is estimated for bound-to-bound transitions.
Why can we neglect the charge correction for the ZPL calculations? For the ZPL, there is a small multipole change between the ground and excited state. This multipole change has a minuscule effect on the ZPL energy. For the h configuration using the FNV correction in both the ground and excited state cause only a 5 meV difference in the ZPL. This energy change is calculated between two local states, whereas the maximum multipole change comes from a local to delocalized transition. An excitation from a local state to the conduction band minimum gives a slightly larger correction of 25 meV. These corrections are still smaller than other uncertainties; hence the charge correction for the ZPL can be neglected.
An analysis of the band gap states for the negatively charged spin-3/2 state plotted in Figure 8 can provide additional information. In this plot, the defect band just above the VBM is sufficiently is most significant at the Γ-point. The IPR is not calculated for the excited state, mainly because it is not needed in the workflow. However, manual testing on the divacancy defect in 4H-SiC does not show any discernible difference in localization between the ground and excited state: the excited state IPR is almost identical compared to the ground state IPR. The ZPL of interest is the transition between this band and the nearest higher energy band. If the band closest to the VBM had not been identified as a defect state, the ZPL would have been a free-to-bound transition. One should be aware that this might happen but as a precaution, the full characterization workflow already calculates the excitations to and from the valence and conduction band edges.
The transition energies with polarization, intensity, and lifetime are displayed in Table I Table I)  The ion relaxation factors ∆R and ∆Q are similar for both configurations. Since these values are small for the silicon vacancy, it suggests that the defect may be a good single photon emitter because of the spectral stability. These values are also close to the values for the NV-center in diamond [27].
The fact that the negatively charged spin-3/2 state is responsible for the ZPL observed in experiment can further be supported by the hyperfine coupling parameters and ZFS. For the hyperfine splitting, the h configuration has slightly larger splitting for both the carbon above and the three carbons below the silicon vacancy than the k configuration. This agrees with the experimental findings in Ref. 18. The theoretical results are about 10 Mhz lower than the experimental values.
The underestimation is due to the choice of the PBE functional, which typically overestimates the delocalization of the orbitals. It has been demonstrated that the HSE functional gives hyperfine splitting results about 5 Mhz higher than the experimental values [18]. The largest hyperfine splitting for silicon atoms should have a multiplicity of 12. However, since the symmetry is off and dynamic Jahn-Teller is neglected, these multiplicities split into two different degeneracies of groupings of 9 and 3 for the h configuration and three different degeneracies of groupings of 6,3, and 3 for the k configuration. Considering the ZFS, the h configuration has a lower ZFS than the k configuration, which also agrees with experiment [18]. Overall, the data produced by the full characterization workflow gives a clear picture of the negatively charged spin-3/2 state of the silicon vacancy.

B. Screening workflow comments
Let us next discuss the accuracy of the results from the screening workflow, described in Appendix B, in comparison to the full characterization workflow. Overall, the ground state properties presented in Table IV agree with the results for the full characterization workflow, shown in  If more accurate calculations are required, higher-order methods can be run by hand or, better yet, constructing additional workflows that use the data present in the database as a starting point.
For example, running the HSE functional on top of the PBE results is an efficient procedure [17].
Overall, ADAQ provides a strong starting point in the search for novel systems based on defects in semiconductors [57].

VIII. CONCLUSION
We have developed ADAQ which consists of automatic workflows for high-throughput calculations of magneto-optical properties of point defects and their complexes in semiconductors. To handle the vast number of possible defects, the following strategy for using ADAQ is proposed.
First, the screening workflow should be employed to consider a large number of defects (∼10 000).
This workflow provides the ZPL and TDM for one potentially interesting transition per defect con- This appendix shows how the defect supercells are generated. Figure 9 shows a schematic picture of the defect generation process. First, the unit cell is analyzed for symmetry, and nonequivalent atom positions are found. The interstitial locations are found by a combination of Wyckoff positions and Voronoi tessellation, similar to Ref. 23. The Voronoi locations, found using Voro++ [77][78][79], are mapped to Wyckoff positions, and the interstitial locations are then added to the unit cell in symmetric order (highest symmetry first). For all interstitial locations, the symmetric copies in the unit cell are also saved. Any new interstitial location must be farther away from the already found interstitial locations with a minimal interstitial-interstitial distance (default: 0.5Å). After all the interstitial locations have been found, the supercell is constructed by copying the unit cell in x, y, and z-direction so each of the lattice vectors meets at least the input supercell size criterion λ (default: 20Å), see sitions are available. Only atoms inside both the blue sphere and the green sphere are potential candidates. Second, the smallest defect-defect distance between a vacancy and an interstitial (1Å) is introduced to avoid that the interstitial relaxes into the vacancy position, thus creating a substitution. Figure 9 d) shows the potential interstitial positions. If the interstitial positions are too close or too far away from the vacancy, they are excluded, leaving two potential interstitial positions for a defect cluster consisting of a vacancy and an interstitial.
For these potential positions, there is a possibility that the defect might have been created before or a symmetric equivalent defect cluster exists. For example, in Figure 9  These symmetric matrices are sorted and saved during the generation. If the same matrix is found, the present defect cluster is neglected. This makes sure that only one of the two defect clusters in Figure 9 d) is created and stored in the database. Also, a unique hash index is calculated from these matrices to keep track of the defect in the database.
The defects can also be generated with extrinsic elements. Here, they can be doped with elements from H through Bi except for the noble gases (He, Ne, Ar, Kr). Above Bi in the periodic table, elements become increasingly radioactive and thus neglected as dopants. Since the noble gases do not form any compound, except Xe, they are also neglected as dopants. One should not expect that the elements with d and f electrons will produce accurate results with the present level of theory, but the defect clusters can be generated never the less. Depending on the charge state of the defect, d or f orbitals might be unoccupied, and some of them may be interesting to look at. In this appendix, a scaled-down version of the full characterization workflow is presented to screen defects. Figure 10 shows different convergence settings used for the ZPL calculations for the divacancy in 4H-SiC. The full characterization workflow runs with the PBE 2 × 2 × 2 settings. The order of the configurations is the same as in experiment, and the absolute values are systematically underestimated by about 0.2 eV [17]. Reducing the k-points to only Γ-point increases the absolute values but the different configurations can no longer be correctly identified, as shown at PBE 1 × 1 × 1. This method still requires both ground and excited state calculations to get these values. However, looking only at the difference between eigenvalues in the ground state, one can estimate the ZPL energy without calculating the excited state. This difference is marked as PBE eigendifference in the plot where the relaxation of the excited state is neglected, and the energy is surprisingly close to experimental values, at least for the divacancy. This agreement with experiment is a fortunate error cancellation that the converged PBE values differ by 0.2 eV, and reducing k-point as well as neglecting the relaxation of the excited state makes up this difference.
One cannot assume that this will be the case for all possible defects. Hence, the eigendifferences are calculated first. If a eigendifference is larger than 0.4 eV, the excited state with Γ point is calculated. This approach is good enough to determine if a defect is bright or dark and identifies a range that contains the ZPL.  Fermi smearing, and the FFT grid is set to 3/2 of the largest wave vector. The following steps have been changed from the ground state workflow, Sec. IV A.

Prerelax
Same as in the ground state workflow but now runs with gamma compiled version.

Relax
The k-point set is reduced to Γ-point only and also runs with gamma compiled version.

Analyze
The same analysis as in the ground state is performed, even though most of the output is not used. When the defect states have been found, the differences between the occupied and unoccupied eigenvalues are calculated, and the lowest difference is estimated to be the ZPL. This works as long as there are no singly occupied, almost degenerate, state present for the defect. To prevent calculating excitations between almost degenerate states, the threshold limit of 0.4 eV is used to exclude any transition below this value. The smallest transition above this limit is called the transition of interest. If no such transition exists, then the workflow proceeds to the cleanup step.

Excited relax
This step calculates the transition of interest found in the analyze step with the same settings as in the relax step above.

Cleanup
If everything has finished correctly, then the cleanup removes the WAVECARs as discussed in Sec. IV C. If the defect does not converge because of a particular charge state is unstable or the starting geometry is too far from a stable position, the job is stopped and cleaned up without analyzing it. These are marked as not converged run and still stored in the database.

Database
Similar data is saved to the database for the screening workflow as for the full characterization workflow. Additional stored data include the ZPL estimate, TDM, and partial density difference.
The ZPL estimate is taken from the lowest eigenvalue difference. The TDM and the partial density difference are calculated for the ground state for all local defect states. The ion relaxation from the input geometry is estimated to show how much the defect relaxed. If a transition of interest has been calculated, the total energy difference, TDM, and partial density difference between the excited and ground state are also saved to the database.

AVAILABILITY
The data presented in this paper is available via figshare [80]. For availability of the ADAQ workflow software, see https://httk.org/adaq/.