Relativistic calculation of indirect NMR spin-spin couplings using the Douglas-Kroll-Hess approximation

We have employed the Douglas-Kroll-Hess approximation to derive the perturbative Hamiltonians involved in the calculation of NMR spin-spin couplings in molecules containing heavy elements. We have applied this two-component quasirelativistic approach using finite perturbation theory in combination with a generalized Kohn-Sham code that includes the spin-orbit interaction self-consistently and works with Hartree-Fock and both pure and hybrid density functionals. We present numerical results for one-bond spin-spin couplings in the series of tetrahydrides CH(4), SiH(4), GeH(4), and SnH(4). Our two-component Hartree-Fock results are in good agreement with four-component Dirac-Hartree-Fock calculations, although a density-functional treatment better reproduces the available experimental data.


I. INTRODUCTION
Nuclear magnetic resonance ͑NMR͒ spectroscopy has been a powerful tool in chemistry since its discovery in 1946. 1,2In isotropic phase, the NMR spectrum is determined mainly by two parameters; the nuclear magnetic shieldings ͑NMSs͒ and the spin-spin couplings ͑SSCs͒.These parameters have been at the center of theoretical research for decades, starting with the pioneer work of Ramsey. 3,4][7] However, this is not the case for NMR parameters involving heavy elements.The main reason for this is that when heavy elements are involved, relativistic effects must be explicitly taken into account.NMR parameters strongly depend on the electronic environment in the vicinity of the nuclei and, thus, they are affected by relativistic effects.[10][11][12][13][14][15] Most of these works are focused on NMSs, while those reporting calculations of SSCs are less abundant.Full fourcomponent calculations were carried out by Visscher et al. in the hydrogen halides HX ͑X = F, Cl, Br, and I͒ 13 and by Enevoldsen et al. in the tetrahydrides series XH 4 ͑X = C, Si, Ge, Sn, and Pb͒ 9 using the Dirac-Hartree-Fock ͑DHF͒ formalism.Calculations based on two-component ͑quasirelativ-istic͒ approximations are more attractive in terms of computational cost.Autschbach and Ziegler employed the zerothorder regular approximation ͑ZORA͒ to calculate SSCs using Slater functions in the framework of density-functional theory ͑DFT͒. 11,16Recently, Filatov and Cremer developed a code for the calculation of SSCs in molecules containing heavy atoms based on the infinite-order regular approximation with modified metric ͑IORAmm͒ and Gaussian basis sets ͑GBS͒. 10n this work we investigate the suitability of the Douglas-Kroll-Hess ͑DKH͒ approach [17][18][19][20] for the calculation of SSCs.2][23] The main attractives of the DKH method are that it is variationally stable and it can be systematically improved by increasing the order of the transformation. 24In particular, the calculation of NMSs was discussed in this context by Baba and Fukui 25 and by Fukuda et al. 26,27 In our derivation of the magnetic operators in the DKH context we follow below steps similar to those of Fukuda et al., 26 although our final expressions for the second-order perturbative Hamiltonians are more compact, simplifying its implementation.The formal theory of electric and magnetic perturbations in the DKH formalism was presented by Dyall. 28However, no previous computational studies on SSCs using the DKH formalism were performed.SSCs critically test the behavior of the wave function near the nuclei, making its calculation very difficult.In this paper, we investigate the performance of the DKH approach within the point-nucleus model taking XH 4 ͑X = C, Si, Ge, and Sn͒ as model compounds.
The organization of this paper is as follows.We first describe in Sec.II the theoretical background related to the standard Ramsey's theory and DKH theory.In Sec.III we outline our implementation and the computational details, and in Sec.IV, we present benchmark results for the XH 4 ͑ X = C, Si, Ge, and Sn͒ series, where accurate four-component Dirac-Hartree-Fock calculations are available. 9

II. THEORY
The indirect NMR spin-spin coupling tensor J͑N , M͒ can be calculated from the second derivative of the electronic energy E͑N , M͒ in the presence of two nuclear magnetic moments N and M ͑atomic units are used throughout this paper͒, where ␥ N and ␥ M are the nuclear gyromagnetic ratios of nuclei N and M, respectively.This energy derivative can be obtained considering the perturbative magnetic Hamiltonians arising from the presence of the nuclear magnetic moments.
In the point-nucleus model, the nuclear magnetic potential A ͑such that the magnetic field B = ٌ ϫ A͒ is given by with r K = r − R K , and c is the speed of light in vacuum.

A. Ramsey's theory
The lowest-order ͑in c −1 ͒ magnetic interaction perturbative Hamiltonian can be obtained by expressing the mechanical momentum in terms of the canonical momentum p =−iٌ as = p + ͑1/c͒A ͑e = −1 a.u.͒ and writing down the kinetic-energy operator of the Pauli Hamiltonian ͑1/2͒ ϫ͑ • p͒ 2 ͑where is the Pauli spin matrix vector͒.In this way it is possible to identify four perturbative Hamiltonians: the Fermi contact ͑FC͒

͑4͒
the paramagnetic spin orbital ͑PSO͒ and the diamagnetic spin orbital ͑DSO͒

͑6͒
The first three hyperfine operators, Eqs.͑3͒-͑5͒, are linear in K ͑also called paramagnetic terms͒, while the last one, Eq. ͑6͒, is linear in K and L ͑diamagnetic term͒.

B. First-order Douglas-Kroll-Hess theory
In the relativistic framework, the magnetic interaction perturbative Hamiltonian is obtained by introducing the mechanical momentum = p + A / c in the Dirac Hamiltonian: where V N represents the nucleus-electron interaction, c is the speed of light, and ␤ and ␣ are the four Dirac matrices.In what follows, we shall consider a general vector potential A, although we are interested in the vector potential of the nuclear magnetic moments.In order to obtain the expression for the operators involved in the calculation of SSCs in the DKH theory, successive DKH transformations are applied to the Dirac Hamiltonian, Eq. ͑7͒.
The first step to decouple the large and small components of H D is to apply the free particle Foldy-Wouthuysen transformation 29 by means of the unitary operator with the following operators diagonal in p space: ͑9͒ where

͑10͒
The unitary transformation U 0 fully decouples the free particle Dirac Hamiltonian.However, in the presence of external fields, the transformed Hamiltonian contains both even ͑i.e., 2 ϫ 2 block-diagonal͒ and odd ͑i.e., off-diagonal͒ terms, where the even terms are given by while the odd terms can be expressed as Restricting the Hamiltonian to the upper 2 ϫ 2 block ͑corresponding to ␤ = +1͒ gives place to the first-order DKH approximation ͑DKH1͒.At this level of approximation, all terms containing V N and A are decoupled.The spin-orbit ͑SO͒ operator appears in E 1 V .The sum E 0 + E 1 V leads to the DKH1 Hamiltonian in the absence of the external magnetic field, while E 1 A is the perturbative Hamiltonian.At this point, it is important to stress that we are interested in the expression for the perturbative Hamiltonian, and therefore only those terms which are at most quadratic in A should be considered.The perturbative Hamiltonian E 1 A is linear in A, and as a consequence, no diamagnetic contribution arises at this level of theory.E 1 A can be splitted into spin-dependent and orbital ͑spin-free͒ contributions using the relation In terms of two-component Pauli matrices, a working equation for the DKH1 magnetic perturbation is obtained: where are the first-order paramagnetic spin-dependent and orbital perturbations, respectively.It is easy to show that in the nonrelativistic limit ͑c → ϱ͒, the kinetic factors K → 1 and R → ͑2c͒ −1 .Therefore, the nonrelativistic limit ͓keeping O͑c −2 ͒ terms͔ of h 1 S corresponds to the FC and SD operators in Ramsey's theory, Eqs.͑3͒ and ͑4͒, while h 1 O maps into the PSO operator, Eq. ͑5͒.
It is worth mentioning that another approach has been suggested by Fukuda et al. to introduce the vector potential A in the DKH approximation. 26This approach consists in replacing the momentum operator p in the field-free DKH expressions by the mechanical momentum p + A and expanding in powers of A / c.However, it was pointed out that up to second order in A / c the same final result is obtained. 26

C. Second-order Douglas-Kroll-Hess theory
In the DKH2 approximation, a second unitary transformation U 1 is introduced: with W satisfying the condition where ͓A , B͔ represents the commutator of A and B. The operator W is defined such that the first-order odd terms in Eq. ͑24͒ cancel, i.e.,

͓W,␤E
Therefore, W is an odd operator and the only remaining odd term in Eq. ͑24͒ is ͓W , E 1 ͔.The even part of H 2 , provided that Eq. ͑25͒ is satisfied, can be expresed as For the present purpose of obtaining a perturbative expansion in A, it is convenient to split W up into two contributions, W A and W V , satisfying the following conditions: 26 The condition given by Eq. ͑27͒ leads to the standard second-order DKH Hamiltonian for the unperturbed system, The second condition, Eq. ͑28͒, yields to operators which are linear and quadratic in the magnetic potential, and which are the perturbative operators.
The reduction from four-component spinors to electronic two-component spinors is carried out by keeping the upper diagonal block of all operators in Eqs.͑30͒ and ͑31͒.As a result, the magnetic perturbative Hamiltonians acting on two-component electronic spinors can be separated into those arising from the DKH1 transformation, Eq. ͑20͒, and those arising from the DKH2 transformation, and where the superscripts ͑1͒ and ͑2͒ refer to the order of the magnetic potential involved in each case, ͓A , B͔ + is the anticommutator of A and B, and the following operators are defined in terms of two-component electronic spinors: and, in p space, and It is interesting to remark that the leading order in c of the second term in the right-hand side of Eq. ͑35͒ is c −2 , while for the first term it is c 0 .Therefore, in the nonrelativistic limit only the first term of O A remains as O͑c 0 ͒, and the O͑c −2 ͒ diamagnetic Ramsey term is recovered in h 2 ͑2͒ .Explicit working expressions for the second-order perturbative operators h 2 ͑1͒ are given in the Appendix.The second-order paramagnetic Hamiltonian h 2 ͑1͒ can be splitted into spindependent and spin-free contributions analogously to h 1 ͑see the Appendix͒.Using the FC operator ͓Eq.͑3͔͒ with a nonrelativistic wave function leads to the evaluation of the wave function at the nuclear positions ͑provided that a pointlike nuclear model is used͒.This is not a problem since the nonrelativistic wave function remains finite at the nuclei.However, the rela-tivistic ͑and quasirelativistic͒ wave function shows a logaritmic divergence at the Coulomb singularities' centers.Therefore, Ramsey's theory cannot be employed in a quasirelativistic framework since the value of the wave function at the nuclei is ill defined.If a Gaussian basis-set expansion is used, this problem will manifest in the lack of basis-set convergence for the FC term.Unlike the ZORA, the regularization of the FC operator in the DKH1 approximation is done only through kinetic factors, while in second order, the terms arising from h 2 ͑1͒ contain products of the nuclear potential V and the vector potential A.

III. IMPLEMENTATION AND COMPUTATIONAL DETAILS
The matrix elements of A K • p and A K ϫ p are obtained using standard integration routines.To evaluate the perturbative Hamiltonians, we proceed in the same way as in the standard DKH transformation.First, the nonrelativistic kinetic-energy matrix ͑in the atomic basis-set expansion͒ is diagonalized to obtain a unitary transformation to a p 2 representation where all the kinetic factors are diagonal, and therefore straightforwardly evaluated.Second, all matrix elements involving the external potentials ͑nuclei-electron interaction and vector potential͒ are transformed to this representation.Third, the matrix elements of the operators in Eqs.͑20͒ and ͑32͒ are evaluated and back transformed to the r representation.The diamagnetic ͑second order in A͒ contribution is kept nonrelativistic.
Relativistic calculations were carried out using the generalized Kohn-Sham ͑GKS͒ and generalized Hartree-Fock ͑GHF͒ approaches.In these schemes, the electronic ground state is described by means of the set of occupied twocomponent spinors.To guarantee the rotational invariance of the total energy in the GKS case, we employed a noncollinear generalization of the spin-density functional energy and potential. 30The first-order response is obtained by including the matrix elements of the perturbative Hamiltonian in the self-consistent-field ͑SCF͒ procedure multiplied by an appropriate constant.In all cases, a tight convergence criterion was imposed.We note in passing that our GKS and GHF codes allow us to handle Hamiltonians with any arbitrary spin dependence.Therefore, only three SCF calculations per perturbed nucleus are needed to obtain the spindependent response, in contrast to six in standard unrestricted codes.Moreover, as spin and orbital perturbed states do not mix at first order, it is possible to obtain the complete linear response ͑spin plus orbital͒ with only three SCF calculations per perturbed nucleus, instead of nine needed in a standard calculation ͑six for spin and three for orbital perturbations͒.However, the lack of analytical linear response 31,32 in our relativistic code is subject to all the wellknown shortcomings of numerical differentiation.
Nonrelativistic calculations were carried out using the coupled-perturbed Kohn-Sham ͑CPKS͒ scheme or coupledperturbed Hartree-Fock ͑CPHF͒.We have checked the accuracy of our numerical differentiation procedure by using the nonrelativistic Hamiltonian in the GHF code and comparing it with the CPHF results.
An interesting feature of our approach is that the SO operator is variationally included in the Hamiltonian.Thus, we can choose to include SO not as a perturbation but fully self-consistently.It is also possible to "turn SO off" to analyze its effect on the calculated SSCs.
Relativistic effects are included in the one-electron Hamiltonian, while the electron-electron Coulomb, the portion of Hartree-Fock exchange, and DFT exchangecorrelation potential remain nonrelativistic.

IV. BENCHMARK RESULTS
All calculations were carried out using modified versions of the GAUSSIAN ͑Ref.33͒ suite of codes.SSC values reported in this paper are referred to the isotopic species with nonzero nuclear magnetic moment and the largest natural abundance, i.e., 1 H, 13 C, 29 Si, 73 Ge, and 119 Sn.In all calculations, nuclear g factors are g͑ 1 H͒ = 5.585 692, g͑ 13  It is known that a good description of the core electrons is needed to achieve reliable results for calculated SSCs, mainly due to the presence of the Fermi contact operator.In the case of the DKH operators, there is an additional basisset dependence since the transformation to p 2 representation is carried out using the same basis-set expansion.Therefore, the election of the basis set is critical for our calculations.
An extensive analysis of the basis-set dependence of nonrelativistic SSCs can be found in Refs.34 and 35.We have first analyzed the basis-set requirements for the calculation of SSCs using the DKH method starting from Dunning's aug-cc-pVTZ ͑Ref.36-38͒ in fully uncontracted form for H, C, Si, and Ge.For Ge, we have progressively added s-type functions at the tight end of these basis sets using a geometric progression with the ratio of the two steepest s-type orbitals.For Sn, we have proceeded in the same way FIG. 1. Dependence of the calculated 1 J͑C,H͒ coupling in CH 4 with the number of additional tight functions added to the fully uncontracted cc-pVTZ basis set.

204112-4
but starting from the universal Gaussian basis set developed by Malli et al., 39 which consists in a fully uncontracted ͑22s18p12d͒ basis set.In Figs. 1 and 2 we show the convergence behavior of calculated 1 J͑X ,H͒ in CH 4 and SiH 4 with the number of s-type functions added ͑for a given n, n steep functions were added to both H and X͒.For 1 J͑C,H͒, the nonrelativistic ͑NR͒ value ͑Ramsey's theory͒ saturates for n = 4, while DKH1 shows a slightly better convergence.On the other hand, DKH2 seems to be converged with five tight s-type functions added.A similar trend is observed for the series 1 J͑Si, H͒, 1 J͑Ge, H͒, and 1 J͑Sn, H͒, although the number of s-type functions that needs to be added to achieve a good basis-set convergence is different in each case.In view of these results, we have decided to include five tight s-type functions for 1 J͑C,H͒, five for 1 J͑Si, H͒, four for 1 J͑Ge, H͒, and six for 1 J͑Sn, H͒ for the benchmark calculations reported in this paper.
In Table I we present our calculated 1 J͑X ,H͒ for the series of XH 4 molecules with Hartree-Fock and different levels of the DKH approximation.For comparison, fourcomponent Dirac-Hartree-Fock and IORAmm results taken from the literature are also shown.In all cases, DKH1 gives a too large relativistic correction, while this trend is reverted by using the DKH2 Hamiltonian.A similar behavior was found in Ref. 40  In Table II we present the ͑Hartree-Fock͒ calculated spin-dependent ͑S͒ and orbital ͑O͒ contributions to 1 J͑X ,H͒ for the series of XH 4 hydrides, as well as the total relativistic correction for these couplings ͑⌬J͒.Comparing the magnitude of spin-dependent and orbital contributions, these 1 J͑X ,H͒ couplings are dominated by the spin-dependent term, mainly due to the fact that they include the FC contribution.Moreover, the relativistic correction is dominated by its spin-dependent contribution.With the exception of CH 4 , where the relativistic correction reported using fourcomponent random-phase approximation ͑4-c RPA͒ is about 0.6 Hz, the DKH2 approximation reproduces the relativistic effects for these couplings.For 1 J͑Sn, H͒, the total DKH2 calculated correction is about −698.7 Hz, in good agreement with the 4-c RPA calculation, −719.6 Hz.
It is known that DFT can deliver SSCs with reasonable accuracy. 31,32,41,42Therefore, it is wise to test the performance of DFT in this case.To this end, we have employed three representative functionals: PBE ͑Ref.43͒ ͑based on the generalized-gradient approximation͒, PBEh ͑Refs.44 and 45͒ ͑hybrid PBE, also referred to as PBE0 and PBE1PBE in the literature͒, and the widely used hybrid B3LYP. 46We have chosen not to use functionals from the local-spin density approximation family since they were shown to perform poorly for nonrelativistic SSCs. 40,42Table III summarizes our DFT results using the DKH2 approximation for the calculation of one-bond 1 J͑X ,H͒ SSCs in the series of tetrahydrides XH 4 .The three functionals employed here tend to correct the    overestimation of the magnitude of the coupling given by the Hartree-Fock approximation, and therefore the agreement with experiment is better.However, the trend is not uniform along the series of molecules studied here: the hybrid B3LYP functional seems to work better for GeH 4 and SnH 4 molecules, while PBE and PBEh reproduce fairly well the onebond coupling in CH 4 .Although a comparison of SSCs calculated with the ZORA method for this series of tetrahydrides would be interesting, we could not find such results in the literature.

V. CONCLUSIONS
We have implemented and tested the DKH approximation for the calculation of NMR spin-spin couplings in the series of tetrahydrides CH 4 , SiH 4 , GeH 4 , and SnH 4 .We show that the basis-set dependence is slightly different than that found for Ramsey's theory, probably due to the transformation to p 2 space.It is worth mentioning that the results are stable and consistent as the Gaussian basis set is enlarged.The first-order DKH1 gives unexpectedly large relativistic corrections to 1 J͑X ,H͒ even for light nuclei, while DKH2 considerably improves the agreement with relativistic corrections given by four-component DHF calculations. 9n our approach, the SO operator is variationally included in the calculation.Although for the molecules tested in this work a perturbative treatment of SO could be adequate, our approach is more suitable for other systems with a stronger SO coupling.The results obtained for the XH 4 series are promising and suggest that DKH2 is a worthy al-ternative for the quasirelativistic calculation of SSCs.However, a deeper analysis of its performance for different couplings should be carried out.Work along these lines is being carried out in our group.

APPENDIX: SECOND-ORDER PERTURBATIVE TERMS
Second-order magnetic perturbation operators are defined in Eqs.͑32͒ and ͑33͒,

͑A1͒
where O V and O A are defined in Eqs.͑34͒ and ͑35͒.The first two terms in Eq. ͑A1͒ refer to second-order paramagnetic contributions, while the last term represents a diamagnetic contribution.
Expanding the commutators and using shorthand notation defined in Eqs.͑36͒ and ͑37͒, h 2 ͑1͒ can be expressed in terms of its spin and orbital contributions: where and M = KR was used for brevity.
for hyperfine coupling tensors.The SO effect is negligible for the one-bond couplings of CH 4 , SiH 4 , and GeH 4 , while for 1 J͑Sn, H͒ in SnH 4 the SO corrections are −21.3 and −13.3 Hz for DKH1 and DKH2, respectively.

FIG. 2 .
FIG. 2. Dependence of the calculated 1 J͑Si, H͒ coupling in SiH 4 with the number of additional tight functions added to the fully uncontracted cc-pVTZ basis set.
This work was supported by DOE Grant No. DE-FG03-01ER15137.Two of the authors ͑J.I.M and M.C.R.A͒ acknowledge support from CONICET ͑PEI 6217 and PIP 2140͒ and UBACyT ͑X222͒.One of the authors ͑J.I.M͒ also acknowledges the hospitality of Rice University.

TABLE I .
Calculated one-bond 1 J͑X ,H͒ SSCs for the series of tetrahydrides XH 4 using the Hartree-Fock approximation.For comparison, fourcomponent random-phase approximation ͑4-c RPA͒ and IORAmm calculations are also included ͑values in Hz͒.

TABLE II .
Calculated spin-dependent ͑S͒ and orbital ͑O͒ contributions to 1 J͑X ,H͒ and total relativistic correction ͑⌬J, defined as J rel − J non-rel ͒ of onebond 1 J͑X ,H͒ SSCs for the series of tetrahydrides XH 4 .In all cases the Hartree-Fock approximation has been employed.All values are in Hz.

TABLE III .
Calculated one-bond 1 J͑X ,H͒ SSCs for the series of tetrahydrides XH 4 using different density functionals.All values were obtained using the DKH2-SO approximation ͑in Hz͒.