Spin–spin interactions in defects in solids from mixed all-electron and pseudopotential first-principles calculations

Understanding the quantum dynamics of spin defects and their coherence properties requires an accurate modeling of spin-spin interaction in solids and molecules, for example by using spin Hamiltonians with parameters obtained from first principles calculations. We present a real-space approach based on density functional theory for the calculation of spin-Hamiltonian parameters, where only selected atoms are treated at the all-electron level, while the rest of the system is described with the pseudopotential approximation. Our approach permits calculations for systems containing more than 1000 atoms, as demonstrated for defects in diamond and silicon carbide. We show that only a small number of atoms surrounding the defect needs to be treated at the all-electron level, in order to obtain an overall all-electron accuracy for hyperfine and zero-field splitting tensors. We also present results for coherence times, computed with the cluster correlation expansion method, highlighting the importance of accurate spin-Hamiltonian parameters for quantitative predictions of spin dynamics.


INTRODUCTION
Spin defects in semiconductors are promising quantum bits (qubits) for quantum information technologies including quantum computation, communication and sensing 1,2 . A prime example of spin defects is the nitrogen-vacancy (NV) center in diamond [3][4][5][6][7] , which can be optically initialized and read-out, and possesses millisecond coherence time even at room temperature. In recent years, much effort has been devoted to realizing spin defects in industrially mature host materials with properties similar or superior to those of diamond NV centers. For instance, several promising spin defects have been identified in silicon carbide, including the divacancy (VV) 8,9 , Cr [10][11][12] , and V impurities 13 . There is also a growing interest in discovering and designing spin qubits in piezo-electric materials such as aluminum nitride 14,15 , in oxides 16 , and in 2D materials 17,18 .
First-principles calculations based on density functional theory (DFT) have played an important role in the discovery and identification of spin defects, in particular in understanding their electronic and thermodynamical properties 19,20 . DFT results have been instrumental in interpreting optical and magnetic measurements, and in predicting atomistic and electronic structures of defects yet to be realized experimentally. In addition, specific spin dynamical properties may be investigated with the aid of DFT calculations, using spin Hamiltonians (SH) with parameters obtained from first principles. For systems with a single effective electron spin, e.g., a spin defect in a semiconductor, the leading terms in the SH are [21][22][23] : (1) where μ B is Bohr magneton; S is the effective electron spin; B is the external magnetic field; I N and γ N are the spin and gyromagnetic ratio of the N th nucleus; g, A, D, and P are rank-2 tensors that characterize the strength of electron Zeeman interaction, hyperfine interaction, zero-field splitting and nuclear quadrupole interaction, respectively. Nuclear spin-spin interactions and the chemical shielding effect in nuclear Zeeman interactions are neglected in Eq. (1).
The SH parameters g, A, D, and P can be determined from first principles electronic structure calculations [24][25][26][27][28][29][30][31][32][33][34][35][36] , the majority of which are based on plane-wave pseudopotential (PW-PP) approaches. DFT calculations of SH parameters using basis sets different from PW have been proposed (e.g., numerical atomic orbitals 37 , linearized augmented plane-wave 38 , linear muffin-tin orbitals 26,39 , and Gaussian orbitals 40 ), but they are often limited to smaller systems than those accessible to PW calculations. In the PW-PP method, PP are used to describe the interaction between valence and core electrons, and single-particle wavefunctions of core electrons in the solid are not explicitly evaluated. All-electron (AE) wavefunctions may be reconstructed, for example using the projected augmented wave (PAW) procedure 41 , and then used to compute the parameters of the SH. Recently, we proposed and benchmarked a realspace AE DFT framework 42 using a finite-element (FE) basis set 43 for accurate predictions of SH parameters in molecules and solids, which does not require any reconstruction of corewavefunctions. This framework allows one to systematically convergence the results of SH parameters with respect to the basis set size and hence to establish robust results to compare with experiments for a chosen level of first principles theory. However, the method is computationally rather demanding and it only permitted the investigation of systems with tens of atoms.
Here we propose a computational scheme where only selected atoms are treated at the AE level, while the rest of the system is described within the PP approximation. We show that in order to obtain accurate SH parameters for spin defects, only a small number of atoms surrounding the defect (of the order of 10) needs to be treated at the AE level. Our approach permits calculations for cells with hundreds of atoms, as shown for the NV center in diamond and the VV in 4H-SiC, for which we used cells with up to 1022 atoms. In addition, using the cluster correlation expansion (CCE) 44 method, we demonstrate the importance of accurate SH parameters for precise predictions of coherence times of spin defects in semiconductors.

RESULTS
We carried out mixed AE-PP DFT calculations of hyperfine constants and zero-field splitting tensors using a finite-element (FE) basis and we compared the results to those obtained with PP calculations using plane waves (PW). A brief discussion on the FE basis is provided in the Supplementary Notes along with systematic convergence trends (Supplementary Table 1). In this work, we applied the mixed AE-PP approach to the negativelycharged NV center in diamond and neutral VV in 4H-SiC (see Fig. 1). Both defects have spin-triplet ground states. For VV we considered the axial kk configuration, where both the carbon and silicon vacancies are located at quasi-cubic sites. In addition, for coherence time calculations we also reported results for VV in 4H-SiC in the basal kh configuration, which exhibits avoided crossing of spin levels due to the lack of axial symmetry. We modeled the diamond and 4H-SiC using cubic and hexagonal supercells, respectively, except for kh VV in 4H-SiC, which we modeled using a orthorhombic supercell. For consistency, both PW and FE calculations were performed using structures relaxed with PW DFT. All calculations were performed with the Perdew-Burke-Ernzerhof exchange-correlational functional 45 . The use of the FE basis sets for SH parameter calculations permits an increase in resolution in the core regions of atoms, where wavefunctions exhibit a highly oscillatory behavior, while a coarser resolution suffices in the valence region. Due to its spatial adaptivity, the FE basis set is perfectly suited to carry out mixed AE-PP calculations. We briefly summarize below our strategy to evaluate hyperfine and zero-field splitting tensors, before presenting our results.
Computational framework for hyperfine and zero-field splitting tensors For each nuclear spin, the hyperfine A N -tensor is composed of an isotropic part A fc originating from the Fermi contact of electrons at the nuclei, and an anisotropic part A sd originating from spin dipolar interactions. A fc and A sd can be evaluated using the electron spin density of the system as: jr À r N j 5 n s ðrÞdr; (3) where n s (r) is the electronic spin density in real space, r N is the position of the nucleus of interest, and ðr À r N Þ a is the a-direction component of r − r N . S is the spin quantum number of the system (S = 0 for a singlet, 1 2 for a doublet, etc.), γ e and γ N are gyromagnetic ratios for electron and nuclei, respectively. We note that, while the Fermi contact term requires an accurate value of the electronic spin density exactly at the nucleus, the spin-dipole term requires an accurate description of the spin density in a region surrounding the nucleus, with the size of the region determined by the compact support of the spin density (i.e., the region where the spin density is non negligible) and by the decay of the kernel in Eq. (3). The spatial adaptivity of the FE basis is important for accurate descriptions of the spin density both in the core and valence regions, thus providing an efficient and systematic means of obtaining converged results, as described in the "Methods" section.
The spin-spin component of the zero-field splitting tensor D is evaluated as the expectation value of the magnetic dipole-dipole operator,~r 2 δ ab À 3rar b r 5 , using a Slater determinant built from Kohn-Sham orbitals 28,42 . The operator is essentially the ab element of the Hessian of the Green's function, Gðr; r 0 Þ ¼ 1 jr À r 0 j ; r is a scalar representing jr À r 0 j;r a is the a-th component of the vector r À r 0 . A straightforward evaluation of the D-tensor in realspace involves computing double integrals that are computationally very demanding. In our recent work 42 , we reformulated the evaluation of the D-tensor in real-space by solving a series of Poisson equations, as detailed in the "Methods" section.
Numerical accuracy of the mixed pseudopotential, all-electron approach In order to validate our approach for the calculation of the SH parameters using a mixed AE-PP approach, we first consider the Fermi contact term of the A-tensor for a NV center in a 3 × 3 × 3 diamond supercell containing 215 atoms, computed using the Γpoint for Brillouin zone sampling. Starting from a case where only the nitrogen atom is treated with an AE description, we gradually increase the number of atoms treated using AE calculations by considering eight cases (or levels), with 1,4,7,13,16,19,22, and 25 atoms near the defect treated at the AE level. As shown in Supplementary Fig. 2, when 16 neighbor atoms are treated at the AE level, which includes all the C atoms with dangling bonds Fig. 1 Spin density contours in two spin qubit systems. a, b Structures and spin densities of the negatively-charged nitrogen-vacancy center in diamond (a) and the axial kk-divacancy in 4H-SiC (b). In both systems, the spin density is localized around the carbon atoms with dangling bonds. By only treating a few atoms at the all-electron (AE) level, and the remaining using the pseudopotential (PP) approximation, an accurate prediction of spin Hamiltonian parameters can be achieved.
K. Ghosh et al. around the N atom, the value of the Fermi contact is converged (see Supplementary Table 2). In fact, even considering only the nitrogen atom and the three dangling bond carbon atoms at the AE level yields a value of the Fermi contact of −2.125 MHz, which is in very close agreement with the value of −2.096 MHz obtained by a full AE calculation. Our results suggest that by only treating a few atoms at the AE level, the mixed AE-PP calculation, henceforth denoted as FE-mixed, is adequate and accurate to obtain the Fermi contact term. Next, we increased the system size to a 4 × 4 × 4 supercell containing 511 atoms, and found that a mixed calculation with the same number of atoms treated at the AE level as in the 215 atom cell, accurately reproduces the Fermi contact term. This indicates that the number of atoms requiring an AE description in a FE-mixed calculation is independent of the system size. We also found that the spin-dipolar term of the A-tensor (cf. Supplementary Fig. 2) is much less sensitive to the number of atoms treated with an AE approach than the Fermi contact.
As previously noted, the evaluation of the D-tensor is computationally more demanding than that of the A-tensor, and a complete AE description becomes intractable even for a system with a few hundred atoms. Thus, to validate our mixed approach, we consider NV center in a 2 × 2 × 2 diamond supercell containing 63 atoms. We performed the mixed calculation with an AE description of the four atoms including the nitrogen atom and the three carbon atoms with dangling bonds. We obtained an excellent agreement for the values of the D-tensor obtained using FE-mixed and FE-AE calculations. Due to the C 3v symmetry of the system, the eigenvalues of the D-tensor We report 3 2 jD 3 j throughout this manuscript. We obtained a value of 2928.31 MHz in a mixed calculation, to be compared to 2939.47 MHz obtained from a full AE description. Thus, we conclude that the SH parameters obtained for the NV center in diamond by using a FE-mixed approach, where only the nitrogen atom and the three carbon atoms with dangling bonds are described at the AE level, are as accurate as those obtained with a calculation where all atoms are treated with an AE approach.
We consider next the VV in 4H-SiC in the kk configuration (referred to as VV-SiC) 46,47 . As in the case of NV-diamond, the VV-SiC has 3 carbon atoms adjacent to a silicon vacancy that contribute to the spin density of the system. Further, VV-SiC also has a silicon atom with 3 dangling bonds, adjacent to the carbon vacancy that give rise to mid-gap states. Based on our results for NV-diamond, we only treat these 6 atoms with an AE description. For validation purposes, we first considered a 4 × 4 × 1 supercell of VV-SiC containing 126 atoms, and we carried out two separate calculations, one using a complete AE description (FE-AE) and the other one using AE descriptions only for the atoms with dangling bonds. We found good agreement for the A-tensors computed with the two methods. Similar to the case of NV-diamond, we conclude that also for VV-SiC, treating a small number of atoms at the AE level (the three carbon and three silicon atoms with dangling bonds) suffices to obtain accurate results. In the case of VV-SiC, a full AE calculation of the D-tensor is not feasible with reasonable computational resources, and we limited our study to mixed calculations only.

Large scale calculations
We now turn to present results for A and D obtained with large supercells, using the FE-mixed method. The importance of large supercells to obtain accurate results of spin defects has been emphasized in several recent papers 9,48 . For the NV in diamond, we consider cubic cells ranging from 2 × 2 × 2 to 4 × 4 × 4 in size, containing 63 to 511 atoms. In the case of VV-SiC, we consider hexagonal cells ranging from 4 × 4 × 1 to 8 × 8 × 2 in size, containing 126 atoms to 1022 atoms.
The A-tensor and D-tensor for the various cell-sizes are shown in Figs. 2, 3, for NV-diamond and VV-SiC, respectively. For comparison, we also report the same parameters obtained using PW-PP DFT calculations. We found notable differences between the results obtained with the PW-PP method and with the FE-mixed approach, which does yield an overall AE accuracy, as shown above. In the case of VV-SiC (cf. Fig. 3), the errors associated to PW-PP values are as large as~20% for the value of the Fermi contact term of the A-tensor. In the case of the D-tensor, we performed PW calculations using pseudowavefunctions obtained with two PPs: GIPAW 49 and ONCV 50 , and found results that agree closely, but differ significantly from those obtained using FE-mixed calculations.
The ability to compute SH parameters for cells of different sizes provides important insights into the cell-size dependence of the results (see Supplementary Tables 3, 4). In order to understand the role of defect-defect interactions, we consider the A-tensor of the dangling bond carbon in NV-diamond computed using 2 × 2 × 2 supercells and 3 × 3 × 3 k-point sampling with that computed using 3 × 3 × 3 supercells and 2 × 2 × 2 k-point sampling. These calculations correspond to the same periodic boundary conditions (Born-von Karman boundary condition) for wavefunctions and thus, any change in the  Table 3).
A-tensor is solely due to cell-size effects arising from defect interactions. The Fermi contact value obtained in the two ways changes from 104.77 to 123.36 MHz. Similar cell-size effects arising from the defect-defect interaction are also evident in the case of VV-SiC (Supplementary Table 4). In the case of the D-tensor, while the cell-size effects are only marginal in NVdiamond, they are substantial for the VV-SiC (cf. Fig. 3). Our results underscore the importance of carrying out large scale calculations with AE accuracy to obtain accurate SH parameters.
Coherence time in weakly coupled nuclear-spin baths: the need for all-electron descriptions Having established an efficient approach to compute SH parameters, we now use the SH to compute dynamical properties of spin defects, in particular coherence times, for which we adopted the CCE method 51 . In CCE calculations, the coherence (ρ is the density matrix of the qubit andŜ þ is the spin raising operator) is approximated as a product of contributions from different nuclear-spin clusters. The coherence time is obtained by fitting the coherence function to a compressed exponential function, L % exp Àðt=T 2 Þ n . The convergence of the results is checked against the size of the spin bath, the number of clusters, and the maximum size of the cluster. Here, the ensemble-averaged coherence function is computed for configurations without nuclear spins with high hyperfine constants (A || > 1 MHz). In our CCE calculations, we used hyperfine tensors predicted from DFT calculations for nuclear spins included in the DFT supercell and adopted the point dipole approximation for nuclear spins at distances larger than those included in the supercell.
In general, there are two ways to control nuclear spins in defect systems. The strongly-coupled nuclear spins (with hyperfine coupling of order~1 MHz) can be directly accessed via radio frequency radiation. Here we define nuclear spins as strongly coupled when their hyperfine parameter is larger than the linewidth of the optically detected magnetic resonance 46 and separate oscillations in the Ramsey sequence due to the nuclear spin are observed 52 . These nuclear spins can be controlled with short gate times but they are highly susceptible to electron spin induced noise.
The second type of nuclear spins, weakly coupled to the electron spin (A ≪ 1 MHz) are controlled by dynamical decoupling schemes 53 . Applying refocusing pulses to the central spin may be used to not only increase the coherence time of the defect but also to isolate and control weakly coupled nuclear spins. These spins provide significantly longer coherence times than stronglycoupled spins, and the number of weakly coupled nuclear spins is not limited by the short distance to the central spin which is required for strong coupling.
The strength of the nuclear-spin induced dephasing mechanism, limiting coherence time T Ã 2 in spin defects is directly related to the A-tensor 54 , which therefore requires accurate calculations. Here, we first estimate the sensitivity of the inhomogeneous dephasing time T Ã 2 and the Hahn-echo coherence time T 2 on the values of the A-tensor by carrying out CCE calculations (Fig. 4); T 2 determines the stability of the qubit under dynamical decoupling schemes. We used two sets of hyperfine coupling computed for 4 × 4 × 4 supercell of the NV center in diamond: one set for which the A-tensor is calculated using the PAW reconstructed spin densities, and a second one based on the A-tensor calculated using FE-AE calculations. Figure 4b shows a histogram of the calculated A || (A z z in the defect reference frame) obtained using the two methods. We found that for large hyperfine coupling values, the relative difference is rather minor compared to the one for the smaller coupling terms. In order to see the impact of these differences on the dephasing of the electron spin, first we compute the ensemble-averaged dephasing time, T Ã 2 , by considering the decay of the coherence function averaged over a set of nuclear-spin configurations. The difference in the ensemble averages dephasing time was found not to be significant (1.35 (FE-AE) μs vs. 1.37 (PAW) μs). We note that the predicted value is close to the generally accepted value of nuclear-spin-limited T Ã 2 in diamond of~1 μs 55 . The dynamical decoupling changes the sensitivity of the qubit to the static noise, and the resulting Hahnecho coherence time T 2 is found to be insensitive to the choice of hyperfine couplings and equal to 0.89 ms, in agreement with previously reported CCE predictions 14 .
Next, we focus on the single-defect dephasing time in the weakly coupled nuclear-spin bath considering the example to NV center in diamond (Fig. 4a). We select nuclear configurations whose Fourier transform of the free induction decay contains only one peak. This procedure ensures that the coherence function of the defects in the chosen subset of nuclear configurations does not contain any oscillations due to the nuclear spins with high hyperfine coupling. Hence the procedure guarantees that the  Table 4).
nuclear baths contain only weakly coupled nuclear spins and we can thus use an exponential decay fit to obtain the value of T Ã 2 . Figure 4c shows the T Ã 2 for the chosen subset of nuclear configurations. In these systems, the difference between PAW and FE-AE results is significant as the dephasing time differs by a factor of 2 for certain nuclear configurations. We note that in general for weakly coupled spins, the hyperfine coupling computed with PAW are higher (as can be seen in Fig. 4b compared to their FE-AE counterpart, leading to significant differences in the predictions of the dephasing time. At the same time the Hahn-echo coherence time (Fig. 4c, right pane) does not depend on the choice of hyperfine couplings, as the main dephasing mechanism under dynamical decoupling is pairwise spin flip of distant spin pairs 14 .
However, in the emerging qubit systems based on the clock transitions at avoided crossings of spin levels [56][57][58] , the sensitivity of the coherence time to the magnetic noise from the nearby nuclear spins is different. As an example, we consider the basal kh divacancy in 4H-SiC (Fig. 4d). Compared to kk-VV, the basal divacancy has a lower symmetry and ZFS tensor contains both axial splitting D 3 and basal splitting E = D 2 − D 1 = 18.4 MHz 59 , leading to emergence of avoided crossing in spin levels at zero magnetic field. Using hyperfine couplings computed with PAW and FE-AE in 6 × 3 × 2 orthorhombic supercell containing 574 atoms (Fig. 4e) we predict the ensemble average T Ã 2 PAW ¼ 0:13 ms, T Ã 2 FE ¼ 0:14 ms and T 2 PAW = 1.04 ms, T Ã 2 FE ¼ 1:11 ms using generalization of CCE method 60 . The PW-PAW results consistently overestimate the hyperfine couplings of the weakly coupled nuclear spins, which leads to notable differences even in the ensemble-averaged coherence times of this system. Both values are smaller than the ones, reported in the previous work (T Ã 2 ¼ 0:16 ms, T 2 = 1.15 ms) 60 where we used larger 9 × 5 × 2 orthorhombic supercell containing 1438 atoms to compute hyperfines with PW-PAW approach, which confirms the importance of larger supercell for precise calculations of hyperfine couplings. Finally, we analyze the coherence times of basal divacancy at clock transition in the single nuclear-spin spatial configurations (Fig. 4f).
We find that both T Ã 2 and T 2 depend significantly on the approximation used to obtain hyperfine couplings of the nearby nuclear spins, leading to significant differences in predicted coherence times.

DISCUSSION
In summary, we presented an efficient approach based on realspace DFT to compute SH parameters. Our approach treats a few selected atoms close to a defect of interest with AE accuracy and uses the PP approximation for the remaining atoms; the method uses a FE basis set and is systematically convergent with respect to the basis set size. We presented results demonstrating that the accuracy obtained in the computation of the SH parameters using the mixed AE-PP approach is commensurate with that of full AE calculations. The mixed approach in conjunction with the spatial adaptivity of the FE basis enabled the computation of hyperfine and zero-field splitting tensors for large systems, containing~1000 atoms, with AE accuracy. Our results revealed significant cell-size effects in the calculations of both the hyperfine and the zero-field splitting interaction, indicating that finite size scaling is important in determining accurate values for these tensors.
We computed hyperfine constants for both strongly and weakly coupled nuclei spins from AE calculations in two systems, NV in diamond and kh-VV in SiC, followed by coherence times estimation in those two systems. This shows that the relative difference between AE and PP predictions of hyperfine tensors for strongly-coupled nuclear spins (those for which A ≥ 1 MHz) is small. However, absolute differences in hyperfine tensors predicted with PW-PP and AE methods for weakly coupled spins (A ≪ 1 MHz), even when similar in magnitude to those found for strongly-coupled spins, may dramatically impact the prediction of T Ã 2 of the NV center, with differences up to a factor of 2 for certain nuclear-spin configurations. In the case of the clock transitions of the basal kh-VV in SiC, the variance in the predictions is even more drastic, with relative difference up to a factor of 4. In this case, the choice of the approximation affects both T 2 and T Ã 2 coherence times. We note that, in addition to coherence time calculations, accurate predictions of zero-field splitting and hyperfine tensors for strongly-coupled nuclear spins are important to identify the atomistic structure of spin defects 19 ; furthermore, accurate predictions of hyperfine tensors for weakly coupled nuclear spins are key for the spatial mapping of experimental multinuclear registers 61 and the prediction of plausible memory units in spin centers 52 .
Finally we note that here we focused on hyperfine coupling and work is in progress to address the accuracy of other parameters entering the definition of the SH. For example, it has been recently shown that the addition of quadrupole splitting for nuclei with spin ≥1 both three-and two-dimensional materials can significantly impact the coherence time of the spin qubits 62 . The calculations of D tensors including spin-orbit interaction are usually computed using AE GTO basis sets for clusters of atoms, rather than conducting calculations in the condensed phase 63 , but recently a PAW formalism has been proposed 64 which would be interesting to compare with FE calculations.
The method introduced here, which may be applied also to the calculations of SH parameters in addition to the hyperfine coupling, represents a substantial progress towards accurate computations of SH parameters of spin defects in large scale condensed and molecular systems, paving the way to highthroughput screening of spin defects for quantum information technologies.

METHODS PAW calculations of SH parameters
PW based DFT calculations are performed with the Quantum ESPRESSO code 65 with a kinetic energy cutoff of 75 Ry. GIPAW PP developed by Ceresoli 49 were used for the calculation of the Aand D-tensor. Moreover, the D-tensor calculations were also carried out using ONCV PP 50 . After solving the Kohn-Sham equations, the A-tensor is evaluated by PAW reconstruction using the GIPAW code 27 , while the D-tensor calculations are performed with the PyZFS code 66 using normalized pseudowavefunctions 9,15,67,68 . In A-tensor calculations, we considered nuclear isotopes 13 C, 14 N, and 29 Si for C, N, and Si, respectively.

FE based calculations of SH parameters
FE based DFT calculations are carried out using the DFT-FE code 43 . The convergence of the SH parameters with respect to the FE mesh discretization parameters is studied on the smallest supercell for each system considered here. Our convergence study was carried out with respect to the FE polynomial order and the FE mesh element size in the vicinity of the AE atoms (cf. Supplementary Table 1 for convergence of A-tensor of nitrogen atom in NV center of a 2 × 2 diamond supercell). Both the polynomial and mesh size requirements are more stringent for the A-tensor calculations compared to the D-tensor ones, as the former depends on the cusps of the spin density at the nucleus. More details related to convergence properties can be found in our previous work 42 .

Computational specifics of FE based A-tensor calculation
The computation of the A-tensor uses the self-consistent spin density obtained from the DFT calculation in the following manner. The spin density is evaluated at the FE quadrature points during a self-consistent DFT calculation (as quadrature values are directly used in the ensuing evaluation of integrals), while the Kohn-Sham wavefunctions are computed at the FE nodal points (cf. 43 for details; also see Supplementary Information for a brief discussion on FE discretization). The Fermi contact term, which requires the value of the spin density at the FE nodes (where the atoms of interest are located), is computed by reevaluating the spindensity at these relevant FE nodes from the Kohn-Sham wavefunctions. The available quadrature point values of the spin density are used for the evaluation of the spin dipolar term, which involves an integral over the spin density (cf. Eq (3)). However, instead of directly carrying out the integral in Eq. (3), we evaluate the integral as the Hessian of the potential (Φ s (r)) resulting from the spin density, where Φ s (r) is computed as a solution of the partial differential equation (PDE), ∇ 2 Φ s (r) = − 4πn s (r) with periodic boundary conditions. This reformulation accounts for the contributions of periodic images, which otherwise would have to be explicitly computed in the direct evaluation of the integral in Eq. (3). We note that in order to ensure spin neutrality in the solution of the Poisson equation with periodic boundary conditions, a uniform and opposite background with spin density (−M/Ω, where M is the net magnetization of the system, and Ω is the volume of the supercell) is added to n s (r). It is straightforward to show that the contribution of such a uniform background on the Hessian of Φ s (r) is exactly zero.

Computational specifics of FE based D-tensor calculation
The calculation of the D-tensor in real-space is cast into the solution of a series of Poisson equations as described below.
ðγ e _Þ 2 X occ: with where ϕ i (r) is the ith single electron wavefunction obtained from DFT, and Λ ij b ðrÞ is the solution of the Poisson equation . The superscript D and E represent the terms corresponding to direct and exchange interactions. χ ij = ± 1 for ith and jth wavefunctions having parallel and antiparallel spins, respectively. The most expensive part of the calculation of the D-tensor is the solution of the NðN þ 1Þ 2 Poisson equations, where N is the number of electrons in the system. These Poisson equations are solved using the same FE discretization used to obtain the Kohn-Sham wavefunctions. Hence the reduction in the degrees of freedom achieved via a spatially adaptive FE discretization also aids in improving the computational efficiency of solving the Poisson equations. In addition, Poisson equations corresponding to different pairs of wavefunctions may be solved in parallel since they are independent from each other. The PDE involved in the D-tensor calculation takes the form ∇ 2 Λ ij b ¼ À4π ∂ðϕ Ã i ðrÞϕ j ðrÞÞ ∂rb . Here, the right hand side of the PDE is computed by first evaluating ϕ Ã i ðrÞϕ j ðrÞ on the FE nodes and then interpolating the gradient to the quadrature points that are subsequently used in the ensuing integrals of the weak formulation of the FE method. The PDEs are solved using the framework for the Poisson problem in the self-consistent cycle of DFT. The potential fields (Λ ij;D b ðrÞ and Λ ij;E b ðrÞ) obtained from the PDE are defined on the FE nodes and are interpolated to quadrature points before carrying out the integrations in Eq. (5). The interpolation errors systematically decrease with increasing FE polynomial order and decreasing FE mesh size.
Computational cost benefit of the mixed framework As we showed above, the number of relevant atoms that need to be treated at the AE level is much smaller (<10) than the total number of atoms in the system, and does not scale with the system size. We note that for systematically convergent calculations, the FE basis functions in AE calculations are~10-fold larger than in PP calculations for the elements considered here. For the benchmark calculations presented in this work, mixed AE-PP calculations reduce the number of basis functions by~80% compared to full AE calculations. This advantage of mixed calculations becomes increasingly important for heavier atoms. Importantly, the DFT-FE code has the ability to treat both AE and PP calculations in the same framework. In DFT-FE, the solution to the Kohn-Sham equations is obtained via the Chebyshev polynomial filtered subspace iteration (ChFSI) procedure 69,70 . The various computational steps in the ChFSI procedure scales as OðMNÞ, OðMN 2 Þ, and