ㄧ個新密度泛函理論推導公式之三維有限元素解

全文

(1)國立交通大學應用數學系碩士論文一個新密度泛函理論推導公式之三維有限元素數值解 3D Finite Element Solution of a New Formulation in Density Functional Theory. 研究生：蔡維哲指導老師：劉晉良教授. 中華民國九十四年六月.

(2) 一個新密度泛函理論推導公式之三維有限元素數值解. 3D Finite Element Solution of a New Formulation in Density Functional Theory. 研究生：蔡維哲. Student：Wei-Che Tsai. 指導教授：劉晉良. Advisor：Jinn-Liang Liu. 國立交通大學應用數學系碩士論文. A Thesis Submitted to Department of Applied Mathematics College of Science National Chiao Tung University in Partial Fulfillment of the Requirements for the Degree of Master in Applied Mathematics June 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年六月.

(3) 一個新密度泛函理論推導公式之三維有限元素數值解指導老師：劉晉良博士. 研究生：蔡維哲. 國立交通大學應用數學系. 摘. 要. 在這一篇論文中有二個部分。在第一個部分裡，我們主要探討電子限制在量子奈米結構中的情況。我們主要學習使用 self-consistence 迭代方式來處理一個非線性偏微分系統，這系統是由普瓦松方程式與薛丁格方程式共同組成。藉由這個非線性偏微分系統來描述電子在量子奈米結構中的變化。在第二個部分裡，我們嘗試使用許正餘教授在密度泛函理論所推導出來一個新的公式，從這新的公式出發，使用有限元素法希望能夠得到電子在原子與分子結構中很好的描述，我們在處理這問題同時，也嘗試學習使用更有效率的數值方法來處理有限元素法整理出來的大型矩陣。. i.

(4) 3D Finite Element Solution of a New Fromulation in Density Functional Theory. Student：Wei-Che Tsai. Advisor：Dr. Jinn-Lieng Liu. Department of Applied Mathematics National Chiao Tung University. Abstract There are two parts in the thesis. In the first part, we discuss that electrons confine in quantum nanostructures. We use self-consistence iteration method for solving the nonlinear PDE system. The system consists of Poisson equation and Schrodinger equation. It describes the electronic changes in quantum nanostructures by the nonlinear PDE system. In the second part, we try to compute a new formulation derived by Prof. Hsu. From the formulation, we hope to get good numerical results by using finite element method. The main goal of this thesis is to simulate electronic properties of well-known atomic models by using our still under development finite element codes for future research on density functional theory.. ii.

(5) 誌謝. 陳之藩先生說:「凡事得之於人者太多，出之於己者太少。」我在碩士生涯中深刻體會到這句話的深意。兩年碩士生涯中，我最要感謝我的指導教授劉晉良教授，劉教授引導我從徬徨到冷靜面對問題，指導我如何做學問，也教導我正確的研究方法與精神，雖然是短短兩年，卻讓我在求學過程中得到最大成長，也因為劉教授，讓我能夠接觸不同領域的學者，開拓我的視野，劉老師謝謝您。我還要特別感謝博士班陳人豪學長，學長用他豐富的研究經驗，帶領我在黑暗中前進，讓我能夠整理出精采的碩士論文，我要特別謝謝陳學長。其次我要感謝我的家人，由於家人的鼓勵與全力支持，我才能無後顧之憂地追求理想，真的非常感謝我的家人。最後，我要感謝我的朋友、同學以及所有曾經給我建議、鼓勵的師長們，因為有你們，我才能夠不斷反省、鞭策自己，因為有你們，讓我可以持續地成長、茁壯，謝謝你們。. iii.

(6) Contents Abstract (in Chinese). i. Abstract (in English). ii. Acknowledgement. iii. Contents. iv. I Electron confinement in quantum nanostructures. 1. 1 Introduction. 1. 2 The Model Problem and Numerical Methods. 1. 2.1 Finite Element Method ………………………………………… 2 2.2 Linear system solvers…………………………………………... 6 2.3 Numerical Results……………………………………………… 8. References. 10. Appendix. 11. II 3D Finite Element Method of Self-consistence Density Functional Theory. 12. 3 Introduction. 12. 4 Literature Survey. 12. 4.1 Theoretical Development………………………………………..12 4.1.1 Hartree approximation…………………………………… .12 4.1.2 Hartree-Fock (HF) approximation………………………...13 4.1.3 Thomas-Fermi (TF) approximation……………………….13 4.1.4 Hohenberg-Kohn theorm………………………………….13 4.1.5 Kohn-Sham equations……………………………………..13. iv.

(7) 4.1.6 Hsu Formulation…………………………………………..14 4.2 Numerical Methods……………………………………………..14. 5 Density Functional Theory. 16. 5.1 Hartree approximation…………………………………………..16 5.2 Hartree-Fock approximation…………………………………….16 5.3 The Kohn-Sham Equations……………………………...............17 5.4 A New Density Functional Theory……………………………...19. 6 Numerical Methods. 20. 6.1 Finite Element Method………………………………………….20 6.2 Numerical Linear Algebra………………………………………28. 7 Numerical Results and Conclusions. 31. 7.1 Hsu’s Equation…………………………………………………..31 7.2 Hartree’s results………………………………………………….32 7.3 Kohn Sham Equation……………………………………………32 7.4 Virial theorem check…………………………………………….33 7.5 3D Finite Element Codes………………………………………..33. References. 36. Appendix. 38. v.

(8) Part I. Electron con…nement in quantum nanostructures 1. Introduction. In the world of quantum, the interesting problems are the problems in threedimension (3D). But the complicated problems are always constructed from the simple problems. Thus, we compute the self-consistent electron states and con…ning potential, V (r; T ), for laterally con…ned cylindrical quantum wires at a temperature T from a numerical solution of the coupled Poisson and Schr oÄdinger (PS) equations. Our purposes are two fold. First, we compute the self-consistent electron density function n(r; T ) and con…ning potential V (r; T ) for laterally con…ned cylindrical quantum wires from a numerical solution of the coupled Poisson and Schr oÄdinger (PS) equations. The second purpose is to contrast the results of the relatively simple ThomasFermi approach with those of the coupled PS theory and thereby assess the validity of the Thomas-Fermi approximation in laterally con…ned nanostructures.. 2. The Model Problem and Numerical Methods. We begin by obtaining the eigenstates of an electron that is bound in two dimensions by a cylindrically symmetric nanostructure. The wave function factors in the usual way, 1 ªl;m;k = p ©l;m (r)eimµ eikz (2.1) L where m = 0; §1; §2; :::is the azimuthal quantum number, l = 1; 2; 3; :::is the radial quantum number, and L is the length of the cylinder. The subband spectrum is then labeled by three quantum numbers, El;m;k = (¹h2=2m¤ )k 2 + ¸l;m , where ¸l;m is the eigenvalue associated with the eigenstate ©l;m of the radial Schr oÄdinger equation, ". ". #. #. ¡¹ h2 1 d d m2 r ¡ ©l;m (r) + V (r)©l;m (r) = ¸l;m ©l;m (r) 2m¤ r dr dr r2. (2.2). and the boundary condition is : ©l;m (R) = 0, at the exterior surface. 1. (2.3).

(9) We assume that the ©l;m have been normalized according to Z. 2¼. R. 0. drr©2l;m (r) = 1.. (2.4). Using the one-electron wave function ªl;m;k , the density function is given by n(r) = N1. X l;m. ©2l;m (r)F¡ 1 [¯(¹ ¡ ¸l;m )] 2. (2.5). where and. N1 ´ (2m¤ =¼¯¹ h2 )1=2 1 F¡ 1 [¯(¹ ¡ ¸l;m )] = p 2 ¼. Z. 0. 1. 1. dxx¡ 2 [1 + exp(x ¡ ¯(¹ ¡ ¸l;m ))]¡1. ,where F¡ 21 is a Fermi-Dirac integral, ¯ is 1=(KB T ), KB T ´ ¹h2 (3¼2 Nd )2=3=2m¤ , and ¹ is the chemical potential(Fermi-level). Equation (2.5) gives the density of the electron gas at a distance r from the central axis. We then obtain the electrostatic potential Á from Poisson’s equation , ¡. 1 d dÁ(r) e [r ] = fNd ¡ n[Á(r)]g , r dr dr ". (2.6). where Nd is the number density of donors and where " is the dielectric constant. The potential energy V of the electron is given by V = ¡eÁ. The boundary condition of equation (2.6) is V (R) ¡ ¹ = 0:7eV. 2.1 2.1.1. Finite Element Method Introduction. The basic idea in the …nite element method is to …nd the solution in in…nite dimensional space by approximating it in a …nite dimensional subspace. Moreover, in the …nite element method, it will often be possible to improve or re…ne the approximate solution by spending more computational e¤ort. 2.1.2. Numerical procedure. In this approach, starting from an initial guess for Á, one repeats solving the Schr oÄdinger equation , computing n(r), and then updating Á(r) , until the iterates converge to a the self-consistent solution. Employing the standard 1D …nite element method, the radial Schr oÄdinger equation (2.2) is equivalent to the matrix eigenvalue problem AÁ = ¸BÁ, where A and B are tridiagonal. 2.

(10) Algorithm) (1)To get an initial guess © by solving Poisson’s eqnation. (2)According to V (r) = ¡e¤ ©(r),we can know the potential energy V of the electron. (3)To Get ©l;m (r) by sloving the SchrÄ odinger equation. RR (4)Eigenvectors normalized by 2¼ 0 drr ©2l;m (r) = 1. (5)Compute n(r). The density function is given by n(r) = N1. X l;m. ©2l;m (r)F¡ 1 [¯(¹ ¡ ¸l;m )] 2. (6)To obtain the electrostatic potential © from Poisson’s equation. (7)Repeating the procedures (1)~(6) until the iterates converge. 2.1.3. Discretization. Finite Element Method for 1D Poisson’s Equation 1 d d©(r) e [r ] ´ fNd ¡ n[©(r)]g r dr " dr Boundary Conditions: V (R) ¡ ¹ = 0:7eV ¡. .. Multiply (2:6) by an aribitrary test function v(x), and integrate over (0; R).Now choose v(x) satisfying v(R) = v(0) = 0, and then we get the form A(©; v) = .. R. ZR 0. Z d© dv V d© ¡ ); F (v) = fvdr ( r dr dr dr 0. =) A(©h ; vh ) = F (vh ), for all vh 2 fv 2 S h ; v(0) = v(R) = 0g.Without considering 8 B.C.s, we have the matrix formulation AU = F . ri¡1 ri 1 ri¡1 1 > > > a = ln(2h + r )( ¡ + ) + ln(2r )(¡ ) + ln(2h + 2r )( ) ii i¡1 i¡1 i > > 2 2 2 > h h h h h > > ri > > + ln(2ri )( 2 ); i = 1:::N < h A=> ri¡1 ri¡1 > ) + ln(2r )( ); i = 2:::N a = ln(2h + 2r )(¡ > i¡1 i¡1 i(i¡1) > > h2 h2 > > > 1 ri¡1 1 ri¡1 > > ) + ln(2ri¡1)( + 2 ); i = 2:::N : a(i¡1)i = ln(2h + 2ri¡1 )(¡ ¡ 2 h h h h 3 2 U1 7 6 6 U2 7 7 6 6 : 7 7 6 7 U =6 6 : 7 7 6 6 : 7 7 6 6 : 7 5 4 UN 3.

(11) F =. ri+1 R. (Fi =. ri¡1. f ©i ¼ f (ri )h; i = 2:::N ¡ 1). h h F1 = f (r1 ); FN = f(rN ) 2 2 Impose B.C., Solve AU = F by the Preconditional Conjugate Gradient Method.. 2.1.4. Finite Element Method for 1D SchrÄ odinger’s Equation m2 ¹h 1 d d [r ] ¡ 2 ]©l;m (r) + V (r)©l;m (r) = ¸l;m ©l;m (r) ¡ ¤[ 2m r dr dr r Boundary condition:©(R) = 0. The process of Finite Element Method:. (1). 1. (1)Choosing an arbitrary test function v(r) and integrate over (0; R). =) ¡. Z. ¹2 h ( 2m¤. R. 0. =). Z R Z R Z R 2 1 d d©(r) m e ©(r)v(r)) + V ©v = ¸ ©v (r )v(r) ¡ r dr dr r2 0 0 0. ¹ 2 v(r) d©(r) R Z R d©(r) v(r) h ¹ 2 2 Z R ©(r)v(r) h ¡ ¤( (r )j0 ¡ )d( )) + m (r 2m r dr dr r 2m¤ r2 0 0 +Ve. =¸. Z. Z. R. 0. R. 0. ©v. ©v. Set the test function v(0) = v(R) = 0 =) h ¹2 2m¤. Z. R. 0. =) h ¹2 2m¤ ). Z. 0. R. dv(r) Z R d©(r) r dr ¡ v(r) h ¹ 2 2 Z R ©(r)v(r) e Z R )( )+ m + V ©v = ¸ ©v (r dr r2 2m¤ r2 0 0 0 d©(r) dv(r) d©(r) v(r) h ¹ 2 2 Z R ©(r)v(r) e ( ¡ )+ ¤ m +V dr dr dr r 2m r2 0 A(. N X. Uj ©j ; vh ) = ¸B(. N X. Z. 0. R. ©v = ¸. Z. 0. R. ©v. Uj ©j ; vh ). j=1. j=1. =) Z Z ¹ 2 Z d© dv Z d© v m2h h ¹2 Z v ~ A(©; v) = ( ¡ )+ ©+V ©v , B(©; v) = ©v 2m¤ 2m¤ dr dr dr r r2. 4.

(12) aij = A(©j ; ©i ); bij = B((©j ; ©i ) We have to only compute ai(i¡1) ; aii ; ai(i+1) ; bi(i¡1) ; bii ; bi(i+1) ,other elements are zero. Example for ai(i+1) : h ¹ 2 Z R d©i+1 (r) d©i (r) d©i+1 (r) ©i (r) ( ¡ ) 2m¤ 0 dr dr dr r ©i+1 (r)©i (r) e Z R +V ©i+1 ©i r2 0. ai(i+1) = A(©i+1 ; ©i ) = +. h ¹2 2 Z R m 2m¤ 0. Linear Transformation:. »+1 r ¡ ri = ri+1 ¡ ri 2. (» + 1)h + ri 2 Linear shape function:Ã1 (») = =) r =. ai(i+1) =. ai(i+1). ai(i+1). 1¡» ; Ã2 (») 2. =. 1+» 2. Z 1 2 h ¹ 2 Z 1 dÃ2 (») dÃ1 (») 2 2 h h Ã1 (») dÃ 2(») ¢ ¢ d») ( ( ¢ ( ) ¢ d») ¡ ¤ 2m ¡1 d» d» h 2 d» (» + 1)h h 2 ¡1 + ri 2 Z 1 Z 1 2 h Ã2 (»)Ã 1 (») h ¹h ¢ d» + Ve Ã2 (»)Ã 1 (») ¢ d» + ¤ m2 (» + 1)h 2m 2 2 ¡1 ¡1 ( + ri )2 2 Z 1 Z 1 1 2 1 1¡» h ¹2 h ¹2 (¡ ¢ d») ¡ ( d») = ¤ ¤ 2m ¡1 4 h 2m ¡1 2 (» + 1)h + 2ri Z 1 Z 1 h 1¡» 1+» h h ¹2 1 ¡ »2 e ¢ d» + V ¢ ¢ d» + ¤ m2 2 2m 2 2 2 2 ¡1 ¡1 ((» + 1)h + 2ri ). h ¹2 1 ¹2 1 Z 1 h 1¡» h ¹2 2 h Z 1 1 ¡ »2 = (¡ ) ¡ d») + m ¢ ( 2m¤ h 2m¤ 2 ¡1 (» + 1)h + 2ri 2m¤ 2 ¡1 ((» + 1)h + 2ri )2 h 4 +Ve ¢ ¢ ( ) 8 3 5.

(13) Using Gaussian Quadrature(Zienkiewicz and Taylor, The FEM,#4th ed,McGrawHill,1989) R1 1 ¡ »2 1¡» d»), ¡1 d». (» + 1)h + 2ri ((» + 1)h + 2ri )2 Hence we can get all values of ai(i+1) ,i = 1; 2; :::; N ¡ 1. compute. R1. ¡1 (. By the same way, the values of ai(i¡1) ; aii ; bi(i¡1) ; bii ; bi(i+1) will be obtained easily. Discretization completed, and the next topic is how to …nd out the e¢cient solver of AU = ¸BU . For large scale matrix,we use Jacobi-Davidson method to compute the eigenpairs.. 2.2. Linear system solvers. 2.2.1. Jacobi-Davidson Method. Algorithm , Jacobi-Davidson method for linear eigenvalue problem:. 2.2.2. Algorithm , JD method for Linear Eigenvalue Problem Conjugate Gradient Method. The pseudocode for the Preconditioned Conjugate Gradient Method is given below. 1. Compute r(0) = b ¡ Ax(0) for some initial guess x(0) . 6.

(14) Begin Loops for i = 1; 2; ::: 2. Solve M z(i¡1) = r(i¡1) T. 3. ½i¡1 = r(i¡1) z (i¡1) 4. If i = 1 p(1) = z (0) Else ¯ i¡1 = ½i¡1=½i¡2 p(i) = z (i¡1) + ¯ i¡1p(i¡1) End If 5. q(i) = Ap(i) T. 6. ®i = ½i¡1 =p(i) q(i) 7. x(i) = x(i¡1) + ®i p(i) 8. r(i) = r (i¡1) ¡ ®i q(i) 9. Check convergence continue if necessary End Loops.. 7.

(15) 2.3. Numerical Results. 8.

(16) 9.

(17) References [1] R.L.Burden and J.D.Faires, Numerical Analysis, Brooks/Cole Publishing Company, 1997. [2] Luscombe and Bouchard and Luban, Electron con…nement in quantum nanostructures:Self-consistent Poisson-Schrodinger theory, Physical Review B,#volume 46,#number 16. [3] D.R.Fokkema, G.L.G.Sleijpen, and H.A.van der Vorst. Jocobi-Davidson style QR and QZ algorithms for the reduction of matrix pencils, SIAM J. Sci. Comput., 20(1):94-125, 1998. [4] J.L.Liu,Lecture Note of Numerical PDE, 2004. [5] O.C.Zienkiewicz and R.L.Taylor, The Finite Element Method: The Basis, Butterworth-Heinemann, 2000.. 10.

(18) APPENDIX. Fermi-Dirac Function. 11.

(19) Part II. 3D Finite Element Method of Self-consistent Density Functional Theory 3. Introduction. Density functional theory is an extremely successful approach for the description of metals, semiconductors, and insulators. The last decade has witnessed a proliferation in methodologies for numerically solving large-scale self-consistent eigenvalue problems. The result of Kohn-Sham equation in density functional theory will be reviewed. We mainly foucs on a new formulation in density functional theory by HSU [HSU,2003] using the 3D …nite element method. The new formulation that drops the exchange-correlation term makes the computation potentially much traceable in physics without any ad hoc assumption. Single atoms (e.g., H,He,Li,Be,B,...) will be included in our test programs. Besides, we will discuss some computations of excited state of single atomic system with considering Singlet and Triplet state. Summary of the above description, some objectives of the thesis are 1.To Implement Kohn-Sham equation in order to review the results of single atomic system. 2.To run programs for new Hsu’s DFT formulation. 3.To show the results of new Hsu’s DFT formulation imposing electrons of Singlet and Triplet state.. 4. Literature Survey. 4.1. Theoretical Development. The key problem in the structure of matter is to solve the SchrÄ odinger equation for a system of N interacting electrons in the external coulombic …eld created by a collection of atomic nuclei. It is a very di¢cult problem in many-body theory and, in fact, the exact solution is known only in the case of the uniform electron gas, for atoms with a small number of electrons and for a few small molecules. It usually gets these solutions by numerical skill. At the analytic level, one always has to resort to approximations. 4.1.1. Hartree approximation. The …rst approximation may be considered the one proposed by Hartree. It consists of postulating that the many-electron wave function can be written as a 12.

(20) simple product of one-electron wave functions. The Hartree approximation treats the electrons as distinguishable particles. 4.1.2. Hartree-Fock (HF) or self-consistent …eld (SCF) approximation. It is to introduce Pauli exclusion principle (Fermi statistics for electrons) by proposing an antisymmetrized many-electron wave function in the form of a Slater determinant. It has been for a long time the way of choice of chemists for calculating the electronic structure of molecules. In fact, it provides a very reasonable picture for atomic systems and also provides a reasonably good description of inter-atomic bonding. 4.1.3. Thomas-Fermi (TF) approximation. Parallel to the development of this line in electronic structure theory, Thomas and Fermi proposed, at about the same time as Hartree (1927-1928), that the full electronic density was the fundamental variable of the many-body problem and derived a di¤erential equation for the density without resorting to one-electron orbitals. The Thomas-Fermi (TF) approximation was actually too crude because it did not include exchange and correlation e¤ects and was also unable to sustain bound states because of the approximation used for the kinetic energy of the electrons. However, it set up the basis for the later development of density functional theory (DFT), which has been the way of choice in electronic structure calculations in condensed matter physics and it also became accepted by the quantum chemistry community because of its computational advantages compared to HF-based methods. 4.1.4. Hohenberg-Kohn theorem. In 1964, Hohenberg and Kohn formulated and proved a theorem that put on solid mathematical grounds the former ideas, which were …rst proposed by Thomas and Fermi. In Hohenberg-Kohn theorem the electronic density determines the external potential, but it is also needed that the density corresponds to some ground state antisymmetric wave function, and this is not always the case. Up till now, both the exact ground state density as well as the Hohenberg-Kohn functional is still unknown, so one cannot make use of the Hohenberg-Kohn theorems to calculate the molecular properties. 4.1.5. Kohn-Sham equations. In 1965, Kohn and Sham proposed the idea of replacing the kinetic energy of the interacting electrons with that of an equivalent non-interacting system, because this latter can be easily calculated. Kohn and Sham introduced a …ctitious system of N non-interacting electrons to be described by a single determinant wave function in N ”orbits.” We have seen that a reasonably good description can be obtained by separating the electrostatic (classical Coulomb energy), exchange and correlation contributions. The biggest di¢culty is to deal with correlation.This 13.

(21) is, in fact, an active …eld of research that has produced signi…cant improvements in the past decade. This correlation term has to tread on an approximate manner. Although there are many di¤erent functions available, almost all of them are derived from the electron density of a uniform electron gas, which can be calculated by means of statistical thermodynamics. 4.1.6. Hsu Formulation. A generic derivation from cluster expansion results in a new DFT formulation without the exchange-correlation term that makes the computation much traceable in physics without ad hoc (e.g., local density spin approximation (LDSA)) assumption.[J.Y.Hsu, PRL 2003]. 4.2. Numerical Methods. The last decade has witnessed a proliferation in methodologies for numerically solving large-scale problems in physics. The rapid growth has been driven because a wide variety of computational methods and numerical algorithms have exploited that physical locality. This review examines one subset of these new computational methods, namely real-space techniques. Real space methods can be loosely categorized as one of three types: …nite di¤erences (FD), …nite elements (FE), or wavelets. All three lead to structured, very sparse matrix representations of the underlying di¤erential equations on meshes in real space. Applications of wavelets in electronic structure calculations have been thoroughly reviewed recently [Arias, 1999]. As implied in the title, the primary focus is on calculations in density functional theory (DFT); real-space methods are in no way limited to DFT, but since DFT calculations comprise a dominant theme in modern electrostatics and electronic structure, the discussion here will mainly be restricted to this particular theoretical approach. The early development of FD and FE methods for solving partial di¤erential equations stemmed from engineering problems involving complex geometries, where analytical approaches were not possible [Strang and Fix, 1973]. Example applications include structural mechanics and ‡uid dynamics in complicated geometries. However, even in the early days of quantum mechanics, attention was paid to FD numerical solutions of the Schrodinger equation [Kimball and Shortley, 1934; Pauling and Wilson, 1935]. Real-space calculations are performed on meshes; these meshes can be as simple as Cartesian grids or can be constructed to conform to the more demanding geometries arising in many applications. Finite-di¤erence representations are most commonly constructed on regular Cartesian grids. They result from a Taylor series expansion of the desired function about the grid points. The advantages of FD methods lie in the simplicity of the representation and resulting ease of implementation in e¢cient solvers. Disadvantages are that the theory is not variational, and it is di¢cult to construct meshes ‡exible enough to conform to the physical geometry of many problems. Finite-element methods, on the other hand, have the advantages of 14.

(22) signi…cantly greater ‡exibility in the construction of the mesh and an underlying variational-type formulation. Other advantages include easier parallel implementation using domain decomposition and possible mesh re…nement in regions where solution changes rapidly, as mentioned earlier. However, the cost of the ‡exibility is an increase in complexity and more di¢culty in the implementation of multiscale or related solution methods [Thomas, 2000].Nevertheless, we mainly discuss the fundamentals of FEM solutions of the self-consistent nonlinear coupled Poisson and Schr oÄdinger (PS) equations in density-functional theory (DFT) in this thesis.. 15.

(23) 5. Density Functional Theory. The Hamiltonian in its original form is very complex: X X X Zk e2 e2 h2 X 2 h2 X 2 X e2 5Rk + 5ri + ¡ ¡ 2m i jri ¡ Rk j 2M k jRk ¡ Rl j i;j jri ¡ rj j i k;l k (2) which involves sums over all electrons/nuclei and their pairs involving kinetic energy and columbic potentials.. H=¡. 5.1. Hartree Approximation. Hartree approximation consists of postulating that the many-electron wave function can be written as a simple product of one-electron wave functions. It seems plausible that it might be useful to start with a wavefunction of the general form ©(r1; r2 ; r3 ; :::; rN ) = '1 (r1) ¤ '2 (r2 ) ¤ ::: ¤ 'N (rN ) which is known as a Hartree P roduct. While this functional form is fairly convenient, it has at least one major shortcoming: it fails to satisfy the antisymmetry principle, which states that a wavefunction describing fermions should be antisymmetric with respect to the interchange of any set of space-spin coordinates.. 5.2. Hartree-Fock Approximation. Hartree-Fock theory is fundamental to much of electronic structure theory. It is the basis of molecular orbital theory, which positulates that each electron’s motion can be described by a single-particle function (orbital) which does not depend explicitly on the instantaneous motions of the other electrons. Hartree-Fock theory was developed to solve the electronic SchrÄ odinger equation that results from the time-independent SchrÄ odinger equation after invoking the Born-Oppenheimer approximation (Neglects motion of nuclei [heavier than electrons]). In atomic units, and with r denoting electronic and R denoting nuclear degrees of freedom, the electronic SchrÄ odinger equation is 2. 3. X ZA ZB X 1 1 X 2 X ZA 4¡ 5 ª(r; R) = Eª(r; R) + + 5i ¡ 2 i r R r AB i>j ij A>B A;i Ai. (5.1). The Hartree approximation treats the electrons as distinguishable particles. A step forward is to introduce Pauli exclusion principle (Fermi statistics for electrons) by proposing an antisymmetrized many-electron wave function in the form of a Slater determinant:. 16.

(24) ¯ ¯ ¯ ' (r1 ; ¾ 1 ) ' (r2 ; ¾2 ) : : : ' (rN ; ¾N ) ¯ 1 1 ¯ ¯ 1 ¯ ¯ ¯ '2 (r1 ; ¾ 1 ) '2 (r2 ; ¾2 ) : : : '2 (rN ; ¾N ) ¯ ¯ ¯ ¯ 1 ¯¯ : : : : ¯ ©(r; R) = SDf'j (ri ; ¾ i )g = p ¯ ¯ : : : : ¯ N! ¯¯ ¯ ¯ ¯ : : : : ¯ ¯ ¯ ' (r1 ; ¾ 1 ) ' (r2 ; ¾2 ) : : : ' (rN ; ¾N ) ¯ N. N. N. Hartree-Fock approximation has been for a long time the way of choice of chemists for calculating the electronic structure. In fact, it provides a very reasonable picture for atomic systems and, although many-body correlations are completely absent, it also provides a reasonably good description of inter-atomic bonding.. 5.3. The Kohn-Sham Equations. The Kohn-Sham self-consistent eigenvalue equations for electronic structure can be written as follows: 1 (¡ 52 +Vef f )ªi (r) = "i ªi (r); (3) 2 where the density-dependent e¤ective potential is Vef f (r) = Ve¡nuc (r) + Ve¡e ([½(r)]) + Vxc ([½(r)]; r), R. 0. (4). 0. ½(r ) ±Exc Ve¡nuc = external potential, Ve¡e = jr¡r , Exc =exchange0 dr , Vxc = ±½ j correlation enery. Local Spin Density Approximation (LSDA) : exchange-correlation energy Exc is a simple known funtion Z 3 6 4=3 4=3 Exc [½(r)] = ¡ ( )1=3 (½" (r) + ½# (r))dr: (5) 4 ¼ The classical electrostatic potential Ve¡nuc (r) + Ve¡e (r) is due to both the electrons and nuclei, and the (in principle) exact exchange-correlation potential Vxc ([½(r)]; r) incorporates all nonclassical e¤ects. The exchange-correlation potential includes a kinetic contribution, since the expectation value of the KohnSham kinetic energy is that for a set of noninteracting electrons moving in the one electron e¤ective potential. The electron density ½(r) is obtained from the occupied orbitals :. ½(r) =. occup: X i. jªi (r)j2 :. (6). Standard electrostatics tells us that doing Ve¡e integral is equivalent to solving Poisson’s equation. Hartree Potential Ve¡e is solution of the Poisson equation: 52 Ve¡e = ¡4¼½(r):. We slove it by Conjugate Gradient Method once a distribution ½ is known. 17. (7).

(25) 5.3.1 1.. Self Consistency 1 (¡ 52 +Ve¡nuc (r)+Ve¡e ([½(r)])+Vxc ([½(r)]; r))ªi (r) = "i ªi (r); i = 1; :::; ioccup: : 2. 2. ½(r) =. occup: X i. jªi (r)j2 :. 3. 52 Ve¡e = ¡4¼½(r): The Self Consistent System can be viewed as a nonlinear eigenvalue problem, because Vxc and Ve¡e both depend on ½. Self-Consistent Iteration : 1. Initial Guess for V ,set Ve¡e = 0 2. 1 Slove (¡ 52 +Veff )ªi (r) = "i ªi (r) 2 3. Calculate new ½(r) =. occup: X i. jªi (r)j2. 4. Show 52 Ve¡e = ¡4¼½(r) for a new Ve¡e 5. Compute Vxc = f[½(r)] , Vnew = Ve¡nuc + Ve¡e + Vxc 6. If jVnew ¡ V j < T OL, STOP. 18.

(26) For example, we want to solve Li atom of con…guration 1s2 2s1 with 2 spin up and 1 spin down orbitals. First, we group all the orbitals into two types, spin up and down and then sum all the spin up orbitals to get ½" and ½# from spin down orbitals. Then put them into the exchange potential for di¤erent spin type. Finally, we solve the KS equations seperately with respect to di¤erent types of spin.We have ½" = jª(1s; ")j2 + jª(2s; ")j2 and ½# = jª(1s; #)j2 . Therefore, the equation we solve is 1 (¡ 52 +Vef f;" )ªi;" = "i;"ªi;" 2. (8). for spin up orbitals and replace " by # for spin down orbitals.We can have the . The correlation functional is the same as exchange form of Vxc;" = ±E[½";½#] ±½" functional if we can …nd or have its functional form. If we can not have the correlation functional form for di¤erent spin types, just use the functional form without considering the spin type.. 5.4. A New Density Functional Theory (J.Y.Hsu, PRL 2003). The paper (J.Y.Hsu, PRL 2003) gives a new derivation of Density Functional Theory from a cluster expansion by truncating the higher-order correlations in one and only one term in the kinetic energy. The new formulation admits excited states and allows self-consistent calculation of the exchange correlation e¤ect without any ad hoc assumptions. The wave function ª is chosen as the product of a single-electron wave function ©, and an N-body correlation function, (9). ª(¡) = ¦N i=1 ©(ri )U(¡). ¡ is the N-particle phase space point equivalent to the expression (r1; r2 ; :::; rn ). The exchange symmetry is imposed on U and on the indistinguishable particles so that each electron is described by the same ©. This gives the density function as follows: n(r) =. Z. 2 2 2 2 ¦N i=2 d¿ i j©(ri )j jU(r1 ; r2 ; :::; rN )©(r)j ´ M(r)j©(r)j ´ jª0 (r)j ; (10). where the index starting from i = 2 is chosen for convenience by imposing the exchange symmetry, and d¿ i is the volume element of ith particle. For the immobile ion I with charge ZI and N representing the number of electrons, the last derived equation is *. +. X 1 ZI (N ¡ 1) 1 "ª0 (r) = ¡ 52 ª0 (r) ¡ ; (11) ª0 (r) + ª0 (r) 2 2 jr ¡ r 0 j I jr ¡ RI j. 19.

(27) where. *. 1 jr ¡ r0 j. +. R. =R. d¿. 0. 0. jª0 (r )j2 0 jr¡r j. d¿ 0 jª0 (r0 )j2. (12). :. The derivation of the density functional theory (DFT) from the cluster expansion corrects the spurious self-interaction energy in the classical DFT, admits the excited states, and has a self-consistent exchange correlation e¤ect. Excited States (Singlet and Triplet States). Triplet States Singlet States For a particle j in an arbitrary orbital, the generalized equation is X X Z 1 ZI 1 jªj 0 (r )j2 "ªj (r) = ¡ 52 ªj (r) ¡ ªj (r) + ªj (r) ªj (r): d¿ 0 2 2 jr ¡ r0 j I jr ¡ RI j j 0 6= j (13) Varying a nonnegation trial function º(r1; r2 ) ´ exp(¯jr1 ¡ r2 j) and holding ª0 constant to minimize the total energy, subject to particale conservation 0. N=. Z. 2. d¿ 1 d¿ 2 jº(r1 ; r2)j ½(r1 )½(r2) ´. Z. d¿ 1d¿ 2jª(r1 ; r2)j2 :. (14). We reduce to the simpli…ed ¯ = 0 version X ZI 1 ! ! ¡ ! "ªs (¡ r1 ) = ¡ r2 ªs (¡ r1 ) ¡ (15) ¡ ! ªs ( r1 ) ¡ ! 2 j r1 ¡ RI j ! ! ! 1 XZ jª¾ (¡ r2 )j2 1 XZ ª¤¾ (¡ r2 )ªs (¡ r2 ) ¡ ! ! + d¿ 2 ¡ ª ( r ) ¨ d¿ ª¾ (¡ r1 ): 2 s 1 ! ¡ ! ¡ ! ¡ ! 2 jr ¡ r j 2 jr ¡ r j ¾6=S. 1. 2. ¾6=S. 1. 2. Example for excited Be(2 1S,1 2S,1 2P) X ZI 1 ª2S (r1 ) "ª2S (r1) = ¡ r2 ª2S (r1 ) ¡ 2 jr1 ¡ R1 j Z Z jª2P (r2 )j2 jª1S (r2 )j2 1 d¿ 2 + d¿ 2 ª2S (r1) + ª2s (r1 ) jr1 ¡ r2 j 2 jr1 ¡ r2 j Z ª¤ (r2)ª2S (r2 ) 1 d¿ 2 2P ª2P (r1 ): ¨ 2 jr1 ¡ r2 j. 20. (16).

(28) 6. NUMERICAL METHODS. 6.1 6.1.1. Finite Element Method Introduction. The Finite Element Method is a numerical method which can be used for the accurate solution of complex engineering problems. The method was …rst developed in 1956 for the analysis of aircraft structural problems. Thereafter, within a decade, the potentialities of the method for the solution of di¤erent types of applied science and engineering problems were recognized. Over the years, the …nite element technique has been so well established that today it is considered to be one of the best methods for solving a wide variety of practical problems e¢ciently. In fact the method has become one of the active research areas for applied mathematicians. One of the main reasons for the popularity of the method in di¤erent …elds of engineering is that once a general computer program is written, it can be used for the solution of any problem simply by changing the input data. Many problems that …nd out its approximate numerical solution to predict the response of physical system subjected to the external in‡uences arise in many areas of engineering, science, and applied mathematics. The Finite Element Method which is a computer-aid mathematical technique is just a powerful method for obtaining approximate numerical solution. Up to now,applications to date have occurred principally in the areas of solid mechanics, heat transfer, ‡uid mechanics, and electromagnetism. New areas of application are continually being discovered, recent ones include solid-state physics and quantum mechanics. 6.1.2. Procedure of Finite Element Method. Discretization The discretization of the domain into subregions is the …rst step in FEM. We choose the tetrahedral element as our basic element shape to partition the domain.. 21.

(29) tetrahedral element By the software (HyperMesh), we can partition the domain.. The meshes of 3D computation domain( r=20 Bohr). 22.

(30) Shape Functions The typical 3-D linear tetrahedral element trial solution can be written e U(r; ®) =. 4 X. ®j Nj (x; y; z). j=1. The special coordinates are introduced de…ned by x = L1 x1 + L2 x2 + L3 x3 + L4 x4. (8.a). y = L1 y1 + L2 y2 + L3 y3 + L4 y4. (8.b). z = L1z1 + L2 z2 + L3 z3 + L4 z4. (8.c). 1 = L1 + L2 + L3 + L4. (8.d). The coordinates of point P [Zienkiewicz and Taylor,2000]. 23.

(31) Sloving Eq(8) gives Li =. ai + bi x + ci y + di z etc. 6Vc ¯ ¯ ¯ 1 x1 y1 z1 ¯ ¯ ¯ ¯1 x y z ¯ ¯ 2 2 2¯ ¯ Vc = ¯¯ ¯ ¯ 1 x3 y3 z3 ¯ ¯1 x y z ¯ 4. and. 4 4. ¯ ¯ ¯ xj yj zj ¯ ¯ ¯ ai = ¯¯¯ xk yk zk ¯¯¯ ¯ x l yl z l ¯ ¯ ¯ ¯ 1 yj z j ¯ ¯ ¯ ¯ bi = ¡ ¯¯¯ 1 yk zk ¯¯ ¯ 1 yl zl ¯ ¯ ¯ ¯x 1 z ¯ j¯ ¯ j ci = ¡ ¯¯¯ xk 1 zk ¯¯ ¯ ¯ xl 1 z l ¯ ¯ ¯ ¯ x j yj 1 ¯ ¯ ¯ ¯ di = ¯¯¯ xk yk 1 ¯¯ ¯ x l yl 1 ¯. The linear shape functions for the linear element are simply Ni = Li ; Nj = Lj ; Nk = Lk ; Nl = Ll where Vc represents the volume of the tetrahedron. 6.1.3. Weighted residual approach. Variational approach requires a knowledge of a variational problem( functional to be extremized or made stationary ) for the given problem. Usually we encounter problems for which the variational principles are not known. Weighted residual approach in the …nite element method can be applied even in these cases.We can derive directly from the governing di¤erential equations of the problem without any need of knowing the ’functional’ by weighted residual approach. For a particle J of one nucleus system in an orbital, the governing equation is written as 1 Z ¡ 52 ªJ (r) ¡ ªJ (r) + ¦J ªJ (r) = "ªJ (r) (17) 2 r 0 2 Z 1 X 0 jªJ 0 (r )j d¿ (18) ¦J = 2 J 0 6= J jr ¡ r0 j. where ªJ is the density function of a a particle J in an orbital that is limited to two electrons with spin polarization to satisfy the Pauli exclusion principle, r and 24.

(32) R is the coordinate from the zero point of the computing space, Z is the number of positive charge of the nucleus, the subscript I means the kind of nucleus . STEP 1. e Using a trial solution U(r; ®) to approximate the density function ª(r), we …rst write down the residual equation for each orbital J : 1 Z e e e RJ (r; ®) = ¡ 52 U(r; ®) ¡ "U(r; ®) (19) ®) ¡ Ue (r; ®) + ¦J U(r; 2 r The typical 3-D element trial solution is always written in the general form P e U(r; ®) = nj=1 ®j Nj (r) , and one residual equation is Z Z Z. RJ (r; ®)Ni (r)dV = 0 ; i = 1; 2; :::; n. (20). where it integrates over one element and n is the number of nodes in an element. =) ¾ Z Z Z ½ Z 1 e e ®) Ni (r)dV = 0 ; i = 1; 2; :::; n ®) ¡ Ue (r; ®) + ¦J Ue (r; ®) ¡ "U(r; ¡ 52 U(r; 2 r (21) STEP 2. Divergence Theorem =) e e e 1 Z Z Z @Ni (r) @ U(r; ®) @Ni (r) @ U(r; ®) @Ni (r) @ U(r; ®) + + dV 2 @x @x @y @y @z @z Z Z Z e Ni (r)U(r; ®) ¡Z dV r Z Z Z + ¦J Ni (r)Ue (r; ®)dV. ¡". 1 = 2. Z Z Z. e Ni (r)U(r; ®)dV. I " e @ U (r; ®). @x. (22). e @ Ue (r; ®) @ U(r; ®) Ni (r)nx dydz + Ni (r)ny dxdz + Ni (r)nz dxdy @y @z. #. where nx , ny and nz are the direction vectors of the outward unit normal to the element boundary. The R.H.S. of above Equation is the boundary term contains the ‡ux must vanish from the system equations for the eigenproblems , so Equation (22) is rewritten as follows: e 1 Z Z Z @Ni (r) @ U(r; ®) @Ni (r) @ Ue (r; ®) @Ni (r) @ Ue (r; ®) + + dV (23) 2 @x @x @y @y @z @z Z Z Z e Ni (r)U(r; ®) dV ¡Z r Z Z Z e + ¦J Ni (r)U(r; ®)dV. ¡". =0. Z Z Z. e Ni (r)U(r; ®)dV. 25.

(33) STEP 3. Equation (23) is rearranged to satisfy the matrix solver as follows : f" [M] + [N]g f®g = 0 where [M ]ij =. Z Z Z. Ni (r)Nj (r)dV. (24) (25). Z Z Z. @Ni (r) @Nj (r) @Ni (r) @Nj (r) @Ni (r) @Nj (r) + + dV (26) @x @x @y @y @z @z Z Z Z Z Z Z Ni (r)Nj (r) +Z dV ¡ ¦J Ni (r)Nj (r)dV r. [N]ij = ¡. 1 2. STEP 4. To transfer the integrals of Equation(24) into a form appropriate for numerical evaluation. [M]ij =. Z Z Z. Ni (r)Nj (r)dV. (27). Vc 10 Vc ( if i 6= j ) = 20 ( if i = j ) =. (28). [N] = [A] + [B] + [C]. Z Z Z @Ni (r) @Nj (r) @Ni (r) @Nj (r) @Ni (r) @Nj (r) 1 + + dV (29) 2 @x @x @y @y @z @z bi bj + ci cj + di dj =¡ 72Vc. [A] = ¡. Z Z Z. Ni (r)Nj (r) dV (Gauss Cubic integration) r N (ai + bi xk + ci yk + di zk )(aj + bj xk + cj yk + dj zk ) Vc X q = wk 36 k=1 x2 + y2 + z 2. [B] = Z. k. k. (30). k. where N is the total numbers of the Gauss points, wk is the weighting factor and the subscript k means that the related values on the Gauss point .. 26.

(34) Numerical integration formula ofZGauss Z Z for tetrahedral element [Zienkiewicz and Taylor, 20 [C] = ¡ ¦j Ni (r)Nj (r)dV (31) STEP 5. We obtain the term ¦J from solving Poisson’s equation 52¦J (r) = ¡4¼jªJ (r)j2. (32). 1 the sphere radius of the domain. (33). Boundary Condition : ¦J (@) =. Do all process of …nite element method again as the same above. Poisson’s equation is rearranged to satisfy the matrix solver as follows : h i. h i. A [X] = B. 27.

(35) =) h i. Z Z Z. @Ni (r) @Nj (r) @Ni (r) @Nj (r) @Ni (r) @Nj (r) + + dV @x @x @y @y @z @z bi bj + ci cj + di dj =¡ 36Vc. A =¡. h i. B = ¡4¼. Z Z Z. jªJ (r)j2Ni (r)dV (Gauss Cubic integration). (34). (35). Deriving the element h iequations h icomplete by above steps. The Conjugate Gradient Method can slove A [X] = B , and it means the ¦j term is known now. We get RRR [C]ij = ¡ ¦j Ni (r)Nj (r)dV to utilize Gauss Cubic integration. All elemental matrix equations are assembled to be algebra system that are then solved by Jacobi-Davidson matrix solver.. 6.2 6.2.1. Numerical Linear Algebra Jacobi-Davidson Method. The Jacobi-Davidson solver can e¢ciently deal with the large-scale sparse eigenvalue matrix equation, which is still a very challenging task even nowadays. One of the advantages in using Jacobi-Davidson algorithm to solve the DFT eigenvalue problem is the feasibility of parallelization in the future, considering the computational demanding of the problem itself. The Jacobi-Davidson method [Voss and Betcke, 2002; Hwang, 2003] is based on a combination of two basic principles. The …rst one is to apply a Galerkin approach for the eigenproblem Ax = ¸x, with respect to some given subspace spanned by an orthonormal basis fv1; :::; vm g. The Galerkin condition is AVm s ¡ µVm s ? fv1 ; :::; vm g where Vm denotes the matrix with columns v1 to vm . This equation has m (m) (m) solutions (µj ; sj ). This m pairs are called the Ritz values and Ritz vectors of A with respect to the subspace spanned by the columns of Vm . These Ritz pairs are approximations for eigenpairs of A, and our goal is to obtain better approximations by a well-chosen expansion of the subspace. Suppose that we have an eigenvector approximation uj for an eigenvector x corresponding to a given eigenvalue ¸. We suggest computing the orthogonal correction t for uj(m). A(u(m) + t) = ¸(u(m) + t) j j (m). with t ? (uj ), t 2 U ? ´ fvjv¤ uj = 0g. The correction equation of Jacobi-Davidson is (A? ¡ µj I)t = ¡rj 28. (36).

(36) where rj = (A ¡ µj I)uj , A? = (I ¡ uj u¤j )A(I ¡ uj u¤j ). The next is to solve the t from Eq:(34) and add t into subspace to expand the search subspace, then iterate again with expansive subspace until convergence. The Jacobi-Davidson method has similar convergence properties as inverse iteration if the correction equation is solved exactly. The details of solving linear eigenvalue problem using Jacobi-Davidson method algorithm are summarized as follows [Wang et al., 2003]: 1. Given A(¸) = ¸A1 + A0 . 2. To choose a random vector V i as the initial subspace. 3. To compute the Galerkin condition as Mi = V ¤ Ai V . 4. To compute the Ritz pairs (µ; s) of (µM1 + M0 )s = 0 and select the desired Ritz pair to be eigenpair with ksk2 = 1. 5. To compute u = V s, and the residual r = A(µ)u. 6. If krk2 < ", ¸ = µ, x = u, Quit. 7. To compute correction term t and orthogonalize t against V , v = t= ktk2. 8. Expand V = [V; v]. 9. To back to process 3 (Restart : use a few of the last Ritz vectors as initial vectors ) and iterate until krk2 < ". 6.2.2. Conjugate Gradient Method. The Conjugate Gradient method is an e¤ective method for symmetric positive de…nite systems. It is the oldest and one of the best known nonstationary methods. (Nonstationary methods di¤er from stationary methods in that the computations involve information that changes at each iteration. Typically , constants are computed by taking inner products of residuals or other vectors arising from the iterative method.) The method proceeds by generating vector sequences of iterates i.e., successive approximations to the solution residuals corresponding to the iterates and search directions used in updating the iterates and residuals. Although the length of these sequences can become large only a small number of vectors needs to be kept in memory. In every iteration of the method two inner products are performed in order to compute update scalars that are de…ned to make the sequences satisfy certain orthogonality conditions. On a symmetric positive de…nite linear system these conditions imply that the distance to the true solution is minimized in some norm.. 29.

(37) The iterates x(i) are updated in each iteration by a multiple ®(i) of the search direction vector p(i): x(i) = x(i¡1) + ®(i)p(i): Correspondingly the residuals r(i) = b ¡ Ax(i) are updated as r(i) = r(i¡1) ¡ ®q(i) where q(i) = Ap(i) T. T. T. The choice ® = ®i = r(i¡1) r(i¡1) =p(i) Ap(i) minimizes r(i) A¡1 r(i) over all possible choices for ®. The search directions are updated using the residuals p(i) = r(i) + ¯ i¡1p(i¡1) T. T. where the choice ¯ i = r(i) r(i) =r (i¡1) r(i¡1) ensures that p(i) and Ap(i¡1) or equivalently, r(i) and r(i¡1) are orthogonal. In fact one can show that this choice of ¯ i makes p(i) and r(i) orthogonal to all previous Ap(j) and r(j) respectively. It uses a preconditioner M; for M = I one obtains the unpreconditioned version of the Conjugate Gradient Algorithm. In that case the algorithm may be further simpli…ed by skipping the ”solve” line and replacing z (i¡1) by r (i¡1).. 30.

(38) 7. Numerical Results and Conclusions. 7.1. Hsu’s Equation. 7.1.1. A.(Basic Atoms):. *. ZI 1 1 1 "ª0 (r) = ¡ 52 ª0(r) ¡ ª0 (r) + ª0 (r) 2 jrj 2 jr ¡ r0 j Atom Z Energy Exp error He 2 78.16 79 1.06% + Li 3 195.79 198 1.11% +2 Be 4 368.80 371 0.59% B+3 5 595.95 600 0.68% C+4 6 877.05 882 0.56% Atom Z Energy Hsu Error Iteration He 2 78.16 78.63 0.60% 4 Li 3 197.87 201.12 1.62% 4 Be 4 386.31 394.24 2.01% 5 B 5 647.71 661.72 2.11% 7 C 6 984.38 1013.63 2.89% 7 N 7 1406.58 1458.69 3.57% 8 O 8 1929.05 2003.00 3.69% 8. +. (37). Time Exp 00:28 79 00:40 203.48 02:08 399.14 03:47 670.96 03:50 1030.08 07:48 1486.02 07:51 2043.74. a. The standard Density Functional approach requires 20 or more self-consistency iterations to reach the ground state, but Hsu’s formulation spends us only 4-8 self-consistency iterations. b. For total energy point of view, Hsu’s formulation provides good numerical results. c. When Z=5~10 , the type of P orbit must separate into Px ,Py ,Pz orbits. d. The contours of 3 orbitals for B are shown in Appendix. 7.1.2. B.(A Simple H2+ molecular model): Z1 Z2 1 ª0 (r) ¡ ª0 (r) "ª0(r) = ¡ 52 ª0 (r) ¡ 2 jr ¡ R1 j jr ¡ R2 j. (38). Z1 Z2 Distance Nucleus energy Total energy Exp (-1,0,0) (1,0,0) 2 13.6eV -16.32eV -16.3eV p (-1,-1,0) (1,1,0) 2p2 9.62eV -15.79eV x (-1,-1,-1) (1,1,1) 2 3 7.85eV -15.21eV x a. During the molecular H2+ computation, the most important advantage is to use XYZ coordinates on locating the nucleus. 31.

(39) b. Because this case only has 1 electron, results are very close to experimental values. c. When the distance of two nucleus is 2 , the energy of H2+ is in the most stable situation. 7.1.3. C.(Be with excited states): X ZI 1 ª2S (r1 ) "ª2S (r1) = ¡ r2 ª2S (r1 ) ¡ 2 jr1 ¡ R1 j Z Z jª2P (r2 )j2 jª1S (r2 )j2 1 d¿ 2 + d¿ 2 ª2S (r1) + ª2s (r1 ) jr1 ¡ r2 j 2 jr1 ¡ r2 j Z ª¤ (r2)ª2S (r2 ) 1 d¿ 2 2P ª2P (r1 ) ¨ 2 jr1 ¡ r2 j. (39). Energy Exp Error Time Be 386.31 399.14eV 3.21% 02:06 Excited Be 381.98 4.33eV 04:32 Triplet State 384.47 1.84eV 2.72eV 32.35% 04:32 Singlet State 379.49 6.82eV 5.28eV 29.17% 04:32 a. Hsu’s new approach Equation admits excited states, but traditional DFT doesn’t. b. We only present ¯ = 0 version of Hsu’s formulation, and the results are still room for improvement.. 7.2. Hartree’s Results:. *. 1 ZI 1 "ª0 (r) = ¡ 52 ª0 (r) ¡ ª0 (r) + (N ¡ 1)ª0(r) 2 jrj jr ¡ r0 j. +. Atom Z Energy Exp Error Iteration Time Hsu He 2 78.04 79 1.22% 6 00:45 78.63 Li 3 203.32 203.48 0.08% 6 01:11 201.12 Be 4 396.67 399.14 0.62% 8 03:14 394.24 B 5 664.83 670.96 0.91% 17 15:00 661.72 a. As Z is more and more large, and it takes more iterations and time. b. Numerical results are very good on computing single atom.. 7.3. Kohn Sham Equation: 1 (¡ 52 +Vef f )ªi (r) = "i ªi (r); 2 32. (40).

(40) Veff = Ve¡nuc + Ve¡e +. ±Exc ; ±½. Z 3 6 4=3 Exc [½(r)] = ¡ ( )1=3 (½4=3 " (r) + ½# (r))dr 4 ¼. Atom Kin Vie Vee Vxc E He 77.49 -182.64 55.17 -23.87 -73.85. (41). LSDA Exp Error -74.09 0.32%. a. We review Kohn Sham Equation’s results in order to compare with Hsu’s results.. 7.4. Virial theorem Check. He Kin V Ratio Hsu 84.97 -161.31 1.898 Hartree 77.39 -155.43 2.008 Kohn Sham 77.48 -151.34 1.953 a. Obviously, Hartree’s result obeyes the virial theorem.. 7.5. 3D Finite Element Codes. a. During our DFT program, the Fortran program is composed of about 2000 lines. volume of the largest element b. ¼ 10000: volume of the smallest element c. For controlling memory e¢ciently, we use the method of random-pack-storage that only records the value of the nonzero entries of matrices. d. Now, we mesh usually 117079 nodes(1740293 nonzero terms) or 282639 nodes(4204067 nonzero terms) on the domain in one single machine(CPU:2.4G ;Ram:1G). e. The Gauss Cubic is a powerful numerical integration that can simplify a complicated integral and obtain a very good approximation. f. Mesh Generation Tool Atom He Ansys HyperMesh Node 18974 14576 Energy 66.57 76.31 We observe that quality of grids in HyperMesh is better than in Ansys, but Ansys tool can extend linear element to quadratic element.. 33.

(41) HyperMesh. 34.

(42) Ansys. 35.

(43) References [1] R.L. Burden and J. D. Faires, Numerical Analysis, Brooks/Cole Publishing Company, 1997. [2] J.Y. Hsu, Derivation of the density functional theory from the cluster expansion, Phys. Rev. Lett., 91:133001, 2003. [3] J.Y.Hsu, C.H.Lin, C.Y.Cheng, Pair correlation e¤ect on the self consistent density functional theory, Phys. Rev. A, 71 052502, 2005. [4] S.Gilbert and F.George, An Analysis of The Finite Element Method, Englewood Cli¤s, NJ/Prentice-Hall, 1973. [5] W.Kohn and Sham, Phys Rev A, 140:1133, 1965. [6] Markus Oppel, Density functional theory,2002. [7] T.L.Beck, Real-space mesh techniques in density functional theory, Reviews of Modern Physics, vol. 72, no. 4, October 2000. [8] G.L.Sleijpen and H.A.van der Vorst. A Jacobin-Davidson iteration method for linear eigenvalue problems. SIAM J. Sci. Matrix Anal. Appl., 17(2):401425, April, 1996. [9] O.C.Zienkiewicz and R.L.Taylor, The Finite Element Method: The Basis, Butterworth-Heinemann, 2000. [10] T.A.Arias, Multiresolution analysis of electronic structure: semicardinal and wavelet bases, Rev. Mod. Phys., 71:267, 1999. [11] T.M.Hwang, Numerical experiences of large-scale eigenvalue problems with Jacobi-Davidson method, 2003. [12] W.Wang, T.M.Hwang, W.W.Lin, and J.L.Liu, Numerical methods for semiconductor heterostructures with band nonparabolicity, J. Comp. Phys., 190:141-158, 2003. [13] J.C.Cuevas, Introduction to Density Functional Theory,www-tfp.physik.unikarlsruhe.de/~cuevas. [14] T.A.Arias, Notes on the ab initio theory of molecules and solids: Density functional theory (DFT). [15] J. Kohano¤ and N.I. Gidopoulos, Density Functional Theory: Basics, New Trends and Applications,Volume 2, Part 5, Chapter 26, pp 532-568 in Handbook of Molecular Physics and Quantum Chemistry. [16] X.M. Tong and S.I.Chu, Density-functional theory with optimized e¤ective potential and self-interaction correction for ground states and autoionizing resonances,volume 55, number 5,Phys.Rev.A. [17] Y.Saad, Numerical Methods in Electronic Structures Calculations, Beijing, Nov.1st, 2004. [18] J.Y.Hsu, Triplet and Singlet in the scDFT. [19] H.A.van der Vorst, Computational Methods for Large Eigenvalue Problems. [20] S.S. RAO, The Finite Element Method In Engineering, 4th, 1982. [21] I.Vasiliev, Serdar, and J.R.Chelikowsky, First-principles density-functional 36.

(44) calaulations for optical spectra of clusters and nanocrystals,Phys.Rev.B, volume 65,115416. [22] J.L.Liu and J.H.Chen, An Introduction to Jacobi-Davidson Method, 2005. [23] A.Goswami, Quantum Mechanics, Wm. C. Brown Publisher, 1992.. 37.

(45) APPENDIX. Atom B eigenvector1. 38.

(46) Atom B eigenvector1. 39.








(54) The ‡ow chart of the FEDFT.. 47.

(55) Jacobi-Davidson algorithm. 48.

(56) the Preconditioned Conjugate Gradient Method. 49.

(57)