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.