A Gaussian Implementation of Yang's Divide-and-Conquer Density

Chapter 5

A Gaussian Implementation of Yang's Divide-and-Conquer Density-Functional Theory Approach

Downloaded by COLUMBIA UNIV on March 16, 2013 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch005

Alain St-Amant Department of Chemistry, University of Ottawa, 10 Marie Curie, Ottawa, Ontario K1N 6N5, Canada

A gaussian implementation of Yang's divide-and-conquer approach to density functional theory has been created. Divide-and-conquer ap proaches to the fits of the density and the exchange-correlation po tentials have been developed. T h e concept of extended buffer space is introduced. It extends the spatial extent of buffer space while keep ing the number of basis functions under control. Tests on dipeptide, tripeptide, and tetrapeptide analogues are presented. T h e results sug gest that we w i l l have to go to much larger systems before our gaussian divide-and-conquer approach outperforms the conventional gaussian density functional approach.

C o m p u t a t i o n a l chemists are increasingly t u r n i n g towards density functional the ory ( D F T ) [1] to perform caclulations previously done by conventional ab initio methods [2]. D F T allows us to treat systems much larger than those currently feasible w i t h i n any correlated post-Hartree-Fock approach. Nevertheless, t r u l y large systems are beyond conventional D F T methods. In theory, we could work directly w i t h the electronic density, a tremendous computational simplification. However, i n practise, we must make use of molecular orbitals, and D F T schemes thus scale between N and N , where Ν is related to system size. Y a n g has proposed a divide-and-conquer ( D A C ) approach offering linear scal ing [3]. E l i m i n a t i n g the need for molecular orbitals that span the entire molecule, it brings us closer to an approach dealing directly w i t h the density. T h e D A C ap proach has been implemented [3] w i t h i n a program similar to D M o l [4]. W e wish to implement a D A C scheme w i t h i n the linear combination of gaussian type or bitals ( L C G T O ) - D F T formalism, originally developed by D u n l a p , Connolly, and Sabin [5], popularized by such programs as DGauss [6] and D e F T [7]. T h e fitting procedures that make L C G T O - D F T so efficient are not found w i t h i n a D M o l - t y p e scheme. W e must address this point w i t h i n our gaussian D A C approach. 2

3

0097-6156/96/0629-0070$15.00/0 © 1996 American Chemical Society In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.

5.

ST-AMANT

71

Yang's Divide-and-Conquer DFT Approach

For the systems of interest to us, we demand a higher level of accuracy t h a n that found i n previous D A C schemes. We would like to see our D A C scheme introduce errors, versus conventional gaussian D F T , no greater t h a n 0.2 kcal m o l " i n relative conformational energies of small peptides. It is important to note that this 0.2 kcal m o l " threshold is not the error relative to experiment, but rather to conventional calculations where the only difference lies i n the fact that the D A C approach is not used. We wish to establish a well-defined D A C protocol that w i l l consistently give us this k i n d of precision so that we can be assured that any major errors are a result of the D F T method and not the D A C implementation. 1

1


T h e K o h n - S h a m Equations In the K o h n - S h a m ( K S ) approach [8] to D F T , the density, p(r), is expressed as the sum of the square m o d u l i of N doubly-occupied (we could readily generalize this to the spin-polarized open shell case) K S molecular orbitals, occ

(1)

p(τ) = 2Σ\Ψ>(τ)\ . ί

»

W i t h i n the K S approach, the electronic energy is partitioned as follows,

Ε [p(τ)\ = T lp(r)) + U [p(r)) + E

xc

(2)

[p(r)}.

Τ \fi(r)] is the kinetic energy of non-interacting electrons,

Τ W)\ =

2

Σ

/ *i(r)^*(r)dr.

(3)

U \p(r)] is simply the classical coulomb energy,

where {ZA} and { R A } are the nuclear charges and coordinates. E [p(r)] is made to contain a l l the remaining effects of exchange and correlation ( X C ) . T h e K S molecular orbitals, {^-(r)}, are obtained by solving the K S equations, XC

Ηφ^τ) = e,^(r),

(5)

where Η is the K S operator. It is given by

where v (r) is the X C potential, xc

(7)

In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.

72

C H E M I C A L APPLICATIONS O F DENSITY-FUNCTIONAL T H E O R Y

T h e t o t a l energy can be expressed i n terms of the K S eigenvalues, {ε,·}, w

~

1

· 'p(rMr') nuclei

/ Κ Γ Μ # +£ [ΗΓ)]+Σ

7 7τ> Α

Σ 4

Κ

·

1

(8)

In the above expression, the s u m of the K S eigenvalues is first corrected for the double-counting of coulomb repulsion. It is then corrected for the fact that H contains v ( r ) . Its contribution is subtracted out and replaced b y E [p(r)]. T h e final t e r m is the t r i v i a l internuclear repulsion t e r m . xc


xc

T h e Conventional Gaussian Density Functional A p p r o a c h In the conventional L C G T O - D F T approach [5, 6], the K S orbitals are expressed as a linear combination of Ν contracted gaussian basis functions,

^) =Σ

Μ

·

(»)

μ T h e coefficient m a t r i x , C , is obtained by solving (10)

H C = SCe, where the m a t r i x elements Η

μν

* ~ = + Χ Ϊ Λ Ζ \

+

and 5

μ ι /

are given by

J 0 ? * '

+

-

V

(

R

)

}

X

-

(

R

)

D

R

and S„„ = / x „ ( r ) , ( r ) d r .

(12)

X l

ε is the vector of K S eigenvalues. C o m b i n i n g equations (1) and (9), p(r) is now Ν

P(r) = 2 E

=

2

Ν N occ t

Σ Σ

μ B a r r i n g any simplification, Η

μν

ν

Σ

^ ,α,Χ (Γ)χ,(Γ). μ

Μ

(13)

χ

w i l l require the evaluation of four-center two-

electron integrals, leading to a formal N scaling. 4

However, p(r) is fit b y a n

auxiliary set of M uncontracted gaussians, {y>*(r)}, M

p(r) « p(r) = Σ «k we must first synthesize />(Ri), v i a equation (13). T h i s scales as N . Since the number of grid points scales linearly w i t h system size, as does N, this procedure over the entire set of grid points also scales as N . Inserting v (r) into equation (11), Η requires the evaluation of three-center overlap integrals, yet another N step. xc

2

3

μι/

xc

3

D i v i d e - a n d - C o n q u e r Density Functional T h e o r y It is apparent from equation (3) that the kinetic energy component of the total energy requires the use of molecular orbitals delocalized over the entire extent of a molecule. Y a n g has recently proposed a D A C approach [3] that brings us closer to true density functional approaches where we work directly w i t h p(r) and eliminate the need for molecular orbitals. A rigorous theoretical justification of the method is given elsewhere [3]. Here we w i l l focus on describing the approach w i t h i n a linear combination of atomic orbitals ( L C A O ) context. A molecule may be divided into a set of N i, chemically intuitive subsystems. For example, the alanine dipeptide analogue, Figure I, may be partitioned into three C H fragments, two N H fragments, two C O fragments, and a C H fragment. 8U

3


74

C H E M I C A L APPLICATIONS O F DENSITY-FUNCTIONAL T H E O R Y

T h e electronic density can be arbitrarily partitioned i n the following fashion,

P(r) =

Σ AO

( ) 18

a

where N

N i,

eub

P (r) = Σ a

Nocc

su

P (r)/»(r) = 2 £

p«(r) £

a

a

a

(19)

Μ*)ΐ,

i

provided that p (r) is a function that obeys, at each point i n space, the constraint Downloaded by COLUMBIA UNIV on March 16, 2013 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch005

a

N.

ub

= 1· (20) or W e can now change the summation over doubly-occupied orbitals i n equation (19) to one over the J V doubly-occupied orbitals and the Nw empty v i r t u a l orbitals, Σ Ρ Ή

r

occ

Nsub

Nocc+N ir v

P(*) = 2 Σ P"( )

(21)

v(e -ei)Mr)\ .

Σ

r

a

2

F

i

In the above expression, 6F is the Fermi energy and η is a Heavyside function equal to one when ε, < SF and zero when ε EFF r o m equation (21), we see that i f we make p (r) go to zero as we move away from the atoms i n subsystem a , we need only know φ%(τ) i n the v i c i n i t y of subsystem a . T h i s can be accomplished w i t h the following definition for p ( r ) , a

a

^

Μ

œ'[pï (\r-RA\)f

_

tm

Σ^Μ^(ΙΓ-ΚΒΙ)]

2

where p (\r — R A | ) is the spherical atomic density of atom A , centered at R A , and the summations i n the numerator and the denominator r u n over a l l the atoms i n subsystem α and over a l l the atoms i n the molecule, respectively. H a v i n g built p (r) to localize the contributions of {V>i(r)}, we begin t o imple ment a D A C approach and approximate the true {^t(r)} by subsystem orbitals, {V>?(r)}. A modified version of equation (10) is solved for each subsystem a , A

tom

a

HC a

a

= SCe, Q

a

(23)

a

where H and S contain the elements of Η and S where both basis functions are part of subsystem a . W i t h C , each φ?(τ) is expanded, v i a equation (9), w i t h i n subsystem a's N basis functions. Plugging the subsystem orbitals into equation (21), p(r) is approximated by a

a

a

Q

p(r) « p(r) = 2 £

p"(r) £

(e

v

F

- ef ) |#»(r)| . 2


(24)

5.

ST-AMANT

75

Yang's Divide-and-Conquer DFT Approach

T o avoid oscillations i n the course of a calculation and achieve self-consistency, i t is necessary to work w i t h a unique value of 6F- W e must thus further approximate p(r) b y replacing the Heavyside function w i t h a F e r m i function

p(r) = 2 Σ p«(r) Σ

^ _ ,

?

|tf (r)|»

)

(25)

where β is an adjustable parameter that makes the F e r m i function approach the Heavyside function as i t is increased. Fortunately, the final results are fairly insensitive to its precise value and we have s i m p l y used the previously suggested value [3]. F i n a l l y , 6F is chosen such that the approximate D A C density, />(r), is normalized to the number of electrons i n the molecule, J V ,


e

/ p(r)dr = 2 Ç Ç

1

+

e

_,

( e F

- ) /IV>f(r)| e ?

(26)

d r = N..

2

T h e divide-and-conquer total energy,

E[p(r)] = 2 £

£

E

- L - e ? )

4/ / n7^

d

r

d

° / ? »

£

r

'

W W ' *

" / ^

r

> ^

d

+ *
Α

(28)

Γ

k

where the M auxiliary basis functions, {jt(r)}, still span the entire molecule. T h e double tilde notation is used to indicate that we are approximating the true L C G T O - D F T p(r) i n two ways: we are adopting a D A C approach to get a n approximate density which we i n t u r n fit. O u r scheme is the following. W i t h i n each subsystem, p (r) is constructed, a

(29)

Σ^Χ,(Γ)

p (r) is fit w i t h i n the auxiliary bases of both the subsystem and buffer atoms, a

M

a

(30)

* ( r ) « ^ ( r ) = £*(r). k

M is the number of auxiliary functions i n the subsystem and its associated buffer. Unfortunately, the fitting procedure can no longer be variational, nor analytical. T h e complexity of the expression for p (r) does not permit an analytical approach. W e must go to a numerical approach where the |r — r ' l " " factor present i n equa tion (15) is no longer feasible. W e thus lose the variational aspect of the fitting procedure. W e must perform a straightforward numerical fit of /> (r), m i n i m i z i n g a

a

1

a

points

Σ

[ P W - ^ R I ) ]

2

(31)

^

on a grid spanning the entire molecule. However, p (r) localizes p (r) to a small region of space, thus e l i m i n a t i n g a large number of points. Once a l l the subsystems fit, we s i m p l y sum the subsystem fit coefficient vectors. T h e i r s u m , a , is then used to construct the fit to the molecule's total density, a

a

T

N.ub

Χ*) = Σ

A Gaussian Implementation of Yang's Divide-and-Conquer Density

Recommend Documents