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 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
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)
The choice ® = ®i = r(i¡1)Tr(i¡1)=p(i)TAp(i) minimizes r(i)TA¡1r(i) over all possi-ble choices for ®. The search directions are updated using the residuals
p(i)= r(i) + ¯i¡1p(i¡1)
where the choice ¯i = r(i)Tr(i)=r(i¡1)Tr(i¡1) ensures that p(i) and Ap(i¡1) or equiv-alently, 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).
7 Numerical Results and Conclusions
Atom Z Energy Exp errorHe 2 78.16 79 1.06%
Li+ 3 195.79 198 1.11%
Be+2 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 Time Exp
He 2 78.16 78.63 0.60% 4 00:28 79
Li 3 197.87 201.12 1.62% 4 00:40 203.48 Be 4 386.31 394.24 2.01% 5 02:08 399.14
B 5 647.71 661.72 2.11% 7 03:47 670.96
C 6 984.38 1013.63 2.89% 7 03:50 1030.08 N 7 1406.58 1458.69 3.57% 8 07:48 1486.02 O 8 1929.05 2003.00 3.69% 8 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):
"ª0(r) =¡1
2 52ª0(r)¡ Z1
jr ¡ R1jª0(r)¡ Z2
jr ¡ R2jª0(r) (38)
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):
"ª2S(r1) =¡1
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. 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+V )ª(r) = "ª(r);
Veff = Ve¡nuc+ Ve¡e+ ±Exc
Atom Kin Vie Vee Vxc E LSDA Exp Error
He 77.49 -182.64 55.17 -23.87 -73.85 -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.
b. volume of the largest element
volume of the smallest element ¼10000:
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
HyperMesh
Ansys
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 expan-sion, 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, Engle-wood 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):401-425, 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.uni-karlsruhe.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 Hand-book 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 res-onances,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
calaulations for optical spectra of clusters and nanocrystals,Phys.Rev.B, vol-ume 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.
APPENDIX
Atom B eigenvector1
Atom B eigenvector1
Atom B eigenvector1
Atom B eigenvector2
Atom B eigenvector2
Atom B eigenvector2
Atom B eigenvector3
Atom B eigenvector3
Atom B eigenvector3
Jacobi-Davidson algorithm
the Preconditioned Conjugate Gradient Method