Two-Center Noninteger-n Overlap, Coulomb, and Kinetic Energy

Subscriber access provided by TULANE UNIVERSITY

Article

Two-Center Noninteger-n Overlap, Coulomb, and KineticEnergy Integrals by Numerical Contour Integration Harris J. Silverstone J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/jp5070159 • Publication Date (Web): 20 Aug 2014 Downloaded from http://pubs.acs.org on September 5, 2014

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 free 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 accessible to all readers and 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 8

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

Two-Center Noninteger-n Overlap, Coulomb, and Kinetic-Energy Integrals by Numerical Contour Integration Harris J. Silverstone⇤ Department of Chemistry, Johns Hopkins University University, Baltimore, MD 21218 E-mail: [email protected]

Abstract

is an integer or not.

Numerical contour integration accurately evaluates two-center overlap, Coulomb, and kineticenergy integrals involving noninteger-n Slatertype orbitals. Keywords: Slater-type orbitals; exponentialtype orbitals; molecular integrals; molecular electronic structure

Notation and formulas A Slater-type orbital (STO) is characterized by parameters n, l, m, and z , and coordinates r: ynlmz (r) ⌘ rn 1 e

Ylm (q , f ).

(1)

Ylm (q , f ) denotes a spherical harmonic with standard Condon-Shortley phase. Below, ya (r) is used as a convenient shorthand for yna la ma za (r). Designate the two-center overlap, Coulomb, and kinetic energy integrals by Sab (R), Cab (R), and Tab (R), defined here by

Introduction Recent renewed interest in noninteger-n Slatertype orbitals, which are asymptotically appropriate for nonrelativistic electronic wave functions, p 2 a 2Z2 ± k and which can match the r behavior of 1–3 relativistic wave functions near r = 0, raises the need for accurate evaluation of molecular integrals over wide ranges of parameters and distances. In a new approach, free from special cases, Ba˘gcı and Hoggan 4 calculated arbitrary noninteger-n overlap integrals to 36-digit accuracy using twodimensional numerical integration. Their idea of combining maximum simplicity and numerical integration is adapted to Fourier-transform-derived formulas, with the numerical integration being along a one-dimensional contour in the complex plane. The same 36-digit accuracy is achieved in significantly shorter computation times. Indeed, numerical contour integration may be a practical technique for other molecular integrals, whether n ⇤ To

zr

Sab (R) ⌘ Cab (R) ⌘ Tab (R) ⌘

Z Z Z

dV ya (r)⇤ yb (r dV1

Z

R)

dV2 r121 ya (r1 )⇤ yb (r2

dV ya (r)⇤

✓

◆

1 2 — yb (r 2

(2) R) (3) R). (4)

The Fourier-transform convolution theorem quickly reduces the multidimensional integrations to one dimension, simplifies the dependence on R, and reveals the essential sameness of these three seemingly different integrals. 5–8 A sketch of the Fourier-transform method follows: 9 Start with the

whom correspondence should be addressed

ACS Paragon Plus Environment

1


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

Fourier transform Fnlmz (k) of an STO, Fnlmz (k) ⌘

Z

dVeik·r ynlmz (r)

fnlz (k) = 4pil

Z • 0

rn+1 e

zr

This finishes the sketch of the basic idea of the Fourier-transform method. The details can be found in Refs. 5–8. When both na and nb are integers, the radial factor (10) can be evaluated analytically, as, for instance, in Ref. 8. When na and nb are not integers, no simple analytical expressions are known, and, in Refs. 6 and 7, series expansions requiring several special cases were developed. The novelty in the present paper is that the radial factors are evaluated by numerical contour integration, and, whether the n’s are integer or not, whether the z ’s are equal or not, there are no special cases.

(5)

= fnlz (k)Ylm (qk , fk )

(6) jl (kr) dr.

(7)

The jl (kr) is a spherical Bessel function, and an explicit formula for fnlz (k) is given in Eq. (21) below. Apply the Fourier-transform inversion and convolution theorems (note the shorthand a for na , la , ma , za , etc.) to Eq. (2). The overlap integral Sab (R) becomes an integral over the Fourier transforms in which the angular integration is straightforward: 3

Sab (R) = (2p) la +lb

=

Z

Â

d 3 k Fa (k)⇤ Fb (k)eik·R (0)

l =|la lb |

dabl (qR , fR )Jabl (R)

Angular factors

(8)

The angular factors arise from the integral of spherical-harmonic triple products and are conveniently expressed in terms of 3 j-symbols (instead of the Condon-Shortley coefficients used previously 8 ), and which are built into, for instance, Mathematica:

(9)

where dabl (qR , fR ) is the result of the angular integration [Eq. (16) below], and where the factor (0) Jabl (R) is defined by the N = 0 case of (N) Jabl (R) =

1 l i 2p 2

Z • 0

Z

k2+2N fa (k) fb (k) jl (kR) dk (10)

[cf. Eq. (17) of Ref. 8]. Evaluating the overlap integral has been transformed to evaluating the one(0) dimensional integral for Jabl (R). Similar manipulations for the Coulomb and kinetic-energy integrals lead to the following unifying formulas, with the distinguishing detail the power of k in Eq. (10): (0) Sab (R) = Iab (R)

(11)

Cab (R) = 4pIab (R)

( 1)

(12)

1 (1) Tab (R) = Iab (R). 2

(13)

(N) Iab (R) =

la +lb

Â

l =|la lb |

(N)

dabl (qR , fR )Jabl (R).

Page 2 of 8

1 2 Ylm (q , f )Ylm (q , f )Yl3 m1 m2 (q , f ) sin q dq df 1 2 p = (2l1 + 1)(2l2 + 1)(2l3 + 1)/4p ✓ ◆✓ ◆ l1 l2 l3 l1 l2 l3 ⇥ 0 0 0 m1 m2 m1 m2 (15)

dabl (qR fR )

p =( 1)la +mb (2l + 1)(2la + 1)(2lb + 1)/4p ✓ ◆✓ ◆ la lb l la lb l ⇥ 0 0 0 ma mb mb ma ⇥Ylmb

ma

(qR , fR ).

(16)

Radial factors (N)

The radial factors Jabl (R) arise from the radial integration in the Fourier-transform convolution theorem. The main point of this paper is to replace Eq. (10) by Eq. (23) below, and then to integrate numerically. Eq. (10) is not ideal for numerical integration because, as k ! •, the spherical Bessel function oscillates trigonometrically, and the integrand goes

(14)


2

Page 3 of 8

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


to 0 only like k to an inverse power, that is, relatively slowly. Two relevant facts — (1) the spherical Bessel function jl (kR) is a sum of eikR and e ikR terms, and (2) the integrand of Eq. (10) is an even function of k (note that la + lb + l must be even) — and a substitution of convenience, k = ix, permit putting Eq. (10) into a more useful initial form: (N) Jabl (R) =

To take advantage of the property that the integrand of Eq. (17) goes to 0 exponentially fast in the right half-plane, deform the integration path from the imaginary axis to a contour L, shown in Figure 1, that starts from • in the first quadrant, passes between 0 and the smaller of za and zb , where the integrand has branch points (or poles), then returns to • in the fourth quadrant. It is an essential point of complex variable theory that such a deformation does not change the value of the integral, but it does reduce dramatically the magnitude of the integrand along the integration path and in this way facilitates numerical integration. When N = 1 and l = la + lb , there is a simple pole at the origin whose contribution to the integral is pi times its residue. The residue plus contour contributions are

( 1)N 4p 2

⇥i

Z

i•

x2+2N fa (ix) fb (ix)Kl (xR) dx

i•

(17)

Kl (z) is a modified spherical Bessel function of the second kind. r 2 Kl (z) = K (z) (18) pz l +1/2 ✓ ◆l z e l 1 d = ( z) (19) z dz z l

=

Â zµ

l 1

µ=0

(2l 2l µ (l

µ)! e µ)!µ!

z

(N)

J abl (R) = dN, 1 dl ,la +lb 4p( 1)la +lb za

⇥ G(na + la + 2)G(nb + lb + 2) (2la + 2lb 1)!! ⇥ (2la + 1)!!(2lb + 1)!!

(20)

⇥ (za R)

(While Eq. (17) does not appear explicitly in Ref. 8, it tracks the discussion between Eqs. (6) and (8) in Ref. 6.) The integrand of Eq. (17) contains the factor e xR [cf. Eq. (20)], which goes to 0 strongly in the right half-plane. An explicit formula for fnlz (ix) (Eq. (13) of Ref. 5), fnlz (ix) = 4p

l

Â

G(n

µ=0

⇣ ⇥ (z + x)µ

n 1

µ + 1) (l + µ)! x 2 2µ (l µ)!µ! ( 1)l

µ

(z

x)µ

( 1)N + i 4p 2

0

4pG(n + l + 2) z (2l + 1)!!

n 1

⌘

Z

L

(zb R)

nb 3/2

lb 1/2

x2+2N fa (ix) fb (ix)Kl (xR) dx (23)

a

b

, Figure 1: Schematic path L for numerical contour integration.

(21)

n l 2

la 1/2

zb

µ 1

reveals that the only singularities of the integrand of Eq. (17) in the right half-plane are poles or branch points on the real axis at x = za and x = zb , depending on whether the n’s are integer or noninteger, and possibly at x = 0 in the Coulomb case. Note that fnlz (ix) ⇠ xl for small x, fnlz (ix) ⇠

na 3/2

Numerical contour integration Eq. (23) is convenient for numerical integration because the integrand vanishes exponentially at the end points of L. For a given integral, the “optimum” path is likely a steepest descent path, along which the magnitude of the integrand decreases most rapidly. The steepest path would be

( x)l . (22)


3


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

0

a

Page 4 of 8

trol the number of significant digits, as is discussed in some detail by Ba˘gcı and Hoggan. 4 Further, not only does Mathematica make the numerical contour integration easy, but it provides some zeroeffort parallelization. E.g., use of (schematically)

b

´

(N)

Figure 2: “Good enough” path contour integration.

L0

ParallelTable[Jabl (R), {l , |la

for numerical

lb |, la + lb , 2}] (25)

instead of (N)

Table[Jabl (R), {l , |la

specific for each integral and time consuming to determine. Fortunately, much simpler, more general paths seem sufficiently accurate (defined here as 36 significant digits) and efficient. In particular, consider a rectangular path L0 as indicated in Figure 2, where the main caveat is to avoid pathologically large values of the integrand near x = za and zb . Such rectangular L0 reproduced to 36-digit accuracy all the integrals reported in Ref. 4 (except for one typo and one probable transcription error).

lb |, la + lb , 2}]

(26)

to parallelize the calculation of the “table” of the (N) Jabl (R) reduced the overall computation time by a factor of 2, and the computation time for higher-l cases by 2.5, on a quad core i7 MacBook Pro. Numerical contour integration appears to work for any (or, at least, reasonable) values of na and nb , integer or noninteger, and any values of za and zb , unequal or equal. There are no special cases. To assemble the molecular integrals, the radial (N) Jabl (R) are multiplied by the angular dabl (qR fR ) computed from Eq. (16), and the products added as in Eq. (14).

Mathematica program Before computational tools like Mathematica, numerical contour integration would have required considerable effort. With current versions of Mathematica, the numerical-integration command NIntegrate can do contour as well as ordinary integration, and it is sufficient just to specify the integrand and contour. Mathematica does the rest, including path truncation, quadrature method, subdivision, and accuracy. Numerical explorations and proof of concept are easily accomplished. For the noninteger-n STO radial factors, there are two definitions to program: fnlz (ix) and Kl (xR), Eqs. (20) and (21). Armed with these definitions, the following instruction in Mathematica carries out the numerical integration along the contour L0 of Figure 2, where x0 ± iy0 are the two corners in L0 ,

Illustrative calculations The recently reported 36-digit-accuracy overlap integrals in Tables I–III of Ba˘gcı and Hoggan, 4 encompassing variety and high accuracy, seem an ideal benchmark. The ynlmz (r) of Eq. (1) needs the normalization factor (2z )n+1/2 Nnz = p , G(2n + 1)

(27)

and, since Ba˘gcı and Hoggan use real Ylm , that accommodation has been made for Tables 1 and 2 here. The basic Mathematica program described above — coded to minimize programming effort, not computational efficiency — calculated the 12 integer-n and 9 noninteger-n overlap integrals to better than 37-digit-accuracy in 19 seconds for the integer-n and 107 seconds for the noninteger-n sets on a 2.3 GHz quad core i7 MacBook Pro running Mathematica version 10. The times quoted are elapsed “SessionTime” — i.e., wall-clock, not CPU times. For integer-n STO’s, the results from numerical

i NIntegrate[ x2 fa (ix) fb (ix)Kl (xR) 4p 2 , {x, •, x0 + iy0 , x0 iy0 , •}] (24) It helps further to specify the NIntegrate options WorkingPrecision, Method, MaxRecursion, and AccuracyGoal, to con-


4

Page 5 of 8

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


Notes and References

contour integration are reported in Table 1 to 37 significant digits. Besides agreeing with Ba˘gcı and Hoggan, 4 the integer-n integrals have been independently checked using analytic formulas. 8 The noninteger-n results are reported in Table 2, and are essentially in agreement with Ba˘gcı and Hoggan. There were two exceptions: there is a typo in the first entry of Table III, 4 which should be “E-01” rather than “E-02”. The other is the extremely small seventh entry of Table I, 4 in which there is disagreement in the 13th digit. Note that this integral is 7.09⇥10 90 , so that its value is of only academic importance; its exact calculation contains a severe numerical cancellation. 10

(1) Drake, G. W. F., Ed. Springer Handbook of Atomic, Molecular, and Optical Physics; Springer: New York, 2006. (2) Grant, I. P. Relativistic Quantum Theory of Atoms and Molecules; Springer: New York, 2007. (3) Ba˘gcı, A.; Hoggan, P. E. Performance of Numerical Approximation on the Calculation of Two-Center Two-Electron Integrals over Non-Integer Slater-Type Orbitals Using Elliptical Coordinates. arXiv:1405.5436 [quant-ph] 2014, 1–33. (4) Ba˘gcı, A.; Hoggan, P. E. Performance of Numerical Approximation on the Calculation of Overlap Integrals with Noninteger Slater-Type Orbitals. Phys. Rev. E 2014, 89, 053307.

Discussion When attacked by the Fourier-transform convolution theorem, STO multicenter integrals typically involve an inverse radial Fourier transform, which is the difficult step. The analytical resolution of two-center noninteger-n integrals can lead to intimidating series expansions or to what are essentially single-center expansions. Ba˘gcı and Hoggan 4 have shown that the single-center type expansions converge slowly compared to numerical integration, and that in practical terms numerical integration is able to attain orders of magnitude higher accuracy, and with integrands that are simple in structure. In this paper, Ba˘gcı and Hoggan’s ideas have been applied to the inverse radial Fourier transform, and the resulting numerical integration is one-dimensional along a contour in the complex plane. Because the integral is onedimensional, the computer time is orders of magnitude shorter than for two-dimensional numerical integration. Because tools like Mathematica have numerical contour integration built-in, programming effort, at least to get benchmark values, is minimal. Programming effort to achieve high efficiency is not expected to be minimal. Beyond two-center noninteger-n STO integrals, the possibility of using numerical contour integration in more complicated STO multicenter integrals is tantalizing.

(5) Silverstone, H. J. On the Evaluation of TwoCenter Overlap and Coulomb Integrals with Noninteger-n Slater-Type Orbitals. J. Chem. Phys. 1966, 45, 4337–4341. (6) Silverstone, H. J. Series Expansion for TwoCenter Noninteger-n Overlap Integrals. J. Chem. Phys. 1967, 46, 4368–4376. (7) Silverstone, H. J. Series Expansion for TwoCenter Noninteger-n Coulomb Integrals. J. Chem. Phys. 1967, 46, 4377–4380. (8) Todd, H. D.; Kay, K. G.; Silverstone, H. J. Unified Treatment of Two-Center Overlap, Coulomb, and Kinetic-Energy Integrals. J. Chem. Phys. 1970, 53, 3951–3956. (9) For more details, see the material from Eqs. (1)–(17) in Ref. 8. (10) The exact value of this integral, obtained from a Mathematica implementation of Ref. 8, has the form, C(A p p B e) 111111e 1001/4 /D, where A, B,C, and D are integers of digit-lengths 53, 53, 79, and 71, respectively, and for which D has no common factor with A, B, or C, and A and B are relatively prime. That there is a cancellation of 43 significant digits in evaluating


5


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

p (A B e) augurs a challenge to accurate numerical calculation.


6

Page 6 of 8

Page 7 of 8

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: Integer-n overlap integrals from numerical contour integration na

la ma za R

nb

lb mb zb R

qR

fR

overlap integral (with real Ylm as in Ref. 4)

0

0

1

0

0

9.98337 28456 63420 62391 22146 22559 01738 7 ⇥ 10 4.42287 76698 82608 80679 54150 24354 52085 1 ⇥ 10

4

8.37290 47190 38628 26960 63231 72822 20235 3 ⇥ 10

40

7.09059 49024 34538 18180 82255 51808 39457 3 ⇥ 10

90

7.23339 14676 87606 82677 55821 58140 80371 9 ⇥ 10

3

1.09293 91789 07866 13706 54526 97582 19446 7 ⇥ 10

1

1

0

0

1 10

1

0

0

8

0

0

5

8

0

0

1 10 1 10

3

2

1

40

3

2

1

10

0

0

6

5

4 190

5

4

3

10

6

5

4

545 4 101 4 1001 4

5

4

3

455 4 99 4 999 4

p 6 p 6

2p 3 2p 3

0

0

0

0

13 12 12 13 12 12 21 10

6

45

25 11

6

40

4

40 5

13 12 12 13 12 12 9

8

6

45

0

0

9009 200

12

8

3

8991 200

p 4

3p 4

3

24

12

4

3

6

0

0

4

3

24

12

4

3

6

25p 36

7p 36

4

4

5

5

4

4

5

0

0

7.15537 44350 12757 24983 22690 59617 20779 1 ⇥ 10

10

3.22601 93043 96471 01300 54876 30236 40989 7 ⇥ 10

9

1.35310 57870 24712 38186 18677 70288 53421 6 ⇥ 10

4

5.38980 68533 81437 73017 27203 24019 14298 9 ⇥ 10

5

9.48379 22083 22556 78538 44190 07653 62745 1 ⇥ 10

2

3.56825 98684 57483 60689 04872 31073 68139 2 ⇥ 10

1

Table 2: Noninteger-n overlap integrals from numerical contour integration na

la ma za R

nb

lb mb zb R qR fR overlap integral (with real Ylm as in Ref. 4)

73 10 32 5 19 5 57 10 18 5 63 10 41 10 103 10 51 10

4

4

3

4

4

1

0

0

1

0

0

0

0

0

0

0

0

0

1

1

1

1

0

0

2

1

2

1

1

4

4

4

4p 3

2

2

2

2

9

0

0

0

0

9

0

3

0

0

4

4

5

11 2 37 10 103 10 51 10

2p 5 2p 3

p

5

15 2 77 50 91 50 14 25 27 2

0

0

27 10 77 25 147 50 26 25 33 2 23 2 15 2

73 10 32 5 11 2 19 5

4

4

5

0

0

1.01734 31495 95668 84009 52107 36427 25566 0 ⇥ 10

1

2.90802 04650 66341 47700 88166 91317 05703 2 ⇥ 10

1

6.49621 73637 32212 98813 20684 64320 73379 8 ⇥ 10

2

2.93541 97236 64768 15362 90672 35811 14085 2 ⇥ 10

2

3.68837 33855 08336 58641 31918 22868 35839 8 ⇥ 10

1

3.12099 12216 53204 22891 71991 67638 91165 6 ⇥ 10

1

8.66889 50632 72588 09962 28315 86072 66015 8 ⇥ 10

1

1.85058 95468 20783 59133 54756 23091 64154 6 ⇥ 10

2

1.52896 89539 18532 09475 62114 33969 06544 6 ⇥ 10

5


7


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

0

zb

za L

Schematic path for numerical contour integration


8

Page 8 of 8

Two-Center Noninteger-n Overlap, Coulomb, and Kinetic Energy

Recommend Documents