Efficient and Accurate Prediction of Nuclear ... - ACS Publications

Subscriber access provided by University of South Dakota

Quantum Electronic Structure

Efficient and accurate prediction of nuclear magnetic resonance shielding tensors with double-hybrid density functional theory Georgi L Stoychev, Alexander A Auer, and Frank Neese J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00624 • Publication Date (Web): 26 Jul 2018 Downloaded from http://pubs.acs.org on July 27, 2018

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 51 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

Efficient and Accurate Prediction of Nuclear Magnetic Resonance Shielding Tensors with Double-Hybrid Density Functional Theory Georgi L. Stoychev, Alexander A. Auer, and Frank Neese∗ Max-Planck-Institut für Kohlenforschung, Mülheim an der Ruhr, Germany

Abstract Analytic calculation of nuclear magnetic resonance chemical shielding tensors, based on gauge-including atomic orbitals, is implemented for double-hybrid density functional theory (DHDFT), using the resolution of the identity (RI) approximation for its second order Møller–Plesset perturbation theory (MP2) correlation contributions. A benchmark set of 15 small molecules, containing 1 H,

13 C, 15 N, 17 O, 19 F,

and

31 P

nuclei, is

used to assess the accuracy of the results in comparison to coupled cluster and empirical equilibrium reference data, as well as to calculations with MP2, Hartree–Fock, and commonly used pure and hybrid density functionals. Investigated are also errors due to basis set incompleteness, the frozen core approximation, different auxiliary basis sets for the RI approximation, and grids used for the chain-of-spheres exchange integral evaluation. The DSD-PBEP86 double-hybrid functional shows the smallest deviations from the reference data with mean absolute relative error in chemical shifts of 1.9 %. This is significantly better than MP2 (4.1 %), spin-component-scaled MP2 (3.9 %), or the best conventional density functional tested, M06L (5.4 %). A protocol (basis sets, grid sizes, etc.) for the efficient and accurate calculation of chemical shifts at the DHDFT level is proposed and shown to be routinely applicable to systems of 100–400 electrons,

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

requiring computation times one to two orders of magnitude longer than for equivalent calculations with conventional (pure or hybrid) density functionals.

1

Introduction

Nuclear magnetic resonance (NMR) chemical shielding constants (CSCs) provide a valuable link between theory and experiment. Their calculation using correlated wave function-based methods and analytic derivative techniques was pioneered by Gauss, first at the level of second order Møller–Plesset perturbation theory (MP2), 1,2 and subsequently at higher orders of perturbation theory and coupled cluster theory. 3–6 MP2 has been shown to offer a good balance between cost and accuracy in the calculated NMR CSCs. 1,2,7–11 Employing a local correlation approach and the resolution of the identity (RI) approximation further reduces the computational effort, allowing MP2 level CSCs to be calculated for systems of a few hundred electrons and a few thousand basis functions within several days of computation time. 12,13 The introduction of linear and sub-linear scaling methodologies by Ochsenfeld and coworkers may make even larger systems accessible, especially if only a few nuclei are of interest. 14–17 A somewhat empirical way to improve MP2 results is to scale the samespin (SS) and opposite-spin (OS) contributions to the energy by different factors in the so-called spin-component-scaled MP2 (SCS-MP2) approach. 18–20 The related spin-oppositescaled MP2 (SOS-MP2) method completely neglects the SS contribution and provides similar accuracy with the added benefit of reduced formal scaling of the computational effort with system size from O (N 5 ) to O (N 4 ). Fitting the scaling parameters in SCS-MP2 (or SOSMP2) to better reproduce CCSD(T) (coupled cluster theory including singles doubles and perturbative triples) reference values is also possible, 21 although this somewhat obscures the physical meaning behind these parameters. A rather more popular approach to the calculation of CSCs is density functional theory (DFT), which has a significantly lower computational cost. The proper treatment of

2


Page 2 of 51

Page 3 of 51 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


magnetic perturbations in the context DFT is the subject of ongoing discussion. 22–26 Most common density functionals (DFs) do not depend on the external magnetic field, and therefore produce exchange and correlation (XC) energies, which are unphysically constant in the presence of such. While it is possible to introduce an explicit magnetic field dependence, 27 the more widely-accepted approach is current-density functional theory (CDFT), pioneered by Vignale and Rasolt, 28,29 who introduced the paramagnetic current-density as an independent variable, although this choice is also subject to debate. 30 CDFT is still actively developed and was recently extended to functionals based on the meta-generalized-gradient approximation (meta-GGA). 26,31–34 However, while the current-density contributions to NMR CSCs can be substantial, 23,26 CDFT results are not necessarily substantially better than those obtained with standard current-independent DFs. 22,34,35 Therefore, for pragmatic reasons the latter are still widely used for the calculation of magnetic properties with different functionals showing varying degrees of accuracy. 11,36,37 While there are DFs specifically optimized for CSC calculations, 38,39 some general application meta-GGAs like TPSS, VS98, and M06L have been shown to be particularly well-suited for NMR chemical shifts. 37,40 According to the “Jacob’s ladder” classification, introduced by Perdew, 41 the highest level (fifth rung) DFT methodologies include a non-local correlation energy contribution by taking into account the virtual molecular orbitals (MOs). One possible approach of this type is double-hybrid DFT (DHDFT), whereby an MP2-like term is added to the total energy. 42,43 Combining this with SCS-MP2 and the empirical dispersion correction (denoted D3BJ in its most popular formulation), introduced by Grimme, 44–46 gives the general formulation of dispersion-corrected spin-component-scaled double-hybrid DFT (DSD-DFT), developed by Kozuch and Martin. 47–49 In extensive benchmark studies of themochemistry, kinetics and noncovalent interactions, double-hybrid density functionals (DHDFs) have been clearly shown to outperform lower-rung DFs, with DSD-BLYP-D3BJ and DSD-PBEP86-D3BJ producing the most accurate energies. 50,51 DSD-DFT functionals also perform better than other common DFs or MP2 in harmonic frequency calculations, for which they were not specifi-

3



cally optimized. 48,49 This suggests they may be sufficiently “universal” to also produce high quality results for other response properties. Note that, as discussed above, the XC functionals used in DHDFT are independent of the magnetic field. While it may be expected that this deficiency is partly offset by the inclusion of Hartree–Fock (HF) exchange and MP2 correlation, this issue is not examined in the current work. The efficient calculation of molecular properties in the framework of DHDFT requires the development and implementation of analytic derivatives. 52–54 An additional complication in the case of NMR CSCs is the gauge origin problem, 55 for which the most widely accepted solution is the use of gauge-including atomic orbitals (GIAOs). 56–60 Hence, the goal of this work is to evaluate the performance of DHDFT for NMR CSC calculations using GIAOs, which have not been reported previously, to the best of our knowledge.

2

Theory

The theory of MP2-level NMR CSC calculations is well known 1,2,61 and adaptations to SCS, 21 local, 12,17 and RI-accelerated local MP2 13 have been developed and implemented. The extension to DHDFT is rather simple and in general presents fewer difficulties than geometric DHDFT derivatives. 52,54 For detailed discussions, the reader is referred to these previous derivations. This section aims to present a concise summary of the working equations, as implemented in the ORCA program, 62,63 in order to provide sufficient context for interpretation of the results.

2.1

Derivative Approach to NMR Chemical Shielding

The chemical shielding tensor of nucleus A can be calculated as the second derivative of the energy with respect to the external magnetic field B and the nuclear magnetic moment MA .

4


Page 4 of 51

Page 5 of 51 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


The following simple expression is obtained:

¯ A = E¯ B,MA = σ

X

¯ B,MA + ¯ µν h D µν

µν

X

B ¯ MA ¯ µν hµν D

(1)

µν

¯ are respectively the method-specific shielding tensor, energy expres¯ and h ¯ D, ¯ E, where σ, sion, density matrix in the atomic orbital (AO) basis, and one-electron part of the Fock matrix and the superscripts denote derivatives. Note that for non-variational methods, the proper definitions of the density matrices must be used. 1,4,6,64–68 The general DSD-DFT energy expression is 47–49

E DSD-DFT = TS + J + Ene + cX EXHF + (1 − cX ) EXDFT + cC ECDFT + cO EOMP2 + cS ESMP2 +s6 ED {z } {z } | | E MP2

E SCF

(2) where the terms are, respectively, the kinetic energy, the electron–electron and electron– nuclear Coulomb energies, exact (HF) exchange, exchange and correlation DF contributions, opposite spin and same spin MP2 energies, and empirical dispersion correction. Depending on the values of the coefficients in eq 2, one may also obtain e.g. “pure” dispersion-corrected DFT (cX = cO = cS = 0), pure HF (cX = 1, cO = cS = cC = s6 = 0), pure MP2 (cX = cO = cS = 1, cC = s6 = 0), as well as simpler DHDFs like B2PLYP (cO = cS = 1 − cC , s6 = 0). These coefficients are assumed to be included in the definitions of E SCF and E MP2 . The DSD-DFT shielding tensor may also be split into SCF and MP2 contributions:

σ DSD-DFT = σ SCF + σ MP2 A A A

(3)

Note that the empirical dispersion correction does not depend on B or MA and therefore does not contribute to the shielding tensor. In ORCA σ SCF is calculated as described in ref A 37 with the only difference that the DFT correlation terms are scaled by cC . Because the usual RI-MP2 energy expression is not variational, it is convenient to obtain the MP2-specific density matrices required in eq 1 using a Lagrangian formulation, which 5



Page 6 of 51

includes the Hylleraas functional, 69 in the generator state matrix formulation, 70,71 and the Brillouin condition Fai = 0:

L=

X hD

Ei

E D X e ij∗ + D0 FT + e ij + Yj∗ Yi† T zai Fai Yj Yi,T T

(4)

ai

i≥j

e are the contravariant where the angle brackets denote a matrix trace, F is the Fock matrix, T MP2 coefficients or “amplitudes” (see eqs 9–10), D0 is the unrelaxed difference density matrix in the MO basis (see eqs 11–12), z is a vector of undetermined coefficients, and Y are the three-index two-electron repulsion integrals, symmetrically orthogonalized in the Coulomb metric (see eqs 6–8). Note that special care needs to be taken regarding complex conjugation, as the magnetic perturbation leads to imaginary matrix elements. After making the Lagrangian stationary with respect to the amplitudes, the unknown coefficients, and rotations among the molecular orbitals (MOs), the MP2 contribution to the shielding tensor takes a form analogous to eq 1: 1

T T B,M D T T M E B,MA B B,T B A σ MP2 = L = cD c h + c D + DU − U D c h A A

(5)

where D is the relaxed difference density matrix, c are the MO coefficients and UB are the orbital rotation coefficients.

2.2

Detailed Working Equations

The three-index RI integrals in eq 4 are expressed as 72–74

i = YaP

X

(ia | Q) V−1/2

QP

= YiPa∗

(6)

Q

Z (pq | Q) = Z VP Q =

−1 ϕp (r1 ) ϕq (r1 ) r12 ηQ (r2 ) dr1 dr2

(7)

−1 ηP (r1 ) r12 ηQ (r2 ) dr1 dr2

(8)

6


Page 7 of 51 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


Throughout this work the indices i, j, k denote occupied, a, b, c – unoccupied, and p, q, r – any p MOs, while P, Q denote auxiliary basis functions. General-index RI-integrals YqP are defined

analogously. Note that for the sake of simpler expressions a symmetric orthogonalization is implied in eq 6, while our implementation actually employs an asymmetric, Cholesky decomposition-based orthogonalization. The covariant and contravariant amplitudes are defined as ij Tab = (εi + εj − εa − εb )−1

X

b = ∆ij YiPa YjP ab

X

b YiPa YjP

(9)

P

P

ij ij ij Teab = (1 + δij )−1 2 (cS + cO ) Tab − 2cS Tba

(10)

where δij is the Kronecker delta symbol. Eqs 9–10 are obtained as in refs 70– 71, employing also the RI-MP2 73,74 and SCS-MP2 18 Ansätze. The unrelaxed MP2 density matrix only has occupied–occupied and virtual–virtual blocks given by

0 Dkj =−

X

D E e ij∗ Tki (1 + δij ) T

(11)

i 0 Dab =

X

e ji∗ Tij + T e ij∗ Tji T

i≥j

(12) ab

The relaxed MP2 density is defined as D = D0 + Z, 64,66,75 where the z-vector, cast into matrix form as Zai = Zia = 21 zai , is obtained by solving the z-vector equations 76 X

(εb − εj ) zbj +

R [Z]ai = 2Xbj

(13)

1 b YaP ΓjaP − R [D0 ]bj 2

(14)

ai

Xbj =

X

YiPj ΓibP −

iP

X aP

where the three-index density matrix is 52,75

ΓiaP =

X

ij j (1 + δij ) Teab YbP

jb

7


(15)


Page 8 of 51

and the Fock response operator for a given density matrix P is 52

R [P]pq =

X

Prs [(rs || pq) + (rs || qp)] + RXC [P]pq

(16)

rs

The antisymmetrized two-electron integrals (in Mulliken (1∗ 1 | 2∗ 2) notation) are defined as

(pq || rs) = 2 (pq | rs) − cX (ps | rq)

(17)

and may be approximated using the RIJK (RI for both Coulomb and exchange integrals), 77 RIJONX (RI for Coulomb only) or RIJCOSX approaches (RI for Coulomb and chain-ofspheres exchange) 78 of which the latter is the most efficient for large systems, due to its more favorable scaling behavior. 37,53,79 The contribution RXC , coming from the XC functional, is given in the Supporting Information and drops out in the case of pure MP2. Taking the magnetic field derivative of the Lagrangian stationarity conditions, we obtain the following expressions for the field-perturbed amplitudes, z-vector, and orbital rotation B . Analogous equations are derived in refs 1, 2, 61, and 13 coefficients Uai

( ij,B Tab =∆ij ab

i X h a,B b,B b YiP YjP + YiPa YjP P

+

X

ij B Tac Fbc + Tcbij Fac

B

c

) i Xh kj B ik B Tab Fkj + Tab Fki −

(18)

k

B (εb − εj ) zbj +

X

B R ZB ai = 2X bj

(19)

ai B

X bj =

i i Xh X h j,B b,B i,B b i i ΓaP YaP + ΓjaP YaP − Γi,B Y + Γ Y bP jP bP jP aP

iP

X 1X 1X 1 B − zaj Fab + zbi FjiB + Dpq (pq || bj)(B) + RXC(B) [D]bj 2 a 2 i 2 pq 1 1X B 1X B Upb R [D]pj + R [D]bp Upj − R UB D bj + R D0B bj − 2 2 p 2 p

8


(20)

Page 9 of 51 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


where

0B Dkj =

X

0B Dab =

X

(1 + δij )

D E D E e ij,B Tki − T e ij Tki,B T

(21)

i

e ji Tij,B + T e ij Tji,B − T e ji,B Tij − T e ij,B Tji T

i≥j

(22) ab

1 B B B = zai = −Zia Zai 2 X eij,B Y j + Teij Y j,B Γi,B = (1 + δ ) T ij aP ab bP bP ab

(23) (24)

jb

ij,B ij,B ij,B Teab = (1 + δij )−1 2 (cS + cO ) Tab − 2cS Tba X R [P]pq = cX Prs [(rp | qs) − (rq | ps)]

(25) (26)

rs

A superscript in parentheses denotes that only derivatives of the basis functions with respect to the perturbation are taken. The full derivative of the RI integrals is

p,B YqP =

X

cB µp =

X

V−1/2

Xh

QP

i B cµp cνq (µν | Q)B + cµp cB − c c (µν | Q) νq µp νq

(27)

µν

Q B cµq Uqp

(28)

q

Note that the auxiliary basis functions are independent of the magnetic field. As pointed out by Loibl et al., 13,80 using GIAO-type fitting functions instead would violate gauge invariance, because the explicit dependence on the gauge origin would not cancel out in the perturbed B three index integrals. The coefficients Uai are computed as solutions to the coupled perturbed

self-consistent field (CPSCF) equations, 81 which we have discussed in detail in ref 37. The B UijB and Uab blocks are also required and can be chosen in different ways, the most common

being 1 (B) UijB = − Sij 2 1 (B) B Uab = − Sab 2 9


(29) (30)


Page 10 of 51

However, a computationally more efficient alternative choice of UijB , resulting in perturbed canonical orbitals, 6,61,82 is presented in the Supporting Information. Expressions for the derivative integrals over GIAOs S(B) , hMA , hB,MA , (pq || jb)(B) and (µν | Q)B in eqs 29, 30, 5, 20, and 27 can be found in ref 37 (eqs 23, 25, 26, 22, and 45, respectively) where also the full derivative of the Fock matrix FB is discussed (eqs 10–16). Thorough discussions of the GIAO approach can also be found in refs 83 and 84. The XC contribution due to the GIAO Ansatz RXC(B) is discussed in the Supporting Information. The relaxed field-perturbed MP2 density is finally completed with the z-vector contribution DB = D0B + ZB

(31)

and the MP2 contribution to the shielding tensor is calculated according to eq 5. Additional care needs to be taken if the frozen core approximation is used. 13,75,85 The necessary changes are discussed in the Supporting Information. If an implicit solvation model is employed, an additional term enters the Fock response contributions – this is also given in the Supporting Information.

2.3

Implementation details

In ORCA the relaxed MP2 density D is evaluated as discussed in ref 52. The algorithm is similar to that of Weigend and Häser, 75 with some notable differences. For detailed descriptions of these implementations, the reader is referred to the original publications. i Here we briefly summarize the main steps: (1) YqP are calculated and stored on disk; (2) the

occupied MOs are split into batches according to available memory; (3) in each batch of i, e ij are calculated (eqs 9 and 10); (4) contributions integrals are read and amplitudes Tij and T to D0 (eqs 11–12), Γi (eq 15), and X (eq 14) are accumulated; (5) after the main loop, the response operator R [D0 ] (eq 16) is constructed and added to X; (6) the z-vector equations (eq 13) are solved and the relaxed density matrix D is stored on disk.

10


Page 11 of 51 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


The expressions for the magnetic field-perturbed quantities D0B , Γi,B , Tij,B (first two B

terms), and X (first four terms) – eqs 21–22, 24, 18, and 20, respectively – are analogous to the respective unperturbed quantities D0 , Γi , Tij , and X – eqs 11–12, 15, 9, and 14, respectively – except that the former contain six times as many contributions (recalling that the magnetic field has three components). Therefore, the algorithm to calculate these terms is also analogous to that for first DHDFT and RI-MP2 derivatives, discussed above and in refs 52 and 75, with six times as many operations (seven times as many in case the unperturbed quantity needs to be recalculated as well – this is the case for Tij and Γi ) and four times higher memory requirements (including the unperturbed quantities), leading to an expected increase of the computation cost by a factor of up to 20–30. Both unperturbed i,B i and YqP and perturbed three index integrals, YqP , are precalculated and stored on disk. The

similarities and differences between the algorithms for the calculation of D and DB are can be seen in Scheme 1. The last two terms in eq 18 require special treatment, which is discussed in the Supporting Information. As explained there, it is usually more efficient to store the unperturbed 0B (eq 21) where amplitudes on disk. Another potential bottleneck are the contributions to Dkj

all Tij and Tij,B are required for a given i. One approach, which is used for the equivalent 0 contributions to Dkj is to keep these amplitudes in memory. However, due to the higher

memory requirements, more batches would be needed, resulting in additional overhead. Alternatively, the amplitudes for the batch can be stored on disk and processed in a second loop over i. This is also produces overhead in the form of disk input/output (I/O) and e ij and T e ij,B and is therefore only necessary recalculation of the contravariant amplitudes T when memory is a limiting factor. B

As is the case for X, 52 the three external index contributions to X (first two terms in

11



i Make and store all YqP

Page 12 of 51

i and Y i,B Make and store all YqP qP

Make list of degenerate orbitals Li for all i B . Complete UB and FB Read D and all Uai

for batch I ∈ {occupied} do

for batch I ∈ {occupied} do

for i ∈ I do

for i ∈ I do

for j ∈ {occupied} do


e ij . Keep Tij in memory Make Tij and T

e ij . Keep Tij in memory Make Tij and T Y, Y B , and Tij contribs. to Tij,B for k ∈ Li ∪ Lj do Make/read Tkj and Tik Tkj and Tik contribs. to Tij,B end for(k) e ij,B . Keep Tij,B in memory Complete Tij,B and T

if j ≤ i then Tij

if j ≤ i then

e ij contribs. to E MP2 and D0 and T

e ij , Tij,B , and T e ij,B contribs. to D0B Tij , T

end if

end if

e ij contribs. to Γi T

e ij and T e ij,B contribs. to Γi and Γi,B T

end for (j)

end for (j)



for k ≤ j do

for k ≤ j do

e ij contribs. to D0 Tik and T

e ij , Tik,B , and T e ij,B contribs. to D0B Tik , T

end for (k)

end for (k)

end for (j)

end for (j)

3-internal Γi contribs. to X

3-internal Γi and Γi,B contribs. to X

end for (i)

end for (i)

3-external Γi contribs. to X

3-external Γi and Γi,B contribs. to X

B

B

end for (I)

end for (I)

Fock response contribs. to X

Fock response contribs. to X

Solve z-vector equations and complete D

Solve first order z-vector equations and complete DB

B

Scheme 1: Algorithms for the calculation of D (left) and DB (right), aligned to highlight the analogous steps.

12


Page 13 of 51 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


eq 20) are evaluated by partially transforming Γi and Γi,B to the AO basis,

ΓiνP =

X

Γi,B νP =

Xh

cνb ΓibP

(32)

b B i cνb Γi,B bP + cνb ΓbP

i

(33)

b

and contracting them with AO-basis three index integrals, generated on the fly: B

X ai ←

X µ

cB µa

X Pν

µ ΓiνP YνP +

X µ

cµa

i X h i,B µ µ,B ΓνP YνP + ΓiνP YνP

(34)

Pν

As in the DHDFT gradient implementation in ORCA, 52 but unlike the original algorithm by Weigend and Häser, 75 Γi and Γi,B for the whole batch are kept in memory. A final point to consider is the evaluation of two-electron integrals in the Fock response terms R [Z] and R ZB in left-hand sides (LHSs) and R [D0 ], (pq || bj)(B) , R D0B , R [D], and R UB D , in the right-hand sides (RHSs) of eqs 13 and 23, as well as those in the perturbed Fock matrix FB . The latter were discussed at length in ref 37 and that discussion largely applies also to the GIAO integrals (pq || bj)(B) . Theoretically it is most consistent to apply to these terms the same approximation (e.g. RIJK or RIJCOSX) that was used in the SCF procedure. However, the RI approximation to the exchange integrals offers no computational advantage when the latter are to be contracted with a density matrix defined in the entire MO space. Therefore, in our implementation the LHS terms are treated with the RIJONX approximation, which allows for better prescreening. An alternative is to store the required (ia | jb) and (ij | ab) integrals on disk, rather than recalculate them at each CPSCF iteration. This is most efficiently done via an RI transformation during the first run through the RI-MP2 program, where the integrals (ia | P ) and (ij | P ) are already available, although (ab | P ) are normally not. Note, however, that the RI-MP2 auxiliary basis set (denoted AuxC) used, rather than the one used for the SCF (denoted AuxJ). In Section 4 we compare the efficiency of both of these approaches.

13



3

Page 14 of 51

Benchmark calculations

In a previous publication we have assembled a benchmark set of molecules for NMR CSC calculations, comprising a total of 34 1 H,

13

C,

15

N,

17

O,

19

F, and

31

P shielding values. 37

The molecules were chosen such that a broad range of values is covered for each nucleus, experimental gas-phase data are available for most CSCs, and the systems are small enough to allow for high-quality CCSD(T) reference calculations. The benchmark set and the reference data are given in Table 1.

14


Page 15 of 51 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 1: Benchmark set from ref 37 and empirical equilibrium values, obtained from previously published gas phase experimental data and theoretically estimated vibrational corrections. All values are in ppm. Molecule 1H

13 C

15 N

17 O

19 F

31 P

PH3 HF H2 O NH3 CH4 (CH3 )2 CO furan (at C2/5) furan (at C3/4) (CH3 )2 CO CO CF4 furan (C2/5) furan (C3/4) (CH3 )2 CO CH4 PN N2 NNO NNO NH3 OF2 (CH3 )2 CO CO furan NNO H2 O F2 OF2 PF3 CF4 HF PN PF3 PH3

CCSD(T)/pS4 29.46 28.82 30.65 31.44 31.39 29.53 24.03 25.02 −10.84 2.56 65.96 47.36 81.67 162.88 199.39 −344.71 −61.16 11.74 106.22 270.40 −446.32 −297.91 −55.42 64.82 198.77 337.63 −192.76 −24.28 231.81 267.58 419.91 51.61 224.80 604.51

Exp.

Vib. Corr.

27.89 28.64 30.05 30.68 30.80

0.33a 0.52a 0.61a 0.63a

Emp. Eq. 28.97 30.57 31.29 31.43

−13.2 0.9 64.4

0.8b 1.0b 1.4b

−12.4 1.9 65.8

157.9 195.0 −349 −61.6 11.3 99.5 264.5 −493.6 −309.1 −62.7 50.3 181.0 323.5 −233.2 −59.7 230.8 258.6 409.6 53.0 222.7 594.5

3.3b 2.6b 5.3a 3.3c 3.1c 6.8c 6.8c 34.0d 3.8d 3.8d 11.6d 8.2d 9.6d 23.6e 23.7e

161.2 197.6 −343.7 −58.3 14.4 106.3 271.3 −459.6 −305.3 −58.9 61.9 189.2 333.1 −209.6 −36.0

7.3e 8.6e 4.4f 2.3f 9.5f

265.9 418.2 57.4 225.0 604.0

a

Exp. Ref. 86 87, 88 87 87 87, 89

90, 91 91 90, 91

90, 90, 92 93 93 93 93 94 94 94 94 94 94 91, 91, 91, 91, 91 92 98 98

91 91

95 96, 97 95 95

Geometry: CCSD(T)/cc-pVTZ; ZPV correction: B3LYP/aug-cc-pCVTZ. 36 Geometry: CCSD(T)/cc-pVTZ; ZPV correction: MP2/cc-pVTZ (force field) + MP2/qz2p (NMR shielding). 7 c Geometry: CCSD(T)/cc-pVQZ; ZPV correction: CCSD(T)/cc-pVTZ (force field) + CCSD(T)/qz2p (NMR shielding). 10 d Geometry: CCSD(T)/cc-pVTZ; ZPV and thermal corrections: MP2/cc-pVTZ (force field) + MP2/qz2p (NMR shielding). 9 e Geometry: CCSD(T)/cc-pVTZ; ZPV and thermal corrections: MP2/cc-pVTZ (force field) + MP2/qz2p (NMR shielding). 8 f Geometry: CCSD(T)/cc-pVQZ; ZPV correction: CCSD(T)/cc-pVTZ (force field) + CCSD(T)/qz3d1f (NMR shielding). 10 b

15



Page 16 of 51

In the following we will assess the accuracy of RI-MP2 and DHDFT CSC calculations by separately examining errors inherent to the method itself (Section 3.2), as well as those due to basis set incompleteness in the orbital (Section 3.3) and auxiliary (Section 3.1) basis sets, the frozen core approximation (Section 3.4), and two-electron integral approximations in the Fock matrix (Section 3.5). All NMR CSC calculations in this work are performed using the pcSseg-n (n = 0−4) basis sets (further denoted pSn for the sake of brevity), 99 a segmented contracted version of the pcS-n basis sets. 100 Both basis set families were optimized for DFT-level NMR calculations but the latter have also been shown to converge rapidly towards the complete basis set (CBS) limit at the MP2-level. 11,101 Canonical MP2 results used as a reference in Section 3.1 were calculated using the CFOUR program. 1,2,102 All other calculations were performed with a development version of ORCA 4.0. 62,63 The SCF equations were converged to a maximum change of density matrix elements of 10−8 . CPSCF and z-vector equations were converged to 10−10 . We use several statistical measures of the deviation of CSCs or relative chemical shifts (RCSs), calculated according to a given protocol (method, basis, etc.) with respect to a given reference (e.g. CCSD(T) data). The mean error (ME), mean absolute error (MAE), mean relative error (MRE), mean absolute relative error (MARE), and maximum errors and relative errors (MaxAE and MaxARE, respectively) are defined as: ref ˜ i = σ i − σi ∆σ |σiref | N 1 X˜ ∆σi MREσ = N i=1

∆σi = σi − σiref N 1 X MEσ = ∆σi N i=1

N N 1 X 1 X ˜ MAEσ = |∆σi | MAREσ = ∆σi N i=1 N i=1 h i ˜ MaxAEσ = max |∆σi | MaxAREσ = max ∆σi

(35)

The same quantities (with an index δ) are defined analogously for the RCSs, which are given 16


Page 17 of 51 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


by δi = σstd − σi , where σstd is the CSC of the given nucleus in a “standard” compound, calculated at the same level of theory. We use CH4 as a standard for 15

N, H2 O for

17

O, HF for

19

F, and PH3 for

31

13

C and 1 H, NH3 for

P.

In the box-and-whiskers plots used throughout the rest of this work boxes show the interquartile range (IQRE or IQRRE), i.e. the range of errors (or relative errors), excluding the top and bottom 25 %; whiskers show the minimum and maximum errors (MinE and MaxE) or relative errors (MinRE and MaxRE); and lines show the median errors (MedE) or relative errors (MedRE) and the ME or MRE. Note that wherever relative errors are calculated, those data with small reference values (less than 10 ppm for CSCs and 1 ppm for RCSs) were excluded from the analysis because they disproportionately skew the results. The excluded nuclei are listed in the figure captions.

3.1

Auxiliary Basis Sets

In order to reliably evaluate the accuracy of DHDFT CSC calculations, we must select suitable AuxC basis sets to make sure that the errors due to the RI-MP2 approximation are negligible. For our benchmarking purposes we have decided on a very conservative threshold of 0.01 and 0.1 % for the MAREσ and MaxAREσ , respectively. In routine applications less stringent thresholds and therefore smaller AuxC basis sets may be sufficient. The chosen pS2, pS3, and pS4 orbital basis sets (OBSs) are of triple-, quadruple-, and quintuple-zeta quality, respectively, and contain tight functions to better describe the core region. Regrettably, no auxiliary basis sets have been specifically optimized for these OBSs. Hence, we have chosen the AuxC basis sets, here denoted cwnC (n = 3, 4, 5), 103 optimized for Peterson’s cc-pwCVXZ (X = T, Q, 5) basis sets, 104 which also contain additional tight functions. Alternatively, we employ the AutoAux scheme to generate large fitting basis sets for each OBS. 105 As well as the one made with the default settings, denoted AA, we generate a near-complete AuxC basis (using the ORCA keywords AutoAuxSize=3 and AutoAuxLMax=true), denoted AA3l, which encompasses the full product space of the 17



atomic OBS. This allows us to validate the correctness of our implementation at the limit of very large AuxC basis sets by comparing our RI-MP2 CSC results to canonical MP2 values. The results of the comparison are presented in Figure 1. The first thing to note is that for very large AuxC sets, the errors are vanishingly small, i.e. the RI-MP2 CSCs converge to the canonical MP2 values. Second, the pS2/cw3C and pS3/cw4C combinations result in errors slightly above our chosen thresholds (MAREσ of 0.03 and 0.02 % and MaxAREσ of 0.17 and 0.13 %, respectively). As mentioned above, these errors are quite acceptable for most applications, especially when compared to the other sources of error (see below). However, in order to accurately estimate these other errors, we have chosen the pS2/cw4C, pS3/cw5C, and pS4/cw5C combinations as default for the rest of this work. It is also worth noting that although not parametrized for CSC calculations, the AutoAux scheme produces sufficiently accurate AuxC basis sets for the purpose, albeit around 1.5 times larger than the cwnC sets of similar quality. In the following sections the RI-MP2 approximation is always used and therefore the RIprefix will be skipped. It is, however, implied by the use of AuxC basis sets denoted in the figures.

3.2

Method Error

In this section we evaluate the deviations, inherent to the methods used, in the calculated CSCs from the CCSD(T) and empirical equilibrium reference data presented in Table 1. In particular, we compare MP2, SCS-MP2, and two DHDFs: the original, and still widely used, B2PLYP and the more recent DSD-PBEP86, using the parameters optimized for DSDPBEP86-D3BJ from ref 49. Note however, that the dispersion correction does not contribute to the CSC. We also include results from our previous study, 37 namely at the HF, B3LYP, TPSS, and M06L levels. The latter two functionals showed the smallest overall deviations in our tests at the SCF (i.e. HF or DFT) level. The pS4 basis set (and the cw5C AuxC set, where applicable) is used for all methods, as it is assumed to be sufficiently close to the CBS 18


Page 18 of 51

Page 19 of 51 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


All nuclei (N = 31) pS2/cw3C (3.O) pS2/cw4C (5.O) pS2/cw5C (5.1) pS2/AA

(7.6)

pS2/AA3l(1O.4) pS3/cw4C (2.5) pS3/cw5C (3.7) pS3/AA

(4.6)

pS3/AA3l (9.6) pS4/cw5C (2.3) pS4/AA

(3.9)

pS4/AA3l (8.5)

0.15 0.10 0.05 0.00 Relative errors in vs canonical MP2 / %

Figure 1: Relative deviations of chemical shielding constants (%), calculated at the RI-MP2 level using different OBS/AuxC combinations, from canonical MP2 results with the same OBS. Excluded: CO, CH3 COCH3 , and F2 O. The average AuxC/OBS size ratios are given in parentheses. Boxes show the IQRREσ , whiskers show the MinREσ and MaxREσ , orange lines show the MedREσ , and green dashed lines show the MREσ . limit. The errors in the CSCs for different groups of nuclei are shown in Figure 2. Because the CCSD(T) results are so close to the empirical equilibrium data, similar conclusions can be drawn from either reference. Therefore, only the CCSD(T) reference is discussed here and a comparison to the empirical equilibrium data is given in the Supporting Information. Out of the SCF methods, the error distribution of HF is the widest, while those of the DFT methods are narrower. The different DFs also show a systematic increase of the CSCs in the order B3LYP < TPSS < M06L. Therefore, overall the M06 functional is most accurate for

13

C,

15

N,

17

O,

19

F, and

31

P CSCs, followed by TPSS and B3LYP, while the order is

reversed for 1 H shielding. MP2 shows a significant improvement over HF and is mostly better than DFT – the systematic deviation is consistently small for all elements and the error distribution is narrower, with the exception of M06L for heteroatom CSCs. Note the particularly good performance of MP2 for carbon CSCs, which has been discussed before. 7,11 SCS-MP2 improves the results for 1 H CSCs but does not appear significantly better for the other nuclei. It is interesting to compare B3LYP, MP2, and B2PLYP: the former two

19



show systematic deviations in the calculated CSCs with opposite signs, while the latter is somewhere in between. While this is not surprising, one can suppose that optimizing the cC parameter can further reduce the systematic error in B2PLYP CSCs. The resulting functional would be analogous to B2K-PLYP and B2T-PLYP, 106 optimized for kinetics and thermochemistry, respectively. At the very least, an analysis similar to that performed in the optimization of the B2GP-PLYP functional, 107 could provide insight into the dependence of the calculated CSCs on the parameters cC and cX (as well as cOS and cSS ). However, such an analysis is outside the scope of the current work. 108 Finally, we note the outstanding performance of the DSD-PBEP86 functional: it produces the narrowest error distributions of all methods tested here and while there is a significant systematic deviation in the carbon CSCs, this is expected to cancel out in RCSs. The effect of error cancellation when calculating RCSs – which are the experimentally measured quantities – can be seen in Figure 3, which shows the relative errors in RCSs for all nuclei. There is still a significant systematic deviation in the B3LYP data, while this is smaller for the two meta-GGA DFs. MP2 is clearly superior to the SCF methods, showing a narrow error distribution and no significant systematic deviation. SCS-MP2 has an even narrower error distribution, however both MP2 and SCS-MP2 have a MaxAREδ > 12 %. Maurer and Ochsenfeld have shown that the accuracy of SCS-MP2 for CSCs can be substantially improved by optimizing the cSS and cOS coefficients. 21 However, the vastly different optimal values for each combination of basis set and target element, makes it impossible to attach any physical interpretation to these parameters, essentially turning SCS-MP2 into a rather expensive semi-empirical method. Over the whole set of RCSs, B2PLYP shows a narrower error distribution that either MP2 or SCS-MP2. However, as for B3LYP, there is a residual systematic error, which does not cancel out. On the other hand, DSD-PBEP86 exhibits the best performance overall with the narrowest error distribution and almost no systematic deviation. The latter may be partly due to the fact that it includes a smaller proportion of DFT correlation (cC = 0.44), compared to B2PLYP (cC = 0.73), although the

20


Page 20 of 51

Page 21 of 51 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


different XC functional is certainly also a factor. In summary, the most accurate meta-GGA functionals explored here – TPSS and M06L – result in MAREδ of 6.4 and 5.4 %, respectively. SCS-MP2 shows a slight improvement over MP2 with a MAREδ of 3.9 vs 4.1 %. By that measure, B2PLYP is not better with MAREδ = 4.3 %. However, DSD-PBEP86 shows remarkable accuracy with MAREδ = 1.9 %, enabling for example the confident assignment of carbon shifts only a few ppm apart. H (N = 8)

C (N = 7)

N,O,F,P (N = 19)

HF/pS4 B3LYP/pS4 TPSS/pS4 M06L/pS4 MP2/pS4/cw5C SCS-MP2/pS4/cw5C B2PLYP/pS4/cw5C DSD-PBEP86/pS4/cw5C

0.5

0.0

0.5

1.0

30 20 10 0 10 Errors in vs CCSD(T)/pS4 / ppm

150

100

50

0

50

Figure 2: Deviations of chemical shielding constants (ppm) for groups of nuclei, calculated using different methods, from CCSD(T)/pS4 results. The number of nuclei in each group is given in parentheses. Boxes show the IQREσ , whiskers show the MinEσ and MaxEσ , orange lines show the MedEσ , and green dashed lines show the MEσ . All nuclei (N = 26) HF/pS4 B3LYP/pS4 TPSS/pS4 M06L/pS4 MP2/pS4/cw5C SCS-MP2/pS4/cw5C B2PLYP/pS4/cw5C DSD-PBEP86/pS4/cw5C

20 0 20 40 Relative errors in vs CCSD(T)/pS4 / %

Figure 3: Relative deviations of chemical shifts (%), calculated using different methods, from CCSD(T)/pS4 results. The number of data points is given in parentheses. Excluded: NH3 and H2 O. Boxes show the IQRREδ , whiskers show the MinREδ and MaxREδ , orange lines show the MedREδ , and green dashed lines show the MREδ .

21


100


3.3

Basis Set Convergence

We have previously shown that for SCF-level methods the basis set incompleteness error for pS2 is already an order of magnitude smaller than the inherent error of the method. 37 However, a slower convergence can be expected for correlated wave function-based methods. Likewise, DHDFs should exhibit a stronger dependence on the basis set size than other DFs, as is the case for energies. 50 In this section we compare the basis set convergence of DSDPBEP86 for CSCs and RCSs and compare it to that of HF, TPSS and MP2, using the pS4 basis set as a reference. The deviations of CSCs for groups of nuclei are presented in Figure 4. Note that the SCF-level results are not quite converged at the pS2 level. However, the deviations are negligible compared to the errors due to the methods. On the other hand, the basis set errors for MP2/pS2 are more than twice as large as for HF/pS2 and fall within the same order of magnitude as the method error, which is due also to the superior accuracy of MP2. As expected, the basis set errors for DSD-PBEP86 are between those for HF and MP2 and also within the same order of magnitude as the method error. Because of the significant systematic deviations in the CSCs, it can be expected that the basis set incompleteness errors for RCSs will be smaller. The latter are shown in Figure 5. Apparently, the systematic errors do cancel out for the SCF methods but not entirely for the perturbative approaches. Therefore it cannot be claimed that either MP2 or DSD-PBEP86 RCSs are sufficiently converged at the triple-zeta (pS2) level. However, due to fortuitous error compensation, the deviations for both approaches with respect to the CCSD(T)/pS4 reference are actually smaller for pS2 than for pS3 with MAREδ of 3.70, 1.41, 4.00, and 1.84 % for MP2/pS2/cw4C, DSD-PBEP86/pS2/cw4C, MP2/pS3/cw5C, and DSD-PBEP86/pS3/cw5C, respectively. This error cancellation has been noted before – e.g. in refs 109, 7, and 11 – and while it is not necessarily reliable, it provides some justification for the use of the smaller pS2 basis set. Note, however, that because SCS-MP2 already underestimates the RCSs at the pS4 level, a reduction of the basis set size to pS2 leads to 22


Page 22 of 51

Page 23 of 51 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


an increase in the total error versus CCSD(T)/pS4 with MAREδ = 5.2 %. H (N = 8)

C (N = 7)

N,O,F,P (N = 19)

HF/pS2 HF/pS3 TPSS/pS2 TPSS/pS3 MP2/pS2/cw4C MP2/pS3/cw5C DSD-PBEP86/pS2/cw4C DSD-PBEP86/pS3/cw5C

0.1

0.0

0.1

0.2

0.3

0

1 2 3 4 Errors in vs pS4 / ppm

10

0

10

Figure 4: Deviations of chemical shielding constants (ppm) for groups of nuclei, calculated using different methods and basis sets, from pS4 results for the same method. The number of nuclei in each group is given in parentheses. Boxes show the IQREσ , whiskers show the MinEσ and MaxEσ , orange lines show the MedEσ , and green dashed lines show the MEσ . All nuclei (N = 26) HF/pS2 HF/pS3 TPSS/pS2 TPSS/pS3 MP2/pS2/cw4C MP2/pS3/cw5C DSD-PBEP86/pS2/cw4C DSD-PBEP86/pS3/cw5C

7.5

5.0 2.5 0.0 2.5 5.0 Relative errors in vs pS4 / %

Figure 5: Relative deviations of chemical shifts (%), calculated using different methods and basis sets, from pS4 results for the same method. The number of data points is given in parentheses. Excluded: NH3 and H2 O. Boxes show the IQRREδ , whiskers show the MinREδ and MaxREδ , orange lines show the MedREδ , and green dashed lines show the MREδ .

3.4

Frozen Core Approximation

The usual justification for using the frozen core (FC) approximation is that the core electron contribution to the correlation energy does not change significantly during chemical bonding and therefore cancels out in the calculation of relative energies (between systems of the same chemical composition). However, this argumentation does not necessary apply to calculated 23



Page 24 of 51

NMR properties for several reasons. First, the chemical shielding tensor is inherently sensitive to changes in the electron density near the nucleus. Second, in the calculation of RCSs the reference system can, and often does, have completely different chemical structure from the studied system and may therefore exhibit very different core correlation effects. While the frozen core approximation is employed for CSC calculations in the local RI-MP2 implementation of Loibl and Schütz, 13 little justification is given for its use. In this section we attempt to gauge the size of the errors introduced by the FC approximation. It is important to note that our test set does not include heavy nuclei or other systems for which core correlation is known to be important. Hence, the FC errors reported here are likely underestimated and should only be used as a guideline. H (N = 8)

C,N,O,F (N = 23)

P (N = 3)

MP2/pS2/cw4C MP2/pS3/cw5C MP2/pS4/cw5C DSD-PBEP86/pS2/cw4C DSD-PBEP86/pS3/cw5C DSD-PBEP86/pS4/cw5C

0.00

0.02

0.04

0.06

0

5 10 15 0 Errors in vs all-electron / ppm

10

20

Figure 6: Deviations of chemical shielding constants (ppm) for groups of nuclei, calculated with the frozen core approximation, from all–electron results for the same method. The number of nuclei in each group is given in parentheses. Boxes show the IQREσ , whiskers show the MinEσ and MaxEσ , orange lines show the MedEσ , and green dashed lines show the MEσ . Figure 6 shows the FC error in the CSCs for MP2 and DSD-PBEP86. Note that unlike other plots in this work, the nuclei are grouped by rows in the periodic table to highlight the larger errors for heavier elements. In addition, hydrogen only has valence electrons, hence for 1 H shielding the FC error is due to the influence of the core electrons from neighboring atoms. Naturally, the deviations are smaller for DSD-PBEP86 because of the overall scaling of the MP2 contribution, hence we will only discuss the MP2 results. One important observation is that all FC errors are positive, i.e. the CSCs are always overestimated in the FC approximation. Such a systematic deviation would cancel out in RCSs, which is a point 24


30

Page 25 of 51 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


in favor of the approximation. However, note also the wide range of errors: depending on the chosen reference, the expected error cancellation might be very different. There is also a clear basis set dependence of the FC error. This is due to the inability of the smaller basis sets, especially pS2, to adequately describe the core region. The pS4 values should therefore be considered the mosts representative. All nuclei (N = 26) MP2/pS2/cw4C MP2/pS3/cw5C MP2/pS4/cw5C DSD-PBEP86/pS2/cw4C DSD-PBEP86/pS3/cw5C DSD-PBEP86/pS4/cw5C

4 3 2 1 0 Relative errors in vs all-electron / %

Figure 7: Relative deviations of chemical shifts (%), calculated with the frozen core approximation, from all–electron results for the same method. The number of data points is given in parentheses. Excluded: NH3 and H2 O. Boxes show the IQRREδ , whiskers show the MinREδ and MaxREδ , orange lines show the MedREδ , and green dashed lines show the MREδ .

25



Page 26 of 51

Table 2: Analysis of the frozen core (FC) errors in the CSC and RCS of 3 1P in PN. All values, except in the last two columns, are in ppm.

σ(PH3 ) σ(PN) δ(PN)

HF

MP2

FC-MP2

MP2 − HF

FC-MP2 − HF

∆FC

∆FC /% MP2

∆FC /% MP2−HF

584 -110 694

609 107 502

616 135 480

26 217 -192

32 246 -213

7 28 -22

1 27 -4

26 13 11

26


Page 27 of 51 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


The RCS data in Figure 5 show that in fact the systematic deviations do not cancel out, at least using the chosen references. This is perhaps best understood by an example: the largest frozen core error in our test set – 28 ppm – is for the phosphorus CSC in PN, which is 27 % of the total shielding constant and 13 % of the MP2 contribution (see Table 2). The FC error for the reference nucleus in PF3 is only 7 ppm (26 % of the MP2 contribution). Therefore the final absolute error in the RCS is 22 ppm, which is only 4 % of the total shift value but 11 % of the MP2 contribution. While an error of 4 % might be considered acceptable in some cases, it is of the same order of magnitude as the method error, thereby increasing the deviation with respect to CCSD(T)/pS4 from 9 to 13 %. Moreover, it is apparent that one cannot rely on error cancellation in the RCSs, as the absolute FC error is not much smaller and the relative error decreases only because the denominator is larger for the RCS than for the CSC (694 vs 107 ppm). While this is the most extreme example in our test set, similar observations can be made for most of the other nuclei. Therefore, we conclude that use of the FC approximation is not justified for NMR RCS calculations, as the resultant errors – MAREδ = 1.268 and 0.457 % for MP2 and DSD-PBEP86, respectively – are of the same order of magnitude as the inherent accuracies of the methods. A final point to make is that while the FC errors are smaller for lighter elements, so are the computational savings gained by the approximation. Conversely, calculations on heavier nuclei could benefit more from freezing the core electrons, but the errors thus introduced would also be greater.

3.5

RIJK/RIJCOSX for the Fock Matrix

In the previous sections, no approximation was used for two-electron integral contributions to the Fock matrix (including the GIAO terms). In this section we apply the RIJK and RIJCOSX approximations to these terms and assess the additional errors thus introduced. Ideally, the latter should be at least an order of magnitude below the method and basis set errors. 27



Page 28 of 51

In our previous work we discussed the theory, implementation, accuracy, and efficiency of these approximations both for regular two-electron integrals, and derivative integrals over GIAOs. 37 We concluded that Weigend’s “universal” Coulomb- and exchange-fitting basis set (denoted def2-JK) 110 is sufficiently flexible to be used with the pS2 and pS3 OBSs for HF and DFT CSC calculations employing either RIJK or RIJCOSX. 37 We also proposed two sets of COSX grid parameters denoted COSX-S and COSX-L. The former uses the g2 angular grid with a radial integration accuracy ε = 4.0 for the SCF and CPSCF RHS contributions, while the latter uses the larger g3 angular grid with ε = 4.0. In both cases a smaller g1 angular grid with ε = 3.3 is used for the CPSCF LHS. For more information about these settings, see the discussion in ref 37. Here we extend these definitions such that the grid used for the CPSCF RHS is also used for the Fock response terms in the RHS of the z-vector equations, while the grid used for the CPSCF LHS is also used for the LHS of the z-vector equations. H (N = 8)

C (N = 7)

N,O,F,P (N = 19)

MP2/pS2/cw4C/RIJCOSX-S MP2/pS2/cw4C/RIJCOSX-L MP2/pS2/cw4C/RIJCOSX-XL MP2/pS2/cw4C/RIJK MP2/pS3/cw5C/RIJCOSX-S MP2/pS3/cw5C/RIJCOSX-L MP2/pS3/cw5C/RIJCOSX-XL MP2/pS3/cw5C/RIJK DSD-PBEP86/pS2/cw4C/RIJCOSX-S DSD-PBEP86/pS2/cw4C/RIJCOSX-L DSD-PBEP86/pS2/cw4C/RIJCOSX-XL DSD-PBEP86/pS2/cw4C/RIJK DSD-PBEP86/pS3/cw5C/RIJCOSX-S DSD-PBEP86/pS3/cw5C/RIJCOSX-L DSD-PBEP86/pS3/cw5C/RIJCOSX-XL DSD-PBEP86/pS3/cw5C/RIJK

0.02

0.01

0.00

0.01

0.0 0.2 0.4 3 Errors in vs exact Fock / ppm

2

1

0

Figure 8: Deviations of chemical shielding constants (ppm) for groups of nuclei, calculated using the RIJK and RIJCOSX approximations for the two-electron Fock matrix contributions, from results using exact two-electron integrals for these terms. The number of nuclei in each group is given in parentheses. Boxes show the IQREσ , whiskers show the MinEσ and MaxEσ , orange lines show the MedEσ , and green dashed lines show the MEσ .

28


1

2

Page 29 of 51 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


The errors for CSCs are shown in Figure 8. The first thing to notice is that the RIJK errors are very small, which confirms that the def2-JK basis set is large enough to be used with pS2 and pS3. On the other hand, RIJCOSX-S errors are an order of magnitude larger. RIJCOSX-L errors are smaller with pS2 but not with pS3, which is unexpected, as we had previously observed that the COSX errors (with a given grid setting) decrease with increasing basis set size. 37 Additional testing revealed that a large part to the error is due to the smaller grid used in the z-vector equations LHSs. Therefore, we propose a third set of grid parameters, denoted RIJCOSX-XL, whereby g3/ε = 4.0 is used for the CPSCF and z-vector equations RHSs, as in RIJCOSX-L, and g2/ε = 4.0 is used for the LHSs. Using these settings, the RIJCOSX errors are roughly of the same magnitude as the RIJK errors. It should be stressed however, that for pS2 the RIJCOSX-S errors are already an order of magnitude below the basis set error and for pS3 the RIJCOSX-L errors are several times smaller then the basis set error, albeit not a whole order of magnitude. Therefore, the RIJCOSX-L settings should be quite sufficient for regular applications and the RIJCOSX-XL settings need only be used when very precise results are required. Because the errors in the CSCs are rather unsystematic, they do not cancel out in RCSs (Figure 9). Hence, the conclusion above applies here as well: RIJCOSX-S may be used with the pS2 basis set, as the resultant errors are far below the basis set error (MAREδ = 0.165 and 0.115 % for MP2 and DSD-PBEP86, respectively), while RIJCOSX-L is more robust overall with MAREδ of 0.039, 0.030, 0.104 and 0.042 % for MP2/pS2/cw4C, DSDPBEP86/pS2/cw4C, MP2/pS3/cw5C, and DSD-PBEP86/pS3/cw5C, respectively. RIJK results in negligible errors (MAREδ < 0.03 % in all cases) and is therefore the preferred approximation for smaller systems, while RIJCOSX should be used for larger calculations due to its more favorable scaling behavior.

3.6

Summary

29



Page 30 of 51

All nuclei (N = 26) MP2/pS2/cw4C/RIJCOSX-S MP2/pS2/cw4C/RIJCOSX-L MP2/pS2/cw4C/RIJCOSX-XL MP2/pS2/cw4C/RIJK MP2/pS3/cw5C/RIJCOSX-S MP2/pS3/cw5C/RIJCOSX-L MP2/pS3/cw5C/RIJCOSX-XL MP2/pS3/cw5C/RIJK DSD-PBEP86/pS2/cw4C/RIJCOSX-S DSD-PBEP86/pS2/cw4C/RIJCOSX-L DSD-PBEP86/pS2/cw4C/RIJCOSX-XL DSD-PBEP86/pS2/cw4C/RIJK DSD-PBEP86/pS3/cw5C/RIJCOSX-S DSD-PBEP86/pS3/cw5C/RIJCOSX-L DSD-PBEP86/pS3/cw5C/RIJCOSX-XL DSD-PBEP86/pS3/cw5C/RIJK

0.8

0.6

0.4 0.2 0.0 0.2 0.4 0.6 Relative errors in vs exact Fock / %

0.8

Figure 9: Relative deviations of chemical shifts (%), calculated using the RIJK and RIJCOSX approximations for the two-electron Fock matrix contributions, from results using exact two-electron integrals for these terms. The number of data points is given in parentheses. Excluded: NH3 and H2 O. Boxes show the IQRREδ , whiskers show the MinREδ and MaxREδ , orange lines show the MedREδ , and green dashed lines show the MREδ .

30


Page 31 of 51 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 3: MAREδ (%) due to different sources of error (excluding H shifts in H2 O and NH3 ). Source of error

Basisa

Methodb

pS4

Frozen

corec

setd

HF

B3LYP

TPSS

M06L

MP2

SCS-MP2

B2PLYP

DSD-PBEP86

10.762

9.048

6.362

5.420

4.118

3.905

4.293

1.914

pS4

1.268

0.457

pS2 pS3

0.766 0.105

0.836 0.133

0.838 0.606

1.860 0.359

1.129 0.228

RIJKe

pS2 pS3

0.023 0.025

0.013 0.019

0.023g 0.025g

0.019 0.026

0.022 0.020

RIJCOSX-Sef

pS2 pS3

0.159 0.082

0.042 0.035

0.165 0.158

0.115 0.071

RIJCOSX-Lef

pS2 pS3

0.033 0.025

0.016 0.018

0.039 0.104

0.030 0.042

RIJCOSX-XLef

pS2 pS3

0.035 0.054

0.031 0.031

RI-MP2

pS2/cw3C pS2/cw4C

0.016h 0.005h

0.009i 0.004i

pS3/cw4C pS3/cw5C

0.025h 0.002h

0.011i 0.001i

pS4/cw5C

0.009h

Basis

DFT DFT DFT DFT DFT DFT

grid grid grid grid grid grid

Totalbk Totalbk

1j 2j 3j 4j 5j 6j

0.275 0.218 0.034 0.028 0.012 0.009

pS2 pS2 pS2 pS2 pS2 pS2 pS2/cw3C pS3/cw4C

11.055 10.832

9.123 9.048

6.271 6.075

0.086 0.076 0.010 0.002 0.001 0.001 3.962 3.385

3.695 4.064

a

5.175 4.014

4.015 4.353

1.420 1.863

The pS2/cw4C, pS3/cw5C, and pS4/cw5C OBS/AuxC combinations were used for the RI-MP2 approximation (where applicable), except where explicitly specified. b Vs CCSD(T)/pS4. c Vs all-electron calculations with the same basis set. d Vs pS4. e Using the def2-JK AuxJ basis, vs the same method/basis with no approximation in the two-electron Fock contributions. f See text for grid settings used. g No exact exchange - RI used only for Coulomb terms. h Vs canonical MP2 results obtained using CFOUR. i Vs AA3l near-complete AuxC basis. j Vs grid 7. 111 k All-electron calculations using RIJCOSX-L (just RIJ for TPSS and M06L) and DFT grid 4.

31



The magnitude of errors, coming from different sources, in the calculated RCSs can be easily compared in Table 3. Largest are the deviations of different methods from the CCSD(T)/pS4 reference: close to the CBS limit the ranking is (from most to least accurate): DSD-PBEP86, SCS-MP2, MP2, B2PLYP, M06L, TPSS, B3LYP, and HF. When the RI-MP2 and RIJCOSX errors are compared to those from the method and basis set, pS2/cw3C, pS3/cw4C, and RIJCOSX-L can be seen as a good balance of cost and accuracy. Although not discussed above, the effect of the DF integration grid was also studied: while the default grid for energy calculations (grid 2) 111 may be somewhat unreliable for RCSs, the values obtained with grids 3 and 4 are an order of magnitude more accurate, hence these grid settings are recommended. Finally, the last two rows of Table 3 show the total error versus CCSD(T)/pS4, when all the relevant approximations (RIJCOSX, RI-MP2, smaller basis sets and DFT grid, but not FC) have been applied. Note that some methods – M06L, MP2, B2PLYP, and DSD-PBEP86 – benefit from error compensation when a smaller basis set is used, while others, such as SCS-MP2, suffer from an accumulation of errors. Therefore, the final ranking at the pS2 level is: DSD-PBEP86, MP2, M06L, B2PLYP, SCS-MP2, TPSS, B3LYP, and HF.

4

Computational cost

In the previous section we have shown the superior accuracy of DHDFs (DSD-PBEP86 in particular) for the computation of NMR RCSs, compared to SCF-level methods and MP2. However, it must be stressed that, although applicable to much larger systems than coupled cluster would be feasible for, these calculations are significantly more time-consuming than the hybrid DFT equivalents, and even more so than pure DFT CSC calculations, where even the iterative solution of the CPSCF equations is not needed. Therefore, in this section we evaluate the performance of our implementation for larger systems. The (all-electron) DSD-PBEP86/pS2/cw3C level of theory was used throughout this section.

32


Page 32 of 51

Page 33 of 51

RIJK(RIJONX) RIJK(RITrafo) RIJCOSX-L

250 Wall-clock time / minutes

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


200 150 100 50 0

2

4

6 8 10 12 Number of carbon atoms

14

Figure 10: Scaling with system size of the total DSD-PBEP86/pS2/cw3C computation time for linear alkane chains (Cn H(2n+2) ) using different approximations for the Fock matrix contributions (the approximations in parentheses are used for the LHS of the CPSCF equations – see text for details). The calculations were performed on 8 Intel Xeon E7-8837 2.67 GHz cores with 8 GB RAM per core. The computational effort is dominated by the evaluation of the MP2 response density DB (a detailed breakdown of contributions to the total computation time is given in the Supporting Information). The latter formally scales as O (N 5 ) with system size and can be up to 20–30 times more expensive than the evaluation of D, as discussed in Section 2.3. However, the approximation used for the two-electron integrals, originating from the Fock operator, also has an effect on the timing. Figure 10 shows the total computation times for DSD-PBEP86 NMR CSC calculations on idealized linear alkane chains (CH3 (CH2 )n−1 H). The calculations using RIJK differ in how the LHSs of the CPSCF and z-vector equations are handled: in one case – denoted RIJK(RIJONX) – RI is only used for the Coulomb part (where applicable) and the exchange part is calculated in an integral-direct fashion, while in the other – denoted RIJK(RITrafo) – the (ia | jb) and (ij | ab) integrals are generated by an RI transformation (note that the AuxC basis is used for this) and stored on disk, as discussed in Section 2.3. The RIJK(RITrafo) option is fastest up to about C10 H22 , after which the RIJCOSX-L approach is faster. RIJK(RIJONX) is slightly faster than RIJCOSXL for the two smallest systems but for C15 H32 it already takes about 44 % more time. For the larger systems, the RIJK(RITrafo) computation times are between those of the other two options. However, this result depends on the speed of disk I/O operations (a RAID 0 33



hard disk configuration was used in this case). Therefore it can be concluded that for very small systems the choice of approximation is largely immaterial, while for larger systems of 100 electrons or more the RIJCOSX-L approximation is recommended.

34


Page 34 of 51

Page 35 of 51 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 4: Comparison of computation times (in minutes) for some medium-sized systems. Natoms , Nel , Nbas , NauxC , NauxJ , and Ngrid denote the number of atoms, electrons, contracted OBS, AuxC, and AuxJ functions, and COSX grid points for the system, respectively. Grid 4 was used for XC functional integration. The calculations were performed on 8 Intel Xeon E5-2687 v4 3.0 GHz cores with 8 GB RAM per core.

Natoms Nel Nbas (pS2) NauxC (cw3C) NauxJ (def2-JK) Ngrid (COSX-L)

Benzene

Caffeine

Coronene

Penicillin

Tweezer

12 42 300 846 558 23218

24 102 644 1854 1242 47803

36 156 1032 3024 2016 74108

41 176 1087 3158 2114 78287

92 374 2520 7296 4852 172573

1.4 2.3

1.7 1.6

10.7 13.5

9.2 4.3 3.0 14.0 108.7

7.6 3.7 4.3 18.4 144.4

63.1 32.2 48.0 668.0 7831.2

5.5 2.6 0.9 27.9 133.8

6.0 3.2 1.6 39.5 185.6

3.8 16.5 139.2 26.2 170.8

3.4 15.6 178.3 35.4 236.0

TPSS SCF (RI) RHS (RI)

0.1 0.1

0.6 0.7

DSD-PBEP86/RIJCOSX-L SCF (RIJCOSX) RHS (RIJCOSX) CPSCF (RIJCOSX) RI-MP2: D RI-MP2: DB

0.4 0.2 0.1 0.3 1.7

2.8 1.3 1.0 2.9 19.0

DSD-PBEP86/RIJK(RITrafo) SCF (RIJK) RHS (RIJK) CPSCF (RITrafo) RI-MP2: D RI-MP2: DB

0.2 0.1 0.0 0.2 1.5

1.4 0.7 0.1 3.5 21.1 Comparison

Total Total Total Total Total a

TPSS PBEP86/RIJCOSX-L DSD-PBEP86/RIJCOSX-L PBEP86/RIJK(RITrafo)a DSD-PBEP86/RIJK(RITrafo)

0.3 0.7 2.7 0.3 2.0

1.3 5.1 27.0 4.1 26.9

Including the RI transformation and storage of (ia | jb) and (ij | ab) integrals.

35


25.4 143.4 8642.7


Page 36 of 51

A further illustration is given in Table 4 for several “real world” systems: the benzene, caffeine, coronene, and penicillin G molecules and the “tweezer” host–guest complex discussed in refs 112, 80, 13, and 37. 113 The latter is the largest system studied here with 374 electrons and 2520 basis functions. Note that, to reduce memory requirements, the perturbed amplitudes for each batch were stored on disk, as discussed in Section 2.3. For comparison, TPSS/pS2 calculations were also performed on these systems. Table 4 provides separate timings for different calculation parts: the SCF solution, assembly of the CPSCF RHS, solution of the CPSCF equations, and calculation of the MP2 relaxed density and response density matrices. Note that for TPSS only the first two steps are necessary. The final lines of the table allow for a quick comparison of the total computation time required for NMR CSC calculations with a pure DF (TPSS), a hybrid DF (here taken as the SCF part of DSD-PBEP86), and a DHDF (DSD-PBEP86). Due to the efficiency of the RIJCOSX approximation, hybrid DFT calculations are consistently only a few times more expensive than pure DFT ones. However, the cost of DHDFT quickly grows to more than an order of magnitude above hybrid DFT, with the largest calculation taking 6 days to complete and requiring 1 TB of disk space. It is clear that in order to apply DHDFT to much larger systems, a local correlation approximation is needed, such as the ones employed by Gauss and Werner, 12 Loibl and Schütz,, 13 and Maurer and Ochsenfeld. 17 We are currently working on the development and implementation of NMR CSC calculations at the level of DLPNOMP2. 114

5

Conclusion

In this work NMR chemical shielding calculations at the DHDFT level using GIAOs and the RI-MP2 approximation were reported for the first time and were shown to be a valuable tool for accurate prediction of NMR chemical shifts. In a benchmark study of 34 1 H, 15

13

C,

N,17 O, 19 F, and 31 P shielding constants, the DSD-PBEP86 functional produced the smallest

36


Page 37 of 51 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


deviation from the CCSD(T) reference data – MAREδ = 1.9 % – which is significantly better than MP2 (MAREδ = 4.1 %) and common meta-GGA DFs like M06L (MAREδ = 5.4 %) and TPSS (MAREδ = 6.4 %). A more extensive evaluation of different DHDFs, as well as explicit optimization of such functionals, could further improve upon this result. It was shown that triple-zeta basis sets (pS2) are not sufficiently close to the CBS limit for perturbative approaches. However, fortuitous error compensation was observed in most cases with the notable exception of SCS-MP2. Errors due to the RI and COSX approximations were studied and found to be negligible if the recommended auxiliary basis sets (def2-JK and cw3C or cw4C) and grid settings (COSX-L) are used. The frozen core approximation, on the other hand, was found to introduce significant deviations, which do not cancel out in RCSs. The DSD-PBEP86/pS2/cw3C/RIJCOSX-L computational protocol was proposed as the best compromise between cost and accuracy, routinely applicable to systems of 100– 400 electrons. All features discussed in this work will be available in the upcoming release of the ORCA program. An implementation of NMR chemical shielding calculations at the DLPNO-MP2 level is currently in development.

Acknowledgement The authors are grateful to the Max Planck Society for generous support of our work. We thank Ute Becker for her assistance with the parallelization of the implemented code. G.L.S. thanks Peter Pinski for the helpful discussions.

Supporting Information Available Discussion of perturbed canonical orbitals and XC functional, frozen core, and implicit solvation terms, as well as additional results, referenced in the main text (PDF). Spreadsheet of all calculated chemical shielding values (XLSX). Coordinates of the large systems used in Section 4 (TXT). 37



Page 38 of 51

References (1) Gauss, J. Calculation of NMR chemical shifts at second-order many-body perturbation theory using gauge-including atomic orbitals. Chem. Phys. Lett. 1992, 191, 614–620. (2) Gauss, J. Effects of electron correlation in the calculation of nuclear magnetic resonance chemical shifts. J. Chem. Phys. 1993, 99, 3629–3643. (3) Gauss, J. GIAO-MBPT(3) and GIAO-SDQ-MBPT(4) calculations of nuclear magnetic shielding constants. Chem. Phys. Lett. 1994, 229, 198–203. (4) Gauss, J.; Stanton, J. F. Coupled-cluster calculations of nuclear magnetic resonance chemical shifts. J. Chem. Phys. 1995, 103, 3561–3577. (5) Gauss, J.; Stanton, J. F. Gauge-invariant calculation of nuclear magnetic shielding constants at the coupled–cluster singles and doubles level. J. Chem. Phys. 1995, 102, 251–253. (6) Gauss, J.; Stanton, J. F. Perturbative treatment of triple excitations in coupled-cluster calculations of nuclear magnetic shielding constants. J. Chem. Phys. 1996, 104, 2574– 2583. (7) Auer, A. A.; Gauss, J.; Stanton, J. F. Quantitative prediction of gas-phase 13 C nuclear magnetic shielding constants. J. Chem. Phys. 2003, 118, 10407–10417. (8) Harding, M. E.; Lenhart, M.; Auer, A. A.; Gauss, J. Quantitative prediction of gasphase

19

F nuclear magnetic shielding constants. J. Chem. Phys. 2008, 128, 244111.

(9) Auer, A. A. Quantitative prediction of gas-phase

17

O nuclear magnetic shielding con-

stants. J. Chem. Phys. 2009, 131, 024116. (10) Prochnow, E.; Auer, A. A. Quantitative prediction of gas-phase

15

magnetic shielding constants. J. Chem. Phys. 2010, 132, 064109. 38


N and

31

P nuclear

Page 39 of 51 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


(11) Flaig, D.; Maurer, M.; Hanni, M.; Braunger, K.; Kick, L.; Thubauville, M.; Ochsenfeld, C. Benchmarking Hydrogen and Carbon NMR Chemical Shifts at HF, DFT, and MP2 Levels. J. Chem. Theory Comput. 2014, 10, 572–578. (12) Gauss, J.; Werner, H.-J. NMR chemical shift calculations within local correlation methods: the GIAO-LMP2 approach. Phys. Chem. Chem. Phys. 2000, 2, 2083–2090. (13) Loibl, S.; Schütz, M. NMR shielding tensors for density fitted local second-order Møller-Plesset perturbation theory using gauge including atomic orbitals. J. Chem. Phys. 2012, 137, 084107. (14) Ochsenfeld, C.; Kussmann, J.; Koziol, F. Ab Initio NMR Spectra for Molecular Systems with a Thousand and More Atoms: A Linear-Scaling Method. Angew. Chemie Int. Ed. 2004, 43, 4485–4489. (15) Kussmann, J.; Ochsenfeld, C. Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within Hartree-Fock and density-functional theory. J. Chem. Phys. 2007, 127, 054103. (16) Beer, M.; Kussmann, J.; Ochsenfeld, C. Nuclei-selected NMR shielding calculations: A sublinear-scaling quantum-chemical method. J. Chem. Phys. 2011, 134, 074102. (17) Maurer, M.; Ochsenfeld, C. A linear- and sublinear-scaling method for calculating NMR shieldings in atomic orbital-based second-order Møller-Plesset perturbation theory. J. Chem. Phys. 2013, 138, 174104. (18) Grimme, S. Improved second-order Møller–Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies. J. Chem. Phys. 2003, 118, 9095–9102. (19) Fink, R. F. Spin-component-scaled Møller–Plesset (SCS-MP) perturbation theory: A

39



generalization of the MP approach with improved properties. J. Chem. Phys. 2010, 133, 174113. (20) Grimme, S.; Goerigk, L.; Fink, R. F. Spin-component-scaled electron correlation methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 886–906. (21) Maurer, M.; Ochsenfeld, C. Spin Component-Scaled Second-Order Møller–Plesset Perturbation Theory for Calculating NMR Shieldings. J. Chem. Theory Comput. 2015, 11, 37–44. (22) Lee, A. M.; Handy, N. C.; Colwell, S. M. The density functional calculation of nuclear shielding constants using London atomic orbitals. J. Chem. Phys. 1995, 103, 10095– 10109. (23) Malkin, V.; Malkina, O.; Salahub, D. Calculations of NMR shielding constants beyond uncoupled density functional theory. IGLO approach. Chem. Phys. Lett. 1993, 204, 87–95. (24) Malkin, V. G.; Malkina, O. L.; Casida, M. E.; Salahub, D. R. Nuclear Magnetic Resonance Shielding Tensors Calculated with a Sum-over-States Density Functional Perturbation Theory. J. Am. Chem. Soc. 1994, 116, 5898–5908. (25) Malkin, V. G.; Malkina, O. L.; Eriksson, L. A.; Salahub, D. R. Theor. Comput. Chem.; 1995; Vol. 2; pp 273–347. (26) Reimann, S.; Ekström, U.; Stopkowicz, S.; Teale, A. M.; Borgoo, A.; Helgaker, T. The importance of current contributions to shielding constants in density-functional theory. Phys. Chem. Chem. Phys. 2015, 17, 18834–18842. (27) Grayce, C. J.; Harris, R. A. Magnetic-field density-functional theory. Phys. Rev. A 1994, 50, 3089–3095.

40


Page 40 of 51

Page 41 of 51 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


(28) Vignale, G.; Rasolt, M. Density-functional theory in strong magnetic fields. Phys. Rev. Lett. 1987, 59, 2360–2363. (29) Vignale, G.; Rasolt, M. Current- and spin-density-functional theory for inhomogeneous electronic systems in strong magnetic fields. Phys. Rev. B 1988, 37, 10685–10696. (30) Tellgren, E. I.; Kvaal, S.; Sagvolden, E.; Ekström, U.; Teale, A. M.; Helgaker, T. Choice of basic variables in current-density-functional theory. Phys. Rev. A 2012, 86, 062506. (31) Becke, A. D. Current density in exchange-correlation functionals: Application to atomic states. J. Chem. Phys. 2002, 117, 6935–6938. (32) Maximoff, S. N.; Scuseria, G. E. Nuclear magnetic resonance shielding tensors calculated with kinetic energy density-dependent exchange-correlation functionals. Chem. Phys. Lett. 2004, 390, 408–412. (33) Tao, J. Explicit inclusion of paramagnetic current density in the exchange-correlation functionals of current-density functional theory. Phys. Rev. B 2005, 71, 205107. (34) Furness, J. W.; Verbeke, J.; Tellgren, E. I.; Stopkowicz, S.; Ekström, U.; Helgaker, T.; Teale, A. M. Current Density Functional Theory Using Meta-Generalized Gradient Exchange-Correlation Functionals. J. Chem. Theory Comput. 2015, 11, 4169–4181. (35) Tellgren, E. I.; Teale, A. M.; Furness, J. W.; Lange, K. K.; Ekström, U.; Helgaker, T. Non-perturbative calculation of molecular magnetic properties within current-density functional theory. J. Chem. Phys. 2014, 140, 034101. (36) Teale, A. M.; Lutnæs, O. B.; Helgaker, T.; Tozer, D. J.; Gauss, J. Benchmarking density-functional theory calculations of NMR shielding constants and spin–rotation constants using accurate coupled-cluster calculations. J. Chem. Phys. 2013, 138, 024111. 41



(37) Stoychev, G. L.; Auer, A. A.; Izsák, R.; Neese, F. Self-Consistent Field Calculation of Nuclear Magnetic Resonance Chemical Shielding Constants Using Gauge-Including Atomic Orbitals and Approximate Two-Electron Integrals. J. Chem. Theory Comput. 2018, 14, 619–637. (38) Keal, T. W.; Tozer, D. J. The exchange-correlation potential in Kohn–Sham nuclear magnetic resonance shielding calculations. J. Chem. Phys. 2003, 119, 3015–3024. (39) Keal, T. W.; Tozer, D. J. A semiempirical generalized gradient approximation exchange-correlation functional. J. Chem. Phys. 2004, 121, 5654–5660. (40) Zhao, Y.; Truhlar, D. G. Improved Description of Nuclear Magnetic Resonance Chemical Shielding Constants Using the M06-L Meta-Generalized-Gradient-Approximation Density Functional. J. Phys. Chem. A 2008, 112, 6794–6799. (41) Perdew, J. P.; Schmidt, K. Jacob’s ladder of density functional approximations for the exchange-correlation energy. AIP Conf. Proc. 2001; pp 1–20. (42) Grimme, S. Semiempirical hybrid density functional with perturbative second-order correlation. J. Chem. Phys. 2006, 124, 034108. (43) Goerigk, L.; Grimme, S. Double-hybrid density functionals. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 576–600. (44) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. (45) 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. (46) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. 42


Page 42 of 51

Page 43 of 51 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


(47) Kozuch, S.; Gruzman, D.; Martin, J. M. L. DSD-BLYP: A General Purpose Double Hybrid Density Functional Including Spin Component Scaling and Dispersion Correction. J. Phys. Chem. C 2010, 114, 20801–20808. (48) Kozuch, S.; Martin, J. M. L. DSD-PBEP86: in search of the best double-hybrid DFT with spin-component scaled MP2 and dispersion corrections. Phys. Chem. Chem. Phys. 2011, 13, 20104. (49) Kozuch, S.; Martin, J. M. L. Spin-component-scaled double hybrids: An extensive search for the best fifth-rung functionals blending DFT and perturbation theory. J. Comput. Chem. 2013, 34, 2327–2344. (50) Goerigk, L.; Grimme, S. Efficient and Accurate Double-Hybrid-Meta-GGA Density Functionals—Evaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 291–309. (51) Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Grimme, S. A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochemistry, kinetics and noncovalent interactions. Phys. Chem. Chem. Phys. 2017, 19, 32184–32215. (52) Neese, F.; Schwabe, T.; Grimme, S. Analytic derivatives for perturbatively corrected “double hybrid” density functionals: Theory, implementation, and applications. J. Chem. Phys. 2007, 126, 124115. (53) Kossmann, S.; Neese, F. Efficient Structure Optimization with Second-Order ManyBody Perturbation Theory: The RIJCOSX-MP2 Method. J. Chem. Theory Comput. 2010, 6, 2325–2338. (54) Kossmann, S. Efficient Novel Approaches for the Calculation of Molecular Response Properties : Second-Order Many-Body Perturbation and Double-Hybrid Density 43



Functional Theory. Ph.D. Dissertation, Rheinischen Friedrich-Wilhelms-Universität Bonn, 2011. (55) Kutzelnigg, W. Ab initio calculation of molecular properties. J. Mol. Struct. THEOCHEM 1989, 202, 11–61. (56) London, F. Théorie quantique des courants interatomiques dans les combinaisons aromatiques. J. Phys. Radium 1937, 8, 397–409. (57) Hameka, H. On the nuclear magnetic shielding in the hydrogen molecule. Mol. Phys. 1958, 1, 203–215. (58) Ditchfield, R. Molecular Orbital Theory of Magnetic Shielding and Magnetic Susceptibility. J. Chem. Phys. 1972, 56, 5688–5691. (59) Ditchfield, R. Self-consistent perturbation theory of diamagnetism. Mol. Phys. 1974, 27, 789–807. (60) 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. (61) 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. (62) Neese, F. The ORCA program system. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 73–78. (63) Neese, F. Software update: the ORCA program system, version 4.0. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, e1327. (64) Salter, E.; Trucks, G. W.; Fitzgerald, G.; Bartlett, R. J. Theory and application of MBPT(3) gradients: The density approach. Chem. Phys. Lett. 1987, 141, 61–70. 44


Page 44 of 51

Page 45 of 51 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


(65) Trucks, G. W.; Salter, E.; Sosa, C.; Bartlett, R. J. Theory and implementation of the MBPT density matrix. An application to one-electron properties. Chem. Phys. Lett. 1988, 147, 359–366. (66) Salter, E. A.; Trucks, G. W.; Bartlett, R. J. Analytic energy derivatives in many-body methods. I. First derivatives. J. Chem. Phys. 1989, 90, 1752–1766. (67) Jørgensen, P.; Helgaker, T. Møller-Plesset energy derivatives. J. Chem. Phys. 1988, 89, 1560–1570. (68) Gauss, J.; Stanton, J. F.; Bartlett, R. J. Coupled-cluster open-shell analytic gradients: Implementation of the direct product decomposition approach in energy gradient calculations. J. Chem. Phys. 1991, 95, 2623–2638. (69) Hylleraas, E. A. Über den Grundterm der Zweielektronenprobleme von H-, He, Li+, Be++ usw. Zeitschrift für Phys. 1930, 65, 209–225. (70) Pulay, P.; Saebo, S.; Meyer, W. An efficient reformulation of the closed-shell selfconsistent electron pair theory. J. Chem. Phys. 1984, 81, 1901–1905. (71) Pulay, P.; Saebø, S. Orbital-invariant formulation and second-order gradient evaluation in Møller-Plesset perturbation theory. Theor. Chim. Acta 1986, 69, 357–368. (72) Vahtras, O.; Almlöf, J.; Feyereisen, M. W. Integral approximations for LCAO-SCF calculations. Chem. Phys. Lett. 1993, 213, 514–518. (73) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of approximate integrals in ab initio theory. An application in MP2 energy calculations. Chem. Phys. Lett. 1993, 208, 359–363. (74) 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.

45



(75) Weigend, F.; Häser, M. RI-MP2: first derivatives and global consistency. Theor. Chem. Accounts Theory, Comput. Model. (Theoretica Chim. Acta) 1997, 97, 331–340. (76) Handy, N. C.; Schaefer, H. F. On the evaluation of analytic energy derivatives for correlated wave functions. J. Chem. Phys. 1984, 81, 5031–5033. (77) Weigend, F. A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. Phys. Chem. Chem. Phys. 2002, 4, 4285–4291. (78) Neese, F.; Wennmohs, F.; Hansen, A.; Becker, U. Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A ’chain-of-spheres’ algorithm for the Hartree-Fock exchange. Chem. Phys. 2009, 356, 98–109. (79) Kossmann, S.; Neese, F. Comparison of two efficient approximate Hartee–Fock approaches. Chem. Phys. Lett. 2009, 481, 240–243. (80) Loibl, S.; Manby, F. R.; Schutz, M. Density fitted, local Hartree–Fock treatment of NMR chemical shifts using London atomic orbitals. Mol. Phys. 2010, 108, 477–485. (81) 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. (82) Lee, T. J.; Rendell, A. P. Analytic gradients for coupled-cluster energies that include noniterative connected triple excitations: Application to cis- and trans-HONO. J. Chem. Phys. 1991, 94, 6229–6236. (83) Helgaker, T.; Jørgensen, P. An electronic Hamiltonian for origin independent calculations of magnetic properties. J. Chem. Phys. 1991, 95, 2595. (84) Gauss, J. In Mod. Methods Algorithms Quantum Chem. Proceedings, Winterschool, 2r1-25 Febr. 2000 Forschungszentrum Jülich, Ger.; Grotendorst, J., Ed.; John von Neumann Institute for Computing: Jülich, Germany, 2000; pp 541–592. 46


Page 46 of 51

Page 47 of 51 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


(85) Lee, T. J.; Racine, S. C.; Rice, J. E.; Rendell, A. P. On the orbital contribution to analytical derivatives of perturbation theory energies. Mol. Phys. 1995, 85, 561–571. (86) Jameson, C. J. Encycl. Magn. Reson.; John Wiley & Sons, Ltd: Chichester, UK, 2007. (87) Raynes, W. T. In Nucl. Magn. Reson. Vol. 7 ; Abraham, R. J., Ed.; Royal Society of Chemistry: Cambridge, 1978; pp 1–25. (88) Hindermann, D. K.; Cornwell, C. D. Vibrational Corrections to the Nuclear-Magnetic Shielding and Spin–Rotation Constants for Hydrogen Fluoride. Shielding Scale for 19 F. J. Chem. Phys. 1968, 48, 4148–4154. (89) Hindman, J. C. Proton Resonance Shift of Water in the Gas and Liquid States. J. Chem. Phys. 1966, 44, 4582–4592. (90) Jameson, A.; Jameson, C. J. Gas-phase

13

C chemical shifts in the zero-pressure limit:

refinements to the absolute shielding scale for

13

C. Chem. Phys. Lett. 1987, 134,

461–466. (91) Sundholm, D.; Gauss, J.; Schäfer, A. Rovibrationally averaged nuclear magnetic shielding tensors calculated at the coupled-cluster level. J. Chem. Phys. 1996, 105, 11051–11059. (92) Raymonda, J.; Klemperer, W. Molecular Beam Electric Resonance Spectrum of 14

31

P

N. J. Chem. Phys. 1971, 55, 232–233.

(93) Jameson, C. J.; Jameson, A. K.; Oppusunggu, D.; Wille, S.; Burrell, P. M.; Mason, J. 15

N nuclear magnetic shielding scale from gas phase studies. J. Chem. Phys. 1981,

74, 81–88. (94) Makulski, W.; Jackowski, K. The 17 O nuclear magnetic shielding scale from gas-phase measurements. J. Mol. Struct. 2003, 651-653, 265–269.

47



(95) Jameson, C. J.; Jameson, A. K.; Burrell, P. M.

19

Page 48 of 51

F nuclear magnetic shielding scale

from gas phase studies. J. Chem. Phys. 1980, 73, 6013–6020. (96) Nebgen, J.; Rose, W.; Metz, F. The

19

F nuclear magnetic resonance spectra of liquid

and gaseous fluorine, oxygen difluoride, and nitrogen trifluoride. J. Mol. Spectrosc. 1966, 20, 72–74. (97) Jameson, C. J.; Jameson, A. K.; Honarbakhsh, J. 19 F nuclear magnetic shielding scale from gas phase studies. II. J. Chem. Phys. 1984, 81, 5266–5267. (98) Jameson, C. J.; De Dios, A.; Keith Jameson, A. Absolute shielding scale for

31

P from

gas-phase NMR studies. Chem. Phys. Lett. 1990, 167, 575–582. (99) Jensen, F. Segmented Contracted Basis Sets Optimized for Nuclear Magnetic Shielding. J. Chem. Theory Comput. 2015, 11, 132–138. (100) Jensen, F. Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods. J. Chem. Theory Comput. 2008, 4, 719–727. (101) Reid, D. M.; Kobayashi, R.; Collins, M. A. Systematic Study of Locally Dense Basis Sets for NMR Shielding Constants. J. Chem. Theory Comput. 2014, 10, 146–152. (102) Stanton, J.; Gauss, J.; Harding, M.; Szalay, P. CFOUR, a quantum chemical program package. http://www.cfour.de. (103) Hättig, C. Optimization of auxiliary basis sets for RI-MP2 and RI-CC2 calculations: core-valence and quintuple-zeta basis sets for H to Ar and QZVPP basis sets for Li to Kr. Phys. Chem. Chem. Phys. 2005, 7, 59–66. (104) Peterson, K. A.; 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.

48


Page 49 of 51 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


(105) Stoychev, G. L.; Auer, A. A.; Neese, F. Automatic Generation of Auxiliary Basis Sets. J. Chem. Theory Comput. 2017, 13, 554–562. (106) Tarnopolsky, A.; Karton, A.; Sertchook, R.; Vuzman, D.; Martin, J. M. L. DoubleHybrid Functionals for Thermochemical Kinetics. J. Phys. Chem. A 2008, 112, 3–8. (107) Karton, A.; Tarnopolsky, A.; Lamère, J.-F.; Schatz, G. C.; Martin, J. M. L. Highly Accurate First-Principles Benchmark Data Sets for the Parametrization and Validation of Density Functional and Other Approximate Methods. Derivation of a Robust, Generally Applicable, Double-Hybrid Functional for Thermochemistry and Thermochemical Kinetics. J. Phys. Chem. A 2008, 112, 12868–12886. (108) We did also test the B2GP-PLYP parametrization and it performs somewhat better than B2PLYP with MAREδ = 3.4 %. The data are available in the Supporting Information. (109) Gauss, J.; Stanton, J. F. Adv. Chem. Phys.; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2002; Vol. 123; pp 355–422. (110) Weigend, F. Hartree–Fock exchange fitting basis sets for H to Rn. J. Comput. Chem. 2008, 29, 167–175. (111) Grids 1–7 employ 50-, 110-, 194-, 302-, 434-, 590-, and 770-point Lebedev angular grids and radial integration parameters of 4.34, 4.34, 4.34, 4.67, 5.01, 5.34, and 5.67, respectively. Default pruning settings in ORCA were used. (112) Brown, S. P.; Schaller, T.; Seelbach, U. P.; Koziol, F.; Ochsenfeld, C.; Klärner, F.G.; Spiess, H. W. Structure and Dynamics of the Host-Guest Complex of a Molecular Tweezer: Coupling Synthesis, Solid-State NMR, and Quantum-Chemical Calculations. Angew. Chemie Int. Ed. 2001, 40, 717–720. (113) Cartesian coordinates for these systems are provided in the Supporting Information. 49



(114) Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. J. Chem. Phys. 2015, 143, 034108.

50


Page 50 of 51

Page 51 of 51 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


Graphical TOC Entry

51


Efficient and Accurate Prediction of Nuclear ... - ACS Publications

Recommend Documents