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.
APPENDIX
Fermi-Dirac Function
Part II
3D Finite Element Method of Self-consistent Density
Functional Theory
3 Introduction
Density functional theory is an extremely successful approach for the descrip-tion 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 func-tional theory by HSU [HSU,2003] using the 3D …nite element method. The new formulation that drops the exchange-correlation term makes the computation po-tentially 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 consider-ing Sconsider-inglet 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
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 calcu-lating 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 sus-tain 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 den-sity 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 exter-nal 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
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 man-ner. 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 trace-able 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 computa-tional 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 func-tional theory (DFT); real-space methods are in no way limited to DFT, but since DFT calculations comprise a dominant theme in modern electrostatics and elec-tronic 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 com-plex 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 [Kim-ball 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
signi…cantly greater ‡exibility in the construction of the mesh and an underlying variational-type formulation. Other advantages include easier parallel implemen-tation 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 multi-scale or related solution methods [Thomas, 2000].Nevertheless, we mainly discuss the fundamentals of FEM solutions of the self-consistent nonlinear coupled Pois-son and Schr Äodinger (PS) equations in density-functional theory (DFT) in this thesis.
5 Density Functional Theory
The Hamiltonian in its original form is very complex:
H =¡h2 which involves sums over all electrons/nuclei and their pairs involving kinetic energy and columbic potentials.
5.1 Hartree Approximation
Hartree approximation consists of postulating that the many-electron wave func-tion can be written as a simple product of one-electron wave funcfunc-tions. 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 antisym-metry principle, which states that a wavefunction describing fermions should be antisymmetric with respect to the interchange of any set of space-spin coordi-nates.
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 ex-plicitly 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 ap-proximation (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 The Hartree approximation treats the electrons as distinguishable particles.
A step forward is to introduce Pauli exclusion principle (Fermi statistics for elec-trons) by proposing an antisymmetrized many-electron wave function in the form of a Slater determinant:
©(r; R) = SDf'j(ri; ¾i)g = 1
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 rea-sonable 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
252+Vef f)ªi(r) = "iªi(r); (3) where the density-dependent e¤ective potential is
Vef f(r) = Ve¡nuc(r) + Ve¡e([½(r)]) + Vxc([½(r)]; r), (4) Ve¡nuc = external potential, Ve¡e = R jr¡r½(r00)jdr0, Vxc = ±E±½xc , Exc =exchange-correlation enery.
Local Spin Density Approximation (LSDA) : exchange-correlation energy Exc
is a simple known funtion
Exc[½(r)] =¡3 4(6
¼)1=3
Z
(½4=3" (r) + ½4=3# (r))dr: (5) 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 poten-tial Vxc([½(r)]; r) incorporates all nonclassical e¤ects. The exchange-correlation potential includes a kinetic contribution, since the expectation value of the Kohn-Sham 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 :
5.3.1 Self Consistency 1.
(¡1
252+Ve¡nuc(r)+Ve¡e([½(r)])+Vxc([½(r)]; r))ªi(r) = "iªi(r); i = 1; :::; ioccup:: 2.
½(r) =
occup:X
i
jªi(r)j2: 3.
52Ve¡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.
Slove (¡1
252+Veff)ªi(r) = "iªi(r) 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
For example, we want to solve Li atom of con…guration 1s22s1 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
252+Vef f;")ªi;"= "i;"ªi;" (8) for spin up orbitals and replace " by # for spin down orbitals.We can have the form of Vxc;" = ±E[½±½";½#]" . The correlation functional is the same as exchange 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,
ª(¡) = ¦Ni=1©(ri)U(¡) (9)
¡ 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
¦Ni=2d¿ij©(ri)j2jU(r1; r2; :::; rN)©(r)j2 ´ M(r)j©(r)j2 ´ jª0(r)j2; (10)
where 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
"ªj(r) =¡1 ª0 constant to minimize the total energy, subject to particale conservation
N =
Z
d¿1d¿2jº(r1; r2)j2½(r1)½(r2)´
Z
d¿1d¿2jª(r1; r2)j2: (14) We reduce to the simpli…ed ¯ = 0 version
"ªs(¡!r1) =¡1
Example for excited Be(2 1S,1 2S,1 2P)
"ª2S(r1) =¡1
6 NUMERICAL METHODS 6.1 Finite Element Method
6.1.1 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 devel-oped 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 ap-plied 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 writ-ten, 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 ar-eas 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.
tetrahedral element
By the software (HyperMesh), we can partition the domain.
The meshes of 3D computation domain( r=20 Bohr)
Shape Functions The typical 3-D linear tetrahedral element trial solution can be written
U(r; ®) =e X4 j=1
®jNj(x; y; z) The special coordinates are introduced de…ned by
x = L1x1+ L2x2+ L3x3+ L4x4 (8.a) y = L1y1+ L2y2+ L3y3+ L4y4 (8.b) z = L1z1+ L2z2+ L3z3 + L4z4 (8.c)
1 = L1+ L2+ L3+ L4 (8.d)
Sloving Eq(8) gives Li = ai+ bix + ciy + diz
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
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
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.
Using a trial solution U(r; ®) to approximate the density function ª(r), wee
…rst write down the residual equation for each orbital J : RJ(r; ®) =¡1
2 52U(r; ®)e ¡ Z
rU (r; ®) + ¦e JU(r; ®)e ¡ "U(r; ®)e (19) The typical 3-D element trial solution is always written in the general form U(r; ®) =e Pnj=1®jNj(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.
=)
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:
1Z Z Z @Ni(r)@U(r; ®)e
+ @Ni(r)@U (r; ®)e
+@Ni(r)@U (r; ®)e
dV (23)
STEP 3.
Equation (23) is rearranged to satisfy the matrix solver as follows :
f" [M] + [N]g f®g = 0 (24)
To transfer the integrals of Equation(24) into a form appropriate for numerical evaluation.
r dV (Gauss Cubic integration) (30)
=Vc
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 .
Numerical integration formula of Gauss for tetrahedral element [Zienkiewicz and Taylor, 20 [C] =¡
Z Z Z
¦jNi(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) Boundary Condition :
=) Deriving the element equations complete by above steps. The Conjugate Gradient Method can slovehAi[X] =hBi, and it means the ¦j term is known now. We get [C]ij =¡R R R ¦jNi(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 Numerical Linear Algebra
6.2.1 Jacobi-Davidson Method
The Jacobi-Davidson solver can e¢ciently deal with the large-scale sparse eigen-value matrix equation, which is still a very challenging task even nowadays. One of the advantages in using Jacobi-Davidson algorithm to solve the DFT eigen-value 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; :::; vmg. The Galerkin condition is
AVms¡ µVms? fv1; :::; vmg
where Vm denotes the matrix with columns v1 to vm. This equation has m solutions (µ(m)j ; s(m)j ). 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 u(m)j .
A(u(m)j + t) = ¸(u(m)j + t) with t ? (u(m)j ), t2 U?´ fvjv¤uj = 0g.
The correction equation of Jacobi-Davidson is
? (36)
where rj = (A¡ µjI)uj, A?= (I ¡ uju¤j)A(I¡ uju¤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 iter-ation if the correction equiter-ation 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¤AiV .
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 meth-ods. (Nonstationary methods di¤er from stationary methods in that the compu-tations involve information that changes at each iteration. Typically , constants are computed by taking inner products of residuals or other vectors arising from
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 meth-ods. (Nonstationary methods di¤er from stationary methods in that the compu-tations involve information that changes at each iteration. Typically , constants are computed by taking inner products of residuals or other vectors arising from