Molecular Integral Evaluation Trygve Helgaker Centre for - - PowerPoint PPT Presentation

Molecular Integral Evaluation

Trygve Helgaker

Centre for Theoretical and Computational Chemistry Department of Chemistry, University of Oslo, Norway

The 11th Sostrup Summer School Quantum Chemistry and Molecular Properties July 4–16 2010

Trygve Helgaker (CTCC, University of Oslo) Molecular Integral Evaluation 11th Sostrup Summer School (2010) 1 / 34

Molecular integral evaluation

◮ one-electron interactions

◮ overlap, multipole-moment, and kinetic-energy integrals ◮ Coulomb attraction integrals

Oab =

• χa(r) ˆ

O(r)χb(r) dr

◮ two-electron interactions

◮ Coulomb integrals

gabcd = χa(r1)χb(r1)χc(r2)χd(r2) r12 dr1 dr2

◮ Coulomb and exchange contributions to Fock/KS matrix

Fab =

• cd
• 2gabcdDcd − gacbdDcd
• ◮ basis functions

◮ primitive Cartesian GTOs

◮ integration schemes

◮ McMurchie–Davidson, Obara–Saika and Rys schemes Trygve Helgaker (CTCC, University of Oslo) Introduction 11th Sostrup Summer School (2010) 2 / 34

Overview

1

Cartesian and Hermite Gaussians

◮ properties of Cartesian Gaussians ◮ properties of Hermite Gaussians 2

Simple one-electron integrals

◮ Gaussian product rule and overlap distributions ◮ overlap distributions expanded in Hermite Gaussians ◮ overlap and kinetic-energy integrals 3

Coulomb integrals

◮ Gaussian electrostatics and the Boys function ◮ one- and two-electron Coulomb integrals over Cartesian Gaussians 4

Sparsity and screening

◮ Cauchy–Schwarz screening 5

Coulomb potential

◮ early density matrix contraction ◮ density fitting Trygve Helgaker (CTCC, University of Oslo) Introduction 11th Sostrup Summer School (2010) 3 / 34

Cartesian Gaussian-type orbitals (GTOs)

◮ We shall consider integration over primitive Cartesian Gaussians centered at A:

Gijk(r, a, A) = xi

Ayj Azk A exp(−ar2 A),

     a > 0

• rbital exponent

rA = r − A

electronic coordinates

i ≥ 0, j ≥ 0, k ≥ 0

quantum numbers

◮ total angular-momentum quantum number ℓ = i + j + k ≥ 0 ◮ Gaussians of a given ℓ constitutes a shell:

s shell: G000, p shell: G100, G010, G001, d shell: G200, G110, G101, G020, G011, G002

◮ Each Gaussian factorizes in the Cartesian directions:

Gijk(a, rA) = Gi(a, xA)Gj(a, yA)Gk(a, zA), Gi(a, xA) = xi

A exp(−ax2 A)

◮ note: this is not true for spherical-harmonic Gaussians nor for Slater-type orbitals

◮ Gaussians satisfy a simple recurrence relation:

xAGi(a, xA) = Gi+1(a, xA)

◮ The differentiation of a Gaussian yields a linear combination of two Gaussians

∂Gi(a, xA) ∂Ax = − ∂Gi(a, xA) ∂x = 2aGi+1(a, xA) − iGi−1(a, xA)

Trygve Helgaker (CTCC, University of Oslo) Cartesian and Hermite Gaussians 11th Sostrup Summer School (2010) 4 / 34

Hermite Gaussians

◮ The Hermite Gaussians are defined as

Λtuv(r, p, P) = ∂ ∂Px t ∂ ∂Py u ∂ ∂Pz v exp

• −pr2

P

• ,

rP = r − P

◮ Like Cartesian Gaussians, they also factorize in the Cartesian directions:

Λt(xP) = (∂/∂Px)t exp

• −px2

P

• ← a Gaussian times a polynomial of degree t

◮ Hermite Gaussians yields the same spherical-harmonic functions as do Cartesian Gaussians

◮ they may therefore be used as basis functions in place of Cartesian functions

◮ We shall only consider their use as intermediates in the evaluation of Gaussian integrals Cartesian product → Gi(xA)Gj(xB) =

i+j

• t=0

E ij

t Λt(xP) ← Hermite expansion

◮ such expansions are useful because of simple integration properties such as

∞

−∞

Λt(x) dx = δt0 π p

◮ McMurchie and Davidson (1978)

◮ We shall make extensive use of the McMurchie–Davidson scheme for integration

Trygve Helgaker (CTCC, University of Oslo) Cartesian and Hermite Gaussians 11th Sostrup Summer School (2010) 5 / 34

Integration over Hermite Gaussians

◮ From the definition of Hermite Gaussians, we have

∞

−∞

Λt(x) dx = ∞

−∞

∂ ∂Px t exp

• −px2

P

• dx

◮ We now change the order of differentiation and integration by Leibniz’ rule:

∞

−∞

Λt(x) dx = ∂ ∂Px t ∞

−∞

exp

• −px2

P

• dx

◮ The basic Gaussian integral is given by

∞

−∞

exp

• −px2

P

• dx =

π p

◮ Since the integral is independent of P, differentiation with respect to P gives zero:

∞

−∞

Λt(x) dx = δt0 π p

◮ only integrals over Hermite s functions do not vanish Trygve Helgaker (CTCC, University of Oslo) Cartesian and Hermite Gaussians 11th Sostrup Summer School (2010) 6 / 34

Hermite recurrence relation

◮ Cartesian Gaussians satisfy the simple recurrence relation

xAGi = Gi+1

◮ The corresponding Hermite recurrence relation is slightly more complicated:

xPΛt = 1 2p Λt+1 + tΛt−1

◮ The proof is simple:

◮ from the definition of Hermite Gaussians, we have

Λt+1 = ∂ ∂Px t ∂ ∂Px exp(−px2

P) = 2p

∂ ∂Px t xPΛ0

◮ inserting the identity

∂ ∂Px t xP = xP ∂ ∂Px t − t ∂ ∂Px t−1 we then obtain Λt+1 = 2p

• xP

∂ ∂Px t − t ∂ ∂Px t−1 Λ0 = 2p (xPΛt − tΛt−1)

Trygve Helgaker (CTCC, University of Oslo) Cartesian and Hermite Gaussians 11th Sostrup Summer School (2010) 7 / 34

Comparison of Cartesian and Hermite Gaussians

Cartesian Gaussians Hermite Gaussians definition Gi = xi

A exp

• −ax2

A

• Λt =

∂t ∂Pt

x exp

• −px2

P

• recurrence

xAGi = Gi+1 xPΛt =

1 2p Λt+1 + tΛt−1

differentiation

∂Gi ∂Ax = 2aGi+1 − iGi−1 ∂Λt ∂Px = Λt+1

3 2 1 1 2 3 2 1 1 2 3 2 1 1 2 3 2 1 1 2

Trygve Helgaker (CTCC, University of Oslo) Cartesian and Hermite Gaussians 11th Sostrup Summer School (2010) 8 / 34

The Gaussian product rule

◮ The most important property of Gaussians is the Gaussian product rule ◮ The product of two Gaussians is another Gaussian centered somewhere on the line

connecting the original Gaussians; its exponent is the sum of the original exponents: exp

• −ax2

A

• exp
• −bx2

B

• =

exp

• −µX 2

AB

• exponential prefactor

exp

• −px2

P

• product Gaussian

where Px = aAx + bBx p ← “center of mass” p = a + b ← total exponent XAB = Ax − Bx, ← relative separation µ = ab a + b ← reduced exponent

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1

◮ The Gaussian product rule greatly simplifies integral evaluation

◮ two-center integrals are reduced to one-center integrals ◮ four-center integrals are reduced to two-center integrals Trygve Helgaker (CTCC, University of Oslo) Gaussian product rule and Hermite expansions 11th Sostrup Summer School (2010) 9 / 34

Overlap distributions

◮ The product of two Cartesian Gaussians is known as an overlap distribution:

Ωij(x) = Gi(x, a, Ax)Gj(x, b, Bx)

◮ The Gaussian product rule reduces two-center integrals to one-center integrals

• Ωij(x) dx = KAB
• xi

Axj B exp

• −px2

P

• dx

◮ the Cartesian monomials still make the integration awkward ◮ we would like to utilize the simple integration properties of Hermite Gaussians

◮ We therefore expand Cartesian overlap distributions in Hermite Gaussians

• verlap distribution → Ωij(x) = i+j

t=0E ij t Λt(xP) ← Hermite Gaussians

◮ note: Ωij(x) is a single Gaussian times a polynomial in x of degree i + j ◮ it may be exactly represented as a linear combination of Λt with 0 ≤ t ≤ i + j

◮ The expansion coefficients may be evaluated recursively:

E i+1,j

t

=

1 2p E ij t−1 + XPAE ij t + (t + 1)E ij t+1

◮ we shall now prove these recurrence relations Trygve Helgaker (CTCC, University of Oslo) Gaussian product rule and Hermite expansions 11th Sostrup Summer School (2010) 10 / 34

Recurrence relations for Hermite expansion coefficients

◮ A straighforward Hermite expansion of Ωi+1,j

Ωi+1,j = KABxi+1

A

xj

B exp

• −pr2

P

• =
• t

E i+1,j

t

Λt

◮ An alternative Hermite expansion of Ωi+1,j

Ωi+1,j = xAΩij = (x − Ax)Ωij = (x − Px)Ωij + (Px − Ax)Ωij = xPΩij + XPAΩij =

• t

E ij

t xPΛt + XPA

• t

E ij

t Λt

=

• t

E ij

t

• 1

2p Λt+1 + tΛt−1 + XPAΛt

• =
• t
• 1

2p E ij t−1 + (t + 1)E ij t+1 + XPAE ij t

• Λt

◮ we have here used the recurrence xPΛt =

1 2p Λt+1 + tΛt−1

◮ A comparison of the two expansions yields the recurrence relations

E i+1,j

t

=

1 2p E ij t−1 + XPAE ij t + (t + 1)E ij t+1

Trygve Helgaker (CTCC, University of Oslo) Gaussian product rule and Hermite expansions 11th Sostrup Summer School (2010) 11 / 34

Overlap integrals

◮ We now have everything we need to evaluate overlap integrals:

Sab = Ga |Gb Ga = Gikm (r, a, A) = Gi(xA)Gk(yA)Gm(zA) Gb = Gjln (r, b, B) = Gj(xB)Gl(yB)Gn(zB)

◮ The overlap integral factorizes in the Cartesian directions:

Sab = SijSklSmn, Sij =

• Gi(xA)
• Gj(xB)
• ◮ Use the Gaussian product rule and Hermite expansion:

Sij =

• Ωij(x) dx =

i+j

• t=0

E ij

t

• Λt(xP) dx =

i+j

• t=0

E ij

t δt0

• π

p = E ij

• π

p

◮ only one term survives!

◮ The total overlap integral is therefore given by

Sab = E ij

0 E kl 0 E mn

π p 3/2

Trygve Helgaker (CTCC, University of Oslo) Simple one-electron integrals 11th Sostrup Summer School (2010) 12 / 34

Dipole-moment integrals

◮ Many other integrals may be evaluated in the same manner ◮ For example, dipole-moment integrals are obtained as

Dij =

• Gi(xA) |xC | Gj(xB)
• =
• Ωij(xP)xC dx

=

i+j

• t=0

E ij

t

• Λt(xP)xcdx

(expand in Hermite Gaussians) =

i+j

• t=0

E ij

t

• [xPΛt(xP) + XPC Λt(xP)] dx

(use xC = xP + XPC ) =

i+j

• t=0

E ij

t 1 2p Λt+1(xP) + XPC Λt(xP) + tΛt−1(xP)

• dx

(use xPΛt =

1 2p Λt+1 + tΛt−1)

= E ij

0 XPC

π p + E ij

1

π p (s functions integration)

◮ The final dipole-moment integral is given by

Dij =

• E ij

1 + XPC E ij

π p

Trygve Helgaker (CTCC, University of Oslo) Simple one-electron integrals 11th Sostrup Summer School (2010) 13 / 34

Kinetic-energy integrals

◮ As a final example of non-Coulomb integrals, we consider the kinetic-energy integrals:

Tab = − 1

2

• Ga
• ∂2

∂x2 + ∂2 ∂y2 + ∂2 ∂z2

• Gb
• Tab = TijSklSmn + SijTklSmn + SijSklTmn

Sij =

• Gi(xA)
• Gj(xB)
• Tij = − 1

2

• Gi(xA)
• ∂2

∂x2

• Gj(xB)
• ◮ Differentiation of Cartesian Gaussians gives

d dx Gj(xB) = −2bGj+1 + jGj−1 d2 dx2 Gj(xB) = 4b2Gj+2 − 2b(2j + 1)Gj + j(j − 1)Gj−2

◮ the first derivative is linear combination of two undifferentiated Gaussians ◮ the second derivative is linear combination of three undifferentiated Gaussians

◮ Kinetic-energy integrals are linear combinations of up to three overlap integrals:

Tij = −2b2Si,j+2 + b(2j + 1)Si,j − 1

2 j(j − 1)Si,j−2

Trygve Helgaker (CTCC, University of Oslo) Simple one-electron integrals 11th Sostrup Summer School (2010) 14 / 34

Coulomb integrals

◮ We shall consider two types of Coulomb integrals:

◮ one-electron nuclear-attraction integrals

• Ga(rA)
• r−1

C

• Gb(rB)
• ◮ two-electron repulsion integrals
• Ga(r1A)Gb(r1B)
• r−1

12

• Gc(r2C )Gd(r2D)
• ◮ Coulomb integrals are not separable in the Cartesian directions

◮ however, they may be reduced to a one-dimensional integral on a finite interval:

Fn(x) = 1 exp

• −xt2

t2ndt ← the Boys function

◮ this makes the evaluation of Gaussian Coulomb integrals relatively simple

◮ We begin by evaluating the one-electron spherical Coulomb integral:

Vp = exp

• −pr2

P

• rC

dr = 2π p F0

• pR2

PC

• ◮ we shall next go on to consider more general Coulomb integrals

Trygve Helgaker (CTCC, University of Oslo) Coulomb integrals 11th Sostrup Summer School (2010) 15 / 34

Coulomb integral over a spherical Gaussian I

◮ We would like to evaluate the three-dimensional Coulomb-potential integral

Vp = exp

• −pr2

P

• rC

dr

1

The presence of r−1

C

is awkward and is avoided by the substitution 1 rC = 1 √π ∞

−∞

exp

• −r2

C t2

dt ← Laplace transform to yield the four-dimensional integral Vp = 1 √π

• exp
• −pr2

P

∞

−∞

exp

• −t2r2

C

• dt dr

2

To prepare for integration over r, we invoke the Gaussian product rule exp

• −pr2

P

• exp
• −t2r2

C

• = exp
• −

pt2 p + t2 R2

CP

• exp
• −
• p + t2

r2

S

• to obtain (where the exact value of S in rS does not matter):

Vp = 1 √π ∞

−∞

• exp
• −(p + t2)r2

S

• dr exp
• −

pt2 p + t2 R2

CP

• dt

– to be continued. . .

Trygve Helgaker (CTCC, University of Oslo) Spherical Coulomb integrals 11th Sostrup Summer School (2010) 16 / 34

Coulomb integral over a spherical Gaussian II

– carried over from previous slide: Vp = 1 √π ∞

−∞

−(p + t2)r2

S

• dr exp
• −

pt2 p + t2 R2

CP

• dt

3

Integration over all space r now yields a one-dimensional integral Vp = 2 √π ∞

• π

p + t2 3/2 exp

• −pR2

CP

t2 p + t2

• dt

4

To introduce a finite integration range, we perform the substitution u2 = t2 p + t2 ⇒ 1 + pt−2 = u−2 ⇒ pt−3 dt = u−3du ⇒ dt = p−1 t2 u2 3/2 du and obtain Vp = 2π p 1 exp

• −pR2

CPu2

du = 2π p F0

• pR2

PC

• ◮ We have here introduced the Boys function

Fn(x) = 1 exp

• −xt2

t2ndt – a 3D integral over all space has been reduced to a 1D integral over [0, 1]

Trygve Helgaker (CTCC, University of Oslo) Spherical Coulomb integrals 11th Sostrup Summer School (2010) 17 / 34

Gaussian electrostatics

◮ The Boys function is related to the error function

erf(x) = 2 √π x exp

• −t2

dt = 2xF0

• x2

√π

◮ Introducing unit charge distributions and the reduced exponent as

ρp(rP) = p π 3/2 exp

• −pr2

P

• ,

α = pq p + q we obtain formulas similar to those of point-charge electrostatics (but damped) ρp (rP) rC dr = erf √pRPC

• RPC

ρp (r1P) ρq (r2Q) r12 dr1dr2 = erf √αRPQ

• RPQ

2 4 6 8 10 0.5 1.0 1.5 Trygve Helgaker (CTCC, University of Oslo) Spherical Coulomb integrals 11th Sostrup Summer School (2010) 18 / 34

The Boys function Fn(x)

◮ The Boys function is central to molecular integral evaluation

◮ it is evaluated by a combination of numerical techniques

◮ Fn(x) > 0 since the integrand is positive:

Fn(x) = 1 exp

• −xt2

t2n dt > 0

◮ Fn(x) is convex and decreasing:

F ′

n(x) = −Fn+1(x) < 0

F ′′

n (x) = Fn+2(x) > 0

2 4 6 8 10 12 0.2 0.4 0.6 0.8 1.0

◮ Evaluation for small and large values:

Fn(x) = 1 2n + 1 +

∞

• k=1

(−x)k k!(2n + 2k + 1) (x small) Fn(x) ≈ ∞ exp

• −xt2

t2ndt = (2n − 1)!! 2n+1

• π

x2n+1 (x large)

◮ Different n-values are related downward recurrence relations:

Fn−1(x) = 2xFn(x) + exp(−x) 2n − 1

Trygve Helgaker (CTCC, University of Oslo) Spherical Coulomb integrals 11th Sostrup Summer School (2010) 19 / 34

Cartesian Coulomb integrals

◮ We have obtained a simple result for one-center spherical Gaussians:

exp

• −pr2

P

• rC

dr = 2π p F0

• pR2

PC

• ◮ We shall now consider general two-center one-electron Coulomb integrals:

Vab =

• Ga
• r−1

C

• Gb
• =

Ωab(r) rC dr

◮ The overlap distribution is expanded in Hermite Gaussians

Ωab(r) = Ωij(x)Ωkl(y)Ωmn(z), Ωij(x) = i+j

t=0 E ij t Λt(xP)

◮ The total overlap distribution may therefore be written in the form

Ωab(r) =

• tuv

E ij

t E kl u E mn v

Λtuv(rP) =

• tuv

E ab

tuvΛtuv(rP)

◮ We now expand two-center Coulomb integrals in one-center Hermite integrals: two-center integral → Vab =

• tuv

E ab

tuv

Λtuv(rP) rC dr ← one-center integrals

◮ Our next task is therefore to evaluate Hermite Coulomb integrals

Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 20 / 34

Coulomb integrals over Hermite Gaussians

◮ Expansion of two-center Coulomb integrals in one-center Hermite integrals

Vab =

• tuv

E ab

tuv

Λtuv(rP) rC dr, exp

• −pr2

P

• rC

dr = 2π p F0

• pR2

PC

• ◮ Changing the order of integration and differentiation by Leibniz’ rule, we obtain

Λtuv(rP) rC dr = ∂t+u+v ∂Pt

x∂Pu y ∂Pv z

exp(−pR2

P)

rC dr = 2π p ∂t+u+vF0(pR2

PC )

∂Pt

x∂Pu y ∂Pv z

= 2π p Rtuv(p, RPC ) where the one-center Hermite Coulomb integral is given by Rtuv(p, RPC ) = ∂t+u+vF0

• pR2

PC

• ∂Pt

x∂Pu y ∂Pv z

◮ The Hermite Coulomb integrals are derivatives of the Boys function

◮ can be obtained by repeated differentiation, using F ′

n(x) = −Fn+1(x)

◮ recursion is simpler and more efficient Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 21 / 34

Evaluation of Hermite Coulomb integrals

◮ To set up recursion for the Hermite integrals, we introduce the auxiliary Hermite integrals

Rn

tuv(p, P) = (−2p)n ∂t+u+vFn

• pR2

PC

• ∂Pt

x∂Pu y ∂Pv z

which include as special cases the source and target integrals Rn

000 = (−2p)Fn → R0 tuv =

∂t ∂Pt

x

∂u ∂Pu

y

∂v ∂Pv

z

F0

◮ The Hermite integrals can now be generated from the recurrence relations:

Rn

t+1,u,v = tRn+1 t−1,u,v + XPC Rn+1 tuv

Rn

t,u+1,v = uRn+1 t,u−1,v + YPC Rn+1 tuv

Rn

t,u,v+1 = vRn+1 t,u,v−1 + ZPC Rn+1 tuv

Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 22 / 34

Summary one-electron Coulomb integral evaluation

1

Calculate Hermite expansion coefficients for overlap distributions by recurrence E 00 = exp

• −µX 2

AB

• E i+1,j

t

=

1 2p E ij t−1 + XPAE ij t + (t + 1)E ij t+1

2

Calculate the Boys function by a variety of numerical techniques Fn(x) = 1 exp

• −xt2

t2ndt (x = pR2

PC )

Fn(x) = 2xFn+1(x) + exp(−x) 2n + 1

3

Calculate Hermite Coulomb integrals by recurrence Rn

000 = (−2p)nFn

Rn

t+1,u,v = XPC Rn+1 tuv + tRn+1 t−1,u,v

4

Obtain the Cartesian Coulomb integrals by expansion in Hermite integrals

• Ga
• r−1

C

• Gb
• =
• tuv

E ij

t E kl u E mn v

R0

tuv

Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 23 / 34

Two-electron Coulomb integrals

◮ Two-electron integrals are treated in the same way as one-electron integrals ◮ Evaluation of one-electron integrals:

• Ga
• r−1

C

• Gb
• =

Ωab rC dr =

• tuv

E ab

tuv

Λtuv rC dr = 2π p

• tuv

E ab

tuvRtuv(p, RPC )

◮ Evaluation of two-electron integrals:

• Ga(a)Gb(1)
• r−1

12

• Gc(2)Gd(2)
• =

Ωab(1)Ωcd(2) r12 dr1dr2 =

• tuv

E ab

tuv

• τνφ

E cd

τνφ

Λtuv(1)Λτνφ(2) r12 dr1dr2 = 2π5/2 pq√p + q

• tuv

E ab

tuv

• τνφ

E cd

τνφ(−1)τ+ν+φRt+τ,u+ν,v+φ(α, RPQ)

◮ Hermit expansion for both electrons, around centers P and Q ◮ integrals depend on RPQ and α = pq/(p + q), with p = a + b and q = c + d Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 24 / 34

Direct SCF and other techniques

◮ The general expression for two-electron integrals is given by

gabcd = 2π5/2 pq√p + q

• tuv

E ab

tuv

• τνφ

E cd

τνφ(−1)τ+ν+φRt+τ,u+ν,v+φ(α, RPQ)

◮ In SCF theories, these integrals make contributions to the Fock/KS matrix:

Fab =

• cd
• 2gabcdDcd
• Coulomb

− gacbdDcd

• exchange
• ◮ in early days, one would write all integrals to disk and read back as required

◮ in direct SCF theories, integrals are calculated as needed ◮ this development (1980) made much larger calculations possible

◮ Many developments have since improved the efficiency of direct SCF

◮ screening of integrals ◮ early contraction with density matrices ◮ density fitting ◮ multipole methods Trygve Helgaker (CTCC, University of Oslo) Cartesian Coulomb integrals 11th Sostrup Summer School (2010) 25 / 34

Screening of overlap integrals

◮ The product of two s functions is by the Gaussian product rule

exp

• −ar2

A

• exp
• −br2

B

• = exp(− ab

a+b R2 AB) exp

• −(a + b)r2

P

• ◮ This gives a simple expression for the overlap integral between two such orbitals:

Sab =

• π

a + b 3/2 exp

• −

ab a + b R2

AB

• ◮ the number of such integrals scales quadratically with system size

◮ however, Sab decreases rapidly with the separation RAB

◮ Let us assume that we may neglect all integrals smaller than 10−k:

|Sab| < 10−k ← insignificant integrals

◮ We may then neglect integrals separated by more than

RAB >

• a−1

minln

• π

2amin 3 102k

• ◮ in a large system, most integrals becomes small and may be neglected

◮ the number of significant integrals increases linearly with system size Trygve Helgaker (CTCC, University of Oslo) Screening of integrals 11th Sostrup Summer School (2010) 26 / 34

Screening of two-electron integrals

◮ A two-electron ssss integral may be written in the form:

gabcd = erf √αRPQ SabScd RPQ

◮ the total number of integrals scales quartically with system size ◮ the number of significant integrals scales quadratically

Sab, Scd → 0 rapidly R−1

PQ → 0

very slowly

◮ Let us decompose the integral into classical and nonclassical parts:

gabcd = SabScd RPQ

classical

− erfc √αRPQ SabScd RPQ

• nonclassical

◮ quadratic scaling of classical part, can be treated by multipole methods ◮ linear scaling of nonclassical part since Sab → 0, Scd → 0, erfc → 0 rapidly

erfc(√αRPQ) ≤ exp(−αR2

PQ)

√παRPO

Trygve Helgaker (CTCC, University of Oslo) Screening of integrals 11th Sostrup Summer School (2010) 27 / 34

The scaling properties of molecular integrals

◮ linear system of up to 16 1s GTOs of unit exponent, separated by 1a0

3 5 7 9 11 13 15 110 3 310 3

nuclear attraction n3 n2 n1

3 5 7 9 11 13 15 110 3 310 3 3 5 7 9 11 13 15 210 4 510 4

electron repulsion n4 n2 n1

3 5 7 9 11 13 15 210 4 510 4 3 5 7 9 11 13 15 110 2 210 2

• verlap

n2 n1

3 5 7 9 11 13 15 110 2 210 2 3 5 7 9 11 13 15 110 2 210 2

kinetic energy n2 n1

3 5 7 9 11 13 15 110 2 210 2 Trygve Helgaker (CTCC, University of Oslo) Screening of integrals 11th Sostrup Summer School (2010) 28 / 34

Integral prescreening

◮ Small integrals (< 10−10) are not needed and should be avoided by prescreening ◮ The two-electron integrals

gabcd = Ωab(1)Ωcd(2) r12 dr1dr2 are the elements of a positive definite matrix with diagonal elements gab,ab ≥ 0

◮ The conditions for an inner product are thus satisfied ◮ The Cauchy–Schwarz inequality yields

• gab,ab
• ≤ √gab,ab√gcd,cd

◮ Precalculate

Gab = √gab,ab and prescreen

• gab,cd
• ≤ GabGcd

Trygve Helgaker (CTCC, University of Oslo) Screening of integrals 11th Sostrup Summer School (2010) 29 / 34

Early contraction with density matrix

◮ In direct SCF, two-electron integrals are only used for Fock/KS-matrix construction

Jab = 2π5/2 pq√p + q

• tuv

E ab

tuv

• τνφ
• cd

E cd

τνφ(−1)τ+ν+φRt+τ,u+ν,v+φ(α, RPQ)Dcd

◮ only integrals that contribute significantly are evaluated ◮ screening is made based on the size of the density matrix

◮ It is not necessary to evaluate significant integrals fully before contraction with Dcd ◮ An early contraction of Dcd with the coefficients E cd

τνφ is more efficient:

K Q

τνφ = (−1)τ+ν+φ cd∈Q

DcdE cd

τνφ

Jab = 2π5/2 pq√p + q

• tuv

E ab

tuv

• Q
• τνφ

K Q

τνφRt+τ,u+ν,v+φ(α, RPQ)

◮ Timings (seconds) for benzene:

cc-pVDZ cc-pVTZ cc-pVQZ late contraction 27 374 3823 early contraction 10 90 687

Trygve Helgaker (CTCC, University of Oslo) Early density-matrix contraction 11th Sostrup Summer School (2010) 30 / 34

Density fitting

◮ Traditionally, the electron density is expanded in orbital products

ρ(r) =

abDab Ωab(r)

◮ the number of terms is n2, where n is the number of AOs

◮ The Coulomb contribution to the Fock/KS matrix is evaluated as

Jab =

• Ωab(r1)
• r−1

12

• ρ(r2)
• ◮ the formal cost of this evaluation is therefore quartic (n4)

◮ Consider now an approximate density

ρ(r), expanded in an auxiliary basis ωα(r):

• ρ(r) =

α cαωα(r),

• ρ(r) ≈ ρ(r)

◮ the size of the auxiliary basis N increases linearly with system size

◮ We may now evaluate the Coulomb contribution approximately as

• Jab =
• Ωab(r1)
• r−1

12

• ρ(r2)
• ◮ the cost of the evaluation is therefore cubic (n2N)

◮ if the

ρ is sufficiently accurate and easy to obtain, this may give large savings

Trygve Helgaker (CTCC, University of Oslo) Density fitting 11th Sostrup Summer School (2010) 31 / 34

Coulomb density fitting

◮ Consider the evaluation of the Coulomb potential

Jab =

cd (ab |cd ) Dcd

◮ We now introduce an auxiliary basis and invoke the resolution of identity:

(ab |cd ) ≈

αβ (ab |α ) (α |β )−1 (β |cd )

◮ in a complete auxiliary basis, the two expressions are identical

◮ We may now calculate the Coulomb potential in two different ways:

Jab =

cd (ab |cd ) Dcd = (ab |ρ )

≈

cd

• αβ (ab |α ) (α |β )−1 (β |cd ) Dcd =

α (ab |α ) cα = (ab |

ρ ) = Jab where the coefficients cα are obtained from the linear sets of equations

• α (β |α ) cα =

cd (β |cd ) Dcd

⇔ (β | ρ ) = (β |ρ )

◮ no four-center integrals, the formal cost is cubic or smaller ◮ note: the solution of the linear equations scales cubically

◮ This particular approach to density fitting is called Coulomb density fitting

Trygve Helgaker (CTCC, University of Oslo) Density fitting 11th Sostrup Summer School (2010) 32 / 34

Robust density fitting

◮ In Coulomb density fitting, the exact and fitted Coulomb matrices are given by

Jab = (ab |ρ ) ,

• Jab = (ab |

ρ ) , (α |ρ ) = (α | ρ )

◮ We may then calculate the exact and density-fitted Coulomb energies as

E =

ab JabDab = (ρ |ρ ) ,

• E =

ab

JabDab = (ρ | ρ )

◮ The true Coulomb energy is an upper bound to the density-fitted energy:

E − E = (ρ |ρ ) − (ρ | ρ ) = (ρ |ρ ) − (ρ | ρ ) − ( ρ |ρ ) + ( ρ | ρ ) = (ρ − ρ |ρ − ρ ) where we have used the relation (˜ ρ|ρ) = (˜ ρ|˜ ρ) which follows from (α |ρ ) = (α | ρ )

◮ the error in the density-fitted energy is quadratic in the error in the density ◮ Coulomb fitting is thus equivalent to minimization of

E − E = (ρ − ρ |ρ − ρ ) ≥ 0

◮ It is possible to determine the approximate density in other ways

◮ quadratic energy error is ensured by using the robust formula

• E = (

ρ |ρ ) + (ρ | ρ ) − ( ρ | ρ ) = E − (ρ − ρ |ρ − ρ )

◮ however, the error in the energy is no longer minimized Trygve Helgaker (CTCC, University of Oslo) Density fitting 11th Sostrup Summer School (2010) 33 / 34

Density fitting: sample calculations

◮ Calculation on benzene

cc-pVDZ cc-pVTZ cc-pVQZ exact density late 27 374 3823 exact density early 10 90 687 fitted density 1 7 34 exact energy −230.09671 −230.17927 −230.19463 density-fitted energy −230.09733 −230.17933 −230.19466

◮ Clearly, very large gains can be achieved with density fitting ◮ Formal scaling is cubical

◮ for large systems, screening yields quadratic scaling for integral evaluation ◮ the cubic cost of solving linear equations remains for large systems

◮ Linear scaling is achieved by boxed density-fitting and fast multipole methods

◮ the density is partitioned into boxes, which are fitted one at a time ◮ fast multipole methods for matrix elements Trygve Helgaker (CTCC, University of Oslo) Density fitting 11th Sostrup Summer School (2010) 34 / 34