Algorithms for Chemical Computations - American Chemical Society

2 Algorithm Design in Computational Quantum

Downloaded by UNIV OF MASSACHUSETTS AMHERST on October 5, 2015 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0046.ch002

Chemistry ERNEST R. DAVIDSON Chemistry Dept., University of Washington, Seattle, WA 98195

Quantum chemistry is a diverse d i s c i p l i n e which uses many different methods to correlate a wide variety of phenomena. In the earliest period of the subject the Schrödinger equation was solved exactly for a few simple model situations. These model solutions were then used to interpret the spectra, kinetics, and thermodynamics of molecules and solids. During this period, accurate solutions for the electronic structure of helium (1) and the hydrogen molecule (2) were obtained i n order to verify that the Schrödinger equation was useful. Most of the e f f o r t , however, was devoted to developing a simple quantum model of electronic structure. Hartree (3) and others developed the self-consistent-field model for the structure of light atoms. For heavier atoms, the Thomas-Fermi model (4) based on t o t a l charge density rather than individual orbitals was used. Models for the electronic structure of polynuclear systems were also developed. Except for metals, where a free electron model of the valence electrons was used, all methods were based on a description of the electronic structure i n terms of atomic o r b i t a l s . Direct numerical solutions of the Hartree-Fock equations were not feasible and the Thomas-Fermi density model gave ridiculous results. Instead, two different models were introduced. The valence bond formulation (5) followed closely the concepts of chemical bonds between atoms which predated quantum theory (and even the discovery of the electron). In this formulation certain reasonable "configurations" were constructed by drawing bonds between unpaired electrons on different atoms. A mathematical function formed from a sum of products of atomic o r b i t a l s was used to represent each configuration. The energy and electronic structure was then 21 In Algorithms for Chemical Computations; Christoffersen, R.; ACS Symposium Series; American Chemical Society: Washington, DC, 1977.


22

ALGORITHMS

FOR

CHEMICAL

COMPUTATIONS

found by the l i n e a r v a r i a t i o n method ( a l s o c a l l e d "resonance" or " c o n f i g u r a t i o n i n t e r a c t i o n " ) . Because of i t s almost one-to-one correspondance w i t h e a r l i e r chemical concepts the valence bond model gained widespread acceptance ( 6 ) . The molecular o r b i t a l model (7) assumed, instead, that the e l e c t r o n s were i n c e r t a i n molecular o r b i t a l s which could be expressed as l i n e a r combinations of atomic o r b i t a l s . C o n f i g u r a t i o n s were then cons t r u c t e d as various ways of arranging e l e c t r o n s i n o r b i t a l s . The molecular o r b i t a l model gave a c l e a r i n t e r p r e t a t i o n of molecular spectra but was l e s s transparent than the valence bond method i n modeling geometrical s t r u c t u r e of molecules (6,8). In almost a l l e a r l y app l i c a t i o n s of valence bond (9) and molecular o r b i t a l (10) models the i n t e g r a l s encountered were too d i f f i c u l t to a c t u a l l y evaluate so e m p i r i c a l values of the i n t e g r a l s were assumed which reproduced the phenomena being studied. With the advent of the stored-program d i g i t a l computer a minor r e v o l u t i o n occurred i n quantum chemistry. The i n t e g r a l s appearing i n the models being used f o r small molecules were a c t u a l l y evaluated and i t became c l e a r that molecules were enormously more complicated than had been a n t i c i p a t e d . The o v e r s i m p l i f i e d valence bond and molecular o r b i t a l methods o f t e n gave q u a l i t a t i v e l y r i d i c u l o u s r e s u l t s when taken l i t e r a l l y (11). As a consequence of these negative r e s u l t s , the f i e l d of ab i n i t i o quantum chemistry developed w i t h the goal of f i n d i n g computer algorithms f o r s o l v i n g the Schrôdinger equation. The prospect of o b t a i n i n g r e l i able r e s u l t s f o r molecular systems not s u s c e p t i b l e t o d i r e c t measurement ( r e p u l s i v e p o t e n t i a l energy surfaces, upper atmosphere free r a d i c a l s , etc.) and c l a r i f y i n g the i n t e r p r e t a t i o n of experimental r e s u l t s which do not f o l l o w simple models a t t r a c t e d i n t e r e s t i n t h i s f i e l d i n s p i t e of the e x t r a o r d i n a r y expense of the approach and the lack of chemical i n s i g h t i n the e a r l y r e s u l t s . In the ab i n i t i o approach the d e s i r e d answers are the experimental observables - s p e c t r a l l i n e p o s i t i o n s , shapes, i n t e n s i t i e s ; s c a t t e r i n g and r e a c t i o n r a t e s ; p o l a r i z a b i l i t i e s and o p t i c a l r o t a r y power; etc. These are t o be obtained from the Schrôdinger equation by numerical methods which are mathematically w e l l - d e f i n e d and i n v o l v e no intermediate parameters not appearing i n the Schrôdinger equation i t s e l f . Usually the Born-Oppenheimer separation of nuclear and e l e c t r o n i c coordinates i s assumed and small terms i n the hamiltonian, such as s p i n - o r b i t coupling, are neglected i n the f i r s t approximation. P e r t u r b a t i o n

In Algorithms for Chemical Computations; Christoffersen, R.; ACS Symposium Series; American Chemical Society: Washington, DC, 1977.


2.

DAVIDSON

Computational Quantum Chemistry

23

theory may be used to c o r r e c t f o r these approximations by coupling e l e c t r o n i c s t a t e s i n the next l e v e l of approximation. Figure 1 o u t l i n e s the r e l a t i o n s h i p between various steps i n the c a l c u l a t i o n of some e x p e r i mental observables. C e n t r a l t o a l l other steps i s the c a l c u l a t i o n of the a d i a b a t i c e l e c t r o n i c wavefunctions f o r a l l s t a t e s of i n t e r e s t . From the wavefunctions one can obtain f i r s t order p r o p e r t i e s and coupling matrix elements f o r e s t i m a t i n g c o r r e c t i o n s due t o coupling of s t a t e s by non-adiabatic o r s p i n - o r b i t e f f e c t s . Methods which by-pass the wavefunction such as Χα or d e n s i t y f u n c t i o n a l models (12) are not yet s u f f i c i e n t l y general to t r e a t t h i s wide c l a s s of chemical problems. Each box i n Figure 1 represents i t s own p e c u l i a r computing problems. The algorithms f o r v a r i o u s steps are at various l e v e l s of s o p h i s t i c a t i o n depending on the r e l a t i v e cost, d i f f i c u l t y , and i n t e r e s t i n the r e s u l t s . The i n i t i a l c a l c u l a t i o n of e l e c t r o n i c wavef u n c t i o n s and energy surfaces have preoccupied quantum chemists f o r t h i r t y years. The c a l c u l a t i o n of a d i a b a t i c s c a t t e r i n g and r e a c t i o n rates has received much a t t e n t i o n i n recent years (JJ3). The accurate c a l c u l a t i o n of v i b r a t i o n a l - r o t a t i o n a l l e v e l s i s nearly as d i f f i c u l t but has received l i t t l e a t t e n t i o n u n t i l very r e c e n t l y . Equally accurate formalisms i n the coupled s t a t e model do not e x i s t because no general a l g o r i t h m e t r i c formalism e x i s t s f o r handling the e l e c t r o n i c part of the problem. No v i b r a t i o n a l - r o t a t i o n a l spectrum has yet been computed from an ab i n i t i o approach t a k i n g f u l l account of BornOppenheimer coupling i n a Jahn-Teller-Renner s i t u a t i o n . Generally speaking the whole area of coupled e l e c t r o n i c s t a t e c a l c u l a t i o n s l a c k s a workable algorithm. F i r s t order p e r t u r b a t i o n theory, while suggestive, i s o f t e n not a q u a n t i t a t i v e t o o l . The r e s t of t h i s paper w i l l deal e x c l u s i v e l y w i t h algorithms f o r c o n s t r u c t i o n of e l e c t r o n i c wavefunctions because these are c e n t r a l t o the o v e r a l l problem. In order t o appreciate the methods used, one must r e c a l l that we are i n t e r e s t e d i n s o l v i n g a p a r t i a l d i f f e r e n t i a l equation eigenvalue problem f o r s e v e r a l wavefunctions at s e v e r a l d i f f e r e n t arrangements of the n u c l e i . This d i f f e r e n t i a l equation i n v o l v e s one- and two-body opera t o r s i n the p o t e n t i a l energy operator and p a r t i a l d e r i v a t i v e s w i t h respect t o 3N coordinates (where Ν i s the number of e l e c t r o n s ) . For benzene, f o r example, there are 12 n u c l e i and 42 e l e c t r o n s . The reasonable a s p i r a t i o n of f i n d i n g the e q u i l i b r i u m geometry and f o r c e constants f o r the f i r s t 10 s t a t e s would i n v o l v e s o l v i n g a p a r t i a l d i f f e r e n t i a l


24

ALGORITHMS FOR CHEMICAL COMPUTATIONS


COUPLED STATE REACTION RATE PERTURBED VIBRATIONALROTATIONAL L E V E L S , JAHNTELLER-RENNER EFFECTS NUCLEAR-ELECTRONIC COUPLED MOTION

GET APPROXIMATE POTENTIAL SURFACES AND E L E C T R O N I C WAVEFUNCTIONS FOR STATES AND GEOMETRIES OF INTEREST USING BORN-OPPENHEIMER APPROXIMATION A N D ONLY COULOMB I N T E R A C T I O N S

NATURAL ORBITALS

ι CHEMICAL INTER•* PRETATION

DENSITY MATRICES

F I R S T ORDER CORRECTIONS TO ENERGY DISTRIBUTION OF SPIN, CHARGE A N D MOMENTUM Figure 1.

CORRECT WAVEFUNCTIONS FOR P E R T U R B A T I O N S (SPIN-ORBIT, EXTERNAL FIELD, RELATIVISTIC, ETC.) WITHIN BORN-OPPENHEIMER APPROXIMATION

— NUCLEAR MOTION ADIABATIC REACTION RATES

POLARIZABILITY

"U" E L E C T R O N I C

TRANSITION RATES AND L I F E T I M E S

VIBRATIONAL AVERAGED PROPERTIES

Flow chart jot ab initio calculions



2.

DAVIDSON

25


equation of t h i s type i n 126 independent v a r i a b l e s . The only reason i t i s p o s s i b l e here i s that (1) the f i x e d f i e l d due to the n u c l e i dominates over the e l e c t r o n - e l e c t r o n r e p u l s i o n so the e l e c t r o n i c motions are u s u a l l y not s t r o n g l y coupled to each other, (2) i t i s impossible f o r a large c o l l e c t i o n of mutuallyr e p u l s i v e p a r t i c l e s t o avoid each other i f they are constrained to remain i n the same region of space, and (3) e l e c t r o n s are i n d i s t i n g u i s h a b l e so the coordinates are permutational equivalent. Hence the antisymmetric independent p a r t i c l e approximation which leads t o a pseudo-separation of v a r i a b l e s i s often a good f i r s t approximat ion. Now consider the resources a v a i l a b l e f o r s o l v i n g t h i s (or a s i m i l a r ) problem i f some government agency decides these r e s u l t s are v i t a l to the n a t i o n a l welfare. I t would then be p o s s i b l e t o spend up to 1 0 hours of CDC7600 time on t h i s problem (about $10,000,000). This w i l l allow about 1 0 a r i t h m e t i c operations ( a d d i t i o n o r ^ m u l t i p l i c a t i o n ) . Also we can assume^that at most 10 words of high speed core memory, 10 words of low^ speed core, 1 0 words of d i s k or drum storage, and 10 words of s e q u e n t i a l tape storage are a v a i l a b l e . By present standards t h i s would be a very large c a l c u l a ^ t i o n s i n c e every member given here i s a f a c t o r of 10 l a r g e r than what i s t y p i c a l l y used. If one wavefunction at one set of nuclear c o o r d i nates were sought by numerical i n t e g r a t i o n using only two p o i n t s i n each coordinate, a g r i d of 2 ^ 10 p o i n t s would be r e q u i r e d . I f spin and antisymmetry are taken i n t o account the s i t u a t i o n i s even worse. Since no two e l e c t r o n s can be at the same point with the same spin at l e a s t Ν p o s i t i o n s must be considered f o r each e l e c t r o n and the minimum g r i d contains 42! = 1 0 points i n 3N space. The only method found so f a r which i s f l e x i b l e enough to y i e l d ground and e x c i t e d s t a t e wavefunctions, t r a n s i t i o n r a t e s and other p r o p e r t i e s i s based on ex panding a l l wavefunctions and operators i n a f i n i t e d i s c r e t e set of b a s i s f u n c t i o n s . That i s , a set of onep a r t i c l e s p i n - o r b i t a l s {φ.}ξ are s e l e c t e d and the wavefunction i s expanded i n S l a t e r determinants based on these o r b i t a l s . A d i r e c t expansion would r e q u i r e w r i t i n g Ψ as 4

1 4

8

1 2 6

3 8

5 1

=1

ψ = ι

Φ

οΦ τ

τ

τ

= det(

dT

and f

Q i j k * • / V-fel2- q +

p

a

/

q

(

R

) ] / a

i

W

R-p—R-q

Hence [ij||k£] - K . . K J d T f ( r ) r - l / d T k

1 2

t

f Î2 d T

f

( r

r

/ sined6/ ^exp(-a r2 ) ο ο 7 r

= /°°cir r Ο

2

t

12

2

χ expi-a^R^)exp(-2a r R cos9) =

33-Γ = Γ

u

2 n

du ο S i m i l a r l y , S l a t e r o r b i t a l s f o r diatomic molecules give i n t e g r a l s of the form (22) 1

e~ u

n

2

( u - ! ) * du

(23) " W "

-

l\

. - " V P V u )

"" 2

2 +

'du


2.

33


DAVIDSON

Rather elaborate r e c u r s i o n r e l a t i o n s can be found f o r a l l these i n t e g r a l s when care i s taken to preserve numerical accuracy. Since u s u a l l y a l l values of η are needed anyway, the intermediate values of η as w e l l as the l a r g e s t value n=N and the smallest n=0 are u s e f u l . For example, the r e c u r s i o n r e l a t i o n (2n+l)F (T) = 2T F Downloaded by UNIV OF MASSACHUSETTS AMHERST on October 5, 2015 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0046.ch002

n

Q + 1

( T ) + e"

T

i s s t a b l e f o r r e c u r r i n g downward on η but i s unstable f o r r e c u r r i n g upwards ( f o r small T/n) because (2n+l)F (T) * e-T. Consequently, e v a l u a t i o n of F ( T ) i n v o l v e s d i f f e r e n t schemes depending on the value of Ν and T. For Τ > N, upward recurrence from F i s p o s s i b l e without l o s s of s i g n i f i c a n t f i g u r e s . For Τ < N, downward recurrence must be used s t a r t i n g from F ( T ) . For most f u n c t i o n s t h i s s i t u a t i o n would r e q u i r e e i t h e r a set of t a b l e s f o r every p o s s i b l e s t a r t i n g value of Ν or e l s e one t a b l e f o r an N* greater than any Ν which can occur followed by downward recurrence from N*. The p a r t i c u l a r f u n c t i o n F d e a l t w i t h here, however, obeys the r e l a t i o n s h i p n

Q

N

é

F

T

T

n< > =

"W >

so the Taylor s e r i e s has the simple form

V > = j W o > / T

T

T

k

k!

0

The convergence r a t e of t h i s s e r i e s i s nearly independent of η ( F k + i / F k « 1 f o r small T) so a t a b l e of F ( T ) at a sequence of i n t e r v a l s of T f o r η from zero N*+K (an i n t e r v a l width of 0.1 r e q u i r e s a Κ of 6 f o r twelve s i g n i f i c a n t f i g u r e s ) s u f f i c e s f o r a l l values of Ν and T. As f o r F , at l a r g e Τ i t i s b e t t e r to evaluate a generalized error function n +

n

n +

Q

Q

Q

G (/T) = n

G

n+1

( / Î )

=

u 2

2 n

e" u du Ο (*+4)G (/T) - T n

N

+

V

T

Hence an e f f i c i e n t algorithm must recognize s e v e r a l ranges of T:



34

ALGORITHMS

Τ = 0

F

n

0 2

lz \ v ) i j i-i...H .J-i...H D i r e c t e v a l u a t i o n i n t h i s form r e q u i r e s d^d2did2 m u l t i p l i c a t i o n . On the other hand the m u l t i p l i c a t i o n Υ = B W ν followed by W , Y r e q u i r e s only z( LL2) I 2 ) ( 2) _ _ _ l 2 2 l l 2 m u l t i p l i c a t i o n s (or d ^ d ] ^ l 2 2 i f the m u l t i p l i c a t i o n s are done i n the opposite order). Figures 3,4 show an o u t l i n e of algorithms f o r the t r i a n g u l a r and rectangular cases f o r matrices small enough t o f i t e n t i r e l y i n t o high speed core. These algorithms are designed with one a d d i t i o n a l p r i n c i p l e i n mind. Namely, the only r e a l v a r i a t i o n between d i f f e r e n t ways of doing matrix m u l t i p l i c a t i o n i s the cost of indexing and amount of s c r a t c h storage used. Double s u b s c r i p t s should u s u a l l y be avoided and as f a r as p o s s i b l e matrix elements should be accessed s e q u e n t i a l l y . For t h i s reason i t i s best t o c a r r y out the rectangular transformation as Y = BÏ W * followed by Β „ =

S

ij •

^kiV^Vk,*

2

1

2

T

/ p v r

d

d

( i

+

d

N

/r

/T

N

r

d

r

d

+

r

-

d

d

d

/P

/ Γ

-Γ Γ -(Γ ) 1

2

N

-(ΓιΓ )

1

2

ΣΆ(Τ2)' Scratch storage i s reduced by using each column of Y as soon as i t i s formed t o do the second m u l t i p l i c a t i o n . The t r i a n g u l a r transformation i s f u r thur complicated by the f a c t that both Β and Β are stored i n a t r i a n g u l a r p a t t e r n which increases the complexity of indexing. Transformation of the two e l e c t r o n i n t e g r a l s i s a much more time consuming step. I f R(ί ,i ,in,i.)is the integral Ί

1

9

Δ

0

4

R(Vi ,i ,i Wf^^ 2

3

4

and R i s the transformed 5

1

1

integral

1

Γ

(

τ

άτ

(^' 2' 3' 4>"/βΐ ^1>*βΐ Ϊ2βΐ ^2>*8ΐ ΐ:2^ 1 2 1

2

3

4


42


ά^ά

INC =

I

2

ICNT = 0


For i = 1... cT|

I

LCNT = 0

For I = 1...d

d

D(£) =

2

l

I A(k+ICNT)*B(k+LCNT) k=l

LCNT = LCNT + d

1

JCNT = 0

ι

For j = Ï...INC i n steps o f cT|

d

x(j) =

2 I D(«*C(X,+JCNT) «,=1

JCNT = JCNT + d

ICNT = ICNT + d

2

]

τ

Figure 3.

Transformation of a real non-symmetric matrix, X = A BC T


2.

DAVIDSON


43

ITRI = 0 ICNT = 0 For i = 1... cT-j LTRI = 0 Downloaded by UNIV OF MASSACHUSETTS AMHERST on October 5, 2015 | http://pubs.acs.org Publication Date: June 1, 1977 | doi: 10.1021/bk-1977-0046.ch002

CLEAR D U ) TO ZERO FOR i = ^...d

]

For ϋ - 1...d I —

1

For k = 1...JM D U ) = D U ) + C(k+ICNT)*B(k+LTRI) D(k) = D(k) + CU+ICNT)*B(k+LTRI)

-

1

_ ( LTRI = LTRI+Jt D{SL) = D{i)

1

+ CU+ICNT)*B(LTRI)

i JCNT = 0 ITR = ITRI+i-1

I

For j = ITRI... ITR d

l

X( j+1 ) = I DU)*CU+JCNT) «,=1 JCNT = JCNT+d 1

ICNT = ICNT+d

1

ITRI = ITRI+i

Figure 4.

Transformation of a real symmetric matrix, X

T

= C BC


44

ALGORITHMS

FOR CHEMICAL

COMPUTATIONS

the R and R are r e l a t e d by a four-index (tensor) trans formation. R(J ,J ,J ,J ) = Σ 1

2

3

.

.

4

. (r


W

ΗV

u

w

r

i

( 4) 4J4

W

—

R ( i l , i 2

'

'

i 3

'

i 4 )

D i r e c t e v a l u a t i o n of t h i s f o u r - f o l d sum would r e q u i r e 4 d d n d2d2d.3d3d d4 m u l t i p l i c a t i o n s t o form a symmetry block of R i n t e g r a l s . By c o n s t r a s t , s e q u e n t i a l oneindex transformations 1

4

j ,i ,i ,i ) = Σ w* 1

2

3

4

Y(J ,J >J J ) 1

2

3 >

I

β

4

W

( r i ) i i J i

(r )i j 0

0

ZU^j^jg,^) = Σ W*

R(i ,i 1

X ( j 0

( r 3 ) i 3 J 3

l

, I

Y(J

2

2 >

, I

3

i ,i ) 3

, I

J ,i

1 )

2

4

4

3 )

)

i ) 4

3

R( J , 3 > J , J ) = I 2

x

3

4

1^

V j i j / ^ r ^ ' V V 4

4 4

r e q u i r e only d-id^d^dj^d^ + d d 2 d 2 d d + d d 2 d d d + d ^ d 2 d d d m u l t i p l i c a t i o n s . These transformations can be organized by t h i n k i n g of R ( i ^ , i , i , i ) f o r f i x e d i o i . as a matrix R( ^3^-4) which i s transformed l i k e a 1

3

4

3

4

1

3

3

4

4

2

3 4

3

4

11^2

one-body operator t o give (i i ) (i i ) Υ -ϊ ) Β W 3

4

+

+

3 3 4

4

( Γ ι

( r 2 >

.

(i i ) If the Y matrices are then reorganized t o give ( J 1 J 2 ) matrices by use of 3

4

Y

(io,i ) 4

1

Δ

ά

4

^ l '

J

2

_(J J ) n

^S'

2

4


DAVIDSON

2.

45


the R i n t e g r a l s can be formed from


U l V

The use of s i x d i f f e r e n t storage patterns f o r the two-electron i n t e g r a l s r e q u i r e s s i x d i f f e r e n t algorithms for c a r r y i n g out the transformation. Only the simplest (RRR) w i l l be presented i n d e t a i l here (.33). Since the number of i n t e g r a l s u s u a l l y exceeds the amount of high speed core a v a i l a b l e (and u s u a l l y low speed core as w e l l ) a transformation using minimum core w i l l be d i s cussed (assuming d i s k i s l a r g e enough t o hold one block of R ) . Suppose the i n t e g r a l s R(±i±2±^±^) are o r i g i n a l l y arranged so that R(1>1), R(2,1),.. appear i n s e q u e n t i a l order on a s e q u e n t i a l f i l e . The range of ( 1 3 1 4 ) can be blocked i n t o d d / n groups of s i z e n (with a smaller group at the end i f needed). Each group of n R matrices can then be transformed by a standard two sub s c r i p t transformation t o leave n Y matrices i n sequen t i a l order ( i n the same space i n core o r i g i n a l l y occupied by the R m a t r i c e s ) . Storage f o r the W matrices 3

4

3 4

3 4

3 4

3 4

and one s c r a t c h region f o r wT -(

3

r

5

4

vR^ " " ^ are needed i n Γχ)-

a d d i t i o n t o the space f o r the R arrays. The 3^3 sub s c r i p t s on each y ( 3 4 ^ array can a l s o be blocked i n t o d d / n 2 blocks of s i z e n and the Y arrays can be w r i t t e n t o d i s k i n blocks of s i z e n-^ ^Y 3 4 dom access f i l e . When a l l R matrices have been t r a n s formed, a block of Ϋ^1^2^ matrices i s e a s i l y formed i n core by reading a l l appropriate pieces from d i s k . The Y arrays can then be transformed by a standard two sub s c r i p t transformation and w r i t t e n t o a s e q u e n t i a l f i l e . This method r e q u i r e s d^d2n words of high-speed core for the i n i t i a l R arrays and d d n^2 words f o r the Y arrays. The intermediate random f i l e contains d d d q d / 1 2 3 4 blocks of s i z e n 3 4 which i s w r i t t e n and read only once. Maximum e f f i c i e n c y u s u a l l y r e q u i r e s making the product n 3 4 l a r g e as p o s s i b l e . Because t h i s i n t e g r a l transformation step involved d^ operations to transform d i n t e g r a l s i t has gained a r e p u t a t i o n as a bottleneck i n c a l c u l a t i o n s . A c t u a l l y , however, u n t i l d i s about 60 the formation of d4 i n t e g r a l s (over con t r a c t e d gaussian o r b i t a l s ) takes longer than the t r a n s 2

I

1

2

I

1

1 2

N

A

S

A

R

A

N

_

2

34

3

4

1

n

N

b

v

N

1 2

N

a

s

1 2

4


2

4

46



formation. For l a r g e r values of d i t i s l i k e l y that a CI matrix of l a r g e dimension w i l l be formed using these i n t e g r a l s ( o r a t h i r d or higher order p e r t u r b a t i o n c a l c u l a t i o n w i l l be done). U s u a l l y these uses of the i n t e g r a l s are more time consuming than t h e i r production so the transformation i s seldom the l i m i t i n g step. Eigenvalue algorithms. Matrix eigenvalue problems a r i s e i n quantum chemistry at both the SCF and CI l e v e l . The Roothaan SCF method r e q u i r e s s o l v i n g a non-ortho gonal eigenvalue problem of the dimension of the b a s i s set on each i t e r a t i o n f o r many of the eigenvalues and eigenvectors. The CI method u s u a l l y r e q u i r e s f i n d i n g the lowest few eigenvalues of a l a r g e matrix i n an orthonormal b a s i s of c o n f i g u r a t i o n s . Several algorithms e x i s t which are s u i t a b l e f o r f i n d i n g a l l of the eigenvalues of any matrix of dimen sion d which can be kept i n c e n t r a l memory. The Jacobi plane r o t a t i o n method i s by f a r the simplest t o program and i s reasonably e f f i c i e n t ( 3 4 ) . As i t i s an i t e r a t i v e method the running time cannot be r i g o r o u s l y defined, but times p r o p o r t i o n a l t o d are expected. Other methods u s u a l l y begin w i t h a n o n - i t e r a t i v e transforma t i o n t o t r i d i a g o n a l form followed by c a l c u l a t i o n of the eigenvalues and eigenvectors and a back transformation to the o r i g i n a l problem ( 3 4 , 3 5 ) . The time r e q u i r e d f o r the transformations i s p r o p o r t i o n a l to d w h i l e the time required to solve the t r i d i a g o n a l problem i s only pro portional to d . The Jacobi method i s g e n e r a l l y slower than these other methods unless the matrix i s nearly diagonal. In SCF c a l c u l a t i o n s one i s faced w i t h the non-orthogonal eigenvalue equation 3

3

2

F C = S C Λ where Λ i s the diagonal matrix of eigenvalues and C i s a matrix of eigenvectors. I f an o r t h o g o n a l i z i n g t r a n s f o r mation W i s known such that W SW=1, then T

T

1

T

W F W W"

C - W

1

S W W"

C A

or F

T

C

1

T

F

1

» W

= C A

where A

F W

and C » W C»


2.

47


DAVIDSON

Usually on the f i r s t i t e r a t i o n of an SCF c a l c u l a t i o n W i s computed by the Schmidt o r t h o g o n a l i z a t i o n method but t h e r e a f t e r W i s chosen t o be the C matrix from the pre vious i t e r a t i o n . This produces an F' matrix which i s nearly diagonal so the Jacobi method becomes q u i t e e f f i c i e n t a f t e r the f i r s t i t e r a t i o n . Further, i n the Jacobi method, F i s d i a g o n a l i z e d by an i t e r a t i v e sequence of simple p l a n e - r o t a t i o n transformations F / N = Χ (η)Ζ' (n)Ht(n V eigenvectors of F can thus be generated aê C = (···(WX( )X(2)*"'X(n)) which avoids the m u l t i p l i c a t i o n of W by C . A disadvantage of the Jacobi method i s that the e r r o r i n the eigenvector i s u s u a l l y p r o p o r t i o n a l t o the square root of the e r r o r i n the eigenvalues. Thus, i n 8 d i g i t a r i t h m e t i c , only 4 f i g u r e s can be obtained i n the eigenvectors. The inverse i t e r a t i o n method of Wilkinson (34) i s a method which gives f u l l accuracy i n the vectors. This method i s based on computing the eigenvector as (X1.-F")C = X where λ i s the eigenvalue and X i s a guess to the eigenvector. Because t h i s method r e q u i r e s s o l v i n g a d i f f e r e n t set of l i n e a r equations f o r each eigenvector i t i s only f e a s i b l e i f F" has an e a s i l y i n v e r t e d form ( s o l v i n g l i n e a r equations i s a d process unless the c o e f f i c i e n t matrix has some s i m p l i f y i n g f e a t u r e ) . I f F" i s t r i d i a g o n a l , then the time f o r each vector i s p r o p o r t i o n a l t o d so the time for d vectors i s p r o p o r t i o n a l t o d . In CI c a l c u l a t i o n s i t i s necessary t o f i n d a few s o l u t i o n s t o the matrix eigenvalue problem f

f

Τ

T

h

e

f i n a l


n + 1

v

x)

3

2

H C = λC 1 5 where Η i s of dimension from 10 t o 10 . For smaller dimensions i t i s most e f f i c i e n t t o use the standard t r i d i a g o n a l i z a t i o n r o u t i n e s . For matrices which are too l a r g e t o f i t i n t o high-speed core, s p e c i a l methods have been developed whose time per eigenvalue i s p r o p o r t i o n a l only t o the number of non-zero matrix elements ( d at most). These methods should be u s e f u l i n other areas of chemistry as w e l l . The f i r s t development i n t h i s area was the Nesbet method (36) f o r f i n d i n g the lowest (or highest) eigen value. This method was reorganized i n t o a b e t t e r algorithm by S h a v i t t (37) and then extended by S h a v i t t , et a l . (3j3) t o f i n d a few non-degenerate eigenvalues. Recently Davidson (39) has combined the fundamental ideas from Nesbet, Lanczos and inverse i t e r a t i o n schemes to form a method which works f o r the f i r s t few eigen values even i f they are degenerate. H i s method, however, 2

American Chemical Society Library 1155 10th St. N. W. In Algorithms for Chemical D. Computations; Christoffersen, R.; Washington, C. 20036 ACS Symposium Series; American Chemical Society: Washington, DC, 1977.

48

A L G O R I T H M S FOR

CHEMICAL

COMPUTATIONS

i n v o l v e s a l i t t l e more input-output than the Nesbet or S h a v i t t methods. The b a s i c concept of the Nesbet-Shavitt method i s based on i t e r a t i v e s e q u e n t i a l o p t i m i z a t i o n of the eigen vector elements. I f the q u a n t i t y p(C)=C HC/C C i s known for some C° and p p ( C ) i s below a l l the diagonal elements of H, then s e q u e n t i a l m i n i m i z a t i o n of p(C) with respect to each element C± [ i . e . s o l v i n g (3p/ac.) =0 and then stepping C? to the new C. before going to the next value of i ] gives T

0 = a

T

0


n o

6C. = C.-C°

= -[(Η-ρ1)0°]./(Η..-ρ)

where ρ = p(C° + e.ôC.) while f o r any value of 6C^ 26C,[(H-p°l)C°], + ( 6 C . ) [ H . . - p ° ] p(C°+e.6C.)-p(C ) i ^ ^ ^ C° C° + (6C ) + 2C°6C 2

0

1

1

±

Nesbet approximated

i

the optimum 6C^ by

6C. = -q^/CH..-ρ ) 0

where q° = [(H-p°l)C°]

JL

while S h a v i t t found 6C. from the s l i g h t l y more exact formula 1

6C.=-2q° / {Η..-ρ° + / ( Η . . - ρ ° ) -4q° [-q°+q(H .-p ) ]/C° C° } . 0

2

T

i:

Both of these formulas can be shown to give monotonie convergence f o r p . More importantly, S h a v i t t showed how use of the hermitian property of H could be used to w r i t e HC as (HÇ). = I H J

C

=

Σ H ji

C J

J

so that E±j and H j ^ d i d not both need to be stored and read from e x t e r n a l s t o r e . S h a v i t t et a l . f u r t h e r



2.


DAVIDSON

49

modified the Nesbet-Shavitt scheme to do excited states by introducing root-shifting and over-relaxation to speed convergences. Their method, however, often f a i l s to converge for nearly degenerate eigenvalues. Davidson introduced a different method for higher eigenvalues which also avoids the need to have the ele ments of H stored i n any particular order. In this method the k£h eigenvector of H for the η·^ iteration i s expanded i n a sequence of orthonormal vectors b i , ΐ=1···η with coefficients found as the k — eigenvector of the small matrix Η with elements b^HBj. Convergence can be obtained for a reasonably small value of η i f the expansion vectors b are chosen appropriately. Davidson defined ç a

( n )

-

k

( n )

l 4 l \ n)

- [H-p(c< )i]c£

n)

and chose b +i as the normalized residual when ξ/ ' was orthogonalized to the preceeding ]>]_··-b. This choice for b +l i s similar to the Nesbet choice (and also to f i r s t order perturbation theory and the inverse i t e r a tion method). By the excited state variation theorem, the k-tfi eigenvalue of H as i t i s sequentially bordered w i l l decrease monotonically to the k-ίΛ eigenvalue of H. Butscher and Kammer (40) have shown how a slight modification of this scheme which tracks on certain large elements of C rather than the index k can find a C with a certain desired pattern of coefficients without prior knowledge of the value of k and without finding any other eigenvectors. n

n

n

Literature Cited 1. 2. 3. 4. 5. 6.

Hylleraas, E.A., Z. Physik. (1930) 65, 209. James, H.M. and Coolidge, A.S., J. Chem. Phys. (1933) 1, 825. Hartree, D.R., Proc. Cambridge P h i l . Soc. (1928) 24, 89. Thomas, L.H. Proc. Cambridge P h i l . Soc. (1927) 23, 542. Slater, J.C., Phys. Rev. (1931) 38, 1109. Pauling, L., J. Am. Chem. Soc. (1931) 53, 1367.


50



7. 8. 9.

Mulliken, R.S., Phys. Rev. (1928) 32, 186. Walsh, A.D., J . Chem. Soc. (1953) 2260. Pauling, L. and Wheland, G.W., J . Chem. Phys. (1933) 1, 362. 10. Pariser, R. and Parr, R.G., J . Chem. Phys. (1953) 21, 466. 11. See for example Karo, A.M. and Allen, L.C., J. Chem. Phys. (1959) 31, 968. 12. Johnson, K.H., Adv. Quantum Chem. (1973) 7, 143. 13. See for example Kouri, D.J., "Energy Structure and Reactivity", Smith, D.W. and McRae, W.B., Eds., John Wiley & Sons, 1973. 14. Boys, S.F., Proc. Roy. Soc. (1950) A200, 542. 15. Hehre, W.J., Stewart, R.F. and Pople, J.A., J . Chem. Phys. (1969) 51, 2657. 16. Dunning, T.H., J . Chem. Phys. (1970) 53, 2823. 17. Shavitt, I., i n "Methods i n Computational Physics", Vol.2, Alder, B., Fernbach, S. and Rotenberg, Μ., Academic Press, 1963. 18. Huzinaga, S., Supp. Prog. Theoretical Phys. (1967) 52. 19. Shipman, L.L., and Christoffersen, R.E., Comp. Phys. Comm. (1971) 2, 201. 20. Elbert, S.T., and Davidson, E.R., J . Comput. Phys. (1974) 16, 391. 21. Whitten, J.L., J . Chem. Phys. (1963) 39, 349. 22. Rüdenberg, Κ., J . Chem. Phys. (1951) 19, 1459. 23. Corbato, F.J., J . Chem. Phys. (1956) 24, 452. 24. Roothaan, C.C.J., Rev. Mod. Phys. (1951) 23, 161. 25. Roothaan, C.C.J. and Bagus, P.S., i n "Methods i n Computational Physics", Vol. 2, Alder, Β., Fernbach, S., and Rotenberg, Μ., Eds., Academic Press, 1963. 26. Guest, M.F., and Saunders, V.R., Mol. Phys. (1974) 28, 219. 27. Hsu, H., Davidson, E.R. and Pitzer, R.M., J. Chem. Phys., J . Chem. Phys. (1976) 65, 609. 28. For a thorough review see Shavitt, I., i n "Modern Theoretical Chemistry", Vol. 2, Schaeffer, H.F. I I I , Ed., Plenum Press, New York, 1976. 29. Møller, Chr., and Plesset, M.S., Phys. Rev. (1934) 44, 618. 30. Elbert, S.T., and Davidson, E.R., Int. J . Quant. Chem. (1974) 8, 857. 31. Davidson, E.R., "Reduced Density Matrices i n Quantum Chemistry, Academic Press, New York, 1976. 32. Davidson, E.R., Int. J . Quant. Chem. (1974) 8, 83. 33. Elbert, S.T., Ab initio Calculations i n Urea, Ph.D. thesis, University of Washington, 1973.


2.

DAVIDSON

34.


35. 36. 37. 38. 39. 40.


51

Wilkinson, J.H., "The Algebraic Eigenvalue Problem", Clarendon Press, Oxford, 1965. Givens, J.W., J. Assoc. Comp. Mach. (1957) 4, 298. Nesbet, R.H., J. Chem. Phys. (1965) 43, 311. Shavitt, I, J. Comput. Phys. (1970) 6, 124. Shavitt, I., Bender, C.F., Pipano, A. and Hosteny, R.P., J. Comput. Phys. (1973) 11, 90. Davidson, E.R., J. Comput. Phys. (1975) 17, 87. Butscher, W. and Kammer, W.E., J. Comput. Phys. (1976) 20, 313.


Algorithms for Chemical Computations - American Chemical Society

Recommend Documents