Citation
DEVELOPMENT OF EFFICIENT ELECTRON CORRELATION METHODS FOR ONE- AND TWO-DIMENSIONAL EXTENDED SYSTEMS AND THEIR APPLICATIONS

Material Information

Title:
DEVELOPMENT OF EFFICIENT ELECTRON CORRELATION METHODS FOR ONE- AND TWO-DIMENSIONAL EXTENDED SYSTEMS AND THEIR APPLICATIONS
Copyright Date:
2008

Subjects

Subjects / Keywords:
Atoms ( jstor )
Conceptual lattices ( jstor )
Electrons ( jstor )
Energy ( jstor )
Geometry ( jstor )
Molecules ( jstor )
Orbitals ( jstor )
Sine function ( jstor )
Spherical harmonics ( jstor )
Symmetry ( jstor )

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright the author. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Embargo Date:
5/4/2003
Resource Identifier:
70786468 ( OCLC )

Downloads

This item is only available as the following downloads:


Full Text

PAGE 1

DEVELOPMENT OF EFFICIENT ELECTR ON CORRELA TION METHODS F OR ONEAND TW O-DIMENSIONAL EXTENDED SYSTEMS AND THEIR APPLICA TIONS By MOTOI TOBIT A A DISSER T A TION PRESENTED TO THE GRADUA TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQUIREMENTS F OR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORID A 2002

PAGE 2

Cop yrigh t 2002 b y Motoi T obita

PAGE 3

A CKNO WLEDGMENTS First, I w ould lik e to thank Prof. Ro dney J. Bartlett for his sup ervision, patience, and time sp en t with me. He alw a ys brings me exciting pro jects and step b y step learning o ccasions. Prof. Sam uel B. T ric k ey has b een of great help for his kno wledge of densit y functional theory and p erio dic system researc h. Discussions ab out analytical deriv ativ e ev aluations of the p erio dic cell parameters w ere esp ecially useful for me to get a solution. Great appreciation go es to Prof. N. Yngv e Ohrn, Prof. Sam uel L. Colgate, and Prof. James D. Winefordner. I learned a great deal of c hemistry through their classes in group theory , thermo dynamics, and sp ectro c hemical analysis. Dr. So Hirata let me use his one-dimensional p olymer co de as I started m y o wn researc h for t w o-dimensional systems. He alw a ys ga v e me clear answ ers for m y questions ranging from science to tec hnical programming. My sp ecial gratitude go es to Dr. Ajith P erera, for his help. I alw a ys felt less tension after talking to him. Mr. Kenneth Wilson and Mr. An thon y Y au, t w o of m y colleagues in the Bartlett researc h group, ha v e b een a great help sho wing me ho w to use computer b etter. Finally , I express a w ord of thanks to m y family for letting me do what I w an t so freely , and m y anc ee, Mari T atsuza w a, for her men tal supp ort. iii

PAGE 4

TABLEOFCONTENTS page A CKN O WLEDGMENT S .. . .......... . .......... . ... . iii LIS T O F T ABLE S ...... . .......... . .......... . ... . vi LIS T O F FIGURE S ..... . .......... . .......... . ... . vii ABSTR A C T ......... . .......... . .......... . ... . ix CHAPTERS 1 GENERA L INT R ODUCTIO N ..... . .......... . ... . 1 2 TH E ANA L YTICA L ENE R G Y GRADIEN T SCHEM E I N THE G A USSIA N BASE D HA R TREEF OC K AN D DENSIT Y FUNCTIONA L THEO R Y F O R EXTENDE D SYSTEM S .... . ... . 5 2. 1 I n tr o ductio n . .......... . .......... . ... . 5 2. 2 Theor y ... . .......... . .......... . ... . 10 2.2. 1 Th e Crystallin e Orbita l Theor y fo r P eri o di c System s . 10 2.2. 2 Th e Multi p ol e Expansio n an d th e F as t Multi p ole Meth o d .......... . .......... . ... . 15 2.2. 3 Analytica l Gradie n t s o f th e T ota l Energ y p e r Uni t Cel l 23 2.2. 4 Gri d W eig h t Deri v ati v e s wit h Res p ec t t o Uni t Cell P arameter s ........ . .......... . ... . 25 2.2. 5 Th e MPE/FM M Energ y Deri v ati v e s ..... . ... . 29 2. 3 Result s ... . .......... . .......... . ... . 30 2. 4 Conclusio n . . .......... . .......... . ... . 36 3 F ORMA L ASPECT S O F A TOMIC-ORBI T A L BASE D COUPLEDCLUSTE R THEO R Y .......... . .......... . ... . 48 3. 1 I n tr o ductio n . .......... . .......... . ... . 48 3. 2 Co nv e n tiona l C C o v ervie w ... . .......... . ... . 50 3. 3 T o ol s Bridgin g Be tw ee n th e M O an d A O Represe n tation s . . 51 3.3. 1 Th e Orthogonali t y Conditio n ......... . ... . 51 3.3. 2 A O t o M O T ransformatio n .......... . ... . 53 3.3. 3 Beh a vio r o f th e A O Creation/annihilatio n O p erator s . 54 3.3. 4 Defnitio n o f th e Referenc e Stat e ....... . ... . 55 3. 4 Th e Groun d Stat e Couple d Cluste r Equation s ... . ... . 55 3. 5 W orkin g F or m ula s ........ . .......... . ... . 60 iv

PAGE 5

3.5. 1 Th e CCS D Groun d Stat e Equation s ..... . ... . 61 3.5. 2 Th e CCS D Gradie n t Equation s ....... . ... . 63 4 A TOMIC-ORBI T A L BASE D COUPLE D CLUSTE R THEO R Y F O R EXTENDE D SYSTEM S ..... . .......... . ... . 67 4. 1 I n tr o ductio n . .......... . .......... . ... . 67 4. 2 F ormalis m . . .......... . .......... . ... . 69 4. 3 Analysi s .. . .......... . .......... . ... . 77 4.3. 1 L o cali t y ......... . .......... . ... . 78 4.3. 2 Sparsi t y ......... . .......... . ... . 79 4.3. 3 Impleme n tatio n ..... . .......... . ... . 79 4.3. 4 Resourc e Manageme n t . . .......... . ... . 80 4. 4 Numerica l Result s ........ . .......... . ... . 85 4. 5 Summar y an d F utur e W or k ... . .......... . ... . 89 APPENDIX...................................94 A ABBREVI A TION S AN D SYMBOL S COMMON L Y USE D TH R OUGH TH E DISSE RTA TIO N ......... . .......... . ... . 94 B ANA L YTICA L GRADIENT S WIT H RESPEC T T O TH E THREE DIMENSIONA L UNI T CEL L P ARAMETER S ....... . ... . 97 C CONNECTIO N AMON G TH E NORMALIZE D SPHERICA L HARMONICS , MU L TIPOL E MOMENT S AN D SOLI D SPHERICAL MOMENT S .... . .......... . .......... . ... . 100 REFERENCE S ........ . .......... . .......... . ... . 103 BIOGRAPHICA L SKETC H . . .......... . .......... . ... . 108 v

PAGE 6

LIS T O F T ABLES T able page 2. 1 Co nv ergenc e o f th e tota l energ y p e r uni t cell , ban d gap , an d analytica l gradie n t s o f th e tota l energ y an d timin g fo r th e gradie n t s e v aluatio n wit h res p ec t t o th e SC F co nv ergenc e criterion , S C F e v aluated a t BN b on d length s o f 1.5 0 A an d 1.3 0 A . Th e tw o-electro n i n tegral cuto parameter , I N T =1 4 an d long-rang e eec t cuto parameter f = 5 ar e fxed . Th e tota l energie s an d gradie n t s ar e gi v e n i n atomic units , ban d gap s ar e i n eV . ........ . .......... . ... . 37 2. 2 Co nv ergenc e o f th e tota l energ y p e r uni t cell , ban d gap , an d analytica l gradie n t s o f th e tota l energ y an d timin g fo r th e gradie n t s e v aluatio n wit h res p ec t t o th e tw o-electro n i n tegra l cuto parameter , I N T e v aluate d a t BN b on d length s o f 1.5 0 A an d 1.3 0 A . Th e SC F conv ergenc e criterio n S C F = 8 an d long-rang e eec t cuto paramete r f =5 ar e fxed . Th e tota l energie s an d gradie n t s ar e gi v e n i n atomi c units, ban d gap s ar e i n eV , an d timing s ar e i n seconds . ....... . ... . 38 2. 3 Co nv ergenc e o f th e tota l energ y p e r uni t cell , ban d gap , an d analytica l gradie n t s o f th e tota l energ y wit h res p ec t t o th e long-rang e Coulo m b eec t cuto paramete r f e v aluate d a t BN b on d length s o f 1.5 0 A an d 1.3 0 A . Th e SC F co nv ergenc e criterio n S C F = 8 an d th e tw o-electron i n tegra l cuto parameter , I N T =1 4 ar e fxed . Th e tota l energie s and gradie n t s ar e gi v e n i n atomi c units , ban d gap s ar e i n eV . ... . ... . 39 2. 4 Th e optimize d geometr y , ban d ga p an d Mulli k e n c harg e o f hexagonal b oro n nitrid e ..... . .......... . .......... . ... . 39 4. 1 Th e n u m b e r o f double s excitatio n amplitude s t o determin e fo r p olyace t ylen e calculation s i n A O-C C an d co nv e n tiona l CC . C O base d amplitude s ar e determine d base d o n 2 0 kp oi n ts . Uni t : millions . . ... . 91 4. 2 Th e MBPT(2)/STO-3 G energ y fo r th e p ol y ace t ylen e c hai n calculated b y th e co nv e n tiona l meth o d a s a functio n o f th e lattic e su m cuto parameter . ...... . .......... . .......... . ... . 91 4. 3 Th e A O-MBPT(2)/STO-3 G energ y p e r uni t cel l fo r th e p ol y ace t ylen e c hai n calculate d a s a functio n o f th e i n tegral s an d amplitud e threshol d ma x= 1 0)Tj/R35 1 Tf 0.04625 Tc 0.828 0 Td(A. Th e cel l summatio n truncatio n an d k p oi n t s fo r the H F par t ar e fxe d a t 12 . EH F= )Tj/R6 1 Tf 0.772 0 Td(75 : 9443069 Eh....... . ... . 91 vi

PAGE 7

LISTOFFIGURES Figure page 2. 1 Th e defnitio n o f th e uni t cel l parameter s a , b , an d r an d th e cel l indice s qxan d qyfo ra genera l tw o-dimensiona l p eri o di c system . ... . 40 2. 2 Th e regio n wher e th e long-rang e Coulo m b i n teraction s ar e computed b y m ulti p ol e expansion . Th e referenc e cel l i s l o cate d a t th e ce n ter. I n thi s particula r fgur e UX= 13 , UY= 14 , TX= 4 an d TY= 5 are c hosen . (a ) th e i n teractio n b e tw ee n th e referenc e cel l an d ea c h cell i n th e shade d regio n i s compute d b y th e MPE . (b ) F o r th e F o c k matri x correctio n an d th e tota l energ y correctio n calculation , th e total co n tributio n i s exactl y th e t wic e th e co n tributio n fro m th e shaded regio n R . ....... . .......... . .......... . ... . 41 2. 3 Illustratio n o f th e FMM . Th e referenc e cel l i s a t th e lef t b otto m corner . I n th eE region , th e Coulo m b i n teraction s ar e calculate d b y the tw o-electro n i n tegrals , whil e a t th e F(i ) (i=0,1,2,.. ) th e Coulo m b interaction s ar e appr o ximate d b y FMM . Th e F(0 ) regio n i s equi v ale n t t o th e MP E wit h m ulti p ole s b ein g defne d fo r cell s o f th e original size , wherea s dista n t cell s ar e regrou p e d i n t o large r cell s a s i n F(1) an d F(2) . ....... . .......... . .......... . ... . 42 2. 4 Th e fou r step s i n constructin ga ne w se t o f m ulti p ole s defne d fo r a consolidate d 3 3 cell . . .......... . .......... . ... . 43 2. 5 Th e pr o cedur e fo r axi s rotatio n t o ma k e th e tw o c o ordinat e systems ce n tere d o n th e x y -plan e h av e a commo n z -axi s ....... . ... . 44 2. 6 Geometr y o f m o de l h ydroge n ruorid e c hai n sta c k . a =4.849 7 A , b =4.0 A , an d r = 9 0 degrees . .......... . .......... . ... . 45 2. 7 A compariso n o f timing s amon g th e m ulti p ol e expansion s o f v arious order s an d th e explici t treatme n t o f long-rang e par t a s a functio n of n u m b e r o f cell s include d i n th e calculation . Timing s ar e gi v e n i n seconds . ......... . .......... . .......... . ... . 46 2. 8 A compariso n o f timing s amon g th e fas t m ulti p ol e meth o d o f v arious order . Timing s ar e gi v e n i n seconds . ... . .......... . ... . 46 vii

PAGE 8

2. 9 L a y ere d structur e fo ra Boron-nitrid e calculation . Bla c k cel l i s the referenc e uni t cell , th e dar k gr a y cell s ar e treate d b y th e explici t calculatio n o f tw o-electro n i n tegrals , th e lighe r gr a y cell s ar e treate d b y th e m ulti p ol e expansio n an d th e lig h tes t gr a y cell s ar e treate d b y the fas t m ultipl e meth o d . . .......... . .......... . ... . 47 2.1 0 Thre e symmetricall y equi v ale n t bu t diere n tl y treate d BN b onds (bro k e n lines) . .... . .......... . .......... . ... . 47 4. 1 Th e dec a y b ehai v o r o f th e maxi m u m v alue s o f th e regula r tw o-electron i n tegral s an d A O-C C transforme d i n tegrals . ......... . ... . 92 4. 2 Timin g p e r iteratio n fo r th e Li H c hai n A O-MBPT(2 ) calculations wit h v ariou s cuto parameter s an d basi s sets . ........ . ... . 93 4. 3 Timin g p e r iteratio n fo r th e Li H c hai n A O-LCC D calculation s with v ariou s cuto parameter s an d basi s sets . . .......... . ... . 93 viii

PAGE 9

Abstract of Dissertation Presen ted to the Graduate Sc ho ol of the Univ ersit y of Florida in P artial F ulllmen t of the Requiremen ts for the Degree of Do ctor of Philosoph y DEVELOPMENT OF EFFICIENT ELECTR ON CORRELA TION METHODS F OR ONEAND TW O-DIMENSIONAL EXTENDED SYSTEMS AND THEIR APPLICA TIONS By Motoi T obita Ma y 2002 Chair: Ro dney J. Bartlett Ma jor Departmen t: Chemistry This dissertation is fo cused on the dev elopmen t of highly accurate electroncorrelation metho ds for oneand t w o-dimensional p erio dic systems. F or onedimensional systems, atomic-orbital based man y-b o dy p erturbation theory and coupled cluster theory are dev elop ed and applied to p oly acet ylene and lithiumh ydride mo del c hain. Use of atomic orbitals instead of con v en tional crystalline (molecular) orbitals enables us to con trol run time and accuracy of the calculation. The gain in the atomic-orbital based coupled cluster metho d originates from the lo calit y and sparsit y of matrices dened in the framew ork of the theory . Pro vided, ecien t sparse matrix-matrix m ultiplication routines, w e obtain go o d estimates of the correlation energy m uc h faster than the con v en tional metho d for large systems. 90 p ercen t of the correlation energy can b e reco v ered v ery quic kly , and 2 to 3 digit accuracy can b e obtained for p olymers with relativ ely simple ix

PAGE 10

unit cells suc h as p oly acet ylene. The formalism b ehind the atomic-orbital based coupled cluster theory is applicable for b oth large molecular systems and p erio dic systems. F ormal asp ects of the atomic-orbital based coupled cluster theory are discussed and corresp ondence b et w een the atomic-orbital based framew ork and crystalline (molecular) orbital based framew ork are sho wn. Tw o-dimensional co de dev elopmen t is based on Hartree-F o c k and densit y functional theories due to the fact that coupled-cluster theory is to o costly . Ecien t inclusion of the Coulom b eects b y the fast m ultip ole metho d and analytical gradien t tec hniques are the core elemen ts that con tribute robustness and computational eciency for t w o-dimensional systems. The fast m ultip ole metho d is an algorithm to include the long-range Coulom b eects for uniform systems with linear-scaling costs for molecular systems and with logarithmic scaling for innite p erio dic systems. The analytical gradien t tec hnique is a p o w erful to ol when optimized geometries or vibrational frequencies are computed. If optim um geometries or vibrational frequencies are required, then analytical gradien ts are for all practical purp ose, a necessit y . x

PAGE 11

CHAPTER 1 GENERAL INTR ODUCTION The end of the 20th cen tury and the b eginning of the 21st cen tury ha v e witnessed the great success of quan tum c hemistry as a practical and cost eectiv e to ol for molecular and materials mo deling. Tw o Nob el prize winners in c hemistry in 1998 from the quan tum c hemistry comm unit y are the solid pro of that quan tum c hemical metho ds are widely accepted and recognized as p o w erful to ols whic h can b enet science and tec hnology comm unities. One of the winners is a ph ysicist who largely con tributed to constructing the fundamen tals of densit y functional theory (DFT). The other is a c hemist whose w ork are fo cused on quan titativ e calculations of relativ ely small molecules using w a v e-function based metho ds. One of the recen t eorts of the quan tum c hemistry comm unit y is to mak e quan tum c hemical to ols, compared to other semi-empirical or molecular mec hanical mo dels, able to treat larger systems. The capabilit y will help us study imp ortan t questions of c hemistry , particularly those related to biology and materials science where the systems of in terest t ypically in v olv e a large or ev en an innite n um b er of atoms. Man y of the surface phenomena, p olymerization pro cesses, c hemical reactions under the presence of catalysts, large molecular clusters, and molecules of biological in terests are all in this category . As sp ecial cases of large systems, there are innite p erio dic systems suc h as p olymers, slabs, and crystals, though p erio dic systems ha v e unique ph ysical features, whic h are not dened in non-p erio dic large systems. Correlated treatmen ts of these p erio dic systems w ere started earlier b y ph ysicists using DFT. Esp ecially , band structure calculations on solids ha v e b een p opular. Also DFT succeeded in describing, qualitativ ely , man y prop erties suc h 1

PAGE 12

2 as stabilities and phase transitions under nite pressure. While w a v e-function metho ds could treat only small systems, they pro vided more quan titativ e answ er. Moreo v er, geometry optimization tec hniques are adv anced with the help of the analytical ev aluation of the energy gradien ts. Also, particular forms of exc hangecorrelation functionals dev elopp ed b y c hemists enhanced the predictabilit y of DFT for molecules. Applications of the metho ds to p erio dic systems ha v e b een dicult o wing to their demanding computational requiremen ts, th us there is a shorter history con tributing in this area from c hemists. T aking these demands and diculties in to accoun t, I dedicate m y eorts to t w o main directions. The rst is the quan tum c hemical treatmen t of t w odimensional p erio dic systems, whic h is discussed in c hapter 2. The second is the highly accurate correlated treatmen t of one-dimensional p erio dic systems. F ormal discussion of the atomic-orbital based coupled cluster theory is discussed in c hapter 3, and applications to innite p olymers are giv en in c hapter 4. The goal in these w ork is to compute quan tities p er unit cell suc h as the total energy and energy gradien ts for eac h atom in the reference cell. In b oth cases, the dicult y lies in ho w to implemen t the metho ds whic h w arran t ecien t executions. Man y parts of the form ulations are similar to their molecular coun terparts, ho w ev er, most of their implemen tations are targeted to w ard small molecules and th us straigh tforw ard applications of suc h algorithms for larger systems, suc h as p olymers, are computationally prohibitiv e. As a result, parts of the form ulations are new and non-trivial, implemen tation is v ery c hallenging, and b oth are in terrelated. In the quan tum-c hemical treatmen t of t w o-dimensional systems, the most c hallenging part of dev eloping an ecien t computer program is the inclusion of the long-range Coulom b eect. Since the explicit ev aluation of t w o-electron in tegrals is computationally in tensiv e, suc h calculations m ust b e conned only within the

PAGE 13

3 sev eral neigh b or cells around the reference unit cell. In the presen t w ork, w e c ho ose that the remaining part b e treated b y the m ultip ole expansion tec hnique and its extension, the fast m ultip ole metho d. The computational cost for the m ultip ole expansion tec hnique for t w o-dimensional systems scales O ( D 2 ) where D is the cuto radius, while the fast m ultip ole metho d enables us to further reduce the scaling to O (log D 2 ). By virtue of the m ultip ole expansion tec hnique, w e can accoun t for a considerable p ortion of the long-range eects ranging through a few micrometers, while w e include the Coulom b in teraction in the region in the length scale of a few millimeters with the fast m ultip ole metho d. Since our reference cell size is in the order of a nanometer, it is reasonable to sa y that it is p ossible to obtain the iden tical n umerical precisions of desired quan tities as the quan tities for true innite systems, via the fast m ultip ole metho d. Another p ortion of m y eorts on t w o-dimensional systems is the implementation of the analytical gradien t of the total energy p er unit cell with resp ect to geometrical parameters including lattice constan ts. F or the N-atom p er cell systems, w e need to p erform 6N+6 times single p oin t calculations (plus and min us displacemen ts for x , y , and z directions for eac h atom, and the three lattice constan ts) to obtain gradien ts for all the atoms and cell parameters for systems with no symmetries. This forces us to study only highly symmetric systems where the n um b er of single p oin t calculations can b e signican tly reduced. The same task is ecien tly done with the analytical gradien t tec hnique, esp ecially when the system con tains a large n um b er of atoms. Again, the analytical ev aluation of the long-range Coulom b in teraction energy deriv ativ es is highly non-trivial. In this w ork, w e presen t the formalism and implemen tation of the analytical gradien ts of the m ultip ole expansion and fast m ultip ole metho d for p erio dic systems using Gaussian basis sets for the rst time.

PAGE 14

4 In the second part, w e will describ e the second-order man y-b o dy p erturbation theory [MBPT(2)] and coupled cluster theory doubles (CCD) form ulated and implemen ted for one-dimensional innite p erio dic systems. The steep cost of the MBPT(2)[ O ( N 5 )] and CCD[ O ( N 6 )] mo dels, with N b eing the n um b er of basis functions in a unit cell, had prev en ted us from p erforming these calculations for p olymers, except for those with v ery small unit cells or basis sets. The problem is o v ercome b y reform ulating the man y-b o dy p erturbation theory and coupled cluster theory in atomic-orbital basis. T aking adv an tage of the lo calit y of atomic orbitals, w e can reduce the scaling of these metho ds to O ( N ) asymptotically in the limit of N ! 1 . This enables us to p erform one of the most accurate calculations for innite systems.

PAGE 15

CHAPTER 2 THE ANAL YTICAL ENER GY GRADIENT SCHEME IN THE GA USSIAN BASED HAR TREE-F OCK AND DENSITY FUNCTIONAL THEOR Y F OR EXTENDED SYSTEMS 2.1 In tro duction One of the in terests in quan tum c hemistry is to pro vide equilibrium geometries and vibrational sp ectra of systems. These prop erties require the rst and second total energy deriv ativ es with resp ect to displacemen ts of the atoms. Therefore, the a v ailabilit y of analytical expressions for the total energy deriv ativ es is of great imp ortance. With analytical rst deriv ativ es, not only is computation time reduced but also n umerical noise is reduced whic h could b e serious in n umerically constructing secondor higher-deriv ativ es. The Hartree-F o c k (HF) analytical energy gradien t and Hessian tec hniques using Gaussian basis functions w ere originally dev elop ed for molecular systems. Deriv ations for these expressions w ere dev elop ed b y sev eral authors [1{3] b et w een the late 1950s and the early 1970s. The HF energy gradien t p er unit cell for one-dimensional (1D) p olymers app eared in 1983 [4, 5]. The analytical gradien t for 1D extended systems using Gaussians w as implemen ted in 1997 [6] and sev eral applications [6{8] ha v e b een p erformed successfully to obtain optimized geometries and vibrational frequencies of p olymers suc h as p oly acet ylene, p olymethineimine, p oly eth ylene, and p olydiacet ylene. Judging from these successes, dev elopmen t of the analytical energy gradien t sc heme for t w oand three-dimensional (2D,3D) p erio dic systems is exp ected to yield signican t computational cost sa vings. Dev elopmen ts in 2D systems in early da ys w ere exclusiv ely done in densit y functional theory (DFT), where band calculations 5

PAGE 16

6 w ere most p opular and energetics follo w ed. Among these, A Gaussian based co de w as rst made b y Calla w a y et al. [9]. Since then there has b een extensiv e amoun t of w ork and sev eral co des p erforming 2D systems with Gaussians exist [10{ 14]. Analytical gradien ts are not implemen ted in an y of ab o v e except the recen t preliminary w ork b y Doll [15] where HF gradien ts for 2D and 3D systems ha v e b een implemen ted excluding those with resp ect to the lattice constan ts. Th us from the p ersp ectiv e of the presen t w ork, the form ulas to calculate analytical gradien ts with resp ect to lattice constan ts and those in the DFT framew ork in higher dimensions are the new con tributions. Extension of the metho d for 1D systems to 2D and 3D in tro duces t w o issues. The rst is the optimization of the unit cell parameters, or the lattice constan ts. In the 1D case w e need to optimize b oth the translational unit length and the atomic p ositions in order to p erform a full geometry optimization, but the expression to ev aluate the deriv ativ e with resp ect to the former [4, 5] is simple. Ho w ev er, in the 2D and 3D cases, w e also need to consider the unit cell angles as w ell as the translational unit lengths, whic h complicates the formalism. Optimizations of the unit cell parameters w ere originally view ed as minimization of stress and force on a unit cell of solids [16, 17]. Nielsen and Martin considered innitesimal p osition c hanges in orbitals and deriv ed an exc hange-correlation functional resp onsible for the isotropic stress. The expression is extended for the surface stress using plane w a v e basis [18], and v ery recen tly the analytic stress tensor for p erio dic systems w as presen ted b y Kudin and Scuseria [19]. If there is no external stress on a system, the optimized geometry should corresp ond to a state with no stress. Th us, the stress tensor w as used for the optimization with resp ect to the unit cell parameters. W e c ho ose to use an alternativ e approac h, namely the direct computation of the energy deriv ativ es with resp ect to the crystallographic lattice constan ts a , b , c , , and r . W e start the deriv ation in terms of in tegrals formed b y atomic

PAGE 17

7 orbitals. Since a crystalline orbital is expanded b y a linear com bination of atomic orbitals, the w orking equations are directly obtained b y considering the c hange in orbitals, in turn, c hange in in tegrals formed b y atomic orbitals. Considering innitesimal c hanges of atomic p ositions, and hence the deriv ativ e in tegrals, yields w orking form ulas for the ev aluation of the analytical gradien t. The gradien ts with resp ect to atomic displacemen ts and to the unit cell parameters enable us to p erform a full analytical geometry optimization. It should b e noted that the analytic ev aluation of the energy gradien ts with resp ect to unit cell parameters do es not require an y additional computation of in tegrals but needs only the in tegral deriv ativ es with resp ect to atomic displacemen ts. Th us, the optimization of the unit cell size and that of the atomic p ositions within the unit cell can b e p erformed at the same time. The second issue that arises when going from 1D to higher dimensions is the treatmen t of the long-range in teraction. The magnitude of the Coulom b in teraction deca ys as 1 =r with r b eing the distance b et w een t w o particles or c harge distributions. Recalling that w e are trying to describ e innite systems, the Coulom b in teractions m ust b e prop erly included up to innit y . Due to the slo wly deca ying b eha vior of the Coulom b in teraction, a truncation of the long-range Coulom b eect often leads to an inaccurate description of innite systems esp ecially when the unit cell has a p ermanen t dip ole. The long-range Coulom b in teraction can b e ecien tly included using the m ultip ole expansion (MPE) tec hnique, pro vided a sophisticated lattice summation tec hnique. A go o d review of the tec hnique is found in Hirsc hfelder et al. [20]. The lattice summation using MPE orignated from Ew ald [21], and extension to electronic strucure w as made b y L owdin [22] where the idea of com bining n uclear c harges and electronic c harges to mak e the total c harge neutral, whic h is required to obtain con v ergence, is in tro duced. Delhalle and co-w ork ers presen ted all the necessary

PAGE 18

8 form ulas to compute the m ultip ole correction for 1D systems [23{ 25] using Cartesian m ultip oles expanded b y Gaussian functions within the HF framew ork. An adv an tage of the Cartesian m ultip oles is their ease with whic h the momen t in tegrals can b e co ded when the maxim um quan tum n um b er l max is small suc h as the dip ole ( l = 1) and quadrup ole ( l = 2). Ho w ev er, as l increases, form ulating the MPE b ecomes exp onen tially complicated. Jacquemin et al. [26] form ulated inclusion of the long-range eect for the total energy p er unit cell and its gradien ts based on a m ultiple T a ylor expansion tec hnique. Using a McMurc hie-Da vidson t yp e in tegral ev aluation sc heme [27], the long-range terms are added in to the n uclear-repulsion energy , o v erlap and kinetic energy in tegrals, electron-n uclear attraction and t w o-electron in tegrals. Alternativ ely , w e c ho ose to use spherical m ultip oles in our m ultip ole expansion. All the n uclear momen ts and electronic momen t in tegrals are ev aluated recursiv ely [28], making it p ossible to easily include large v alues of l 's. F ormally , the cost scaling of the MPE is O ( D n ), where D is the cuto distance within whic h the Coulom b in teractions are ev aluated b y MPE and n is the dimension of the system. In 1D systems, b ecause all the cells are linearly aligned, the summation o v er all the cells is sho wn to b e the Riemann-zeta function. This reduces the cost scaling of the 1D MPE from O ( D ) to O (1) b y computing the Riemann zeta function in adv ance. In higher dimensions, w e are unable to use the Riemann-zeta function as the innite summation tec hnique. Th us, a straigh tforw ard application of MPE scales nominal O ( D 2 ) and O ( D 3 ). Apparen tly w e need to c ho ose an innite summation tec hnique to rapidly include the long-range eects. Tw o ma jor approac hes are Ew ald summation tec hnique [21] and its v arian ts and fast m ultip ole metho d [29{ 32] (FMM). Detailed theoretical and n umerical comparisons of these can b e found elsewhere [33], but summarizing in short, b oth are capable of including the long-range con tributions rapidly . One of their dierences is that the Ew ald based approac hes ev aluate the truely

PAGE 19

9 innite eects b y ev aluating the long-range eects in k -space, while the FMM in tro duces a cuto in real space. Though the cuto distance can b e c hosen arbitrary large and th us n umerical results w e obtain w ould not b e aected, the metho d lac ks the mathmatical rigor to include the truely innite eects. On the other hand, for the same reason, FMM has a strong connection to treat large but nite systems. A generalized Ew ald tec hnique for 2D p erio dic systems w as rstly implemen ted b y Birk enheuer et al. [34] whic h w as an exten tion of w ork b y Min tmire et al. [35] within the lo cal densit y appro ximation with a tting of the Exc hange p oten tial and b y Bo ettger and T ric k ey [36]. The fast m ultip ole metho d w as initially dev elop ed b y Greengrad and Rokhlin [29] for systems with p oin t c harge, and exten tion to p erio dic systems follo w ed [30]. The metho d w as in tro duced to electronic structure theory b y White et al. [31] and also extended to p erio dic systems [32]. W e c ho ose to use the fast m ultip ole metho d due to the strong connection with large molecular systems. The idea of the FMM is that w e group m ultip oles of adjacen t distan t cells together as the distance D b ecomes larger to reduce computational cost. This grouping can b e done recursiv ely and reduces the scaling to O (log D n ). In molecular cases, w e need to consider all the pairs of cells to compute the Coulom b eects. The maxim um FMM algorithm can ac hiev e is linear scaling. Ho w ev er, in extended systems, one cell of the pair is alw a ys the reference cell, th us the cost of ev aluating the Coulom b eects decreases as the distance b et w een cells increases, ac hieving logarithmic scaling. Using FMM, a calculation is carried out in a m ulti-la y ered structure: the sev eral nearest neigh b oring cells around the reference cell are treated explicitly using the regular oneand t w o-electron in tegrals, and the cells just b ey ond the nearest neigh b ors are treated b y MPE. Outside of the MPE region, w e p erform FMM recursiv ely . The b oundary b et w een the explicit and MPE region can b e c hosen based on the cuto distance, and th us w e ha v e con trol o v er accuracy . W e c ho ose the sizes of the

PAGE 20

10 in teracting cells to b e 3 n , 9 n times that of the reference cell. In this sc heme w e could include the long-range Coulom b eects up to microto milli-meter region around the reference cell. This w ould giv e almost the same n umerical result to a true innite system calculation. The total energy and analytical gradien ts for 2D HF and DFT are implemen ted in to the quan tum c hemistry co de pol ymer2.0 [37]. An o v erview of the crystalline orbital theory is giv en in the next section follo w ed b y the analytical gradien t expressions with an in-depth discussion of the long-range in teraction treatmen t and the problems sp ecic to p erio dic DFT. The gain in CPU time to accurately include the long-range eects b y MPE and FMM is signican t, whic h is n umerically sho wn b y a h ydrogen-ruoride c hain 2D stac k mo del system. Also, the p erformance of the analytical gradien t and its con v ergence prop erties with resp ect to v arious adjustable parameters is in v estigated. The in v estigation and geometry optimizations are p erformed on a hexagonal b oron-nitride innite sheet. In the app endices, mathematical details of the FMM algorithm and the form ulas required for 3D extension are presen ted. 2.2 Theory 2.2.1 The Crystalline Orbital Theory for P erio dic Systems The form ulas for computing the analytical total energy gradien ts within the framew ork of Gaussian-function-based Hartree-F o c k and densit y functional theories are presen ted in this section. As a prerequisite for the energy gradien t deriv ation, the total energy expression is briery review ed. Detailed deriv ations ha v e b een giv en elsewhere [1{6]. Some emphasis is placed on 2D systems, but most parts of the discussion are for innite p erio dic systems in general. The total energy p er

PAGE 21

11 unit cell in p erio dic systems is expressed as E = X ; X q P ( q ) H ( q ) + 1 2 X ; X ; X q ; r ; s P ( q ) P ( s r ) ( ( 0 ) ( q ) j ( r ) ( s ) ) + m 1 Z cell f [ ; r ] d r m 2 4 X ; X ; X q ; r ; s P ( q ) P ( s r ) ( ( 0 ) ( r ) j ( q ) ( s ) ) + E NR + ( E LR ) ; (2.1) where ; ; and are atomic orbital lab els, and q ; r ; and s are cell v ector lab els. The p osition of the cell q is giv en as q = q x a + q y b : (2.2) with a and b b eing the t w o lattice v ectors. The dierence b et w een the molecular orbital theory and the crystalline orbital theory lies in summing o v er the unit cells. The dimension of the v ector is iden tical to the dimension of the p erio dicit y , i.e. , 0 for molecules, 1 for p olymers, 2 for slabs, and 3 for crystals. This is illustrated in Figure 2.1 for a 2D system. The reference unit cell is denoted as ( q x , q y )=(0,0)= 0 . m 1 and m 2 are mixing co ecien ts of the exc hange-correlation energy and the exact single-determinan t exc hange energy . f is an exc hange-correlation functional in a general form. E NR denotes the n uclear-n uclear repulsion energy p er unit cell whic h is the sum of the classical Coulom b energies generated b y n uclear c harges. E LR is in tro duced as a result of truncating the lattice summation and can b e ev aluated b y MPE [23{ 25] or FMM [29], whic h will b e discussed in detail later. The one-electron in tegrals in Eq.(2.1) are dened as H ( q ) = Z ( 0 ) ( r ) 1 2 r 2 ( q ) ( r ) d r X A X s Z ( 0 ) ( r ) Z A j r R ( s ) A j ( q ) ( r ) d r ; (2.3)

PAGE 22

12 where Z A is the c harge of a n ucleus A at p osition R ( s ) A and ( q ) is a Gaussian t yp e atomic basis function lo cated at cell q . The elemen ts of the densit y matrix P ( q ) are dened as P ( q ) = 2 K occ: X j B Z X k c [ k ] j c [ k ] j exp ( i k q ) ; (2.4) where K is the n um b er of Brillouin zone sampling p oin ts. The sampling in 2D is done with an ev enly spaced mesh in the, in general, paralellogram Brillouin zone. The t w o-electron in tegrals are dened as ( ( 0 ) ( q ) j ( r ) ( s ) ) = Z Z ( 0 ) ( r 1 ) ( q ) ( r 1 ) 1 r 12 ( r ) ( r 2 ) ( s ) ( r 2 ) d r 1 d r 2 ; (2.5) where the 1 =r 12 is the distance b et w een the cen ters of c harge distributions made b y electron 1 and 2. Ultimately w e ha v e the k -dep enden t self-consisten t eld equations to solv e X ( F [ k ] " [ k ] i S [ k ] ) c [ k ] i = 0 ; (2.6) where the F o c k and o v erlap matrix elemen ts in realand recipro cal spaces are giv en b y F [ k ] = X q F ( q ) exp ( i k q ) ; (2.7) S [ k ] = X q S ( q ) exp ( i k q ) ; (2.8) F ( q ) = H ( q ) + X s ; r X P ( s r ) [( ( 0 ) ( q ) j ( r ) ( s ) ) m 2 2 ( ( 0 ) ( r ) j ( q ) ( s ) )] + m 1 X ( q ) ; (2.9) S ( q ) = Z ( 0 ) ( r ) ( q ) ( r ) d r : (2.10) The exc hange-correlation energy , in general, is expressed as a functional of the densit y and the densit y gradien ts. The exc hange-correlation in tegral matrix

PAGE 23

13 elemen t X ( q ) is giv en [38, 39] b y X ( q ) = Z ( 0 ) ( r ) v X C ( q ) ( r ) + r ( 0 ) ( r ) g X C ( q ) ( r ) + ( 0 ) ( r ) g X C r ( q ) ( r ) d r ; (2.11) with v X C = @ f @ ; (2.12) g X C = @ f @ r + 1 2 @ f @ r r ; (2.13) r = r r ; (2.14) r = r r ; (2.15) = 1 2 X q X P ( q ) ( 0 ) ( q ) ; (2.16) = + ; (2.17) where and denote spin, and and r are the densit y and densit y gradien t in v arian ts. As men tioned, a unique feature of the p erio dic system form ulations in con trast to the molecular coun terpart is the presence of the lattice sums i.e. , the summation o v er q in Eq.(2.1). In principle, the summation m ust extend to innit y for accurate results. Ho w ev er, in practice w e truncate the innite lattice summations with resp ect to q , r and s after nite terms. In this w ork, w e emplo y the so-called \Nam ur cuto criteria" [24, 37] for this purp ose. This criterion in tro duces t w o parameters T and U , and the t w o-electron in tegrals are ev aluated only when the indices on the in tegrals satisfy the follo wing conditions: T i q i T i ; (2.18) U i r i U i ; (2.19) r i T i s i r i + T i ; (2.20)

PAGE 24

14 U i s i q i U i (2.21) i = x; y ; z : (2.22) It is not costly to obtain con v erged results with resp ect to the T parameter, since the t w o-electron in tegrals deca y exp onen tially with increasing the distance b et w een the atomic orbitals ( 0 ) and ( q ) and that b et w een ( r ) and ( s ) . In con trast, the rate of con v ergence of the total energy with resp ect to U is m uc h lo w er since the t w o-electron in tegrals ha v e asymptotic 1 =r b eha vior with r b eing the distance b et w een the cen ters of the t w o c harge distributions. This b eha vior also applies to the n uclear-electron attraction in tegrals and the n uclear-n uclear repulsion energy term. This slo w con v ergence problem has b een successfully o v ercome using MPE [23{ 25] or FMM [29]. The E LR term in Eq.(2.1) is in tro duced to accoun t for the trancated lattice sum. Tw o-dimensional MPE and FMM are discussed in depth and 3D extension is briery addressed in the follo wing subsection. W e emplo y the recursion relation of Obara and Saik a [40] for the computation of the t w o-electron in tegrals and the in tegral deriv ativ es. This requires explicit ev aluation of only ( ss j ss ) t yp e in tegrals to generate in tegrals of higher angular momen tum pro vided the atomic orbital cen ter and the Gaussian exp onen t are iden tical. Since the storage requiremen ts for the in tegrals exceed sev eral gigab ytes with a mo dest basis set in calculations of p erio dic systems with no symmetry , most suc h calculations ha v e to b e p erformed b y in tegral-direct metho ds [41, 42]. In a direct algorithm, the t w o-electron in tegrals are calculated at eac h SCF iteration, th us completely a v oiding the storage. Ev aluation of `nearly zero' in tegrals is a v oided b y a prescreening tec hnique [42]. Only when j A B j 10 I N T (2.23)

PAGE 25

15 is satised, do es computation of ( j ) pro ceed, where A = j ( j ) j 1 2 (2.24) B is similarly dened using and . I N T is a user-dened in tegral cuto parameter. 2.2.2 The Multip ole Expansion and the F ast Multip ole Metho d As discussed in the previous subsection, the innite summation o v er the cell indices is the ma jor limitation in CO calculations. w e use the FMM to rapidly ev aluate the long range Coulom b in teractions. The accuracy of the FMM dep ends on three factors. The rst is the order of the m ultip ole expansion truncation l max . F or l max = 1 , w e w ould obtain the exactly the same answ er as the exact result for a giv en in teraction. The second is c hoice of the b oundaries dening the explicit, MPE and FMM regions. And nally , the ph ysical cuto of the FMM regions. Note that all of the three are adjustable making the calculations to ha v e arbitrary precision. W e compute the long-range in teractions as the sum of the in teractions b et w een c harge distributions in the reference cell and in all the distan t cells. The m -th comp onen t of the 2 l th electronic momen ts expressed in the spherical co ordinates whic h are asso ciated with the orbital pro duct, all n uclei or the total c harge as N ( q )( l ;m ) , N ( l ;m ) n uc and N ( l ;m ) tot resp ectiv ely , are dened as N ( q )( l ;m ) = Z ( 0 ) ( r ) r l P m l (cos ) e im ( l + j m j )! ( q ) ( r ) d r h j N ( l ;m ) j i ; (2.25) N ( l ;m ) n uc = # atoms X A =1 Z A r l A P m l (cos A ) e im A ( l + j m j )! # atoms X A =1 Z A N ( l ;m ) A ; (2.26) N ( l ;m ) ele = X q X P ( q ) N ( q )( l ;m ) ; N ( l ;m ) tot = N ( l ;m ) n uc + N ( l ;m ) ele ; (2.27)

PAGE 26

16 with X q = T X q x = T T X q y = T ; (2.28) and P m l is the asso ciated Legendre function. First, let us consider the Coulom b energy p er unit cell con tributed from one distan t cell. T aking adv an tage of the translational symmetry in p erio dic systems, a c harge distribution at the reference cell and one at a distan t cell are iden tical. Consequen tly , the Coulom b in teraction energy p er unit cell con tributed from a distance cell, E (t w o ), is giv en as E (t w o) = X l ;l 0 ;m C ( l ; l 0 ; m ) N ( l ;m ) tot N ( l 0 ;m ) tot r l + l 0 +1 ; (2.29) where r is the distance b et w een the cen ters of the t w o c harge distributions and C ( l ; l 0 ; m ) is a constan t factor. The form ula is sligh tly dieren t from that usually used in molecular linear scaling metho ds using m ultip ole expansion tec hniques [43{45], where the denominator is incorp orated in one of the m ultip oles, as w ell as a c hange in the prefactor, yielding summations tak en only o v er l and m . Instead, t w o dieren t m ultip oles ha v e to b e ev aluated as a trade-o. Owing to the translational symmetry in p erio dic systems, w e nd that the ab o v e form ula has more computational adv an tage for calculating the long-range in teractions for extended systems. The form ula (2.29) assumes that t w o m ultip oles are dened b y a common z -axis and parallel x and y -axes. F or 1D systems, this condition is easily satised b y c ho osing the z -axis and the translational v ector to b e parallel. Ho w ev er, in 2D and 3D systems, the m ultip oles m ust b e rotated so that they are dened at the co ordinate system with the z -axis in common. This requiremen t is a trade-o to increase the eciency of the algorithm. The gain is that w e compute the n uclear and electronic momen t in tegrals only once taking adv an tage of the translational symmetry of the p erio dic systems. The in teractions b et w een m ultip oles at the

PAGE 27

17 reference cell and a distance cell are dieren t and can b e ev aluated b y computing momen t in tegrals and constructing m ultip oles for eac h cell. Instead, computing only one set of m ultip oles and rotating them as needed yields the same result whic h con tributes to eciency . A similar idea w as used elsewhere for a 3D crystal case [46]. Th us, ev aluating Eq.(2.29) generally requires rotations of the m ultip oles b eforehand. Denoting the prop erly rotated m ultip oles b y an o v erbar, w e obtain the total long-range Coulom b energy on the reference cell caused b y all the other distan t cells, E LR , are obtained b y simply summing o v er all the in teraction energies b et w een the reference cell and a distance cell, E(t w o). E MPE LR = X q E (t w o) ( q x ; q y ) = X q 1 X l =0 1 X l 0 =0 min ( l ;l 0 ) X m = min ( l ;l 0 ) C ( l ; l 0 ; m ) N ( l ;m ) tot (0 ; 0) N ( l 0 ;m ) tot ( q x ; q y ) r l + l 0 +1 ; (2.30) with X q = ( U x X q x = U x U y X q y = U y T x X q x = T x T y X q y = T y ) ; (2.31) C ( l ; l 0 ; m ) = ( l + l 0 )!( 1) m f ( 1) l + ( 1) l 0 g ; (2.32) r = q ( q x a + q y b cos r ) 2 + ( q y b sin r ) 2 : (2.33) where q is the gra y region illustrated in Figure 2.2(a). Because of the translational symmetry , the actual summation range can b e halv ed as in Figure 2.2(b) (the shaded area) while m ultiplying b y a factor of 2 on the ab o v e energy expression.

PAGE 28

18 The m ultip oles can b e re-written only with real quan tities as 1 X l =0 1 X l 0 =0 min ( l ;l 0 ) X m = min ( l ;l 0 ) N ( l ;m ) N ( l 0 ;m ) = 1 X l =0 1 X l 0 =1 N +( l ; 0) N +( l 0 ; 0) +2 1 X l =0 1 X l 0 =1 min ( l ;l 0 ) X m =1 N +( l ;m ) N +( l 0 ;m ) + N ( l ;m ) N ( l 0 ;m ) ; (2.34) for n uclear, electronic, and total quan tities. N + and N are the real and imaginary parts of N , resp ectiv ely . Since the summation o v er q x and q y has to span a large n um b er of cells, due to the 1 =r deca y rate of the Coulom b in teraction, the rotation of the m ultip ole axis and ev aluation of the energy for eac h distan t cell are the computationally demanding pro cedures. The FMM is a computationally attractiv e sc heme. The basic idea is that m ultip oles of distan t cells can b e group ed together if the distance b et w een the reference cell and distan t cells is large enough. A sc hematic description of the FMM treatmen t is giv en in Figure 2.3. As w e see in the gure, recursiv ely executing the FMM pro cedure con tributes to the p erformance of the algorithm. The zeroth-order of FMM recursion, denoted as F(0)=M (M for MPE and F for FMM), is the same as the MPE, as describ ed ab o v e except that only the distan t cells in the rst tier [F(0) region in the Figure 2.3] are considered. The rst recursion F(1) and later consist of three steps. The rst step is the reconstruction of m ultip oles for a larger cell using the set of m ultip oles dened for the cell of the previous recursion. In this w ork, w e c ho ose to group 9 (3 3) cells of the previous recursion and mak e one large cell in ev ery recursion. The redening of the m ultip oles in the larger cell is done b y shifting the origins of the m ultip oles to the common origin as illustrated in Figure 2.4 and adding them together. Let the origin shifting matrix b e M x y ( r ), where the r is the v ector connecting the old

PAGE 29

19 and new origins, then the large m ultip oles are dened b y the so-called \c hildren to paren t transformation" M 0 ( l ;m ) = 9cells X i l X =0 X = M ( ; ) M m l ( r i ) : (2.35) The origin shifting matrix for r = ( r ; ; ) is dened b y M m l ( r ) = r l P m l (cos ) e im ; (2.36) M m l ( r ) = M m l ; (2.37) and the new m ultip oles M ( l ;m ) are constructed b y M ( l ;m ) = ( N +( l ;m ) i N ( l ;m ) )( l + j m j )! for m 0 ; (2.38) M ( l ;m ) = ( N +( l ;m ) + i N ( l ;m ) )( l + j m j )! for m < 0 : (2.39) The summation o v er i for the nine cells can b e decomp osed in to v e. The rst of the v e is the cen ter cell, where the origins of the small and large cells coincide, th us the origin shifting matrix is the unit matrix. The other four corresp ond to step 1 to 4 in Figure 2.4. Due to the symmetry of the origin shifting matrix for v ectors r and r , there is a computational adv an tage to treat the t w o cells at a time. The second step of the three is the rotation of the axis to ha v e t w o m ultip oles dened with a common z -axis. This is p erformed b y con v ersion of the m ultip oles in to the spherical harmonics, the actual rotation, and nally bac k con v ersion of the rotated spherical harmonics in to m ultip oles. Since the spherical harmonics ha v e rotational in v ariance, w e ha v e to rst con v ert the m ultip oles in to the spherical

PAGE 30

20 harmonics b y r l Y ( l ;m ) = ( 1) m s (2 l + 1) 4 ( l j m j )! ( l + j m j )! ( 1) m < ( M ( l ;m ) ) + i ( 1) m +1 = ( M ( l ;m ) ) for m 0 ; (2.40) (2.41) r l Y ( l ;m ) = ( 1) m s (2 l + 1) 4 ( l j m j )! ( l + j m j )! M ( l ;m ) for m < 0 : (2.42) where M (or M 0 ) is the m ultip ole dened for the original size cell or a larger size cell, and Y (or Y 0 ) is the corresp onding spherical harmonic. The actual rotation is giv en b y Y n l ( ; ) = [ D l ( R )] T 1 Y m l ( ; ) ; (2.43) where Y is the rotated spherical harmonics b y the rotation matrix D . The explicit form of the rotation matrix is giv en in app endix C, but the rotation w e consider is limited due to the t w o-dimensional p erio dicit y as illustrated in Figure 2.5. A t rst, the m ultip oles at the reference cell and the m ultip oles at a distan t ( q x ; q y ) are dened at the xed and common Cartesian co ordinates as in Figure 2.5(a). In order to mak e the t w o m ultip oles ha v e a common z -axis, w e rotate the axis b y the angle along the z -axis to the coun ter-clo c kwise direction (Figure 2.5(b)). Then, w e tip the z -axis of the systems to w ard the common new x -axis b y 90 degree. These t w o steps complete the rotation. Though a general rotation is describ ed b y angles ; , and r , in the 2D case, and r are xed to 90 and 0 degrees, resp ectiv ely .

PAGE 31

21 The nal step is the bac k con v ersion of the spherical harmonics to the real m ultip oles. The op eration is dened b y N +( l ;m ) = ( 1) m s 4 (2 l + 1)( l + j m j )!( l j m j )! < ( r l Y ( l ; m ) + ( 1) m r l Y ( l ;m ) ) ; (2.44) N ( l ;m ) = ( 1) m s 4 (2 l + 1)( l + j m j )!( l j m j )! = ( r l Y ( l ; m ) ( 1) m r l Y ( l ;m ) ) : (2.45) Again, the quan tit y can b e dened in the original size cell or a larger size cell ( N 0 + ; N 0 made b y Y 0 ). As w e can see, the complex m ultip oles M , the real m ultip oles N + and N , and the spherical harmonics Y are closely in ter-related. Their in ter-connections are presen ted in app endix C. After prop er rotation is p erformed, quan tities suc h as the long-range energy and F o c k matrix con tributions are ev aluated. The FMM form ulas to ev aluate the long-range Coulom b in teraction energy and the F o c k matrix con tributions are E FMM LR = 2 f max X f =0 X q f C ( l ; l 0 ; m ) r l + l 0 +1 1 X l =1 1 X l 0 =1 N +( l ; 0) tot (0 ; 0) N 0 +( l 0 ; 0) tot ( f ) +2 1 X l =1 1 X l 0 =1 min ( l ;l 0 ) X m =1 N +( l ;m ) tot (0 ; 0) N 0 +( l 0 ;m ) tot ( f ) + N ( l ;m ) tot (0 ; 0) N 0 ( l 0 ;m ) tot ( f ) ; (2.46) F (LR) ( q ) = 2 f max X f =0 X q f C ( l ; l 0 ; m ) r l + l 0 +1 1 X l =0 1 X l 0 =1 N +( q )( l ; 0) (0 ; 0) N 0 +( l 0 ; 0) tot ( f ) +2 1 X l =0 1 X l 0 =1 min ( l ;l 0 ) X m =1 N +( q )( l ;m ) (0 ; 0) N 0 +( l 0 ;m ) tot ( f ) + N ( q )( l ;m ) (0 ; 0) N 0 ( l 0 ;m ) tot ( f ) : (2.47) with N 0 m ultip oles dened in a 3 f 3 f size cell. The summation index q f stands for F(f ) region in Figure 2.3. The regions for higher f are dened similarly .

PAGE 32

22 Summarizing the pro cedure of the FMM, (I) N ( l ;m ) ! M ( l ;m ) Generate the real m ultip oles N +( l ;m ) and N ( l ;m ) using the recursion form ula b y P erez-Jord a and Y ang [28] and con v ert them in to the complex form, M ( l ;m ) Then, en ter the FMM recursion lo op o v er f . (I I) M ( l ;m ) ! M 0 ( l ;m ) The c hildren to paren t transformation is p erformed to obtain the m ultip ole dened in 9 times larger cell. (I I I) M 0 ( l ;m ) ! Y 0 ( l ;m ) Con v ert the m ultip oles in to the spherical harmonics to allo w a rotaion of axes. (IV) Y 0 ( l ;m ) ! Y 0 ( l ;m ) The rotation to allo w t w o co ordinate systems of m ultip oles ha v e the common z -axis and parallel x and y -axes. (V) Y 0 ( l ;m ) ! N 0 ( l ;m ) ; Con v ert the rotated spherical harmonics in to the real m ultip oles. (VI) Ev aluation of the long-range energy and F o c k matrix con tribution. This ends the FMM lo op started in (I I). No w, w e rep eat with larger cells. F or 1D p olymers, steps (I I I), (IV), and (V) are skipp ed since w e ma y alw a ys c ho ose a system to b e aligned on the z -axis. In order to p erform the same pro cedure for 3D systems, t w o sligh t mo dications are required. First, rotations of the m ultip oles m ust b e p erformed in a general direction (unlik e just in the x y plane for 2D systema). Second, the FMM cells increase b y a factor of 27 (3 3 3) at eac h iteration. F or the most general case, 13 steps of the origin shifting op erations are required.

PAGE 33

23 2.2.3 Analytical Gradien ts of the T otal Energy p er Unit Cell The analytical expression of the total energy gradien ts with resp ect to displacemen t of an atom A in a unit cell in direction Q i ( i = x; y ; z ) is giv en b y @ E @ Q A i = X ; X q P ( q ) @ H ( q ) @ Q A i + 1 2 X ; X ; X q ; r ; s P ( q ) P ( s r ) @ @ Q A i ( ( 0 ) ( q ) j ( r ) ( s ) ) + m 1 X ; X q P ( q ) Z f ( 0 ) v X C @ ( q ) @ Q A i + @ ( 0 ) @ Q A i v X C ( q ) + r ( 0 ) g X C @ ( q ) @ Q A i + @ r ( 0 ) @ Q A i g X C ( q ) + ( 0 ) g X C @ r ( q ) @ Q A i + @ ( 0 ) @ Q A i g X C r ( q ) g d r m 2 4 X ; X ; X q ; r ; s P ( q ) P ( s r ) @ @ Q A i ( ( 0 ) ( r ) j ( q ) ( s ) ) X ; X q W ( q ) @ S ( q ) @ Q A i + @ E NR @ Q A i + ( @ E LR @ Q A i ) ; (2.48) where W is the energy-w eigh ted densit y matrix, whose elemen ts are dened as W ( q ) = 2 K o cc : X j BZ X k " [ k ] j c [ k ] j c [ k ] j exp ( i k q ) : (2.49) Those who are in terested in a more detailed deriv ation of the HF/KS analytical energy gradien t are referred to Hirata and Iw ata [6]. In calculations of 2D p erio dic systems, w e need to compute the energy gradien ts with resp ect to the lattice constan ts a; b , and r , the lengths of t w o lattice v ectors, and the angle formed b y the t w o lattice v ectors, as w ell as gradien ts with resp ect to eac h atomic displacemen t. W e can align the v ector a to lie parallel to the x -axis and restrict the v ector b to b e in the x y plane, th us making the angle r coplanar with the x y plane. Under this picture, let the cen ter of the atom A in the reference cell as ( Q A (0 ; 0) X ; Q A (0 ; 0) Y ; Q A (0 ; 0) Z ), then the cen ter of a basis set Q at the unit cell ( q x , q y ) on

PAGE 34

24 atom A ma y b e written as Q A ( q x ;q y ) X = q x a + q y b cos r + Q A (0 ; 0) X ; (2.50) Q A ( q x ;q y ) Y = q y b sin r + Q A (0 ; 0) Y ; (2.51) Q A ( q x ;q y ) Z = Q A (0 ; 0) Z ; (2.52) for a general 2D case. By virtue of the fact that eac h of the @ E @ a , @ E @ b , @ E @ r is dep enden t on @ E @ Q x , @ E @ Q y , and @ E @ Q z , using the c hain rule, w e can write the deriv ativ es with resp ect to the unit cell lengths and the angle as @ E @ a = X A X q x X q y @ Q A ( q x ;q y ) X @ a @ E @ Q A ( q x ;q y ) X + @ Q A ( q x ;q y ) Y @ a @ E @ Q A ( q x ;q y ) Y + @ Q A ( q x ;q y ) Z @ a @ E @ Q A ( q x ;q y ) Z ! = X A X q x X q y q x @ E @ Q A ( q x ;q y ) X ; (2.53) @ E @ b = X A X q x X q y @ Q A ( q x ;q y ) X @ b @ E @ Q A ( q x ;q y ) X + @ Q A ( q x ;q y ) Y @ b @ E @ Q A ( q x ;q y ) Y + @ Q A ( q x ;q y ) Z @ b @ E @ Q A ( q x ;q y ) Z ! = X A X q x X q y q y cos r @ E @ Q A ( q x ;q y ) X + q y sin r @ E @ Q A ( q x ;q y ) Y ! ; (2.54) @ E @ r = X A X q x X q y @ Q A ( q x ;q y ) X @ r @ E @ Q A ( q x ;q y ) X + @ Q A ( q x ;q y ) Y @ r @ E @ Q A ( q x ;q y ) Y + @ Q A ( q x ;q y ) Z @ r @ E @ Q A ( q x ;q y ) Z ! = X A X q x X q y q y b sin r @ E @ Q A ( q x ;q y ) X + q y b cos r @ E @ Q A ( q x ;q y ) Y ! ; (2.55)

PAGE 35

25 The deriv ativ e expression with resp ect to the cell parameter a is iden tical to that for 1D systems [4, 5]. A similar but sligh tly more complicated pro cedure w ould yield the cell parameter deriv ativ e expressions for 3D systems. The form ulas are presen ted in app endix B. With the form ulas giv en ab o v e, the force with resp ect to the unit cell parameters is accum ulated as the force with resp ect to atomic displacemen ts is ev aluated. Hence, the analytical gradien ts with resp ect to the lattice constan ts a; b; c; ; , and r can b e directly obtained with marginal additional computational costs when analytical gradien ts with resp ect to atomic displacemen ts are computed. 2.2.4 Grid W eigh t Deriv ativ es with Resp ect to Unit Cell P arameters When w e ev aluate the analytical gradien ts with resp ect to cell parameters in grid-based DFT for p erio dic systems, w e need to compute one more term in addition to Eq.(2.48), whic h is the deriv ativ e of the grid w eigh t. The follo wing deriv ation for the grid deriv ativ e is based on Bec k e's original n umerical quadrature [38] and its extension b y T reutler and Ahlric hs [47]. The notation used here closely follo ws theirs. The exc hange correlation energy p er unit cell is giv en as E X C = Z cel l f ( r ) d r cel l X i X N p i ( r N ) f ( r N ) ; (2.56) where p i ( r N ) is the w eigh t of the grid and f ( r N ) is the v alue ev aluated at r N . The grid w eigh t expression at r for i -atom cen tered grid is giv en as p i ( r ) = P ( 0 ) i ( r ) P n; q P ( q ) n ( r ) : (2.57)

PAGE 36

26 Straigh tforw ardly , the grid w eigh t deriv ativ e with resp ect to a unit cell parameter p = ( a; b; c; ; ; r ) is @ p i ( r ) @ p = ( P ( 0 ) i ( r )) 0 P n; q P ( q ) n ( r ) P ( 0 ) i ( r )( P n; q P ( q ) n ( r )) 0 ( P n; q P ( q ) n ( r )) 2 : (2.58) Th us, for a giv en P ( q ) i ( r ), w e need to compute its deriv ativ es. P ( q ) i ( r ) itself is expressed as P ( q ) i ( r ) = j ( q 0 ) 6 = i ( 0 ) s k ( ij ) ; (2.59) where s k ( ij ) = 1 2 [1 g k ( ij )] ; (2.60) g k ( ij ) = h f h [ ::h | {z } k times ( ij )] g ; (2.61) h ( ij ) = 3 2 ij 1 2 3 ij ; (2.62) ij = ij + a ij (1 2 ij ) ; (2.63) ij = r i r j R ij ; (2.64) a ij = u ij u 2 ij 1 ; (2.65) u ij = 1 + 1 ; (2.66) = s R i R j : (2.67) Th us the deriv ativ e of P ( q ) i ( r ) is constructed b y s k and the deriv ativ e of s k . The deriv ativ e for eac h s k ( ij ) is computed b y the c hain rule @ s k ( ij ) @ p = @ s k ( ij ) @ g 3 @ g 3 @ g 2 @ g 2 @ g 1 @ g 1 @ ij @ ij @ ij @ ij @ p ; (2.68) where r i and r j are the lengths of atom i and j from the grid p oin t r = ( x e ; y e ; z e ), resp ectiv ely , R ij is the relativ e distance b et w een atoms i and j , and R i and R j

PAGE 37

27 denote the Bragg-Slater radius [48] of atom i and j . r i = p ( x i x e ) 2 + ( y i y e ) 2 + ( z i z e ) 2 (2.69) R ij = q ( x i x j ) 2 + ( y i y j ) 2 + ( z i z j ) 2 (2.70) Eac h term in Eq.(2.68) is giv en as @ s k ( ij ) @ g 3 = @ [ 1 2 (1 g 3 )] @ g 3 = 1 2 ; (2.71) @ g 3 @ g 2 = @ ( 3 2 g 2 1 2 g 3 2 ) @ g 2 = 3 2 3 2 g 2 2 ; (2.72) @ g 2 @ g 1 = @ ( 3 2 g 1 1 2 g 3 1 ) @ g 1 = 3 2 3 2 g 2 1 ; (2.73) @ g 1 @ ij = @ ( 3 2 ij 1 2 3 ij ) @ ij = 3 2 3 2 2 ij ; (2.74) @ ij @ ij = 1 2 a ij ij ; (2.75) and @ ij @ m = @ ( r i r j R ij ) @ m = @ ( r i r j R ij ) @ x i @ x i @ p + @ ( r i r j R ij ) @ y i @ y i @ p + @ ( r i r j R ij ) @ z i @ z i @ p + @ ( r i r j R ij ) @ x j @ x j @ p + @ ( r i r j R ij ) @ y j @ y j @ p + @ ( r i r j R ij ) @ z j @ z j @ p + @ ( r i r j R ij ) @ x e @ x e @ p + @ ( r i r j R ij ) @ y e @ y e @ p + @ ( r i r j R ij ) @ z e @ z e @ p : (2.76)

PAGE 38

28 The deriv ativ e terms with resp ect to l (= x; y ; z ) can b e computed with @ ( r i r j R ij ) @ l i = 1 R ij l i l r r i 1 R 3 ij ( l i l j )( r i r j ) ; (2.77) @ ( r i r j R ij ) @ l j = 1 R ij l j l r r j + 1 R 3 ij ( l i l j )( r i r j ) ; (2.78) @ ( r i r j R ij ) @ l e = 1 R ij l i l e r i l j l e r j : (2.79) T erms @ l @ p ( l = x; y ; z ) and ( p = a; b; r ) are giv en in Eq.(2.53), (2.54), and (2.55) for 2D systems and in app endix B for 3D systems. This completes the deriv ation of the grid w eigh t deriv ativ es with resp ect to the lattice constan ts in p erio dic systems. It should b e noted that the grid deriv ativ e term is an essen tial quan tit y when the gradien t with resp ect to unit cell parameters is computed regardless of the neness of the grid. It is rep orted that the grid w eigh t deriv ativ e term has a negligible eect when the grid is ne enough for molecular systems [49]. This is true for the gradien ts with resp ect to atomic displacemen ts but not to the unit cell parameters. The reason wh y the term is essen tial for the gradien ts with resp ect to unit cell parameters lies in the atomic partitioning of the grid. The grid b oundary dep ends on the unit cell parameters, and innitesimal c hanges to the unit cell parameters c hanges the p osition and w eigh t of the grid p oin ts in suc h a w a y that the grid deriv ativ es giv e a nite con tribution to the gradien t ev en in the limit of an innitely ne grid. The c hange in the grid w eigh ts w ould b e signican t esp ecially near the b oundary b et w een the reference unit cell and the neigh b oring cell, and hence the grid deriv ativ e can alternativ ely b e ev aluated as the surface in tegral at the b oundary , where the surface in tegral is transformed in to a three dimensional space in tegral b y virtue of Gauss theorem [7].

PAGE 39

29 2.2.5 The MPE/FMM Energy Deriv ativ es The deriv ativ e of the m ultip ole expansion energy with resp ect to a displacemen t of an atom A in the Q i (i= x; y ; z ) direction is deriv ed b y taking the rst deriv ativ e of Eq.(2.46) with resp ect to Q A i , @ E FMM LR @ Q A i = 2 f max X f =0 X q f " 1 X l =0 1 X l 0 =1 C ( l ; l 0 ; 0) r l + l 0 +1 @ N +( l ; 0) tot @ Q A i N 0 +( l 0 ; 0) tot + N +( l ; 0) tot @ N 0 +( l 0 ; 0) tot @ Q A i +2 1 X l =0 1 X l 0 =1 min ( l ;l 0 ) X m =1 C ( l ; l 0 ; m ) r l + l 0 +1 @ N +( l ;m ) tot @ Q A i N 0 +( l 0 ;m ) tot + N +( l ;m ) tot @ N 0 +( l 0 ;m ) tot @ Q A i + @ N ( l ;m ) tot @ Q A i N 0 ( l 0 ;m ) tot + N ( l ;m ) tot @ N 0 ( l 0 ;m ) tot @ Q A i # : (2.80) Quan tities w e need are @ N ( l ;m ) tot @ Q A i N 0 ( l 0 ;m ) tot + N ( l ;m ) tot @ N 0 ( l 0 ;m ) tot @ Q A i = @ ( N ( l ;m ) n uc + N ( l ;m ) ele ) @ Q A i N 0 ( l 0 ;m ) tot + N ( l ;m ) tot @ ( N 0 ( l ;m ) n uc + N 0 ( l ;m ) ele ) @ Q A i ; (2.81) where @ N ( l ;m ) n uc @ Q A i = # atoms X N =1 @ ( Z A N ( l ;m ) A ) @ Q A i = # atoms X N =1 @ ( Z A N ( l ;m ) A ) @ Q A i AN = Z Q @ N ( l ;m ) Q @ Q A i ; (2.82) and @ N ( l ;m ) ele @ Q A i = X q X P ( q ) @ h ( 0 ) j N ( l ;m ) j ( q ) i + @ ( P ( q ) ) N ( q )( l ;m ) @ Q A i = X q X P ( q ) h @ ( 0 ) @ Q A i j N ( l ;m ) j ( q ) i A + h ( 0 ) j N ( l ;m ) j @ ( q ) @ Q A i i A + @ ( P ( q ) ) N ( q )( l ;m ) @ Q A i # ; (2.83)

PAGE 40

30 A = 8 > > < > > : 1 if the A O is cen tered on the atom A ; 0 otherwise : The last term is included in the energy w eigh ted densit y matrix term in Eq.(2.48). Th us, the new quan tities are required to ev aluate Eq.(2.80): the n uclear momen t deriv ativ es and the momen t in tegral deriv ativ es. The pro cedure to ev aluate Eq.(2.80) is similar to that for the long-range energy con tribution. F or deriv ativ es with resp ect to unit cell parameters, w e need additional terms as w ell as the terms in (2.53),(2.54), and (2.55) due to the dep endence of the denominator r l + l 0 +1 on the cell parameters. The additional term is written as @ E FMM LR @ k = X R C 0 ( l ; l 0 ; m ) r l + l 0 +3 l max X l =1 l max X l 0 =1 N +( l 0 ; 0) tot N 0 +( l 0 ; 0) tot g ( p ) +2 l max X l =1 l max X l 0 =1 min ( l ;l 0 ) X m =1 ( N +( l 0 ;m ) tot N 0 +( l 0 ;m ) tot + N ( l 0 ;m ) tot N 0 ( l 0 ;m ) tot ) g ( p ) ; (2.84) where p is the lattice parameters a , b , or r for 2D systems and C 0 ( l ; l 0 ; m ) = C ( l ; l 0 ; m ) ( l + l 0 + m ) ; (2.85) g ( a ) = q 2 x a + q x q y b cos r ; (2.86) g ( b ) = q x q y b cos r + q 2 y b; (2.87) g ( r ) = q x q y ab sin r : (2.88) 2.3 Results In this section, w e discuss the n umerical p erformance of the metho ds. First, w e analyze the n umerical p erformance of treating the long-range in teractions, then w e discuss the analytical gradien ts. A t w o-dimensional h ydrogen-ruoride

PAGE 41

31 sheet is c hosen as a mo del 2D system as sho wn in Figure 2.6. The STO-3G basis set is used. Brillouin zone sampling is p erformed on an 8 8 mesh. Accuracy is conrmed with a calculation with a ner mesh. The CPU time to include the long-range in teractions b y the MPE and FMM is measured. Figure 2.7 sho ws plots of the n um b er of cells included in the MPE calculations v ersus CPU time needed for one SCF iteration with or without MPE. When w e compare the plot without MPE, lab eled b y \Exact" in the legend of the gure, and the plots with MPE, w e observ e that the slop es are dieren t. The cost for ev aluating the long-range in teraction b y MPE is prop ortional to the n um b er of cells included. Since the calculation is of a 2D system, the cost scaling observ ed agrees with the formal scaling of O ( D 2 ). Changing the order of the m ultip ole expansion, L max from 2 to 4,6, or 8 aects the prefactor of the cost, but the scaling order is main tained since the four sets are parallel. Figure 2.8 sho ws the timing information obtained b y FMM. The gure sho ws a dramatic decrease in computational cost. Using FMM, with L max =2, w e can include more than a million times larger n um b ers of distan t cells b y just doubling the computational eort. The scaling prop ert y is similar for dieren t L max . F or L max =8, the increase in the computation time to include a million times larger n um b ers of distan t cells is ab out a factor of three. Another result w e should p oin t out is the gain b y going from MPE to FMM. Comparing Figure 2.7 and Figure 2.8, w e see that, in this particular case, w e already ac hiev e a sp eed up b y a factor of t w o orders of magnitude. Since their cost scaling is dieren t, the gain of FMM o v er the MPE w ould b e increasingly larger as w e include more long-range eects. W e applied our metho ds not only to the mo del h ydrogen ruoride 2D stac k but also to an exp erimen tally kno wn system. Boron nitride is one of the b est studied t w o-dimensional p erio dic insulators. Con v ergence b eha vior of the analytical gradien ts is discussed with resp ect to three parameters: SCF con v ergence

PAGE 42

32 (10 S C F ), long range Coulom b truncation b y FMM ( f ), and in tegral ev aluation cuto (10 I N T ). It is customary in molecular calculations to tigh ten the SCF con v ergence criteria b ecause it is computationally inexp ensiv e. Ho w ev er, in p erio dic systems, esp ecially in 2D and 3D, ev en Hartree-F o c k calculations can b e v ery time-consuming and in man y cases the t w o-electron in tegrals ha v e to b e calculated b y the direct metho d [41, 42]. Th us, reducing the n um b er of SCF iterations b y setting lo w er con v ergence criteria has a certain computational adv an tage. The truncation parameter of the explicit ev aluation of the Coulom b in teractions and the in tegral cuto parameter also signican tly aect the run time of the calculations. The long-range Coulom b truncation parameter denes the b oundary b et w een the exact, y et computationally in tensiv e, region and the appro ximate, y et computationally ecien t, region. The in tegral cuto parameter has to b e c hosen carefully . Setting a cuto parameter to o aggressiv ely w ould destro y the balance b et w een the electron-electron repulsion term and the other terms, whic h ma y lead to erroneous results. The main fo cus of these in v estigations is to nd the most aggressiv e parameters whic h do not sacrice the accuracy w e need. Calculations are p erformed with a hexagonal geometry , the B-N distances 1.50 A and 1.30 A. Kno wing that the exp erimen tal equilibrium b ond distance is 1.45 A [50], the ev aluation of analytical gradien ts at the t w o geometries allo ws us to see the con v ergence of analytical gradien ts at near and far p oin ts from the equilibrium. Brillouin zone sampling is done with a 10 10 mesh for all the calculation using the STO-3G basis set. Calculations with the 3-21G basis set is p erformed with a 12 12 mesh. T able 2.1 sho ws the total energy p er unit cell, band gap (the dierence b et w een the highest o ccupied crystalline orbital and the lo w est uno ccupied crystalline orbital energies), and the energy gradien t p er unit cell with resp ect to x -direction of the b oron atom and with resp ect to the translational p erio d a and the unit cell angle r as w ell as the n um b er of SCF iterations needed to meet the

PAGE 43

33 dened con v ergence criterion. The SCF calculation is considered to b e con v erged when the largest of the dierences b et w een the curren t and the previous densit y matrix elemen ts go es b elo w the user-dened con v ergence criterion 10 S C F . S C F is v aried from 4 to 8 with an incremen t of 1. The total energy con v erges up to 10 5 Hartree with S C F =7, while the band gap con v erges m uc h faster, whic h agrees with w ell-kno wn b eha vior. It is observ ed that con v ergence of the energy gradien t with resp ect to atomic displacemen t is faster than the gradien t with resp ect to lattice parameters. This can b e understo o d from Eqs.(2.53) to (2.55). Since the energy gradien t with resp ect to the unit cell parameters are calculated using the gradien t with resp ect to atomic displacemen t m ultiplied b y factors related to the lattice constans and the cell p osition, the error in unit cell parameters gradien ts are amplied b y these factors. In order to obtain four decimal accuracy w e only need S C F =6. Once all the gradien ts go b elo w 1 : 0 10 3 , w e consider the geometries to b e w ell con v erged for most cases. Based on this fact, four decimal accuracy in gradien ts ma y b e enough, if the SCF pro cedure is the computational b ottlenec k. T able 2.2 sho ws the con v ergence with resp ect to the in tegral cuto parameter I N T , where only the t w o-electron in tegrals whose absolute v alue is larger than 10 I N T are considered in the calculation. The general trend is the same as that observ ed in T able 2.1. The band gap con v erges rapidly and the gradien t with resp ect to the atomic displacemen t con v erges faster than that with resp ect to lattice constan ts. The timings to ev aluate the analytical gradien ts are giv en, whic h are measured on an IBM RS/6000 Mo del 397, 160MHz with 1GB of memory . The CPU time required for ev aluating analytical gradien ts is a strong function of the in tegral cuto. Though, I N T = 4 is to o aggressiv e to b e meaningful, setting a mo dest cuto, I N T = 8, con tributes time sa ving without sacricing results.

PAGE 44

34 T able 2.3 sho ws the con v ergence with resp ect to the long-range Coulom b cuto parameter f . The parameter is v aried from 0 to 5 as w ell as `None'. `None' means that the calculations are p erformed only within the dark grey region in Figure 2.9, f =0 includes the ligh t grey region b y MPE and f =1,2,.. includes the outer region b y FMM for f times recursiv ely . W e see that the energy gradien ts do not con v erge at all when w e do not include enough of the long-range eects. The cost ev aluating the long-range eect b y FMM is negligible compared to the cost for calculating the short-range region. Th us, w e suggest that w e alw a ys c ho ose a large enough cuto parameter. f =5 ma y b e enough to incorp orate the imp ortan t part of the long-range con tributions. W e p erformed geometry optimizations with and without including the longrange Coulom b part as a test. In Figure 2.10, the righ t b ottom cell surrounded b y b old lines is the reference cell. A full geometry optimization determines the B-N b ond length, t w o translational p erio ds a and b , and the angle r . When the long-range eects are neglected w e obtain the B-N b ond length 1.428 A, a =4.707 A, b =4.708 A, and r =58.3 degrees. Based on these geometrical parameters t w o other B-N b onds in the Figure (brok en lines across the cells) are calculated to b e 1.444 A for b oth b onds. The truncation of the long-range eects sligh tly shifts the minim um of the p oten tial energy surface from the hexagonal symmetry p oin t. Note that abruptly cutting o the Coulom b eect leads to man y other problems. F or example, in the hexagonal b oron-nitride system, the Mullik en c harge on the 2 p x t yp e basis function and the 2 p y t yp e basis function should b e the same. Monitoring this c harge helps estimate ho w go o d our long-range treatmen t is. Once the long-range eects are included b y FMM, w e obtain the iden tical b ond lengths for all the B-N b onds, and also the condition a = b = p 3( B N ) and r =60 degrees are satised.

PAGE 45

35 Understanding these issues related to the innite summation of the Coulom b eect, geometry optimizations are p erformed with the hexagonal geometrical constrain t. The n umerical study sho w ed us that the energy gradien ts with resp ect to atomic displacemen t con v erge faster than those for lattice constan ts. Th us, calculating energy gradien ts only for the B-N b ond and generating the lattice constan ts based on the B-N b ond length will b e more ecien t using the relationship a = b = p 3( B N ). T able 2.4 sho ws the optimized B-N b ond length, band gap and Mullik en c harge on eac h atom. The B-N b ond length is consisten t b et w een HF and B3L YP with b oth basis sets, STO-3G and 3-21G, and is also in go o d agreemen t to the exp erimen tal b ond length [50] 1.45 A. The band gap, ho w ev er, sho ws v ery dieren t b eha vior b et w een the HF and B3L YP results. The fact that HF o v erestimates the band gap and that the correlation eect decreases the HF band gap is w ell kno wn, therefore the result is not surprising. The prior results for the band gap using a semi-empirical theory [51], the one-particle man y-b o dy Green-F unction tec hnique [52] and the Laplace-transformed second-order p erturbation theory [53] qualitativ ely coincide with our result. Since the latter t w o w orks presen t only STO-3G basis set results, this pap er is the rst to presen t some insigh t in to the basis set eects on the band gap at this lev el of correlation treatmen t. Although the 3-21G basis is not large, it giv es quan titativ e results m uc h closer to large basis set calculations suc h as 6-31G** for a one-dimensional p olymer [8], while the STO-3G basis set often giv es erroneous results. The B3L YP/3-21G result should b e the highest lev el at presen t at least from quan tum c hemical approac hes for extended system. Compare with solid state application [54, 55] where a large basis set is used. The B3L YP/3-21G result agrees to the exp erimen tal band gap [56] within 0.5 eV. But, this ma y b e fortuitous due to the w ell kno wn tendency of DFT metho ds to underestimate the band gap. The basis set eects and dierences

PAGE 46

36 b et w een exp erimen ts and calculations has to b e b etter calibrated. Sensitivit y of the band gap to the basis set size is p oin ted out b y V ra ck o et al. Also a 0.4 eV dierence in calculated band gap for the mono-la y er and crystal is rep orted [54]. The Mullik en c harges on eac h atom are also sensitiv e to the c hoice of basis sets. The results b y the STO-3G and 3-21G basis sets dier as m uc h as ab out 0.5. 2.4 Conclusion W e ha v e giv en a formalism for the analytical total energy gradien t for Hartree-F o c k and densit y functional theories for p erio dic systems, with some emphasis on t w o-dimensional systems. The analytical gradien ts of the total energy are implemen ted including those with resp ect to lattice constan ts. The long-range Coulom b in teraction energy and its gradien ts are treated ecien tly b y FMM. Numerical examples conrm that the long-range Coulom b treatmen t is an essen tial quan tit y for the total energy . W e also sho w that the long-range eects are an imp ortan t piece of energy gradien ts. The analytical gradien ts including fast computation of the long-range eects are nessesary set up in order to p erform geometry optimizations and harmonic frequency calculations practically . Numerical studies of con v ergence with resp ect to sev eral parameters sho w that S C F =6, I N T =8, and f =5 ma y b e used to obtain useful results without sacricing accuracy .

PAGE 47

37 T able 2.1: Con v ergence of the total energy p er unit cell, band gap, and analytical gradien ts of the total energy and timing for the gradien ts ev aluation with resp ect to the SCF con v ergence criterion, S C F ev aluated at B-N b ond lengths of 1.50 A and 1.30 A. The t w o-electron in tegral cuto parameter, I N T =14 and long-range eect cuto parameter f =5 are xed. The total energies and gradien ts are giv en in atomic units, band gaps are in eV. B-N b ond length: 1.50 A SCF E H F Band Gap @ E @ Q B x @ E @ a @ E @ r Iteration 4 -78.28557 13.41 -0.0035 0.0651 -0,2506 6 5 -78.28273 13.46 -0.0035 0.0645 -0.2488 9 6 -78.28218 13.46 -0.0035 0.0645 -0.2489 11 7 -78.28219 13.46 -0.0035 0.0645 -0.2489 13 8 -78.28219 13.46 -0.0035 0.0645 -0.2489 15 B-N b ond length: 1.30 A SCF E H F Band Gap @ E @ Q B x @ E @ a @ E @ r Iteration 4 -78.22905 14.74 -0.0058 -0.1220 0.2149 5 5 -78.22507 14.74 -0.0058 -0.1203 0.2108 7 6 -78.22487 14.75 -0.0058 -0.1204 0.2110 11 7 -78.22487 14.75 -0.0058 -0.1204 0.2110 13 8 -78.22487 14.75 -0.0058 -0.1204 0.2110 14

PAGE 48

38 T able 2.2: Con v ergence of the total energy p er unit cell, band gap, and analytical gradien ts of the total energy and timing for the gradien ts ev aluation with resp ect to the t w o-electron in tegral cuto parameter, I N T ev aluated at B-N b ond lengths of 1.50 A and 1.30 A. The SCF con v ergence criterion S C F =8 and long-range effect cuto parameter f =5 are xed. The total energies and gradien ts are giv en in atomic units, band gaps are in eV, and timings are in seconds. B-N b ond length: 1.50 A INT E H F Band Gap @ E @ Q B x @ E @ a @ E @ r Timing 4 -79.16506 5.46 -0.0173 0.4949 -1.3814 2800 5 -78.27210 13.46 -0.0028 0.0650 -0.2481 4000 6 -78.28217 13.62 -0.0034 0.0641 -0.2474 5400 7 -78.28223 13.46 -0.0035 0.0648 -0.2496 7000 8 -78.28221 13.46 -0.0035 0.0645 -0.2490 8900 9 -78.28219 13.46 -0.0035 0.0645 -0.2489 11000 10 -78.28219 13.46 -0.0035 0.0645 -0.2489 13300 B-N b ond length: 1.30 A INT E H F Band Gap @ E @ Q B x @ E @ a @ E @ r Timing 4 Did not con v erge 5 -78.23276 16.87 -0.0042 -0.1000 0.1615 6400 6 -78.22725 14.73 -0.0058 -0.1165 0.2013 8900 7 -78.22480 14.75 -0.0058 -0.1205 0.2144 11700 8 -78.22488 14.74 -0.0058 -0.1204 0.2100 15000 9 -78.22487 14.75 -0.0058 -0.1204 0.2100 18500 10 -78.22487 14.75 -0.0058 -0.1204 0.2100 22700

PAGE 49

39 T able 2.3: Con v ergence of the total energy p er unit cell, band gap, and analytical gradien ts of the total energy with resp ect to the long-range Coulom b eect cuto parameter f ev aluated at B-N b ond lengths of 1.50 A and 1.30 A. The SCF con v ergence criterion S C F =8 and the t w o-electron in tegral cuto parameter, I N T =14 are xed. The total energies and gradien ts are giv en in atomic units, band gaps are in eV. B-N b ond length: 1.50 A f E H F Band Gap @ E @ Q B x @ E @ a @ E @ r None -78.27728 13.20 -0.0004 0.0641 -0.2526 0 -78.28008 13.35 -0.0021 0.0645 -0.2496 1 -78.28146 13.42 -0.0030 0.0645 -0.2491 2 -78.28195 13.48 -0.0033 0.0645 -0.2490 3 -78.28212 13.46 -0.0035 0.0645 -0.2490 4 -78.28217 13.46 -0.0035 0.0645 -0.2490 5 -78.28219 13.46 -0.0035 0.0645 -0.2490 B-N b ond length: 1.30 A f E H F Band Gap @ E @ Q B x @ E @ a @ E @ r None -78.21721 14.34 -0.0003 -0.1213 0.2072 0 -78.22139 14.56 -0.0030 -0.1206 0.2109 1 -78.22364 14.68 -0.0048 -0.1205 0.2111 2 -78.22461 14.72 -0.0055 -0.1205 0.2111 3 -78.22474 14.74 -0.0057 -0.1204 0.2111 4 -78.22484 14.74 -0.0058 -0.1204 0.2111 5 -78.22487 14.75 -0.0058 -0.1204 0.2110 T able 2.4: The optimized geometry , band gap and Mullik en c harge of hexagonal b oron nitride Metho d B-N Band Gap Charge N Charge B HF/STO-3G 1.44 13.8 7.554 4.446 B3L YP/STO-3G 1.46 5.7 7.444 4.556 B3L YP/3-21G 1.46 6.3 7.936 4.064 Exp erimen t 1.45 [50] 5.8 [56]

PAGE 50

40 X Y ¡ a b (q x ,q y )=(0,0) (1,0) (2,0) (0,1) (1,1) Figure 2.1: The denition of the unit cell parameters a , b , and r and the cell indices q x and q y for a general t w o-dimensional p erio dic system.

PAGE 51

41 0 4 13 -4 -13 0 5 -5 -14 14 0 4 13 -4 -13 0 5 -5 -14 14 R (a) (b) Figure 2.2: The region where the long-range Coulom b in teractions are computed b y m ultip ole expansion. The reference cell is lo cated at the cen ter. In this particular gure U X = 13, U Y = 14, T X = 4 and T Y = 5 are c hosen. (a) the in teraction b et w een the reference cell and eac h cell in the shaded region is computed b y the MPE. (b) F or the F o c k matrix correction and the total energy correction calculation, the total con tribution is exactly the t wice the con tribution from the shaded region R.

PAGE 52

42 F(2) F(1) M=F(0) E Figure 2.3: Illustration of the FMM. The reference cell is at the left b ottom corner. In the E region, the Coulom b in teractions are calculated b y the t w o-electron in tegrals, while at the F(i) (i=0,1,2,..) the Coulom b in teractions are appro ximated b y FMM. The F(0) region is equiv alen t to the MPE with m ultip oles b eing dened for cells of the original size, whereas distan t cells are regroup ed in to larger cells as in F(1) and F(2).

PAGE 53

43 Original Step 1. Step 2. Step 3. Step 4. Final Figure 2.4: The four steps in constructing a new set of m ultip oles dened for a consolidated 3 3 cell.

PAGE 54

44 (a) x x y y (c) z'' z'' y''= y' y''= y' x'' (b) x' x' y' y' a z'= z z z z'= z x'' Figure 2.5: The pro cedure for axis rotation to mak e the t w o co ordinate systems cen tered on the xy -plane ha v e a common z -axis

PAGE 55

45 F H F H F H F H F H F H F H F H .............. ..............b a ¡ Figure 2.6: Geometry of mo del h ydrogen ruoride c hain stac k. a =4.8497 A, b =4.0 A, and r = 90 degrees.

PAGE 56

46 Performance of MPE Number of cells includedTiming per SCF iteration / sec 10 2 10 3 10 4 10 5 10 1 10 2 10 3 10 4 10 5 L max =2 L max =4 L max =6 L max =8 Exact Figure 2.7: A comparison of timings among the m ultip ole expansions of v arious orders and the explicit treatmen t of long-range part as a function of n um b er of cells included in the calculation. Timings are giv en in seconds.Timing per SCF iteration / secPerformance of FMM Number of cells included 10 4 10 5 10 6 10 7 10 8 10 9 10 10 10 11 10 1 10 2 10 3 10 4 L max =2 L max =4 L max =6 L max =8 Figure 2.8: A comparison of timings among the fast m ultip ole metho d of v arious order. Timings are giv en in seconds.

PAGE 57

47 Figure 2.9: La y ered structure for a Boron-nitride calculation. Blac k cell is the reference unit cell, the dark gra y cells are treated b y the explicit calculation of t w o-electron in tegrals, the ligher gra y cells are treated b y the m ultip ole expansion and the ligh test gra y cells are treated b y the fast m ultiple metho d. B N N B B N B Figure 2.10: Three symmetrically equiv alen t but dieren tly treated B-N b onds (brok en lines).

PAGE 58

CHAPTER 3 F ORMAL ASPECTS OF A TOMIC-ORBIT AL BASED COUPLED-CLUSTER THEOR Y 3.1 In tro duction Man y-b o dy p erturbation theory (MBPT) and coupled-cluster theory (CC) ha v e pro v en to b e useful to ols due to their ecien t and p o w erful inclusion of electron correlation. The second-order MBPT [MBPT(2)] and CC singles, doubles, and p erturbativ e triples, CCSD(T) [57, 58] corrections are often our metho ds of c hoice in order to obtain quan titativ e results and the n um b er of n umerical applications using these metho ds is almost coun tless. The formal construction of these metho dologies suc h as CCD [59{61], CCSD [62], CCSDT [63], and CCSDTQ [64] are all dev elop ed based on an exp onen tial ansatz. Their p erturbativ e [57, 58] and non-p erturbativ e [65] appro ximations, excited states [66{ 68], and analytical gradien ts [69{71] all deriv e from an exp onen tial ansatz as w ell. Although the ground state, excited state and analytical gradien t formalisms lo ok v ery dieren t at a glance, the actual computations are v ery similar. All the theories are presen ted as manipulation of in tegrals and CC amplitudes, and the expressions can b e readily deriv ed with the aid of the diagrammatic tec hnique [72]. In practice, computations are p erformed as matrix-matrix m ultiplications. Size-extensivit y [73] is one of the features in MBPT/CC theory , and is essen tial to treat large systems suc h as p olymers. Its rigorous denition is the absence of \unlink ed diagrams" whic h guaran tees correct scaling with the n um b er of electrons [74]. Size-extensiv e electron-correlation metho ds for large systems including p olymers are often v ery exp ensiv e since the MBPT/CC computational 48

PAGE 59

49 scaling is giv en as p o w ers of system size. F or example, MBPT(2) has an O ( N 5 ) in tegral transformation step, and CCSD(T) has the iterativ e O ( N 6 ) and one-time O ( N 7 ) computational step with N b eing the n um b er of basis functions. Th us, the MBPT/CC pro cedure b ecomes impractical for large systems due to the large N . F or large systems, it will b e more ecien t to use atomic-orbital (A O) based metho ds rather than the con v en tional MO based metho ds describ ed ab o v e. The A O based metho ds ha v e O ( N ) scalings with a m uc h larger prefactor than the corresp onding con v en tional metho ds for the large system limit due to the shortrange nature of the electron correlation. F or small molecules, the con v en tional metho ds are more ecien t due to the small prefactor, but there is a cross-o v er p oin t where the A O based metho ds b ecome more ecien t. Dev elopmen ts in linear scaling metho ds are v ery activ e starting from Hartree-F o c k and densit y functional theories as w ell as MBPT(2) and CC for molecules and p olymers. Regardless of the large amoun t of w ork in A O based metho ds, their metho dological dev elopmen ts are based on the corresp onding molecular orbital (MO) treatmen t. Starting from the already established MO based form ulation, using the MO to A O transformation form ula on the nal equation yields a new set of nal equations in A O basis. This is the most exp edien t w a y to obtain the w orking equations, with giv en MO equations. W e are able to obtain the iden tical A O-based equations b y using the MO to A O transformation at an earlier stage. The creation and annihilation op erators, zeroth-order Hamiltonian, excitation amplitudes and the reference and excited congurations, whic h are giv en b y a Slater determinan t, can b e written as a linear com bination of atomic-orbital based quan tities. The A O-based equations can b e deriv ed based up on suc h quan tities. Ho w ev er, the A O-based quan tities ha v e MO indices due to the fact that the F ermi v acuum state is dened only in the MO basis. Th us, the deriv ation pro cess is h ybrid b et w een MO and A O, and only the

PAGE 60

50 resuling equations are giv en solely in A O. In this pap er, A O based equations are deriv ed based on this approac h. This w ork adds to the formal framew ork of A O based MBPT/CC theory and elucidates similarities and dierences b et w een the MO and A O based formalism. It is also p ossible to write do wn the zeroth-order Hamiltonian and excitation amplitudes solely based on A O basis. Their denitions are giv en, ho w ev er, the deriv ation of CC equations based on these denitions are highly complicated and length y b ecause w e need to start from the ph ysical v acuum state rather than the F ermi v acuum state. 3.2 Con v en tional CC o v erview Though, w e are trying to deriv e the A O-based CC equations, it is useful to review the con v en tional CC to b etter understand the connection b et w een the t w o approac hes. The coupled cluster w a v e function is written as j C C i = e ^ T j 0 i ; (3.1) where T is excitation op erator, H N is the normal ordered zeroth-order Hamiltonian and 0 is the zeroth-order w a v e function. The HF orbital is a p opular c hoice as a zeroth-order description of the system. The excitation op erator is the sum of 1-, 2-, ... up to n-electron excitation op erators ^ T = ^ T 1 + ^ T 2 + + ^ T n ; (3.2) with n b eing the n um b er of electrons in the system. The n -electron excitation op erator is dened b y ^ T n = 1 ( n !) 2 X ij::: ab::: t ab::: ij::: f ^ a y ^ b y ::: ^ j ^ i g ; (3.3) where t is a set of co ecien ts to determine, called excitation amplitudes.

PAGE 61

51 The normal ordered Hamiltonian is giv en as H N = X pq f p q f ^ p y ^ q g + 1 4 X pq r s w pq r s f ^ p y ^ q y ^ s ^ r g ; (3.4) where f and w are oneand t w o-electron in tegrals constructed b y the A O to MO transformation of the F o c k matrix and HF t w o-electron in tegrals. Note that lab els i; j; k ; l ; ::: represen t o ccupied orbitals, a; b; c; d; ::: represen t virtual orbitals, and p; q ; r ; s; ::: represen t an y orbitals b y con v en tion. The equations to calculate the correlation energy and excitation amplitudes are obtained b y pro jecting the CC w a v efunction on to zeroth-order w a v efunction and sets of excitation congurations, resp ectiv ely E cor r = h 0 j H N j 0 i ; (3.5) 0 = h m j H N j 0 i ; (3.6) where m is a set of a i and ab ij for the single-, double-excitations and higher lev el excitations are dened in a similar manner. 3.3 T o ols Bridging Bet w een the MO and A O Represen tations Some useful connections b et w een the MO and A O represen tations when w e deriv e the A O based coupled cluster equations are presen ted. In addition to the MO lab el con v en tion, greek letters are used to represen t A O lab els. 3.3.1 The Orthogonalit y Condition The orthogonalit y of the MO basis is expressed in terms of A O as pq = h p j q i ; = X ; h j ( c p ) c q j i ;

PAGE 62

52 = X c p h j i c q ; = X c p S c q : (3.7) Here, S is the o v erlap matrix and c is the MO co ecien t. Dening c q = X S c q ; (3.8) then pq = X c p c q : (3.9) Multiplying the Eq.(3.7) b y ( c p r ) from the left, b y ( c q ) from the righ t, and summing o v er p and q yields X pq ( c p r ) pq ( c q ) = X pq X ( c p r ) c p S c q ( c q ) : (3.10) Note that the righ t hand side has non-zero v alue when p = q = i or p = q = a due to the delta function. Dening the densit y matrix P and its virtual orbital and an y orbital v arian ts Q , and R as P = X i c i c i ; (3.11) Q = X a c a c a ; (3.12) R = X p c p c p ; (3.13) the ab o v e relationships ma y b e simplied as X P r S P = P r for p = q = i; (3.14) X Q r S Q = Q r for p = q = a; (3.15) X R r S R = 0 otherwise : (3.16)

PAGE 63

53 The last case includes ( p = i; q = j ), ( p = i; q = a ), and ( p = a; q = b ). If w e c ho ose to use the orthogonalized A O basis, whic h ma y b e obtained b y rotating the orbitals, then w e ha v e a simpler orthogonalit y condition pq = X c p c q : (3.17) In this case, Eqs.(3.14) to (3.16) analogs are X P r P = P r for p = q = i; (3.18) X Q r Q = Q r for p = q = a; (3.19) X R r R = 0 otherwise : (3.20) 3.3.2 A O to MO T ransformation The A O to MO transformation form ula for a general MO quan tit y X and an A O quan tit y Y is X pq ::: r s::: = X :::::: c p c q :::c r c s :::Y ::: ::: ; (3.21) Multiplying b oth sides of the ab o v e equation b y P pq r s c p c q c r c s ::: and using Eq.(3.9) yields X pq r s c p c q c r c s :::X pq ::: r s::: = X :::::: X pq r s ( c p c p )( c q c q ) ::: ( c r c r )( c s c s ) :::Y ::: ::: = Y ::: ::: : (3.22) This is the MO to A O transformation in the general form. Creation/annihilation op erators are dened in a similar manner to the A O to MO transformation. ^ i = X c i ^ ; (3.23)

PAGE 64

54 ^ i y = X c i ^ y ; (3.24) ^ a = X c a ^ ; (3.25) ^ a y = X c a ^ y : (3.26) The righ t hand form needs to b e considered only when w e w an t to eliminate all the MO indices throughout the deriv ation of A OCC equations. Otherwise, the creation/annihilation op erators are alw a ys giv en as their original form. 3.3.3 Beha vior of the A O Creation/annihilation Op erators There is a ma jor dierence b et w een the MO creation/annihilation op erators and their A O analog, that is the orthogonalit y of MO and non-orthogonalit y of A O. The con traction of a general creation and an annihilation op erator denoted b y underline o v er the t w o op erators is giv en as ^ p ^ q y = ^ p ^ q y ( ^ q y ^ p ) = [ ^ p; ^ q y ] + = pq ; (3.27) whereas in terms of A O, the delta function is replaced b y the o v erlap matrix. ^ ^ y = ^ ^ y ( ^ y ^ ) = [ ^ ; ^ y ] + = S : (3.28) When w e form ulate CC equations in terms of MO, it is more con v enien t to in tro duce the F ermi v acuum as the reference state and re-dene the creation and annihilation op erators in the so-called, hole-particle form. In the A O-based form ulation, w e ha v e no suc h distinction b et w een holes and particles.

PAGE 65

55 3.3.4 Denition of the Reference State The F ermi v acuum state is dened in the MO represen tation as j 0 i = ^ i y 1 ^ i y 2 ^ i y n j 0 i = j i 1 i 2 i n i ; (3.29) where j 0 i is the ph ysical v acuum and n is the n um b er of electrons in the system. The F ermi v acuum state cannot b e dened solely in the A O basis, b ecause w e require a F ermi lev el that distinguishes b et w een particle and holes. Hence w e use the MO to A O transformation. Here, w e review imp ortan t prop erties related to the F ermi v acuum state in the MO formalism. When w e consider the normal ordering relativ e to the F ermi v acuum in the MO form ulation, all the virtual creation and o ccupied annihilation op erators are left of all the virtual annihilation and o ccupied creation op erators. The con tractions of op erator strings relativ e to the F ermi v acuum, denoted b y o v erbars across the op erators, are ^ i y ^ j = ij ; (3.30) ^ i ^ j y = 0 ; (3.31) ^ a y ^ b = 0 ; (3.32) ^ a ^ b y = ab : (3.33) 3.4 The Ground State Coupled Cluster Equations Using the basic to ols in tro duced in previous subsections, the A O based CC equations are deriv ed. The zeroth-order Hamiltonian and n-electron excitation amplitudes in the A O basis are obtained b y rst transforming those in the MO

PAGE 66

56 basis and p erforming the summation o v er MO indices, H AO N = X X pq c p c q f f ( X c p ^ y )( X r c q ^ r ) g + 1 4 X X pq r s c p c q c r c s w f ( X c p ^ y )( X c q ^ y )( X r c s r ^ r )( X c r ^ ) g = X X r R R r f f ^ y ^ r g + 1 4 X X r R R R R r w f ^ y ^ y ^ r ^ g ; (3.34) and ^ T AO n = 1 ( n !) 2 X ::: ::: X ij::: ab::: c i c j :::c a c b :::t ::: ::: f ( X c a ^ y )( X c b ^ y ) ::: ( X c j ^ )( X o c i o ^ o ) g = 1 ( n !) 2 X ::: ::: X ::: o ::: P o P :::Q Q :::t ::: ::: f ^ y ^ y ::: ^ ^ o g : (3.35) where P , Q , and R dened in Eqs.(3.12) to (3.13) are used. In order to deriv e the CC equations with these pure A O expressions, w e need to in tro duce MO excitations and transform them to A O's as ab o v e. E C C D = h 0 j H AO N (1 + ^ T AO 2 ) j 0 i = h 0 j H AO N ^ T AO 2 j 0 i = h 0 j ^ i n ^ i n 1 ^ i 2 ^ i 1 H AO N ^ T AO 2 ^ i y 1 ^ i y 2 ^ i y n 1 ^ i y n j 0 i : (3.36) No w w e mo v e the creation op erators to the left. Though sw apping t w o adjacen t creation/annihilation op erators in tro duces a phase factor of (-1), mo ving the ^ i y 1 to the left of H AO N ^ T AO 2 giv es a phase factor of 1, since the op eration in v olv es an ev en n um b er of in terc hanges. Using the fact ^ i ^ i y = 1 ; (3.37)

PAGE 67

57 mo ving all the creation op erators to the left of H AO N ^ T AO 2 simplies the expression as E C C D = h 0 j H AO N ^ T AO 2 j 0 i : (3.38) Substituting the Eqs.(3.34) and (3.35) in to the ab o v e expression yields, E C C D = 1 4 h 0 j [ X X r R R r f f ^ y ^ r g + 1 4 X X r R R R R r w f ^ y ^ y ^ r ^ g ] [ X X o P o P Q Q t f ^ y ^ y ^ ^ o g ] j 0 i : (3.39) The rst term in the square brac k et v anishes b ecause there are t w o op erators in the f term, whereas there are four in the amplitude part, th us there is no w a y of making external con tractions for all the op erators. Wic k's theorem [72], whic h sa ys the v acuum exp ection v alue of an y normal pro duct with con tractions is zero unless all the op erators are con tracted, leads us to the conclusion that the rst term is zero. Th us, E C C D = 1 16 h 0 j [ X X r R R R R r w f ^ y ^ y ^ r ^ g ] [ X X o P o P Q Q t f ^ y ^ y ^ ^ o g ] j 0 i : (3.40) Lo oking at the con traction part, w e observ e that ^ y and ^ y m ust b e con tracting with ^ and ^ o and that ^ r and ^ m ust b e con tracting with ^ y and ^ y . Th us, there are four p ossible com binations whic h w ould giv e non-zero con tributions. T aking a particular com bination of con traction ^ y and ^ o , ^ y and ^ , ^ r and ^ y , and ^ and ^ y ,

PAGE 68

58 eac h con traction yields an o v erlap matrix b y virtue of Eq.(3.28). Th us, E C C D = 1 16 X X r X X o h 0 j [ R R R R r w ]( S o S S r S +t hr ee other s )[ P o P Q Q t ] j 0 i : (3.41) Recalling the prop erties deriv ed from the orthogonal condition, Eqs.(3.14) to (3.16), w e ha v e non-zero con tributions only when the rst t w o R 's are P 's and the later R 's are Q 's. Reorganizing terms as X o P S o P o = P ; (3.42) X P S P = P ; (3.43) X r Q r S r Q = Q ; (3.44) X Q S Q = Q ; (3.45) then, w e nally get the energy expression as E C C D = 1 16 X X ( Q Q P P + t hr ee other s ) w t = 1 4 X t ; (3.46) where = X Q Q P P w : (3.47) This is a form of A O-based energy expression. Note that w e obtain exactly the same expression for the orthogonalized A O case. Amplitude equations can b e deriv ed in a similar manner. Here w e only emphasize the diference from the energy

PAGE 69

59 expression deriv ation. The doubles amplitude equations are written as 0 = h ab ij j H N (1 + ^ T 2 + 1 2 ^ T 2 2 ) j 0 i = h 0 j X ij ab ( c i )( c j )( c a )( c b ! ) ^ i n ^ i n 1 ^ i i +1 ^ a ^ i i 1 ^ i j +1 ^ b ^ i j 1 ^ i 2 ^ i 1 H N (1 + ^ T 2 + 1 2 ^ T 2 2 ) ^ i y 1 ^ i y 2 ^ i i ^ i j ^ i y n 1 ^ i y n j 0 i : (3.48) Mo ving the creation op erators to the left lea v es four op erators uncanceled due to the double excitation. Renaming ^ i i as ^ i and ^ i j as ^ j , and expand them using Eqs.(3.23) to (3.26), w e ha v e 0 = h 0 j X ij ab ( c i )( c j )( c a )( c b ! )( X c i ^ y )( X c j ^ y )( X c a ^ )( X c b ^ ) H N (1 + ^ T 2 + 1 2 ^ T 2 2 ) j 0 i = X h 0 j P P Q Q ! f ^ y ^ y ^ ^ g H N (1 + ^ T 2 + 1 2 ^ T 2 2 ) j 0 i : (3.49) Though more complicated, the remaining part has to follo w the same pro cedure as the total energy deriv ation: substituting H N , T 2 , and T 2 2 b y the A O expressions, Eqs.(3.34) and (3.35), considering all the com binations of con tractions of A O creation and annihilation op erators, pic king up non-zero con tributions P S P and QS Q and nally reorganizing terms. There can b e man y forms of energy expressions b y com bining P 's and Q 's in dieren t w a ys and also b y inserting the resolution of the iden tit y , Eq.(3.9), in the expression. Scuseria ga v e t w o v ersions of suc h A O expressions [84]. One of them features R 3 deca y prop ert y of in tegrals as a function of separation b et w een c harge distributions. He obtained this particular form of the CC equation b y p erforming the MO to A O transformation on nal equations deriv ed in the MO basis. Here, w e deriv e the same form ula within an A O-based framew ork from inception. The complexit y comes from the facts that w e ha v e to constract P 's and Q 's b y taking

PAGE 70

60 one c from Eq.(3.34) and the other from Eq.(3.35) and that w e insert four of resolution of iden tit y t yp e terms, Eq.(3.9). Ho w ev er, after understanding the A O-based equation in the ab o v e deriv ation, an exp edien t w a y to deriv e w orking form ulas is simply tak e the nal equations form ulated in MO, and transform them or simply replace the MO lab els b y corresp onding A O lab els and lift all the o ccupied and virtual orbital restrictions in summations. 3.5 W orking F orm ulas W e ha v e already seen that there are man y similarities b et w een the MOand A O-based CC form ulations. The same pro cedure can b e applied for CC equations b ey ond the ground state energy and amplitudes suc h as for analytical deriv ativ es and EOM-CC equations. When w e ev aluate the analytical deriv ativ e expressions, w e need deriv ativ es of the Hamiltonian, deriv ativ es of the oneand t w o-electron in tegrals and deriv ativ es of the MO co ecien ts as w ell as deriv ativ es of the excitation amplitudes. In tro ducing n-equations, whic h parallels the -equation [75, 76] in the MO represen tation, reduces the n um b er of linear equations to solv e from 3 N to one [77] where acoun ts for the rst order resp onse of the molecular system to an external p erturbation. Giv en the con v erged excitation amplitudes, w e solv e the n equations. The energy deriv ativ es with resp ect to atomic displacemen ts are written b y the rstand second-order reduced densit y matrices, whic h are dened b y the excitation amplitudes and ! amplitudes, and the deriv ativ e in tegrals. The deriv ativ e in tegrals consist of deriv ativ es with resp ect to the co ecien ts (P or Q) and deriv ativ es with resp ect to the bare t w o-electron in tegrals v . The P and Q co ecien ts are written in terms of the MO co ecien ts, th us their deriv ativ es are obtained b y solving the coupled p erturb ed HF equations.

PAGE 71

61 In this section, w orking CCSD, EOM-CCSD and CCSD gradien t [78] expressions are presen ted in terms of A O using the R 3 deca ying in tegrals. All the form ulas are in a quasi-linearized form [73, 76] ensuring p erformance and ease of co ding. 3.5.1 The CCSD Ground State Equations The CCSD energy E C C S D = X x f + 1 4 X % : (3.50) The single excitation amplitudes: 1 equation Z = x f + X v F X o F + X x F X C 1 2 X H : (3.51) The double excitation amplitudes: 2 equation Z = D + P ( ) X ( v F 1 2 X x F ) P ( ) X ( o F + 1 2 X x F ) + 1 2 X % A + 1 2 X % B + P ( ) P ( ) X ( E E ) + P ( ) X I P ( ) X J : (3.52) The transformed A O-CC one-electron in tegrals o f ; v f ; and x f are o F = X Q f Q ; (3.53) v F r = X P r f P ; (3.54) x F = X P f Q : (3.55)

PAGE 72

62 The transformed A O-CC t w o-electron in tegrals A; B ; C ; D ; E ; G; H ; I ; J ; K ; L; M and are dened as X = 2 x x (3.56) ( X = A; B ; C ; D ; E ; G; H ; I ; J ; K ; L; M ; ; x = a; b; c; d; e; g ; h; i; j; k ; l ; m; ) ; and the small letter in tergrals are a = X P P v P P ; (3.57) b = X Q Q v Q Q ; (3.58) c = X P Q v P Q ; (3.59) d = X P P v Q Q ; (3.60) e = X P Q v Q P ; (3.61) g = X P Q v Q Q ; (3.62) h = X P P v Q P ; (3.63) i = X Q Q v Q P ; (3.64) j = X P Q v P P ; (3.65) k = X Q P v Q Q ; (3.66) l = X P P v P Q ; (3.67) m = X Q P v P Q ; (3.68) = X P P v Q Q : (3.69) Note that for the orthogonalized A O basis P = P and Q = Q . Th us some of the in tegrals ab o v e b ecome iden tical.

PAGE 73

63 The one-elecron in termediates o F ; v F ; and x F are o F = (1 ) o f + 1 2 X x f + X L + 1 2 X ~ % ; (3.70) v F = (1 ) v f 1 2 X x f + X G 1 2 X ~ % ; (3.71) x F = x f + X : (3.72) The t w o-electron in termediates A ; B ; and E are A = A + P ( ) X L + 1 4 X % ; (3.73) B = B P ( ) X K + 1 4 X % ; (3.74) E = E + X G X H X ( 1 2 + ) : (3.75) Eectiv e t w o-particle excitation amplitudes are dened as ~ % = + 1 2 ( ) ; (3.76) % = + : (3.77) The denominator terms are dened as Z = o f v f ; (3.78) Z = o f + o f v f v f : (3.79) 3.5.2 The CCSD Gradien t Equations The single excitation amplitudes: n 1 equation Z ! = x F + X ! v ~ F X ! o ~ F + 1 2 X ! ~ I

PAGE 74

64 + X ! ~ M 1 2 X ! ~ J X v G K X o G L + X ( X v G X o G ) : (3.80) The double excitation amplitudes: n 2 equation Z ! = D + P ( ) X ! v ~ F P ( ) X ! o ~ F + 1 2 X ! ~ B + 1 2 X ! ~ A + P ( ) P ( ) X ! ~ M + P ( ) X ( v G + X ! ) + P ( ) X ( o G + X ! ) + P ( ) P ( ) ! x F + P ( ) X ! K P ( ) X ! L : (3.81) One-electron in termediates v ~ F = v F 1 2 X x F ; (3.82) o ~ F = o F 1 2 X x F : (3.83) Tw o-electron in termediates ~ A = A + 1 4 X % ; (3.84) ~ B = B + 1 4 X % ; (3.85) ~ M = M + 1 2 X ; (3.86) ~ J = J X x F X o o ~ A o + 1 2 X G % P ( ) X ~ ~ E + P ( ) X o L o o ; (3.87)

PAGE 75

65 ~ I = I + X x F + X ~ B + 1 2 X o H o % o P ( ) X ~ ~ E + P ( ) X K ; (3.88) with the double tilde in termediate giv en as ~ ~ E = E X : (3.89) Also for computational reasons, another set of one-electron in termediates are dened as o G = 1 2 X ! ; (3.90) v G = 1 2 X ! : (3.91) When a quan tit y is dened as a function of and ! amplitudes suc h as G in termediates, the tensor notation cannot b e applied, th us all the indices are written as subscripts. Giv en b oth sets of the and ! amplitudes, the CCSD energy gradien t with resp ect to an external p erturbation is expressed as dE d = X o D @ o F @ + X v D @ v F @ + X AB C D I L X X X ( ; ) @ X @ ; (3.92) with the one-particle densit y matrices o D = 1 4 P + ( ) X ! 1 2 P + ( ) X ! ; (3.93) v D = 1 4 P + ( ) X ! + 1 2 P + ( ) X ! ; (3.94) and the t w o-particle densit y matrices A ( ; r ) = 1 8 P + ( ; r ) A r ; (3.95)

PAGE 76

66 B ( ; ) = 1 8 P + ( ; ) B ; (3.96) C ( ; ) = 1 4 P + ( ; ) C + 1 8 P + ( ; ) X ! 1 8 P + ( ; ) ! ; (3.97) D ( ; ) = 1 8 % + 1 8 ! + 1 16 X % A 1 8 P ( ) X % ( o G + X ! ) 1 8 P ( ) X % ( v G X ! ) 1 8 P ( ) P ( ) X ( + 2 )( X C ! + ! ) + 3 2 P ( ) P ( ) X ! ! ; (3.98) I ( ; ) = 1 8 X ! % + 1 8 X ! 1 8 X B + 1 4 P ( ) X C 1 8 P ( ) v G ; (3.99) L ( ; r ) = 1 8 X ! r % 1 8 X r ! + 1 8 X A r + 1 4 P ( ) X C r 1 8 P ( ) o G r ; (3.100) and in termediates A = 1 2 X % ! ; (3.101) B = 1 2 X % ! ; (3.102) C = 1 2 X ! : (3.103)

PAGE 77

CHAPTER 4 A TOMIC-ORBIT AL BASED COUPLED CLUSTER THEOR Y F OR EXTENDED SYSTEMS 4.1 In tro duction Recen t success in molecular ab initio quan tum c hemistry has triggered emerging demands for highly accurate electron correlation metho ds for large systems. Ho w ev er, b ecause of the so called \p o w er-la w dep endence" of the electron correlation metho ds, applications of highly correlated metho ds to large systems ha v e b een computationally prohibitiv e. F or example, the cost scaling for CCSD is iterativ e O ( N 6 ), and for CCSD(T) is iterativ e O ( N 6 ) plus non-iterativ e O ( N 7 ). Recen t eorts on linear scaling metho ds mak e it p ossible, in principle, to p erform correlated calculations suc h as second-order man y-b o dy p erturbation theory [MBPT(2)] [79{ 81] and coupled cluster theory (CC) [82{84] for suc h large systems. A form ulation, taking adv an tage of the lo calit y of atomic orbitals, w as giv en b y Scuseria and Ay ala [84] for the MBPT(2) and CC singles and doubles (CCSD) for large molecules. A linear scaling non-iterativ e triples correction w as giv en b y Sc h utz [83] in the follo wing y ear. Both form ulations ac hiev e linear scaling, O ( N ), at the large molecule limit with a larger prefactor than the con v en tional metho ds. There m ust b e a critical system size where the computation time b y con v en tional metho ds and b y A O-based metho ds coincide, but the critical system size is still under discussion. Extended systems, suc h as p olymers, slabs and crystals, are sp ecial cases of large systems in the sense that w e need to explicitly treat a large n um b ers of atoms, though there are distinctions from large but nite systems, and innite 67

PAGE 78

68 but non-p erio dic systems. T raditionally electron correlation in innite p erio dic systems is treated b y densit y functional theory , while more accurate results ma y b e obtained b y the CC theory . F or the treatmen t of innite p olymers, v arian ts of linear scaling metho ds are b eing dev elop ed but the common idea of the linear scaling correlated metho ds is to use the lo calit y of the electron correlation eects frequen tly imp osing some arbitrary cutos based up on pro conceiv ed c hemical lo calit y argumen ts. W e pursue other approac hes to b e able to treat an y systems without an y c hemical in tuitions. W a v e-function based electron correlation metho ds for p olymer ha v e b een pioneered b y Suhai [85{87] with his w ork on the fourthorder p erturbation theory . Sun and Bartlett [88{90] obtained the quasi-particle energy , x-ra y and ultra-violet photo-electron sp ectra b y MBPT(2), while analytical gradien ts for MBPT(2) w as dev elop ed b y Hirata and Iw ata [91]. F orner and co-w ork ers [92] form ulated and implemen ted CC doubles (CCD) and linearized CCD (LCCD). Reinhardt [93] considered coupled-electron-pair appro ximations (CEP A) from a dressed-conguration-in teraction p oin t of view, and Dolg and co-w ork ers [94{96] demonstrated the eectiv eness of the so-called incremen tal approac h for CC for extended systems. V ery recen tly , CCSD [97] has b een implemen ted for extended systems in our researc h group. The ultimate question in implemen ting highly correlated metho ds for extended systems is ho w w e ac hiev e the rapid execution and accurate results. This, in turn, raises the question, what kind of orbitals should w e c ho ose in the form ulation. W e form ulate and implemen t MBPT(2), LCCD, and CCD for onedimensional extended systems using atomic-orbital based tec hniques. Our previous implemen tation [97] is based on crystalline orbitals. These t w o form ulations are compared formally and n umerically . Comparison of the n um b er of amplitudes to determine, the deca y rate of in tegrals, and a relationship b et w een the cell summation cuto parameter and the correlation energy retriev al are sho wn. In this

PAGE 79

69 pap er, the formalism for the atomic orbital based CCD and its appro ximation, LCCD and MBPT(2) are presen ted initially . Then the in tegral transformations for innite p olymers are discussed. F ormal and n umerical comparisons b et w een the crystalline-orbital-based and the atomic-orbital-based implemen tation follo ws. The lithium h ydride c hain and tr ansp oly acet ylene are c hosen to sho w the strength and w eakness of our implemen tation compared with the con v en tional CC implemen tation for p olymers. 4.2 F ormalism In this section, the A O-CC form ulation for one-dimensional p olymers is presen ted. Extensions to t w oand three-dimensional systems are straigh tforw ard. F or the higher-dimensional systems w e ha v e to ha v e v ector quan tities k and Q instead of scalar. The notation used here closely follo ws the A O-CC form ulation for molecules b y Scuseria and Ay ala [84]. F or extended system sp ecic parts, the follo wing v e rules are applied: (I) The same subscript and sup erscript cancel when summed. (I I) The star refers to the complex conjugate and its op eration in terc hanges the subscript and sup erscript of the starred quan tit y . (I I I) The square brac k ets are used to refer to k -space cell indices, while the paren theses are used for real-space indices. (IV) p 1 ; p 2 ; and p 3 are cell indices determined b y the Nam ur cuto criterion for the bare t w o-electron in tegrals, while q 1 ; q 2 ; and q 3 are another set of parameters for the transformed A O-CC in tegrals. (V) In tegrals and amplitudes ha v e translational symmetry in the direction of the p erio dicit y . By virtue of the symmetry , w e need to kno w only N 4 Q 3 n um b ers of in tegrals and amplitudes instead of N 4 Q 4 . F or a general quan tit y Y , the translational symmetry is dened as Y ( Q 1 ) ( Q 2 ) ( Q 3 ) ( Q 4 ) = Y ( Q 1 + X ) ( Q 2 + X ) ( Q 3 + X ) ( Q 4 + X ) , where X

PAGE 80

70 can b e an y in teger. Th us w e can alw a ys c ho ose the upp er left cell indices to b e zero b y c ho osing X = Q 1, This prop ert y is used extensiv ely through the deriv ation of the A O-CC form ulas for extended systems. F ollo wing the pro cedure discussed in Chapter 3, w e can obtain man y forms of A O based expression. One of them is E cor r = 1 K X ij ab X k 1 k 2 k 3 k 4 w i [ k 1 ] j [ k 2 ] a [ k 3 ] b [ k 4 ] t a [ k 3 ] b [ k 4 ] i [ k 1 ] j [ k 2 ] = 1 K X ij ab X k 1 k 2 k 3 k 4 w i [ k 1 ] j [ k 2 ] a [ k 3 ] b [ k 4 ] t a [ k 3 ] b [ k 4 ] i [ k 1 ] j [ k 2 ] ( 1 K X [ c i [ k 1 ] ] c i [ k 1 ] ) ( 1 K X [ c j [ k 2 ] ] c j [ k 2 ] )( 1 K X c a [ k 3 ] [ c a [ k 3 ] ] )( 1 K X c b [ k 4 ] [ c b [ k 4 ] ] ) = 1 K X X k 1 k 2 k 3 k 4 w [ k 1 ] [ k 2 ] [ k 3 ] [ k 4 ] t [ k 3 ] [ k 4 ] [ k 1 ] [ k 2 ] ; (4.1) where w [ k 1 ] [ k 2 ] [ k 3 ] [ k 4 ] = 1 K 2 X ij ab c i [ k 1 ] c j [ k 2 ] [ c a [ k 3 ] ] [ c b [ k 4 ] ] w a [ k 1 ] b [ k 2 ] i [ k 3 ] j [ k 4 ] ; (4.2) t [ k 3 ] [ k 4 ] [ k 1 ] [ k 2 ] = 1 K 2 X ij ab [ c i [ k 1 ] ] [ c j [ k 2 ] ] c a [ k 3 ] c b [ k 4 ] t i [ k 3 ] j [ k 4 ] a [ k 1 ] b [ k 2 ] : (4.3) This oers a straigh tforw ard transformation of the con v en tional CC to A O-CC form ulation but do es not exploit sparsit y b ey ond the t w o-electron in tegrals transformation. W e can alternativ ely manipulate the correlation energy as follo wing using b oth the resolution of iden tit y and MO to A O transformation. E cor r = 1 K X ij ab X k 1 k 2 k 3 k 4 w i [ k 1 ] j [ k 2 ] a [ k 3 ] b [ k 4 ] t a [ k 3 ] b [ k 4 ] i [ k 1 ] j [ k 2 ]

PAGE 81

71 = 1 K X ij ab X k 1 k 2 k 3 k 4 X Q 1 Q 2 Q 3 Q 4 X r c i [ k 1 ] c j [ k 2 ] [ c a [ k 3 ] r ] [ c b [ k 4 ] ] exp ( i ( k 1 Q 1 + k 2 Q 2 k 3 Q 3 k 4 Q 4 ) a ) w [ Q 1 ] [ Q 2 ] r [ Q 3 ] [ Q 4 ] c a [ k 3 ] r c b [ k 4 ] [ c i [ k 1 ] ] [ c j [ k 2 ] ] exp ( i ( k 1 Q 1 k 2 Q 2 + k 3 Q 3 + k 4 Q 4 ) a ) t r [ Q 3 ] [ Q 4 ] [ Q 1 ] [ Q 2 ] ( 1 K X c i [ k 1 ] [ c i [ k 1 ] ] )( 1 K X c j [ k 2 ] [ c j [ k 2 ] ] )( 1 K X [ c a [ k 3 ] ] c a [ k 3 ] ) ( 1 K X [ c b [ k 4 ] ] c b [ k 4 ] ) = X X Q 1 Q 2 Q 3 Q 4 ( Q 1 ) ( Q 2 ) ( Q 3 ) ( Q 4 ) ( Q 3 ) ( Q 4 ) ( Q 1 ) ( Q 2 ) ; (4.4) where ( Q 1 ) ( Q 2 ) ( Q 3 ) ( Q 4 ) = 1 K 2 X k 1 k 2 k 3 k 4 X ij ab X r [ c i [ k 1 ] ] c i [ k 1 ] [ c j [ k 2 ] ] c j [ k 2 ] c a [ k 3 ] [ c a [ k 3 ] r ] c b [ k 4 ] [ c b [ k 4 ] ] exp ( i ( k 1 Q 1 + k 2 Q 2 k 3 Q 3 k 4 Q 4 ) a ) w ( P 1 ) ( P 2 ) r ( P 3 ) ( P 4 ) ; (4.5) ( Q 3 ) ( Q 4 ) ( Q 1 ) ( Q 2 ) = 1 K 2 X k 1 k 2 k 3 k 4 X ij ab X r c i [ k 1 ] [ c i [ k 1 ] ] c j [ k 2 ] [ c j [ k 2 ] ] [ c a [ k 3 ] ] c a [ k 3 ] r [ c b [ k 4 ] ] c b [ k 4 ] exp ( i ( k 1 Q 1 k 2 Q 2 + k 3 Q 3 + k 4 Q 4 ) a ) t r ( Q 3 ) ( Q 4 ) ( Q 1 ) ( Q 2 ) ; (4.6) t r ( Q 3 ) ( Q 4 ) ( Q 1 ) ( Q 2 ) = 1 K t r ( Q 3 ) ( Q 4 ) ( Q 1 ) ( Q 2 ) : (4.7) Note the c pro ducts compared to Eqs.(4.2) and (4.3). In this form ulation, the in tegrals and amplitudes are giv en in real-space, in con trast with the k -space form ulations for the CO-based metho d and the A O-based metho d constructed using only the resolution of iden tit y . T aking adv an tage of the translational symmetry , E cor r = X X q 1 q 2 q 3 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( q 1 ) ( q 3 + q 2 ) (0) ( q 3 ) = X X q 1 q 2 q 3 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) (0) ( q 3 + q 2 q 1 ) ( q 1 ) ( q 3 q 1 ) : (4.8) This allo ws us to dene the an ti-symmetrized in tegrals and amplitudes as (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = 2 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) (0) ( q 3 ) ( q 3 + q 2 ) ( q 1 ) ; (4.9)

PAGE 82

72 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = 2 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( q 3 ) (0) ( q 1 ) ( q 3 + q 2 ) = 2 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) (0) ( q 3 ) ( q 1 q 3 ) ( q 2 ) : (4.10) This form ulation has an adv an tage o v er the rst form ulation due to the faster deca y prop ert y of the in tegrals and amplitudes [84], th us for the rest of this rep ort w e fo cus on the second form ula. Emplo ying these transformations on eac h in tegrals and amplitudes sho wing up in the con v en tional crystal orbitals amplitude equation, w e w ould obtain the amplitude equations in the A O basis. First, w e lo ok at in tegrals. Dening the densit y and its asso ciated quan tities P (0) ( q ) = 1 K X j; k c j [ k ] [ c j [ k ] ] exp ( ik q a ) ; (4.11) Q (0) ( q ) = 1 K X b; k c b [ k ] [ c b [ k ] ] exp ( ik q a ) ; (4.12) P (0) ( q ) = 1 K X j; k c j [ k ] [ c j [ k ] ] exp ( ik q a ) ; (4.13) Q (0) ( q ) = 1 K X b; k [ c b [ k ] ] c b [ k ] exp ( ik q a ) ; (4.14) One-electron in tegrals in A O basis are v F (0) ( q 1 ) = X X q q Q (0) ( q ) f ( q ) ( q ) Q ( q ) ( q 1 ) = X X q q Q (0) ( q ) f (0) ( q q ) Q (0) ( q 1 q ) ; (4.15) o F (0) ( q 1 ) = X X q q P ( q ) ( q 1 ) f ( q ) ( q ) P (0) ( q ) = X X q q P (0) ( q 1 q ) f (0) ( q q ) P (0) ( q ) : (4.16) Similarly one of the t w o-electron in tegrals in the A O basis are ( Q 1 ) ( Q 2 ) ( Q 3 ) ( Q 4 ) = X r X P P P r P P ( Q 1 ) ( P ) P ( Q 2 ) ( P ) v ( P ) ( P ) r ( P r ) ( P ) Q r ( P r ) ( Q 3 ) Q ( P ) ( Q 4 )

PAGE 83

73 = X r X P P P r P P (0) ( P Q 1 ) P (0) ( P Q 2 ) v (0) ( P P ) r ( P r P ) ( P P ) Q r (0) ( Q 3 P r ) Q (0) ( Q 4 P ) : (4.17) Note that all the in tegrals and amplitudes dened in the curren t formalism are real quan tities, in con trast to the crystalline orbital based formalism where suc h quan tities are complex. Redening the ab o v e v ariables as q 1 = Q 3 Q 1 ; q 2 = Q 4 Q 2 ; q 3 = Q 2 Q 1 ; (4.18) p 1 = P r P ; p 2 = P P ; p 3 = P P ; p = p Q 1 ; (4.19) (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X r X p 1 p 2 p 3 p P (0) ( p ) P (0) ( q 3 + p 3 + p ) v (0) ( p 3 ) r ( p 1 ) ( p 3+ p 2) Q r (0) ( q 1 p 1 p ) Q (0) ( q 3 + q 2 p 3 p 2 p ) : (4.20) The transformation ab o v e can b e done with an N 5 ( P ; Q ) 4 computational cost b y p erforming a quarter of the transformations at a time. Let the a quarter-, a half-, and three-quarter-transformed in tegrals b e Q 1, Q 2, and Q 3, resp ectiv ely , Q 1 (0) ( p 3 ) r ( p 1 ) ( q 3 + q 2 ) = X X ( p 3 + p 2 ) v (0) ( p 3 ) r ( p 1 ) ( p 3+ p 2) Q ( p 3 + p 2 ) ( q 3 + q 2 ) = X X ( p 3 + p 2 ) v (0) ( p 3 ) r ( p 1 ) ( p 3+ p 2) Q (0) ( q 3 + q 2 p 3 p 2 ) ; (4.21) Q 2 (0) ( q 3 ) r ( p 1 ) ( q 3 + q 2 ) = X X p 3 Q 1 (0) ( p 3 ) r ( p 1 ) ( q 3+ q 2) P ( q 3 ) ( p 3 ) = X X p 3 Q 1 (0) ( p 3 ) r ( p 1 ) ( q 3+ q 2) P (0) ( p 3 q 3 ) ; (4.22) Q 3 (0) ( q 3 p ) ( q 1 ) ( q 3 + q 2 ) = X r X p 1 Q 2 (0) ( q 3 ) r ( p 1 ) ( q 3+ q 2) Q r ( p 1 ) ( q 1 ) = X r X p 1 Q 2 (0) ( q 3 ) r ( p 1 ) ( q 3+ q 2) Q r (0) ( q 1 p 1 ) ; (4.23)

PAGE 84

74 (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X X p Q 3 (0) ( q 3 p ) ( q 1 p ) ( q 3+ q 2 p ) P ( p ) (0) = X Q 3 (0) ( q 3 p ) ( q 1 p ) ( q 3+ q 2 p ) P (0) ( p ) : (4.24) When j p 3 j > C 2 , the t w o-electron in tegral v w ould b e w ell appro ximated b y the m ultip ole expansion tec hnique [23{ 25]. The bare t w o-electron in tegral is expanded using the m ultip oles as v (0) ( p 3 ) r ( p 1 ) ( p 3+ p 2) = X l X l 0 X m ( 1) l 0 + m ( l + l 0 )! 1 ( a p 3 ) l + l 0 +1 Z d r 1 (0) ( r 1 ) r l 1 P m l (cos 1 ) e im 1 ( l + j m j )! r ( p 1 ) ( r 1 ) Z d r 2 (0) ( r 2 ) r l 0 2 P m l 0 (cos 2 ) e im 2 ( l 0 + j m j )! ( p 2 ) ( r 2 ) X l X l 0 X m ( 1) l 0 + m ( l + l 0 )! 1 ( a p 3 ) l + l 0 +1 h (0) j N l ;m j r ( p 1 ) ih (0) j N l 0 ;m j ( p 2 ) i ; (4.25) where a is the unit cell length. No w the in tegral can b e giv en as (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X r X p 1 p 2 p 3 p P (0) ( p ) P (0) ( q 3 + p 3 + p ) X l X l 0 X m ( 1) l 0 + m ( l + l 0 )! ( a p 3 ) l + l 0 +1 h (0) j N l ;m j r ( p 1 ) ih (0) j N l 0 ;m j ( p 2 ) i Q r (0) ( q 1 p 1 p ) Q (0) ( q 3 + q 2 p 3 p 2 p ) : (4.26) Dening new quan tities to simplify the and the other in tegrals U (0) ( q 1 ) ( p; l ; m ) = X r X p 1 P (0) ( p ) h (0) j N l ;m j r ( p 1 ) i Q r (0) ( q 1 p 1 p ) ; (4.27) V (0) ( q 1 ) ( p; l ; m ) = X r X p 1 P (0) ( p ) h (0) j N l ;m j r ( p 1 ) i P r (0) ( q 1 p 1 p ) ; (4.28) W (0) ( q 1 ) ( p; l ; m ) = X r X p 1 Q (0) ( p ) h (0) j N l ;m j r ( p 1 ) i Q r (0) ( q 1 p 1 p ) ; (4.29) X (0) ( q 1 ) ( p; l ; m ) = X r X p 1 P (0) ( p ) h (0) j N l ;m j r ( p 1 ) i Q r (0) ( q 1 p 1 p ) ; (4.30) Y (0) ( q 1 ) ( p; l ; m ) = X r X p 1 Q (0) ( p ) h (0) j N l ;m j r ( p 1 ) i P r (0) ( q 1 p 1 p ) : (4.31)

PAGE 85

75 then, (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X p 3 p X l ;l 0 ;m ( 1) l 0 + m ( l + l 0 )! ( a p 3 ) l + l 0 +1 U (0) ( q 1 ) ( p; l ; m ) U (0) ( q 2 + p 3 ) ( p; l 0 ; m ) (4.32) Note that the summation o v er p 3 is tak en o v er only a limited range since U deca ys rapidly with rep ect to cell index. The m ultip ole expansion tec hnique reduces the in tegral transformation scaling from N 5 Q 4 to ( N 4 Q 3 LM + N 4 Q 2 L 2 M ). The short-range part still has to b e done with the N 5 ( P ; Q ) 4 algorithm b ecause the m ultip ole expansion tec hinique requires distribution of t w o w ell separated c harges. Rep eating the same pro cedure, the other transformed and t w o-electron in termediates are giv en as a (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X " X p 1 p 2 p 3 p P (0) ( q 1 p 1 p ) P " (0) ( q 3 + q 2 p 3 p 2 p ) v (0) ( p 3 ) ( p 1 ) " ( p 3+ p 2) P (0) ( p ) P (0) ( p 3 q 3 + p ) ; (4.33) b (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X " X p 1 p 2 p 3 p Q (0) ( q 1 p 1 p ) Q " (0) ( q 3 + q 2 p 3 p 2 p ) v (0) ( p 3 ) ( p 1 ) " ( p 3+ p 2) Q (0) ( p ) Q (0) ( p 3 q 3 + p ) ; (4.34) c r (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X " X p 1 p 2 p 3 p P (0) ( q 1 p 1 p ) Q " (0) ( q 3 + q 2 p 3 p 2 p ) v (0) ( p 3 ) ( p 1 ) " ( p 3+ p 2) P r (0) ( p ) Q (0) ( p 3 q 3 + p ) ; (4.35) d (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X " X p 1 p 2 p 3 p P (0) ( q 1 p 1 p ) P " (0) ( q 3 + q 2 p 3 p 2 p ) v (0) ( p 3 ) ( p 1 ) " ( p 3+ p 2) Q (0) ( p ) Q (0) ( p 3 q 3 + p ) ; (4.36) e r (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = X " X p 1 p 2 p 3 p P (0) ( q 1 p 1 p ) Q " (0) ( q 3 + q 2 p 3 p 2 p ) v (0) ( p 3 ) ( p 1 ) " ( p 3+ p 2) Q r (0) ( p ) P (0) ( p 3 q 3 + p ) ; (4.37) This is done b y the quarter b y quarter transformations as is describ ed using the in tegral. The long-range part, to o, can b e done using the m ultip ole expansion

PAGE 86

76 tec hnique where the U U term in Eq.(4.32) is replaced b y V V for a , W W for b , V W for c , X X for d , and X Y for e in tegrals. The one-electron in termediates are F (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = n (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( X 6 = X q v F (0) ( q ) ( q ) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) X 6 = X q o F ( q ) ( q 3 + q 2 ) ( q 3 ) (0) ( q ) ( q 1 ) ) = X 6 = X q v F (0) ( q ) ( q ) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) + X 6 = X q v F ( q 3 ) ( q ) ( q ) (0) ( q 3 + q 2 ) ( q 1 ) X 6 = X q o F ( q ) ( q 1 ) (0) ( q 3 ) ( q ) ( q 3 + q 2 ) X 6 = X q o F ( q ) ( q 3 + q 2 ) ( q 3 ) (0) ( q ) ( q 1 ) = X 6 = X q v F (0) ( q ) (0) ( q 3 q ) ( q 1 q ) ( q 3 + q 2 q ) + X 6 = X q v F (0) ( q q 3 ) (0) ( q ) ( q 3 + q 2 q ) ( q 1 q ) X 6 = X q o F (0) ( q 1 q ) (0) ( q 3 ) ( q ) ( q 3 + q 2 ) X 6 = X q o F (0) ( q 3 + q 2 q ) (0) ( q 3 ) ( q q 3 ) ( q 1 q 3 ) = X 6 = X q v F (0) ( q ) (0) ( q 3 q ) ( q 1 q ) ( q 3 + q 2 q ) + X 6 = X q v F (0) ( q ) (0) ( q q 3 ) ( q 2 q ) ( q q 3 + q 1 ) X 6 = X q o F (0) ( q ) (0) ( q 3 ) ( q 1 q ) ( q 3 + q 2 ) X 6 = X q o F (0) ( q ) (0) ( q 3 ) ( q 2 q ) ( q 3 + q 1 ) : (4.38) The t w o-electron in termediates are j (0) r ( q r ) ( q 1 ) ( q ) = e (0) r ( q r ) ( q 1 ) ( q ) 1 2 ( q ) ( q ) ( q 1 ) ( q ) r ( q r ) ( q ) ( q ) ( q ) + 1 2 (0) ( q ) ( q 1 ) ( q ) r ( q r ) ( q ) ( q ) ( q ) = e (0) r ( q r ) ( q 1 ) ( q ) 1 2 (0) ( q q ) ( q 1 q ) ( q q ) r (0) ( q q r ) ( q q r ) ( q q r ) + 1 2 (0) ( q ) ( q 1 ) ( q ) r (0) ( q q r ) ( q q r ) ( q q r ) ; (4.39) k r ( q r ) (0) ( q 1 ) ( q ) = k r (0) ( q r ) ( q 1 q r ) ( q q r ) = c r ( q r ) (0) ( q 1 ) ( q ) 1 2 ( q ) ( q ) ( q 1 ) ( q ) r ( q r ) ( q ) ( q ) ( q ) = c r (0) ( q r ) ( q 1 q r ) ( q q r ) 1 2 (0) ( q q ) ( q 1 q ) ( q q ) r (0) ( q q r ) ( q q r ) ( q q r ) : (4.40) Finally the A O-based doubles amplitude equation for extended systems is (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = d (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) + F (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( ) + R (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( ) o F ( q 1 ) ( q 1 ) + o F ( q 3 + q 2 ) ( q 3 + q 2 ) v F (0) (0) v F ( q 3 ) ( q 3 )

PAGE 87

77 = d (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) + F (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( ) + R (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) ( ) o F (0) (0) + o F (0) (0) v F (0) (0) v F (0) (0) ; (4.41) with R (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) = a ( q ) ( q ) ( q 1 ) ( q 3 + q 2 ) (0) ( q 3 ) ( q ) ( q ) + ( q ) ( q ) ( q 1 ) ( q 3 + q 2 ) b (0) ( q 3 ) ( q ) ( q ) + ( q ) ( q ) ( q 1 ) ( q 3 + q 2 ) r ( q r ) ( q ) ( q ) ( q ) (0) ( q 3 ) r ( q r ) ( q ) +n (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) [ (0) ( q ) r ( q r ) ( q ) r ( q r ) ( q ) ( q ) ( q ) ( q ) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) r ( q r ) ( q ) ( q 1 ) ( q ) ( q ) ( q ) r ( q r ) ( q ) (0) ( q 3 ) ( q ) ( q 3 + q 2 ) + j (0) r ( q r ) ( q 1 ) ( q ) ( q ) ( q 3 ) r ( q r ) ( q 3 + q 2 ) k r ( q r ) (0) ( q 1 ) ( q ) ( q ) ( q 3 ) r ( q r ) ( q 3 + q 2 ) k r ( q r ) ( q 3 ) ( q 1 ) ( q ) (0) ( q ) r ( q r ) ( q 3 + q 2 ) ] = a (0) ( q q ) ( q 1 q ) ( q 3 + q 2 q ) (0) ( q 3 ) ( q ) ( q ) + (0) ( q q ) ( q 1 q ) ( q 3 + q 2 q ) b (0) ( q 3 ) ( q ) ( q ) + (0) ( q q ) ( q 1 q ) ( q 3 + q 2 q ) r (0) ( q q r ) ( q q r ) ( q q r ) (0) ( q 3 ) r ( q r ) ( q ) +n (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) [ (0) ( q ) r ( q r ) ( q ) r (0) ( q q r ) ( q q r ) ( q q r ) (0) ( q 3 q ) ( q 1 q ) ( q 3 + q 2 q ) r (0) ( q q r ) ( q 1 q r ) ( q q r ) (0) ( q q r ) r ( q r q r ) ( q q r ) (0) ( q 3 ) ( q ) ( q 3 + q 2 ) + j (0) r ( q r ) ( q 1 ) ( q ) (0) ( q 3 q ) r ( q r q ) ( q 3 + q 2 q ) k r (0) ( q r ) ( q 1 q r ) ( q q r ) (0) ( q 3 q ) r ( q r q ) ( q 3 + q 2 q ) k r (0) ( q 3 q r ) ( q 1 q r ) ( q q r ) (0) ( q ) r ( q r ) ( q 3 + q 2 ) ] : (4.42) Solving for giv es the A O-CCD amplitudes for extended systems. Neglecting three-term m ultiplications from the Eq.(4.42) and in termediates therein will lead to the A O-LCCD appro ximation and R = 0 will lead to A O-MBPT(2). Also, setting all the cell-related indices to zero w ould reduce to the molecular formalism. 4.3 Analysis When w e estimate p oten tial p erformance of the con v en tional CC and A OCC metho ds, lo oking at the n um b er of non-zero in tegrals/amplitudes is a go o d measure. The n um b er of A O-based in tegrals/amplitudes is considerably larger than the con v en tional CC in tegrals/amplitudes for small systems, ho w ev er, man y of the A O-based in tegrals/amplitudes b ecome negligible as the system

PAGE 88

78 b ecomes larger. Let us in tro duce t w o terminologies \lo calit y" and \sparsit y" whic h con tribute to reducing the n um b er of non-zero in tegrals/amplitudes in the A O-based sc heme. The \lo calit y" stands for the c haracteristics that the A O-based in tegrals/amplitudes deca y faster than the con v en tional in tegrals/amplitudes as a function of c harge separation. This enables us to treat only the limited ph ysical space (n um b er of unit cells included in a calculation) to reco v er the correlation energy . The \sparsit y" stands for the c haracteristics that man y of the in tegrals/amplitudes in the lo cal space are negligible. Here, the meaning of negligible dep ends on what kind of accuracy a user w an ts. In turn, w e are able to con trol accuracy of calculations b y setting this single parameter. 4.3.1 Lo calit y In the Blo c h oribital based coupled cluster form ulation [97], the n um b er of the t-amplitudes are giv en b y ( no ) 2 ( nv ) 2 (2 k max ) 3 ; (4.43) where ( no ) is the n um b er of o ccupied crystal orbitals in a unit cell, ( nv ) is the n um b er of virtual crystal orbitals in a unit cell and (2 k max ) is the n um b er of sampling p oin ts in k -space, whic h can b e set to the same v alue as the n um b er of nearest neigh b or cells included in a calculation, c 2 . While in the A O-based form ulation, the n um b er of amplitudes is giv en b y ( no + nv ) 4 (2 c 0 1 ) 2 (2 c 0 2 ) ; (4.44) where ( c 0 1 ) and ( c 0 2 ) are the Nam ur-lik e cuto parameters for the A O-based quan tities. Considering the fact that the con v en tional CC in tegrals/amplitudes deca y as r 1 , while the A O-based in tegrals/amplitudes deca y as r 3 , implies c 0 2 << k max . c 0 1 is the cuto parameter whic h has the exp onen tial deca y prop ert y .

PAGE 89

79 Though the cuto v alue c 1 dep ends on the degree of electron densit y spread, in practice c 0 1 < c 0 2 for insurators and semi-conductors giving an o v erall lo calit y factor, whic h is the ratio of (2 c 0 1 ) 2 (2 c 0 2 ) to (2 k max ) 3 b eing less than 1. Note that in molecular calculations, there are no cell related parameters c 0 1 ; c 0 2 and k , ho w ev er the same principle can b e applied on the basis of distance b et w een the electron 1 and electron 2. 4.3.2 Sparsit y The A O-CC in tegrals/amplitudes ha v e another preferable feature, sparsit y , whic h the con v en tional CC do es not ha v e. In the con v en tional CC w e need to calculate the transformed in tegrals suc h as v I [0] J [ k 3 ] A [ k 1 ] B [ k 3 + k 2 ] for all the I ; J ; A; B ; k 1 ; k 2 ; and k 3 regardless of the system size. While in A O-CC, the in tegrals and amplitudes deca y rapidly as a function of distance and as a function of the orbital exp onen t as w ell as symmetries imp osed b y the geom try of the system. Th us, it is essen tial to compute only elemen ts whic h are ab o v e the user dened threshold for a desired accuracy to ac hiev e b etter p erformance than the con v en tional CC metho ds. 4.3.3 Implemen tation Summarizing the last t w o subsections, lo calit y is a distance-based cuto whic h reduces the innite cell summation o v er cells to a nite one, and the sparsit y is an absolute v alue based cuto whic h pic ks up only the n umerically imp ortan t in tegrals and amplitudes. After p erfoming these screening pro cesses using lo calit y and sparsit y , if there are a smaller n um b er of surviving in tegrals and amplitudes than those in the con v en tional CC metho ds, w e w ould exp ect that the A O-CC could excute faster than con v en tional CC. A prescription to do the screening follo ws. First, determine the lo calit y related cuto parameters c 0 1 and c 0 2 . c 0 1 is

PAGE 90

80 dened as the maxim um q 1 , where the largest elemen t of j d (0) (0) ( q 1 ) ( q 1 ) j is ab o v e the threshold v alue. Analogously c 0 2 is dened as the maxim um q 3 , where the largest elemen t of j d (0) ( q 3 ) (0) ( q 3 ) j is ab o v e the threshold. F or a quan tit y d (0) ( q 3 ) ( q 1 ) ( q 3 + q 2 ) , the distance b et w een and and the one b et w een and should b e treated equally , th us the condition j q 2 + q 3 q 1 j c 0 2 (4.45) ensures balance in the cuto sc heme. This is kno wn as the Nam ur cuto criterion [23{ 25]. In CC calculations, construction of R or the con traction of the CC in tegrals and amplitudes, is the ma jor computational b ottlenec k. As seen in the form ulation, the in tegrals and amplitudes are coupled in complicated manner b y cell indices as w ell as A O indices. Certainly , reducing the total n um b er of ( q 1 ; q 2 ; q 3 ) com binations greatly con tributes to the sp eed-up. Th us, at the CC in tegrals transformation stage, when only a small fraction of in tegrals is ab o v e threshold for a giv en ( q 1 ; q 2 ; q 3 ), w e can c ho ose to neglect all of the in tegrals for the ( q 1 ; q 2 ; q 3 ). All asp ects of the lo calit y and sparsit y will b e revisited in a later section including the in v estigation of relationship b et w een the cuto parameters and n umerical results. 4.3.4 Resource Managemen t The prededing discussion ab out the lo calit y and sparsit y is rather formal, ho w ev er w e need more practical lev el discussions to actually p erform the calculations. In this subsection, ho w w e ac hiev e the most ecien t use of a giv en disk space and memory is discussed. In order to p erform A O-CC calculations for p olymers, disk and memory requiremen ts are quite sev ere, th us w e are frequen tly forced to pa y more CPU time as a trade o to reduce suc h requiremen ts except for the simplest systems.

PAGE 91

81 Let us rst discuss the data storage. F or the A O-CC in tegrals it is essen tial to store only elemen ts whic h are ab o v e the threshold v alue, otherwise w e ha v e to store a h uge n um b er of n umerically negligible quan tities whic h will exhaust disk space. Therefore, at eac h quarter of the in tegral transformation step, only the v alues ab o v e the threshold are stored on disk as a sequen tial le. If needed, the les are compressed b y a system call after the writing is complete. Reading a compressed le tak es longer than reading an uncompressed le for the same amoun t of data, ho w ev er, w e could exp ect ab out t wice to ten times the compression ratio, b y exp erience, dep ending on the con ten t of a le. The in tegrals a; b; c; d; and e app ear only once in the equation, th us they can b e pre-sorted in the order of the maxim um v ectorizing eciency . In tegrals d can b e calculated on the ry if w e w an t to a v oid data storage completely , since w e need the in tegral only once p er iteration. Ho w ev er, other quan tities are needed more than once in one CC iteration, th us the direct ev aluation of those w ould b e enormously costly . Also k eeping the half-transformed in tegrals on disk and p erforming the other half on the ry w ould roughly sa v e 50 p ercen t of in tegral storage space compared with the fully transformed in tegrals. In tegrals a and b , also in termediate j and k ( e and c in the case of LCCD) app ear in the equation more than once when the p erm utation is explicitly considered. Th us, there is no uniquely optim um w a y of sorting these in tegrals. If the disk space p ermits, w e could mak e a cop y of the in tegrals sorted in dieren t orders to ac hiev e the b est p erformance for eac h term, otherwise the ev aluation of some terms requires m ultiple readings of an in tegral le. In relativ ely small molecular CC calculations, it is often p ossible to store a set of iden tical in tegrals in dieren t orders so that all the con traction could b e w ell v ectorized. In other w ords, b etter CPU timing is obtained for the molecular calculations at the cost of additional disk space. In p olymer CC, the situation is often opp osite. The

PAGE 92

82 disk space requiremen t is a b ottlenec k, th us w e ha v e to pa y more CPU time to manage to run calculations within the giv en disk space. amplitudes sho w up ev erywhere in the equation and they need to b e read in man y dieren t orders. Th us, b eing able to k eep all the amplitudes on memory w ould b e ideal, but otherwise w e k eep the most up-to-date set of theta amplitudes in a direct access le indexed b y a set of cell indices ( q 1 ; q 2 ; q 3 ). If needed, again, w e could compress the le, though the dra wbac k of uncompressing the theta le ev ery time they are read could b e substan tial. When the N 4 dimension of the matrix is dense for a giv en ( q 1 ; q 2 ; q 3 ), then storing in a direct access le is the b est w a y . It w ould b e more rexible to k eep the dense parts of the amplitudes in a direct access le and sparse parts of theta as an indexed sequen tial le. A question is where is the p erformance crossing p oin t b et w een dense and sparse. A large p ortion of the disk storage comes from sets of amplitudes and their elemen t b y elemen t dierences storage needed for the DI IS [98] amplitude con v ergence acceleration. These les are needed only when the DI IS extrap olation is p erformed, th us they can b e k ept compressed except when an writing op eration is in action. Since w e need at most t w o les to b e op ened at the same time, it is p ossible to ha v e all the other les compressed to sa v e disk storage. Whether w e should store the sets of amplitudes in direct access les or sequen tial les is another question. Though, the size of an uncompressed le is m uc h larger than the size of a direct access le, b ecause w e store zero elemen ts in the le, w e exp ect a higher compression eciency since man y parts of the les can b e simply sequences of zeros. It is our exp erience that compressed sequen tial les and compressed direct access les do not dier so m uc h in size for the same n um b er of non-zero elemen ts, ev en though the sizes of uncompressed les are substan tially dieren t. Therefore, if w e ha v e enough disk space a v ailable to hold t w o uncompressed sets of direct access les for the DI IS op eration, this w ould b e the c hoice since

PAGE 93

83 w orking on a direct access le is a fully v ectorized pro cess. If w e use sequen tial les, since w e store only elemen ts whic h are ab o v e the threshold, the step to calculate a dierence of elemen ts in the n -th and n + 1-th iterations and the step to calculate extrap olated requires a one-b y-one index matc hing ev en though w e store the amplitudes in a sp ecic order, b ecause, in general, the n um b er of surviving amplitudes c hanges in eac h iteration. Also w e need to ha v e double sets of dierence amplitudes les if w e c ho ose to use the sequen tial les b ecause w e need to kno w a dieren t part of the same le at the same time at the DI IS matrix construction part. Secondly , the amoun t of a v ailable memory also con tributes to the o v erall p erformance of the co de. The in tegral transformation part w ould b e ecien tly executed if N 4 Q 2 size of memory is a v ailable since eac h of the transformation steps giv en in Eq.(4.21)-(4.24) shares one q (or p ) index in common on the left and righ t hand side of the equation. Similarly , one of the -amplitude con traction terms a will run with maxim um p erformance with N 4 Q 2 since q 3 in R and is in common, while to ev aluate the term needs N 4 Q 3 order of memory for the b est p erformance. In realit y , if there is only N 4 memory a v ailable, these terms ha v e to b e ev aluated part b y part b y reading the same in tegrals m ultiple times. Th us, applying a dieren t algorithm based on the memory requiremen t and a v ailable memory w ould lead to the b est p erformance within a giv en computer resource. Considering those factors whic h con trol the actual run time of calculations, the curren t implemen tation tak es a v ailable disk space and memory as input and determines the optim um algorithm, taking full adv an tage of the maxim um computational resource. As for memory , N 4 Q 2 and N 4 v ersions are implemen ted for b oth the in tegral transformations and con traction step. CPU time requiremen ts for b oth algorithms are roughly inden tical, ho w ev er N 4 requires Q 2 times more disk reading op erations. F or disk space managemen t, w e rst p erform the in tegral

PAGE 94

84 transformations. Based on the size of the transformed in tegral les, the size of the in termediates and DI IS related les are estimated. Based on the actual disk space a v ailable and the estimated total amoun t of uncompressed les, les can b e c hosen to b e compressed to satisfy the requiremen t. The order of priorities for the le compression is the follo wing. First of all, the DI IS related les are the targets of compression since they are used only when the DI IS extrap olation is p erformed. If compressing the DI IS les do es not pro vide enough disk space the in tegral le d is compressed and then c , e , a , b , , j , k in this order. The bare t w o-electron in tegral le can b e deleted if the A O-CC energy calculation is the nal ob jectiv e of the calculation. Ho w ev er, if one p erforms an excited state calculation or analytical gradien t calculation w e still need to k eep the t w o-electron in tegral le. Also, the an ti-symmetrized amplitudes sho w up sev eral times in the equation, this can b e constracted on the ry using the original amplitudes, . Of course, the amplitude le itself can b e compressed, ho w ev er, this is not recommended since sho ws up ev erywhere in the equation, and th us compressing this le leads to a signican t slo wdo wn. The nal disk space gain could b e ac hiev ed b y the h ybrid storage of in to a direct access le and a sequen tial le based on the sparsit y as a function of set ( q 1 ; q 2 ; q 3 ). It ma y b e ideal if it is p ossible to mak e the CC calculation quic kly con v erged without using DI IS in order to sa v e disk space, ho w ev er DI IS seems an essen tial to ol to mak e the calculation con v erge. Within all the foregoing computational constrain ts, the p erformance of the calculations is fully determined b y ho w fast w e can con tract the in tegrals and amplitudes. Since w e k eep only the in tegrals ab o v e a user-dened threshold, lo calit y and sparsit y are already tak en in to accoun t for the in tegrals. The amplitudes are read in as a N 4 size of sup er-matrix including zeros when they are stored in a direct access le, and read one b y one for the requested cell and A O indices when they are stored in a sequen tial le. There m ust b e a crossing

PAGE 95

85 p oin t in whic h the con traction b ecomes faster for one b y one m ultiplication of a corresp onding amplitude and an in tegral than the v ectorized ro w, whic h includes man y zeros times the in tegral. In v estigating this crossing p oin ts leads to fast computation of sparse matrix m ultiplication. 4.4 Numerical Results The A O-CC is implemen ted in the quan tum c hemistry co de POL YMER 2.0 [37]. T r ansp oly acet ylene and the LiH c hain are c hosen to see the v arious asp ects of the A O-CC metho d. Standard basis sets STO-3G [99],6-31G [100],631G** [101] are used in this w ork. W e use the B3L YP/6-31G* optimized geometry [7] for all the calculations presen ted, where the geometrical parameters are C = C 1.369 A, C C 1.426 A, C H 1.091 A, 6 C C C 124.5 degrees and 6 C C H 118.3 degrees. The T able 4.1 sho ws the n um b er of doubles amplitudes in v olv ed in p oly acet ylene c hain calculations with v arious basis sets in the A O-based and CO-based sc hemes. Some data sho wn in paren theses are rough estimates b y analogy with the trend observ ed in the h ydrogen and LiH c hain (not presen ted). The actual n um b er of the surviving amplitudes will b e kno wn only at run time in the A O sc heme. The CO-based metho d, on the other hand, pro vides the exact n um b er of amplitudes once system and basis set are c hosen due to the delo calized nature of CO. In v oking the dropp ed core appro ximation [102], where the excitations from the lo w energy o ccupied MOs to high energy virtual MOs are neglected, yields greater eciency as is w ell kno wn. The n um b er of amplitudes in the A O-based sc heme is con trolled b y the amplitude/in tegral cuto parameter. The n um b er of surviving amplitudes increase v ery rapidly as the cuto parameter is c hanged from 10 4 to 10 5 and 10 6 . Comparing the n um b er of amplitudes in the A O-based sc heme and CO-based sc heme for the amplitude cuto range c hosen, the A O sc heme has

PAGE 96

86 few er amplitudes to determine whic h indicates a computational adv an tage in the A O sc heme o v er the CO sc heme. Ho w ev er, comparing the A O-based result with the 10 6 cuto and the CO-based result with the dropp ed core appro ximation, the n um b er of amplitudes increases more rapidly in the A O-based sc heme than the CO-based sc heme as w e increase the basis set size. With the 6-31G** basis set, the n um b er of amplitudes in the A O sc heme is larger than that of the CO sc heme. This is due to the ratio of o ccupied orbitals and virtual orbitals. In order to simplify the discussion, let us consider a molecular analog of Eq.(4.43) and Eq.(4.44), where the ratio of the n um b er of amplitudes in the A O sc heme to the MO sc heme is giv en b y ratio = ( no + nv ) 4 no 2 nv 2 f sp : (4.46) F or (p oly)acet ylene case, no = 7 and nv = 5 ; 15 ; 27 ; and 33 for STO-3G, 6-31G, 6-31G* and 6-31G** basis set, resp ectiv ely . f sp has a strong correlation with the cuto parameter. The constan t part of the ratio ( no + nv ) 4 no 2 nv 2 tak es the smallest n um b er, 16, when no = nv , and as the dierence b et w een the no and nv b ecomes larger, the constan t part also b ecomes larger. F or this system, the constan t part for the four basis sets are 16.9, 21.2, 37.4 and 48.0, resp ectiv ely . This is the reason the n um b er of A O-based amplitudes increases more rapidly with resp ect to the basis set size than the n um b er of CO-based amplitudes. This fact is somewhat unfortunate b ecause an STO-3G basis set usually giv e the b est ratio, though w e w ould prefer to use larger basis sets. Ho w ev er, w e ha v e con trol o v er the n um b er of amplitudes, whic h will aect the timing and accuracy of calculations. This is a great feature of the A O-based sc heme. T able 4.2 presen ts the Hartree-F o c k energy and correlation energy p er unit cell calculated b y the con v en tional MBPT(2) metho d as a function of lattice summation cuto for p oly acet ylene using the STO-3G basis set. Since p oly acet ylene

PAGE 97

87 do es not ha v e a p ermanen t dip ole, the long-range eect is not included in the calculation. The Hartree-F o c k energy con v erges v ery rapidly b elo w 10 5 Hartree, while the correlation energy sho ws a slo w er con v ergence b eha vior and the fth decimal deviates as w e increase the lattice sum cuto parameter. In T able 4.3, w e sho w the calculated energy using A O-MBPT(2) for the same system. The Hartree-F o c k calculations are p erformed with the t w o-electron in tegral cuto and the k -p oin t sampling parameter b eing 12 whose energy is guaran teed to b e con v erged to the fth decimal as sho wn in T able 4.2. The correlation energy p er unit cell is ev aluated as a function of amplitude/in tegral cuto threshold. The rst thing to note is that the calculated correlation energy is accurate up to ten times the amplitude/in tegral cuto threshold. This relationship w as observ ed in preliminary computations of the correlation energy for the h ydrogen c hain as w ell. Though more cases should b e in v estigated for the general relationship b et w een the accuracy of the correlation energy and amplitude/in tegral cuto, this is p erhaps applicable, more or less, b ecause the correlation energy is a pro duct of in tegrals and amplitudes, so the error in the energy is caused b y the error in the in tegrals and amplitudes. In the presen t case, the accuracy is sho wn to b e ten times of the cuto. The second observ ation is that the A O-based form ulation leads to faster con v ergence in the correlation energy . Comparing T ables 4.2 and 4.3, the correlation energy calculated b y the con v en tional MBPT(2) deviates at the fth decimal, whereas the energy b y A O-MBPT(2) deviates at the sixth decimal. This is due to the faster deca ying prop erties of the in tegrals denition in our curren t formalism, and faster con v ergence o v er the con v en tional metho d should o ccur univ ersally regardless of the system. Setting the amplitude/in tegral cuto parameter determines the lattice summation cuto C2' in ternally . Note that in the A O-based formalism, the lattice sum for Hartree-F o c k and that for a correlated calculation are t w o separate quan tities. F rom a dieren t viewp oin t, lo oking at Eq.(4.20), the

PAGE 98

88 left hand side is the in tegrals used in the A O-CC calculation and has the cell lab el q , while the bare t w o-electron in tegrals v app ear on the righ t hand size and ha v e the cell lab el p . Th us, they are t w o separated quan tities. Ho w ev er, it is b etter to dene the cuto based on the absolute v alue of amplitude/in tegral rather than the cell cuto. The n um b er of cells, of course, c hanges as a function of system size, ho w ev er cuto criteria based on the absolute v alue of amplitude/in tegral can b e applied to systems of an y size and, as discussed, w e w ould ha v e b etter con trol o v er the accuracy of the correlation energy . The separation of the lattice sum in Hartree-F o c k and correlated calculations leads to an in teresting p oin t. W e can roughly estimate the correlation energy within a v ery short time, whic h w as not p ossible in the con v en tional CC metho d. In the case of p oly acet ylene with the STO-3G basis set, there are enormous dierences for the n um b er of cells included b et w een the cuto parameter 10 3 result and 10 4 result. The result with the cuto parameter A=3 is esp ecially striking. The corresp onding C2' for A=3 is only 2, implying the cell summation calculations w ould b e m uc h faster than the same calculation with A=4, where C2' is 12. The faster con v ergence b eha vior in the A O-based sc heme than CO-based sc heme can b e seen from a more fundamen tal p oin t of view. Figure 4.1 sho ws the deca y b eha vior of the A O-CC in tegral d and bare t w o-electron in tegral, whic h is used to construct the CO based CC in tegrals. The plot is a maxim um for the in tegral v alues with resp ect to distance of cell separation. W e observ e that the A O-CC in tegrals deca y more rapidly , and also they ha v e m uc h smaller absolute v alues. The rapid deca y of the A O-CC in tegrals guaran tees the faster con v ergence b eha vior of the A O sc heme with resp ect to the n um b er of cells included. Figure 4.2 and 4.3 sho ws the timing information obtained for the LiH c hain with v arying basis sets and cuto parameters for A O-MBPT(2) and A O-LCCD.

PAGE 99

89 The timing is obtained on a IBM RS/6000 42P mo del 270, 375MHz and 1GB of memory . In the MBPT(2) case, the 10 4 cuto and 10 6 cuto yield ab out an order of magnitude dierence in the computation time, while in LCCD the dierence is as m uc h as t w o orders of magnitudes. MBPT(2) calculations seem fast enough since the largest calculation with 6-31G** tak es ab out 100 second p er iteration. LCCD calculations are more demanding. Since w e store the in tegral j (= c ) and k (= e ) only in one w a y , not all the terms gain the maxim um v ectorized eciency . Ho w ev er, sorting the in tegrals les for all the p erm utations to obtain the maxim um eciency requires more disk space and is prohibitiv e. This resource constrain t w orsens a dierence b et w een the eciency of A O-MBPT(2) and A OLCCD compared to what it should b e. Better lo calization than the presen t A O-based quan tities, and b etter in tegralamplitude con traction sc hemes will b e k eys to ac hiev e b etter p erformance as w ell as exploiting p oin t and space group symmetries for future w ork. 4.5 Summary and F uture W ork W e ha v e implemen ted A O-based MBPT(2) and CCD for innite p erio dic 1D systems. The strengths and w eaknesses are compared b oth formally and n umerically with the CO-based metho ds. The A O-based sc heme has strong correlation among the in tegrals/amplitudes cuto parameter, accuracy of the result and execution time. F or the mE h target accuracy , whic h ma y b e obtained b y setting the cuto to 10 3 to 10 4 , the n um b er of amplitudes in the A O sc heme is exp ected to b e m uc h smaller than those in the CO sc heme. By reducing the n um b er of imp ortan t in tegrals/amplitudes in the A O-based sc heme o v er the CObased sc heme, calculations whic h w ere to o large to run are no w p ossible to run. Also, a v ery go o d estimate of the correlation energy can b e p erformed in shorter time. Applications of A O-CC and A O-MBPT(2) to mo del o xygen p olymers,

PAGE 100

90 p oly acet ylene, p olymethineimine and in v estigation the energetics and innite order eects o v er MBPT(2) on these p olymers are planed in the future.

PAGE 101

91 T able 4.1: The n um b er of doubles excitation amplitudes to determine for p olyacet ylene calculations in A O-CC and con v en tional CC. CO based amplitudes are determined based on 20 k-p oin ts. Unit : millions. A O Based CO Based 10 4 10 5 10 6 F ull Dropp ed core STO-3G 0.1 0.5 2 19 5 6-31G 1.1 ( 7 ) ( 25) 88 34 6-31G* (4 ) (25 ) (100) 286 125 6-31G** (7 ) (45 ) (200) 427 192 T able 4.2: The MBPT(2)/STO-3G energy for the p oly acet ylene c hain calculated b y the con v en tional metho d as a function of the lattice sum cuto parameter. C2=K E H F E M B P T (2) E T O T AL 6 -75.9442437 -0.1210051 -76.0652488 8 -75.9442973 -0.1209820 -76.0652793 10 -75.9443056 -0.1209730 -76.0652787 12 -75.9443069 -0.1209690 -76.0652748 16 -75.9443071 -0.1209628 -76.0652698 20 -75.9443071 -0.1209623 -76.0652674 40 -75.9443071 -0.1209570 -76.0652641 80 -75.9443072 -0.1209561 -76.0652633 T able 4.3: The A O-MBPT(2)/STO-3G energy p er unit cell for the p oly acet ylene c hain calculated as a function of the in tegrals and amplitude threshold max = 10 A . The cell summation truncation and k -p oin ts for the HF part are xed at 12. E H F = 75 : 9443069 E h A C2' E M B P T (2) E T O T AL 3 2 -0.12 -76.06 4 12 -0.120 -76.064 5 13 -0.1209 -76.0652 6 16 -0.12092 -76.06523 7 18 -0.120927 -76.065234 8 20 -0.1209275 -76.0652345

PAGE 102

92 Integrals decay behavior Unit cellMaximum integral value 0 2 4 6 8 10 10 -4 10 -3 10 -2 10 -1 10 0 10 1 AO-CC integral Two-electron integral Figure 4.1: The deca y b ehaiv or of the maxim um v alues of the regular t w o-electron in tegrals and A O-CC transformed in tegrals.

PAGE 103

93 AO-MBPT(2) LiH chain calculation timingTime per iteration (sec) 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -4 10 -5 10 -6 STO-3G 6-31G 6-31G* 6-31G** Figure 4.2: Timing p er iteration for the LiH c hain A O-MBPT(2) calculations with v arious cuto parameters and basis sets. AO-LCCD LiH chain calculation timingTime per iteration (sec) 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 -4 10 -5 10 -6 STO-3G 6-31G 6-31G* 6-31G** Figure 4.3: Timing p er iteration for the LiH c hain A O-LCCD calculations with v arious cuto parameters and basis sets.

PAGE 104

APPENDIX A ABBREVIA TIONS AND SYMBOLS COMMONL Y USED THR OUGH THE DISSER T A TION ; ; ; ; etc: A tomic orbital lab els. i; j; k ; l ; etc: Occupied molecular orbital lab els. a; b; c; d; etc: Virtual molecular orbital lab els. ; Spin lab els. a; b; c T ranslational lengths of a unit cell. (Distinction b et w een these and the virtual orbital lab els should b e clear) ; ; r Angles of a unit cell. (Distinction b et w een these and the spin lab els should b e clear) ( q ) ; ( p ) Real space cell lab els. The dimension of the ( q ) is 0 for molecule, 1 for p olymers, 2 for slabs and 3 for crystals. Eac h comp onen t is giv en as ( q ) = f q x ; q y ; q z g . [ k ] Recipro cal space sampling p oin t lab els. The dimension of the [ k ] is 0 for molecule, 1 for p olymer, 2 for slab and 3 for crystal. Eac h comp onen t is giv en as [ k ] = f k x ; k y ; k z g . Q A general direction x; y ; or z . A A atomic lab el. S ( 0 ) ( q ) An o v erlap matrix elemen t in real space. F ( 0 ) ( q ) A F o c k matrix elemen t in real space. T ( 0 ) ( q ) A kinetic matrix elemen t in real space. N ( 0 ) ( q ) A Nuclear-electron attraction matrix elemen t in real space. C ( 0 ) ( q ) A Coulom b matrix elemen t in real space. 94

PAGE 105

95 X ( 0 ) ( q ) A HF exc hange matrix elemen t in real space. Y ( 0 ) ( q ) A DFT exc hange-correlation matrix elemen t in real space. M ( 0 ) ( q ) A m ultiple-expansion or fast-m ultip ole correction matrix elemen t in real space. W ( 0 ) ( q ) An energy w eigh ted densit y matrix elemen t in real space. P ( 0 ) ( q ) A densit y matrix elemen t in real space. " ( 0 ) i [ k ] A crystalline orbital energy in recipro cal space. c ( 0 ) i [ k ] A crystalline oribital co ecien t in recipro cal space. l ; m The orbital and magnetic quan tum n um b ers. N ( 0 ) ( q ) ( l ; m ) An electronic m ultip ole matrix elemen t. N nuc ( l ; m ) A n uclear m ultip ole matrix elemen t. N tot ( l ; m ) A total m ultip ole matrix elemen t. f The exc hange-correlation functional. v xc The exc hange-correlation p oten tial. g xc The exc hange-correlation p oten tial gradien ts. m 1 ; m 2 Mixing co ecien ts b et w een HF exc hange and DFT exc hangecorrelation. ( m 1 + m 2 = 1). An atomic basis function. r The p osition of an electron. R The p osition of a n ucleus. T ; U The Nam ur cuto parameters. C The co ecien ts used in MPE and FMM. p Grid w eigh t used in densit y functional calculations. O ( n ) Scaling factor b eing n . T The excitation op erator in MO and CO basis. The excitation op erator in A O basis. t Excitation amplitudes in MO and CO basis.

PAGE 106

96 Excitation amplitudes in A O basis. D Cuto radius of the long-range Coulom b term. N The n um b er of basis functions in a system.

PAGE 107

APPENDIX B ANAL YTICAL GRADIENTS WITH RESPECT TO THE THREE DIMENSIONAL UNIT CELL P ARAMETERS The analytical energy gradien t expressions with resp ect to the unit cell parameters for the most general 3D lattice, triclinic, (a 6 = b 6 = c, 6 = 6 = r ) are deriv ed in terms of those with resp ect to eac h atom in x , y and z directions. The form ulae mak e it p ossible to optimize an y kind of p erio dic system geometries. Deriv ativ e expressions for higher-symmetry 3D systems as w ell as 1D and 2D form ulae can b e extracted from the presen ted form ulae. Here, again, ~ a is assumed on x -axis, ~ b is on x y plane. Angle is formed b y ~ b and ~ c , and angle is formed b y ~ c and ~ a . The form ulae are based on the righ t-handed system. The co ordinates of the origin of the cell ( q x ; q y ; q z ) are Q A ( q x ;q y ;q z ) X = q x a + q y b cos r + q z c cos + Q A (0 ; 0 ; 0) X ; Q A ( q x ;q y ;q z ) Y = q y b sin r + q z c cos cos cos r sin r + Q A (0 ; 0 ; 0) Y ; Q A ( q x ;q y ;q z ) Z = q z c 1 sin r q cos 2 cos 2 + sin 2 r + 2 cos cos cos r + Q A (0 ; 0 ; 0) Z ; Using the c hain rule, @ E @ a = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ a @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ a @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ a @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q x @ E @ Q A ( q x ;q y ;q z ) X ; 97

PAGE 108

98 @ E @ b = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ b @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ b @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ b @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q y cos r @ E @ Q A ( q x ;q y ;q z ) X + q y sin r @ E @ Q A ( q x ;q y ;q z ) Y ! ; @ E @ c = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ c @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ c @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ c @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q z cos @ E @ Q A ( q x ;q y ;q z ) X + q z cos cos cos r sin r @ E @ Q A ( q x ;q y ;q z ) Y + q z 1 sin r q cos 2 cos 2 + sin 2 r + 2 cos cos cos r @ E @ Q A ( q x ;q y ;q z ) Z ! ; @ E @ = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q z c sin sin r @ E @ Q A ( q x ;q y ;q z ) Y +2 q z c sin sin r cos cos cos r p cos 2 cos 2 + sin 2 r + 2 cos cos cos r @ E @ Q A ( q x ;q y ;q z ) Z ! ; @ E @ = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q z c sin @ E @ Q A ( q x ;q y ;q z ) X + q z c sin sin r cos r @ @ E @ Q A ( q x ;q y ;q z ) Y +2 q z c sin sin r cos cos cos r p cos 2 cos 2 + sin 2 r + 2 cos cos cos r @ E @ Q A ( q x ;q y ;q z ) Z ! ;

PAGE 109

99 @ E @ r = X A X q x X q y X q z @ Q A ( q x ;q y ;q z ) X @ r @ E @ Q A ( q x ;q y ;q z ) X + @ Q A ( q x ;q y ;q z ) Y @ r @ E @ Q A ( q x ;q y ;q z ) Y + @ Q A ( q x ;q y ;q z ) Z @ r @ E @ Q A ( q x ;q y ;q z ) Z ! = X A X q x X q y X q z q y b sin r @ E @ Q A ( q x ;q y ;q z ) X + ( q y b cos + q z c cos sin 2 r ) @ E @ Q A ( q x ;q y ;q z ) Y + q z c ( cos r cos cos sin r p cos 2 cos 2 + sin 2 r + 2 cos cos cos r q cos 2 cos 2 + sin 2 r + 2 cos cos cos r cos r sin 2 r @ E @ Q A ( q x ;q y ;q z ) Z ! :

PAGE 110

APPENDIX C CONNECTION AMONG THE NORMALIZED SPHERICAL HARMONICS, MUL TIPOLE MOMENTS AND SOLID SPHERICAL MOMENTS The normalized spherical harmonics, m ultip ole momen ts and solid spherical harmonics are denoted as Y m l , M m l and N m l . They are dened in terms of the asso ciated Legendre function P m l , Y m l = i j m j m s (2 l + 1) 4 ( l j m j )! ( l + j m j )! P m l (cos ) e im ; (C.1) Y m l = ( 1) m Y m l ; (C.2) M m l = r l P m l (cos ) e im ; (C.3) M m l = M m l ; (C.4) N m l = r l ( l + j m j )! P m l (cos ) e im ; (C.5) N m l = N m l : (C.6) They are in terconnected b y the follo wing form ula. r l Y m l = ( 1) m r ( l + j m j )!( l j m j )!(2 l + 1) 4 N m l ; (C.7) and M m l = ( 1) m s 4 (2 l + 1) ( l + j m j )! ( l j m j )! r l Y m l : (C.8) Multip ole in tegrals and n uclear m ultip ole momen ts are calculated in terms of N rst, then they are con v erted in to Y to p erform the axis rotation. Another transformation from Y to M follo ws since the t w o-cen ter m ultip ole expansion is 100

PAGE 111

101 form ulated in terms of M . Rotation of normalized spherical harmonics are dened b y Y m l ( ; ) = l X m = l D l ( R ) nm Y n l ( ; ) (C.9) ! [ D l ( R )] T Y n l ( ; ) ; (C.10) Y n l ( ; ) = [ D l ( R )] T 1 Y m l ( ; ) : (C.11) where D l ( R ) nm = e inr d l ( ) nm e im ; (C.12) d l ( ) n;m = ( 1) m s ( l n )!( l m )! ( l + n )!( l + m )! l X s =1 ( 1) s ( l + s )! ( l s )!( s m )!( s n )! sin 2 s n m 2 cos n + m 2 for n > 0 ; m > 0 ; d l ( ) n; m = ( 1) l + n d l ( )) n;m ; (C.13) d l ( ) n;m = ( 1) l + m d l ( )) n;m ; (C.14) d l ( ) n; m = ( 1) n + m d l ( )) n;m : (C.15) When n and/or m is zero, D l ( R ) n; 0 = ( 1) n r 4 2 l + 1 Y n l ( ; r ) ; (C.16) D l ( R ) 0 ;m = r 4 2 l + 1 Y m l ( ; ) : (C.17) A new co ordinate system is dened b y the old co ordinate system and angle ; and r . First, the original co ordinate system is rotated b y an angle , ab out the z -axis in the coun ter-clo c k wise direction (so that the x -axis mo v es to w ard the y -axis). The rotated co ordinate system is then tipp ed b y an angle ab out the new y-axis so that the z -axis mo v es in the x z plane in the direction corresp onding to the p ositiv e x -direction. The tipp ed co ordinate system is no w rotated b y an angle r , ab out the new z -axis in the coun terclo c kwise direction.

PAGE 112

102 In the t w o-dimensional extended systems with the x y plane p erio dicit y , m ultip ole-m ultip ole in teractions are calculated as follo wing. First, compute m ultip ole momen ts in the original co ordinate system (Fig. 2.5-a). Since in teraction b et w een m ultip oles are dened when t w o m ultip oles are aligned on the common z -axis, co ordinate systems m ust b e rotated. First rotation is ab out the z -axis, so that w e will ha v e a common x -axis (Fig. 2.5-b). Then the co ordinate system is tipp ed 90 degrees ab out the new y -axis. This mak es the z -axis common to b oth m ultip oles (Fig. 2.5-c). = 90 degrees and r = 0 ; 180 degrees ( e imr is the unit matrix) constrain ts come from the 2D p erio dicit y .

PAGE 113

REFERENCES [1] S. Bratoz, Collo q. In t. Cen tre Natl. Rec h. Sci. P aris 82, 287 (1958). [2] D. M. Bishop and M. Randic, J. Chem. Ph ys. 44, 2480 (1966). [3] R. Mo ccia, Chem. Ph ys. Lett. 5, 260 (1970). [4] H. T eramae, T. Y amab e, C. Satok o, and A. Imam ura, Chem. Ph ys. Lett. 101, 149 (1983). [5] H. T eramae, T. Y amab e, and A. Imam ura, J. Chem. Ph ys. 81, 3564 (1984). [6] S. Hirata and S. Iw ata, J. Ph ys. Chem. A 102, 8426 (1998). [7] S. Hirata and S. Iw ata, J. Chem. Ph ys. 107, 10075 (1997). [8] M. T obita, S. Hirata, and R. J. Bartlett, J. Chem. Ph ys. 114, 9130 (2001). [9] C. S. W ang and J. Calla w a y , Comput. Ph ys. Comm un. 14, 327 (1978). [10] J. A. Applebaum and D. R. Hamann, Solid State Comm un. 27, 881-83 (1978). [11] J. W. Min tmire and J. R. Sabin, In t. J. Quan tum. Chem. Symp. 14, 707 (1980). [12] J. C. Bo ettger and S. B. T ric k ey , J. Ph ys. F 16, 693 (1986). [13] C. Pisani, R. Do v esi and C. Ro etti, Hartree-F o c k ab-initio treatmen t of crystalline systems, Lecture Notes in Chemistry 48 (Springer V erlag, Heidelb erg, 1988). [14] J. C. Bo ettger, In t. J. Quan tum Chem. Symp. 29 , 197 (1995). [15] K. Doll, Comput. Ph ys. Comm un, 137, 74, (2001). [16] O. H. Nielsen and R. M. Martin, Ph ys. Rev. B 32, 3780 (1985). [17] O. H. Nielsen and R. M. Martin, Ph ys. Rev. B 32, 3792 (1985). [18] P . J. F eib elman, Ph ys. Rev. B 44, 3916 (1991). [19] K. N. Kudin and G. E. Scuseria, Ph ys. Rev.B 61 , 5141 (2000). 103

PAGE 114

104 [20] J. O. Hirsc hfelder, C. F. Curtiss and R. B. Bird, Molecular Theory of Gases and Liquids (John Wiley and Sons, New Y ork, 1954). [21] P . Ew ald, Ann. Ph ys. 64, 253 (1921); see also S. W. DeLeeu w, J. W. P erram and E. R. Smith, Pro c. Ro y . So c. Lond. A373, 27 (1980). [22] P .-O. L owdin, Philosop. Magaz. Suppl. 5, 1, (1956). [23] J. Delhalle, J. M. Andr e, Ch. Demanet, and J. L. Br edas, Chem. Ph ys. Lett. 54, 186 (1978). [24] L. Piela and J. Delhalle, In t. J. Quan tum Chem. 13, 605 (1978). [25] J. Delhalle and L. Piela, Ph ys. Rev. B 22, 6254 (1980). [26] D. Jacquemin, J.-M. Andr e and B. Champagne, J. Chem. Ph ys. 5306, (1999). [27] L. E. McMurc hie and E. R. Da vidson, J. Comput. Ph ys. 26, 218, (1978). [28] J. M. P erez-Jord a and W. Y ang, J. Chem. Ph ys. 104, 8003 (1996). [29] L. Greengard and V. Rokhlin, J. Comput. Ph ys. 73, 325 (1987). [30] K. E. Sc hmidt and M. A. Lee, J. Stat. Ph ys. 63, 1223 (1991). [31] C. A. White, B. G. Johnson, P . M. W. Gill and M. Head-Gordon, Chem. Ph ys. Lett. 230, 8, (1994). [32] M. Challacom b e, C. A. White and M. Head-Gordon, J. Chem. Ph ys. 107, 10131, (1997). [33] E. L. P ollo c k and J. Glosli, Comput. Ph ys. Comm un. 95, 93, (1996). [34] U. Birk enheuer, J. C. Bo ettger, and N. R osc h, J. Chem. Ph ys. 100, 6826, (1994). [35] J. W. Min tmire, J. R. Sabin, and S. B. T ric k ey , 26, 1743 (1982). [36] J. C. Bo ettger and S. B. T ric k ey , Ph ys. Rev. B 32, 1356 (1985). [37] S. Hirata, M. T obita, M. T asumi, H. T orii, S. Iw ata, M. Head-Gordon, and R. J. Bartlett, P olymer v ersion 2.0 (2002). [38] A. D. Bec k e, J. Chem. Ph ys. 88, 2547 (1988). [39] B. G. Johnson, P . M. W. Gill, and J. A. P ople, J. Chem. Ph ys. 98, 5612 (1993). [40] S. Obara and A. Saik a, J. Chem. Ph ys. 84, 3963 (1986). [41] J. Alml of, K. F aegri, Jr., and K. Korsell, J. Comput. Chem. 2, 385 (1982).

PAGE 115

105 [42] M. H aser and R. Ahlric hs, J. Comput. Chem. 10, 104 (1989). [43] C. A. White and M. Head-Gordon, J. Chem. Ph ys. 101, 6593 (1994). [44] R. Kuttec h, E. Apr a and J. Nic hols, Chem. Ph ys. Lett. 238, 173 (1995). [45] E. Sc h w egler, M. Challacom b e and M. Head-Gordon, J. Chem. Ph ys. 109, 8764 (1998). [46] K. Kurki-Suonio, Israel J. Chem. 16, 115 (1977). [47] O. T reutler and R. Ahlric hs, J. Chem. Ph ys. 102, 346 (1995). [48] J. C. Slater, J. Chem. Ph ys. 41, 3199 (1964). [49] B. G. Johnson and M. J. F risc h, Chem. Ph ys. Lett. 216, 133 (1993). [50] R. S. P ease, Acta Cryst. 5, 356 (1962). [51] M. V ra ck o, C. M. Liegener and J. Ladik, In t. J. Quan tum Chem. 37, 241 (1990). [52] K. Koba y ashi, T. Sano and Y. J. I'Ha y a, Chem. Ph ys. Lett. 219, 53 (1994). [53] P . Y. Ay ala, K. N. Kudin and G. E. Scuseria, J. Chem. Ph ys. 115, 9698 (2001). [54] X. Blase, A. Rubio, S. G. Louie, and M. L. Cohen. Ph ys. Rev. B 51, 6868 (1995). [55] S. B. T ric k ey , R. J. Mathar, and J. C. Bo ettger in Pro c. Third UNAMCra y Sup ercomputing Conference: Computational Chemistry edited b y G. Cisneros, J. A. Cogordan, M. Castro, and C. W ang. (W orld Scien tic, Singap ore, 1997), 239. [56] A. Zunger, A. Katzir and A. Halp erin, Ph ys. Rev. B13, 5560 (1976). [57] K. Ragha v ac hari, G. W. T ruc ks, J. A. P ople and M Head-Gordon, Chem. Ph ys. Lett. 157, 479 (1989). [58] R. J. Bartlett, J. D. W atts, S. A. Kuc harski and J. Noga, Chem. Ph ys. Lett. 165, 513 (1990). [59] F. Co ester, Nucl. Ph ys. 1, 421 (1958). [60] J. Cizek, J. Chem. Ph ys. 45 , 4256 (1966). [61] R. J. Bartlett and G. D. Purvis, In t. J. Quatum Chem. 14, 561, (1978). [62] G. D. Purvis and R. J. Bartlett, J. Chem. Ph ys. 76, 7918 (1982). [63] J. Noga and R. J. Bartlett, J. Chem. Ph ys. 86, 7041 (1987).

PAGE 116

106 [64] S. A. Kuc harski and R. J. Bartlett, J. Chem. Ph ys. 97, 4281 (1992). [65] J. Noga, R. J. Bartlett and M. Urban, Chem. Ph ys. Lett. 134, 126 (1987). [66] H. J. Monkhorst, In t. J. Quan tum Chem. (Symp.) 11, 421 (1977). [67] J. F. Stan ton and R. J. Bartlett, J. Chem. Ph ys. 98, 7029 (1993). [68] D. C. Comeau and R. J. Bartlett, Chem. Ph ys. Lett. 207, 414 (1993). [69] L. Adamo wicz, W. D. Laidig and R. J. Bartlett, In t. J. Quatum Chem. 18, 245 (1984). [70] E. A. Salter, G. W. T ruc ks and R. J. Bartlett, J. Chem. Ph ys. 90, 1752 (1988). [71] E. A. Salter and R. J. Bartlett, J. Chem. Ph ys. 90, 1767 (1988). [72] I. Sha vitt and R. J. Bartlett, Man y-Bo dy Metho ds in Quan tum Chemistry , in preparation. [73] R. J. Bartlett in Mo dern Electronic Structure Theory , P art I edited b y D. R. Y ark on y . (W orld Scien tic, Singap ore, 1995) 1047. [74] R. J. Bartlett and G. D. Purvis I I I, In t. J. Quan tum Chem. Symp. 14, 561 (1978). [75] G. Fitzgerald, R. J. Harrison, and R. J. Bartlett, J. Chem. Ph ys. 85, 5143 (1986). [76] J. Gauss, J. F. Stan ton, and R. J. Bartlett, J. Chem. Ph ys. 95, 2623 (1991). [77] L. Adamo wics, W. D. Laidig, and R. J. Bartlett, In t. J. Quan tum Chem. Symp. 18, 245 (1984). [78] E. A. Salter, G. W. T ruc ks, and R. J. Bartlett, J. Chem. Ph ys. 90, 1752 (1989). [79] Marco H aser, Theor. Chim Acta, 87, 147, (1993). [80] P . Y. Ay ala and G. E. Scuseria, J. Chem. Ph ys. 110, 3660, (1999). [81] M. Sc h utz, G. Hetzer and H.-J. W erner, J. Chem. Ph ys. 111, 5691 (1999). [82] C. Hample and H.-J. W erner, J. Chem. Ph ys. 104, 6286 (1996). [83] M. Sc h utz, J. Chem. Ph ys. 113, 9986 (2000). [84] G. E. Scuseria and P . Y. Ay ala, J. Chem. Ph ys. 111, 8330 (1999). [85] S. Suhai, Ph ys. Rev. B 50, 14791 (1994).

PAGE 117

107 [86] S. Suhai, J. Chem. Ph ys. 101, 9766 (1994). [87] S. Suhai, Ph ys. Rev. B 51, 16553 (1995). [88] J.-Q. Sun and R. J. Bartlett, Ph ys. Rev. Lett. 77, 3669 (1996). [89] J.-Q. Sun and R. J. Bartlett, J. Chem. Ph ys. 108, 301 (1998). [90] J.-Q. Sun and R. J. Bartlett, J. Chem. Ph ys. 104, 8553 (1996). [91] S. Hirata and S. Iw ata, J. Chem. Ph ys. 109, 4147(1998). [92] W. F orner, R. Knab, J. C zek, and J. Ladik, J. Chem. Ph ys. 106, 10248 (1997). [93] P . Reinhardt, Theor. Chem. Acc. 104, 426 (2000). [94] M. Y u, S. Kalv o da, and M. Dolg, Chem. Ph ys. 224, 121 (1997). [95] A. Ab durahman, A. Sh ukla, and M. Dolg, J. Chem. Ph ys. 112, 4801 (2000). [96] A. Abrurahman, A. Sh ukla, and M. Dolg, Chem. Ph ys. 257, 301 (2000). [97] S. Hirata, I. Grab o wski, M. T obita, and R. J. Bartlett, Chem. Ph ys. Lett. 345, 475 (2001). [98] P . Pula y , Chem. Ph ys.Lett. 73, 393 (1980). [99] W. J. Hehre, R. F. Stew art, and J. A. P ople, J. Chem. Ph ys. 51, 2657 (1969). [100] W. J. Hehre, R. Ditc held, and J. A. P ople, J. Chem. Ph ys. 56, 2257 (1972). [101] P . C. Hariharan and J. A. P ople, Theor. Chim. Acta 28, 213 (1973). [102] A. P . Rendell and T. J. Lee, J. Chem. Ph ys. 94, 6219 (1991).

PAGE 118

BIOGRAPHICAL SKETCH Motoi T obita w as b orn in Kohzaki-mac hi, Chiba-pref., Japan, whic h is an eastern suburban part of T oky o. While he w as in high-sc ho ol, the rst Japanese astronaut w en t in to space. The news stim ulated him enough to pursue his career as an astronaut. Since he w as more in terested in scien tic researc h than w orking on space v ehicles, he decided to ma jor in materials science and engineering at W aseda Univ ersit y , T oky o. In his senior thesis researc h, he w as in v olv ed in a theoretical prediction of molecular materials using quan tum c hemistry . His target system w as h ydrogen b onding b et w een DNA base pairs. Through the course of the senior thesis w ork, he realized ho w m uc h computer based c hemistry p oten tially could help designing new materials and molecules and he w an ted to do more in quan tum c hemistry . After earning B. Eng. in materials engineering he had sev eral dieren t dreams. Three of them w ere (1) to study at the highest lev el institution in quan tum c hemistry , (2) to get a Ph.D degree, and (3), to master English. The second and third are requiremen ts to apply to b e an astronaut in his coun try . Coming to the Univ ersit y of Florida as a graduate studen t seemed to satisfy all three. A t the Univ ersit y of Florida, he w as a c hemistry graduate studen t under the sup ervision of Prof. R. Bartlett. Cho osing his ma jor as c hemical ph ysics, he had to sp end a lot of time studying c hemistry and ph ysics. Sub jects lik e quan tum mec hanics, inorganic c hemistry and analytical c hemistry w ere totally new to the B. Eng. bac kground studen t. Nev ertheless, surviving the rst y ear completing all of his course w ork successfully and b eing ev aluated v ery w ell b y his studen ts when he taugh t c hemistry lab with his clumsy English, ga v e him a lot of condence to 108

PAGE 119

109 nish the Ph.D program. In his researc h, Prof. Bartlett alw a ys considered his bac kground and in terests. One of his researc h topics w as high energy densit y materials as a p ossible ro c k et prop ellan t. Another, and one of the main topics of his dissertation, is t w o-dimensional materials suc h as surfaces or thin-lms. Through his study , he learned a lot ab out n umerical computer programming and its application. He wishes to use his scien tic kno wledge in space science and engineering after completing his Ph.D program, and also wishes to b e an astronaut as a scien tist if a c hance comes to him. A t this momen t, the c hance to apply as an astronaut has not o ccurred. In his leisure time, he w as m uc h in v olv ed in pla ying comp etitiv e sp orts, esp ecially tennis. He pla y ed as a mem b er of the Univ ersit y of Florida tennis club team for three y ears, and also pla y ed in the UST A (United States T ennis Asso ciation) league, whic h included a n um b er of times tra v eling to comp ete against other tennis teams. Also, he enjo y ed b eing a mem b er of the In ternational Gourmet Asso ciation and v olun teered as treasurer for t w o semesters. He met a great v ariet y of p eople through these activities.

PAGE 120

I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Ro dney J. Bartlett, Chair Graduate Researc h Professor of Chemistry I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Sam uel B. T ric k ey Professor of Ph ysics I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. N. Yngv e Ohrn Professor of Chemistry I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. Sam uel O. Colgate Professor of Chemistry I certify that I ha v e read this study and that in m y opinion it conforms to acceptable standards of sc holarly presen tation and is fully adequate, in scop e and qualit y , as a dissertation for the degree of Do ctor of Philosoph y. James D. Winefordner Graduate Researc h Professor of Chemistry

PAGE 121

This dissertation w as submitted to the Graduate F acult y of the College of Lib eral Arts and Sciences and to the Graduate Sc ho ol and w as accepted as partial fulllmen t of the requiremen ts for the degree of Do ctor of Philosoph y. Ma y 2002 Dean, College of Lib eral Arts and Sciences Dean, Graduate Sc ho ol

PAGE 122

DEVELOPMENT OF EFFICIENT ELECTR ON CORRELA TION METHODS F OR ONEAND TW O-DIMENSIONAL EXTENDED SYSTEMS AND THEIR APPLICA TIONS Motoi T obita (352) -392-6365 Departmen t of Chemistry Chair: Ro dney J. Bartlett Degree: Do ctor of Philosoph y Graduation Date: Ma y 2002 One of the imp ortan t questions for theoretical c hemists/ph ysicists/materials scien tists is that is it p ossible to in v en t new molecules or materials using the computer. The answ er is turning more and more to `Y es' as the ev olution of computer hardw are, soft w are and theories underling the soft w are. T rying to mak e new molecules or materials exp erimen tally tak es so man y times of trial and error. The step is v ery often time-consuming and exp ensiv e, and sometimes v ery dangerous suc h as explosiv e or p oisonous or can b e imp ossible to p erform. The computer can certainly do this preliminary trial and error pro cess m uc h faster, safer and c heap er. F or example, w e can no w tell p eople that \this semi-conductor materials starts conducting curren t at the 0.5 eV of external v oltage applied" or that \a c hemical reaction ABC ! AB + C releases so m uc h energy whic h is more p o w erful than dynamite" without p erforming an exp erimen t but using only the computer. W e only need to kno w what atoms and molecules materials consist of. Th us, our computational to ol are univ ersal to an y kind of molecules and materials, in principle. W e can do ev en more. Predicting crystal structure under sup er high pressure and/or temp erature ma y b e one example, switc hing one of the atoms in the DNA molecule with a dieren t atom ma y b e another. These kind of sim ulations on non-existing molecules/materials are the p o w er of computer aided molecule/materials design. W e could mak e suc h highly accurate predictions only for atoms and molecules consisting of up to only sev eral atoms un til tens of y ears ago, ho w ev er recen t rapid adv ances of computer and soft w are dev elopmen t, w e can

PAGE 123

no w conduct high accuracy calculation for large molecules or materials. What I mean b y \high accuracy" is that the underling metho dology for computation uses quan tum mec hanics rather than classical mec hanics. This enables us to understand what is going on at the atomic size scale. My dissertation topic is theory and soft w are dev elopmen t p erforming computer aided materials design, with some emphasis on large molecules, p olymers, and surfaces using the highly accurate quan tum mec hanical treatmen ts that include \electron-correlation". With the theory and soft w are I ha v e dev elop ed, it is no w p ossible to mak e most accurate prediction of prop erties of materials for relativ ely simple p olymers and surfaces.