for the calculation of NMR shielding tensors using density-fitting

Nuclei-selected atomic-orbital response-theory formulation

for the calculation of NMR shielding tensors using density-fitting

Chandan Kumar,^1,a)Thomas Kjærgaard,^2,b)Trygve Helgaker,¹and Heike Fliegl^1,c)

1Centre for Theoretical and Computational Chemistry (CTCC), Department of Chemistry, University of Oslo, P.O. Box 1033, Blindern, 0315 Oslo, Norway

2Department of Chemistry, Aarhus University, Langelandsgade 140, 8000 Aarhus C, Denmark

(Received 16 September 2016; accepted 1 December 2016; published online 20 December 2016) An atomic orbital density matrix based response formulation of the nuclei-selected approach of Beer, Kussmann, and Ochsenfeld [J. Chem. Phys.134, 074102 (2011)] to calculate nuclear magnetic resonance (NMR) shielding tensors has been developed and implemented into LSDalton allowing for a simultaneous solution of the response equations, which significantly improves the performance. The response formulation to calculate nuclei-selected NMR shielding tensors can be used together with the density-fitting approximation that allows efficient calculation of Coulomb integrals. It is shown that using density-fitting does not lead to a significant loss in accuracy for both the nuclei-selected and the conventional ways to calculate NMR shielding constants and should thus be used for applications with LSDalton.Published by AIP Publishing.[http://dx.doi.org/10.1063/1.4972212]

I. INTRODUCTION

The main drawback of standard quantum-chemical meth- ods such as Hartree–Fock^1–3 and Kohn–Sham⁴ theories is their formal cubic (N³) scaling behavior with molecular size N or with respect to the number of atomic basis functions.

Thus, calculations for large molecular systems of more than a thousand atoms are out of reach using conventional imple- mentations. To overcome this bottleneck, the development of linear-scaling methods has become an attractive alterna- tive.^5–8 Available linear-scaling methods efficiently compute not only the molecular energy but also response properties such as nuclear magnetic resonance (NMR) chemical shifts.^7–15 Illustrative reviews can be found in Refs.5and6.

A nucleus in a molecule is affected differently by an exter- nal magnetic field depending on its local electronic environ- ment. This effect is studied in NMR spectroscopy, where infor- mation about the local environment of the nuclei is extracted from the measured shielding or deshielding of resonance fre- quencies.^16–18Theoretically, the chemical shift of a nucleus is calculated as the difference of the NMR shielding constant of the nucleus of interest relative to that of a reference compound.

This value can directly be compared with experiment, neglect- ing temperature and solvent effects. The shielding constant is calculated as one third of the trace of the NMR shield- ing tensor, which is the key response quantity that must be computed.

The NMR shielding tensor is the second derivative of the molecular electronic energy with respect to the external mag- netic fieldBand the nuclear magnetic momentm_jof nucleus j. Since the NMR shielding tensor is bilinear in the pertur- bations, one may eliminate from the expressions (using the

a)chandan.kumar@kjemi.uio.no b)tkjaergaard@chem.au.dk c)heike.fliegl@kjemi.uio.no

2n + 1 rule) either all terms that involve the response of the wave function to the external magnetic fieldBor all terms that involve the response to the nuclear magnetic momentsm_j.^19,20 Traditionally, the second route has been followed since there are only three components of the magnetic fieldBbut poten- tially a large number of nuclear dipole momentsm_j^19–21that need to be taken into account. More detailed information on the derivation as well as on historical aspects of the NMR shield- ing tensor onab initio and density functional theory (DFT) level following the traditional route can be found, for exam- ple, in Refs.22-30and references therein. Application aspects have been discussed in a recent review by Lodewyk and co- workers.³¹Locality is not exploited in the traditional approach.

Thus, the calculation of the shielding tensor of one particular nucleus in a large system implies the calculation of unneces- sary information that cannot be avoided. In practice, this means that, albeit often only one particular nucleus is of interest, the shielding tensors for all nuclei are calculated. In cases where the NMR shielding constants of only a few among several hun- dred nuclei are required, it is beneficial to use a method that is able to target the nuclei of interest exclusively and avoid unnecessary calculations.¹¹

Such a method was recently presented by Beer, Kuss- mann, and Ochsenfeld, who calculate directly the response of the molecular system to selected nuclear magnetic moments m_j in sublinear time, making it possible to calculate shield- ing constants of very large compounds.¹⁴ Their approach has been implemented at the Hartree–Fock, Kohn–Sham, and second-order Møller–Plesset levels of theory.^14,15Beer and co- workers demonstrated that their alternative approach scales better and therefore becomes faster for very large systems, with more than a thousand atoms. Moreover, sublinear scal- ing is achieved by a sophisticated truncation of the London or gauge-including atomic-orbital (GIAO) expansion.^32–35 The key quantity computed in their formulation is the first-order perturbed density matrix with respect to the magnetic moment

0021-9606/2016/145(23)/234108/9/$30.00 145, 234108-1 Published by AIP Publishing.

of the selected nucleus, for which the coupled-perturbed self-consistent-field (CPSCF) equations are solved.^13,14 In pure Kohn–Sham theory, the CPSCF equations become uncou- pled, improving efficiency of the calculations. The uncoupled Kohn–Sham scheme for calculating NMR shielding tensors introduced by Malkin et al.,³⁶ extended to London orbitals by Leeet al.,³⁷ and implemented in the Dalton program by Helgakeret al.³⁸is also part of LSDalton and is used in this work. Inspired by the approach of Beer and co-workers,¹⁴we solve for nuclear-moment response rather than field responses but use instead the atomic-orbital (AO) density matrix formu- lation of Larsenet al.^9,10as implemented in LSDalton. In this work we present a new way to derive the final equations used in Ref.14in an atomic orbital basis including density-fitting and local density-fitting, which has not been tried before. The per- turbed density matrices are obtained using the general response solver developed by Corianiet al.⁸ The response framework and the general response solver allow for all requested response equations to be solved simultaneously, significantly improving performance. In contrast to our work, Beer and co-workers¹⁴ extended the CPSCF equations to suit their needs. Instead using our general response framework, the response solver can be used directly^39,40 as it is done in the present work.

The GIAO basis-set truncation described in Ref.14to achieve sublinear scaling is not included in the present implementa- tion. Instead, the density-fitting or resolution-of-the-identity approximation^41–48for efficient computation of Coulomb inte- grals has been implemented and tested for both the standard and the alternative approach. The present implementation allows for the calculation of nuclei-selected shielding tensors on DFT and Hartree–Fock level.

The paper is organized as follows. In SectionII, the theory of nuclei-selected NMR approach implemented in LSDalton is outlined, starting from the Kohn–Sham energy expression.

Computational details are given in SectionIIIand some appli- cations are presented and discussed in SectionIV. SectionV contains some concluding remarks.

II. THEORY

The NMR shielding tensor elements can be written as second derivatives of the Kohn–Sham energy with respect to a Cartesian component of the external magnetic fieldB_iand a Cartesian component of the magnetic momentm_jbelonging to a nucleus in the molecule,

ij= d²E

dBidm_{j Bi=0,mj=0}. (1) In the following, we discuss the Kohn–Sham energy and den- sity matrix before considering the evaluation of the NMR shielding tensor.

A. The Kohn–Sham energy and density matrix

We consider a molecule ofNelectrons and with nuclei of chargeZ_Aat positionsRAand wish to express the Kohn–Sham energy in terms of a density matrix expanded in a basis of AOs

µ(x). Introducing an anti-Hermitian matrixXthat contains our variational parameters, a valid Kohn–Sham AO density matrix can be expressed as⁴⁹

D(X)=exp ( XS)D0exp (SX), (2) whereSis the AO overlap matrix with elements

S_µ⌫ =

⌅ ⇤

µ(x) _⌫(x) dx (3)

and where the reference density matrixD0satisfies the usual symmetry, rank, and idempotency conditions of a valid density matrix.¹⁰In terms of such an AO density matrix, the Kohn–

Sham energy is given by E(X)= trD(X)h+1

2trD(X)G(D(X))+E_xc(⇢)+h_nuc, (4) where the one-electron Hamiltonian integrals are given by

h_µ =

⌅

µ(x)* ,

1

2r² X

A

Z_A

|R_A r|+

- (x) dx, (5) the two-electron Coulomb interactions and exchange interac- tions (weighted by ) are described by the matrix elements

G(D)_µ =X

⇢

(g_{µ ⇢} g_{µ ⇢} )D _⇢, (6)

g_{µ ⇢} =

⌅⌅ µ(x1) (x1) _⇢(x2) (x2)

|x₁ x₂| dx1dx2, (7)

andE_xc(⇢) is the exchange-correlation functional expressed in terms of the electron density

⇢(x)=X

µ⌫

µ(x)D_µ⌫(X) _⌫(x). (8) The final term in the energy expression hnuc is an additive nuclear–nuclear repulsion term.

Let us consider the conditions that must be imposed on the reference matrixD0to be an optimal density matrix. Per- forming an asymmetric Baker–Campbell–Hausdorff (BCH) expansion,⁴⁹we obtain

D(X)=D0+[D0,X]_S+1

2[[D0,X]_S,X]_S+· · ·, (9)

where we have introduced the overlap-dependent S commuta- tor [D,X]_S=DSX XSD. Setting the first-order variation in the Kohn–Sham energy to zero, we obtain the conditions

E=trF D(X)= trF[D, X]_S=tr (FD0S SD0F) X=0, (10) which must hold for each X. The variational conditions can therefore be written as

FD₀S=SD₀F (11) which uniquely determines the optimal Kohn–Sham density matrixD0.

B. The Kohn–Sham NMR shielding tensor

We evaluate Eq.(1)by differentiating the energy expres- sion in Eq.(4)and noting thath_nucdoes not depend onB_iand m_j, we obtain

ij =trD0h^Bi,mj+¹₂ trD0G^Bi,mj(D0)+trD^Bi,mj[h+G(D0)]

+E^B_xc^i,mj+trD^Bi⇥h^mj+G^mj(D₀)⇤ +trD^mj f

h^Bi +G^Bi(D₀)g

+trD^mjG(D^Bi), (12)

where superscripts B_i and m_j denote differentiation with respect to B_i andm_j, respectively, atB_i =m_j = 0, while D0

is the density matrix atB_i=m_j= 0. In deriving the expression, we have used the symmetry trD^BG(D^mj)=trD^mjG(D^B). The derivative of the exchange-correlation energy is given by

E_xc^Bi,mj=trF_xcD^Bi,mj+trG_xc(D^Bi)D^mj+trF^B_xcⁱD^mj (13) in terms of the matrices

(Fxc)_µ⌫=

⌅ E_xc(⇢)

⇢(r) ^µ(r) ⌫(r) dr, (14)

Gxc(M)_µ⌫=X

⇢

M_⇢

⌅⌅ 2E_xc(⇢)

⇢(r) ⇢(r⁰) ^µ(r) _⌫(r)

⇥ ^⇢(r⁰) (r⁰) drdr⁰, (15) andF^B_xcⁱ given in Eq. (95) of Ref.50. We note that the first and second functional derivatives of Exc are the exchange- correlation potential and kernel, respectively. By the anti- symmetry of D^Bi, we have Gxc(D^Bi) = 0. Introducing the Kohn–Sham matrix

F=h+G(D0)+Fxc (16) and noting that G^mj(D₀)=G^Bi,mj(D₀)=0 but G^Bi(D₀),0 since the two-electron integrals are independent of the nuclear magnetic moments but depend on the magnetic field, we arrive at the following expression for the shielding tensor:

ij=trD₀h^Bi,mj+tr (D^Bih^mj +D^mjh^Bi)+trD^Bi,mjF

+trD^mjF^B_xcⁱ +trD^mjG^Bi(D0)+trD^mjG(D^Bi). (17) In pure Kohn–Sham theory, we have = 0 in Eq.(6) and G(D^Bi) then vanishes.^36–38

C. Derivatives of the density matrix

In this section, we discuss how the perturbed den- sity matrices needed to compute the NMR shield- ing tensor are constructed, following the procedure of Refs.9and10.

Differentiating the BCH expansionD(X) in Eq.(9)with respect to a general perturbation parameteraorbatX=0, we obtain the first and second orders¹⁰

D^a =D^a₀+[D₀,X]^a_S, (18) D^ab=D^ab₀ +[D₀,X]^ab_S +1

2[[D₀,X]_S,X]^ab_S

=P_ab⇣

D^ab₀ +[D₀,X^ab]_S+2[D^a₀,X^b]_S +2[D₀,X^a]_S^b+[[D₀,X^a]S,X^b]_S⌘

, (19)

where we have introduced the permutation operator P_abM^aN^b=¹₂(M^aN^b + M^bN^a). The derivatives of D₀ are obtained from the idempotency relationD0SD0=D0as¹⁰

D^a₀= D0S^aD0, (20) D^ab₀ = D0S^abD0+2PHP_abD0S^aD0S^bD0, (21)

whereP_H which extracts the Hermitian part of a matrix A asP_HA=¹₂(A+A^†). The first-order perturbed density matrix with respect toB_iandm_jthen becomes

D^Bi =D^B₀ⁱ +[D0,X^Bi]_S= D0S^BiD0+[D0,X^Bi]_S, (22) D^mj =D^m₀^j +[D0,X^mj]_S=[D0,X^mj]_S, (23) where we used the fact that the overlap matrix S does not depend onm_j.

Writing the NMR shielding tensor in terms of response vectors we obtain

ij =trD0h^Bi,mj trD0S^BiD0h^mj +tr⇣

[D₀,X^Bi]Sh^mj+[D₀,X^mj]Sh^Bi⌘ +trD^Bi,mjF+tr [D0,X^mj]_S⇣

G^Bi(D0)+F^B_xcⁱ⌘ trD₀S^BiD₀G([D₀,X^mj]_S)

+tr [D₀,X^Bi]SG([D₀,X^mj]S). (24) This expression can be further simplified by considering that ij is bilinear in the perturbations, implying that we may eliminate either the terms that involve the response of the density matrix to the external magnetic field X^Bi or the terms containing the response of the density matrix to the nuclear magnetic moment X^mj.^10,20 In the conven- tional formulation, the terms containingX^mj are eliminated, yielding²⁰

ij =trD₀h^Bi,mj trD₀S^BiD₀h^mj+tr [D₀,X^Bi]_Sh^mj. (25) Alternatively, eliminating the terms containingX^Bi, we find

ij=trD0h^Bi,mj trD0S^BiD0h^mj +tr ([D^B₀ⁱ,X^mj]_S+[D₀,X^mj]_SBi)F +tr [D₀,X^mj]S

⇣h^Bi+G^Bi(D₀)+F^B_xcⁱ⌘

trD0S^BiD0G([D0,X^mj]_S). (26) From Eqs. (25) and (26), we see that the expression for the NMR shielding tensor incorporates two different types of response matrices: X^Bi and X^mj, respectively. Both can be obtained using the response equation solver described in Refs. 8 and 12. In Sec. II D, we discuss how the density matrix derivatives are obtained from the response equations.

D. Response equations

To obtainD^Bi andD^mj, we must evaluateX^Bi andX^mj, respectively. Differentiating the variational condition given in Eq.(11), we obtain

P_A(F^aD₀S+FD^aS+FD₀S^a)=0, (27) wherePAA=¹₂(A A^†). Introducing the first derivative of the Kohn–Sham matrix

F^a=h^a+G^a(D0)+G(D^a)+F_xc^a(D0) (28) and the notation

K(A)=PA[G(A)D0S+FAS], (29) we obtain

K(D^a)= P_A⇥h^aD₀S+G^a(D₀)D₀S+F^a_xcDS+FD₀S^a⇤. (30)

Introducing the expression for the perturbed densityD^a=D^a₀ +[D0,X^a]_S, we can rearrange Eq.(30)in the manner

K([D₀,X^a]_S)= P_A[h^aD₀S+G^a(D₀)D₀S+F^a_xcD₀S +FD0S^a+K(D^a₀)]. (31) Introducing the notation

E^[2]X=K([D0,X]_S) (32) and noting thatG^a(D0),F^a_xc,S^a, andD^a₀vanish whena=mj

but not whena=B_i, we arrive at the linear response equations E^[2]X^mj = P_A h^mjD0S , (33) E^[2]X^Bi = P_A⇣

h^BiD0S+G^Bi(D0)D0S+F^B_xcⁱD0S +FD₀S^Bi+K(D^B₀ⁱ)⌘

. (34)

These equations are solved iteratively, applyingE^[2]to a trail vectorbat each iteration,^8,51

E^[2]b=G([D0,b]_S)D0S SD0G([D0,b]_S)

+F[D0,b]_SS S[D0,b]_SF. (35) The response equations Eqs.(33)and(34)are solved using an iterative preconditioned conjugate gradient subspace algo- rithm, see Refs. 8, 50, and 52. The computationally most demanding step is the linear transformation Eq.(29) of the trial vectors.

To obtain the linear transformation, the matrixG([D0,bi]_S) must be calculated from the non-symmetric matrix [D0,bi]_S, wherebiis theith trial vector. Because of the non-symmetry, this step is slightly more expensive than the evaluation of the corresponding Kohn–Sham matrix. The two-electron integrals cannot be stored in memory but must be recalculated at each iteration, which is the computationally most expensive step of the procedure. However, since the integrals are independent of the perturbation, they can at each iteration be used for several Kohn–Sham matrices. Furthermore, the local expression for the right-hand side of response equation Eq.(33)ensures that the response solver can use sparse trial vectors, which means that the linear transformation Eq. (35) can exploit density matrix based integral screening techniques such as LinK,^53,54 allowing several Kohn–Sham matrices to be constructed with limited additional cost compared with a single Fock matrix construction. Consequently, the cost of the calculation scales sublinearly with the number of Kohn–Sham matrices and it becomes beneficial to solve equations for several perturba- tions simultaneously. Moreover, several vectors are then added

to the reduced space at each iteration, reducing the number of iteration needed to obtain convergence. Further benefits fol- low from a linear dependence between the solution vectors. In Ref.55, a similar strategy was used when solving for several frequencies simultaneously.

In principle, it is possible to construct approximate h^Bi where London orbitals are only used in the vicinity of the nucleus of interest, as in the GIAO truncation scheme. This would impose additional sparsity in the right-hand side, improving convergence and reducing cost.

E. Density-fitting

In applying Coulomb density fitting^41,42,45 in the eval- uation of NMR shielding constants, it is worth noting that density-fitting can be used straightforwardly by replacing the Coulomb matrix of some general matrixM

J_µ (M)=X

⇢

g_{µ ⇢} M _⇢ (36)

with the robustly density-fitted alternative^56,57 J˜_µ (M)=X

⇢

M _⇢*.

, X

P

(µ |P)C_⇢^P +X

Q

C^Q_µ (Q|⇢ ) X

PQ

C^Q_µ (Q|P)C^P_⇢ +/

-

, (37)

where (µ |P) is a general three-center two-electron repulsion integral

(µ |P)=

⌅⌅

µ(x₁) (x₁) 1

|x1 x2| ^P(x₂) dx₁dx₂, (38)

and where the fitting coefficients C_⇢^P =X

Q

(P|Q) ¹(Q|µ ) (39) are obtained by minimizing the self-interaction energy errors as available in LSDalton⁵⁶

µ ⇢ =X

Q

⇣(µ |⇢ ) C_µ^P (P|⇢ ) (µ |Q)C^Q_⇢

+C_µ^P (P|Q)C_⇢^Q ⌘

. (40)

Such a straightforward replacement is possible since the density-fitted Coulomb matrix depends onMin the same way as the non density-fitted Coulomb matrix,

FIG. 1. Structure of ↵-D-glucose (a) and amylose₆₄(b). The nucleus whose NMR shielding constant has been calcu- lated is labeled H1.

FIG. 2. Comparison of total B3LYP/6-31G* wall times in seconds for the nuclei-selected and standard NMR approach obtained via thesimultaneous solution of the response equations for nucleus H1.

@J(D(X))

@D(X) _D(X)=D0M=J(M), (41)

@J(D(X))˜

@D(X) D(X)=D0

M=J(M).˜ (42) The magnetic derivative density-fitted Coulomb matrix becomes

J˜^B_µⁱ(M)=X

⇢

M _⇢2666664 X

P

(µ |P)^BiC^P_⇢ +X

Q

C_µ^Q(Q|⇢ )^Bi

+X

P

*.

,

(µ |P) X

Q

C_µ^Q(Q|P)+/

-

@C^P_⇢

@Bi

+X

Q

@C^Q_µ

@B_i *

,(Q|⇢ ) X

P

(Q|P)C^P_⇢ + - 3777 775 . .

Note that Ref.47pointed out that, if GIAOs are used for auxil- iary fitting functions, they break gauge-origin invariance. We therefore do not use GIAOs for the auxiliary fitting functions,

J˜^B_µⁱ (M)=X

P

((µ )^Bi|P)X

⇢

C^P_⇢ M _⇢

+X

Q

X

⇢

C_µ^Q(Q|(⇢ )^Bi)M ⇢ (43)

making explicit the field independence of the auxiliary basis set.

In this case, M is the symmetric density matrix which reduces to

J˜^B_µⁱ (D)=X

P

((µ )^Bi|P)X

⇢

C^P_⇢ D _⇢. (44) Recently, Hollmann et al. presented a three-center integral screening technique (SQV`)⁵⁸that may be extended to mag- netic derivative three-center integrals to yield an improved efficiency of the density-fitting method.

While density fitting offers significant speed-ups at little loss of accuracy, the density-fitted Coulomb matrix construc- tion is not linear scaling due to the non-local fitting coefficients.

The computational time to solve the fitting equations in the Coulomb metric scales cubically with system size.

Baerendset al.⁵⁹fitted the electron density in the overlap metric, giving errors one order of magnitude greater than in the Coulomb metric. The coefficients obtained in the overlap metric decay more or less exponentially with distance, whereas the coefficients obtained in the Coulomb metric decay more slowly at long distances.

Using the overlap operator as the metric ˆw = (r1 r2) we may write the Coulomb matrix as

J˜_µ⌫(D)=

✓C_µ⌫^↵(↵|⇢ )+(µ⌫| )C⇢ +C_µ⌫^↵ (↵| )C⇢

C¯^↵_µ⌫h↵|wˆ|⇢ i+C¯^↵_µ⌫h↵|wˆ| iC_⇢ hµ⌫|wˆ| iC_⇢ +C^↵_µ⌫h↵|wˆ| iC_⇢ ◆

D_⇢ , (45) where the new local fitting coefficients and multipliers are determined from the linear equations

C¯^↵_µ⌫h↵|wˆ| i (µ⌫| )+C_µ⌫^↵ (↵| )=0, (46)

h↵|wˆ| iC¯_⇢ (↵|⇢ )+(↵| )C_⇢ =0. (47) Exploiting the fact that the auxiliary fitting functions have no magnetic dependence and that the density matrix is symmetric, the magnetic derivative Coulomb matrix becomes

J˜^B_µⁱ (D)=X

P

((µ )^Bi|P)X

⇢

C_⇢^P D _⇢ X

P

h(µ )^Bi|wˆ|PiX

⇢

C¯^P_⇢ D _⇢. (48)

TABLE I. Comparison of total wall times in seconds for the calculated B3LYP/6-31G* NMR shielding constants for↵-D-glucose and a set of amylose chains obtained via both thesimultaneousandnon-simultaneoussolutions of the response equations. The number of atoms is given in parentheses.

Nuclei-selected Standard NMR

Molecule Simultaneous Non-simultaneous Speed-up Simultaneous Non-simultaneous Speed-up

↵-D-glucose (24) 27 69 2.6 115 153 1.3

Amylose2(45) 112 198 1.8 627 744 1.2

Amylose₄(87) 344 586 1.7 253 6 288 7 1.1

Amylose8(171) 971 162 8 1.7 619 4 687 3 1.1

Amylose₁₆(339) 269 1 457 0 1.7 155 80 173 37 1.1

Amylose₃₂(675) 931 1 142 54 1.5 581 93 753 96 1.3

Amylose48(1011) 202 97 281 50 1.4 129 155 162 875 1.3

Amylose₆₄(1347) 380 54 460 18 1.2 204 193 264 887 1.3

TABLE II. Calculated B3LYP/6-31G*—percentage of non-zero numbers—of the right-hand side for ↵-D- glucose and a set of amylose chains obtained via thesimultaneoussolution of the response equations. The total number of regular and primitive regular basis functions is given as well.

Number of basis functions Nuclei-selected Standard NMR Molecule Regular Primitive regular Minimum Maximum Average Minimum Maximum Average

↵-D-glucose 192 372 99.9 99.9 99.9 100.0 100.0 100.0

Amylose₂ 366 709 99.9 99.9 99.9 99.9 100.0 99.9

Amylose4 714 138 3 93.9 95.0 94.4 99.9 99.9 99.9

Amylose₈ 141 0 273 1 54.7 56.5 55.5 85.8 87.7 87.0

Amylose16 280 2 542 7 30.7 34.8 33.2 53.8 55.5 54.9

Amylose₃₂ 558 6 108 19 9.2 12.5 11.4 29.5 31.6 30.9

Amylose48 837 0 162 11 4.9 7.6 6.7 20.3 22.2 21.6

Amylose₆₄ 111 54 216 03 3.0 5.2 4.5 15.4 17.2 16.6

III. COMPUTATIONAL DETAILS

The approach for the calculation of nuclei-selected NMR shielding tensors as described above has been imple- mented in a local version of the LSDalton program pack- age, including density-fitting and local density-fitting for the efficient calculation of Coulomb integrals as available in LSDalton.^60,61 Our implementation has been tested using a set of amylose chains, whose structures have been taken from Ref. 39. These systems have been used to test the linear scaling behavior in previous work.¹⁴ The molecular struc- ture of↵-D-glucose is depicted in Figure1; chains contain- ing up to 1347 atoms have been considered in the present study as also indicated in the figure. The notation amylosen

indicates that n amylose units are connected. The nucleus selected for the NMR shielding calculation is labeled H1 and highlighted.

The NMR shielding tensor calculations have been per- formed at the Kohn–Sham level of theory, using the hybrid Becke-three-parameter exchange functional combined with the Lee–Yang–Parr correlation functional (B3LYP)^62,63 and the Keal–Tozer KT3⁶⁴ functional, that has been developed specifically for the prediction of NMR shielding constants.

In all calculations, the 6-31G* basis set was employed.⁶⁵The performance of density-fitting was tested at the non-hybrid KT3/6-31G* level of theory with the recommended def2-SVP auxiliary basis.⁶⁶The threshold for the response solver was set to 10 ⁵times the square root of the number of electrons in the system.

IV. RESULTS AND DISCUSSION

All calculations have been performed for↵-D-glucose and a set of amylose chains with up to 1347 atoms, calculating in each case the shielding constant of H1, see Fig.1. Good agree- ment is observed between shielding constants calculated in the conventional manner and in the nuclei-selected manner (all values in the range 29.33–27.41 ppm), but the computational cost is significantly reduced for the nuclei-selected approach as shown in Figure2. A comparison of the numerical values of the NMR shielding constants for both approaches is given in Table S1 of thesupplementary material. In each case, the three response equations were solved in a simultaneous and separate manner, giving the same results.

In general, the time determining step is not the NMR step but the self-consistent-field (SCF) step. For example, for amylose16at the B3LYP/6-31G* level of theory, the SCF step takes about 6300 seconds on one node with sixteen central processing units (CPUs), whereas the nuclei-selected response solver takes about 2700 s. In the following, we discuss the NMR step in more detail.

The wall times for the NMR calculation reported in TableIand visualized in Figure2are obtained on one node with 16 CPUs. The simultaneous solution of the response equations is in all cases faster than the non-simultaneous solution. The speed-up is greater for the nuclei-selected approach due to the locality of the perturbation, which means that the solution of one response vector benefits from the trial vectors generated for another response. The speed-up for the simultaneous solu- tion arises mainly from the reduced overhead of calculating the three trial-vector transformations at once; it decreases for larger systems because of increased memory requirements.

In Figure 2, we have plotted the wall time against the number of atoms for the nuclei-selected and the conventional calculations of the H1 NMR shielding constant. The nuclei- selected approach is in all cases faster than the conventional one, decreasing much less steeply with increasing number of atoms.

FIG. 3. Comparison of the average of the percentage of non-zero num- bers in the right-hand side for the nuclei-selected and conventional NMR approaches employing thesimultaneous solutionof the response equations at the B3LYP/6-31G* level.

TABLE III. Calculated B3LYP/6-31G*—percentage of non-zero numbers—of trial vectors for the response equation (Equation(32)) for a set of amylose chains obtained via thesimultaneoussolution of the response equations. In addition, the number of iterations until the response convergence threshold of 10 ⁵is reached is given. The minimum, maximum, and average are reported for all the three components of the trial vector.

Nuclei-selected Standard NMR

Molecule Number of iteration Minimum Maximum Average Number of iteration Minimum Maximum Average

↵-D-glucose 4 100.0 100.0 100.0 5 100.0 100.0 100.0

Amylose₂ 4 99.9 99.9 99.9 6 99.9 99.9 99.9

Amylose4 4 98.7 98.9 98.9 6 99.8 99.8 99.8

Amylose₈ 4 75.8 77.2 77.3 6 88.7 91.9 90.7

Amylose16 4 39.1 45.8 42.9 6 56.6 62.0 58.6

Amylose₃₂ 4 17.6 22.7 20.9 7 31.0 32.6 32.0

Amylose48 4 10.8 15.1 13.3 7 22.2 23.6 22.7

Amylose₆₄ 4 7.6 11.0 9.6 7 15.9 16.8 16.5

The reduced computational time for the nuclei-selected NMR approach can be understood from the sparsity (here, the percentage of zero numbers) of the right-hand side of the response equations, which is much higher in the nuclei- selected approach, see TableIIand Figure3. As expected, the sparsity of the right-hand side increases with increasing system size.

In TableIII, the percentage of non-zero numbers of the response trial vectors is reported, in addition to the number of iterations needed for convergence, using the simultaneous solution. The minimum, maximum, and average matrix spar- sities are listed for all the three components of the trial vector.

The standard NMR equations take more iterations to converge than in the nuclei-selected approach—for example, amylose48

requires seven iterations with the standard approach but only four with the nuclei-selected approach. This is not surprising, bearing in mind that the trial vectors as well as the right- hand side of the response equations are more sparse in the nuclei-selected approach, by the locality of the nuclear mag- netic moment. The residual norm of the response trial vectors using the simultaneous solution is reported in Table S4 of the supplementary material.

Coulomb density-fitting, which has been employed for both the SCF and the NMR steps, affects the calcu- lated shielding constant by less than 1 ppb (parts per billion), see Table S8 of the supplementary material. In TableIV, the first (DF) column gives density-fitting timings for the nuclei-selected approach, respectively, whereas the

TABLE IV. Computation time in seconds for NMR shielding constants on B3LYP/6-31G*/def2-SVP level for↵-D-glucose and a set of amylose chains obtained via thesimultaneoussolution of the response equations. The label DF indicates that density-fitting has been employed.

Nuclei-selected Standard NMR

Molecule DF DF

↵-D-glucose 22 95

Amylose₂ 78 526

Amylose4 221 203 5

Amylose₈ 498 595 2

Amylose₁₆ 1209 145 51

Amylose32 3247 350 33

Amylose₄₈ 7412 713 36

FIG. 4. Comparison of total B3LYP/6-31G* wall times in seconds for the nuclei-selected NMR approach with and without density-fitting using the simultaneoussolution of the response equations for nucleus H1.

last column gives density-fitting timings for a conventional calculation. Comparing with Table I, we see that density- fitting in all cases reduces timings significantly. The same conclusion is drawn from Figure4, where timings with and without density-fitting are compared. Local density-fitting is available as a pilot implementation. Albeit local density- fitting has a larger error for the energy as compared to tradi- tional density-fitting, the calculated NMR shielding constants are with a maximum deviation of 0.04 ppm not strongly influenced.

TABLE V. Computation time in seconds for NMR shielding constants on KT3/6-31G*/def2-SVP level for↵-D-glucose and a set of amylose chains obtained via thesimultaneoussolution of the response equations.

Nuclei-selected

Molecule DF Without DF

↵-D-glucose 5 16

Amylose₂ 13 36

Amylose4 31 95

Amylose₈ 86 281

Amylose₁₆ 252 911

Amylose32 959 3263

Amylose₄₈ 2244 7182

Finally, we consider the non-hybrid KT3 functional, designed for accurate shielding calculations. As for the B3LYP functional, density-fitting significantly reduces timings for the NMR step (see TableV), while affecting the calculated shield- ing by less than 1 ppb (see Table S9 of the supplementary material). The KT3 and B3LYP shieldings differ by 0.3 ppm.

Since the Hessian is diagonal for non-hybrid functionals, the response equations are solved in one iteration.

V. SUMMARY AND CONCLUSIONS

The theory for the calculation of nuclei-selected NMR shielding tensors using London orbitals and molecular response theory has been developed and implemented in LSDalton, using Coulomb density-fitting and a simultaneous solution of response equations. A comparison of computa- tional timings between the non-simultaneous and simultane- ous solution of the response equations is presented showing that the simultaneous solution is superior and is thus recom- mended for future usage. As expected, density-fitting reduces computation time dramatically, while affecting the calculated shieldings by less than 1 parts per billion. The accuracy of the NMR shielding constants obtained with local density- fitting was with a maximum deviation of 0.04 ppm almost the same as with traditional density-fitting. The present imple- mentation does not include explicit field dependence of the auxiliary functions. A significantly improved but not linear- scaling behavior is observed for amylose chains containing up to 1347 atoms at the B3LYP/6-31G* level of theory.

For the calculation of a single shielding constant in a large molecule, a significant reduction in computing cost is observed in the nuclei-selected approach making the latter promising for future NMR studies including solvent effects⁶⁷and molecular dynamics.

SUPPLEMENTARY MATERIAL

Seesupplementary materialfor extra tables.

ACKNOWLEDGMENTS

This work was supported by the Norwegian Research Council through the CoE Centre for Theoretical and Com- putational Chemistry (CTCC) Grant Nos. 179568/V30 and 231571/F20. This work has received support from the Nor- wegian Supercomputing Program (NOTUR) through a grant of computer time (Grant No. NN4654K). The research lead- ing to these results has received funding from the European Research Council under the European Union’s Seventh Frame- work Programme (FP/2007-2013)/ERC Grant Agreement No.

291371.

1D. R. Hartree,Math. Proc. Cambridge Philos. Soc.24, 89 (1928).

2V. Fock,Z. Phys.61, 126 (1930).

3J. C. Slater,Phys. Rev.81, 385 (1951).

4W. Kohn and L. J. Sham,Phys. Rev.140, A1133 (1965).

5C. Ochsenfeld, J. Kussmann, and D. S. Lambrecht,Rev. Comput. Chem.

23, 1 (2007).

6J. Kussmann, M. Beer, and C. Ochsenfeld,Wiley Interdiscip. Rev.: Comput.

Mol. Sci.3, 614 (2013).

7P. Sałek, S. Høst, L. Thøgersen, P. Jørgensen, P. Manninen, J. Olsen, B. Jans´ık, S. Reine, F. Pawłowski, E. Tellgren, T. Helgaker, and S. Coriani, J. Chem. Phys.126, 114110 (2007).

8S. Coriani, S. Høst, B. Jans´ık, L. Thøgersen, J. Olsen, P. Jørgensen, S. Reine, F. Pawłowski, T. Helgaker, and P. Sałek,J. Chem. Phys.126, 154108 (2007).

9H. Larsen, P. Joørgensen, J. Olsen, and T. Helgaker,J. Chem. Phys.113, 8908 (2000).

10H. Larsen, T. Helgaker, J. Olsen, and P. Jørgensen,J. Chem. Phys.115, 10344 (2001).

11C. Ochsenfeld, J. Kussmann, and F. Koziol,Angew. Chem., Int. Ed.43, 4485–4489 (2004).

12K. Kristensen, P. Jørgensen, A. J. Thorvaldsen, and T. Helgaker,J. Chem.

Phys.129, 214103 (2008).

13M. Beer and C. Ochsenfeld,J. Chem. Phys.128, 221102 (2008).

14M. Beer, J. Kussmann, and C. Ochsenfeld,J. Chem. Phys.134, 074102 (2011).

15M. Maurer and C. Ochsenfeld,J. Chem. Phys.138, 174104 (2013).

16R. R. Ernst and W. A. Anderson,Rev. Sci. Instrum.37, 93 (1966).

17R. R. Ernst,J. Chem. Phys.45, 3845 (1966).

18H. Kessler, M. Gehrke, and C. Griesinger,Angew. Chem., Int. Ed.27, 490 (1988).

19T. Helgaker and P. Jørgensen,J. Chem. Phys.95, 2595 (1991).

20T. Helgaker, M. Jaszu´nski, and K. Ruud,Chem. Rev.99, 293 (1999).

21O. Vahtras, H. Ågren, P. Jørgensen, H. J. A. Jensen, S. B. Padkjær, and T. Helgaker,J. Chem. Phys.96, 6120 (1992).

22P. Pyykk¨o,Theor. Chem. Acc.103, 214 (2000).

23G. Schreckenbach and T. Ziegler,J. Phys. Chem.99, 606 (1995).

24J. Gauss,Chem. Phys. Lett.191, 614 (1992).

25J. Gauss and J. F. Stanton,Adv. Chem. Phys.123, 355 (2003).

26M. Kaupp, M. B¨uhl, and V. G. Malkin,Calculation of NMR and EPR Parameters(Wiley-VCH, Berlin, 2004).

27J. Autschbach,Principles and Applications of Density Functional Theory in Inorganic Chemistry I(Springer, Berlin, Heidelberg, 2004), pp. 1–48.

28J. Gauss, “Molecular properties,” inModern Methods and Algorithms of Quantum Chemistry, edited by J. Grotendorst (John von Neumann Institute for Computing, Julich, 2000), Vol. 3, pp. 541–592.

29M. H¨aser, R. Ahlrichs, H. P. Baron, P. Weis, and H. Horn,Theor. Chim.

Acta83, 455 (1992).

30M. Schindler and W. Kutzelnigg,J. Chem. Phys.76, 1919 (1982).

31M. W. Lodewyk, M. R. Siebert, and D. J. Tantillo,Chem. Rev.112, 1839 (2012).

32F. London,J. Phys. Radium8, 397 (1937).

33H. F. Hameka,Mol. Phys.1, 203 (1958).

34R. Ditchfield,Mol. Phys.27, 789 (1974).

35K. Wolinski, J. F. Hinton, and P. Pulay,J. Am. Chem. Soc.112, 8251 (1990).

36V. G. Malkin, O. L. Malkina, and D. R. Salahub,Chem. Phys. Lett.204, 80 (1993).

37A. M. Lee, N. C. Handy, and S. M. Colwell,J. Chem. Phys.103, 10095 (1995).

38T. Helgaker, P. J. Wilson, R. D. Amos, and N. C. Handy,J. Chem. Phys.

113, 2983 (2000).

39J. Kussmann and C. Ochsenfeld,J. Chem. Phys.127, 054103 (2007).

40C. Ochsenfeld and M. Head-Gordon,Chem. Phys. Lett.270, 399 (1997).

41B. I. Dunlap, J. W. D. Connolly, and J. R. Sabin,J. Chem. Phys.71, 3396 (1979).

42O. Vahtras, J. Alml¨of, and M. Feyereisen,Chem. Pyhs. Lett.213, 514 (1993).

43R. A. Kendall and H. A. Fr¨uchtl,Theor. Chem. Acc.97, 158 (1997).

44M. Feyereisen, G. Fitzgerald, and A. Komornicki,Chem. Phys. Lett.208, 359 (1993).

45K. Eichkorn, O. Treutler, H. ¨Ohm, M. H¨aser, and R. Ahlrichs,Chem. Phys.

Lett.240, 283 (1995).

46F. Weigend,Phys. Chem. Chem. Phys.4, 4285 (2002).

47S. Loibl, F. R. Manby, and M. Sch¨utz,Mol. Phys.108, 477 (2010).

48S. Loibl and M. Sch¨utz,J. Chem. Phys.137, 084107 (2012).

49T. Helgaker, P. Jørgensen, and J. Olsen,Molecular Electronic-Structure Theory(Wiley, Chichester, 2000).

50T. Kjærgaard, P. Jørgensen, A. J. Thorvaldsen, P. Sałek, and S. Coriani,J.

Chem. Theory Comput.5, 1997 (2009).

51T. Kjærgaard, P. Jørgensen, J. Olsen, S. Coriani, and T. Helgaker,J. Chem.

Phys.129, 054106 (2008).

52J. Kauczor, P. Jørgensen, and P. Norman,J. Chem. Theory Comput.7, 1610 (2011).

53C. Ochsenfeld, C. A. White, and M. Head-Gordon,J. Chem. Phys.109, 1663 (1998).

54C. Ochsenfeld,Chem. Phys. Lett.327, 216 (2000).

55J. Kauczor and P. Norman,J. Chem. Theory Comput.10, 2449 (2014).

56S. Reine, E. Tellgren, A. Krapp, T. Kjærgaard, T. Helgaker, B. Jansik, S. Høst, and P. Salek,J. Chem. Phys.129, 104101 (2008).

57B. I. Dunlap,J. Mol. Struct.: THEOCHEM501–502, 221 (2000).

58D. S. Hollman, H. F. Schaefer, and E. F. Valeev,J. Chem. Phys.142, 154106 (2015).

59E. J. Baerends, D. E. Ellis, and P. Ros,Chem. Phys.2, 41 (1973).

60K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O.

Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, E. K. Dalskov, U. Ekstr¨om, T. Enevoldsen, J. J. Eriksen, P. Ettenhuber, B. Fern´andez, L. Ferrighi, H. Fliegl, L. Frediani, K. Hald, A. Halkier, C. H¨attig, H. Heiberg, T.

Helgaker, A. C. Hennum, H. Hettema, E. Hjertenæs, S. Høst, I.-M. Høyvik, M. F. Iozzi, B. Jans´ık, H. J. Aa. Jensen, D. Jonsson, P. Jørgensen, J. Kauczor, S. Kirpekar, T. Kjærgaard, W. Klopper, S. Knecht, R. Kobayashi, H. Koch, J. Kongsted, A. Krapp, K. Kristensen, A. Ligabue, O. B. Lutnæs, J. I. Melo, K. V. Mikkelsen, R. H. Myhre, C. Neiss, C. B. Nielsen, P. Norman, J. Olsen,

J. M. H. Olsen, A. Osted, M. J. Packer, F. Pawlowski, T. B. Pedersen, P.

F. Provasi, S. Reine, Z. Rinkevicius, T. A. Ruden, K. Ruud, V. V. Rybkin, P. Sałek, C. C. M. Samson, A. S. de Mer´as, T. Saue, S. P. A. Sauer, B. Schim- melpfennig, K. Sneskov, A. H. Steindal, K. O. Sylvester-Hvid, P. R. Taylor, A. M. Teale, E. I. Tellgren, D. P. Tew, A. J. Thorvaldsen, L. Thøgersen, O. Vahtras, M. A. Watson, D. J. D. Wilson, M. Ziolkowski, and H. Ågren, Wiley Interdiscip. Rev.: Comput. Mol. Sci.4, 269 (2014).

61A linear scaling molecular electronic structure program, Release Dal- ton2015.1, LSDalton, 2015, seehttp://daltonprogram.org.

62A. D. Becke,J. Chem. Phys.98, 5648 (1993).

63C. Lee, W. Yang, and R. G. Parr,Phys. Rev. B37, 785 (1988).

64T. W. Keal and D. J. Tozer,J. Chem. Phys.121, 5654 (2004).

65W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys.56, 2257 (1972).

66F. Weigend,Phys. Chem. Chem. Phys.8, 1057 (2006).

67D. Flaig, M. Beer, and C. Ochsenfeld,J. Chem. Theory Comput.8, 2260 (2012).