NMR Shielding Tensors and Chemical Shifts in Scalar-Relativistic

Jan 8, 2019 - An efficient formulation of scalar-relativistic NMR shielding tensors based on (one-electron) spin-free exact two-component theory (X2C)...
3 downloads 0 Views 1MB Size
Subscriber access provided by Iowa State University | Library

Quantum Electronic Structure

NMR Shielding Tensors and Chemical Shifts in ScalarRelativistic Local Exact Two-Component Theory Yannick Joel Franzke, and Florian Weigend J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b01084 • Publication Date (Web): 08 Jan 2019 Downloaded from http://pubs.acs.org on January 12, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

NMR Shielding Tensors and Chemical Shifts in Scalar-Relativistic Local Exact Two-Component Theory Yannick J. Franzke∗,† and Florian Weigend∗,‡ †Institute of Physical Chemistry, Karlsruhe Institute of Technology (KIT), Kaiserstraße 12, 76131 Karlsruhe, Germany ‡Institute of Nanotechnology, Karlsruhe Institute of Technology (KIT), Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany E-mail: [email protected]; [email protected] Abstract An efficient formulation of scalar-relativistic NMR shielding tensors based on (oneelectron) spin-free exact two-component theory (X2C) is presented. It utilizes the diagonal local approximation to the unitary decoupling matrix (DLU), which we recently applied to analytical derivatives [J. Chem. Phys. 148, 104110 (2018)]. This allows for routine calculations of large molecules with heavy atoms. Here, the computation times of the NMR shielding tensors of all nuclei formally scale cubically with the size of the system, while memory demands increase quadratically. Efficiency is demonstrated for heavy-element clusters and organometallic complexes with more than 120 atoms and 3200 contracted basis functions. The accuracy of the DLU scheme is evaluated based on

13 C, 17 O, 29 Si, 73 Ge, 119 Sn, 129 Xe, 183 W

and

207 Pb

NMR shielding constants

and chemical shifts using different basis sets. The finite nucleus model is available throughout.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

INTRODUCTION

Today, nuclear magnetic resonance (NMR) spectroscopy is established as probably the most valuable and common analysis among X-ray structure determination. In contrast to the latter, NMR spectroscopy can be carried out in all phases without the need to obtain a single crystal. While the spectra of organic compounds can be easily interpreted based on the large amount of data collected so far, 1 the spectra of inorganic - especially heavy element containing - systems may require theoretical studies to determine the structure and to interpret the corresponding spectra. It is well known that relativistic effects have to be taken into account in heavy element chemistry, even for a qualitative correct description of the electronic structure. 2–13 The most economic way to treat relativistic effects are effective core potentials (ECPs) or pseudo potentials. 14 However, the NMR shielding tensor depends on the density in the core region and then its calculation for heavy elements requires allelectron relativistic treatments. We further note that ECPs are not gauge-invariant without modifications. 15,16 A prominent relativistic theory currently used in quantum chemical calculations is the four-component (4c) approach. This theory has also been applied to NMR shielding tensors by several authors. 17–28 Four-component theory is regarded to be highly accurate but its application to chemical problems may be restricted due to the high demands on computational resources. In the last few decades, two-component (2c) theories were established as a somewhat cheaper alternative. These decouple the electronic and positronic states. Most of the decoupling schemes are based on a unitary transformation 29 to block diagonalize the (one-electron) Hamiltonian. One of the most prominent examples is the second-order Douglas-Kroll-Hess (DKH2) theory 30–32 and its generalization to higher orders 33–36 up to infinite order. 37–42 NMR shielding constants were presented at the DKH2 level more than ten years ago, 43,44 however, problems concerning gauge-invariance arise due to the evaluation of the terms in the momentum space. This has been resolved just recently. 45 Another well established Hamiltonian is the zeroth-order regular approximation (ZORA), 46–48 which 2

ACS Paragon Plus Environment

Page 2 of 62

Page 3 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

has also been applied to NMR shielding constants. 49 It should be noted that the original ZORA formulation depends on the gauge of the electrical potential. This can be resolved with minor modifications. 50,51 While the chemical shifts within ZORA are often sufficiently accurate, the trace of the shielding tensor is not. 52 We further mention the two-component linear response elimination of small component formalism (LRESC) discussed by Aucar et al. in Ref. 53. The currently most promising approach is the (one-electron) exact two-component theory (X2C), in which the Hamiltonian is constructed in matrix form and not in an operator formalism. 54–62 X2C is closely related to the normalized elimination of the small component (NESC), 63–67 however, the decoupling is carried out in one single step. Analytical derivative theory for the calculation of properties then can be applied in rather straightforward manner, 68 which is a major advantage compared to DKH. From a conceptual point of view infinite-order DKH and X2C are superior to ZORA, yet they require the block diagonalization of the parent four-component Dirac matrix in the uncontracted or primitive basis, which can become a bottleneck. This is even more pronounced in case of efficient density functional theory (DFT) or Hartree-Fock (HF) calculations using the resolution of the identity (RI) approximation, 69–73 where the two-electron part is comparably cheap. We have shown that the application of the diagonal local approximation to the unitary transformation (DLU) 74 changes this dramatically. 75,76 In this way calculations with more than 100 heavy atoms on a single CPU with reasonable computational resources become feasible. The effort for the X2C step of Ag147 (one-component variant, polarized double zeta basis) is reduced from more than one day to less than 10 minutes. 75 Similarly, the time for the gradient calculations of Ag13 is reduced from 20 minutes to 1.5 minutes. 76 This allows for applications to mid-size and large molecules of chemical interest. 77 The theoretical formulation of NMR shieldings in X2C has already been derived, 78–82 but routine applications were not possible so far due to the high computational demands. Therefore, we apply the DLU scheme to NMR shielding tensors in this work. We restrict ourselves

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to the scalar-relativistic approach here. This approach is appropriate for organometallic compounds featuring alkynyl compounds, where the main relativistic effect is due to the scalar-relativistic contraction of the s-shells. 83 A scalar-relativistic treatment is sufficient in cases where bonding between a light and a heavy element is dominated by a p-orbital. 84 The effect of spin-orbit coupling on the NMR shielding constants of light-elements has been thoroughly investigated in Ref. 85. An efficient implementation of the two-component variant requires major modifications of our recently published non-relativistic algorithm employing the (multipole-accelerated) resolution of the identity approximation. 86 Thus, we will present that in future works. This paper is organized as follows: In Sec. 2 we give the necessary equations of our implementation, which we describe in Sec. 3. Accuracy and efficiency as well as applications to large chemical systems of the DLU-X2C approach are demonstrated in Sec. 4. Our main results are summarized in Sec. 5, where the conclusions are drawn.

2

THEORY

As previously, 75,76 matrices in the space of spin-free basis functions {λµ } are indicated as M, while matrices in the basis of 2c functions {φµ } are written as M. The 2c spinor basis functions are chosen as direct product of the scalar basis functions with spin functions, {λµ } ⊗ {α, β}. M refers to the corresponding matrices in the space of 4c functions. M refers to the corresponding matrices in the space of 4c functions. A split notation for large (L) and small (S) components is used. Moreover, atomic units and Einstein’s summation convention are used throughout this work except in the DLU equations. We review the Dirac equation in the presence of a magnetic field in 2.1 and apply the X2C transformation in 2.2 and the DLU approximation in 2.3. The calculation of NMR shielding tensors is described in subsection 2.4.

4

ACS Paragon Plus Environment

Page 4 of 62

Page 5 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

2.1

Dirac Equation in the Presence of a Magnetic Field

ˆ is introduced To account for external magnetic fields, the generalized momentum operator ~π according to ˆ = p~ˆ + 1 A. ~ ~π c

(1)

This operator not only includes the linear momentum operator p~ˆ but also the vector potential ~ The vector potential consists of a contribution stemming from the external magnetic field A. ~ A ~ 10 , and from the magnetic moments m ~ 01 : B, ~ I of nucleus I, A O I ~ = A ~ 10 + A ~ 01 = A ~ 10 + A

N X

~ I01 , A

(2)

I

~ O10 A ~ I01 A

1~ ~ O, B × ~rO , ~rO = ~r − R = 2 Z ~I) wI (R ~ ~ dR = −m ~ I × ∇GI , GI = ~ |~r − R|

(3) (4)

~ O and R ~ I denote the position of the gauge origin and the corresponding nucleus, where R ~ acts respectively. wI is a weight function describing the size and shape of the nucleus, and ∇ upon the electronic coordinates. 87 In the point-charge model wI is given by Dirac’s delta distribution, δ(RI ). We also employ the Gaussian charge distribution with the exponential prefactor η taken from Ref. 88,

wI (RI ) =

 η 3/2 π

 exp −η(R − RI )2 .

(5)

Within the Gaussian model for the finite nuclei the vector potential can be further simplified according to Ref. 89 using the lower incomplete gamma function P ( 21 , ηrI2 ), ~ 01 = m ~ I GI = m ~I A ~I ×∇ ~I ×∇ I

5



P ( 21 , ηrI2 ) rI

ACS Paragon Plus Environment

 .

(6)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 62

~ I forms the gradient with respect to the position of nucleus I. This facilitates the Then, ∇ integral evaluation. The point-charge model is obtained by substituting GI = 1/rI or letting η rise to infinity. The one-electron Dirac equation is formulated within a restricted magnetic balance condition 17,18,90 ensuring the correct non-relativistic limit and numerical stability. Thus, only the vector potential of the external field is used in the basis set expansion of the small component. This is done with gauge-including atomic orbitals 91,92 (GIAOs) to account for the gauge-origin invariance and reads

L ψi = cLµi e−iΛµO |φµ i ,   1 ~ 10 ˆ ~σ · p~ + c Aµ S ψi |φµ i , = cSµi e−iΛµO 2c

(7) (8)

here ~σ denotes the vector of Pauli spin matrices. The application of GIAOs leads to a complex ~ µO = R ~ µ −R ~ O . The basis set expansion where R     1 ~ 10 1 ~ 10 19 −iΛµO −iΛµO ˆ ˆ =e p~ + c Aµ . The of the small components uses the identitiy, p~ + c AO e

phase factor including ΛµO =

1 ~ ~ (RµO ×~r)· B, 2c

one-electron Dirac equation in the presence of a magnetic field reads  V  Π









 CL− CL+   S =  02 ( 4c12 W − T ) CS− CS+ Π







02  CL− CL+  − 02   .  1 C C 0  T S− S+ 2 + 2 2c

(9)

Here, the positive sign indicates the chemically relevant electronic states, while the negative sign refers to the positronic solutions. The electronic solutions obey the relation 2 + c2 > 0, where c2 equals the electronic rest energy in atomic units. The matrix representation of the

6

ACS Paragon Plus Environment

Page 7 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

overlap integrals and the scalar potential Vˆ in the basis of 2c functions are block diagonal 



S 0  iΛ S =   , Sµν = hλµ |e µν |λν i , 0 S   V 0  iΛ ˆ V =   , Vµν = hλµ |e µν V|λ νi . 0 V

(10)

(11)

The generalized momentum matrix and the small-small block of the metric contain the Pauli matrices in ~σ = (σx , σy , σz ) due to the restricted magnetic balance condition. So, these matrices are not block diagonal in the presence of a magnetic field. The matrices read 1 ~ 10  1 ~ 10 1 ~ 01  1 hφµ |eiΛµν ~σ · (p~ˆ + A ~σ · (p~ˆ + A ) |φν i , ν + AI ) 2 c c c ν 1 1 ~ 10 1 ~ 01  1 ~ 10  = hφµ |eiΛµν ~σ · (p~ˆ + A ~σ · (p~ˆ + A + AI ) |φν i , ν ) 2 c c ν c   1 1 10 10 iΛµν ~ ) Vˆ ~σ · (p~ˆ + A ~ ) |φν i , = hφµ |e ~σ · (p~ˆ + A c ν c ν 1 1 ~ 10  1 ~ 10  = hφµ |eiΛµν ~σ · (p~ˆ + A ~σ · (p~ˆ + A ) |φν i . ν ) 2 c c ν

Π †µν =

(12)

Π µν

(13)

Wµν Tµν

(14) (15)

~ ν10 indicates the position of the gauge origin, whereas The subscript of the vector potential A ~ 01 denotes the atom center of the basis function and the corresponding the subscript of A I magnetic moment. To arrive at a scalar-relativistic approximation, the explicit spin-dependent parts are separated from the spin-independent or spin-free ones by applying the Dirac identity. The spin-free part is again block diagonal in the basis of the 2c functions. Thus, the problem can be reduced to an one-component approach. 80

7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.2

Page 8 of 62

Block Diagonalization of the One-Electron Dirac Hamiltonian in One Step

The central idea of the scalar-relativistic X2C ansatz is to obtain the matrix representation ˆ operator in a single step by diagonalization of the spin-free Dirac matrix and using of the X the positive energy solutions CS+ and CL+ X = CS+ (CL+ ).−1

(16)

ˆ operator is used to calculate the unitary transformation for block diagonalizing the The X Dirac Hamiltonian 29 according to Heully et al. 93 





†  1 −X  R 0  U =  , X 1 0 R’

(17)

where the unitary matrix is written as a product of a decoupling and a renormalization operation. Here, R refers to the electronic and R’ to the positronic states, respectively. Similar to gradient calculations the transformation to an orthogonal basis 58,75 is omitted in this work since the diagonalization is not costly compared to solving the first- and secondorder response equations. The electrons-only Hamiltonian is then formulated based on the NESC matrix L 56,60,63 and a renormalization matrix R 61

h = R† LR

(18)

with †



 1 W−T X 4c2

(19)

˜ = S + 1 X† TX. S1/2 , S 2c2

(20)

L=V+X Π+Π X+X





and R=S

−1/2



S

−1/2 ˜

SS

−1/2

−1/2

8

ACS Paragon Plus Environment

Page 9 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

The renormalization matrix allows for straightforward incorporation of the all-electron relativistic one-electron Hamiltonian into the existing non-relativistic code since it transforms the relativistic metric to the non-relativistic one.

2.3

Diagonal Local Approximation to the Unitary Decoupling Matrix

The relativistic decoupling is carried out in the (linear-independent) primitive spherical atomic orbital basis (AO basis). Hence, the dimension of the matrix, which has to be diagonalized, is demanding for larger systems. The diagonalization scales like N 3 , where N measures the size of the system. In low-order-scaling approaches such as DFT, the decoupling will thus dominate the computation times for large systems and restrict the applicability of the approach. Local approximations are therefore introduced to reduce the dimension and to accelerate algebraic operations. Within the framework of atom-centered basis functions, the Dirac matrix is divided into atomic blocks. Matrix elements with the same atom-center in the bra and ket basis functions belong to the atomic diagonal blocks, while the others make up the atomic off-diagonal blocks. The inclusion of the latter is crucial to achieve high accuracy. 94 The diagonal local approximation to the unitary decoupling matrix (DLU) by Peng and Reiher 74 achieves high accuracy and reduces the scaling of the X2C step from N 3 to N 2 . The unitary matrix is approximated as the direct sum of the atomic blocks

U LL =

M

LL UAA =

A

U SL =

M

M

RAA ,

(21)

XAA RAA .

(22)

A SL UAA =

A

M A

The decoupling has to be carried out only for the atomic diagonal blocks, which greatly facilitates the calculation. In terms of Eq. (18), the full decoupling matrix in the NESC matrix is approximated as direct sum of the atomic matrices including the full potential. 9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 62

This is combined with a diagonal local approximation to the renormalization matrix. The atomic diagonal blocks of the X2C-Hamiltonian then read † h+ AA = RAA LAA RAA ,

(23)

which can be also written as

h+ AA

=

R†AA

    1 † † † VAA + ΠAA XAA + XAA ΠAA + XAA WAA − TAA XAA RAA . 4c2

(24)

The atomic off-diagonal blocks are obtained by matrix multiplications according to † h+ AB = RAA LAB RBB

(25)

or

h+ AB

=

R†AA

    1 † † † WAB − TAB − XBB RBB . (26) VAB + ΠAB XBB + XAA ΠAB + XAA 4c2

This scheme was called DLU(All) in Ref. 74. For large molecular systems, the atomic offdiagonal block will become more demanding than the atomic diagonal blocks. However, further approximations to accelerate the evaluation of the off-diagonal blocks like the DLU(NR) or DLU(BP) schemes of Ref. 74 do not come with the desired accuracy. In addition, Peng and Reiher proposed an effective linear scaling technique, DLU(NB). 74 Here, the relativistic corrections of the off-diagonal blocks are only calculated for the nearest neighbors, whereas the non-relativistic approximation is used for the other blocks. However, only a limited amount of data has been provided yet. We carried out exploratory energy and gradient calculations on the organometallic iridum complexes of Ref. 76 with the DLU(NB) scheme based on the van-der-Waals radii. The convergence behavior turned out to be unfavorable and the error increases with the number of atoms. Thus, we discuss only the DLU(All)

10

ACS Paragon Plus Environment

Page 11 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

concept herein. We note that diffuse or augmented basis sets increase the error of the DLU scheme due to the local approximation to the renormalization matrix: The DLU errror of the cohesive energy for Ag13 (158 kJ/mol) is increased from 0.20 kJ/mol to 0.40 kJ/mol when adding a diffuse (2s2p2d) set. 75 In case of structure optimizations a similar picture is revealed. Here, the structure of the astatine molecule of Ref. 76 with an augmented basis set (x2c-SVPall2c 95 with additional diffuse functions taken from aug-cc-pVDZ_PP 96 ) shows a DLU error of 0.3 pm, which is still small, in particular compared to the effect of the augmentation itself, which leads to 3.0 pm longer bonds. Thus, the DLU scheme is well balanced in terms of cost and accuracy.

2.4

Analytical Derivatives and NMR Shielding Tensors within OneElectron Exact Two-Component Theory

I The NMR shielding tensor σαβ (α, β = x, y, z) of nucleus I is defined as the derivative of

the energy with respect to the magnetic field and the magnetic moment

I σβα

 =

∂ 2E ∂Bα ∂mI,β



 = tr P 0

∂ 2h ∂Bα ∂mI,β



 + tr

0

∂P ∂h ∂Bα ∂mI,β

 .

(27)

0

The notation (...)0 indicates the limit Bα = 0, mI,β = 0 and P is the one-particle density matrix. The corresponding derivative, also called perturbed density matrix, is obtained via the coupled-perturbed Hartree-Fock (CPHF) 97,98 or Z-vector 99 equations. To solve these equations, the derivative of the one-electron Hamiltonian with respect to the magnetic field is required to construct the derivative of the Fock matrix. In case of non-hybrid density functionals and the uncoupled approximation, no iteration steps are necessary in DFT calculations. 100,101 Techniques for efficient use of GIAOs have been presented in the groups of Pulay 102 and Ahlrichs. 103

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.4.1

Page 12 of 62

Perturbed Density Contribution

The perturbed density contribution requires the first derivatives of the one-electron Hamiltonian with respect to the magnetic field and with respect to the magnetic moment of the corresponding nucleus. The first one is needed to solve the CPHF equations for the perturbed density matrix. The scalar-realtivistic derivatives of the matrices in Eq. (9) with respect to the magnetic field read 80 

∂Vµν Bα



∂Π†µν Bα

!



=

i ~ µν × ~r)α V|νi ˆ , hµ|(R 2c

(28)

=

i ~ µν × ~r)α pˆ2 |νi + 1 hµ|(~r × p~ˆ)α |νi , hµ|(R 4c 2c

(29)

0

0

∂Wµν i ~ µν × ~r)α p~ˆ · Vˆ p~ˆ|νi + 1 hµ|(~r × p~ˆ)α Vˆ + Vˆ (~r × p~ˆ)α |νi , (30) = hµ|(R Bα 2c 2c 0 !   † ∂Πµν ∂Tµν , (31) = Bα 0 Bα 0   i ∂Sµν ~ µν × ~r)α |νi . = hµ|(R (32) Bα 0 2c

Only the small-large and large-small block of the Dirac matrix shows non-vanishing moment derivatives, ∂Π†µν mI,β

! = 0

  1 ~ I GI ) × p~ˆ |νi . hµ| (∇ 2c β

(33)

All first derivatives here are purely imaginary. Thus, complex algebra can be avoided as ~ and m discussed in Ref. 103. In the following, we use λ for the perturbation (B ~ I ) and drop the bracket (...)0 indicating the corresponding limit, 

∂M ∂λ



= M.λ

(34)

0

Note that the unperturbed quantities are evaluated in the same limit. This leads to the usual expressions for the spin-free Dirac matrix without the presence of a magnetic field, which are real and symmetric. 12

ACS Paragon Plus Environment

Page 13 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

To obtain the derivative of Eq. (18), the derivative of the NESC matrix is needed. This reads λ

λ

†,λ



λ

†,λ



λ

†,λ



L = V +Π X+Π X +X Π+X Π +X     1 λ 1 † λ † +X W −T X+X W − T Xλ . 42 42

 1 W−T X 42

(35)

The crucial part is the derivative of the decoupling matrix Xλ since this requires the perturbed coefficients of the large and small component, which are obtained by solving first-order response equations. Ways to calculate the derivative of the decoupling matrix are discussed in the appendix of this work (see Eqs. (57) and (58)). The derivative of the renormalization matrix Rλ can be calculated based on the Sylvester matrix equation as outlined by Cheng and Gauss, 104 and Zou et al.. 105 Here we used the generalized eigenvalue decomposition method of the latter authors. We have implemented both methods and found that the generalized eigenvalue decomposition method is faster and also easier to apply to second derivatives. Application of the DLU scheme to first derivatives was presented in Ref. 76. The diagonal blocks of the Hamiltonian are based on the algorithm above with † † λ λ hλAA = R†,λ AA LAA RAA + RAA LAA RAA + RAA LAA RAA .

(36)

The off-diagonal blocks are obtained by matrix multiplications according to † † λ λ hλAB = R†,λ AA LAB RBB + RAA LAB RBB + RAA LAB RBB ,

13

ACS Paragon Plus Environment

(37)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 62

with † †,λ † λ λ LλAB = VλAB + Π†,λ AB XBB + ΠAB XBB + XAA ΠAB + XAA ΠAB +     1 1 λ † †,λ λ + XAA WAB − TAB XBB + XAA W − TAB XBB 42 42 AB   1 † + XAA WAB − TAB XλBB . 42

(38)

For the evaluation of the perturbed Hamiltonian, only the atomic diagonal elements of the decoupling and the renormalization matrix are required (see Eqs. (52) and (58)). For this, the atomic diagonal blocks of the Hamiltonian have to be diagonalized. This is done in the linear-independent atomic orbital basis, which is obtained for each block by transforming the integrals to the eigenbasis of the overlap matrix. Here, the eigenvectors associated with eigenvalues less than 10−14 are removed. The perturbed overlap and scalar potential matrix do not contribute to the derivative of the atomic diagonal block with respect to the magnetic field since these matrices depend on the distance, Rµν , which then vanishes. 2.4.2

Unperturbed Density Contribution

The unperturbed density contribution involves the second derivative of the Hamiltonian. Only the small-large and large-small block of the Dirac matrix contain non-zero integrals ∂ 2 Π†µν ∂Bα mI,β

! = 0

  1 ~ ~ rν )β |νi hµ|δ ( ∇G ) · ~ r I ν − (∇GI )α (~ αβ 4c2   i ~ µν × ~r)α (∇G ~ I ) × p~ˆ |νi . + 2 hµ|(R 4c β

(39)

Note that these integrals are real. In non-relativistic quantum chemistry the first term is the diamagnetic contribution, whereas the second is denoted paramagnetic unperturbed density contribution to the NMR shielding tensor. Again, we drop the notation (...)0 for the zero

14

ACS Paragon Plus Environment

Page 15 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

limit in the following and use the notation 

∂ 2h ∂Bα ∂mI,β



= h.Bα ,mI,β

(40)

0

Differentiating the X2C-Hamiltonian in Eq. (18) yields hBα ,mI,β =

  R†,Bα ,mI,β L R + R† L RBα ,mI,β + R†,Bα L RmI,β + R†,mI,β L RBα   (41) + R†,Bα LmI,β R + R† LmI,β RBα + R†,mI,β LBα R + R† LBα RmI,β  + R† LBα ,mI,β R .

The second derivative of the NESC matrix reads LBα ,mI,β = Π†,Bα ,mI,β X + Π†,Bα XmI,β + Π†,mI,β XBα + Π† XBα ,mI,β + X†,Bα ,mI,β Π + X†,Bα ΠmI,β + X†,mI,β ΠBα + X† ΠBα ,mI,β  + X†,Bα ,mI,β (W − T) X + X†,mI,β WBα − TBα X

(42)

+ X†,Bα (W − T) XmI,β + X†,mI,β (W − T) XBα  + X† WBα − TBα XmI,β + X† (W − T) XBα ,mI,β . The reader is referred to the appendix, where our calculations of the second derivatives of the decoupling matrix based on second-order response equations is described (see Eqs. (60) and (61)). To calculate the corresponding derivatives of the renormalization matrix, the Sylvester equation  −1 Bα ,mI,β ˜ S RBα ,mI,β R + R RBα mI,β = S − RBα RmI,β − RmI,β RBα

(43)

is solved analogously to that for the first derivatives. 81,105 The above equations are used to get the second derivatives of the atomic diagonal blocks

15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 62

in the DLU procedure. The off-diagonal blocks are then obtained as B ,mI,β

α hAB

†,B ,mI,β

= RAAα

m

m

†,Bα I,β I,β α LAB RBB + R†,B AA LAB RBB + RAA LAB RBB

†,m

B ,mI,β

†,m

m

α † α + RAA I,β LB AB RBB + RAA LAB

m

I,β α RBB + R†AA LB AB RBB

B ,mI,β

α † † I,β Bα α + RAA I,β LAB RB BB + RAA LAB RBB + RAA LAB RBB

(44) .

Again, all the equations can be evaluated without complex algebra utilizing ordinary matrix multiplications. The second derivatives involve terms with first derivatives of the magnetic field and the magnetic moments. This makes the separation of paramagnetic and diamagnetic contributions as done in the non-relativistic limit problematic. Thus, we only distinguish between the perturbed and the unperturbed density contribution.

3

IMPLEMENTATION

We have implemented the scalar-relativisitc DLU-X2C NMR shielding tensors in the MPSHIFT module 86,103,106–108 of TURBOMOLE. 109–111 The necessary additional integrals for the small-small block and the finite nucleus model are implemented using Gauss-Hermite and Gauss-Rys quadrature. 112–114 We have also reimplemented parts of the other one-electron integrals to increase the efficiency of the code. The finite nucleus model based on a Gaussian ~ and charge distribution 88 is available for the scalar potential V and the vector potential A the corresponding derivatives based on the ideas of Hennum et al. in Ref. 89. Linear algebra operations are carried out with the corresponding LAPACK routines. 115 An OpenMP version is available. The X2C energies 75 and gradients 76 were also extended to an OpenMP implementation. The integrals are evaluated in the primitive Cartesian atomic orbital (pCAO) basis with screening criteria applied and then being transformed to the primitive spherical atomic orbital (pAO) basis. The atomic blocks are constructed in the pAO basis. To facilitate the diagonalization and the actual X2C steps, the atomic diagonal blocks are transformed to a 16

ACS Paragon Plus Environment

Page 17 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

linear-independent pAO basis by diagonalizing the overlap matrix. The calculation of the atomic diagonal and off-diagonal blocks follows the scheme outlined in Ref. 76. The calculation of the unperturbed density contribution is more demanding than the perturbed density contribution due to the second-order response equations, as the solution of the second-order response equations (Eqs. (60) and (61) in the appendix) requires a higher number of matrix multiplications. Both formally scale like N 3 when calculating the NMR shielding tensors of all atoms, N being a measure for the size of the system. The number of atomic diagonal blocks scales linearly, while that of the off-diagonal blocks scales quadratically. The derivative of the Hamiltonian with respect to the magnetic moments of all nuclei yields another factor of N , so the above scaling is established. In contrast, the full (non-DLU) X2C algorithm scales like N 4 due to the N 3 scaling algebraic operations. The effort can be reduced by selecting only certain atoms whose NMR shielding tensors have to be calculated. If only a single nucleus is of interest, this results in an N 2 (DLU-X2C) and N 3 (X2C) scaling behavior. This is the same scaling as for the one-electron Hamiltonian in energy calculations. The derivative of the Hamiltonian with respect to the magnetic field scales quadratically (DLU-X2C) and cubically (X2C), respectively. Memory demands increase quadratically and no quantities have to be stored on disk. The reduced memory requirements due to the atomic blocks in the DLU scheme are crucial, since otherwise the application of relativistic all-electron theory to large molecules may be restricted by memory demands. In case of symmetry, only the NMR shielding tensors of symmetry-non-redundant nuclei are evaluated and symmetry-redundant integrals are skipped. However, symmetry is presently not exploited in the X2C response equations. Due to the magnetic field, matrix elements between different irreducible representations have to be calculated. This may lead to rectangular matrices. Thus, symmetry is only fully exploited after the relativistic decoupling, as described in Ref 107. Contributions from the conductor-like screening model (COSMO) 116 and additional point

17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 62

or Gaussian smeared charges are calculated without picture-change correction.

4

Assessment of Accuracy and Efficiency

The gauge-origin and translation invariance of the DLU-X2C Hamiltonian for NMR shifts is demonstrated in subsection 4.1. The accuracy of the local approximation is assessed with different basis sets and various molecules with 207

13

C,

17

O,

29

Si,

73

Ge,

119

Sn,

129

Xe,

183

W and

Pb NMR shifts in 4.2. Therein, the DLU error is first determined for organometallic com-

pounds of the fourth main-group elements. We discuss the chemical shifts of xenon fluroides at various levels of theory complementing the data of Ref. 80. The all-electron relativistic ansatz is compared to ECPs for

17

O NMR shifts of transition-metal oxo compounds. In

subsection 4.3, the efficiency is shown for large metal clusters to estimate the scaling of the NMR calculations with the system size. Organometallic complexes are briefly considered to show the practical applicability of the implementation to systems with more than 100 atoms. In this section, the NMR tensors of all atoms of the molecule are calculated. As discussed in the last section, the scaling can be reduced by considering only the nuclei of interest and their NMR shifts.

4.1

Demonstration of Gauge-Origin and Translational Invariance

Due to the employed GIAOs, the working equations do not depend on the gauge-origin. However, a reference to the origin of the coordinate system is present. Therefore, the perturbed and unperturbed density contribution of a nucleus depends on its point in space and may change under translation just as in non-relativistic GIAO calculations. Nevertheless, the isotropic shielding constant needs to be invariant under translation. Gauge-origin and translation invariance were studied for a set of eight molecules. The structures of HF, HI, WO, NH3 and BiH3 were taken from Ref. 117 and the structures of the xenon fluorides XeF2 , XeF4 and XeF6 from Ref. 80 (X2C-CCSD(T)). The molecules were translated along

18

ACS Paragon Plus Environment

Page 19 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

the z axis by 120 bohr. The xenon fluorides were additionally rotated about the z axis by 60 degree. The HF method was applied together with the uncontracted Dyall triple zeta basis plus 5s5p valence correlating functions 118–122 for I, Xe, W, Bi and an uncontracted cc-pVTZ basis 123 for H, O, F. Following the suggestion of Dyall, 118 the most diffuse primitive set of d functions from the SCF set of the basis was replaced with the 2d correlating functions in cases where these are similar (here I, Xe). The finite nucleus model was selected and a SCF convergence threshold of 10−8 Eh was chosen. The results are given in Tab. 1. The largest change under translation is found for xenon tetrafluoride. Here the isotropic shielding constant changes by 0.14 ppm. The effect is still very small considering the absolute value of 485.8 ppm and the size of the basis set. The variations under translation and/or rotation are very similar for the full X2C Hamiltonian, thus the DLU approximation does not introduce any error concerning gauge-origin and translation invariance. Table 1: Isotropic NMR shielding constants in ppm for various molecules using the DLUX2C/HF method with Dyall-VTZ/cc-pVTZ(unc) basis. LA refers to the lighter atom and HA to the heavier atom, respectively. ∆ Isotropic denotes the difference of the isotropic shielding constant when translating the molecule by 120 a.u. along the z axis. Additionally, the xenon compounds were rotated about the z axis by 60 degree. LA Isotropic HF HI WO NH3 BiH3 XeF2 XeF4 XeF6

27.8 30.4 −372.6 31.3 27.4 455.2 272.3 158.5

∆ LA Isotropic HA Isotropic ∆ HA Isotropic 9 × 10−5 −3 × 10−5 3 × 10−4 3 × 10−6 −7 × 10−5 −5 × 10−6 −6 × 10−3 3 × 10−3

19

409.9 4489.2 36 282.2 261.9 7131.4 1982.9 485.8 1009.9

ACS Paragon Plus Environment

−8 × 10−4 6 × 10−3 2 × 10−3 8 × 10−4 8 × 10−2 −1.0 × 10−4 1.4 × 10−1 −9 × 10−3

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Mean absolute error (MAE) of the DLU-X2C Hamiltonian in 13 C, 29 Si, 73 Ge, 119 Sn and 207 Pb NMR shits for a small test set of organometallic compounds. Max. Error refers to the maximal absolute error while range denotes the span of the calculated shifts. Errors and range of the NMR spectra are given in ppm. MAE Max. Error Range 13

C Si 73 Ge 119 Sn 207 Pb 29

4.2

0.01 0.06 0.07 0.13 0.17

0.02 0.12 0.09 0.19 0.2

109 86 56 301 601

Accuracy of the Local Approximation and Further Studies

4.2.1

Organometallic Main-Group Compounds

First, 13 C, 29 Si, 73 Ge, 119 Sn and 207 Pb NMR shifts based on a small test set of organometallic compounds 83 were analyzed to assess the accuracy of the DLU approximation. The shifts were calculated with the reference compounds SiMe4 (for 13 C, 29 Si), GeMe4 (for 73 Ge), SnMe4 (for

119

Sn) and PbMe4 (for

207

Pb). The structures of these were optimized with the spin-

orbit DLU-X2C Hamiltonian 75,76 at the PBE0 124,125 level of theory (grid 4a) using the x2cTZVPPall-2c basis set 95 employing the RI-J approximation and the finite nucleus model. These grids feature an increased number of radial points for all-electron calculations of heavy elements and lead to an improved convergence behavior and accuracy in the self-consistentfield procedure (convergence criterion 10−8 Eh ). Moreover, Grimme’s dispersion correction D3 126 was included. The structures of the test set were taken from Ref. 83, however, the molecules were moved so that the coordinate axes become principial axes of inertia. The NMR calculations were carried out using the PBE0 functional and Dyall’s core-valence triplezeta basis set for Ge, Sn and Pb, 118–121 and the cc-pCVTZ basis for C and Si, 127,128 and the cc-pVTZ basis for H. 123 All basis sets were employed in decontracted form. To reduce linear dependency of the Dyall core-valence triple-zeta basis, the two most diffuse d functions from the SCF set of Ge were removed. All NMR shielding constants and chemical shifts are given in the Supporting Information, we only discuss the statistical errors here. The mean absolute 20

ACS Paragon Plus Environment

Page 20 of 62

Page 21 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

errors of each type of nucleus are presented in Tab. 2. The maximum DLU error is 0.2 ppm, while the spectra cover a range of 86 ppm (29 Si) to 601 ppm (207 Pb). Thus, the error of the DLU approximation is negligible. NMR shielding tensors are strongly affected by the geometrical parameters. The above procedure has assessed the DLU error of NMR shielding tensors itself. Usually the DLU approximation is applied throughout. This means that the structure is also optimized using the DLU-X2C-Hamiltonian. We have recently shown that the DLU error introduced into the structure is very small. 76 Nevertheless, the combined error is discussed in the following sections.

4.2.2 129

Xenon Fluorides

Xe NMR shieldings were studied by Cheng and Gauss in wavefunction based theory. 80

These compounds serve as a further example to asses the accuracy of the DLU approximation, also compared to the size of the impact of relativity and to the effect of different density functionals. The results obtained at the spin-free X2C-CCSD(T) level provide a scalar-relativistic benchmark value for these systems in wavefunction based theory. The deviation between the experimental results and CCSD(T) can mainly be addressed to spin-orbit coupling. 80 Note that the experimental shifts were measured with respect to neat XeOF4 in the gas-phase 129 while theoretical studies use the xenon atom as reference. Both references are common for

129

Xe shifts. The two scales can be converted to each other by the shift

of xenon gas at infinite dilution (−5460 ppm) with respect to neat XeOF4 . 130 This value is obtained from the experimental

129

Xe shift in solution by extrapolation.

In our study, the BP86, 131,132 PBE, 124 KT3, 133 TPSS, 134 B3LYP, 135–137 PBE0 124,125 and the TPSSh 138 functional were selected to cover the most common types of functionals together with fine grids (grid 4a) for numerical integration of the exchange-correlation potential. The uncontracted relativistic atomic natural orbital basis sets, 139 ANO-RCC(unc), were used together with strict SCF thresholds of 10−9 Eh . Structures were taken from Ref.

21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 62

Table 3: 129 Xe NMR chemical shifts of xenon fluorides obtained at various levels of theory. Chemical shifts are given in ppm with respect to the isotropic shielding of the Xe atom. Method 1 uses structures from Ref. 80 at the X2C-CCSD(T)-ANO-RCC(unc) level. The ANO-RCC(unc) basis set was used in the NMR shielding calculations of that method. Method 2 indicates NMR calculations based on the geometry optimized at the corresponding level with the x2c-TZVPall-2c bases employed in the relativistic calculations and the def2-TZVP/TZVPall basis in the non-relativistic ones (NR). EXP refers to the experimental value in the gas phase taken from Ref. 129. The CCSD(T) results were taken from Ref. 80. XeF2 Level

XeF4

XeF6

Method

NR

X2C DLU

NR

X2C DLU

NR

X2C DLU

HF BP86 PBE KT3 TPSS B3LYP PBE0 TPSSh CCSD(T)

1 1 1 1 1 1 1 1 1

3551 3406 3383 3282 3284 3546 3412 3307 3238

4018 3713 3693 3594 3594 3897 3759 3632 3564

4018 3713 3693 3594 3594 3897 3759 3632 –

4901 5616 5577 5385 5412 5577 5360 5344 4982

5504 6185 6147 5950 5979 6147 5931 5911 5509

5504 6185 6147 5950 5979 6176 5931 5911 –

4289 5227 5183 5032 5090 5145 4991 5028 4641

4858 5903 5856 5699 5762 5819 5646 5694 5258

4858 5903 5856 5699 5762 5819 5645 5693 –

HF BP86 PBE KT3 TPSS B3LYP PBE0 TPSSh MP2

2 2 2 2 2 2 2 2 2

3165 3148 3113 2986 3040 3288 3135 3060 2986

3797 3654 3622 3502 3538 3840 3677 3574 3503

3797 3654 3622 3502 3538 3840 3677 3574 3503

4429 5533 5488 5161 5298 5338 5064 5144 4759

5299 6380 6328 5997 6126 6223 5919 5982 5677

5298 6380 6328 5998 6127 6223 5919 5982 5677

4127 5713 5661 5293 5536 5461 5148 5346 4992

4670 6351 6290 5918 6153 6102 5751 5954 5664

4670 6352 6291 5919 6154 6103 5752 5955 5665

EXP

3386

5623

5425

80. Therein, the octahedral structure was employed for XeF6 . The finite nucleus model was employed in the relativistic calculations here. The results are given in Tab. 3. The shielding constants are given in the Supporting Information of this work. The DLU approximation shows a mean absolute error with respect to full X2C of 0.18 ppm, while the range of the NMR shifts of the three molecules covers more than 2000 ppm and the scalar-relativistic effect leads to an increase by ca. 300 ppm for XeF2 , 550 ppm for XeF4 and 700 ppm for XeF6 . The error of the DLU scheme is thus negligible. The experimental

22

ACS Paragon Plus Environment

Page 23 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

shifts increase from XeF2 to XeF6 by about 2000 ppm and to XeF4 by about 1800 ppm. The non-relativistic Hamiltonian does not accurately reproduce this trend. In particular the increase from XeF2 to XeF6 is significantly underestimated by roughly 300 ppm. This holds for all quantum chemical methods used herein. Taking scalar-relativistic effects into account greatly improves the agreement with the experimental findings. Pure and hybrid density functionals yield very reasonable results. The trend is well reproduced by TPSSh. All functionals overestimate the chemical shifts compared to the CCSD(T) and the experimental results. It should be noted that CCSD(T) does not accurately describe the trend of the

129

Xe shifts. The calculated shifts increase by 1700 ppm from XeF2 to XeF6 and by

1250 ppm from XeF2 to XeF4 . Despite the good performance of DFT, it should be noted that DFT calculations are usually carried out using a double or triple zeta basis and the corresponding structure optimization where error cancellation might even lead to a better agreement with the experimental results. In addition, we carried out structure optimizations and relativistic NMR calculations with the x2c-TZVPall-2c bases. 95 In the non-relativistic calculations the def2-TZVP basis 140 was used for fluorine and the TZVPall 141 for xenon. Here, the corresponding symmetry constraints for XeF2 (D∞h ), XeF4 (D4h ) and XeF6 (Oh ) were employed. Furthermore, MP2 calculations were carried out. All structures are given and compared to the experimental results 142,143 in the Supporting Information. The structure of XeF6 has been ambitiously discussed in the literature, see Ref. 144 for a good overview. Both the octahedral and the C3v structure are minima of the potential energy surface. The latter is likely to be the global minimum in the gas phase. 144 However, the energy difference is very small and the results depend on relativistic, electron correlation and basis set effects. 145–149 For instance, we found that at the HF level the octahedral structure is a saddle point according to numerical frequency calculations, this has also been observed before. 148 Stretching along the corresponding normal mode leads to a C2v transition state which is lower in energy. The shallow potential energy surfaces requires a more sophisticated treatment than HF. The

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

structures obtained at the DFT levels are minima. In the NMR experiments only one type of fluorine atoms is detected. 129 Thus, fast structure conversion is likely to occur and the calculations herein only consider the octahedral structure. The structure optimization at the corresponding level of theory improves the calculated values in some cases compared to the experimental ones as outlined in the second part of the table. The DLU error is again negligible (0.5 ppm). The PBE0 functional reproduces the trend of the xenon fluorides very well while overestimating all shifts by roughly 300 ppm. This is in line with our previous findings that hybrid functionals perform well for main group elements. 86 Note that MP2 leads to very similar results for XeF4 and XeF6 . This also holds for the BP86 and the TPSSh functional. TPSS does not even qualitatively reproduce the trend. Moreover, in case of comparably small basis sets optimized for DFT calculations such as the x2c-type bases additional modifications might be advantageous to improve the basis set convergence and the accuracy. Basis-set extensions for magnetic properties were presented even in non-relativistic quamtum chemistry. 150 Thus, the two-component x2c-type basis set was employed here to ensure more flexibility. This one supplies especially more ptype functions. However, further improvements of this basis set for NMR calculations is in progress in our group and will be reported in the future.

4.2.3

Transition-Metal Oxo Compounds

Scalar-relativistic effects are of crucial importance for the

17

O shifts of transition-metal oxo

2− complexes. 84 Here, we compare different methods to treat these effects in CrO2− 4 , MoO4 , − − − WO2− 4 , MnO4 , TcO4 , ReO4 , FeO4 , RuO4 , OsO4 . The structures of these molecules were

optimized at the DFT level (grid 4a, SCF convergence criterion 10−9 Eh ). The functionals BP86, TPSS and PBE0 were chosen to represent the most common steps of the functional ladder. Non-relativistic calculations were carried out with the def2-TZVP basis 140 for oxygen and the TZVPalls2 basis 141,152 for the transition metals. Relativistic calculations with ECPs were performed using the def2-TZVP basis with Wood-Boring (WB) ECPs 153 and the dhf-

24

ACS Paragon Plus Environment

Page 24 of 62

Page 25 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 4: Comparison of 17 O NMR shifts calculated at various levels of theory and experimental values. The experimental values are taken from Ref. 151. NR denotes the non-relativistic calculations, while WB-ECP and DF-ECP refer to the calculations employing Wood-Boring and Dirac-Fock ECPs for the metal atoms. The TPSS results are given in the Supporting Information. BP86

CrO2− 4 MoO2− 4 WO2− 4 MnO− 4 TcO− 4 ReO− 4 FeO4 RuO4 OsO4

PBE0

NR

WB-ECP

DF-ECP

X2C

DLU

NR

WB-ECP

DF-ECP

X2C

DLU

EXP

782 625 567 1077 860 761 1513 1198 1033

– 554 452 – 736 608 – 1034 868

– 587 458 – 764 595 – 1046 798

765 565 451 1050 769 604 1465 1054 806

765 565 451 1050 769 604 1465 1054 806

845 637 572 1250 919 795 2059 1412 1157

– 561 453 – 784 620 – 1162 877

– 594 458 – 800 605 – 1178 842

824 570 450 1210 803 615 1955 1185 855

824 570 450 1210 803 615 1956 1185 855

799 504 384 1183 713 533 – 1083 760

TZVP 154 basis with Dirac-Fock (DF) ECPs, 155,156 while all-electron relativistic calculations utilized the x2c-TZVPall-2c basis together with the finite nucleus model and the X2C or 2− DLU-X2C Hamiltonian. COSMO was employed to treat the counter ions for CrO2− 4 , MoO4 , − − − WO2− 4 , MnO4 , TcO4 , ReO4 . The obtained structures were confirmed to be minima by

numerical frequency calculations. The deviation between the experimental bond lengths given in Ref. 84 and the calculated ones is below 2.5 pm throughout. FeO4 has been studied theoretically using various methods. 157,158 The tetrahedral structure featuring FeVIII is lower in energy than the C2v structure with FeVI . 158 Thus, we only considered the Td geometry here. To the best of our knowledge, no NMR data exist. The

17

O shifts were calculated

with respect to H2 O in the gas phase. Thus, the experimental data were converted with the experimental vapor/liquid shift of 36 ppm. 159 The calculated NMR shifts are shown in Tab. 4. The results at the BP86 and TPSS level of theory are similar. Hence, the TPSS results are only given in the Supporting Information. The mean absolute error of the DLU scheme amounts to 0.12 ppm, which is negligible considering the range of the

17

O shifts given herein (ca. 800 ppm). The non-relativistic

Hamiltonian and the X2C/DLU-X2C Hamiltonian lead to similar results only in the 3d

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

compounds. For the other systems scalar-relativistic effects have to be treated. Typically, the deviation between the results of X2C/DLU-X2C Hamiltonian and the DF-ECPs is less than 10 ppm. This does not came as a surprise since the basis of both is the four-component Dirac equation. These methods reproduces the experimental NMR shifts within a range of roughly 120 ppm. Previous studies 84 reported a range of ca. 200 ppm solely based on computational methods. The impact of the relativistic effects depends on the functional revealing an interplay between electron correlation and relativistic effects. At first glance, the values obtained at the BP86 level are in better agreement with the experimental ones. However, the PBE0 functional overestimates all shifts and the corresponding error can be corrected systematically. We assume a linear relation between the experimental shifts, δexp , and the calculated shifts, δcalc , described by δexp = a + b · δcalc . A linear regression with the DLU-X2C shifts yields the values a = −76.2 and b = 1.01 with a coefficient of determination of R2 = 0.989. Using these we obtain the following corrected DLU-X2C shifts for PBE0: 756 ppm (CrO2− 4 ), 2− − − − 499 ppm (MoO2− 4 ), 379 ppm (WO4 ), 1146 ppm (MnO4 ), 735 ppm (TcO4 ), 545 ppm (ReO4 ),

1899 ppm (FeO4 ), 1121 ppm (RuO4 ) and 787 ppm (OsO4 ), which is in good agreement with the experimental findings.

4.2.4

Tungsten Compounds

While ECPs can be used to model scalar-relativistic effects originating from the neighboring atoms, an all-electron relativistic Hamiltonian is needed for the shifts of the heavy-element itself. Thus, exemplary tungsten compounds are considered in this section. In view of the results of the transition-metal compounds in Sec. 4.2.3, the same functionals (BP86, TPSS and PBE0) were selected (grid 4a, SCF threshold 10−9 Eh ). The results for BP86 and TPSS are again very similar so that the results for TPSS are given in the Supporting Information. To compare with, non-relativistic calculations were performed with the TZVPalls2 basis for tungsten 152 and the def2-TZVP basis for the light elements. The x2c-TZVPall-2c bases were

26

ACS Paragon Plus Environment

Page 26 of 62

Page 27 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 5: Comparison of non-relativistic and relativistic calculated 183 W shifts (with respect to WO2− 4 ) with experimental data taken from Refs. 160–163. For details on the computational methods see text. BP86 NR WO3 S2− WO2 S2− 2 WOS2− 3 WS2− 4 W(CO)6 WF6 WCl6

X2C

PBE0 DLU

NR

X2C

DLU

EXP

723 736 736 781 775 775 841 1563 1609 1608 1734 1728 1727 1787 2454 2558 2558 2769 2788 2787 2760 3389 3576 3575 3875 3937 3936 3769 −3802 −3828 −3832 −3681 −3816 −3820 −3446 −454 −584 −583 −830 −926 −925 −1121 2544 2226 2226 2522 2086 2086 2181

used in the X2C and DLU-X2C calculations with the finite nucleus model. Again, COSMO was employed to simulate the counter ions for the anions. Structures were optimized at each level of theory and confirmed to be minima by numerical frequency calculations. The NMR shifts are listed in Tab. 5 while the structures and shielding constants are given in the Supporting Information. Again, the mean absolute error of the DLU approximation is very small, 0.8 ppm, compared to the range of the NMR shifts (ca. 7700 ppm). The largest error is found for W(CO)6 with 4 ppm, for other molecules the error is typically below 1 ppm. The treatment of scalar-relativistic effects generally results in a good agreement with the experimental data. A non-relativistic approach is clearly insufficient for WS2− 4 and WCl6 at the BP86 and TPSS level. Here, relativisitc effects amount to 200 and 300 ppm, respectively. The NMR shifts of WF6 and W(CO)6 are the only ones which are not reproduced within a range of ca. 200 ppm. Taking spin-orbit coupling into account does not significantly alter that as the shifts are dominated by the paramagnetic contribution, which contains the perturbed density matrix. 164 The incorporation of HF exchange in the PBE0 functional improves the NMR shifts of WF6 by ca. 300 ppm. When using the PBE0 functional, relativistic effects are important for WF6 and WCl6 while they do not significantly alter the NMR shifts for the other molecules. The PBE0 functional yields better results than the BP86 functional compared to the experiment. However, it does not change the situation for W(CO)6 . The 27

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

− Figure 1: Octahedral silver cluster anions: Ag− 13 (a) and Ag55 (b). These are studied to show the computational efficiency. 1521 and 6435 primitive basis functions are employed, respectively.

deviation from the experimental result still amounts to roughly 400 ppm (ca. 10% of the NMR shift). The deviation from the calculated structure to the experimental one 165 is only 2 pm. The discrepancy between experiment and theory is thus likely to originate from an insufficient treatment of electron correlation, vibrational or environmental effects in the NMR calculation. MP2 does not improve the results but leads to a deterioration. Thus, high-level coupled-cluster calculations may be needed to improve the agreement with the experimental findings. Correcting the NMR shifts of the DLU-X2C-PBE0 method by means of a linear regression as done in Sec. 4.2.3 yields a = 71.5 and b = 0.955 with a coefficient of determination of R2 = 0.996. The corrected NMR shifts are: 812 ppm (WO3 S2− ), 1721 ppm 2− 2− (WO2 S2− 2 ), 2733 ppm (WOS3 ), 3830 ppm (WS4 ), −3576 ppm (W(CO)6 ), −812 ppm (WF6 )

and 2063 ppm (WCl6 ). This is in good agreement with the experimental results.

4.3

Comparison of Computation Times

For an illustrative comparison of the X2C and the DLU-X2C Hamiltonian, we consider the − octahedral silver cluster anions Ag− 13 and Ag55 shown in Fig. 1 at the BP86/x2c-SVPall − level within the RI-J approximation. The computation times for Ag− 13 and Ag55 are given in

28

ACS Paragon Plus Environment

Page 28 of 62

Page 29 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

R Xeon R Processor E5-2687W v2, 3.4 GHz) Table 6: CPU times (single CPU of type Intel for the important steps in a X2C/DLU-X2C NMR shielding calculation of silver cluster anions at the BP86/x2c-SVPall level of theory and organometallic iridum complexes, OLED 1 and 2, at the PBE0/x2c-TZVPPall-2c level. The RI-J approximation is used in the Coulomb integrals and an SCF threshold of 10−8 Eh was applied. The silver cluster Ag− 13 features 1521 primitive basis functions and 1703 RI-auxiliary basis functions, while Ag− uses 55 6435 and 7205 functions, respectively. The clusters are of octahedral symmetry, however, it is not exploited in the calculations. 2319 and 4557 primitive basis functions are employed in the OLED calculations. Note that 9 iterations were needed to solve the CPHF equations. The first-order response with respect to the magnetic moment of the heaviest atom is required in each iteration to check the convergence. The difference between the total computation time and the sum of the steps given below is due to the exchange integrals. Times are given in minutes.

Ag− 13

Ag− 55

OLED 1

OLED 2

Operation

X2C

DLU

X2C

DLU

X2C

DLU

X2C

DLU

1st Order Response B 2nd Order Response RI-J and DFT grid CPHF 1st Order Response m

1.4 101.6 3.0 – 15.3

0.1 3.9 3.0 – 0.6

104.4 31 697.0 80.4 – 4911.7

4.4 355.9 81.0 – 40.6

5.1 1640.7 17.4 846.9 680.9

1.2 89.7 16.8 217.5 30.8

37.1 28 509.1 64.5 8104.4 7397.7

6.4 977.2 79.5 988.1 212.2

Total

121.1

7.6

36 794.8

483.0

2836.6

628.7

38 019.2

3393.1

Tab. 6. For Ag− 13 1521 primitive basis functions are used (689 contracted) while 6435 (2915 contracted) are employed for Ag− 55 . In case of the X2C approach, the computational effort is clearly dominated by the second-order response equations - even for Ag− 13 (83% of the total CPU time). This changes drastically upon the application of the local approximation. Then the one-electron and the Coulomb/DFT part are balanced and the second-order response roughly takes half of the computation time. This allows for application to Ag− 55 with reasonable computational effort for which the NMR shielding calculation still takes roughly 8 hours on a single CPU. Most of the time (73%) is spent in the DLU-X2C part, but this rather reflects the high efficiency of the two-electron part for the density functionals and the RI-J for the Coulomb part. Illustrative calculations on octahedral Ag− 147 show a computation time of less than 8 days and thus revealing the cubic scaling behavior. Again the DLU-X2C part dominates. Matters are different as soon as HF exchange is required then the exchange inte-

29

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

Figure 2: Molecular structure of the organometallic iridum complexes, denoted as OLED 1 (a) and OLED 2 (b) in Ref. 76. OLED 1 consists of 63 atoms while the other complex is build up of 139 atoms. grals and the CPHF procedure will be the most time-consuming step due to the efficiency of the DLU approach. We carried out NMR shift calculations (PBE0, x2c-TZVPPall-2c, RI-J) for the organometallic iridum complexes of Ref. 76, which are displayed in Fig. 2. 2319 (1576 contracted) basis functions were employed for OLED 1 while OLED 2 utilizes 4557 (3229 contracted) functions. The computation time is given in the last columns of Tab. 6. The DLU-X2C calculations take 10 (OLED 1) and 57 hours (OLED 2) on a single CPU of type R Xeon R Processor E5-2687W v2, 3.4 GHz. Here, the second-order response equation Intel

amounts to 14 % (OLED 1) and 29 % (OLED 2), respectively. The CPHF procedure (9 iterations) is more demanding: 48 % and 40 % of the total computation time. The non DLU algorithm reveals a different picture. The second-order X2C response clearly dominates the computation time.

5

SUMMARY AND OUTLOOK

We have presented an implementation of one-electron spin-free exact two-component theory (X2C) for NMR shielding tensors. Efficiency is ensured by the application of the diagonal 30

ACS Paragon Plus Environment

Page 30 of 62

Page 31 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

local approximation to the unitary transformation (DLU) to second derivatives in the spirit of our previous work on geometry gradients. This approximation reduces the formal scaling of the one-electron all-electron relativistic approach from N 4 to N 3 (N measures the system size). This allows for routine applications to large molecules with more than 120 atoms or 3200 contracted basis functions. For systems of this size, the computational effort is dominated by the DLU-X2C part as long as pure density functionals are used together with the resolution-of-the-identity fitting approximation for the Coulomb term (RI-J). For hybrid functionals, which require the calculation of Hartree-Fock exchange and the iterative solution of the coupled-perturbed Hartree-Fock (CPHF) equations, the two-electron part is more consuming than the DLU-X2C part. The accuracy of the approximation was shown for xenon fluorides and various transitionmetal compounds as well as for other organometallic main-group systems with different basis sets. It is thus desirable to extend the implementation herein to a two-component framework utilizing the RI-J approximation to allow for an efficient treatment of spin-orbit effects within a density functional theory framework.

APPENDIX We review and discuss different ways to obtain the derivatives of the coefficients C. We compare a semi-canonical and a non-canonical ansatz for the first- and second-order response equations. The resulting equations are similar to the CPHF equations, 97,98 but can be solved in a non-iterative manner.

First Derivatives The one-electron Dirac equation in a matrix form reads

DC = MCE.

31

ACS Paragon Plus Environment

(45)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 62

The form is the same as the Roothaan-Hall equations 166,167 in HF. In the CPHF framework, the perturbed coefficients Cλ are given as linear-combinations of the unperturbed ones 



λ CL−



λ CS−





λ CL+ 

CL− CL+  λ = U . λ CS+ CS− CS+

(46)

Obviously, the perturbation mixes coefficients of the positive and negative energy states. To find suitable expressions for the matrix Uλ , the terms of the Dirac equation are expanded in powers of the external perturbation λ. This is done in the basis of the unperturbed ˜ = C† DC = E. The energy matrix E is diagonal and the normalization solutions, denoted as D ˜ = C† MC = I, holds. The first-order response equation then becomes condition, M ˜λ − M ˜ λE = E ˜ λ + Uλ E − EUλ . D

(47)

Using the commutator, the expression can be further simplified according to Yoshizawa et al. 81 yielding   λ   λ ˜ −M ˜ λ E − Eλ . U ,E = D

(48)

This equation can be used to calculate the elements of the matrix Uλ . However, for diagonal elements or near-degenerate and degenerate states it will become numerically unstable. Expressions for these elements are derived from the normalization condition and its expansion

˜ λ + Uλ = 0. U†,λ + M

(49)

Assuming canonicity, i.e. Eλ to be diagonal, Eqs. (48) and (49) result in ˜λ − M ˜ λ Eqq D pq pq ∀p 6= q, Eqq − Epp 1 ˜λ = − M . 2 pp

Uλpq =

(50)

Uλpp

(51)

32

ACS Paragon Plus Environment

Page 33 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Note that p and q run over positive and negative energy states. The accuracy of this ansatz relies on the threshold for near-degenerate states. A threshold of 10−4 a.u. was suggested by Cheng and Gauss 68 and also used in our implementation of energy gradients. 76 This threshold was validated by comparison with numerical derivatives. We refer to this ansatz as the semi-canonical one. The derivative of the decoupling matrix X is then obtained according to Cheng and Gauss 104 as

 −1 λ λ CL+ . − XCL+ X λ = CS+

(52)

This relation follows directly from the definition of the decoupling matrix. Similar to the CPHF equations 99 a different formulation without a threshold and the canonicity assumption is possible. A corresponding approach for the derivative of the decoupling matrix has been given by Zou and co-workers. 105 The matrix Uλ is partitioned into different blocks according to the negative and positive energy solutions: 

 λ U−−

Uλ = 

λ U−+ 

λ λ U++ U+−



(53)

The perturbed positive energy coefficients of the large and the small component now contain λ the negative energy coefficients due to the positronic-electronic block U−+ :

λ λ λ CL+ = CL− U−+ + CL+ U++ ,

(54)

λ λ λ = CS− U−+ + CS+ U++ . CS+

(55)

Inserting this into Eq. (52) and using the normalization of the large component on the metric S˜ † CL+ S˜ CL+ = I,

33

ACS Paragon Plus Environment

(56)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 62

the derivative of the decoupling matrix becomes † λ ˜ X λ = (CS− − XCL− ) U−+ CL+ S.

(57)

The elements of the matrix U−+ read λ )kl = (U−+

˜ λ )kl (E++ )ll ˜ λ )kl − (M (D −+ −+ . (E++ )ll − (E−− )kk

(58)

The indices k and l run over the states of the corresponding subblock of the matrix in the λ = 02 is assumed in line with Ref. 99. No problem concerning four-component basis. E−+

numerical instability arises and only one of the four blocks of the matrix Uλ needs to be evaluated. Equations for the positronic-electronic blocks are given in Ref. 105. We have implemented both approaches for geometry gradients. With the given threshold for the semi-canonical solver the results are the same up to the given screening and convergence criteria. However, the non-canonical ansatz is more efficient since the number of required matrix multiplications is smaller and not all blocks are needed. We have thus adopted it in our implementation. The speed-up in case of a DLU-X2C calculation is rather small since the atomic off-diagonal blocks are more demanding. Within the non-canonical ansatz symmetry can now be fully exploited in X2C structure optimizations.

Second Derivatives Second-order response equations to calculate the corresponding derivatives of the coefficients were given by Cheng and Gauss, 168 however, without detailed information about numerical stability or the threshold for near-degeneracy regarding the generalization of Eq. (50) to second derivatives. Since a non-canonical algorithm is desirable due to higher efficiency, we pursue that route and do not consider their approach here in detail. In the spirit of Zou et

34

ACS Paragon Plus Environment

Page 35 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

al. we introduce the ansatz

  ˜ κλ Cκλ = C Uκλ + Uκ Uλ = CU

(59)

˜ κλ can again be partitioned into for the second derivatives of the coefficients. The matrix U the four blocks according to Eq. (53). The corresponding derivative of the decoupling matrix follows as ˜ κλ − U κ U λ − U λ U κ X κλ = (CS− − XCL− )[U −+ −+ ++ −+ ++ † κ λ − U−+ CL+ S˜ CL− U−+

(60)

† † κ λ ˜ S. S˜ CL− U−+ ] CL+ − U−+ CL+

Differentiating Eq. (48) and using the relations from above, the second-order response equation can be written as 81,105 h

i ˜ κλ , E = D ˜ κλ − M ˜ κλ E − Eκλ U + Uκ [Uλ , E] + Uλ [Uκ , E]

(61)

˜λ − M ˜ λ E) + (D ˜λ − M ˜ λ E)Uκ + U (D †,κ

˜κ − M ˜ κ E) + (D ˜κ − M ˜ κ E)Uλ . + U†,λ (D To evaluate the second derivative of the decoupling matrix in Eq. (60) only the positronicelectronic block of this matrix is needed. However, all blocks of the first-order matrices Uλ and Uκ have to be calculated to solve the second-order response equations. In atomic NMR shielding calculations Yoshizawa et al. 81 distinguished between positronic and electronic spinors, where the latter are further partitioned into occupied and unoccupied spinors. This scheme is not easy to apply in case of molecular calculations using contracted basis sets since the X2C part is done in the primitive spherical atomic orbital basis. In the spirit of Ref. 169 and 99, and to ensure numerical stability we make the following assumptions for the

35

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 62

first-order quantities (the same assumptions are made for the perturbation κ)

λ E−+ = 02 ,

(62)

†,λ λ U−− = U−− ,

(63)

†,λ λ . U++ = U++

(64)

This leads to numerically stable expressions for the positronic-positronic and electronicelectronic block 1 ˜λ λ U−− = − M , 2 −− 1 ˜λ λ U++ . = − M 2 ++

(65) (66)

˜ κλ is Using these expressions the positronic-electronic block of the second-order matrix U fully defined: ˜ κλ )kl (U −+

 1 ˜ κλ )kl − (M ˜ κλ )kl (E++ )ll (D = −+ −+ (E++ )ll − (E−− )kk κ λ κ λ λ κ λ κ + U−− U−+ + U−+ U++ + U−− U−+ + U−+ U++

 kl

(E++ )ll

 κ λ κ λ λ κ λ κ − U−− E−− U−+ + U−+ E++ U++ + U−− E−− U−+ + U−+ E++ U++ kl h    i †,κ †,κ ˜ λ ) − (M ˜ λ )(E++ ) + U+− ˜ λ ) − (M ˜ λ )(E++ ) + U−− (D (D −+ −+ ++ ++ kl h    i ˜λ −M ˜ λ E−− U κ + D ˜λ −M ˜ λ E−+ U κ + D −− −− −+ −+ −+ ++ kl    h i †,λ †,λ κ κ κ κ ˜ ˜ ˜ ˜ + U−− (D−+ ) − (M−+ )(E++ ) + U+− (D++ ) − (M++ )(E++ ) kl h    i  ˜κ −M ˜ κ E−− U λ + D ˜κ −M ˜ κ E−+ U λ + D −− −− −+ −+ −+ ++

(67)

kl

We further compared the results of the semi-canonical and the non-canonical approach. Similar to first derivatives, the latter is the method of choice in terms of stability and efficiency.

36

ACS Paragon Plus Environment

Page 37 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Author Information Corresponding Authors *E-mail: [email protected] and [email protected] ORCID Yannick J. Franzke: 0000-0002-8344-113X Florian Weigend: 0000-0001-5060-1689 Funding Y. J. F. acknowledges financial support from Fonds der Chemischen Industrie, Kekulé scholarship. Notes The authors declare no competing financial interest.

Acknowledgement The authors thank Robert Treß for test calculations of the all-electron grids for the DFT part. The authors further thank Patrik Pollak for supplying the non-relativistic all-electron basis sets (TZVPalls2) for tungsten, rhenium and osmium.

Supporting Information Available The Supporting Information is available on the ACS Publications website at DOI: XX/xxxxx • Optimized molecular structures (Cartesian coordinates) and detailed tables of the results of the

13

C,

17

O,

29

Si,

73

Ge,

119

Sn,

129

Xe,

183

W and

207

Pb NMR calculations, and

data for TOC entry (PDF). This material is available free of charge via the Internet at http://pubs.acs.org/.

37

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Friebolin, H. Basic One- and Two-Dimensional NMR Spectroscopy, 5th ed.; WileyVCH: Weinheim, 2010. (2) Dyall, K. G.; Fægri Jr., K. Introduction to Relativistic Quantum Chemistry; Oxford University Press: New York, 2007. (3) Reiher, M.; Wolf., A. Relativistic Quantum Chemistry - The Fundamental Theory of Molecular Science, 2nd ed.; Wiley-VCH: Weinheim, 2015. (4) Liu, W., Ed. Handbook of Relativistic Quantum Chemistry; Springer: Berlin, Heidelberg, 2017. (5) Dolg, M., Ed. Computational Methods in Lanthanide and Actinide Chemistry; John Wiley & Sons: Chichester, 2015. (6) Pyykkö, P. Relativistic effects in structural chemistry. Chem. Rev. 1988, 88, 563–594. (7) Pyykkö, P. Theoretical chemistry of gold. III. Chem. Soc. Rev. 2008, 37, 1967–1997. (8) Saue, T. Relativistic Hamiltonians for Chemistry: A Primer. ChemPhysChem 2011, 12, 3077–3094. (9) Peng, D.; Reiher, M. Exact decoupling of the relativistic Fock operator. Theor. Chem. Acc. 2012, 131, 1081. (10) Pyykkö, P. The Physics behind Chemistry and the Periodic Table. Chem. Rev. 2012, 112, 371–384. (11) Liu, W. Ideas of relativistic quantum chemistry. Mol. Phys. 2010, 108, 1679–1706. (12) Liu, W. Advances in relativistic molecular quantum mechanics. Phys. Rep. 2014, 537, 59 – 89.

38

ACS Paragon Plus Environment

Page 38 of 62

Page 39 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(13) Liu, W. Big picture of relativistic molecular quantum mechanics. Natl. Sci. Rev. 2016, 3, 204–221. (14) Dolg, M.; Cao, X. Relativistic Pseudopotentials: Their Development and Scope of Applications. Chem. Rev. 2012, 112, 403–480. (15) van Wüllen, C. On the use of effective core potentials in the calculation of magnetic properties, such as magnetizabilites and magnetic shieldings. J. Chem. Phys. 2012, 136, 114110. (16) Reiter, K.; Kühn, M.; Weigend, F. Vibrational circular dichroism spectra for large molecules and molecules with heavy elements. J. Chem. Phys. 2017, 146, 054102. (17) Aucar, G. A.; Saue, T.; Visscher, L.; Jensen, H. J. A. On the origin and contribution of the diamagnetic term in four-component relativistic calculations of magnetic properties. J. Chem. Phys. 1999, 110, 6208–6218. (18) Komorovský, S.; Repiský, M.; Malkina, O. L.; Malkin, V. G.; Malkin Ondík, I.; Kaupp, M. A fully relativistic method for calculation of nuclear magnetic shielding tensors with a restricted magnetically balanced basis in the framework of the matrix Dirac-Kohn-Sham equation. J. Chem. Phys. 2008, 128, 104101. (19) Komorovský, S.; Repiský, M.; Malkina, O. L.; Malkin, V. G. Fully relativistic calculations of NMR shielding tensors using restricted magnetically balanced basis and gauge including atomic orbitals. J. Chem. Phys. 2010, 132, 154101. (20) Komorovský, S.; Repiský, M.; Ruud, K.; Malkina, O. L.; Malkin, V. G. FourComponent Relativistic Density Functional Theory Calculations of NMR Shielding Tensors for Paramagnetic Systems. J. Phys. Chem. A 2013, 117, 14209–14219. (21) Iliaš, M.; Saue, T.; Enevoldsen, T.; Jensen, H. J. A. Gauge origin independent calcu-

39

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

lations of nuclear magnetic shieldings in relativistic four-component theory. J. Chem. Phys. 2009, 131, 124119. (22) Olejniczak, M.; Bast, R.; Saue, T.; Pecul, M. A simple scheme for magnetic balance in four-component relativistic Kohn-Sham calculations of nuclear magnetic resonance shielding constants in a Gaussian basis. J. Chem. Phys. 2012, 136, 014108. (23) Olejniczak, M.; Bast, R.; Saue, T.; Pecul, M. Erratum: “A simple scheme for magnetic balance in four-component relativistic Kohn–Sham calculations of nuclear magnetic resonance shielding constants in a Gaussian basis” [J. Chem. Phys. 136, 014108 (2012)]. J. Chem. Phys. 2012, 136, 239902. (24) Xiao, Y.; Peng, D.; Liu, W. Four-component relativistic theory for nuclear magnetic shielding constants: The orbital decomposition approach. J. Chem. Phys. 2007, 126, 081101. (25) Xiao, Y.; Liu, W.; Cheng, L.; Peng, D. Four-component relativistic theory for nuclear magnetic shielding constants: Critical assessments of different approaches. J. Chem. Phys. 2007, 126, 214101. (26) Cheng, L.; Xiao, Y.; Liu, W. Four-component relativistic theory for NMR parameters: Unified formulation and numerical assessment of different approaches. J. Chem. Phys. 2009, 130, 144102. (27) Cheng, L.; Xiao, Y.; Liu, W. Four-component relativistic theory for nuclear magnetic shielding: Magnetically balanced gauge-including atomic orbitals. J. Chem. Phys. 2009, 131, 244113. (28) Xiao, Y.; Sun, Q.; Liu, W. Fully relativistic theories and methods for NMR parameters. Theor. Chem. Acc. 2012, 131, 1080.

40

ACS Paragon Plus Environment

Page 40 of 62

Page 41 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(29) Foldy, L. L.; Wouthuysen, S. A. On the Dirac Theory of Spin 1/2 Particles and Its Non-Relativistic Limit. Phys. Rev. 1950, 78, 29–36. (30) Douglas, M.; Kroll, N. M. Quantum electrodynamical corrections to the fine structure of helium. Ann. Phys. (NY) 1974, 82, 89 – 155. (31) Hess, B. A. Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators. Phys. Rev. A 1986, 33, 3742–3748. (32) Jansen, G.; Hess, B. A. Revision of the Douglas-Kroll transformation. Phys. Rev. A 1989, 39, 6016–6017. (33) Nakajima, T.; Hirao, K. Numerical illustration of third-order Douglas-Kroll method: atomic and molecular properties of superheavy element 112. Chem. Phys. Lett. 2000, 329, 511 – 516. (34) Nakajima, T.; Hirao, K. The higher-order Douglas-Kroll transformation. J. Chem. Phys. 2000, 113, 7786–7789. (35) Wolf, A.; Reiher, M.; Hess, B. A. The generalized Douglas-Kroll transformation. J. Chem. Phys. 2002, 117, 9215–9226. (36) van Wüllen, C. Relation between different variants of the generalized Douglas-Kroll transformation through sixth order. J. Chem. Phys. 2004, 120, 7307–7313. (37) Reiher, M.; Wolf, A. Exact decoupling of the Dirac Hamiltonian. I. General theory. J. Chem. Phys. 2004, 121, 2037–2047. (38) Reiher, M.; Wolf, A. Exact decoupling of the Dirac Hamiltonian. II. The generalized Douglas-Kroll-Hess transformation up to arbitrary order. J. Chem. Phys. 2004, 121, 10945–10956.

41

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(39) Wolf, A.; Reiher, M. Exact decoupling of the Dirac Hamiltonian. III. Molecular properties. J. Chem. Phys. 2006, 124, 064102. (40) Wolf, A.; Reiher, M. Exact decoupling of the Dirac Hamiltonian. IV. Automated evaluation of molecular properties within the Douglas-Kroll-Hess theory up to arbitrary order. J. Chem. Phys. 2006, 124, 064103. (41) Reiher, M.; Wolf, A. Regular no-pair Dirac operators: Numerical study of the convergence of high-order Douglas-Kroll-Hess transformations. Phys. Lett. A 2007, 360, 603 – 607. (42) Peng, D.; Hirao, K. An arbitrary order Douglas-Kroll method with polynomial cost. J. Chem. Phys. 2009, 130, 044102. (43) Fukuda, R.; Hada, M.; Nakatsuji, H. Quasirelativistic theory for the magnetic shielding constant. I. Formulation of Douglas-Kroll-Hess transformation for the magnetic field and its application to atomic systems. J. Chem. Phys. 2003, 118, 1015–1026. (44) Fukuda, R.; Hada, M.; Nakatsuji, H. Quasirelativistic theory for magnetic shielding constants. II. Gauge-including atomic orbitals and applications to molecules. J. Chem. Phys. 2003, 118, 1027–1035. (45) Hayami, M.; Seino, J.; Nakai, H. Gauge-origin independent formalism of twocomponent relativistic framework based on unitary transformation in nuclear magnetic shielding constant. J. Chem. Phys. 2018, 148, 114109. (46) Chang, C.; Pelissier, M.; Durand, P. Regular Two-Component Pauli-Like Effective Hamiltonians in Dirac Theory. Phys. Scr. 1986, 34, 394–404. (47) van Lenthe, E.; Baerends, E. J.; Snijders, J. G. Relativistic regular two-component Hamiltonians. J. Phys. Chem. 1993, 99, 4597–4610.

42

ACS Paragon Plus Environment

Page 42 of 62

Page 43 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(48) van Lenthe, E.; Baerends, E. J.; Snijders, J. G. Relativistic total energy using regular approximations. J. Chem. Phys. 1994, 101, 9783–9792. (49) Wolff, S. K.; Ziegler, T.; van Lenthe, E.; Baerends, E. J. Density functional calculations of nuclear magnetic shieldings using the zeroth-order regular approximation (ZORA) for relativistic effects: ZORA nuclear magnetic resonance. J. Chem. Phys. 1999, 110, 7689–7698. (50) van Leeuwen, R.; van Lenthe, E.; Baerends, E. J.; Snijders, J. G. Exact solutions of regular approximate relativistic wave equations for hydrogen-like atoms. J. Chem. Phys. 1994, 101, 1272–1281. (51) Filatov, M.; Cremer, D. A gauge-independent zeroth-order regular approximation to the exact relativistic Hamiltonian-Formulation and applications. J. Chem. Phys. 2005, 122, 044104. (52) Autschbach, J. The role of the exchange-correlation response kernel and scaling corrections in relativistic density functional nuclear magnetic shielding calculations with the zeroth-order regular approximation. Mol. Phys. 2013, 111, 2544–2554. (53) Aucar, G. A.; Melo, J. I.; Aucar, I. A.; Maldonado, A. F. Foundations of the LRESC model for response properties and some applications. Int. J. Quantum Chem. 2018, 118, e25487. (54) Kutzelnigg, W.; Liu, W. Quasirelativistic theory equivalent to fully relativistic theory. J. Chem. Phys. 2005, 123, 241102. (55) Kutzelnigg, W.; Liu, W. Quasirelativistic theory I. Theory in terms of a quasirelativistic operator. Mol. Phys. 2006, 104, 2225–2240. (56) Liu, W.; Peng, D. Infinite-order quasirelativistic density functional method based on the exact matrix quasirelativistic theory. J. Chem. Phys. 2006, 125, 044102. 43

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(57) Liu, W.; Peng, D. Erratum: “Infinite-order quasirelativistic density functional method based on the exact matrix quasirelativistic theory” [J. Chem. Phys. 125, 044102 (2006)]. J. Chem. Phys. 2006, 125, 149901. (58) Iliaš, M.; Saue, T. An infinite-order two-component relativistic Hamiltonian by a simple one-step transformation. J. Chem. Phys. 2007, 126, 064102. (59) Liu, W.; Kutzelnigg, W. Quasirelativistic theory. II. Theory at matrix level. J. Chem. Phys. 2007, 126, 114107. (60) Peng, D.; Liu, W.; Xiao, Y.; Cheng, L. Making four- and two-component relativistic density functional methods fully equivalent based on the idea of “from atoms to molecule”. J. Chem. Phys. 2007, 127, 104106. (61) Liu, W.; Peng, D. Exact two-component Hamiltonians revisited. J. Chem. Phys. 2009, 131, 031104. (62) Sikkema, J.; Visscher, L.; Saue, T.; Iliaš, M. The molecular mean-field approach for correlated relativistic calculations. J. Chem. Phys. 2009, 131, 124116. (63) Dyall, K. G. Interfacing relativistic and nonrelativistic methods. I. Normalized elimination of the small component in the modified Dirac equation. J. Chem. Phys. 1997, 106, 9618–9626. (64) Dyall, K. G. Interfacing relativistic and nonrelativistic methods. II. Investigation of a low-order approximation. J. Chem. Phys. 1998, 109, 4201–4208. (65) Dyall, K. G.; Enevoldsen, T. Interfacing relativistic and nonrelativistic methods. III. Atomic 4-spinor expansions and integral approximations. J. Chem. Phys. 1999, 111, 10000–10007. (66) Dyall, K. G. Interfacing relativistic and nonrelativistic methods. IV. One- and twoelectron scalar approximations. J. Chem. Phys. 2001, 115, 9136–9143. 44

ACS Paragon Plus Environment

Page 44 of 62

Page 45 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(67) Cremer, D.; Zou, W.; Filatov, M. Dirac-exact relativistic methods: the normalized elimination of the small component method. WIREs Comput. Mol. Sci. 2014, 4, 436– 467. (68) Cheng, L.; Stopkowicz, S.; Gauss, J. Analytic energy derivatives in relativistic quantum chemistry. Int. J. Quantum Chem. 2014, 114, 1108–1127. (69) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Integral Approximations for LCAO-SCF Calculations. Chem. Phys. Lett. 1993, 213, 514–518. (70) Eichkorn, K.; Treutler, O.; Öhm, H.; Häser, M.; Ahlrichs, R. Auxiliary Basis Sets to Approximate Coulomb Potentials. Chem. Phys. Lett. 1995, 242, 283–290. (71) Eichkorn K.,; Weigend, F.; Treutler, O.; Ahlrichs, R. Auxiliary Basis Sets for Main Row Atoms and Transition Metals and Their Use to Approximate Coulomb Potentials. Theor. Chem. Acc. 1997, 97, 119–124. (72) Weigend, F. Accurate Coulomb-fitting basis sets for H to Rn. Phys. Chem. Chem. Phys. 2006, 8, 1057–1065. (73) Weigend, F.; Kattannek, M.; Ahlrichs, R. Approximated electron repulsion integrals: Cholesky decomposition versus resolution of the identity methods. J. Chem. Phys. 2009, 130, 164106. (74) Peng, D.; Reiher, M. Local relativistic exact decoupling. J. Chem. Phys. 2012, 136, 244108. (75) Peng, D.; Middendorf, N.; Weigend, F.; Reiher, M. An efficient implementation of two-component relativistic exact-decoupling methods for large molecules. J. Chem. Phys. 2013, 138, 184105. (76) Franzke, Y. J.; Middendorf, N.; Weigend, F. Efficient implementation of one- and

45

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

two-component analytical energy gradients in exact two-component theory. J. Chem. Phys. 2018, 148, 104410. (77) Lichtenberger, N.; Franzke, Y. J.; Massa, W.; Weigend, F.; Dehnen, S. The Identity of “Ternary” A/Tl/Pb or K/Tl/Bi Solid Mixtures and Binary Zintl Anions Isolated From Their Solutions. Chem. Eur. J. 2018, 24, 12022–12030. (78) Sun, Q.; Liu, W.; Xiao, Y.; Cheng, L. Exact two-component relativistic theory for nuclear magnetic resonance parameters. J. Chem. Phys. 2009, 131, 081101. (79) Sun, Q.; Xiao, Y.; Liu, W. Exact two-component relativistic theory for NMR parameters: General formulation and pilot application. J. Chem. Phys. 2012, 137, 174105. (80) Cheng, L.; Gauss, J.; Stanton, J. F. Treatment of scalar-relativistic effects on nuclear magnetic shieldings using a spin-free exact-two-component approach. J. Chem. Phys. 2013, 139, 054105. (81) Yoshizawa, T.; Zou, W.; Cremer, D. Calculations of atomic magnetic nuclear shielding constants based on the two-component normalized elimination of the small component method. J. Chem. Phys. 2017, 146, 134109. (82) Yoshizawa, T.; Hada, M. Calculations of nuclear magnetic shielding constants based on the exact two-component relativistic method. J. Chem. Phys. 2017, 147, 154104. (83) Demissie, T. B. Relativistic effects on the NMR parameters of Si, Ge, Sn, and Pb alkynyl compounds: Scalar versus spin-orbit effects. J. Chem. Phys. 2017, 147, 174301. (84) Kaupp, M.; Malkin, V. G.; Malkina, O. L.; Salahub, D. R. Scalar Relativistic Effects on 17O NMR Chemical Shifts in Transition-Metal Oxo Complexes. An ab Initio ECP/DFT Study. J. Am. Chem. Soc. 1995, 117, 1851–1852.

46

ACS Paragon Plus Environment

Page 46 of 62

Page 47 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(85) Vícha, J.; Komorovský, S.; Repiský, M.; Marek, R.; Straka, M. Relativistic Spin-Orbit Heavy Atom on the Light Atom NMR Chemical Shifts: General Trends Across the Periodic Table Explained. J. Chem. Theory Comput. 2018, 14, 3025–3039. (86) Reiter, K.; Mack, F.; Weigend, F. Calculation of Magnetic Shielding Constants with meta-GGA Functionals Employing the Multipole-Accelerated Resolution of the Identity: Implementation and Assessment of Accuracy and Efficiency. J. Chem.Theory Comput. 2018, 14, 191–197. (87) Bohr, A.; Weisskopf, V. F. The Influence of Nuclear Structure on the Hyperfine Structure of Heavy Elements. Phys. Rev. 1950, 77, 94–98. (88) Visscher, L.; Dyall, K. Dirac-Fock atomic electronic structure calculations using different nuclear charge distributions. At. Data Nucl. Data Tables 1997, 67, 207 – 224. (89) Hennum, A. C.; Klopper, W.; Helgaker, T. Direct perturbation theory of magnetic properties and relativistic corrections for the point nuclear and Gaussian nuclear models. J. Chem. Phys. 2001, 115, 7356–7363. (90) Kutzelnigg, W. Diamagnetism in relativistic theory. Phys. Rev. A 2003, 67, 032109. (91) London, F. Théorie quantique des courants interatomiques dans les combinaisons aromatiques. J. Phys. Radium 1937, 8, 397–409. (92) Ditchfield, R. Self-consistent perturbation theory of diamagnetism. Mol. Phys. 1974, 27, 789–807. (93) Heully, J.-L.; Lindgren, I.; Lindroth, E.; Lundqvist, S.; Mårtensson-Pendrill, A.-M. Diagonalisation of the Dirac Hamiltonian as a basis for a relativistic many-body procedure. J. Phys. B: At. Mol. Phys. 1986, 19, 2799–2815. (94) Gagliardi, L.; Handy, N. C.; Ioannou, A. G.; Skylaris, C.-K.; Spencer, S.; Willetts, A.;

47

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Simper, A. M. A two-centre implementation of the Douglas-Kroll transformation in relativistic calculations. Chem. Phys. Lett. 1998, 283, 187 – 193. (95) Pollak, P.; Weigend, F. Segmented Contracted Error-Consistent Basis Sets of Doubleand Triple-ζ Valence Quality for One- and Two-Component Relativistic All-Electron Calculations. J. Chem. Theory Comput. 2017, 13, 3696–3705. (96) Peterson, K. A. Systematically convergent basis sets with relativistic pseudopotentials. I. Correlation consistent basis sets for the post-d group 13-15 elements. J. Chem. Phys. 2003, 119, 11099–11112. (97) Stevens, R. M.; Pitzer, R. M.; Lipscomb, W. N. Perturbed Hartree-Fock Calculations. I. Magnetic Susceptibility and Shielding in the LiH Molecule. J. Chem. Phys. 1963, 38, 550–560. (98) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Derivative studies in hartreefock and Møller-plesset theories. Int. J. Quantum Chem. 1979, 16, 225–241. (99) Handy, N. C.; Schaefer, H. F. On the evaluation of analytic energy derivatives for correlated wave functions. J. Chem. Phys. 1984, 81, 5031–5033. (100) Malkin, V.; Malkina, O.; Salahub, D. Calculations of NMR shielding constants by uncoupled density functional theory. Chem. Phys. Lett. 1993, 204, 80 – 86. (101) Schreckenbach, G.; Ziegler, T. Calculation of NMR Shielding Tensors Using GaugeIncluding Atomic Orbitals and Modern Density Functional Theory. J. Chem. Phys. 1995, 99, 606–611. (102) Wolinski, K.; Hinton, J. F.; Pulay, P. Efficient implementation of the gaugeindependent atomic orbital method for NMR chemical shift calculations. J. Am. Chem. Soc. 1990, 112, 8251–8260.

48

ACS Paragon Plus Environment

Page 48 of 62

Page 49 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(103) Häser, M.; Ahlrichs, R.; Baron, H. P.; Weis, P.; Horn, H. Direct computation of secondorder scf properties of large molecules on workstation computers with an application to large carbon clusters. Theor. Chim. Acta 1992, 83, 455–470. (104) Cheng, L.; Gauss, J. Analytic energy gradients for the spin-free exact two-component theory using an exact block diagonalization for the one-electron Dirac Hamiltonian. J. Chem. Phys. 2011, 135, 084114. (105) Zou, W.; Filatov, M.; Cremer, D. Development, Implementation, and Application of an Analytic Second Derivative Formalism for the Normalized Elimination of the Small Component Method. J. Chem. Theory Comput. 2012, 8, 2617–2629. (106) Kollwitz, M.; Gauss, J. A direct implementation of the GIAO-MBPT(2) method for calculating NMR chemical shifts. Application to the naphthalenium and anthracenium ions. Chem. Phys. Lett. 1996, 260, 639 – 646. (107) Kollwitz, M.; Häser, M.; Gauss, J. Non-Abelian point group symmetry in direct second-order many-body perturbation theory calculations of NMR chemical shifts. J. Chem. Phys. 1998, 108, 8295–8301. (108) Huniar, U. Berechnung der chemischen Verschiebung der NMR mit Methoden der Dichtefunktionaltheorie (DFT). Diploma Thesis, University of Karlsruhe, Germany, 1999. (109) Local version of TURBOMOLE V7.3 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. (110) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165 – 169.

49

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(111) Furche, F.; Ahlrichs, R.; Hättig, C.; Klopper, W.; Sierka, M.; Weigend, F. Turbomole. WIREs Comput. Mol. Sci. 2014, 4, 91–100. (112) King, H.; Dupuis, M. Numerical integration using rys polynomials. J. Comput. Phys. 1976, 21, 144–165. (113) Dupuis, M.; Rys, J.; King, H. F. Evaluation of molecular integrals over Gaussian basis functions. J. Chem. Phys. 1976, 65, 111–116. (114) Rys, J.; Dupuis, M.; King, H. F. Computation of electron repulsion integrals using the rys quadrature method. J. Comput. Chem. 1983, 4, 154–157. (115) Anderson, E.; Bai, Z.; Bischof, C.; Blackford, S.; Demmel, J.; Dongarra, J.; Du Croz, J.; Greenbaum, A.; Hammarling, S.; McKenney, A.; Sorensen, D. LAPACK Users’ Guide, 3rd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1999. (116) Schäfer, A.; Klamt, A.; Sattel, D.; Lohrenz, J. C. W.; Eckert, F. COSMO Implementation in TURBOMOLE: Extension of an efficient quantum chemical code towards liquid systems. Phys. Chem. Chem. Phys. 2000, 2, 2187–2193. (117) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143 – 152. (118) Dyall, K. G. Relativistic and nonrelativistic finite nucleus optimized triple-zeta basis sets for the 4p, 5p and 6p elements. Theor. Chem. Acc. 2002, 108, 335–340. (119) Dyall, K. G. Erratum: Relativistic and nonrelativistic finite nucleus optimized triple zeta basis sets for the 4p, 5p and 6p elements. Theor. Chem. Acc. 2003, 109, 284–284. (120) Dyall, K. G. Relativistic Quadruple-Zeta and Revised Triple-Zeta and Double-Zeta Basis Sets for the 4p, 5p, and 6p Elements. Theor. Chem. Acc. 2006, 115, 441–447.

50

ACS Paragon Plus Environment

Page 50 of 62

Page 51 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(121) Dyall, K. G. Relativistic double-zeta, triple-zeta, and quadruple-zeta basis sets for the 5d elements Hf–Hg. Theor. Chem. Acc. 2004, 112, 403–409. (122) Basis sets available from the Dirac web site, http://dirac.chem.sdu.dk. (123) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. (124) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (125) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. (126) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (127) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. V. Core-valence basis sets for boron through neon. J. Chem. Phys. 1995, 103, 4572–4585. (128) Peterson, K. E.; Dunning, T. H. Accurate correlation consistent basis sets for molecular core-valence correlation effects: The second row atoms Al-Ar, and the first row atoms B-Ne revisited. J. Chem. Phys. 2002, 117, 10548–10560. (129) Gerken, M.; Hazendonk, P.; Nieboer, J.; Schrobilgen, G. J. NMR spectroscopic study of xenon fluorides in the gas phase and of XeF2 in the solid state. J. Fluor. Chem. 2004, 125, 1163 – 1168. (130) Jameson, C. In Multinuclear NMR; Mason, J., Ed.; Springer US: New York, 1987; Chapter 18, pp 463–477.

51

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(131) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (132) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 33, 8822–8824. (133) Keal, T. W.; Tozer, D. J. A semiempirical generalized gradient approximation exchange-correlation functional. J. Chem. Phys. 2004, 121, 5654–5660. (134) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 146401. (135) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (136) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (137) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627. (138) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes. J. Chem. Phys. 2003, 119, 12129–12137. (139) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. Main Group Atoms and Dimers Studied with a New Relativistic ANO Basis Set. J. Phys. Chem. A 2004, 108, 2851–2858. (140) Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and

52

ACS Paragon Plus Environment

Page 52 of 62

Page 53 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. (141) Ahlrichs, R.; May, K. Contracted all-electron Gaussian basis sets for atoms Rb to Xe. Phys. Chem. Chem. Phys. 2000, 2, 943–945. (142) Bürger, H.; Kuna, R.; Ma, S.; Breidung, J.; Thiel, W. The vibrational spectra of krypton and xenon difluoride: High-resolution infrared studies and ab initio calculations. J. Chem. Phys. 1994, 101, 1–14. (143) Bürger, H.; Ma, S.; Breidung, J.; Thiel, W. Ab initio calculations and high resolution infrared investigation on XeF4. J. Chem. Phys. 1996, 104, 4945–4953. (144) Gawrilow, M.; Beckers, H.; Riedel, S.; Cheng, L. Matrix-Isolation and QuantumChemical Analysis of the C3v Conformer of XeF6, XeOF4, and Their Acetonitrile Adducts. J. Phys. Chem. A 2018, 122, 119–129. (145) Crawford, T. D.; Springer, K. W.; Schaefer, H. F. A contribution to the understanding of the structure of xenon hexafluoride. J. Chem. Phys. 1995, 102, 3307–3311. (146) Kaupp, M.; van Wüllen, C.; Franke, R.; Schmitz, F.; Kutzelnigg, W. The Structure of XeF6 and of Compounds Isoelectronic with It. A Challenge to Computational Chemistry and to the Qualitative Theory of the Chemical Bond. J. Am. Chem. Soc. 1996, 118, 11939–11950. (147) Dixon, D. A.; de Jong, W. A.; Peterson, K. A.; Christe, K. O.; Schrobilgen, G. J. Heats of Formation of Xenon Fluorides and the Fluxionality of XeF6 from High Level Electronic Structure Calculations. J. Am. Chem. Soc. 2005, 127, 8627–8634. (148) Peterson, K. A.; Dixon, D. A.; Stoll, H. The Use of Explicitly Correlated Methods on XeF6 Predicts a C3v Minimum with a Sterically Active, Free Valence Electron Pair on Xe. J. Phys. Chem. A 2012, 116, 9777–9782. 53

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(149) Cheng, L.; Gauss, J.; Stanton, J. F. Relativistic coupled-cluster calculations on XeF6: Delicate interplay between electron-correlation and basis-set effects. J. Chem. Phys. 2015, 142, 224309. (150) Jensen, F. Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods. J. Chem. Theory Comput. 2008, 4, 719–727. (151) Figgis, B. N.; Kidd, R. G.; Nyholm, R. S. Oxygen-17 nuclear magnetic resonance of inorganic compounds. Proc. Royal Soc. Lond. A 1962, 269, 469–480. (152) Private communication with Patrik Pollak, the TZVPalls2 basis sets for W, Re and Os are given in the Supporting Information. (153) Andrae, D.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H. Energy-adjusted ab initio pseudopotentials for the second and third row transition elements. Theor. Chim. Acta 1990, 77, 123–141. (154) Weigend, F.; Baldes, A. Segmented contracted basis sets for one- and two-component Dirac-Fock effective core potentials. J. Chem. Phys. 2010, 133, 174102. (155) Peterson, K. A.; Figgen, D.; Dolg, M.; Stoll, H. Energy-consistent relativistic pseudopotentials and correlation consistent basis sets for the 4d elements Y-Pd. J. Chem. Phys. 2007, 126, 124101. (156) Figgen, D.; Peterson, K. A.; Dolg, M.; Stoll, H. Energy-consistent pseudopotentials and correlation consistent basis sets for the 5d elements Hf-Pt. J. Chem. Phys. 2009, 130, 164108. (157) Tran, V. T.; Hendrickx, M. F. A. Description of the geometric and electronic structures responsible for the photoelectron spectrum of FeO4-. J. Chem. Phys. 2011, 135, 094505.

54

ACS Paragon Plus Environment

Page 54 of 62

Page 55 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(158) Huang, W.; Pyykkö, P.; Li, J. Is Octavalent Pu(VIII) Possible? Mapping the Plutonium Oxyfluoride Series PuOnF8-2n (n = 0–4). Inorg. Chem. 2015, 54, 8825–8831. (159) McFarlane, W.; McFarlane, H. C. E. In Multinuclear NMR; Mason, J., Ed.; Springer US: New York, 1987; Chapter 14, pp 403–416. (160) Mann, B. In Annual Reports on NMR Spectroscopy; Webb, G., Ed.; Academic Press, 1991; Vol. 23; pp 141 – 207. (161) Banck, J.; Schwenk, A. 183Tungsten NMR studies. Z. Phys. B 1975, 20, 75–80. (162) Gheller, S. F.; Hambley, T. W.; Rodgers, J. R.; Brownlee, R. T. C.; O’Connor, M. J.; Snow, M. R.; Wedd, A. G. Synthesis and characterization of complexes of thiomolybdates and thiotungstates with copper(I) and silver(I) cyanides, including molybdenum-95 and tungsten-183 NMR properties and the crystal and molecular structures of (n-Pr4N)2[(CN)CuS2MoS2], (n-Pr4N)2[(CN)AgS2WS2], and (Ph4As)2[(CN)CuS2MoS2Cu(CN)].H2O. Inorg. Chem. 1984, 23, 2519–2528. (163) Keiter, R. L.; Velde, D. G. V. Direct observation of 183W NMR of W(CO)6 and phosphine derivatives. J. Organomet. Chem. 1983, 258, C34 – C36. (164) Rodriguez-Fortea, A.; Alemany, P.; Ziegler, T. Density Functional Calculations of NMR Chemical Shifts with the Inclusion of Spin-Orbit Coupling in Tungsten and Lead Compounds. J. Phys. Chem. A 1999, 103, 8288–8294. (165) Elschenbroich, C. Organometallchemie, 6th ed.; Vieweg+Teubner: Wiesbaden, 2008; p 330. (166) Roothaan, C. C. J. New Developments in Molecular Orbital Theory. Rev. Mod. Phys. 1951, 23, 69–89. (167) Hall, G. G. The Molecular Orbital Theory of Chemical Valency. VIII. A Method of Calculating Ionization Potentials. Proc. Roy. Soc. Lond. A 1951, 205, 541–552. 55

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(168) Cheng, L.; Gauss, J. Analytic second derivatives for the spin-free exact two-component theory. J. Chem. Phys. 2011, 135, 244104. (169) O’Shea, S. F.; Santry, D. P. Non-empirical molecular orbital theory of the electronic structure of molecular crystals. Theor. Chim. Acta 1975, 37, 1–16.

56

ACS Paragon Plus Environment

Page 56 of 62

Page 57 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Graphical TOC Entry

57

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

ACS Paragon Plus Environment

Page 58 of 62

Page 59 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

ACS Paragon Plus Environment

Page 60 of 62

Page 61 of 62 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

89x36mm (150 x 150 DPI)

ACS Paragon Plus Environment

Page 62 of 62