All-electron density functional calculations for electron and nuclear spin interactions in molecules and solids

TL;DR

This paper presents a first-principles approach using all-electron density functional theory with finite element methods to accurately compute spin interaction parameters in molecules and solids, including hyperfine and quadrupole tensors.

Contribution

It introduces a systematic all-electron FE-based DFT method for calculating spin Hamiltonian parameters, providing a benchmark for accuracy and convergence.

Findings

FE calculations can be systematically converged for spin parameters

Results agree well with Gaussian and plane-wave methods

Establishes reference values for spin Hamiltonian parameters

Abstract

The interaction between electronic and nuclear spins in the presence of external magnetic fields can be described by a spin Hamiltonian, with parameters obtained from first principles, electronic structure calculations. We describe an approach to compute these parameters, applicable to both molecules and solids, which is based on Density Functional Theory (DFT) and real-space, all-electron calculations using finite elements (FE). We report results for hyperfine tensors, zero field splitting tensors (spin-spin component) and nuclear quadrupole tensors of a series of molecules and of the nitrogen-vacancy center in diamond. We compare our results with those of calculations using Gaussian orbitals and plane-wave basis sets, and we discuss their numerical accuracy. We show that calculations based on FE can be systematically converged with respect to the basis set, thus allowing one to establish reference values for the spin Hamiltonian parameters, at a given level of DFT.

Figures

Tables

Table 1. Table 1: Isotropic hyperfine tensor (see Eq. 2 ) (MHz) computed by DFT with finite-element (FE), Gaussian-type orbital (GTO) and plane-wave (PW) basis sets. For PW calculations three different treatments of the core-relaxation effect are considered, which include Slater exchange (Slater-X), exchange (X), and exchange+correlation (XC) in the perturbative potential for the calculation of spin densities at the core region.

System	Atom	FE	GTO(EPR-III)	GTO(IGLO-III)	PW(Slater-X)	PW(X)	PW(XC)	Exp
CN ($S=\frac{1}{2}$)	${}_{}{}^{13}\text{C}$	504.21	500.50	509.63	539.55	536.04	566.57	588 [Easley and Weltner Jr, 1970]
	${}_{}{}^{14}\text{N}$	-12.81	-12.47	-12.25	-14.87	-15.15	-12.43	-13 [Easley and Weltner Jr, 1970]
BO ($S=\frac{1}{2}$)	${}_{}{}^{11}\text{B}$	1007.71	998.17	1002.34	983.68	980.02	1009.31	1027 [Tanimoto et al., 1986]
	${}_{}{}^{17}\text{O}$	-7.34	-7.18	-7.26	-7.83	-8.13	-7.24
AlO ($S=\frac{1}{2}$)	${}_{}{}^{27}\text{Al}$	590.80		646.95	564.18	560.09	626.93	766 [Knight Jr and Weltner Jr, 1971]
	${}_{}{}^{17}\text{O}$	11.47		12.20	-15.22	-14.78	-22.99
NH ($S=1$)	${}_{}{}^{14}\text{N}$	11.27	10.20	9.77	24.24	22.00	33.60	20 [Weltner, 1989]
	${}_{}{}^1\text{H}$	-53.52	-53.13	-47.74	-51.60	-51.60	-51.60	-70 [Weltner, 1989]
Diamond NV ($S=1$)	${}_{}{}^{14}\text{N}$	-2.32			-2.60	-2.60	-2.56	2.23 [He et al., 1993], -2.51 [Felton et al., 2009], -2.53 [Yavkin et al., 2016]
	DB ${}_{}{}^{13}\text{C}$	98.72			100.27	99.05	108.51	146.7 [Felton et al., 2009]

Table 2. Table 2: Spin dipolar hyperfine tensor (see Eq. 3 ) (MHz) computed by DFT with finite-element (FE), Gaussian-type orbital (GTO) and plane-wave (PW) basis sets. The eigenvalue with the largest absolute value is shown.

System	Atom	FE	GTO(EPR-III)	GTO(IGLO-III)	PW	Exp
CN ($S=\frac{1}{2}$)	${}_{}{}^{13}\text{C}$	115.33	118.47	117.43	124.20	89.9 [Easley and Weltner Jr, 1970]
	${}_{}{}^{14}\text{N}$	44.51	42.62	42.40	45.25	30.8 [Easley and Weltner Jr, 1970]
BO ($S=\frac{1}{2}$)	${}_{}{}^{11}\text{B}$	53.71	53.38	53.76	55.37	54.254 [Tanimoto et al., 1986]
	${}_{}{}^{17}\text{O}$	-47.83	-46.47	-45.97	-51.55
AlO ($S=\frac{1}{2}$)	${}_{}{}^{27}\text{Al}$	114.23		111.67	112.51	106 [Knight Jr and Weltner Jr, 1971]
	${}_{}{}^{17}\text{O}$	-122.65		-116.42	-127.57
NH ($S=1$)	${}_{}{}^{14}\text{N}$	-47.87	-45.82	-46.01	-49.59	-46 [Weltner, 1989]
	${}_{}{}^1\text{H}$	58.08	58.50	59.92	58.02	60 [Weltner, 1989]
Diamond NV ($S=1$)	${}_{}{}^{14}\text{N}$	-0.07			-0.05	-0.13 [He et al., 1993], 0.37 [Felton et al., 2009], 0.33 [Yavkin et al., 2016]
	DB ${}_{}{}^{13}\text{C}$	54.87			58.34	52.9 [Felton et al., 2009]

Table 3. Table 3: The spin-spin component of the zero-field splitting tensor (see Eq. 4 ) ( cm − 1 superscript cm 1 \text{cm}^{-1} ) computed by DFT with finite-element (FE), Gaussian-type orbital (GTO) and plane-wave (PW) basis sets. Scalar parameters D = 3 2 ​ D 33 𝐷 3 2 subscript 𝐷 33 D=\frac{3}{2}D_{33} are reported. Scalar parameters E = 1 2 ​ ( D 22 − D 11 ) 𝐸 1 2 subscript 𝐷 22 subscript 𝐷 11 E=\frac{1}{2}(D_{22}-D_{11}) , if non-zero, are reported in brackets.

System	FE	GTO(EPR-III)	GTO(IGLO-III)	PW(GIPAW)	PW(ONCV)	Exp
${\text{O}}_2$ ($S=1$)	1.894	1.893	1.843	1.642	1.695	3.960 [Huber, 1979]a
${\text{C}\text{H}}_2$ ($S=1$)	0.894 (0.051)	0.895 (0.052)	0.895 (0.052)	0.896 (0.052)	0.908 (0.049)	0.760 (0.068) [Wasserman et al., 1971]
NH ($S=1$)	1.861	1.862	1.857	1.795	1.815	1.860 [Dixon, 1959]
${\text{C}}_5{\text{H}}_5^{+}$ ($S=1$)	0.123	0.120	0.110	0.123	0.123	0.187 [Saunders et al., 1973]
Diamond NV ($S=1$)	0.100			0.101	0.100	0.096 [Loubser and van Wyk, 1978]

Table 4. Table 4: Quadrupole coupling constants e ​ Q ​ V 33 𝑒 𝑄 subscript 𝑉 33 eQV_{33} (see Eq. 13 ) (MHz) computed by DFT with finite-element (FE), Gaussian-type orbital (GTO) and plane-wave (PW) basis sets. The scalar parameters η = | V 22 − V 11 V 33 | 𝜂 subscript 𝑉 22 subscript 𝑉 11 subscript 𝑉 33 \eta=\left|\frac{V_{22}-V_{11}}{V_{33}}\right| , if non-zero, are reported in brackets.

System	Atom	FE	GTO(EPR-III)	GTO(IGLO-III)	PW	Exp
HCN ($S=0$)	${}_{}{}^2\text{H}$	0.20	0.21	0.22	0.20
	${}_{}{}^{11}\text{C}$	-2.95	-2.87	-2.77	-3.73
	${}_{}{}^{14}\text{N}$	-4.72	-4.85	-4.79	-5.07	-4.02 [Latosińska, 2003]
NCCN ($S=0$)	${}_{}{}^{14}\text{N}$	-4.72	-4.80	-4.75	-4.99	-4.27 [Latosińska, 2003]
	${}_{}{}^{11}\text{C}$	-1.84	-1.95	-1.76	-2.72
${\text{N}}_2$ ($S=0$)	${}_{}{}^{14}\text{N}$	-5.08	-5.52	-5.46	-5.77	-4.65 [Latosińska, 2003]
${\text{H}}_2\text{O}$ ($S=0$)	${}_{}{}^2\text{H}$	0.31 (0.12)	0.31 (0.13)	0.32 (0.13)	0.30 (0.13)	0.31 (0.14) [Verhoeven et al., 1969]
	${}_{}{}^{17}\text{O}$	10.10 (0.81)	10.69 (0.74)	10.90 (0.74)	10.85 (0.72)	10.17 (0.75) [Verhoeven et al., 1969]
Diamond NV ($S=1$)	${}_{}{}^{14}\text{N}$	-7.30			-7.47	-6.68 [Felton et al., 2009]a
	DB ${}_{}{}^{11}\text{C}$	1.96 (0.03)			-1.77 (0.10)

Table 5. Table 5: The principle values of A 𝐴 A (top, in MHz), D 𝐷 D (center, in cm − 1 superscript cm 1 \text{cm}^{-1} ), and V 𝑉 V (bottom, in a.u.) tensors for the NV center computed by FE-based DFT using different finite-element polynomial degrees.

	14N				DB 13C
	Adip11	Adip22	Adip33	Aiso	Adip11	Adip22	Adip33	Aiso
FE3	0.033	0.033	-0.066	-2.362	-27.303	-27.592	54.896	100.492
FE4	0.034	0.034	-0.067	-2.307	-27.199	-27.654	54.853	99.708
FE5	0.035	0.035	-0.070	-2.319	-27.189	-27.664	54.854	99.016
FE6	0.035	0.035	-0.070	-2.316	-27.171	-27.696	54.867	98.721

Equations

$$H={\mu }_B\mathbf{B}\cdot \mathbf{g}\cdot \mathbf{S}+\sum _N{\gamma }_N\mathbf{B}\cdot {\mathbf{I}}_N+\sum _N\mathbf{S}\cdot {\mathbf{A}}_N\cdot {\mathbf{I}}_N+\mathbf{S}\cdot \mathbf{D}\cdot \mathbf{S}+\sum _N{\mathbf{I}}_N\cdot {\mathbf{P}}_N\cdot {\mathbf{I}}_N$$

$$A^{\text{iso}}=-\frac{1}{3S}{\mu }_0{\gamma }_e{\gamma }_N{\hslash }^2{n}_s({\mathbf{r}}_N),$$

$$A_{ab}^{\text{dip}}=\frac{1}{2S}\frac{{\mu }_0}{4\pi }{\gamma }_e{\gamma }_N{\hslash }^2\int \frac{|\mathbf{r}-{\mathbf{r}}_N{|}^2{\delta }_{ab}-3(\mathbf{r}-{\mathbf{r}}_N{)}_a(\mathbf{r}-{\mathbf{r}}_N{)}_b}{|\mathbf{r}-{\mathbf{r}}_N{|}^5}{n}_s(\mathbf{r})d\mathbf{r},$$

$$D_{ab}=\frac{1}{2S(2S-1)}\frac{{\mu }_0}{4\pi }({\gamma }_e\hslash {)}^2\left[\sum _{i<j}^{occ.}{\chi }_{ij}\int \int {\Phi }_{ij}^{\ast }(\text{r},{\text{r}}^{{}^{\prime }})\frac{{\tilde{r}}^2{\delta }_{ab}-3{\tilde{r}}_a{\tilde{r}}_b}{r^5}{\Phi }_{ij}(\text{r},{\text{r}}^{{}^{\prime }})d\text{r}d{\text{r}}^{{}^{\prime }}\right],$$

$$M_{ab}^{ij,D}=\int \int {\varphi }_i(\text{r}){\varphi }_j({\text{r}}^{{}^{\prime }})\frac{{\partial }^{2}G(\text{r},{\text{r}}^{{}^{\prime }})}{\partial r_a\partial r_b^{{}^{\prime }}}{\varphi }_i^{\ast }(\text{r}){\varphi }_j^{\ast }({\text{r}}^{{}^{\prime }})d\text{r}d{\text{r}}^{{}^{\prime }},$$

$$M_{ab}^{ij,E}=\int \int {\varphi }_i(\text{r}){\varphi }_j({\text{r}}^{{}^{\prime }})\frac{{\partial }^{2}G(\text{r},{\text{r}}^{{}^{\prime }})}{\partial r_a\partial r_b^{{}^{\prime }}}{\varphi }_i^{\ast }({\text{r}}^{{}^{\prime }}){\varphi }_j^{\ast }(\text{r})d\text{r}d{\text{r}}^{{}^{\prime }}.$$

$$M_{ab}^{ij,D}=\int \int \frac{\partial ({\varphi }_i(\text{r}){\varphi }_i^{\ast }(\text{r}))}{\partial r_a}G(\text{r},{\text{r}}^{{}^{\prime }})\frac{\partial ({\varphi }_j({\text{r}}^{{}^{\prime }}){\varphi }_j^{\ast }({\text{r}}^{{}^{\prime }}))}{\partial r_b^{{}^{\prime }}}d\text{r}d{\text{r}}^{{}^{\prime }},$$

$$M_{ab}^{ij,E}=\int \int \frac{\partial ({\varphi }_i(\text{r}){\varphi }_j^{\ast }(\text{r}))}{\partial r_a}G(\text{r},{\text{r}}^{{}^{\prime }})\frac{\partial ({\varphi }_i^{\ast }({\text{r}}^{{}^{\prime }}){\varphi }_j({\text{r}}^{{}^{\prime }}))}{\partial r_b^{{}^{\prime }}}d\text{r}d{\text{r}}^{{}^{\prime }}.$$

$$M_{ab}^{ij,D}=\int \frac{\partial ({\varphi }_i(\text{r}){\varphi }_i^{\ast }(\text{r}))}{\partial r_a}{\Lambda }_b^{jj,D}(\text{r})d\text{r}$$

$$M_{ab}^{ij,E}=\int \frac{\partial ({\varphi }_i(\text{r}){\varphi }_j^{\ast }(\text{r}))}{\partial r_a}{\Lambda }_b^{ij,E}(\text{r})d\text{r},$$

$$D_{ab}=\frac{1}{2S(2S-1)}\frac{{\mu }_0}{4\pi }({\gamma }_e\hslash {)}^2\sum _{i<j}^{occ.}{\chi }_{ij}(M_{ab}^{ij,D}-M_{ab}^{ij,E}).$$

$$\begin{array}{rl}{H}_Q& =\mathbf{I}\cdot \mathbf{P}\cdot \mathbf{I}\\ & =\frac{eQ}{6I(2I-1)}\sum _{a,b}{V}_{ab}\left[\frac{3}{2}(I_a{I}_b+I_b{I}_a)-{\delta }_{ab}I(I+1)\right],\end{array}$$

$$\begin{split}V_{ab}&=[\nabla_{a}\nabla_{b}V(\bm{r})]|_{\bm{r}=\bm{r}_{N}}\\ &=\left\{\nabla_{a}\nabla_{b}\left[-\int d\bm{r}^{\prime}\frac{n(\bm{r}^{\prime})}{|\bm{r}-\bm{r}^{\prime}|}+\sum_{I\neq N}\frac{Z_{I}}{|\bm{r}-\bm{r}_{I}|}\right]\right\}\Biggr{|}_{\bm{r}=\bm{r}_{N}}\end{split}$$

Full text

I Introduction

Electron spins in molecules, nanostructures and solids are important resources in many areas including spintronics Žutić et al. (2004) and quantum information science Weber et al. (2010). For instance, high-spin magnetic molecules can be utilized as single-molecule magnets and are promising platforms for next-generation data storage devices Guo et al. (2018); in the solid state, spin-carrying deep centers in semiconductors can serve as quantum bits for quantum information processing Koehl et al. (2015). In order to understand the physical properties of electron spins in molecules and solids, one needs to describe the interaction of electron and nuclear spins, in the presence of external electromagnetic fields. Such a description may be achieved by using spin Hamiltonians, with parameters derived from experiments or from calculations. For systems with a single effective electron spin, the leading terms in the spin Hamiltonian are Schweiger and Jeschke (2001); Harriman (2013); Abragam and Bleaney (2013):

[TABLE]

where $\mu_{B}$ is Bohr magneton; $\bm{S}$ is the effective electron spin; $\bm{B}$ is the external magnetic field; $\bm{I}_{N}$ and $\gamma_{N}$ are the spin and gyromagnetic ratio of the $N^{\text{th}}$ nucleus; $\bm{g}$, $\bm{A}$, $\bm{D}$, and $\bm{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 spin Hamiltonian parameters $\bm{g}$, $\bm{A}$, $\bm{D}$ and $\bm{P}$ may be obtained by electron paramagnetic resonance (EPR), nuclear quadrupole resonance (NQR) and related spectroscopic techniques Weil and Bolton (2007). Theoretically their values can be determined by first-principles electronic structure calculations, which also provide important information complementary to experiments. For example, in the case of spin defects in solids often times the atomistic structure and charge state of the defect are not straightforward to determine, experimentally. Comparing the computed spin Hamiltonian parameters for candidate structures and charge states with experimental results is a useful means to identify the properties of the defect. In addition, first-principles calculations can provide insights into the structure-property relations of molecules and spin defects, thus facilitating the rational design of molecules and materials with desirable spin properties. Finally, by simulating spin systems under external perturbations such as mechanical strain or applied electromagnetic fields, one can obtain valuable information and guidance for the experimental manipulation of electron spins Falk et al. (2014); Whiteley et al. .

Therefore, in order to devise predictive computational strategies, the development of robust methods for the calculation of spin Hamiltonian parameters is an important task. In spite of important progress in the fields of materials science Van de Walle and Blöchl (1993); Blügel et al. (1987); Pickard and Mauri (2001, 2002); Bahramy et al. (2007); Rayson and Briddon (2008); Bodrog and Gali (2013); Biktagirov et al. (2018) and quantum chemistry Olsen et al. (2002); Latosińska (2003); Sinnecker and Neese (2006); Neese (2005); Reviakine et al. (2006); Kossmann et al. (2007); Neese (2007), there is not yet a general and well established computational protocol that can reliably predict various spin Hamiltonian parameters with high accuracy for broad classes of systems. At present, the method most often adopted for spin Hamiltonian parameter calculations is Density Functional Theory (DFT). While calculations using ab initio wavefunction-based methods have also appeared in the literature Sayfutyarova and Chan (2018); Sugisaki et al. (2009), they have so far been limited to relatively small molecular systems due to their high computational cost. To solve the Kohn-Sham equations in DFT, single particle electronic wavefunctions are usually represented using basis sets, with Gaussian-type orbitals (GTO) and plane-waves (PW) being among the most popular choices for molecular and extended systems, respectively. In PW-based DFT calculations, pseudopotentials are employed and the electronic wavefunctions near the nuclei are not explicitly evaluated, and one generally needs to perform a so-called projected augmented wave (PAW) reconstruction Blöchl (1994) to extract all-electron wavefunctions for the calculation of certain spin Hamiltonian parameters. Besides PW, there are studies exploring other basis sets including numerical atomic orbitals Kadantsev and Ziegler (2008), linearized augmented plane-wave Schwarz and Blaha (2003), linear muffin-tin orbitals Daalderop et al. (1996); Overhof and Gerstmann (2004), and GTO Dovesi et al. (2018) for the calculation of $A$-tensors and $V$-tensors (electric field gradient tensor, see Section II.3) for solids.

In this work we present calculations of spin Hamiltonian parameters carried out, for the first time, using a real-space finite-element (FE) formulation of DFT Motamarri et al. (2013); dft . The FE basis is a piece-wise continuous polynomial basis Brenner and Scott (2009) that allows for systematic convergence of calculations with increasing polynomial order and decreasing element size. An important attribute of the FE basis is its spatial adaptivity that can provide increased resolution in specific regions of interest in real space, while using coarser descriptions elsewhere. In the present context, the FE basis can be chosen to have higher resolution in the core region to accurately describe the highly oscillatory nature of the single particle wavefunctions, and a coarser resolution far from the core where the orbitals are smoother. Further, FE-based calculations can be performed with either open or periodic boundary conditions, and therefore molecular and extended systems can be treated on an equal footing. There are several advantages in using FE-based DFT calculations for computing spin Hamiltonian parameters. The cusp of the wavefunctions near the nuclei can be more efficiently represented than with GTO basis sets, and this is an important requisite to compute quantities such as the Fermi contact component of the $A$-tensor. In addition, FE-based calculations can be systematically converged with respect to the basis set size in a more straightforward manner than GTO-based calculations, and they do not mandate the use of pseudopotentials and PAW reconstructions, as required when using PWs.

The rest of the paper is organized as follows. In Section II we present the formalism for computing the spin Hamiltonian parameters in Eq. 1, except those involving spin-orbit coupling contributions. Specifically, we consider the isotropic (Fermi contact) and the spin dipolar contribution to the hyperfine $A$-tensor, the spin-spin component of the zero-field splitting $D$-tensor, and the nuclear quadrupole $P$-tensor. The calculation of the $g$-tensor and the spin-orbit contributions to the $A$- and $D$-tensors will be subjects for future studies. In Section III we discuss our results for both molecules and solids, and Section IV concludes the paper.

II Real space computation of spin Hamiltonian parameters

The adaptivity of the FE basis to accurately and efficiently describe the all-electron Kohn-Sham orbitals is determined by the size of finite elements and the order of the polynomial basis functions. The FE mesh can be chosen adaptively (h-refinement), so as to assign smaller elements to regions requiring higher resolution (e.g. around the nuclei) and coarser elements elsewhere. In addition, for a given mesh, the order of the polynomial basis functions can be chosen (p-refinement) to provide a high order function approximation. While both refinement approaches can be used to obtain systematic convergence, h-refinement is suitable for realizing spatial adaptivity, and p-refinement provides the flexibility to realize higher accuracy over the complete simulation domain.

II.1 A-tensor

The isotropic (Fermi contact) and the spin dipolar component of the $A$-tensor are given by:

[TABLE]

[TABLE]

where $a,b=x,y,z$, S is the effective electron spin ($S=0$ for a singlet, $\frac{1}{2}$ for a doublet, etc.); $n_{s}$ is the electron spin density; $\bm{r}_{N}$ is the position of the nucleus; $(\bm{r}-\bm{r}_{N})_{a}$ is the $a$-direction component of $\bm{r}-\bm{r}_{N}$; $\gamma_{e}$ and $\gamma_{N}$ are gyromagnetic ratios for electron and nuclei, respectively. $\gamma_{N}$ for various nuclear isotopes can be obtained from standard databases, e.g. Ref. Web, .

As can be seen from Eq. 2, the isotropic (Fermi contact) component of the $A$-tensor exhibits a strong dependence on the electron spin density at the nuclei. An all-electron A-tensor calculation in real-space requires very refined finite elements near the nuclei to accurately compute the electron spin density. The spatial adaptivity of the finite element mesh (h refinement) is hence extremely useful here. On the other hand, the dipolar component of the A-tensor involves an integration with $\frac{1}{r^{5}}$ and $\frac{1}{r^{3}}$ kernels. This requires high accuracy in the electronic spin density within a certain region surrounding the nuclei, which can be systematically improved through the p-refinement.

II.2 D-tensor

The spin-spin component of the $D$-tensor evaluated using the Kohn-Sham wavefunctions, is given by McWeeny (1961); Rayson and Briddon (2008)

[TABLE]

where the summation is over all pairs of occupied orbitals, and $\Phi_{ij}(\textbf{r},\textbf{r}^{{}^{\prime}})$ are $2\times 2$ determinants formed from orbitals $\phi_{i}$ and $\phi_{j}$, \Phi_{ij}(\textbf{r},\textbf{r}^{{}^{\prime}})=\frac{1}{\sqrt{2}}\Big{[}\phi_{i}(\textbf{r})\phi_{j}(\textbf{r}^{{}^{\prime}})-\phi_{i}(\textbf{r}^{{}^{\prime}})\phi_{j}(\textbf{r})]; $\chi_{ij}=\pm 1$ for parallel and antiparallel spins respectively; $\tilde{r}$ is a scalar representing $|\textbf{r}-\textbf{r}^{{}^{\prime}}|$; $\tilde{r}_{a}$ represents the $a$-direction component of the vector $\textbf{r}-\textbf{r}^{{}^{\prime}}$. The operator $\frac{\tilde{r}^{2}\delta_{ab}-3\tilde{r}_{a}\tilde{r}_{b}}{r^{5}}$ is the $ab$ element of the Hessian of the Green’s function of $-\frac{1}{4\pi}\nabla^{2}$, i.e. $G(\textbf{r},\textbf{r}^{{}^{\prime}})=\frac{1}{|\textbf{r}-\textbf{r}^{{}^{\prime}}|}$. Since the operator, $\frac{\partial^{2}G(\textbf{r},\textbf{r}^{{}^{\prime}})}{\partial r_{a}\partial r^{{}^{\prime}}_{b}}$, is invariant under particle exchange, the real-space integrals in Eq. 4 can be split into direct ($M^{ij,D}_{ab}$) and exchange terms ($M^{ij,E}_{ab}$) given by

[TABLE]

and

[TABLE]

Equation 5 and Eq. 6 can be rewritten as

[TABLE]

and

[TABLE]

While the equivalence of Eqs. 5-6 with Eqs. 7-8 is trivial to see for molecular systems (using integration by parts), showing the equivalence for periodic systems requires a more complex manipulation (see Supplemental Material).

In order to evaluate the double integrals in Eq. 7 and Eq. 8, we note that the kernel of extended interactions is the Green’s function of $-\frac{1}{4\pi}\nabla^{2}$, and we take recourse to the solution of the Poisson equation. Thus, we obtain,

[TABLE]

and

[TABLE]

where $\nabla^{2}\Lambda^{jj,D}_{b}(\textbf{r})=-4\pi\frac{\partial(\phi_{j}(\textbf{r})\phi^{*}_{j}(\textbf{r}))}{\partial r_{b}}$ and $\nabla^{2}\Lambda^{ij,E}_{b}(\textbf{r})=-4\pi\frac{\partial(\phi^{*}_{i}(\textbf{r})\phi_{j}(\textbf{r}))}{\partial r_{b}}$. Thus, finally, the $D$-tensor can be expressed as

[TABLE]

The computationally expensive part of the D-tensor calculation involves the solution of Poisson problems, which are solved on the same FE mesh that represents the KS wavefunctions. However, this computation is embarrassingly parallel over the pairs of orbitals $\phi_{i}$ and $\phi_{j}$. We note that, unlike the $A$-tensor, the dipole-dipole integral entering the $D$-tensor expression (Eq. 4) does not explicitly depend on the nuclear coordinates, and thus we expect the $D$-tensor to be less sensitive to the cusps in the spin density at the nuclei. Therefore, a p-refinement is ideal to systematically improve the accuracy in the calculation of $D$.

II.3 Electric Field Gradient Tensor

The nuclear quadrupole interaction $P$-tensor is directly related to the electric field gradient (EFG) V-tensor. We denote the nuclear quadrupole moment by $Q$ and the quantum number ($a$ component) of the nuclear spin as $I$ ($I_{a}$); the nuclear quadrupole Hamiltonian is given by Schweiger and Jeschke (2001)

[TABLE]

where the EFG $V$-tensor is the second derivative of the electrostatic potential at the nucleus:

[TABLE]

Here $n$ is the electron density (defined as positive), and $Z_{I}$ and $\bm{r}_{I}$ are the charge and position of the $I^{\text{th}}$ nucleus in the system, respectively.

Calculation of the nuclear contribution to the $V$-tensor (second term in Eq. 13) is trivial, and only requires the knowledge of the nuclear charges and the respective positions of the nuclei. We note that the electronic contribution to the $V$-tensor is given by the Hessian of the electrostatic potential. To this end, from a converged self-consistent DFT calculation, we extract the Hartree potential and compute the Hessian at the FE quadrature points. By construction, every nucleus is on an FE node in the FE mesh. Thus, the value of the Hessian at each nucleus is obtained via a projection of the quadrature point values to nodal value. As the $V$-tensor involves point-wise second-order derivatives, a careful convergence study of both h and p refinement is required.

III Results and discussions

We carried out calculations of spin Hamiltonian parameters for a series of molecules/radicals and the nitrogen-vacancy (NV) center in diamond. For the calculation of the NV center, the -1 charge state was considered, which is the most relevant charge state for NV-based quantum information processing. A 64-atom supercell of diamond and $\Gamma$-point sampling of the Brillouin zone were used. In the following discussion of $A$ and $V$-tensors of the NV center, we focus on the nitrogen atom and the three carbon atoms with dangling bonds (DB). All calculations were performed with the PBE functional Perdew et al. (1996). When treating charged systems we included a neutralizing jellium background. All structures were optimized with plane-wave DFT using the QUANTUM ESPRESSO code Giannozzi et al. (2009) and the same structures were used for all-electron calculations.

All-electron FE calculations were performed with the DFT-FE code dft using adaptive real-space meshes. The tensor elements were converged with respect to the FE basis through $h$ and $p$ refinements, within 1-2 MHz for the $A$-tensor, $5\times 10^{-4}\,\text{cm}^{-1}$ for the $D$-tensor and 0.05 a.u. for the $V$-tensor of molecules. Convergence of the spin Hamiltonian parameters for the NV center with respect to the FE basis is presented later in the discussion.

In order to verify our FE results, we also performed PW-based calculations for all systems and GTO-based calculations for molecules. PW calculations of the $A$ and $V$ tensors were carried out with the GIPAW code using the GIPAW pseudopotentials (PP) Cer . PW calculations of the $D$-tensor were conducted with two different PP: GIPAW and ONCV Schlipf and Gygi (2015). We followed the numerical method in Ref.  Rayson and Briddon, 2008 to evaluate Eq. 4 in reciprocal space, using normalized pseudo-wavefunctions Ivády et al. (2014); Falk et al. (2014); Seo et al. (2017); Whiteley et al. (without PAW reconstructions) from the QUANTUM ESPRESSO code. A kinetic energy cutoff of 200 Ry was used for PW calculations of molecules; for the NV center we used 100 Ry for computational efficiency. GTO calculations of $A$, $D$ and $P$ tensors were carried out with the ORCA code Neese (2012). Two Gaussian basis sets were considered: EPR-III Rega et al. (1996) and IGLO-III Kutzelnigg et al. (1990), both of which are designed for an accurate representation of core electrons. We also tested a series of general-purpose basis sets from Dunning and co-workers (cc-pVDZ, cc-pVTZ, cc-pVQZ and cc-pV5Z) Dunning Jr (1989), but we found that the values of $A$, $D$ and $P$ tensors converge poorly as a function of the basis set and the poor convergence prevented any meaningful extrapolation to the complete basis set limit. We present GTO results obtained with cc- basis sets in the Supplemental Material.

Table-I and Table-II show the isotropic (Fermi contact) and spin dipolar component of the $A$-tensor for several molecules (CN, BO, AlO, NH) and the NV center. Due to the symmetry of the systems considered here, the dipolar $A$-tensor has only one independent component (except for DB carbons in the NV center). Denoting the principal values of the dipolar $A$-tensor as $A^{\text{dip}}_{11},A^{\text{dip}}_{22},A^{\text{dip}}_{33}$ ($|A^{\text{dip}}_{11}|=|A^{\text{dip}}_{22}|=\frac{1}{2}|A^{\text{dip}}_{33}|$), we show $A^{\text{dip}}_{33}$ in Table-II. In PW calculations we tested three different treatments of core relaxation (Slater exchange-only, exchange-only and exchange-correlation) implemented in the GIPAW code Bahramy et al. (2007). Experimental values are also shown in the Tables for reference. We note that all of the results presented here, in addition to numerical errors which are quantified and discussed in detail below, suffer from systematic errors introduced by the use of a specific, approximate exchange-correlation functional, the PBE functional. A previous study has shown that more advanced functionals, such as certain hybrid and meta-GGA functionals, may improve the agreement with experiments, relative to GGA functionals, for the $A$-tensor of small radicals and transition metal complexes Kossmann et al. (2007). However, there is yet no consensus on which functional is the most accurate one, in general, for the calculation of the $A$-tensor or other spin Hamiltonian parameters.

We found that in general, GTO results obtained with EPR-III and IGLO-III basis sets are similar, with a mean absolute deviation (MAD) of 3.2 (0.6) MHz for $A^{\text{iso}}$ ($A^{\text{dip}}$) for the systems considered here. FE and GTO results agree well: the MAD between FE and GTO@EPR-III results is 2.5 (1.5) MHz for $A^{\text{iso}}$ ($A^{\text{dip}}$). However, for the Al atom in AlO, FE and GTO@IGLO-III yield different values of $A^{\text{iso}}$ by 56 MHz (9%). We expect the difference to originate from inaccuracies of the IGLO-III basis set used in GTO calculation; for example, we found that GTO calculations using different cc- basis sets yield large variations, between 580 to 520 MHz, for the $A^{\text{iso}}$ value of Al (see Supplemental Materials). Overall, the agreement between FE and GTO results serves as a verification of our FE implementation for the calculation of the $A$-tensor. We note that EPR-III and IGLO-III sets are specialized GTO basis designed for spin Hamiltonian parameter calculations, and they are not available for all elements (for instance, an EPR-III basis set for Al is not available). FE-based calculations, on the other hand, can be performed for any element in the periodic table and the results can be systematically converged with respect to the basis set.

We found that PW calculations agree well with all-electron FE and GTO calculations for $A^{\text{dip}}$, while they deviate slightly for $A^{\text{iso}}$. For $A^{\text{dip}}$, the MAD between FE and PW results is 2.7 MHz, while the MAD for $A^{\text{iso}}$ ranges from 13-17 MHz depending on the treatment of core relaxation in PW calculations. Notably, in the case of the AlO molecule, PW calculations predicted a different sign for the $A^{\text{iso}}$ of the O atom compared to all-electron FE and GTO calculations.

PW and FE calculations for the NV center yielded qualitatively similar values for $A^{\text{iso}}$ and $A^{\text{dip}}$ for both nitrogen and DB carbons. The larger value of $A^{\text{iso}}$ compared to $A^{\text{dip}}$ for the nitrogen atom reveals a strong s character of the spin density on the nitrogen. The spin density on the DB carbons has instead a significant p-type contribution as revealed by the comparable values of $A^{\text{iso}}$ and $A^{\text{dip}}$. There is a sizable difference between DFT results and experimental values for the $A^{\text{iso}}$ of DB carbons (30%), which might be due to the use of a small (64-atom) supercell for the NV center.

In Table-III we present the computed zero field splitting $D$-tensor for several spin-triplet molecules/radicals (O2, CH2, NH, C5H5+) as well as for the NV center. We report the scalar parameter $D=\frac{3}{2}D_{33}$, where $D_{11},D_{22},D_{33}$ are principal values of the $D$-tensor such that $|D_{11}|\leq|D_{22}|\leq|D_{33}|$. For low symmetry systems such as the CH2 carbene, we additionally report the scalar parameter $E=\frac{1}{2}(D_{11}-D_{22})$.

Overall, GTO results show a weak dependence on the basis set, with a MAD of 0.016 $\text{cm}^{-1}$ between values obtained with EPR-III and IGLO-III basis sets. PW calculations show a weak dependence on the chosen pseudopotential, with a MAD of 0.017 $\text{cm}^{-1}$ between ONCV and GIPAW results. Similar to the case of the $A$-tensor, GTO and FE results agree well, with a MAD of 0.001 $\text{cm}^{-1}$ between FE and GTO@EPR-III values. Due to the use of pseudo-wavefunctions for the evaluation of Eq. 4 and the lack of PAW reconstruction, PW results deviate from all-electron ones, with a MAD of 0.064 $\text{cm}^{-1}$ between FE and PW@GIPAW values. For the case of the NV center, results from FE, PW and experiments appear to be in good agreement.

Table-IV summarizes the electric field gradient $V$-tensor for several closed-shell molecules (HCN, NCCN, N2, H2O) and for the NV center. Following the convention of the NQR spectroscopy literature, we report the quadrupole coupling constants $eQV_{33}$, where $V_{11},V_{22},V_{33}$ are principal values of the $V$-tensor such that $|V_{11}|\leq|V_{22}|\leq|V_{33}|$. When computing the $eQV_{33}$ values in Table-IV we considered isotopes with non-zero nuclear quadrupole moment $Q$, and the values of $Q$ are obtained from Ref. Web, . For low symmetry systems, we additionally report $\eta=\left|\frac{V_{22}-V_{11}}{V_{33}}\right|$.

Unlike the $A$ tensor, which depends on charge density differences, the $V$-tensor depends on the absolute value of the charge density and thus it is more sensitive to the details of the electronic structure. Differences are indeed observed for GTO calculations with different basis sets (MAD = 0.09 MHz), as well as between GTO and FE calculations (MAD = 0.18 MHz between FE and GTO@EPR-III). PW results significantly deviate from all-electron GTO and FE results, with a MAD of 0.76 MHz between FE and PW values. In the case of the NV center, PW and FE yield similar nuclear quadrupole coupling for nitrogen, in qualitative agreement with experiment, while for DB carbons the predicted $V_{33}$ values using PW and FE have opposite signs.

Finally, to demonstrate the convergence of the FE results with respect to the basis set, in Table-V we show A-, D- and V-tensors for the NV center computed with different FE polynomial degrees. We denote calculations with $n^{\text{th}}$-order polynomials as FE$n$. For the A-tensor and V-tensor calculations a mesh size of 0.1 Bohr was used surrounding the nuclei, while for the D-tensor calculation the mesh size was 0.5 Bohr. We see in Table-V that our results for the A-tensor are well converged at the FE6 level, as indicated by the small difference (less than $3\%$) between FE5 and FE6 results. Similarly, D-tensor values are well converged at the FE5 level. The numerical value of the $V$ tensor is sensitive to the details of the electronic wavefunctions around the nuclei, as mentioned previously, and its convergence is indeed more challenging compared to that of the A- and D-tensors. We performed V-tensor calculations with polynomial degrees up to 7. At the FE7 level, most of the computed V tensor elements are converged within $10\%$, based on asymptotic estimates obtained by power law extrapolations.

IV Conclusions

In this work, we presented an approach to compute spin Hamiltonian parameters based on DFT, which uses all-electron calculations and finite element basis sets to solve the Kohn-Sham equations. The approach can be applied to both solids and molecules and offers the important advantage of straightforward convergence of the calculations with respect to the basis set, which can be systematically achieved by refinement of the finite element basis.

We reported calculations of the Fermi contact and dipolar component of the $A$-tensor, the spin-spin component of the $D$-tensor and the nuclear quadrupole $P$-tensor for several molecules and for the NV center in diamond. We presented detailed comparisons of results obtained using FE, GTO and PW basis. For molecules, we showed that all-electron results obtained with FE basis sets are in good agreement with those obtained with GTO basis sets.

The approach introduced in our work represents the first step towards building a robust protocol for the first-principles prediction of various spin Hamiltonian parameters based on finite element density functional theory. There are multiple prospects of future work in this direction, both in terms of the level of physics and computational efficiency. It is important to extend the current formalism to include relativistic effects since proper treatment of scalar relativistic effects will be crucial for accurate calculations of spin Hamiltonian parameters of heavy elements. The ability to include spin-orbit coupling effects will also allow for the computation of additional spin Hamiltonian parameters, including the $g$-tensor and the spin-orbit component of the $A$ and $D$-tensor. Further, it would be interesting to develop and test more advanced density functionals, such as meta-GGAs and hybrid functionals, and to establish which functional performs better, compared to experiments. With regards to the computational efficiency, the FE basis functions can be enriched using compactly supported precomputed enrichment functions Kanungo and Gavini (2017), which will drastically reduce the computational cost, while providing systematic convergence. Finally, we plan to utilize a combination of all-electron and pseudopotential based calculation under the same framework, where certain atoms of interest are treated at all-electron level and other atoms are treated using pseudopotential approximation, which will enable the computation of spin Hamiltonian parameters in systems involving thousands of atoms.

We thank P. Motamarri for useful discussions on various aspects of the numerical implementation. This work was supported by MICCoM as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division through Argonne National Laboratory, under contract number DE-AC02-06CH11357. V.G. also gratefully acknowledges the support of the Department of Energy, Office of Basic Energy Sciences, through grant number DE-SC0017380, under the auspices of which the computational framework for FE-based all-electron DFT calculations was developed. K.G. acknowledges the computational resources from the University of Michigan through the Flux computing platform. H. M. acknowledges computational resources from the University of Chicago Research Computing Center.

K.G. and H.M. have contributed equally to this work.