An Efficient Method for Calculating Effective Core Potential Integrals

Subscriber access provided by Stockholm University Library

Article

An Efficient Method for Calculating Effective Core Potential Integrals Simon C. McKenzie, Evgeny Epifanovsky, Giuseppe M. J. Barca, Andrew T.B. Gilbert, and Peter M.W. Gill J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.7b12679 • Publication Date (Web): 21 Feb 2018 Downloaded from http://pubs.acs.org on February 26, 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.

The Journal of Physical Chemistry A 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 31 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 Journal of Physical Chemistry

An Efficient Method for Calculating Effective Core Potential Integrals. Simon C. McKenzie,† Evgeny Epifanovsky,‡ Giuseppe M.J. Barca,† Andrew T.B. Gilbert,† and Peter M.W. Gill∗,† †Research School of Chemistry, Australian National University, ACT 2601, Australia ‡Q-Chem Inc., 6601 Owens Dr, Pleasanton, CA 94588, USA E-mail: [email protected]

Abstract Effective core potential (ECP) integrals are amongst the most difficult one-electron integrals to calculate due to the projection operators. The radial part of these operators may include r0 , r−1 and r−2 terms. For the r0 terms, we exploit a simple analytic expression for the fundamental projected integral to derive new recurrence relations and upper bounds for ECP integrals. For the r−1 and r−2 terms, we present a reconstruction method that replaces these terms by a sum of r0 terms and show that the resulting errors are chemically insignificant for a range of molecular properties. The new algorithm is available in Q-Chem 5.0 and is significantly faster than the ECP implementations in Q-Chem 4.4, GAMESS (US) and Dalton 2016.

1

Introduction

Applying ab initio quantum chemical methods to heavy elements encounters two main challenges: large numbers of electrons and increasingly significant relativistic effects. Effective 1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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

core potentials (ECPs) partially address both problems by explicitly modelling only the valence electrons and by incorporating scalar-relativistic corrections. 1 The downside to using ECPs is that they introduce unprojected and projected one-electron integrals, the latter being difficult to evaluate. Several authors have developed methods to compute these problematic projected integrals. Barthelat et al. 2 derived analytic expressions by using the direct differentiation approach of Boys. 3 However, the resulting expressions rapidly become unwieldy, limiting their approach to low orders of angular momentum on the projector and basis set. The method of McMurchie and Davidson 4 is perhaps the most widely used due to its reliability and generality. They factor the integrals into angular and radial parts and treat the latter via asymptotic and power series expansions, which are relatively expensive to evaluate. In separate works, Kolar, 5 and Bode and Gordon 6 developed recurrence relations (RRs) for evaluating components of the radial part of the projected integral. However, a significant number of radial components are required for each projected integral making the method expensive. Pelissier et al. 7 expressed the projected integrals as linear combinations of overlap integral products but this approach introduces linear dependencies that lead to numerical instabilities. Skylaris et al. 8 and Flores-Moreno et al. 9 introduced adaptive quadratures to calculate the projected integrals but both the cost and accuracy of their methods can be difficult to predict. Hu et al. calculated ECP integrals over solid harmonic rather than cartesian Gaussian-type orbitals, computing the radial part via RRs. 10 Whilst the efficient calculation of integrals is desirable, it is even more important to avoid calculating negligible integrals. This can be achieved by using upper bounds but, curiously, there has been little work on developing rigorous upper bounds for either projected or unprojected ECP integrals. Formally, the total number of these integrals grows as O(M N 2 ), where M is the number of ECP centers and N is the number of basis functions. However, unless the basis functions and ECP are spatially close, the resulting integral will be negligible and, as a result, the number of significant ECP integrals grows as only O(M ). Song et al. recently

2


Page 2 of 31

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


derived an upper bound, but failed to demonstrate this ideal O(M ) scaling. Moreover, the computation of their bound required an irksome minimization over a numerical grid. 11 Shaw and Hill developed a rigorous bound for the radial part of ECP integrals but, once again, failed to demonstrate the ideal O(M ) scaling. 12 In the following Section we develop efficient recurrence relations for calculating both projected and unprojected ECP integrals. Our approach relies on the vast simplifications that occur when the potentials contain only r0 terms, in particular for the fundamental projected integral. 2 In Section 3 we present a powerful upper bound for these integrals that achieves O(M ) scaling and, in Section 4, we show how potentials with radial powers other than r0 can be handled via a reconstruction process. Finally, in Section 5 we outline an algorithm that couples these three ideas together, and in Section 6 we present comparative timing results for a slab of platinum atoms.

2

Recurrence Relations

The unnormalized primitive Gaussian-type orbital

|a] = a(r) = (x − Ax )ax (y − Ay )ay (z − Az )az exp(−α|r − A|2 )

(1)

is defined by its exponent α, center A = (Ax , Ay , Az ), angular momentum vector a = (ax , ay , az ) and total angular momentum a = ax + ay + az . A primitive shell, |a], is written in non-bold type and comprises all Gaussian basis functions of the same total angular momentum, exponent and center. A contracted Gaussian-type orbital is a linear combination of primitives

|ai =

Ka X

Di |a]i

(2)

i=1

where Ka is the degree of contraction and Di is a contraction coefficient. Without loss of 3



Page 4 of 31

generality, we can assume that the ECP is centered at the origin and has the form

U(r) = UL (r) +

L−1 X ` X

|Y`m i U` (r) hY`m |

(3)

`=0 m=−`

where UL (r) and U` (r) are radial potentials that have the form

U (r) =

K X

Dk rnk exp(−ηk r2 )

(4)

k=1

ηk is the ECP exponent and

P

m

|Y`m i hY`m | is the spherical harmonic projector of angular

momentum `. Most commonly, nk = −2, −1 or 0 and L rarely exceeds 5. For the remainder of the paper we are concerned with primitive integrals and, therefore, drop the contraction index k, i.e. n ≡ nk . Primitive matrix elements of this ECP in the Gaussian basis are constructed from unprojected integrals Z [a|b]L =

a(r) rn exp(−ηr2 ) b(r) dr

(5)

X

(6)

and projected integrals Z [a|b]` = 4π

∞

r2+n exp(−ηr2 )

0

ha|Y`m i hY`m |bi dr

m

A class of integrals is written in non-bold type, i.e. either [a|b]L or [a|b]` , and comprises all integrals over the basis functions in shells |a] and |b]. These are built up from the momentumless fundamental integrals [0|0]L and [0|0]` . 13 Contracted analogues of these quantities are written using angle brackets.

2.1

Unprojected Integrals

The trivial identity (t − Bt ) = (t − At ) + (At − Bt ) ,

4


(7)

Page 5 of 31 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 t ∈ {x, y, z}, yields the horizontal recurrence relation (HRR) 14

ha|b + 1t iL = ha + 1t |biL + (At − Bt ) ha|biL

(8)

where 1t is the unit 3-vector in the tth direction. The HRR shifts angular momentum from center A to center B and, therefore, we need only a vertical recurrence relation (VRR) that builds angular momentum on center A. The unprojected fundamental integral is the three-center overlap Z

rn exp(−ηr2 − α|r − A|2 − β|r − B|2 ) dr n+3 n 3 2 − n+3 2 2 2 = 2π Γ ζ 2 exp(−αA − βB + ζR ) 1 F1 − , , −ζR 2 2 2

[0|0]L =

(9)

where ζ = α + β + η and R = |R| = |(αA + βB)/ζ| are the exponent and center of the product gaussian and 1 F1 is the confluent hypergeometric function. 15 We have derived VRRs for general n, but we have not pursued these because they are much more complicated than the n = 0 case, where the hypergeometric reduces to unity, the integral factorizes, and the subsequent recursion simplifies dramatically. For n = 0, one finds [0|0]L = (π/ζ)3/2 exp(−αA2 − βB 2 + ζR2 )

(10)

and the higher [a|0]L satisfy the two-term VRR

[a + 1t |0]L = (Rt − At )[a|0]L +

at [a − 1t |0]L 2ζ

(11)

which is a simplification of the McMurchie-Davidson one-center RR for two electron repulsion integrals. 16 A near-optimal strategy for its use is known. 17

5



2.2

Page 6 of 31

Projected Integrals

Unfortunately, the presence of the spherical harmonic projectors means the HRR cannot be used with projected ECP integrals. Hence, we require a VRR to build angular momentum on center A as well as a subsequent VRR for building angular momentum on center B. The angular integrations in [0|0]` are elementary and, upon performing these, one obtains 4 2

2

[0|0]` = exp(−αA − βB )(2` + 1)P`

A·B AB

Z

∞

4π

r2+n exp(−ζr2 )i` (2αAr) i` (2βBr) dr

0

(12) where P` is a Legendre polynomial and i` is a modified spherical Bessel function of the first kind. 15 For most values of n, the radial integral in (12) leads to awkward combinations of confluent hypergeometric functions but, in the n = 0 case, it can be shown that [0|0]` = G i` (T ) (2` + 1)P`

A·B AB

(13)

where G = (π/ζ)3/2 exp(−α0 A2 − β 0 B 2 ) T =

2αβAB , ζ

(14) (15)

α0 = α(ζ − α)/ζ and β 0 = β(ζ − β)/ζ. This result was discovered by Barthelat et al. 2 in 1977, but has not received the attention that it deserves, except for Hu et al. 10 It is convenient to generalize (13) by defining the three-index auxiliary fundamental (m,p,q)

[0|0]`

= G A2p B 2m Im+p+q+` L`

(16)

In = T −n in (T )

(17)

where

6


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


n

Ln = (2n + 1)T Pn

A·B AB

.

(18)

In and Ln satisfy the two-term RRs

In−1 = (2n + 1)In + T 2 In+1

Ln+1

2n + 3 = T n+1

A·B n Ln − T Ln−1 AB 2n − 1

(19) (20)

and the fundamental integral (16) enjoys the derivative property 2αβ 2 At pAt α−ζ ∂ (m,p,q) (m,p+1,q) (m,p−1,q+1) (m,p,q) At [0|0]` + [0|0]` + [0|0]` [0|0]` = 2 2α∂At ζ ζ α " b(`−1)/2c 2k X (2` + 1)β 2αβ (m+k,p+k,q+1) + Bt [0|0]`−2k−1 ζ ζ k=0 # 2k+1 b(`−2)/2c X 2αβ (m+k+1,p+k,q+1) [0|0]`−2k−2 (21) − At ζ k=0 Using (16) and (21), Ahlrichs’ method 18 yields the 6 + 2` term VRR (m,p,q)

[a + 1t |b]`

α−ζ αAt 2β 2 pAt (m,p,q) (m+1,p,q) (m,p−1,q+1) At [a|b]` + [a|b]` + [a|b]` ζ ζ ζ α at at β 2 p at (m,p,q) (m+1,p,q) (m,p−1,q+1) + [a − 1t |b]` + 2 [a − 1t |b]` + 2 [a − 1t |b]` 2ζ ζ 2α " b(`−1)/2c 2k X (2` + 1)β 2αβ (m+k,p+k,q+1) + Bt [a|b]`−2k−1 ζ ζ k=0 2k+1 b(`−2)/2c X 2αβ (m+k+1,p+k,q+1) − At [a|b]`−2k−2 ζ k=0 2k b(`−1)/2c bt X 2αβ (m+k,p+k,q+1) + [a|b − 1t ]`−2k−1 2β k=0 ζ # 2k+1 b(`−2)/2c 2αβ at X (m+k+1,p+k,q+1) − [a − 1t |b]`−2k−2 (22) 2α k=0 ζ =

which is a key result of this Paper. This can be used, together with the analogous RR

7



Page 8 of 31

for building momentum on center B, to construct the target class of projected integrals (0,0,0)

[a|b]` ≡ [a|b]`

.

We note that these RRs relate integrals over different projectors; e.g. a d-projected class is constructed from d-, p- and s-projected classes. It is also worth noting that, although our RR has three auxiliary indices, the index decrements in the third and sixth terms cause many of the intermediates to disappear, unlike other multi-index RRs. 19

3

Upper Bounds

If we can show that an ECP integral is smaller than an upper bound that, in turn, is smaller than a cutoff threshold τ , the ECP integral can be safely ignored. It follows that the ideal upper bound is both strong (i.e. similar in magnitude to the integral) and simple (i.e. much cheaper than the integral). 20 For a given ECP and shell ha|, it is possible to form a two-center upper bound hB2 i for all possible ECP classes ha|bi and this can be used to eliminate many ECP-shell pairs. Similarly, for a given ECP, shell ha| and shell |bi, one can form a (tighter) three-center upper bound hB3 i for the particular class ha|bi and this can be used to eliminate many ECP-shell-shell triples. This two-layer strategy is analogous to the shell-pair and shell-quartet screening used for electron repulsion integrals. 13 We derive our two- and three-center upper bounds using shell-bounding gaussians (SBGs) 21 and this allows us to exploit the simple analytic expressions for n = 0 fundamental integrals.

3.1

Shell-Bounding Gaussians

An SBG ˚ a is an s-type gaussian that bounds all of the functions in a primitive shell |a], i.e.

˚ a(r) ≥ ||a]| ∀|a] ∈ |a]

8


(23)

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


If the shell is unnormalized, one can show that

˚ a(r) = Na exp(−˜ α |r − A|2 ) a/2 a Na = 2eασa

(24) (25)

where α ˜ = (1 − σa )α is an effective exponent and σa is an arbitrary parameter whose value can be chosen to minimize the resulting upper bound. 21 Because SBGs bound an entire shell of basis functions, they lead naturally to more efficient class bounds rather than integral bounds. In the following Sections, we present key formulae for a variety of upper bounds. Additional details can be found in Appendix A.

3.2 3.2.1

Unprojected Integrals Two-Center Bounds

The largest integral |[a|b]L | = max |[a|b]L |

|a] ∈ |a], |b] ∈ |b]

(26)

in an unprojected integral class possesses the two-center class bound

|[a|b]L | ≤ [B2 ]L = [˚ a]L Q

(27)

α ˜ η A2 [˚ a]L = Na exp − α ˜+η

(28)

where

Q=

(ˆ a + 3)aˆ+3 (π/3)3 (2e)aˆ α ˇ aˆ+3

1/2

a ˆ is the largest angular momentum and α ˇ is the smallest exponent in the basis set.

9


(29)


Page 10 of 31

The optimal σ minimizes [˚ a]L and exactly satisfies ∂[˚ a]L =0. ∂σ

(30)

However, for computational efficiency, we use the inexpensive but accurate approximation

σa =

a(α + η)2 . 2α(η 2 A2 + a(α + η))

(31)

Only shells for which [˚ a]L ≥ τ /Q need to be retained for subsequent processing and, in this way, most of the shells that are far from the ECP are eliminated.

3.2.2

Three-Center Bounds

An unprojected integral class possesses the tighter three-center class bound

|[a|b]L | ≤ [B3 ]L = Na Nb [˚ a˚ b]L

(32)

˜ 3/2 exp −˜ ˜ 2 + ζ˜R ˜2 [˚ a˚ b]L = (π/ζ) αA2 − βB

(33)

where

˜ ˜ Ideally, σa and σb would be obtained by ˜ = |(˜ ˜ = |R| ζ˜ = α ˜ + β˜ + η and R αA + βB)/ ζ|. minimising [B3 ]L . However, the resulting expression is too expensive to be used as a bound, so we use the values obtained from the two-center bound (i.e. eq (31)). The primitive [B3 ]L bound can be naively extended to obtain a contracted bound for the unprojected integrals

| ha|biL | ≤

Ka X i

|Di |

Kb X

|Dj |

j

KL X

|Dk |([B3 ]L )ijk ,

(34)

k

but this requires O(K 3 ) computational work. An O(K) contracted upper bound can be

10


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


obtained by following the approach of Barca et al. 21

b | ha|biL | ≤ hB3 iL = h˚ ai h˚ bi h˚ a˚ biL

(35)

where h˚ ai =

Ka X

|Di |Nai

(36)

i

KL b ˇ˜ 2 X |D |(π/η )3/2 exp ζˇ˜ R ˇ˜ 2 ˇ˜ A2 − βB h˚ a˚ biL = exp −α k k k k

(37)

k

ˇ˜ | = ˇ˜ = |R ˇ˜ + βˇ˜ + ηk and R ˇ˜ is the smallest effective exponent in the contracted shell, ζˇ˜k = α α k k ˇ˜ ˇ˜ A + βB)/ ζˇ˜k |. The optimal σa in Nai of (36) is calculated using the smallest exponent in |(α the contracted ECP.

3.3

Projected Integrals

To construct upper bounds for projected integral classes, we begin by bounding the projector by its absolute value and, subsequently, by a scaled s-projector

|Y`m i ≤ ||Y`m |i ≤

3.3.1

p

(2` + 1) |Y00 i .

(38)

Two-Center Bounds

A projected integral class possesses the two-center class bound

|[a|b]` | ≤ [B2 ]` = (2` + 1)2 [˚ a]` Q

where [˚ a]` ≡ [˚ a]L is given by eq (28) and Q is defined in eq (29).

11


(39)


3.3.2

Page 12 of 31

Three-Center Bounds

A projected integral class possesses the tighter three-center class bound

|[a|b]` | ≤ [B3 ]` = (2` + 1)2 Na Nb [˚ a˚ b]`

(40)

˜ 2+α [˚ a˚ b]` = (π/η)3/2 exp −˜ αA2 − βB ˜ 2 A2 /ζ˜ + β˜2 B 2 /ζ˜ f (T )

(41)

where

f (T ) =

   exp(T )/(2T ) T > 1   sinh(1)

(42)

T ≤1

˜ 2˜ αβAB T˜ = ζ˜

(43)

and we have made use of the inequality i0 (x) ≤ exp(x)/(2x). The σa and σb that appear in Na and Nb are determined independently by considering the interaction of each SBG with the ECP projector exp(−ηr2 ) |Y00 i and, because Y00 is spherically symmetric, they are given by (31). As in the unprojected case, a contracted three-center bound for the projected class is

c |ha|bi` | ≤ hB3 i` = (2` + 1)2 h˚ ai h˚ bi h˚ a˚ bi`

(44)

where K` c ˇ˜ 2 X |D |(π/η )3/2 exp α ˇ˜2 B 2 /ζ˜ f (Tˇ˜ ) 2 2 ˜ ˇ˜ A2 − βB ˇ h˚ a˚ bi` = exp −α ˜ A / ζ + β k k k k k

(45)

k

ˇ˜ ˇ˜ βAB 2α ˇ ˜ Tk = ζˇ˜k

12


(46)

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


3.4

Performance

To investigate the performance of the new upper bounds, we considered slabs of platinum atoms, a single unit cell thick and of varying sizes, (Figure 1) with the SRSC 22–28 ECP and basis set. As can be seen in Figure 2, the number of ECP centers, M , increases linearly

Figure 1: Pt36 slab constructed from face-centered unit cells. with the size of the Pt slab and the total number of ECP classes grows as O(M N 2 ) where N is the number of basis functions. Asymptotically, the number of significant ECP classes increases only linearly, and this ideal O(M ) scaling can already be observed for Pt256 , which has a maximum dimension of 27.5 ˚ A and 104 basis functions. The three-center screening overestimates the number of significant classes by a factor of approximately four, but it still has the correct O(M ) scaling. The two-center screening is weaker and overestimates the number of significant classes by a factor of about 25. It also predicts O(M ) significant classes, although it approaches this scaling more slowly.

4

ECP Reconstruction

The efficient method to evaluate ECP projected integrals outlined in Sections 2 and 3 applies only to radial potentials U (r) where n = 0. Whilst this is the case for some ECPs, such as those in the Stuttgart-Cologne family, 22–28 many ECPs have radial potentials with 13



1013

y ∝ M3.00 y ∝ M2.08 y ∝ M1.19 y ∝ M1.09 y ∝ M1.07

● No Screening ▲ Shell-Pair Screening ◆ ECP-Shell Screening

Number of Significant ECP Classes

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

■ ECP-Shell-Shell Screening 10

11

Page 14 of 31

▼ Ideal

●

● ▲ ● ●

109 ● ▲ ◆ ■

● 107

105 10

● ▲ ◆ ■ ▼

▲ ◆ ■ ▼

●

●

▲ ◆

● ▲ ◆

■

■

▼

▼

▲

▲

◆

◆

◆

■

■

■

▼

▼

▼

◆

▲ ◆

■

■

▼

▲

▼

▼

50

100 Number of ECP Centers (M)

500

1000

Figure 2: Scaling of significant ECP classes with increasing size of the Pt slab. The SRSC ECP was used and the threshold set to 10−8 . n = −2 and −1 terms, such as LANL2DZ, 29–31 HWMB, 29 SBKJC, 32–34 CRENBL 35–40 and CRENBS. 36,37,40,41 Examples of these ECP radial potentials for the d-projector on a Pt atom are shown in Figure 3. The significant variation between these radial potentials at small r is striking, and highlights the non-uniqueness of ECPs. 42 Despite this variation, all of these ECPs produce similar chemistry, indicating the resulting valence orbitals are insensitive to the behaviour of the ECP potential at small r. 43 Furthermore, the performance of the SRLC and SRSC ECPs, which only include n = 0 terms, indicates that the inclusion of n = −2 and −1 terms is by no means essential. In this Section we propose a method for reconstructing the n = −2 and −1 radial potentials using only n = 0 terms to allow our recurrence relations and bounds to be applied more generally.

14


Page 15 of 31 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


r 2 U2 (r) SRSC 14

LANL2DZ CRENBS

12

SBKJC

10 8 6 4 2 0.0

0.2

0.4

0.6

0.8

1.0

r (a.u.)

Figure 3: SRSC, LANL2DZ, CRENBS and SBKJC radial potentials for the d-projector, U2 (r), on Pt.

4.1

Method

To model a generic n = −2 term by a sum of three n = 0 terms, we minimized the integral Z

∞

" r−2 exp(−r2 ) −

Z2 = 0

3 X

#2 w(r)r2 dr

ck exp(−λk r2 )

(47)

k=1

and chose the weight function w(r) = r4 , which gave the best fit to the chemically important regions of those we tested. The optimal scale factors that resulted were λ1 = 1.5, λ2 = 5.4 and λ3 = 30.5. To model a generic n = −1 term by a sum of three n = 0 terms, we minimized the integral Z Z1 = 0

∞

" r−1 exp(−r2 ) −

3 X

#2 dk exp(−µk r2 )

w(r)r2 dr

(48)

k=1

with the same weight function. The optimal scale factors that resulted were µ1 = 1.2, µ2 = 2.8 and µ3 = 11.6. Armed with these results, we can reconstruct an ECP using a three-step protocol: 1. For each n = −2 term with exponent η, form n = 0 terms with exponents λ1 η, λ2 η, 15



λ3 η. 2. For each n = −1 term with exponent η, form n = 0 terms with exponents µ1 η, µ2 η, µ3 η. 3. Combine the n = 0 terms from steps 1 and 2 with any n = 0 terms in the original potential and use these as a basis to least-squares fit the original radial potential with the same weight function. The fitting process is straightforward and performed at run-time to allow any ECP with n = −2 and n = −1 terms to be reconstructed. Linear dependence is uncommon but, to guard against numerical problems, the least-squares equations are solved using the MoorePenrose pseudo-inverse.

4.2

Reconstructed ECP Performance

Figure 4 shows an example of a reconstructed SBKJC ECP using the above procedure. The original SBKJC d-projector for Pt has one n = 0 term, and one n = −2 term. The latter yields three new n = 0 terms, resulting in a four-fold contracted potential. We choose to show an SBKJC ECP because this family of ECPs is challenging to reconstruct. This is due to their small degree of contraction, K, which results in less flexibility for the fit. For r ≤ 0.1, the reconstructed ECP is qualitatively incorrect but this region is chemically unimportant. For r > 0.1, the reconstruction is very accurate and the potentials become essentially indistinguishable. The differences are smaller than those between potentials from different ECPs, as shown in Figure 3. To test further the effects of the reconstruction procedure, we performed energy, geometry and frequency calculations on several diatomic molecules using the original and reconstructed ECPs. The reconstruction changes the total energy of an atom by an average of 2 mEh. However, because this arises from the inner core region, it cancels almost perfectly when relative energies are calculated. 16


Page 16 of 31

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


r 2 U2 (r) 14 SBKJC Reconstructed SBKJC

12

10

8

6

4

2

0.2

0.4

0.6

0.8

1.0

r (a.u.)

Figure 4: Comparison of the SBKJC radial potential for the d-projector U2 (r) on Pt with its reconstructed form.

Table 1: Absolute differences in dissociation energies (∆De / kJ mol−1 ), equilibrium bond lengths (∆re / pm) and vibrational frequencies (∆ωe / cm−1 ) between original and reconstructed LANL2DZ, SBKJC and CRENBL ECPs

TaO PtO AuCl TlF

LANL2DZ ∆De ∆re ∆ωe 0.01 0.001 0.10 0.07 0.004 0.23 0.02 0.003 0.01 0.00 0.000 0.00

∆De 1.65 0.03 0.64 0.56

17

SBKJC ∆re 0.064 0.352 0.417 0.110

∆ωe 0.68 3.57 1.51 0.21


CRENBL ∆De ∆re ∆ωe 0.01 0.001 0.01 0.03 0.012 0.15 0.10 0.026 0.08 0.02 0.002 0.06


All calculations were run at the HF level of theory and used the 3-21G basis set on non-metal atoms. The results are shown in Table 1. The mean deviations between the original and reconstructed ECPs for dissociation energies, equilibrium bond distances and fundamental vibrational frequencies are 0.3 kJ mol−1 , 0.08 pm and 0.6 cm−1 , respectively. These deviations are orders of magnitude smaller than the differences between different ECPs, and are of the same order of magnitude as those introduced by the quadrature methods used by Flores-Moreno et al. 9 In this sense, we claim that the ECP reconstructions are accurate.

5

Computational Details

The present method has been implemented within Q-Chem’s new integrals library, libqints, replacing the previous ECP implementation. 44 libqints provides a general framework to both compute and validate molecular integrals and, because validation is a key part of our approach, we have provided suitable reference data in the supplementary material. We first compute and store the basis constant Q, and the h˚ ai and h˚ aiL upper bound factors for each contracted shell. These are computed over batches of similar basis function and ECP projector angular momentum pairs (a, `) to allow for vectorization and efficient memory usage. 13 As with the calculation of the upper bound factors, the ECP integrals are also processed in batches of similar shell and projector angular momenta, (a, `, b). For each batch, we loop over the ECP centers, shells |ai and shells |bi in the batch. Inside these loops the ECP integral screener forms a list of significant ECP-shell pairs using either the hB2 iL or hB2 i` bounds for unprojected and projected integrals, respectively. These significant ECP-shell pairs are combined to form a list of ECP-shell-shell triples, which are further screened using either the hB3 iL or hB3 i` bounds. The advantage of this approach is the bulk of insignificant integrals are screened by the less expensive ECP-shell

18


Page 18 of 31

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


screener before proceeding to the stronger, but more expensive, ECP-shell-shell screener. The surviving ECP-shell-shell triples are ordered so that a ≥ b in preparation for the VRRs. The computation of the ECP integrals is divided into unprojected and projected cases. For the unprojected integrals, we first construct the required fundamental integrals (10), then build angular momentum on center A using the VRR (11), then contract the integrals and finally shift angular momentum from center A to B using the HRR (8). The construction of the projected fundamental integrals and the application of the VRR (22) is relatively complicated. At runtime, we create drivers that determine the small set of intermediates required to build the desired integrals from the auxiliary fundamental integrals (16). A diagram of the [d|0]d driver is shown in Figure 5 and shows how these intermediates are combined to form the target integral class. These drivers are created for each class, (a, `, b), of projected integrals and our implementation can handle arbitrary ECP projector and basis function angular momenta. We compute the factors (G, Ln , Im ) for the fundamental integrals separately, with the Ln and Im factors being computed via the RRs (19) and (20). These factors are combined to build the range of required fundamental projected integrals, (16). The projected VRR (22) is executed for each step of the driver, first building angular momentum on center A, and then on center B using the corresponding VRR. Finally, the unprojected and projected integrals from each ECP center are accumulated into the target ECP matrix element.

6

Results and Discussion

To evaluate the efficiency of the present ECP integral method, timings were carried out on the same Pt slab systems from Section 3 using the new method as implemented in QChem 5.0, the old ECP method in Q-Chem 4.4, 44 GAMESS (US) 6,45 and Dalton 2016. 46 All calculations were performed on a 2.6 GHz Intel Xeon E5-2690 v4 platform using a single CPU core.

19



[d|0]d (0,0,0)

[0|0]d (0,0,0)

[0|0]d (1,0,0)

[p|0]d (0,0,0)

[p|0]s(1,0,1)

[p|0]d (1,0,0)

[0|0]s(1,0,1)

[0|0]p (0,0,1)

[0|0]s(2,0,1)

[p|0]p (0,0,1)

[0|0]d (2,0,0)

[0|0]p (1,0,1)

[0|0]s(0,0,2)

Figure 5: A diagram showing how auxiliary classes are combined to form a [d|0]d class of integrals. Colors indicate projector momentum.

106

105

t ∝ M2.70 ▼ Present t ∝ M1.68

● Dalton

● 104 CPU Time (s)

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

Page 20 of 31

●

● ●

1000

▼ ▼

100

●

▼

▼

▼ ▼

10 ▼ 1 10

20

50 100 Number of ECP Centers (M)

200

500

Figure 6: ECP integral timings for Pt slabs using the SRSC ECP and an integral threshold of 10−8 .

20


Page 21 of 31

Figure 6 shows timing data for Pt slabs with SRSC ECPs. Q-Chem 4.4 and GAMESS (US) cannot handle the g-projectors in this ECP, so only Dalton 2016 and the present method are shown. It is immediately obvious that the present method is around two orders of magnitude faster than Dalton across all system sizes considered. Dalton does not employ any screening and therefore calculates the full O(M N 2 ) set of ECP integrals, as demonstrated by the asymptotic scaling of its timings O(M 2.70 ). Whilst the pair/triple screener predicts the ideal O(M ) number of integrals, the present method does not achieve this ideal asymptotic scaling in timings O(N 1.68 ) due to the small O(M N ) cost of the screener and the O(M N 2 ) overhead of looping over the batches. 105

104 CPU Time (s)

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


t ∝ M2.22 ▲ GAMESS (US) t ∝ M2.82 t ∝ M2.04 ■ Q-Chem 4.4 t ∝ M1.56 ▼ Present

● Dalton

▲ ●

● ▲

■

● ▲ ■

● 1000

▲ ■

● ▲

100

10

▼ ■ ▼

●

■

▲

▼

■

▼

▼

▼

■ ▼ 1 10

20

50 100 Number of ECP Centers (M)

200

500

Figure 7: ECP integral timings for Pt slabs using the SBKJC ECP. The equations show the results of linear regressions performed on the last few data points. Calculations were performed using the original SBKJC ECP except for the present method, which used the reconstructed SBKJC ECP. The integral threshold was set to 10−8 . Figure 7 shows timing data for the same Pt slabs except this time using the SBKJC ECP, which includes a maximum of f -projectors on Pt. Q-Chem 4.4, GAMESS (US) and Dalton are able to calculate ECP integrals where n 6= 0 and, therefore, the original Pt SBKJC ECP, with K = 1–3 was used. The present method used the reconstructed Pt SBKJC ECP with K 0 = 3–5. Despite the greater number of ECP primitives, the present method 21



is significantly faster than Q-Chem 4.4, GAMESS (US) and Dalton across all system sizes considered. GAMESS (US) and Dalton employ no screening, where as Q-Chem 4.4 uses shellpair screening and the present method uses pair/triple screening. The asymptotic scaling of the timings for each method reflects these screening strategies.

7

Conclusion

We have presented an efficient method for computing ECP integrals over Gaussian basis functions. The key to this new method is restricting the ECP radial potentials, U (r), to contain only r0 terms. This allows for a simple analytic form for the otherwise difficult projected fundamental integral. Using this analytic form, we have derived the first RRs in the literature to directly build angular momentum on the Gaussian basis functions. These RRs can handle arbitrary ECP projector and basis function angular momenta. We have also developed powerful upper bounds to screen ECP-shell pairs and ECP-shellshell triples using SBGs. 21 These new bounds are able to reduce the number of ECP integrals from formally cubic scaling to the ideal linear scaling. In order to make our method more widely applicable, we presented a procedure that converts ECPs containing r−2 and r−1 terms into reconstructed ECPs that contain only r0 terms. The discrepancies introduced by this reconstruction were shown to be chemically insignificant and similar in magnitude to those introduced by quadratures used in other ECP methods. 9 Our method has been implemented in the Q-Chem 5.0 package, and this was used to show our approach significantly reduces the cost of calculating ECP integrals compared to Q-Chem 4.4, GAMESS (US) and Dalton 2016. In light of the advantages of our method, we recommend only r0 terms be included in future ECPs.

22


Page 22 of 31

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


Acknowledgement S.C.M. thanks the Westpac Bicentennial Foundation for a Future Leaders scholarship and the Australian Government for an Australian Government Research Training Program scholarship. P.M.W.G. thanks the Australian Research Council for funding (Grants No. DP140104071 and DP160100246) and the NCI National Facility for generous grants of supercomputer time.

8

Supporting Information

The supporting information contains reference ECP integral values for a variety of classes.

9

Appendix A

To derive our upper bounds, we exploit the absolute value inequality (AVI) Z Z f (r) dr ≤ |f (r)| dr

(49)

Z Z 1/q 1/p Z p q f (r)g(r) dr ≤ |g(r)| dr |f (r)| dr

(50)

and Hölder’s inequality 15

where p−1 + q −1 = 1.

9.1

Unprojected Integrals

Using the AVI and SBGs, we can bound a primitive unprojected integral by

|[a|b]L | ≤ [˚ a| exp(−ηr2 )|˚ b]

23


(51)


9.1.1

Page 24 of 31

Two-center Bounds

Applying Hölder’s inequality to (51), and taking the limit as p → ∞ and q → 0, gives the unprojected upper bound |[a|b]L | ≤ [˚ a]L [˚ b]

(52)

˜ 3/2 . [˚ b] = Nb (π/β)

(53)

where [˚ a]L is given by (28) and

The σb that minimizes [˚ b] is σb =

b . b+3

(54)

Maximizing [˚ b] with respect to b and β yields the Q factor (29) which is independent of |b].

9.1.2

Three-center Bounds

Using the analogous fundamental unprojected integral (10) solution, (51) gives the threecenter bound (32) for an unprojected class.

9.2

Projected Integrals

The AVI, SBGs and the inequality |Y`m | ≤ Z

p

(2` + 1)Y00 yield the projected integral bound

r2 exp(−ηr2 )

0

Z ≤ 4π

` X

∞

|[a|b]` | ≤ 4π

|ha|Y`m i| |hY`m |bi| dr

m=−` ` X

∞ 2

2

r exp(−ηr ) 0

h|a|||Y`m |i h|Y`m |||b|i dr

m=−`

≤ 4π (2` + 1)

2

Z

∞

r2 exp(−ηr2 ) h˚ a|Y00 i hY00 |˚ bi dr

0

24


(55)

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


9.2.1

Two-center Bounds

Applying Hölder’s inequality to (55) gives the following projected upper bound

|[a|b]` | ≤ (2` + 1)2 [˚ a]L [˚ b]

(56)

where [˚ a]L and [˚ b] are defined as in the unprojected two-center bound derivation, (28) and (53), respectively. Following the same procedure leads to the two-center bound for a primitive class of projected integrals (39).

9.2.2

Three-center Bounds

Applying the simple projected fundamental integral solution (13) to (55) gives the threecenter bound (40) for a projected class

˜ 2+α ˜ i0 (T ) [˚ a˚ b]` = (π/η)3/2 exp(−˜ αA2 − βB ˜ 2 A2 /ζ˜ + β˜2 B 2 /ζ)

(57)

where T is defined in (43). Exploiting the inequality

i0 (x) ≤

   exp(x)/2x x > 1   sinh(1)

(58)

x≤1

allows us to bound (57) and avoid evaluating the expensive modified spherical Bessel function, as seen in (41).

References (1) Cundari, T. R.; Benson, M. T.; Lutz, M. L.; Sommerer, S. O. Effective core potential approaches to the chemistry of the heavier elements. Rev. Comp. Chem. 2007, 8, 145– 202.

25



(2) Barthelat, J. C.; Durand, P.; Serafini, A. Non-empirical pseudopotentials for molecular calculations. Mol. Phys. 1977, 33, 159–180. (3) Boys, S. F. Electronic wave functions. I. A general method of calculation for the stationary states of any molecular system. Proc. Roy. Soc. London A 1950, 200, 542–554. (4) McMurchie, L. E.; Davidson, E. R. Calculation of integrals over ab initio pseudopotentials. J. Comp. Phys. 1981, 44, 289–301. (5) Kolar, M. Pseudopotential matrix elements in the Gaussian basis. Comp. Phys. Commun. 1981, 23, 275–286. (6) Bode, B. M.; Gordon, M. S. Fast computation of analytical second derivatives with effective core potentials: Application to Si8C12, Ge8C12, and Sn8C12. J. Chem. Phys. 1999, 111, 8778–8784. (7) Pelissier, M.; Komiha, N.; Daudey, J. P. One-center expansion for pseudopotential matrix elements. J. Comput. Chem. 1988, 9, 298–302. (8) Skylaris, C.-K.; Gagliardi, L.; Handy, N. C.; Ioannou, A. G.; Spencer, S.; Willetts, A.; Simper, A. M. An efficient method for calculating effective core potential integrals which involve projection operators. Chem. Phys. Lett. 1998, 296, 445–451. (9) Flores-Moreno, R.; Alvarez-Mendez, R. J.; Vela, A.; Koster, A. M. Half-numerical evaluation of pseudopotential integrals. J. Comput. Chem. 2006, 27, 1009–1019. (10) Hu, A.; Chan, N. W.; Dunlap, B. I. Orbital angular momentum eigenfunctions for fast and numerically stable evaluations of closed-form pseudopotential matrix elements. J. Chem. Phys. 2017, 147, 074102. (11) Song, C.; Wang, L.-P.; Sachse, T.; Preis, J.; Presselt, M.; Martinez, T. J. Efficient implementation of effective core potential integrals and gradients on graphical processing units. J. Chem. Phys. 2015, 143, 14114. 26


Page 26 of 31

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


(12) Shaw, R. A.; Hill, J. G. Prescreening and efficiency in the evaluation of integrals over ab initio effective core potentials. J. Chem. Phys. 2017, 147, 074108. (13) Gill, P. M. W. Molecular integrals over gaussian basis functions. Adv. Quantum Chem. 1994, 25, 141–205. (14) Head-Gordon, M.; Pople, J. A. A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations. J. Chem. Phys. 1988, 89, 5777–5786. (15) Olver, F. W. J., Lozier, D. W., Boisvert, R. F., Clark, C. W., Eds. NIST Handbook of Mathematical Functions; Cambridge University Press: New York, 2010. (16) McMurchie, L. E.; Davidson, E. R. One- and two-electron integrals over Cartesian Gaussian functions. J. Comp. Phys. 1978, 26, 218–231. (17) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. Exact and approximate solutions to the McMurchie-Davidson tree-search problem. Int. J. Quantum Chem. 1991, 40, 809–827. (18) Ahlrichs, R. A simple algebraic derivation of the Obara-Saika scheme for general twoelectron interaction potentials. Phys. Chem. Chem. Phys. 2006, 8, 3072–3077. (19) Barca, G. M. J.; Loos, P.-F.; Gill, P. M. W. Many-electron integrals over Gaussian basis functions. I. Recurrence relations for three-electron integrals. J. Chem. Theory Comput. 2016, 12, 1735–1740. (20) Gill, P. M. W.; Johnson, B. G.; Pople, J. A. A simple yet powerful upper bound for Coulomb integrals. Chem. Phys. Lett. 1994, 217, 65–68. (21) Barca, G. M. J.; Loos, P.-F.; Gill, P. M. W. Many-electron integrals over Gaussian basis functions. II. Upper bounds for two- and three-electron integrals. in preparation. (22) Fuentealba, P.; Preuss, H.; Stoll, H.; Von Szentpaly, L. A proper account of corepolarization with pseudopotentials: single valence-electron alkali compounds. Chem. Phys. Lett. 1982, 89, 418–422. 27



(23) Fuentealba, P.; von Szentpaly, L.; Preuss, H.; Stoll, H. Pseudopotential calculations for alkaline-earth atoms. J. Phys. B 1985, 18, 1287. (24) Igel-Mann, G.; Stoll, H.; Preuss, H. Pseudopotentials for main group elements (IIIa through VIIa). Mol. Phys. 1988, 65, 1321–1328. (25) Kuchle, W.; Dolg, M.; Stoll, H.; Preuss, H. Ab initio pseudopotentials for Hg through Rn. Mol. Phys. 1991, 74, 1245–1263. (26) Bergner, A.; Dolg, M.; Kuchle, W.; Stoll, H.; Preuss, H. Ab initio energy-adjusted pseudopotentials for elements of groups 13-17. Mol. Phys. 1993, 80, 1431–1441. (27) Dolg, M.; Stoll, H.; Preuss, H.; Pitzer, R. M. Relativistic and correlation effects for element 105 (hahnium, Ha): a comparative study of M and MO (M = Nb, Ta, Ha) using energy-adjusted ab initio pseudopotentials. J. Phys. Chem. 1993, 97, 5852–5859. (28) Kaupp, M.; Schleyer, P. v. R.; Stoll, H.; Preuss, H. Pseudopotential approaches to Ca, Sr, and Ba hydrides. Why are some alkaline earth MX2 compounds bent? J. Chem. Phys. 1991, 94, 1360–1366. (29) Hay, P. J.; Wadt, W. R. Ab initio effective core potentials for molecular calculations. Potentials for K to Au including the outermost core orbitals. J. Chem. Phys. 1985, 82, 299–310. (30) Hay, P. J.; Wadt, W. R. Ab initio effective core potentials for molecular calculations. Potentials for the transition metal atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270–283. (31) Wadt, W. R.; Hay, P. J. Ab initio effective core potentials for molecular calculations. Potentials for main group elements Na to Bi. J. Chem. Phys. 1985, 82, 284–298. (32) Cundari, T. R.; Stevens, W. J. Effective core potential methods for the lanthanides. J. Chem. Phys. 1993, 98, 5555–5565.

28


Page 28 of 31

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


(33) Stevens, W. J.; Basch, H.; Krauss, M. Compact effective potentials and efficient sharedexponent basis sets for the first- and second-row atoms. J. Chem. Phys. 1984, 81, 6026–6033. (34) Stevens, W. J.; Krauss, M.; Basch, H.; Jasien, P. G. Relativistic compact effective potentials and efficient, shared-exponent basis sets for the third-, fourth-, and fifth-row atoms. Can. J. Chem. 1992, 70, 612–630. (35) Fernandez Pacios, L.; Christiansen, P. A. Ab initio relativistic effective potentials with spin-orbit operators. I. Li through Ar. J. Chem. Phys. 1985, 82, 2664–2671. (36) Hurley, M. M.; Pacios, L. F.; Christiansen, P. A.; Ross, R. B.; Ermler, W. C. Ab initio relativistic effective potentials with spin-orbit operators. II. K through Kr. J. Chem. Phys. 1986, 84, 6840–6853. (37) LaJohn, L. A.; Christiansen, P. A.; Ross, R. B.; Atashroo, T.; Ermler, W. C. Ab initio relativistic effective potentials with spin-orbit operators. III. Rb through Xe. J. Chem. Phys. 1987, 87, 2812–2824. (38) Ross, R. B.; Gayen, S.; Ermler, W. C. Ab initio relativistic effective potentials with spin-orbit operators. V. Ce through Lu. J. Chem. Phys. 1994, 100, 8145–8155. (39) Ermler, W. C.; Ross, R. B.; Christiansen, P. A. Ab initio relativistic effective potentials with spin-orbit operators. VI. Fr through Pu. Int. J. Quantum Chem. 1991, 40, 829– 846. (40) Nash, C. S.; Bursten, B. E.; Ermler, W. C. Ab initio relativistic effective potentials with spin-orbit operators. VII. Am through element 118. J. Chem. Phys. 1997, 106, 5133–5142. (41) Ross, R. B.; Powers, J. M.; Atashroo, T.; Ermler, W. C.; LaJohn, L. A.; Chris-

29



tiansen, P. A. Ab initio relativistic effective potentials with spin-orbit operators. IV. Cs through Rn. J. Chem. Phys. 1990, 93, 6654–6670. (42) Cohen, M. H.; Heine, V. Cancellation of kinetic and potential energy in atoms, molecules, and solids. Phys. Rev. 1961, 122, 1821–1826. (43) Dolg, M.; Cao, X. Relativistic pseudopotentials: their development and scope of applications. Chem. Rev. 2012, 112, 403–480. (44) Shao, Y.; Zhengting, G.; Epifanovsky,E.; Gilbert, A.T.B.; Wormit, M.; Kussmann, J.; Lange, A.W.; Behn, A.; Deng, J.; Feng, X. et al. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 2015, 113, 184– 215. (45) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. et al. General atomic and molecular electronic structure system. J. Comput. Chem. 1993, 14, 1347–1363. (46) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P. et al. The Dalton quantum chemistry program system. WIREs Comput. Mol. Sci. 2014, 4, 269–284.

30


Page 30 of 31

Page 31 of 31 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 [d|0]d (0,0,0)

[0|0]d (0,0,0)

[0|0]d (1,0,0)

[p|0]d (0,0,0)

[p|0]s(1,0,1)

[p|0]d (1,0,0)

[0|0]s(1,0,1)

[0|0]p (0,0,1)

[0|0]s(2,0,1)

31

[p|0]p (0,0,1)

[0|0]d (2,0,0)

[0|0]p (1,0,1)


[0|0]s(0,0,2)

An Efficient Method for Calculating Effective Core Potential Integrals

Recommend Documents