Accurate finite element method for atomic calculations
based on density functional theory and HartreeFock method
Taisuke Ozaki (JAIST) and Masayuki Toyoda (JAIST)
Copyright 2011 Elsevier B.V. This article is provided by the author for the reader's personal use only.
Any other use requires prior permission of the author and Elsevier B.V.
NOTICE: This is the author's version of a work accepted for publication by Elsevier. Changes resulting
from the publishing process, including peer review, editing, corrections, structural formatting
and other quality control mechanisms, may not be reflected in this document. Changes may have been
made to this work since it was submitted for publication. A definitive version is subsequently
published in Computer Physics Communications 182, 12451252 (2011),
DOI:10.1016/j.cpc.2011.02.010 and may be found at
here.
Abstract:
An accurate finite element method is developed for atomic calculations based
on density functional theory (DFT) within local density approximation (LDA) and
HartreeFock (HF) method.
The radial wave functions are expanded by cubic Hermite spline functions
on a uniform mesh for
, and all the associated integrals are
analytically evaluated in conjunction with fitting procedures of the hartree
and the exchangecorrelation potentials to the same cubic Hermite spline
functions using a set of recurrence formulas. The total energy of atoms
systematically converges from above, and the error algebraically decays
as the mesh spacing decreases.
When the mesh spacing
is taken to be
bohr
, the total
energy for an atom of atomic number
can be calculated within error of
hartree for both the LDA and HF methods.
The equal applicability of the method to DFT and the HF method with
a similar degree of high accuracy enables the method to be a reliable platform
for development of new functionals in DFT such as hybrid functionals.
The electronic structure calculation of atoms[1,2] is
one of most fundamental bases in not only understanding electronic structures
of molecules and solids, but also developing efficient and accurate electronic
structure methods. For the latter case, it is indispensable to distinguish
the intrinsic error produced by the theoretical framework itself
from that caused by the other numerical problems such as incompleteness
of basis set and inaccurate numerical integration,[3]
which will be referred to as nonintrinsic error hereafter.
If the nonintrinsic error is extremely small, the validation of
electronic structure methods can be very precisely performed, which
may highlight strength and weakness of each method without suffering from
any appreciable numerical error.
By the recent advance of the electronic structure methods, requirement
for the intrinsic error has been approaching chemical
accuracy (1mHartree).[4,5,6]
One can imagine that the requirement for the intrinsic error should be
smaller than the chemical accuracy for further development of
accurate electronic structure methods.
A potential approach which can minimize the nonintrinsic error
is a finite element (FE) method in which wave functions are expressed by
a linear combination of piecewise polynomial functions.
[7,8,10,9,11,12,13,14,15,16]
Since the approach based on the FE method can be regarded as a traditional
basis set method, once a highly accurate FE method is established,
it is apparent for the FE method to be widely used as a platform for
development of new functionals in density functional
theory (DFT)[17,18] and post HartreFock (HF) methods.
The equal applicability of the method to DFT and the HF method with a similar
degree of high accuracy is also highly important because hybrid methods, in which DFT
and wave function theories such as the HF method are unified in a single framework,
have recently attracted much attention.[6,19,20,21,22,23]
Although a wide variety of FE methods have been already
developed for atomic calculations so far,[7,8,10,9,11,12]
the equal applicability of the method to DFs and wave function
theories has not been clearly established. In DFT calculations numerical integrations
have to be essentially employed due to nonanalytic nature of exchangecorrelation
functionals, while two electron repulsion integrals can be calculated analytically
in wave function theories such as the HF method. The difference in evaluating matrix
elements for exchangecorrelation potentials limits the equal applicability of
each method to DFs and wave function theories with a similar degree of high accuracy.
In this paper, in order to establish a basis for development of
electronic structure methods which may consist of a DF and a wave function
theory, we present an accurate finite element method for atomic calculations
based on not only DFT, but also the HF method. To address the equal applicability of the
method to DFs and wave function theories,
we use one of the simplest piecewise polynomial functions, i.e., cubic Hermite
spline functions, as basis functions, and show that the use of the cubic
Hermite spline functions allows one to analytically evaluate all of integrals
involved in conjunction with fitting procedures of the hartree and the
exchangecorrelation potentials to the same cubic Hermite spline functions.
As a result of the analytic evaluation of matrix elements, it is found that
the nonintrinsic error in the DFT and HF calculations can be systematically
reduced by only decreasing
the mesh spacing, and that eventually the FE method can provide numerically exact
solutions within machine precision for all the atoms (Z=1103) in the periodic table.
This paper is organized as follows:
In Sec. II the theoretical framework of the FE method using the cubic Hermite
spline functions is presented for DFT within local density approximation (LDA)
and the HF method. In Sec. III numerical accuracy of the FE method is demonstrated
by the total energy energy calculations of all the atoms (Z=1103) in the periodic table.
In Sec. IV we summarize the FE method, and discuss a future possible
application of the method in development of a new functional in DFT.
Throughout the paper, we use atomic units, where
.
Let us start to write the one particle wave function
as



(1) 
where and are a spherical harmonic and a radial wave function,
respectively. By introducing change of variables ,[24]
which allows us to eliminate the cusp of wave functions at the origin in
coordinate, and assuming a spherical potential , the radial
Schrödinger equation is given by



(2) 
with



(3) 
where the operators and are defined by
and the potential is the sum of three contributions:



(6) 
The first term is the attractive potential of the nucleus with atomic number ,
and the second and third are
the hartree and exchange correlation potentials, respectively.
Although we restrict ourselves to the spherical potential in the paper,
the FE method we discuss may be generalized to the nonspherical case by combining
with the Slater method.[2]
We now expand the radial function using cubic Hermite spline functions
as basis functions[28],
which are placed on a regular mesh with spacing in coordinate,
as follows:



(7) 
where the spacing is determined by
, and
and , shown in Fig. 1(a), are defined by
The spline function satisfies a set of conditions that
, , and
,
while satisfies , , and
.
The superscript for in Eq. (7) means that the origin of
is shifted to , where is given by .
Figure 1:
(Color online)
(a) Cubic Hermite spline functions and defined by Eqs. (8) and (9).
(b) Two sorts of regular meshes and in coordinate, which are
used for fitting of the hartree and exchangecorrelation potentials.

The accuracy in description of can be systematically controlled by increasing
as shown later. Employing Eq. (7) for Eq. (2) readily gives the following generalized
eigenvalue problem:



(10) 
where and are the hamiltonian and overlap matrices, respectively,
and these matrices becomes septuple diagonal in LDA, since basis functions
interact with basis functions in the nearest neighbor
sites .
All the matrix elements, except for that of the exchangecorrelation potential
such as LDA,[25,26] can be analytically evaluated.
For example, all the nonzero matrix elements for the operator are
found to be
offsite elements
onsite elements for
onsite elements for
It turns out that the matrix elements depend on the site as a result that the extent
of region spanned by the basis function in coordinate varies depending on the site .
Also it is noted that four onsite matrix elements at have to be
calculated by taking account of a fact that the basis functions span only
the region ranging from to in coordinate, resulting in Eqs. (18)(20).
As well as , the matrix elements for , , and
and the overlap matrix elements can be analytically evaluated, and all the analytic
formulas and MATHEMATICA codes for generating them are provided in
Secs. S112 of the supplemental material.
Although the matrix elements for can be analytically evaluated as shown
in the section of HF method, for the LDA calculation we present
an alternative numerical
method which is much faster than the analytic counterpart, while keeping numerical accuracy.
Thus, one may see that in LDA the remaining matrix elements which are numerically evaluated
are only those for the hartree and exchangecorrelation potentials in the FE method.
Here we show that even the matrix elements for the hartree and exchangecorrelation
potentials can be very accurately evaluated by utilizing the same cubic spline functions
in the FE method.
Let us introduce two sorts of regular meshes with spacing in coordinate,
that is, one is as defined before and the other
,
as shown in Fig. 1(b). We evaluate the charge density on the two regular meshes,
where, by recalling that the basis functions are strictly localized in real space,
one can compute the charge density by considering only contributions of the neighboring sites at
each mesh point, and fit the charge density on the meshes to the following function:



(21) 
where and are uniquely determined by a set of recurrence formulas:
The recurrence formulas are derived by noting that only is nonzero at
as shown in Fig. 1(b) and that is only the unknown parameter if we start to
fit
to the function Eq. (21) from , which results in
that Eq. (23) has to be recursively solved starting from .
Once is expanded in terms of the Hermite spline functions, by considering the
spherical charge density distribution, we can analytically
evaluate the hartree potential as
where the integrals are analytically evaluated due to the integrands by simple
polynomial functions. As well,
on the other mesh are
calculated in the same way as for
. The explicit formulas of the
integrals can be found in Secs. S13 and 14 of the supplemental material.
Using the similar way to the expansion of charge density,
the set of
and
can be fitted to
where the coefficients and are uniquely determined
by the similar recurrence formulas to Eqs. (22) and (23).
Once
and
are fitted to Eq. (25),
the calculation of matrix elements for is straightforward as follows:




(23) 
The integrals involved survive only if , or for ,
and or for
, and there are 24 and 16 nonzero
elements for the former and the latter, respectively, which can be easily evaluated
in analytic formulas.
The other combinations give always zero elements due to the strictly
localized spline functions. All the analytic formulas and MATHEMATICA codes
for generating them are provided in Secs. S1521 of the supplemental material.
The matrix elements for the exchangecorrelation potential can be calculated by just
replacing with in Eqs. (25) and (26) after
and
are calculated by LDA.
Thus, from the above derivation we see that all the matrix elements required in the
FE method within LDA are evaluated without introducing numerical integration which
can be a serious source of numerical error.
It is also pointed out that the extension of the method to generalized gradient
approximation (GGA)[27] has no difficulty. We only have to perform the evaluation of
by GGA instead of LDA.
Since the resultant hamiltonian and overlap matrices are septuple diagonal,
the eigenvalues and eigenstates can be efficiently calculated by a combination of
a shiftandinverse Lanczos method and a shiftinvert method,[29]
which are used to estimate approximate eigenstates and to refine the approximate
eigenstates, respectively.
To calculate approximate eigenstates by the shiftandinverse Lanczos method,
the generalized eigenvalue problem Eq. (10) is now transformed to



(24) 
where the transformed hamiltonian , eigenvalues , and eigenvectors
are given by
The shift
for the eigenvalues is taken to be an approximate lowest
eigenvalue for each angular momentum , which can be found from the results at
the previous selfconsistent field (SCF) step.
Then, the Lanczos iteration for Eq. (27) is performed by the following recurrence formulas:
where the initial vector
is generated by random numbers, but
normalized with respect to .
The multiplication between a matrix and a vector,
and
, can be calculated by making use of
the LU and Cholesky factorizations, respectively, in O() operations.
The recursion level of 100200 is used to obtain sufficient convergence.
By diagonalizing the tridiagonal matrix of which diagonal elements are and
the offdiagonal elements are , and backtransforming
and by Eqs. (29)(30), one can obtain a set of approximate
eigenvalues
and eigenvectors starting from
the lowest state.
The approximate eigenvalues
and eigenvectors
obtained by the shiftandinverse Lanczos method can be further refined
by the following shiftinvert method:
where one of the approximate vectors, corresponding to an occupied state,
is chosen as the initial vector
. The iteration is
repeated until a condition
is satisfied.
Only less than 10 iterations are enough to achieve the condition for
each eigenstate. It is also noted that the shiftandinverse Lanczos
method can be skipped as the SCF iteration converges, since the eigenvectors
found at the previous SCF step are good approximation for the shiftinvert
method at the next SCF step, which allows us to accelerate the calculation.
The difference between the LDA and HF methods lies in treatment of the exchange potential.
In the HF method, the following potential is used instead of Eq. (6)

(35) 
where is the nonlocal exchange operator which acts as

(36) 
One can analytically evaluate matrix elements for the nonlocal exchange potential
thanks to the simple polynomial of the Hermite spline functions. In addition,
the matrix elements for the hartree potential are also analytically evaluated
to keep consistency in evaluating the matrix elements for the hartree and exchange
potentials in the HF method, while the matrix elements are evaluated in
an alternative way in the LDA calculation.
To evaluate the exchange potential with the one particle wave functions Eq. (1),
the Coulomb operator is expanded in terms of the spherical harmonics



(37) 
where and are the greater and lesser of and ,
respectively.
By using the expansion of the radial function Eq. (7) along with Eq.
40, the exchange energy is computed as follows:
where is the spherically averaged density matrix for each
channel

(39) 
and
is dependent matrix elements of the exchange potential
in a representation with the Hermite spline basis functions.
Similarly, the Hartree energy is computed as



(40) 
with independent matrix elements
.
After performing the integration in the angular coordinates,
the matrix elements of and , which are
parts of the hamiltonian matrix elements, are obtained as

(41) 
and
where the coefficient
is
and the quantity denoted by a closed bracket is the twoelectron
integrals
The factor of a half in Eq. 45 appears because
degenerate spin configurations are assumed here.
Note also in Eq. 45 that the summation over is
truncated due to the fact that the coefficient
is always zero when or
.
The integral 47 is invariant under the
following rotations of indices
Due to this invariance, the ordering among the mesh indices

(46) 
can be assumed without losing generality.
Furthermore, by considering that the integral has a nonzero value
only when and specify the same mesh point or the neighboring
points with each other, and so are and , it is also possible
to assume

(47) 
To summarize above, one can safely assume that

(48) 
and

(49) 
This is validated because the integral is always zero whenever Eqs.
51 and 52 cannot be satisfied
simultaneously under any rotation in Eq. 48.
Then, within this assumption, the integration ranges for and
in Eq. 47 are



(50) 



(51) 
These ranges may and may not overlap with each other.
The simpler case is when or, equivalently,
and thus they have no overlap.
In this case, the integral is given as



(52) 
where

(53) 
Note that the integration range is bounded to be nonnegative.
Therefore, the actual lower limit of the above integral is
, although we denote it as for simplicity.
Readers must keep in mind that the similar notations are also used
in the following equations.
For the other case where or, equivalently,
, the ranges have overlap with each
other for

(54) 
Then, the integral is given as
where
All the integrals , , and in Eqs. 56,
59, 60 and 61, respectively, can be
analytically solved.
For example, the integral with are:
When ,
When ,
When ,
It turns out that the combinations of the mesh indices which gives
nonzero values of are classified into the following cases
and otherwise is zero.
All the analytic formulas for , , , and and MATHEMATICA codes for generating
them are provided in Sec. S22 of the supplemental material.
Since the exchange term is nonlocal, the hamiltonian is not septuple diagonal.
Therefore, one cannot apply the same techniques as the LDA calculations to solve
the generalized eigenvalue problem. However, by noting that the septuple diagonal
elements in the hamiltonian are still dominant, the refinement procesure can be
generalized even for the dense hamiltonian matrix in the HF method. To derive
the generalized method, we first divide and into
and
,
where is the septuple diagonal hamiltonian and is a reference energy.
By putting these expressions into Eq. (10), one can derive the following equation:
Based on the equation, the shiftinvert method of Eqs. (35)(37) can be generalized as follows:
where one of the approximate vectors for an occupied state is chosen as the initial vector
. As well as the shiftinvert method by Eqs. (35)(37), one can
achieve sufficient convergence by only less than 10 iterations by Eqs. (73)(77).
The matrix multiplication with
is performed
by making use of the LU factorization in O() operations as in the LDA calculation.
We employ the conventional scheme for the diagonalization in the initial stage
of the SCF calculation, and switch it to the above shiftinvert method after several
SCF iterations, which accelerates the calculation, since a few eigenstates only
have to be evaluated in the scheme.
The numerical accuracy of the solution for the Schrödinger equation can be
evaluated by the virial theorem. If the solution is exact, the virial theorem
rigorously holds. Thus, the numerical deviation in the virial theorem is
a measure of examining numerical error of the solution.
Considering that the correlation energy in LDA includes a part of kinetic energy,
the virial theorem for LDA is defined by



(75) 
where and are the conventional kinetic
and potential energies in LDA. The contribution from the correlation
energy to the kinetic energy,
, is given by



(76) 
with the definition of an energy density



(77) 
where
is the correlation energy density, and
.
The expressions, Eqs. (78)(80), for the virial theorem can be derived
as shown in Ref. [30] by using the generalized procedure
by Slater.[31]
On the other hand, the virial theorem simply holds in the HF method
without any correction.
Table 1:
The virial theorem and the total energy in hartree calculated
by the LDA and HF methods for the ground state of a helium atom.
of 10 bohr is used as the maximum range of coordinate.
The bold font means that the number is exact.
Meshes 


Total energy 


LDA 

10 
0.023649134368226 

2.709633955981526 
20 
0.000766923670941 

2.835007522850960 
40 
0.000040592237419 

2.834807057468094 
80 
0.000000722804997 

2.834835146173011 
160 
0.000000011744304 

2.834835616626474 
320 
0.000000000190544 

2.834835623943877 
640 
0.000000000003354 

2.834835624053601 
1280 
0.000000000000075 

2.834835624055133 


HF 

10 
0.043770457887893 

2.847096711441144 
20 
0.000890014514660 

2.861255456882009 
40 
0.000024238408751 

2.861671203112089 
80 
0.000000459034967 

2.861679838221988 
160 
0.000000007541142 

2.861679993078111 
320 
0.000000000118509 

2.861679995572653 
640 
0.000000000001846 

2.861679995611623 
1280 
0.000000000000028 

2.861679995612229 
In Table I we show the convergence of the virial theorem and the total energy
for the LDA and HF calculations of a helium atom as a function of the number of meshes,
where of 10 bohr is used as the maximum range of coordinate
for all the cases. The analytic functional form parametrized by Vosko, Wilk, and
Nusair[26] is used for the LDA calculations.
All the calculations in the study are performed using long double of 80 bit
which has 19 significant digits decimally.
It is found that the errors in the virial theorem and the total energy for
the both cases algebraically decay as the number of meshes increases.
Also we see that the order of the errors in the virial theorem and the total energy
are almost equivalent to each other, which supports that the evaluation of the virial
theorem can be a measure of checking the accuracy of the total energy.
Using 1280 mesh points, corresponding to bohr,
the total energy is calculated with accuracy of 14 digits for the LDA
and HF calculations.[32]
It is worth mentioning that the total energy converges from
above as the number of meshes increases for both the LDA and HF calculations, which
indicates that the method can be regarded as a variational method in practice.
We further verify the accuracy of the method by applying the FE method
to all the atoms (Z=1103) in the periodic table within LDA and
a series of rare gas atoms within the HF method,[33]
where the spherical charge density distribution is assumed for the
nonspin polarized ground state electronic configuration in the LDA calculations.
Figure 2 shows the absolute value of the virial theorem,
,
as a function of atomic numbers. By considering that the eigenenergies of hydrogen
like atoms scales as , the mesh spacing are taken to be inversely
proportional to so that the bare Coulomb potential at
can be proportional to .
The error in the virial theorem for the LDA calculations with the mesh spacing
of is
and
hartree
for hydrogen and lawrencium atoms, respectively,
and the errors of the other cases are in between those of the two atoms, which
suggests that using the mesh spacing the total energy for LDA can be computed
within error of hartree for all the atoms in the periodic table.
In addition to
,
we also check a variant of the virial theorem
(not shown), and
find that
is about for all the atoms, indicating that
the relative error is almost constant and the number of accurate digits of the
total energy is almost equivalent to each other, which is 1314 digits in the case
of the mesh spacing of . It is also confirmed that all the calculated
results by LDA are consistent with the results by Kotochigova et al.[34]
The error in the virial theorem in the HF calculations with the mesh spacing of
varies in a similar way to that of the LDA calculations.
The same analysis as the LDA case implies that the total
energy of the HF calculations can be obtained within error of hartree
and that the number of accurate digits of the total energy is 1112 digits for all
the rare gas atoms in the case of the mesh spacing of
.
The LDA calculations with the coarse mesh spacing are also presented for comparison,
showing that the error in the HF calculation is comparable to the corresponding LDA
calculation if the same grid spacing is used. The comparison leads to another
conclusion that the numerical fitting of exchangecorrelation and hartree potentials
in the LDA calculations is not a source of numerical errors. In both the LDA and HF
calculations the error in the total energy mainly comes from expansion of
the wave functions by the finite basis functions.
Figure:
(Color online)
The absolute value (hartree) of the virial theorem,
,
as a function of atomic numbers calculated by the LDA and HF methods,
and corresponding curves of error estimated by Eq. (85).
The nonspin polarized ground state electronic configuration is considered for
all the cases, and
the spherical charge density distribution is assumed in the LDA calculations.
The mesh spacing is taken to be ,
,
bohr.
of 10 bohr is used as the maximum range of coordinate
for all the cases.

We now derive a formula which gives the absolute error in the total energy calculated
by the FE method.
In this derivation it is assumed that a large is used so that
the truncation of tail of wave functions cannot be a source of numerical error.
The calculation in Fig. 2 suggests that of 10 bohr is large
enough to avoid the error for all the elements (Z=1103) in the periodic table.
From the two observations that the error decreases algebraically as the mesh spacing
decreases in Table I, and that the number of accurate digits of the total energy is
nearly constant when the mesh spacing is taken to be inversely proportional to
, the number of accurate digits of the total energy,
, may be written as



(78) 
where and are parameters to be fitted, and is the mesh spacing as
discussed before. Even if we have the same number of accurate digits of the
total energy for different elements, the absolute error in the total energy depends
on the absolute magnitude of the total energy. Therefore, as the next step, let us roughly
estimate the total energy of atoms.
Suppose that the total energy can be estimated by the sum of eigenvalues of hydrogen like atoms.
Then, the total energy of an atom in which electrons fully occupy all the states up to the
principal quantum number is obtained by



(79) 
In the case, the atomic number is calculated by
which gives
as increases. Substituting the asymptotic form
for Eq. (82), we have



(81) 
Noting that the number of accurate decimal places is given by
,
and using Eqs. (81) and (84), one can derive the following expression to estimate the absolute
error, , in the total energy,
where , , and are found to be 6, 2, and 0.3 by fitting Eq. (85) to the
data shown in Fig. 2. As shown in Fig. 2 the estimated error by Eq. (85) fits well with
that in the whole range we study.
Strictly speaking, the expression of Eq. (85) estimates the error in the the virial theorem,
since the fitting is done using the virial theorem
. However, the expression
can be used to estimate the error in the total energy due to the near equivalence of
the error between the virial theorem and the total energy.
The expression implies that the error in the total energy approximately scales as
or when or are fixed.
As well as the convergence property of the total energy, the eigenvalue of state
and the charge density at the origin in the ground state of a helium atom are shown as
a function of the number of meshes in Table II. Although these quantities slowly converge
compared to the total energy shown in Table I, it can be seen that the systematic
convergence produces highly accurate results for both the LDA and HF calculations.
Table 2:
The eigenvalue (hartree) of 1s state and charge density at the origin in the
ground state of a helium atom.
of 10 bohr is used as the maximum range of coordinate.
The bold font means that the number is exact.
Meshes 
Eigenvalue of 

at the origin 


LDA 

10 
0.471608882934828 

4.5260991473798461147 
20 
0.570859970667500 

3.5721180605164599170 
40 
0.570412462140008 

3.5237334305245833802 
80 
0.570424629988695 

3.5258950648066678425 
160 
0.570424750048491 

3.5267674469147037484 
320 
0.570424730001095 

3.5268447250668697409 
640 
0.570424724176759 

3.5268499263432255490 
1280 
0.570424722705590 

3.5268502577049796584 


HF 

10 
0.917364393558815 

4.6752954006527842707 
20 
0.917885224607414 

3.6462035714690724058 
40 
0.917953892393717 

3.5931129688148539896 
80 
0.917955530938917 

3.5949601670283128317 
160 
0.917955562327678 

3.5958316993332024767 
320 
0.917955562848009 

3.5959123561764878905 
640 
0.917955562856210 

3.5959178873966144104 
1280 
0.917955562856337 

3.5959182401606711429 
We have developed an accurate FE method for atomic calculations based
on DFT within LDA and the HartreeFock method. Cubic Hermite spline
functions on a uniform mesh for are used as basis functions to expand the radial
wave functions. The numerical integrations being a source of numerical error can be
avoided due to the simplicity of the cubic Hermite spline functions, and all the
associated integrals are analytically evaluated in conjunction with fitting of
the hartree and exchangecorrelation potentials to the same cubic Hermite spline functions.
By taking account of the localized nature of the basis functions in real space,
the generalized eigenvalue problem is efficiently solved using a generalized shiftinvert
iterative method for not only the LDA but also HF calculations.
The numerical calculations show that the convergence is systematically
controlled by the mesh spacing and in practice numerically exact solutions can be
obtained within machine precision for all the elements (Z=1103) in the periodic table.
The convergence of the total energy from above implies that the FE method can
be regarded as a variational scheme with respect to
the Hermite spline functions for both the LDA and HF methods.
Based on the virial theorem and an intuitive analysis the absolute error in the total energy
is estimated to be proportional to .
The absolute error is less than hartree for the LDA and HF calculations when
the mesh spacing is taken to be
bohr.
Since the FE method can provide high quality numerical solutions
with a similar degree of accuracy for both DFT and the HF method, this can be
utilized as a platform, which is free from
the nonintrinsic numerical error in the validation of newly developed methods,
for development of new functionals in DFT such as hybrid functionals.
Along this line, we have been trying to develop a hybrid exchange hole model
based on the FE method, which will be discussed elsewhere.
ACKNOWLEDGMENT
The authors were partly supported by the Fujitsu lab., the Nissan Motor Co., Ltd.,
Nippon Sheet Glass Co., Ltd., and the Next Generation Super Computing Project,
Nanoscience Program, MEXT, Japan.

 1

J.C. Slater,
Quantum Theory of Atomic Structure, Vol. 1, McGrawHill, New York, 1960.
 2

J.C. Slater,
Quantum Theory of Atomic Structure, Vol. 2, McGrawHill, New York, 1960.
 3

D. Feller, K.A. Peterson, and T.D. Crawford,
J. Chem. Phys. 124 (2006) 054107.
 4

R.M. Martin,
Electronic Structure: Basic Theory and Practical Methods,
Cambridge University Press, New York, 2008.
 5

T. Helgaker, P. Jorgensen, and J. Olsen,
Molecular ElectronicStructure Theory,
Wiley, 2000.
 6

C.E. Dykstra, G. Frenking, K.S. Kim, and G.E. Scuseria,
Theory and Applications of Computational Chemistry: The First 40 Years
(A Volume of Technical and Historical Perspectives),
Elsevier, Amsterdam (2005).
 7

B.W. Shore,
J. Chem. Phys. 58 (1973) 3855.
 8

J.L. Gázquez and H.J. Silverstone,
J. Chem. Phys. 67 (1977) 1887.
 9

C.F. Fischer and W. Guo,
J. Comp. Phys. 90 (1990) 486.
 10

H.T. Jeng and C.S. Hsue,
Chin. J. Phys. 35 (1997) 215.
 11

Z. Romanowski,
Modelling Simul. Mater. Sci. Eng. 16, 015003 (2008);
ibid. 17 (2009) 045001.
 12

D. Engel, M. Klews, and G. Wunner,
Comp. Phys. Comm. 180 (2009) 302.
 13

S.R. White, J.W. Wilkins, and M.P. Teter,
Phys. Rev. B 39 (1989) 5819.
 14

E. Tsuchida and M. Tsukada, Phys. Rev. B 54 (1996) 7602.
 15

J.E. Pask and P. Sterne,
Modelling Simul. Mater. Sci. Eng. 13 (2005) R71.
 16

E.J. Bylaska, M. Holst, and J.H. Weare,
J. Chem. Theory Comput. 5 (2009) 937.
 17

P. Hohenberg and W. Kohn, Phys. Rev. 136 (1964) B864.
 18

W. Kohn and L. J. Sham, Phys. Rev. 140 (1965) A1133.
 19

A.D. Becke, J. Chem. Phys. 98 (1993) 1372.
 20

T. Leininger, H. Stoll, H.J. Werner, and A. Savin, Chem. Phys. Lett. 275 (1997) 151.
 21

H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115 (2001) 3540.
 22

J. Heyd, G.E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118 (2003) 8207.
 23

T.M. Henderson, A.F. Izmaylov, G. Scalmani, and G.E. Scuseria,
J. Chem. Phys. 131 (2009) 044108.
 24

In addition to the change of variables, , we tried two other cases,
and , and found that the latter two cases lead to
complication in the formalism.
 25

D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45 (1980) 566.
 26

S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, (1980) 1200;
S.H. Vosko and L. Wilk, Phys. Rev. B 22, (1980) 3812.
 27

J.P. Perdew, K. Burke, and M. Ernzerhof,
Phys. Rev. Lett. 77 (1996) 3865.
 28

In addition to the cubic Hermite spline functions, we investigated
the convergence of the quintic Hermite spline functions, and found that
the convergence ratio with respect to the total number of basis functions is
comparable to the cubic case. Thus, it is concluded that the cubic Hermite spline
functions is an optimum choice with respect to the convergence ratio and
the simplicity of analytic expressions derived for matrix elements.
 29

Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst,
Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide,
Society for Industrial Mathematics (1987).
 30

F.W. Averill and G.S. Painter,
Phys. Rev. B 24 (1981) 6795.
 31

J.C. Slater,
J. Chem. Phys. 57 (1972) 2389.
 32

The total energy we obtain corresponds to that in the nonrelativistic HF limit,
which lacks the correlation energy of 0.042044 hartree
compared to the exact energy of 2.903724 hartree.
 33

The program code, ADPACK, and the calculated results are
available on a web site (http://www.openmxsquare.org/).
 34

S. Kotochigova, Z.H. Levine, E.L. Shirley, M.D. Stiles, and C.W. Clark,
Phys. Rev. A 55, 191 (1997); a web site of the database
(http://www.nist.gov/physlab/data/dftdata/)
20110301