Efficient projector expansion in LCAO electronic structure calculations: Ver. 1.0
The original has been published in Phys. Rev. B 72, 045121 (2005).

Abstract

A projector expansion method is presented for an efficient and accurate implementation of the first-principles electronic structure calculations using pseudo potentials and atomic basis functions. By expressing the rapidly varying local potential in the vicinity of nuclei by a separable projector expansion, the difficulty involved in the grid integration using the regular real-space grid is remarkably reduced without increasing the computational effort. To illustrate the capability, it is shown that the proposed 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.

I. INTRODUCTION

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($N$) 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.

II. PROJECTOR EXPANSION

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 $\{{\chi }_{i\alpha }\}$,[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 $E_{\mathrm{tot}}$ in the DFT is given by the sum of the kinetic energy $E_{\mathrm{kin}}$, the electron-core Coulomb energy $E_{\mathrm{ec}}$, the electron-electron Coulomb energy $E_{\mathrm{ee}}$, the exchange-correlation energy $E_{\mathrm{xc}}$, and the core-core Coulomb energy $E_{\mathrm{cc}}$ between pseudo core charges $Z_i$ and $Z_j$ as follows:

$E_{\mathrm{tot}}=E_{\mathrm{kin}}+E_{\mathrm{ec}}+E_{\mathrm{ee}}+E_{\mathrm{xc}}+E_{\mathrm{cc}}$ (1)

with

$E_{\mathrm{kin}}=2{\displaystyle \sum _{i\alpha ,j\beta }}{\rho }_{i\alpha ,j\beta }{\displaystyle \int d(𝐫){\chi }_{i\alpha }\widehat{T}{\chi }_{j\beta }},$ (2)

$E_{\mathrm{ec}}$ $=$ $E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ec}}^{(\mathrm{NL})}$ (3)

$=$ ${\displaystyle \int d(𝐫)n(𝐫)\sum _k{V}_{\mathrm{L},k}(𝐫-{𝐑}_k)}+{\displaystyle \int d(𝐫)n(𝐫)\sum _k{V}_{\mathrm{NL},k}(𝐫-{𝐑}_k)},$

$E_{\mathrm{ee}}$ $=$ ${\displaystyle \frac{1}{2}}{\displaystyle \int \int d(𝐫)d({𝐫}^{\prime })\frac{n(𝐫)n({𝐫}^{\prime })}{|𝐫-{𝐫}^{\prime }|}},$ (4)

$E_{\mathrm{xc}}={\displaystyle \int d(𝐫)n{ϵ}_{\mathrm{xc}}(n)},$ (5)

$E_{\mathrm{cc}}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j}}{\displaystyle \frac{Z_i{Z}_j}{|{𝐑}_i-{𝐑}_j|}},$ (6)

where $i$ and $\alpha$ are the site and basis function indices, respectively. ${\rho }_{i\alpha ,j\beta }$ is a density matrix associated with two basis functions ${\chi }_{i\alpha }$ and ${\chi }_{j\beta }$, and is defined by ${\sum }_{\nu }\mathrm{\Theta }({\epsilon }_{\nu }-\mu )c_{i\alpha ,\nu }{c}_{j\beta ,\nu }$ with the LCAO coefficient $c_{i\alpha ,\nu }$, the one-particle eigenenergy ${\epsilon }_{\nu }$, the chemical potential $\mu$, and a step function $\mathrm{\Theta }(x)$. $n$ is the electron density defined by $2{\sum }_{i\alpha ,j\beta }{\rho }_{i\alpha ,j\beta }{\chi }_{i\alpha }{\chi }_{j\beta }$. $V_{\mathrm{L},k}$ and $V_{\mathrm{NL},k}$ are the local part and non-local part in the norm-conserving pseudo potential of atom $k$, respectively. The factor 2 in $E_{\mathrm{kin}}$ and $n$ is for the spin multiplicity. In this formulation, we consider only non spin-polarized case, non-Bloch expression of one-particle wave functions, and the LDA for simplicity, but the extensions of the below description to the spin polarized case, the Bloch expression, and a generalized gradient approximation (GGA) [24] to the exchange-correlation term are straightforward. Although the expression Eq. (1) is commonly used, it can be transformed into a more tractable form without any approximation.[1, 2] Reorganizing the sum of $E_{\mathrm{ec}}^{(\mathrm{L})}$, $E_{\mathrm{ee}}$ and $E_{\mathrm{cc}}$ by introducing an atomic electron density $n_i^{(a)}$, we can rewrite the sum $E_{\mathrm{ec0}}+E_{\mathrm{ee}}+E_{\mathrm{cc}}$ by that of two majority contributions $E_{\mathrm{na}}$ and $E_{\mathrm{ecc}}$ consisting of short-range terms and a minority contribution $E_{\delta \mathrm{ee}}$ consisting of a long-range term as follows:

$E_{\mathrm{ec}}^{(\mathrm{L})}+E_{\mathrm{ee}}+E_{\mathrm{cc}}=E_{\mathrm{na}}+E_{\mathrm{scc}}+E_{\delta \mathrm{ee}}$ (7)

with

$E_{\mathrm{na}}={\displaystyle \int d(𝐫)n(𝐫)\sum _i{V}_{\mathrm{na},i}(𝐫-{𝐑}_i)},$ (8)

$E_{\mathrm{scc}}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j}}\left[{\displaystyle \frac{Z_i{Z}_j}{|{𝐑}_i-{𝐑}_j|}}-{\displaystyle \int 𝑑𝐫n_i^{(a)}(𝐫-{𝐑}_i)V_{\mathrm{H},j}^{(a)}(𝐫-{𝐑}_j)}\right],$ (9)

$E_{\delta \mathrm{ee}}={\displaystyle \frac{1}{2}}{\displaystyle \int 𝑑𝐫\delta n(𝐫)\delta V_{\mathrm{H}}(𝐫)},$ (10)

where a difference electron density $\delta n(𝐫)$ is defined by the difference between the electron density $n$ and the sum, $n^{(a)}$, of atomic electron densities $n_i^{(a)}$ as:

$\delta n(𝐫)$ $=$ $n(𝐫)-n^{(a)}(𝐫),$ (11)

$=$ $n(𝐫)-{\displaystyle \sum _i}{n}_i^{(a)}(𝐫).$

Through this paper, when the site index is dropped in the designation for a quantity specifiable with the site index, it means the inclusion of the summation over the site index except for the symbolic use. The atomic electron density $n_i^{(a)}$ is evaluated from pseudo wave functions under the confinement potential in this study.[4, 5] So, $\delta n$ is considerably smaller than $n^{(a)}$. Also, it should be noted that $n_i^{(a)}$ is finite only within the confinement radius $r_{c,i}$. From the sum of the local part of pseudo potential $V_{\mathrm{L},i}$ and a Hartree potential $V_{\mathrm{H},i}^{(\mathrm{a})}$ associated with the atomic electron density, a neutral atom potential $V_{\mathrm{na},i}$ is defined by $V_{\mathrm{na},i}=V_{\mathrm{L},i}+V_{\mathrm{H},i}^{(a)}$. The neutral atom potential is spherical because of the spherical atomic electron density, and becomes zero beyond the confinement radius $r_{c,i}$ due to Gauss’s law, indicating that it is a short-range spherical potential. Since the electron density $n$ can be evaluated the short-range quantities $\chi$ and $\rho$, we see that the neutral atom potential energy $E_{\mathrm{na}}$ consists of short-range terms. In addition to this, the screened core-core Coulomb energy $E_{\mathrm{scc}}$ is also evaluated by taking account of only the neighboring atoms, since the second term in Eq. (9) becomes exactly equivalent to the first term when $r_{c,i}+r_{c,j}\le |{𝐑}_i-{𝐑}_j|$, and thereby the long-range terms in Eq. (9) vanish. Therefore, the first and second terms being the majority contributions in the right side of Eq. (7) are calculated by making use of information of adjacent atoms in an O($N$) operation. Since the second term in the right side of Eq.  (9) is a function of only the distance between two atoms, it is accurately evaluated using a finer real-space grid, and can be tabulated as a function of the distance between two atoms at the first stage of the calculation. So, the evaluation of $E_{\mathrm{scc}}$ is far from the difficulty involved in the real-space grid integration. The difference electron-electron Coulomb energy $E_{\delta \mathrm{ee}}$ is a minority contribution to the total energy, but is constructed by a long-range term $\delta V_{\mathrm{H}}$, since the difference Hartree potential $\delta V_{\mathrm{H}}$ is associated with the every difference electron density $\delta n(𝐫)$ in the real space. It is worth mentioning that the Ewald summation over the core-core Coulomb energy is not required by the reorganization for parts of the total energy, which is one of advantages in the reorganized expression for the total energy. Considering the variation of the total energy with respect to the LCAO coefficient $c_{i\alpha ,\nu }$, we have a well-known generalized eigenvalue equation $H{c}_{\nu }={\epsilon }_{\nu }S{c}_{\nu }$, where $H$ and $S$ are a Hamiltonian matrix defined by $⟨{\chi }_{i\alpha }|\widehat{H}|{\chi }_{j\beta }⟩$ and an overlap matrix defined by $⟨{\chi }_{i\alpha }|{\chi }_{j\beta }⟩$, respectively. In the Kohn-Sham (KS) Hamiltonian $\widehat{H}\equiv \widehat{T}+V_{\mathrm{eff}}$, the effective potential $V_{\mathrm{eff}}$ is given by

$V_{\mathrm{eff}}={\displaystyle \sum _k}{V}_{\mathrm{NL},k}+{\displaystyle \sum _k}{V}_{\mathrm{na},k}+\delta V_{\mathrm{H}}+V_{\mathrm{xc}},$ (12)

where $V_{\mathrm{xc}}$ is the exchange-correlation potential. In the matrix generalized eigenvalue equation, the matrix elements for the overlap matrix, the kinetic operator, and the non-local part of pseudo potential $V_{\mathrm{NL},k}$ can be evaluated by two-center integrals. Since $V_{\mathrm{NL},k}$ is expressed by a separable projector expansion,[25, 26] the matrix elements are reduced to a product of two-center integrals. The resultant two-center integrals in the evaluation of these matrix elements can be very accurately evaluated using a finer grid in the momentum space within almost the same computational time because they only have to be evaluated once before the SCF loop.[1, 2]

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 $V_{\mathrm{na},k}$, $\delta V_{\mathrm{H}}$, and $V_{\mathrm{xc}}$ 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 $V_{\mathrm{na},k}$ is difficult to be accurately calculated within a modest fineness of the real-space grid. In Fig. 1 shows potentials $V_{\mathrm{na}}$, $\delta V_{\mathrm{H}}$, and $V_{\mathrm{xc}}$ in Eq. (12) along the bond axis in a carbon monoxide (CO). We can see that the neutral atom potential $V_{\mathrm{na}}$ rapidly varies in the vicinity of nuclei, while the other potentials $\delta V_{\mathrm{H}}$ and $V_{\mathrm{xc}}$ smoothly vary. The comparison suggests that a finer grid is required to accurately evaluate the matrix elements for the neutral atom potential $V_{\mathrm{na}}$, and that the calculations are not converged when a modest fineness of the real-space grid is used as shown later on.

Figure 1: Potentials $V_{\mathrm{na}}$, $\delta V_{\mathrm{H}}$, and $V_{\mathrm{xc}}$ in Eq. (12) along the bond axis in a carbon monoxide (CO).

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 $V_{\mathrm{na}}$. The neutral atom potential $V_{\mathrm{na},i}$ is spherical and is defined within the finite range determined by the cutoff radius $r_{c,i}$ of the confinement potential. Therefore, the potential $V_{\mathrm{na},i}$ can be expressed by a projector expansion as follows:

${\widehat{V}}_{\mathrm{na},k}={\displaystyle \sum _{lm}^{L_{\max }}}{\displaystyle \sum _{\zeta }^{N_{\mathrm{rad}}}}|V_{\mathrm{na},k}{\overline{R}}_{l\zeta }{Y}_{lm}⟩{\displaystyle \frac{1}{c_{l\zeta }}}⟨Y_{lm}{\overline{R}}_{l\zeta }{V}_{\mathrm{na},k}|,$ (13)

where a set of radial functions $\{{\overline{R}}_{l\zeta }\}$ is an orthonormal set defined by a norm $\int r^2𝑑rR{V}_{\mathrm{na},k}{R}^{\prime }$ for radial functions $R$ and $R^{\prime }$, and is calculated by the following Gram-Schmidt orthogonalization:[26]

${\overline{R}}_{l\zeta }$ $=$ $R_{l\zeta }-{\displaystyle \sum _{\eta }^{\zeta -1}}{\overline{R}}_{l\eta }{\displaystyle \frac{1}{c_{l\eta }}}{\displaystyle \int r^2𝑑r{\overline{R}}_{l\eta }{V}_{\mathrm{na},k}{R}_{l\zeta }},$ (14)

$c_{l\zeta }$ $=$ ${\displaystyle \int r^2𝑑r{\overline{R}}_{l\zeta }{V}_{\mathrm{na},k}{\overline{R}}_{l\zeta }}.$ (15)

The radial function $R_{l\zeta }$ used in this study is pseudo wave functions for both the ground and excited states under the same confinement potential as used in the calculation of the atomic electron density $n_i^{(a)}$.[4, 5] It is expected that the projector expansion rapidly converges with respect to the summation over $R$, since radial functions $R_{l\zeta }$ forms an orthonormal set within the finite range of $V_{\mathrm{na},k}$. The projector expansion for the neutral atom potential becomes exact when the summation over the angular and the radial parts goes to infinity. The most important feature in the projector expansion defined by Eq. (13) is that the neutral atom potential is expressed by a separable form, and thereby we only have to evaluate the two-center integrals to construct the matrix elements for the neutral atom potential once before the SCF loop as well as the other two-center integrals. As discussed above, the two-center integrals can be accurately evaluated in the momentum space instead of the real space. Therefore, we can avoid the grid integration using the real-space grid for the rapidly varying neutral atom potential in the vicinity of nuclei, indicating that the difficulty involved in the grid integration is almost suppressed. It should be noted that the computational time by the introduction of the projector expansion is not increased within the SCF loop, since the basis functions and the projectors in Eq. (13) are Fourier-transformed once at the first stage of calculations, and the two-center integrals between them are evaluated before the SCF loop. Furthermore, the computational demand before the SCF loop is not large. In fact, the comparison of the total computational time shows that the computational time by the projector expansion method is only 13 % longer than that by the conventional method with the same real-space grid fineness as used in the projector expansion method in a molecular dynamics (MD) simulation (100 MD steps) of a C${}_{60}$ molecule using 4 processors of Pentium 4. In addition to this, these calculations can be easily parallelized on a parallel computer due to no data-communication in their parallel computation.

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 $R_{l\zeta }$ introduced in Eq. (14) are used for the basis functions ${\chi }_{j\beta }$($\equiv Y_{l^{\prime }{m}^{\prime }}{R}_{l^{\prime }{\zeta }^{\prime }}$) in the LCAO method, the matrix elements $⟨{\chi }_{i\alpha }|{\widehat{V}}_{\mathrm{na},k}|{\chi }_{j\beta }⟩$ for the neutral atom potential are evaluated in case of $i=k$ and/or $j=k$ without depending on the convergence of the summation over the angular and radial parts. Since $R_{l^{\prime }{\zeta }^{\prime }}={\sum }_{\zeta }{a}_{\zeta }^{(l^{\prime }{\zeta }^{\prime })}{\overline{R}}_{l^{\prime }\zeta }$, e.g., in case of $j=k$, we have

${\widehat{V}}_{\mathrm{na},j}|{\chi }_{j\beta }⟩$ $=$ ${\displaystyle \sum _{\zeta }}{a}_{\zeta }^{(l^{\prime }{\zeta }^{\prime })}|V_{\mathrm{na},j}{\overline{R}}_{l^{\prime }\zeta }{Y}_{l^{\prime }{m}^{\prime }}⟩$ (16)

$=$ $V_{\mathrm{na},j}|{\chi }_{j\beta }⟩.$

Therefore, the one-center integral ($i=j=k$) and the two-center integral ($i=k$ or $j=k$) are highly accurately evaluated in the evaluation of $⟨{\chi }_{i\alpha }|{\widehat{V}}_{\mathrm{na},k}|{\chi }_{j\beta }⟩$, where $⟨{\chi }_{i\alpha }|{\widehat{V}}_{\mathrm{na},k}|{\chi }_{j\beta }⟩$ is referred to as integral to distinguish it from one- and two-center integrals in the projector expansion defined by Eq. (13).

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 ($i=j\ne k$). In this case, the neutral atom potential $V_{\mathrm{na},k}$ and the product of $R_{l\zeta }$ and $R_{l^{\prime }\eta }$ 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 $n_{\mathrm{pcc}}$ 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 $E_{\mathrm{xc}}$ and $E_{\delta \mathrm{ee}}$ are evaluated using the real-space grid integration. Consequently, it is not trivial whether the force associated with $E_{\mathrm{xc}}$ and $E_{\delta \mathrm{ee}}$ 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.

III. NUMERICAL RESULTS

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.

Figure 2: Convergence properties of of the total energy in a water molecule as a function of (a) the max L, $L_{\max }$, and (b) the number of radial projectors, $N_{\mathrm{rad}}$, in the summation of Eq. (13). In the calculations (a) and (b), $N_{\mathrm{rad}}$ and $L_{\max }$ were fixed at 4 and 6. The cutoff energy of 201 (Ryd) was used for the grid integration using the real-space grid. Non-Projector means that the real-space grid integration is employed for construction of the matrix elements for the neutral atom potential.

In Fig. 2 shows the convergence properties of the total energy in a water molecule as a function of the max $L$, $L_{\max }$, and the number of radial projectors, $N_{\mathrm{rad}}$, 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 $L_{\max }$ and $N_{\mathrm{rad}}$ approach 6 and 4, respectively. In this case, we have to use at least 2 for $L_{\max }$ and $N_{\mathrm{rad}}$ 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: $L_{\max }=L_{\max }^{\mathrm{basis}}+4$ and $N_{\mathrm{rad}}=4$, where $L_{\max }^{\mathrm{basis}}$ is the maximum angular momentum quantum number of basis functions. From several test calculations, we confirm that the values for $L_{\max }$ and $N_{\mathrm{rad}}$ 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 $V_{\mathrm{na}}$, and is removed by replacing it by the projector expansion.

Figure 3: Total energy (Hartree) of a water dimer as a function of the distance between two oxygen atoms calculated by the projector expansion method (Projector) and the grid integration using the regular real-space grid (Non-Projector) for the neutral atom potential. The value in parenthesis means the cutoff energy for the grid integration using the regular real-space grid. Each energy curve was shifted by adding a constant for ease of comparison.

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.

TABLE. I. Dependence of optimized structure (Å and degrees) and dipole moment (Debye) of a water molecule on the initial position relative to the regular real-space grid. The initial Cartesian coordinates (Å) are given by (a) O=(0.0,0.0,0.0), H${}_1$=(0.76,0.59,0.0), H${}_2$=(-0.76,0.59,0.0). The other initial Cartesian coordinates are generated by rotating the coordinates (a) by (b) $R_x(30)$, (c) $R_y(40)R_x(30)$, and (d) $R_z(50)R_y(40)R_x(30)$, where $R_x(30)$ means a rotational matrix which rotates the coordinate by 30 degrees on the x-axis. The cutoff energy of 177 (Ryd) was used for the grid integration using the regular real-space grid. The experimental values are taken from Ref. [29, 30].

Projector
Initial Geo.	$r_{\mathrm{OH}}$	$\mathrm{\angle }(\mathrm{HOH})$	${\mu }_{\mathrm{B}}$
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.	$r_{\mathrm{OH}}$	$\mathrm{\angle }(\mathrm{HOH})$	${\mu }_{\mathrm{B}}$
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

IV. CONCLUSIONS

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.

ACKNOWLEDGMENT

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.

Appendix A

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 $E_{\delta \mathrm{ee}}$ and $E_{\mathrm{xc}}$ 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 $E_{\delta \mathrm{ee}}$ and the exchange-correlation energy $E_{\mathrm{xc}}$ can be discretized using the regular real-space grid as follows:

$E_{\delta \mathrm{ee}}$ $=$ ${\displaystyle \frac{1}{2}}\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩})\delta V_{\mathrm{H}}({𝐫}_{𝐩}),$ (17)

$E_{\mathrm{xc}}=\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}n({𝐫}_{𝐩}){ϵ}_{\mathrm{xc}}(n({𝐫}_{𝐩})),$ (18)

where $𝐩$ is the index in a vector form for specifying the position in the regular real-space grid, and $\mathrm{\Delta }V$ is the volume per grid. Then, the derivative of $E_{\delta \mathrm{ee}}$ with respect to the atomic coordinate ${𝐑}_k$ is given by

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial {𝐑}_k}}={\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {𝐑}_k}}{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n({𝐫}_{𝐩})}}+{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n^{\mathrm{a}}({𝐫}_{𝐩})}{\partial {𝐑}_k}}{\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n^{\mathrm{a}}({𝐫}_{𝐩})}}.$ (19)

In this study, $\delta V_{\mathrm{H}}({𝐫}_{𝐩})$ is defined through the fast Fourier transform (FFT) of $\delta n(𝐫)$ by

$\delta V_{\mathrm{H}}({𝐫}_{𝐩})$ $=$ ${\displaystyle \sum _{𝐆}}\delta {\tilde{V}}_{\mathrm{H}}(𝐆){\mathrm{e}}^{i𝐆\cdot 𝐫},$ (20)

$=$ ${\displaystyle \frac{4\pi }{N_{\mathrm{rsg}}}}{\displaystyle \sum _{𝐆}}{\displaystyle \frac{1}{{|𝐆|}^2}}{\displaystyle \sum _{𝐩}}\delta n({𝐫}_{𝐩}){\mathrm{e}}^{i𝐆\cdot (𝐫-{𝐫}_{𝐩})},$

where $N_{\mathrm{rsg}}$ is the number of grid points in the regular real-space grid. Considering this expression of $\delta V_{\mathrm{H}}(𝐫)$, we obtain an explicit compact expression for the derivative of $E_{\delta \mathrm{ee}}$ with respect to $n({𝐫}_{𝐩})$ as follows:

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n({𝐫}_{𝐩})}}$ $=$ ${\displaystyle \frac{1}{2}}\mathrm{\Delta }V\{\delta V_{\mathrm{H}}({𝐫}_{𝐩})+{\displaystyle \sum _{𝐪}}\delta n({𝐫}_{𝐪}){\displaystyle \frac{\partial \delta V_{\mathrm{H}}({𝐫}_{𝐪})}{\partial n({𝐫}_{𝐩})}}\},$ (21)

$=$ ${\displaystyle \frac{1}{2}}\mathrm{\Delta }V\{\delta V_{\mathrm{H}}({𝐫}_{𝐩})+{\displaystyle \frac{4\pi }{N_{\mathrm{rsg}}}}{\displaystyle \sum _{𝐆}}{\displaystyle \frac{1}{{|𝐆|}^2}}{\displaystyle \sum _{𝐪}}\delta n({𝐫}_{𝐪}){\mathrm{e}}^{i𝐆\cdot ({𝐫}_{𝐪}-{𝐫}_{𝐩})}\},$

$=$ $\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐩}).$

Similarly, the derivative of $E_{\delta \mathrm{ee}}$ with respect to $n^{\mathrm{a}}({𝐫}_{𝐩})$ is given by

${\displaystyle \frac{\partial E_{\delta \mathrm{ee}}}{\partial n^{\mathrm{a}}({𝐫}_{𝐩})}}$ $=$ $-{\displaystyle \frac{1}{2}}\mathrm{\Delta }V\{\delta V_{\mathrm{H}}({𝐫}_{𝐩})-{\displaystyle \sum _{𝐪}}\delta n({𝐫}_{𝐪}){\displaystyle \frac{\partial \delta V_{\mathrm{H}}({𝐫}_{𝐪})}{\partial n^{\mathrm{a}}({𝐫}_{𝐩})}}\},$ (22)

$=$ $-{\displaystyle \frac{1}{2}}\mathrm{\Delta }V\{\delta V_{\mathrm{H}}({𝐫}_{𝐩})+{\displaystyle \frac{4\pi }{N_{\mathrm{rsg}}}}{\displaystyle \sum _{𝐆}}{\displaystyle \frac{1}{{|𝐆|}^2}}{\displaystyle \sum _{𝐪}}\delta n({𝐫}_{𝐪}){\mathrm{e}}^{i𝐆\cdot ({𝐫}_{𝐪}-{𝐫}_{𝐩})}\},$

$=$ $-\mathrm{\Delta }V\delta V_{\mathrm{H}}({𝐫}_{𝐩}).$

The derivative of $n({𝐫}_{𝐩})$ with respect to ${𝐑}_k$ is found by considering the definition of $n(𝐫)$ as follows:

${\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {𝐑}_k}}$ $=$ ${\displaystyle \sum _{i\alpha ,j\beta }}{\displaystyle \sum _{\nu }}\{{\displaystyle \frac{\partial c_{i\alpha ,\nu }^{*}}{\partial {𝐑}_k}}{c}_{j\beta ,\nu }{\chi }_{i\alpha }(𝐫){\chi }_{j\beta }(𝐫)+c_{i\alpha ,\nu }^{*}{\displaystyle \frac{\partial c_{j\beta ,\nu }}{\partial {𝐑}_k}}{\chi }_{i\alpha }({𝐫}_{𝐩}){\chi }_{j\beta }({𝐫}_{𝐩})\}$ (23)

$+2{\displaystyle \sum _{\alpha ,j\beta }}{\rho }_{k\alpha ,j\beta }{\displaystyle \frac{\partial {\chi }_{k\alpha }({𝐫}_{𝐩})}{\partial {𝐑}_k}}{\chi }_{j\beta }({𝐫}_{𝐩}).$

The derivative of the LCAO coefficient $c$ with respect to ${𝐑}_k$ can be transformed to the derivative of the overlap matrix by taking account of the orthonormality relation in one-particle wave functions as usually made in the LCAO method with the non-orthogonal basis set.[1] Since the overlap matrix element is a two-center integral, the matrix element and its derivatives are evaluated in the momentum space as well as the other two-center integrals. Also, the derivative of the basis function $\chi$ with respect to ${𝐑}_k$ is analytically evaluated using cubic splines because of the use of numerical basis function. The derivative of $n^{\mathrm{a}}({𝐫}_{𝐩})$ with respect to ${𝐑}_k$ is formed by only the contribution from the same atomic site $k$ due to the independent grid position from the atomic positions as follows:

${\displaystyle \frac{\partial n^{\mathrm{a}}({𝐫}_{𝐩})}{\partial {𝐑}_k}}$ $=$ ${\displaystyle \frac{\partial n_k^{\mathrm{a}}({𝐫}_{𝐩})}{\partial {𝐑}_k}}.$ (24)

The derivative of $E_{\mathrm{xc}}$ with respect to the atomic coordinate ${𝐑}_k$ is easily evaluated by

${\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial {𝐑}_k}}$ $=$ ${\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {𝐑}_k}}{\displaystyle \frac{\partial E_{\mathrm{xc}}}{\partial n({𝐫}_{𝐩})}},$ (25)

$=$ $\mathrm{\Delta }V{\displaystyle \sum _{𝐩}}{\displaystyle \frac{\partial n({𝐫}_{𝐩})}{\partial {𝐑}_k}}{v}_{\mathrm{xc}}(n({𝐫}_{𝐩})).$

Even for the GGA, the derivative of $E_{\mathrm{xc}}$ with respect to the atomic coordinate ${𝐑}_k$ is explicitly expressed, while the GGA is implemented by a finite difference scheme in our implementation. As a result, the derivative of the discretized $E_{\delta \mathrm{ee}}$ and $E_{\mathrm{xc}}$ are analytically evaluated. In addition, considering the localized basis functions $\chi$, the evaluation of Eqs. (19) and (25) can be quite efficiently performed in an O($N$) operation as well as the other contributions to the force.

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 $E_{\delta \mathrm{ee}}$ and $E_{\mathrm{xc}}$ 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.