Taisuke Ozaki (JAIST) and Hiori Kino (NIMS)
The electronic structure method based on a linear combination of numerical
atomic local basis orbitals (LCAO)[1,2,3,4,5] and
pseudo potentials [6,7] is a possible way of extending the
applicability of a density functional theory (DFT) [8] to
large-scale systems, since the generalized eigenvalue problem with
the resultant sparse matrices can be solved by O()
methods.[9,10,11,12,13,14] The pioneering work by
Sankey et al.[1] and its generalization to a fully
self-consistent field (SCF) calculation by the SIESTA group [2]
enable us to deal with large-scale systems within a fully self-consistent DFT.
However, in general, one of difficult problems in the implementation
of these LCAO methods is how the matrix elements for the effective potential
are accurately evaluated within a modest
computational effort.[2,3,15,16,17,18]
Although this numerical integration for the construction of matrix
elements is a highly technical issue, serious problems enough to
lose the validity of calculations are often elicited when we
calculate weakly interacting systems such as hydrogen bonding
molecules and include deep semicore states in pseudo potentials.[2]
Due to the non-integrable form of the effective potential in the
conventional DFT, the numerical integration is indispensable to
construct the matrix elements, even if analytic basis functions
such as Gaussian functions are used. While a scheme to avoid the
numerical integration is to fit the potential into integrable
analytic functions,[17,20] however, it can be regarded
as just a replacement of the numerical integration with the numerical fitting.
In the LCAO methods, usually, two kinds of numerical grids are used for the numerical integration for the construction of matrix elements.[2,3,15,16,17,18,19] One of them is a grid decomposed into radial and angular parts such as Gauss-Legendre and Lebedev grids which are combined with a partitioning scheme to decompose a multicenter integral into one-center integrals.[3,15,16,17,18,19] The other is a regular real-space grid similar to that used in conventional plane-wave DFT methods.[2] Although the former possesses a benefit that the total energy is independent of the rotation of system unlike the latter, the latter is superior to the former in terms of computational efficiency. Since all the basis functions use the same regular real-space grid, once the value of basis function on grid is calculated and stored in the computational memory before the SCF calculation, there is no need of the recalculation of the value of basis function on grid during the SCF calculation. Therefore, in the latter, the matrix elements are quite efficiently constructed in an element-by-element fashion for non-zero values of basis functions on grid stored in the computational memory. On the other hand, in the former, it is difficult to store the value of basis function on grid because of the requirement of the extensive memory size, since the quadrature grid is determined for each set of two atoms.
Although the regular real-space grid possesses an advantage in terms of computational efficiency, however, serious problems remain in the grid integration using the regular real-space grid, since the total energy depends on the relative position between the system and the grid. [2,21,22] A typical illustration is that the optimized structure depends on the initial arrangement relative to the grid position. This dependence causes a severe problem when weakly interacting systems such as hydrogen bonding molecules are calculated. In such a system a highly fine grid is needed to obtain a convergent result, which demands great computational resources. The same problem involved in the real-space grid integration is also reported in other real-space methods,[21,22] which implies that the difficulty involved in the grid integration is a major concern in the efficient and accurate implementation of real-space methods. Thus, our aim is to establish a method which addresses the regular real-space grid for the simplicity and overcomes the difficulty involved in the grid integration within a modest computational effort. In this paper, we present a projector expansion method in which a rapidly varying local potential in the vicinity of nuclei is expressed by a separable projector expansion, and thereby the contribution of the rapidly varying local potential to matrix elements is accurately calculated in the momentum space instead of the real space. Since the real-space integration is applied to only slowly varying parts in the effective potential, as a result, the difficulty involved in the grid integration is remarkably reduced. This paper is organized as follows. In Sec. II, a projector expansion method is presented in order to overcome the difficulty involved in the grid integration in the first-principles electronic structure calculations using atomic basis functions. In Sec. III, as illustrations of the accuracy of the proposed method, it is shown that the projector expansion method significantly suppresses not only a spurious oscillation in the energy curve for the atomic displacement involved in a weak interaction such as hydrogen bonding, but also the dependence of optimized structure on relative position to the real-space grid in the geometry optimization within a modest grid fineness. In Sec. IV, we summarize the projector expansion method and the capability. Finally, in the Appendix, the analytic evaluation of force on atom in the LCAO method using the projector expansion method is briefly described.
Let us start our formulation from the total energy expression
of the DFT [8] to elucidate the difficulty involved in
the grid integration. In the framework of atomic basis functions
,[1,2,4,5] norm-conserving
pseudo potentials,[6,7] and a local density
approximation (LDA) [23] to the exchange-correlation term,
the total energy
in the DFT is given by the sum of
the kinetic energy
, the electron-core Coulomb energy
,
the electron-electron Coulomb energy
, the exchange-correlation
energy
, and the core-core Coulomb energy
between
pseudo core charges
and
as follows:
![]() |
(1) |
![]() |
(2) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
(3) |
![]() |
![]() |
![]() |
(4) |
![]() |
(5) |
![]() |
(6) |
![]() |
(7) |
![]() |
(8) |
![]() |
(9) |
![]() |
(10) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
(11) |
![]() |
(12) |
Rather than the matrix elements expressible by two-center integrals,
the difficulty involved in the grid integration comes from the remaining
contributions in the effective potential. Although the matrix elements
for the remaining potentials
,
,
and
has been evaluated by the grid integration using
the regular real-space grid,[2] it can be pointed out
that the matrix element for the neutral atom potential
is difficult to be accurately calculated within a modest fineness of the
real-space grid. In Fig. 1 shows potentials
,
,
and
in Eq. (12) along the bond axis in a carbon monoxide (CO).
We can see that the neutral atom potential
rapidly varies
in the vicinity of nuclei, while the other potentials
and
smoothly vary. The comparison suggests that a finer grid
is required to accurately evaluate the matrix elements for the neutral
atom potential
, and that the calculations are not converged
when a modest fineness of the real-space grid is used as shown later on.
Thus, in order to suppress the difficulty involved in the real-space
grid integration, we propose a projector expansion method to accurately
evaluate the matrix elements for the neutral atom potential .
The neutral atom potential
is spherical and is defined
within the finite range determined by the cutoff radius
of
the confinement potential. Therefore, the potential
can be expressed by a projector expansion as follows:
![]() |
(13) |
![]() |
![]() |
![]() |
(14) |
![]() |
![]() |
![]() |
(15) |
There exists another advantage of the projector expansion method. The choice of the local part in the separable pseudo potential is an important factor to enhance the transferability of the pseudo potential. However, choosing a deep local part has tended to be avoided because finer real-space mesh is necessary to calculate converged results in the real-space grid integration. In this projector expansion method, there is no difficulty in the choice of the deep local part. Thus, it is possible not only to choose the deep local part, but also to include deep semicore states in the pseudo potential without increasing computational effort.
If the same radial functions introduced in Eq. (14) are
used for the basis functions
(
)
in the LCAO method, the matrix elements
for the neutral atom potential are evaluated in case of
and/or
without depending on the convergence of the summation over the angular and
radial parts.
Since
,
e.g., in case of
, we have
![]() |
![]() |
![]() |
|
![]() |
![]() |
(16) |
Although we are able to complete our formulation for the projector
expansion here, a different scheme is furthermore introduced to
accurately calculate the two-center integral ().
In this case, the neutral atom potential
and the
product of
and
are Fourier-transformed
at the first stage of the calculation, and therefore the
two-center integral between the neutral atom potential
and the product of two basis functions is accurately evaluated
in the momentum space apart from the projector expansion by Eq. (13).
As a result, all the one- and two-center integrals, which are
large components in the matrix elements, are calculated with a high
degree of precision without depending on the convergence of the
summation over the angular and radial parts. The projector expansion
by Eq. (13) is applied only for the three-center integrals
which are small components in the matrix elements. Thus, the error
associated with the convergence of the summation in the projector
expansion is significantly suppressed.
There is another origin of the difficulty involved in the
real-space grid integration. The partial core correction (PCC) charge
density is often used in Eq. (5) to take into
account the non-linearity in the exchange-correlation terms.
In this case, if a highly localized and large PCC charge density
in the vicinity of nuclei is employed, the difficulty associated
with Eq. (5) appears in the real-space grid integration.
However, the difficulty can be avoidable by using a modest
PCC charge density. In addition, we can avoid the inclusion of
the large PCC charge density by using the pseudo potential including
the semicore states, since it is relatively easy to include the
semicore states in the projector expansion method as discussed above.
Here, we would like to comment on the evaluation of force on atom
in the LCAO method using the projector expansion scheme.
The majority contributions to the total energy are expressed by
the two-center integrals in the projector expansion method.
So, the contribution to the force associated with these parts can
be analytically calculated by differentiating the two-center
integrals evaluated in the momentum space without any difficulty.
On the other hand, the remaining and
are evaluated using the real-space grid integration.
Consequently, it is not trivial whether the force associated
with
and
is analytically appraisable
or not. However, in fact, it is easy to analytically evaluate
the force on atom in the regular real-space scheme coupled with
the projector expansion method. The details are given in the Appendix.
In this section, we show numerical results on the convergence of the summation over the angular and radial parts in Eq. (13) and demonstrate the capability of the projector expansion method for the suppression of the difficulty involved in the real-space grid integration.
![]() |
In Fig. 2 shows the convergence properties of the total energy in
a water molecule as a function of the max ,
, and
the number of radial projectors,
, in the summation
in Eq. (13). In this study, an optimized double valence plus
polarization function (DVP) was used as a basis set.[4,5]
To replace the deep core potential into a shallow potential,
a norm-conserving pseudo potential [7] was used in a separable
form with multiple projectors.[26] The cutoff radii of pseudo
potentials are listed in Table I of Ref.[5]. A generalized
gradient approximation (GGA) [24] is used for the
exchange-correlation without the nonlinear partial core correction.
The cutoff energy of 3600 (Ryd) and the Gauss-Legendre grid [27]
of 128 were used for the evaluation of the two-center integrals (integrals)
in the momentum space. The cutoff energy for the real-space grid integration
is given in the caption of the Figures and Table. All calculations were
performed by our DFT code, OpenMX.[28] From the convergence
properties, we see that the fully convergent result is achieved
when
and
approach 6 and 4, respectively.
In this case, we have to use at least 2 for
and
to accurately evaluate the one- and two-center integrals,
since the double valence plus polarization function is used as
a basis set. If the values are less than 2, the accuracy of
the one- and two-center integrals is not assured.
Based on the convergence properties, we generally use the
following values:
and
,
where
is the maximum angular momentum
quantum number of basis functions. From several test calculations,
we confirm that the values for
and
are
enough to achieve the fully convergent results for systems including
other elements, while the results are not shown.
To illustrate the accuracy of the projector expansion method,
we show a comparison of the total energies of a water dimer molecule,
which is a typical hydrogen bonding system with a weak interaction,
calculated by the projector expansion method (Projector) and the
grid integration using the regular real-space grid (Non-Projector)
in Fig. 3. The total energies are calculated as a function of the
distance between two oxygen atoms in water molecules, each of which
possesses the experimental geometry. The dihedral angle between the
water molecules is fixed at 180 degrees as depicted in Fig. 3.
We see that the total energy calculated by the grid integration
using the regular real-space grid with a lower cutoff energy
considerably oscillates, while the oscillation can be suppressed
using the high cutoff energy of 1000 Ryd.
This spurious oscillation is due to the inaccurate calculation
of the matrix elements for the neutral atom potential around
the tail of basis functions in the real-space grid integration.
On the other hand, the projector expansion method entirely suppresses
the spurious oscillation of the total energy curve even in the lowest
cutoff energy, thus indicating that the contribution
around the tail of basis functions to the matrix elements
associated with the neutral atom potential are accurately calculated.
The calculated energy curve by the projector method is completely
consistent with that by the non-projector method using 1000 Ryd.
It should be noted that the error in the calculation of the total
charge is only 1.0e-5 even in the non-projector method using 200 Ryd.
Nevertheless, the spurious oscillation appears in amplitude of
0.001 Hartree in the non-projector method. Our analysis clearly
shows that the spurious oscillation comes from the deep neutral
atom potential , and is removed by replacing it
by the projector expansion.
![]() |
As shown above, the spurious oscillation appears in the energy curve calculated by the real-space grid integration. Therefore, it would be considered that the optimized structure is trapped in a local minimum near the initial geometry if the geometry optimization is performed using the real-space grid integration, which could be a serious problem in the calculations of weakly interacting systems. In contrast, the difficulty is overcome in the projector expansion method. To study this dependence of the optimized structure on the initial geometry, we optimize the structure of a water molecule with various initial positions relative to the regular real-space grid. The optimized structure by the real-space grid integration strongly depends on the initial geometry relative to the regular real-space grid as shown in Table I, which indicates that the structure is trapped in a local minimum near the initial geometry. On the other hand, the optimized structure by the projector expansion method is almost independent of the initial geometry relative to the regular real-space grid, while the difference in the bond angle is 0.4 degrees at a maximum. These two clearly illustrate that the difficulty involved in the grid integration is successfully suppressed by introducing the projector expansion for the rapidly varying neutral atom potential in the vicinity of nuclei.
Projector | |||
Initial Geo. | ![]() |
![]() |
![]() |
a | 0.990, 0.990 | 104.8 | 1.95 |
b | 0.990, 0.990 | 104.5 | 1.95 |
c | 0.990, 0.990 | 104.8 | 1.94 |
d | 0.990, 0.990 | 104.9 | 1.95 |
Non-Projector | |||
Initial Geo. | ![]() |
![]() |
![]() |
a | 0.989, 0.989 | 101.2 | 2.00 |
b | 0.982, 0.982 | 102.0 | 1.98 |
c | 0.993, 0.993 | 101.4 | 1.99 |
d | 1.004, 0.993 | 107.0 | 1.92 |
Expt. | 0.957 | 104.5 | 1.86 |
In the real-space implementation of the DFT using atomic basis functions, one of difficult problems is how the matrix elements for the effective potential are accurately evaluated in an efficient numerical scheme within a modest computation effort. Especially, this problem becomes obvious when weakly interacting systems such as hydrogen bonding molecules are calculated. In this context, to suppress the difficulty involved in the regular real-space grid integration, we have developed a projector expansion method which expresses the rapidly varying neutral atom potential in the vicinity of nuclei by a separable projector. The neutral atom potential being the origin of the difficulty involved in the real-space grid integration is factorized by the separable projectors, and the resultant two-center integrals are accurately evaluated in the momentum space instead of the real space. To illustrate the accuracy of the projector expansion method, we have shown that the spurious oscillation in the energy curve of a water dimer molecule is entirely suppressed and that the optimized structure is almost independent of the initial geometry relative to the regular real-space grid within a modest grid fineness. Thus, we conclude that the projector expansion method is a highly useful technique for the accurate and efficient implementation of the DFT based on the LCAO method.
This work was partly supported by NEDO under the Nanotechnology Materials Program, Research and Development for Applying Advanced Computational Science and Technology of Japan Science and Technology Corporation (ACT-JST), and NAREGI Nanoscience Project and Special Coordination Funds of the Ministry of Education, Culture, Sports, Science and Technology, Japan.
In this Appendix, we show an explicit expression for the analytic
evaluation of the force on atom in the projector expansion method.
On the parts of the total energy expressed by the two-center
integrals (integrals), its contributions to the force are
analytically calculated by differentiating the two-center
integrals (integrals) evaluated in the momentum space.
Then, each term in the summation over the discretized radial
grid in the momentum space are analytically differentiating
using cubic splines without any difficulty.
Therefore, we focus on the analytic differentiation of the
remaining
and
with respect
to atomic coordinate. So, this derivation is a general formulation
for the analytic force evaluation in the LCAO method based on the
grid integration using the regular real-space grid.
The difference electron-electron Coulomb energy
and the exchange-correlation energy
can be discretized
using the regular real-space grid as follows:
![]() |
![]() |
![]() |
(17) |
![]() |
(18) |
![]() |
(19) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
(20) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
||
![]() |
![]() |
(21) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
||
![]() |
![]() |
(22) |
![]() |
![]() |
![]() |
|
![]() |
(23) |
![]() |
![]() |
![]() |
(24) |
![]() |
![]() |
![]() |
|
![]() |
![]() |
(25) |
The force on atom in the LCAO methods based on the numerical
grid integration has been often evaluated by appealing the
Hellman-Feymann theorem with the Pulay correction [3] rather
than strictly differentiating the discretized energy expression.
However, our derivation shows that it is possible to evaluate
the analytic force consistent with the total energy in an efficient way,
even though both the numerical basis functions and the numerical
grid integration are introduced. A defect of the analytic forces
derived here is that they do not fulfill the law of action and reaction.
The defect comes from the independent determination of the regular
real-space grid from the atomic positions. Since the grid position
is fixed in the absolute Cartesian coordinate, the calculated force
suffers from 'a pinning effect'. Nevertheless, the effect is relatively
small since only the minor contributions
and
to the total energy are evaluated using the real-space grid integration.
A way of recovering the law of action and reaction is to introduce
an adaptive coordinate [31] which is dependent on the atomic
positions, while the simplicity in the regular real-space grid is lost.
A study toward this direction will be a future work.