Efficient O($N$) divide-conquer method with localized single-particle natural orbitals: Ver. 1.0
The original has been published in Phys. Rev. B 98, 245137 (2018).

Abstract

An efficient O($N$) divide-conquer (DC) method based on localized single-particle natural orbitals (LNOs) is presented for large-scale density functional theories (DFT) calculations of gapped and metallic systems. The LNOs are non-iteratively calculated by a low-rank approximation via a local eigendecomposition of a projection operator for the occupied space. Introducing LNOs to represent the long range region of a truncated cluster reduces the computational cost of the DC method while keeping computational accuracy. A series of benchmark calculations and high parallel efficiency in a multilevel parallelization clearly demonstrate that the O($N$) method enables us to perform large-scale simulations for a wide variety of materials including metals with sufficient accuracy in accordance with development of massively parallel computers.

1 INTRODUCTION

First-principles electronic structure calculations based on the density functional theories (DFT) [1, 2] have been playing a versatile role in a wide variety of materials sciences to deeply understand physical and chemical properties of existing materials and even to design novel materials having a desired property before actual experiments [3, 4, 5, 6]. In recent years more complicated materials with secondary structures such as heterointerfaces [7, 8] and dislocations [9] have been becoming a scope of application by DFT calculations. Since these complicated structures cannot be easily modeled by a small unit cell, development of efficient DFT methods in accordance with development of massively parallel computers is crucial to realize such large-scale DFT calculations. Among efficient DFT methods [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 26, 25, 27, 28, 29, 30, 31, 32, 33, 34], O($N$) methods whose computational cost scales linearly as a function of the number of atoms have enabled us to extend the applicability of DFT to large-scale systems [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Nevertheless applications of the O($N$) methods to metallic systems have been still limited because of the fundamental difficulty of a truncation scheme in real space, which is an idea commonly adopted in most of the O($N$) methods [12, 13], in realizing the O($N$) methods [14] as discussed in the Appendix. A theoretically proper approach to go beyond the truncation scheme is to take account of the contribution from the external region beyond the truncated region via a self energy in Green function formalisms [26, 27, 28, 29, 30]. Another straightforward approach is to use a relatively large cutoff radius in the truncation scheme in order to reach a sufficient accuracy. The latter approach might be suited to the divide-conquer (DC) method [31, 32] among the O($N$) methods proposed so far in the following twofold aspects: (i) There is a way to reduce the computational cost in the framework of the DC method by lowering the dimension of matrices with an introduction of a Krylov subspace [33, 34]. (ii) The computational time of the DC method can be reduced by a massive parallelization, since the calculations in the DC method are performed nearly independently for each atom [35, 36]. With the two aspects the improvement of the DC method can be a promising direction to develop an accurate, efficient, and robust O($N$) method applicable to not only insulators and semi-conductors, but also metals. Along this line, the DC method based on the Krylov subspace has been proven to be an efficient and accurate O($N$) method by a wide variety of applications such as dynamics of Li ions in a lithium ion battery [37, 38, 39] and structure optimization of semicoherent heterointerfaces in steel [7, 8]. However, there exist drawbacks in the generation of the Krylov subspace [34]. Since the Krylov subspace is generated at the first self-consistent field (SCF) step and kept unchanged during the subsequent SCF calculation, calculated quantities such as the electron density and total energy depend on the initial guess of electron density or Hamiltonian matrix elements. If the Krylov subspace is regenerated every SCF step to avoid the dependency on the initial guess, the computational efficiency must be largely degraded. In addition the iterative calculations in the generation of the Krylov subspace tend to suffer from numerical round-off error, leading to an uncontrollable behavior in the SCF calculation if the Krylov subspace is regenerated every SCF step [34]. Therefore, a more robust approach needs to be developed to improve the DC method, which overcomes the drawbacks inherent in the DC method based on the Krylov subspace.

In this paper, we focus on localized single-particle natural orbitals (LNOs) calculated by a low-rank approximation to perform a coarse graining of basis functions, and apply the LNOs to the DC method, which will be referred to as the DC-LNO method hereafter, to reduce the dimension of matrices without sacrificing the accuracy. The DC-LNO method overcomes the drawbacks in the DC method based on the Krylov subspace, while taking account of reduction of the prefactor in the computational cost and a simple algorithm with less communication leading to a high parallel efficiency. A series of benchmark calculations clearly demonstrates that the DC-LNO method is an accurate, efficient, and robust O($N$) method applicable to not only insulators and semi-conductors, but also metals in step with recent development of massively parallel computers.

The paper is organized as follows: In the section II we propose a method to generate LNOs, and present the DC-LNO method as an extension of the DC method. In the section III the implementations of the method are discussed in detail. In the section IV a series of benchmark calculations is presented. In the section V we summarize the theory of the DC-LNO method and numerical aspects.

2 THEORY

2.1 General

We consider an extension of a divide-conquer (DC) approach [31, 32] by introducing a coarse graining of basis functions as shown in Fig. 1. For each atom the Kohn-Sham (KS) Hamiltonian and overlap matrices are truncated within a given cutoff radius, and the resultant truncated cluster problem is solved atom by atom, leading to the O($N$) scaling in the computational cost. As expected the error of truncation can be systematically reduced in exchange for the increase of computational cost as the cutoff radius increases. The number of atoms in a truncated cluster exceeds 300 atoms in many cases in order to attain a sufficient accuracy as discussed later on. To reduce the computational cost we introduce a coarse graining of basis functions that the original basis functions, pseudo-atomic orbitals (PAOs) in our case [40, 41], are replaced by localized single-particle natural orbitals (LNOs) in the long range (yellow) region to represent the Hamiltonian and overlap matrices, while the PAO functions remain unchanged in the short range (orange) region. In the following subsections a method of generating LNOs and a DC method with LNOs are discussed in detail.

Figure 1: Truncation of a system in the DC method with LNOs. The short range (orange) and long range (yellow) regions are represented by PAOs and LNOs, respectively.

2.2 Generation of LNOs

We present a method of generating LNOs based on a low-rank approximation via a local eigendecomposition of a projection operator. The method might be applicable to any local basis functions, and not limited to the application to the DC method we discuss in the paper. Even in the conventional O($N^3$) calculations, the LNOs can be easily obtained by a non-iterative calculation using the density matrix and overlap matrix elements. Since a smaller number of LNOs well reproduce dispersion of occupied bands, they can be an alternative basis set for a compact representation for the Hamiltonian and overlap matrices. Though in this subsection we present a method of calculating LNOs in a general form starting from Bloch functions to clarify a mathematical basis of LNOs, it should be clearly noticed in this place that the density matrix is directly calculated in the DC-LNO method without calculating the Bloch functions. The general formulation presented here provides a theoretical basis of three step algorithm which will be discussed in the end of the subsection.

Under the Born-von Karman boundary condition we expand a Kohn-Sham (KS) orbital ${\varphi }_{𝐤\mu }$, indexed with a k-vector $𝐤$ in the first Brillouin zone and band index $\mu$, using PAOs $\chi$ [33, 34], being a real function, as

$|{\varphi }_{𝐤\mu }⟩={\displaystyle \frac{1}{\sqrt{N_{\mathrm{BC}}}}}{\displaystyle \sum _{𝐑}}{\mathrm{e}}^{\mathrm{i}𝐤\cdot 𝐑}{\displaystyle \sum _{i\alpha }}{c}_{𝐤\mu ,i\alpha }|{\chi }_{𝐑i\alpha }⟩,$ (1)

where $𝐑$, $N_{\mathrm{BC}}$, and $c$ are a lattice vector, the number of cells in the boundary condition, and linear combination of pseudo-atomic orbital (LCPAO) coefficients, respectively. It is also noted that $⟨𝐫|{\varphi }_{\mu }^{(𝐤)}⟩\equiv {\varphi }_{\mu }^{(𝐤)}(𝐫)$ and $⟨𝐫|{\chi }_{𝐑i\alpha }⟩\equiv {\chi }_{i\alpha }(𝐫-{\tau }_i-𝐑)$, where $i$ and $\alpha$ are atomic and orbital indices, respectively, and ${\tau }_i$ is the position of atom $i$. We assume that $M_i$ PAO functions are allocated to atom $i$. Throughout the paper we do not consider the spin dependency on the formulation for sake of simplicity, but the generalization is straightforward. By considering overlap matrix elements $S_{𝐑i\alpha ,{𝐑}^{\prime }j\beta }\equiv ⟨{\chi }_{𝐑i\alpha }|{\chi }_{{𝐑}^{\prime }j\beta }⟩$, one can introduce two alternative localized orbitals:

$|{\underset{¯}{\chi }}_{𝐑i\alpha }⟩$ $=$ ${\displaystyle \sum _{{𝐑}^{\prime }j\beta }}|{\chi }_{{𝐑}^{\prime }j\beta }⟩S_{{𝐑}^{\prime }j\beta ,𝐑i\alpha }^{-1/2},$ (2)

$|{\tilde{\chi }}_{𝐑i\alpha }⟩$ $=$ ${\displaystyle \sum _{{𝐑}^{\prime }j\beta }}|{\chi }_{{𝐑}^{\prime }j\beta }⟩S_{{𝐑}^{\prime }j\beta ,𝐑i\alpha }^{-1},$ (3)

where $S^{-1/2}$ and $S^{-1}$ are calculated from an overlap matrix $S$ for a supercell consisting of $N_{\mathrm{BC}}$ primitive cells in the Born-von Karman boundary condition. Hereafter, $\underset{¯}{\chi }$ and $\tilde{\chi }$ will be referred to as Löwdin [42] and dual orbitals [33], respectively. It is noted that we have the following relations:

$⟨{\underset{¯}{\chi }}_{𝐑i\alpha }|{\underset{¯}{\chi }}_{{𝐑}^{\prime }j\beta }⟩$ $=$ ${\delta }_{{\mathrm{𝐑𝐑}}^{\prime }}{\delta }_{ij}{\delta }_{\alpha \beta },$ (4)

$⟨{\chi }_{𝐑i\alpha }|{\tilde{\chi }}_{{𝐑}^{\prime }j\beta }⟩$ $=$ $⟨{\tilde{\chi }}_{𝐑i\alpha }|{\chi }_{{𝐑}^{\prime }j\beta }⟩={\delta }_{{\mathrm{𝐑𝐑}}^{\prime }}{\delta }_{ij}{\delta }_{\alpha \beta },$ (5)

and that the identity operator $\widehat{I}$ can be expressed in Bloch functions $\varphi$, Löwdin orbitals $\underset{¯}{\chi }$, or PAOs $\chi$ and dual orbitals $\tilde{\chi }$ as

$\widehat{I}$ $=$ ${\displaystyle \sum _{𝐤\mu }}|{\varphi }_{𝐤\mu }⟩⟨{\varphi }_{𝐤\mu }|,$ (6)

$=$ ${\displaystyle \sum _{𝐑i\alpha }}|{\underset{¯}{\chi }}_{𝐑i\alpha }⟩⟨{\underset{¯}{\chi }}_{𝐑i\alpha }|,$

$=$ ${\displaystyle \sum _{𝐑i\alpha }}|{\chi }_{𝐑i\alpha }⟩⟨{\tilde{\chi }}_{𝐑i\alpha }|={\displaystyle \sum _{𝐑i\alpha }}|{\tilde{\chi }}_{𝐑i\alpha }⟩⟨{\chi }_{𝐑i\alpha }|.$

We now define a projection operator for the occupied space by

$\widehat{P}={\displaystyle \sum _{𝐤\mu }}|{\varphi }_{𝐤\mu }⟩f({\epsilon }_{𝐤\mu })⟨{\varphi }_{𝐤\mu }|,$ (7)

where $f$ and $\epsilon$ are the Fermi-Dirac function and an eigenvalue of the KS equation, respectively. By using the second line of Eq. (6), one obtains an alternative expression of the projection operator as

$\widehat{P}={\displaystyle \sum _{𝐑i\alpha ,{𝐑}^{\prime }j\beta }}|{\underset{¯}{\chi }}_{𝐑i\alpha }⟩{\underset{¯}{\rho }}_{𝐑i\alpha ,{𝐑}^{\prime }j\beta }⟨{\underset{¯}{\chi }}_{{𝐑}^{\prime }j\beta }|$ (8)

with

${\underset{¯}{\rho }}_{𝐑i\alpha ,{𝐑}^{\prime }j\beta }$ $=$ ${\displaystyle \sum _{𝐤\mu }}⟨{\underset{¯}{\chi }}_{𝐑i\alpha }|{\varphi }_{𝐤\mu }⟩f({\epsilon }_{𝐤\mu })⟨{\varphi }_{𝐤\mu }|{\underset{¯}{\chi }}_{{𝐑}^{\prime }j\beta }⟩,$

$=$ ${\displaystyle \frac{1}{N_{\mathrm{BC}}}}{\displaystyle \sum _{𝐤\mu }}{\mathrm{e}}^{\mathrm{i}𝐤\cdot \left(𝐑-{𝐑}^{\prime }\right)}f({\epsilon }_{𝐤\mu })b_{𝐤\mu ,i\alpha }{b}_{𝐤\mu ,j\beta }^{*},$

where ${𝐛}_{𝐤\mu }=S^{1/2}{𝐜}_{𝐤\mu }$, and ${𝐜}_{𝐤\mu }$ is a column vector whose elements are LCPAO coefficients $\{c_{𝐤\mu ,i\alpha }\}$. As well, one can derive another expression of the projection operator by applying the third line of Eq. (6) to Eq. (7) as

$\widehat{P}={\displaystyle \sum _{𝐑i\alpha ,{𝐑}^{\prime }j\beta }}|{\chi }_{𝐑i\alpha }⟩{\rho }_{𝐑i\alpha ,{𝐑}^{\prime }j\beta }⟨{\chi }_{{𝐑}^{\prime }j\beta }|$ (10)

with

${\rho }_{𝐑i\alpha ,{𝐑}^{\prime }j\beta }$ $=$ ${\displaystyle \sum _{𝐤\mu }}⟨{\tilde{\chi }}_{𝐑i\alpha }|{\varphi }_{𝐤\mu }⟩f({\epsilon }_{𝐤\mu })⟨{\varphi }_{𝐤\mu }|{\tilde{\chi }}_{{𝐑}^{\prime }j\beta }⟩,$

$=$ ${\displaystyle \frac{1}{N_{\mathrm{BC}}}}{\displaystyle \sum _{𝐤\mu }}{\mathrm{e}}^{\mathrm{i}𝐤\cdot \left(𝐑-{𝐑}^{\prime }\right)}f({\epsilon }_{𝐤\mu })c_{𝐤\mu ,i\alpha }{c}_{𝐤\mu ,j\beta }^{*}.$

It is noted that $\underset{¯}{\rho }$ is related to $\rho$ by the following relation:

$\underset{¯}{\rho }=S^{1/2}\rho S^{1/2}.$ (12)

Remembering that the number of electrons $N_{\mathrm{ele}}$ in the supercell consisting of $N_{\mathrm{BC}}$ primitive cells can be obtained by the trace of $\widehat{P}$, and that we have two alternative expressions Eqs. (8) and (10) for $\widehat{P}$, one has the following expressions for $N_{\mathrm{ele}}$:

$N_{\mathrm{ele}}$ $=$ $2\mathrm{t}\mathrm{r}\left[\widehat{P}\right],$ (13)

$=$ $2{\displaystyle \sum _{𝐑i\alpha }}⟨{\underset{¯}{\chi }}_{𝐑i\alpha }|\widehat{P}|{\underset{¯}{\chi }}_{𝐑i\alpha }⟩=2\mathrm{t}\mathrm{r}\left[S^{1/2}\rho S^{1/2}\right],$

$=$ $2\mathrm{t}\mathrm{r}\left[\rho S\right]=2{\displaystyle \sum _{𝐑i\alpha }}⟨{\tilde{\chi }}_{𝐑i\alpha }|\widehat{P}|{\chi }_{𝐑i\alpha }⟩,$

where the factor of 2 is due to spin degeneracy. Since each term in the summation over $𝐑$ equally contributes to $N_{\mathrm{ele}}$, we have $N_{\mathrm{ele}}=N_{\mathrm{BC}}{N}_{\mathrm{ele}}^{(\mathrm{𝟎})}$, where $N_{\mathrm{ele}}^{(\mathrm{𝟎})}=2{\sum }_{i\alpha }⟨{\tilde{\chi }}_{\mathrm{𝟎}i\alpha }|\widehat{P}|{\chi }_{\mathrm{𝟎}i\alpha }⟩$. Thus, it is enough to consider $N_{\mathrm{ele}}^{(\mathrm{𝟎})}$ instead of $N_{\mathrm{ele}}$ for further discussion. By introducing a notation for a subset of orbitals $\{\chi \}$ and $\{\tilde{\chi }\}$ with

$|{\chi }_{𝐑i})$ $=$ $(|{\chi }_{𝐑i1}⟩,|{\chi }_{𝐑i2}⟩,\mathrm{\cdots },|{\chi }_{𝐑i{M}_i}⟩),$ (14)

$|{\tilde{\chi }}_{𝐑i})$ $=$ $(|{\tilde{\chi }}_{𝐑i1}⟩,|{\tilde{\chi }}_{𝐑i2}⟩,\mathrm{\cdots },|{\tilde{\chi }}_{𝐑i{M}_i}⟩),$ (15)

one can write $N_{\mathrm{ele}}^{(\mathrm{𝟎})}$ as

$N_{\mathrm{ele}}^{(\mathrm{𝟎})}=2{\displaystyle \sum _i}{\mathrm{tr}}_{\mathrm{𝟎}i}\left[({\tilde{\chi }}_{\mathrm{𝟎}i}|\widehat{P}|{\chi }_{\mathrm{𝟎}i})\right]=2{\displaystyle \sum _i}{\mathrm{tr}}_{\mathrm{𝟎}i}\left[{\mathrm{\Lambda }}_{\mathrm{𝟎}i}\right],$ (16)

where ${\mathrm{tr}}_{\mathrm{𝟎}i}$ means a partial trace over orbitals associated with an atom $i$ in the central cell with $𝐑=\mathrm{𝟎}$, and ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ is defined by

${\mathrm{\Lambda }}_{\mathrm{𝟎}i}={\displaystyle \sum _{𝐑j}}{\rho }_{\mathrm{𝟎}i,𝐑j}{S}_{𝐑j,\mathrm{𝟎}i}$ (17)

with definition of block elements:

${\rho }_{𝐑i,{𝐑}^{\prime }j}$ $=$ $({\tilde{\chi }}_{𝐑i}|\widehat{P}|{\tilde{\chi }}_{{𝐑}^{\prime }j}),$ (18)

$S_{𝐑i,{𝐑}^{\prime }j}$ $=$ $({\chi }_{𝐑i}|{\chi }_{{𝐑}^{\prime }j}).$ (19)

These block elements ${\rho }_{𝐑i,{𝐑}^{\prime }j}$ and $S_{𝐑i,{𝐑}^{\prime }j}$ are $M_i\times M_j$ matrices, where $M_i$ and $M_j$ are the number of PAO functions allocated to atoms $i$ and $j$, respectively. Therefore, ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ defined by Eq. (17) is a $M_i\times M_i$ matrix. It should be emphasized that Eq. (16) giving $N_{\mathrm{ele}}^{(\mathrm{𝟎})}$ by the sum of the partial trace is an important relation in calculating LNOs, since it shows that a local similarity transformation on an atomic site $i$ does not change $N_{\mathrm{ele}}^{(\mathrm{𝟎})}$ because of a property of the trace. Noting that ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ is nonsymmetric, we consider a general eigendecomposition of a nonsymmetric matrix for ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ among similarity transformations as

$V_{\mathrm{𝟎}i}^{-1}{\mathrm{\Lambda }}_{\mathrm{𝟎}i}{V}_{\mathrm{𝟎}i}={\lambda }_{\mathrm{𝟎}i},$ (20)

where ${\lambda }_{\mathrm{𝟎}i}$ is a diagonal matrix having eigenvalues $\{{\lambda }_{\mathrm{𝟎}i\gamma }\}$ of ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ as diagonal elements. Since an eigenvalue ${\lambda }_{\mathrm{𝟎}i\gamma }$ of ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ gives the population for the corresponding eigenstate of ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$, one can distinguish LNOs spanning the occupied space from others among all the eigenstates of ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ by monitoring the eigenvalues. To see the idea more clearly, we define an operator by

${\widehat{\mathrm{\Lambda }}}_{\mathrm{𝟎}i}={\displaystyle \sum _{\gamma }}|v_{\mathrm{𝟎}i\gamma }⟩{\lambda }_{\mathrm{𝟎}i\gamma }⟨{\tilde{v}}_{\mathrm{𝟎}i\gamma }|,$ (21)

where $|v_{\mathrm{𝟎}i\gamma }⟩$ is the $\gamma$-th column vector of $V_{\mathrm{𝟎}i}$, and $⟨{\tilde{v}}_{\mathrm{𝟎}i\gamma }|$ is the $\gamma$-th row vector of $V_{\mathrm{𝟎}i}^{-1}$. Note that $⟨{\tilde{v}}_{\mathrm{𝟎}i\gamma }|$ is the dual orbital of $|v_{\mathrm{𝟎}i\gamma }⟩$, and $⟨{\tilde{v}}_{\mathrm{𝟎}i\gamma }|v_{\mathrm{𝟎}i\eta }⟩={\delta }_{\gamma \eta }$. It is easy to confirm that ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ and ${\lambda }_{\mathrm{𝟎}i}$ in a matrix form can be obtained by representing ${\widehat{\mathrm{\Lambda }}}_{\mathrm{𝟎}i}$ with $\{{\chi }_{\mathrm{𝟎}i\alpha }\}$ and $\{{\tilde{\chi }}_{\mathrm{𝟎}i\alpha }\}$, and with $\{v_{\mathrm{𝟎}i\gamma }\}$ and $\{{\tilde{v}}_{\mathrm{𝟎}i\gamma }\}$, respectively. Thus, we have ${\mathrm{tr}}_{\mathrm{𝟎}i}[{\widehat{\mathrm{\Lambda }}}_{\mathrm{𝟎}i}]={\mathrm{tr}}_{i\mathrm{𝟎}}[{\mathrm{\Lambda }}_{\mathrm{𝟎}i}]$. If the eigenvalue ${\lambda }_{\mathrm{𝟎}i\gamma }$ is nearly zero, the contribution of the corresponding eigenstate $|v_{\mathrm{𝟎}i\gamma }⟩$ is negligible in the summation of Eq. (21). Therefore, ${\widehat{\mathrm{\Lambda }}}_{\mathrm{𝟎}i}$ can be approximated by excluding eigenstates whose eigenvalues are less than a threshold value ${\lambda }_{\mathrm{th}}$ that we will discuss later on. The treatment can be regarded as a low-rank approximation [43]. Then, using Eqs. (16) and (21) we can approximate the projection operator $\widehat{P}$ defined by Eq. (7) as

$\widehat{P}\simeq {\displaystyle \sum _{𝐑i}}{\displaystyle \sum _{\gamma }^{{\lambda }_{\mathrm{th}}\le {\lambda }_{𝐑i\gamma }}}|v_{𝐑i\gamma }⟩{\lambda }_{𝐑i\gamma }⟨{\tilde{v}}_{𝐑i\gamma }|,$ (22)

where terms satisfying a condition ${\lambda }_{\mathrm{th}}\le {\lambda }_{𝐑i\gamma }$ are taken into account in the summation over $\gamma$. If all the terms are included, Eq. (22) becomes equivalent to Eq. (7). We now define our LNOs by $\{v\}$ whose eigenvalues are larger than or equal to the threshold value ${\lambda }_{\mathrm{th}}$. By comparing Eq. (20) with Eq. (16), one has $V_{\mathrm{𝟎}i}^{-1}{\mathrm{\Lambda }}_{\mathrm{𝟎}i}{V}_{\mathrm{𝟎}i}=V_{\mathrm{𝟎}i}^{-1}({\tilde{\chi }}_{\mathrm{𝟎}i}|\widehat{P}|{\chi }_{\mathrm{𝟎}i})V_{\mathrm{𝟎}i}$. Thus, LNOs are defined by $\{v\}$, and it should be noted that $\{\tilde{v}\}$ are the corresponding dual orbitals. The approximate formula Eq. (22) for $\widehat{P}$ implies that a set of orbitals $\{v\}$, LNOs, well spans the occupied space, while the number of orbitals is reduced compared to the original PAOs. It is worth noting that LNOs can be independently calculated for each atom by a non-iterative calculation via Eqs. (17) and (20). The fact makes the method of generating LNOs very efficient, and also guarantees that the resultant LNOs associated with an atom $i$ are expressed by a linear combination of PAOs allocated to only the atom $i$ [44].

The computational procedure to generate LNOs is summarized as follows: (i) Calculation of $\mathrm{\Lambda }$. For each atom $i$ the matrix ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ is calculated by Eq. (17), where the summation over $𝐑$ is limited within a finite range because of the locality of PAOs in real space. (ii) Diagonalization of $\mathrm{\Lambda }$. Since the matrix ${\mathrm{\Lambda }}_{\mathrm{𝟎}i}$ is nonsymmetric, the diagonalization in Eq. (20) is performed by a generalized eigenvalue solver for a nonsymmetric matrix such as DGEEV in LAPACK [48]. (iii) Selection of $v$. Eigenvectors $\{v_{\mathrm{𝟎}i\gamma }\}$ whose eigenvalues are larger than or equal to the threshold value ${\lambda }_{\mathrm{th}}$ are selected as LNOs. Only the overlap and density matrices are required to calculate LNOs through the steps (i)-(iii) above. Therefore, either conventional O($N^3$) methods or O($N$) methods can be employed as eigenvalue solver as long as they generate the density matrix. As for LNOs other than those in the central cell with $𝐑=\mathrm{𝟎}$, it is apparent from the derivation that one can obtain $|v_{𝐑i\gamma }⟩$ by parallel translation of $|v_{\mathrm{𝟎}i\gamma }⟩$ with the lattice vector $𝐑$. The method can also be extended to choose another energy window. In the projection operator $\widehat{P}$ defined by Eq. (7), the Fermi-Dirac function is introduced to choose the occupied space. However, one can choose a proper energy window in the definition of the projection operator $\widehat{P}$ depending on what we discuss, which enables us to focus on specific bands such as localized $d$-bands near the Fermi level. In this sense, LNOs can be utilized like Wannier functions (WFs) [49, 50], while WFs are obtained through a unitary transformation of Bloch functions rather than the low-rank approximation, and they are orthonormal each other unlike LNOs. It is also worth pointing out that our method shares the basic idea based on the projection with other methods such as the quasiatomic orbitals scheme [51, 52, 53], a projection method [54], and a method via selected columns of the density matrix [55, 56].

It might be possible for the method we present in the paper to be applied for other localized basis functions such as finite element methods [57, 58] and finite difference methods [59, 60, 61]. In those cases one may introduce spatial partitioning methods such as Voronoi tessellation to decompose basis functions rather than focusing on basis functions on a single grid. A similar procedure can be applied for a set of partitioned basis functions.

2.3 DC method with LNOs

Here we consider an extension of the DC method [31, 32] using LNOs discussed in the previous subsection. Our theoretical basis to formulate the O($N$) DC method is that the total energy and atomic forces in the KS framework can be calculated by using electron density $n(𝐫)$, density matrix $\rho$, and energy density matrix $e$ defined by

$n(𝐫)$ $=$ ${\displaystyle \sum _i}\left(2{\displaystyle \sum _{\alpha ,𝐑j\beta }}{\rho }_{\mathrm{𝟎}i\alpha ,𝐑j\beta }{\chi }_{\mathrm{𝟎}i\alpha }(𝐫){\chi }_{𝐑j\beta }(𝐫)\right),$ (23)

$=$ ${\displaystyle \sum _i}{n}_i(𝐫),$

$\rho =-{\displaystyle \frac{1}{\pi }}\mathrm{Im}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}G(E+i{0}^{+})f(E)𝑑E,$ (24)

and

$e=-{\displaystyle \frac{1}{\pi }}\mathrm{Im}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}G(E+i{0}^{+})f(E)E𝑑E,$ (25)

where $G$ is the Green function defined by $G(Z)\equiv {(ZS-H)}^{-1}$ with the overlap matrix $S$ and KS matrix $H$, and the factor of 2 in Eq. (23) is due to spin degeneracy. It is remarked that forces on atoms in the DC method are not calculated variationally, but evaluated by using the formula derived by assuming that numerically exact KS wave functions are available as discussed in Ref. [34]. The DC method calculates the Green function $G(Z)$ approximately by introducing the truncation scheme as shown in Fig. 1. The KS matrix of the truncated cluster for atom $i$ are constructed using PAOs and LNOs as follows:

(26)

where the top left and bottom right blocks correspond to the short range (orange) region represented by PAOs, and the long range (yellow) region represented by LNOs, respectively, as shown in Fig. 1. The top right and the bottom left block consist of the hopping matrix elements bridging the two regions, and they are represented by both PAOs and LNOs. As well, the same structure is found for the overlap matrix. Noting that the computational bottleneck is mainly governed by the eigenvalue problem for the truncated clusters, and that the matrix size can be reduced by introducing LNOs compared to the conventional DC method, one can expect a considerable reduction of the computational cost as the size of the long range region increases. The idea of reducing the matrix dimension by introducing an effective representation of Hamiltonian is similar to that in the O($N$) Krylov subspace method [34] and the absolutely localized molecular orbitals (ALMO) method [23]. By solving the eigenvalue problem $H^{(i)}{c}_{\mu }^{(i)}={\epsilon }_{\mu }^{(i)}{S}^{(i)}{c}_{\mu }^{(i)}$ for the truncated cluster of atom $i$, we calculate matrix elements associated with the atom $i$ for the Green function as [62]

$G_{\mathrm{𝟎}i\alpha ,𝐑j\beta }^{(i)}(Z)={\displaystyle \sum _{\mu }}{\displaystyle \frac{c_{\mu ,\mathrm{𝟎}i\alpha }^{(i)}{(c_{\mu ,𝐑j\beta }^{(i)})}^{*}}{Z-{\epsilon }_{\mu }^{(i)}}}.$ (27)

Matrix elements ${\rho }_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$ and $e_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$ associated with the atom $i$ can be analytically calculated by inserting Eq. (27) into Eqs. (24) and (25), respectively. We only have to calculate ${\rho }_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$ and $e_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$ only if $S_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$ is non zero, since the other elements do not contribute to the total energy and forces on atoms in case of semi-local functionals such as local density approximations (LDA) [63, 64] and generalized gradient approximations (GGA) [65]. Then, $n_i(𝐫)$ in Eq. (23) is easily computed from ${\rho }_{\mathrm{𝟎}i\alpha ,𝐑j\beta }$. By applying the procedure for all the atoms in a system, all the necessary information to calculate the total energy and forces on atoms are obtained.

The overall procedure of the DC method with LNOs is summarized as follows: (i) Calculation of LNOs. LNOs are calculated by Eq. (20) for all atoms in the central cell with $𝐑=\mathrm{𝟎}$. At every SCF step LNOs are updated, leading to self-consistent determination of LNOs. (ii) Construction of $H^{(i)}$ and $S^{(i)}$. For each atom $i$ the KS and overlap matrices for a truncated cluster associated with the atom $i$ are constructed by Eq. (26). (iii) Diagonalization of $H^{(i)}{c}_{\mu }^{(i)}={\epsilon }_{\mu }^{(i)}{S}^{(i)}{c}_{\mu }^{(i)}$ by making use of a parallel eigenvalue solver. (iv) Finding a common chemical potential to conserve the total number of electrons in the system using Eq. (27), as discussed in Ref. [34] in detail. (v) Calculation of density matrix $\rho$ and energy density matrix $e$ using Eqs. (24), (25), and (27). (vi) Calculation of electron density $n$ using Eq. (23).

3 IMPLEMENTATIONS

We have implemented the DC-LNO method into the OpenMX DFT software package [66] which is based on norm-conserving pseudopotentials (PPs) [67, 68] and optimized pseudo-atomic orbitals (PAOs) [40, 41] as basis set. All the benchmark calculations were performed with a computational condition of a production level. The basis functions used are C6.0-s2p2d1, Si7.0-s2p2d1, Ti7.0-s2p2d1, O6.0-s2p2d1, Li8.0-s3p2, Al7.0-s2p2d1, and Fe5.5-s3p2d2 for carbon, silicon, titanium, oxygen, lithium, aluminum, and iron, respectively, where in the abbreviation of basis functions such as C6.0-s2p2d1, C stands for the atomic symbol, 6.0 the cutoff radius (Bohr) in the generation by the confinement scheme, and s2p2d1 means the employment of two, two, and one optimized radial functions for the $s$-, $p$-, and $d$-orbitals, respectively. The radial functions were optimized by a variational optimization method [40]. These basis functions we used can be regarded as double zeta plus polarization basis sets if we follow the terminology of Gaussian basis functions. As valence electrons in the PPs we included $2s$ and $2p$, $3s$ and $3p$, $3s$, $3p$, $3d$, and $4s$, $2s$ and $2p$, $1s$, $2s$, and $2p$, $3s$ and $3p$, and $3s$, $3p$, $3d$, and $4s$ states for carbon, silicon, titanium, oxygen, lithium, aluminum, and iron, respectively. All the PPs and PAOs we used in the study were taken from the database (2013) in the OpenMX website [66], which were benchmarked by the delta gauge method [69]. Real space grid techniques are used for the numerical integrations and the solution of the Poisson equation using FFT with the energy cutoff of 300 Ryd [70]. We used a generalized gradient approximation (GGA) proposed by Perdew, Burke, and Ernzerhof to the exchange-correlation functional [65]. An electronic temperature of 300 K is used to count the number of electrons by the Fermi-Dirac function for all the systems we considered.

The short and long range regions depicted in Fig. 1 are determined as follows: (i) We first pick up atoms in a sphere with a given cutoff radius $r_{\mathrm{L}}$. (ii) Among the atoms selected by the step (i) we distinguish the first neighboring atoms (FNAs) having non-zero overlap with the central atom in terms of basis functions, and remaining atoms other than FNAs are called the second neighboring atoms (SNAs), where the number of FNAs and SNAs are $N_{\mathrm{F}}$ and $N_{\mathrm{S}}$, respectively. (iii) The short range region is determined by adjusting a cutoff radius $r_{\mathrm{S}}$ so that the number of atoms in a sphere with a radius of $r_{\mathrm{S}}$ can be as close as possible to $N_{\mathrm{F}}+\kappa N_{\mathrm{S}}$, where the parameter $\kappa$ can vary from 0 to 1, and we will discuss the choice of $\kappa$ later on. (iv) The long range region consists of remaining atoms other than atoms selected by the step (iii). If we assign FNAs to atoms in the short range region, the total energy does not converge to the numerically exact one calculated by the conventional diagonalization method even if the cutoff radius $r_{\mathrm{L}}$ increases systematically. This is because the error with the low-rank approximation by Eq. (22) keeps increasing as the cutoff radius increases. To avoid the situation, we add a buffer region consisting of about $\kappa N_{\mathrm{S}}$ atoms as described by the step (iii) above, which guarantees the convergence of the total energy and other quantities as a function of the cutoff radius $r_{\mathrm{L}}$. Throughout the study, we used $\kappa$ of $\frac{3}{10}$ for all the systems by taking the accuracy into account more than the efficiency, and did not adjust the parameter, while a smaller value, which well balances both the accuracy and efficiency, can be employed for some systems.

The way of parallelization for the DC-LNO method on parallel computers will be discussed together with its benchmark calculations later on.

4 NUMERICAL RESULTS

4.1 Band dispersions by LNOs

In order to investigate to what extent LNOs can span occupied spaces, we compare band dispersions of gapped and metallic systems calculated with PAOs and LNOs. Figure 2(a) and (b) show band dispersions of diamond and silicon calculated by a conventional O($N^3$) diagonalization method with PAOs and LNOs. For both the cases the SCF calculations were performed by using PAOs. For the case of LNOs the band dispersions were calculated with the LNOs after the SCF calculations with PAOs. The number of LNOs per atom is 4 for both carbon and silicon atoms.

Figure 2: Band dispersions of (a) diamond and (b) silicon in the diamond structure with experimental lattice constants (3.567 and 5.430 Å) calculated by PAOs (black) and LNOs (red). A conventional O($N^3$) method was used for the diagonalization, where the number of k-points for the Brillouin zone sampling is $71\times 71\times 71$ for both the cases. In the case of LNOs, the SCF calculations were performed by using PAOs, and after determining the SC electron density, LNOs were used to calculate the band dispersion. The numbers of PAOs and LNOs per atom are shown in the parenthesis.

It is found that in both the cases the band dispersions of occupied space are well reproduced with LNOs compared to those calculated by PAOs, while a large difference can be seen in conduction bands between PAOs and LNOs as expected. The good agreement between PAOs and LNOs in describing the occupied bands implies that the low-rank approximation by Eq. (22) is reasonably valid. As shown in Table 1 we see that the first four eigenvalues of the matrix $\mathrm{\Lambda }$ are actually dominant for both diamond and silicon, justifying the low-rank approximation. The largest and the next three eigenvalues correspond to a $s$-orbital and $p$-like orbitals deformed by contribution of $d$-orbitals, respectively. It is also noted that the eigenvalues can be negative, which is related to a negative value of Mulliken populations for delocalized orbitals [71, 72]. As well as the gapped systems, similar calculations were performed for metals, lithium in the body centered cubic (BCC) structure and aluminum in the face centered cubic (FCC) structure as shown in Figs. 3 and 4, respectively. The band dispersions calculated with the minimal LNOs are reasonably compared to those by PAOs, while the use of the five and nine LNOs for Li and Al atoms fully reproduce the band dispersions including conduction bands as shown in Figs. 3(b) and 4(b), respectively. One can confirm again in Table 1 that eigenvalues for the minimal LNOs are dominant even for metals, while the magnitude of the subsequent eigenvalues is relatively large compared to those of the gapped systems. Thus, we conclude that LNOs can be regarded as a compact basis set spanning well the occupied space for both gapped and metallic systems.

Figure 3: Band dispersions of BCC lithium with an experimental lattice constant of 3.491 Å calculated by PAOs (black) and (a) two and (b) five LNOs (red). The number of k-points for the Brillouin zone sampling is $101\times 101\times 101$. The other details are the same as in the caption of Fig. 2.

Figure 4: Band dispersions of FCC aluminum with an experimental lattice constant of 4.050 Å calculated by PAOs (black) and (a) four and (b) nine LNOs (red). The number of k-points for the Brillouin zone sampling is $111\times 111\times 111$. The other details are the same as in the caption of Fig. 2.

Table 1: Eigenvalues $\lambda$ of the matrix $\mathrm{\Lambda }$ for diamond, silicon, BCC lithium, and FCC aluminum. The corresponding eigenvectors were used as LNOs to calculate the band dispersions shown in Figs. 2-4.

	Diamond	Si	Li	Al
${\lambda }_1$	0.514	0.639	0.999	0.551
${\lambda }_2$	0.483	0.430	0.277	0.253
${\lambda }_3$	0.483	0.430	0.075	0.253
${\lambda }_4$	0.483	0.430	0.075	0.253
${\lambda }_5$	0.012	0.018	0.074	0.049
${\lambda }_6$	0.012	0.018	0.001	0.049
${\lambda }_7$	0.012	0.018	0.001	0.049
${\lambda }_8$	0.011	0.012	0.001	0.033
${\lambda }_9$	0.011	0.012	-0.002	0.033
${\lambda }_{10}$	-0.001	0.002	-	-0.002
${\lambda }_{11}$	-0.001	0.002	-	-0.002
${\lambda }_{12}$	-0.001	0.002	-	-0.002
${\lambda }_{13}$	-0.018	-0.011	-	-0.015

4.2 Total energies by DC-LNO

As a first step of validation of the DC-LNO method with respect to computational accuracy and efficiency, we show in Fig. 5 the absolute error in the total energy of gapped and metallic systems calculated by the DC and the DC-LNO methods, where the reference energies were calculated by the conventional O($N^3$) method with dense k-points for the Brillouin zone sampling as given in the caption of Figs. 2-4. For all the cases the threshold value ${\lambda }_{\mathrm{th}}$ in Eq. (22) was set to be 0.1 which gives us the minimal LNOs corresponding to orbitals of valence electrons. In the gapped systems, diamond, silicon, and rutile TiO${}_2$, the absolute error decreases almost exponentially as the number of atoms in a truncated cluster increases. The overall behaviors of the error between the DC and DC-LNO methods are similar to each other, while the error by the DC-LNO method is accidentally much smaller than that by the DC method in the case of large truncated clusters of diamond. As well as the gapped systems the absolute error for metals calculated by the DC-LNO method decreases as increasing the number of atoms in a truncated cluster in a similar way to the DC method. The relatively large oscillating behavior observed in the metals might be related to long range characteristics of the off-diagonal Green functions as discussed in the Appendix. For all the cases including metals it is found that a truncated cluster including about 300 atoms is required to attain the millihartree accuracy corresponding to the error less than a few millihartree/atom in the total energy. From the comparison between the DC and DC-LNO methods, we see that the computational accuracy does not degrade largely even if basis functions for atoms in the long range region are approximated by LNOs, and that thereby the computational accuracy can be controlled mainly by the size of truncated cluster just like for the DC method. Since controlling only the single parameter allows us to balance the computational accuracy and efficiency, it is expected that the feature makes the DC-LNO method easy to use for a wide variety of applications.

Figure 5: Absolute error in the total energy (Hatree/atom) for (a) diamond, (b) silicon in the diamond structure, (c) rutile TiO${}_2$, (d) BCC lithium, (e) FCC aluminum, and (f) BCC iron as a function of the number of atoms in a truncated cluster calculated by the DC and DC-LNO methods. The experimental lattice constants were used for all the cases.

4.3 Computational time

Since the total number of basis functions to represent the Hamiltonian of the truncated cluster is reduced by introducing LNOs while keeping the accuracy, it is expected that the computational time can be substantially reduced. In the whole procedure of the DC-LNO method the calculation of LNOs and construction of Hamiltonian and overlap matrices occupy a small fraction of the whole computational time, typically less than $10\%$, and thereby the computational time is mainly governed by solving of the eigenvalue problem $H^{(i)}{c}_{\mu }^{(i)}={\epsilon }_{\mu }^{(i)}{S}^{(i)}{c}_{\mu }^{(i)}$ for each atom $i$. Noting that the computational time to solve the eigenvalue problem scales as the third power of the dimension of the matrices, the ratio of computational time between the DC-LNO and DC methods for elemental systems might be estimated by

${\displaystyle \frac{t_{\mathrm{DC}-\mathrm{LNO}}}{t_{\mathrm{DC}}}}={\displaystyle \frac{{\left[M_{\mathrm{P}}(N_{\mathrm{F}}+\kappa N_{\mathrm{S}})+M_{\mathrm{L}}(1-\kappa )N_{\mathrm{S}}\right]}^3}{M_{\mathrm{P}}^3{(N_{\mathrm{F}}+N_{\mathrm{S}})}^3}},$ (28)

where $M_{\mathrm{P}}$ and $M_{\mathrm{L}}$ are the number of PAOs and LNOs associated with each atom, respectively, and $\kappa$ is a factor, which is fixed to $\frac{3}{10}$ in this study, to control the size of the buffer region as discussed in the section ’IMPLEMENTATIONS’. For example, if the cutoff radius $r_{\mathrm{L}}$ is set to be $8.7$ Å in the diamond structure with the experimental lattice constant of 3.567 Å, $N_{\mathrm{F}}$ and $N_{\mathrm{S}}$ are found to be 167 and 298, and the resultant number of atoms in the short (long) range region becomes 275 (190). Then, $t_{\mathrm{DC}-\mathrm{LNO}}/t_{\mathrm{DC}}$ can be estimated to be 0.37 in the case that the number of PAOs and LNOs per atom is 13 and 4, which implies that the computational time of the DC-LNO method becomes about one third of that calculated by the DC method. Figure 6 shows actual timing results of the DC and DC-LNO methods for (a) diamond and (b) BCC lithium. We see that the actual $t_{\mathrm{DC}-\mathrm{LNO}}/t_{\mathrm{DC}}$ for the diamond case is 0.40 in the case that the number of atoms is 465 (=167+298) as shown in Fig. 6(a), which is well compared to the estimated value of 0.37. As indicated by Eq. (28) and in Figs. 6(a) and (b), it is concluded that the DC-LNO method becomes much faster than the DC method as the size of the truncated cluster increases.

Figure 6: Comparison of computational time of the diagonalization part per molecular dynamics (MD) step between the DC and DC-LNO methods for (a) diamond and (b) BCC lithium. Both the calculations were performed for the primitive cell using 14 MPI processes per atom on Intel Xeon CPUs E5-2690v4 @ 2.60 GHz.

In Fig. 7 we further show timing results of the DC-LNO method as a function of the number of atoms in the unit cell of Si crystal. It is confirmed from the linear increase of the computational time that the DC-LNO method is a linear scaling approach in practice. As a comparison the computational time is also shown for the conventional diagonalization method [73]. The crossing point between the two methods in the computational time is located around 100 atoms when 1280 CPU cores is used. Over 100 atoms the DC-LNO method is much faster than the conventional diagonalization method. The reason why the crossing point is located at such a small number of atoms is partly due to a better parallel efficiency of the DC-LNO method as discussed in the next subsection.

Figure 7: Computational time (sec.) of the diagonalization part per 10 SCF steps for Si in the diamond structure. In the DC-LNO method the cutoff radius $r_{\mathrm{L}}$ of $11.3$ Å was used, resulting in the truncated cluster containing 293 atoms. In the O($N^3$) method the k-grids of $2\times 2\times 2$ were used for the Brillouin zone sampling in all the cases. All the calculations were performed using 1280 MPI processes on a cluster machine consisting of Intel Xeon CPUs E5-2680v2 @ 2.80 GHz connected with Infiniband 4X FDR @ 56 Gbps.

4.4 Parallelization

To minimize the computational time on massively parallel computers we introduce a multilevel parallelization using message passing interface (MPI). In our implementation there are three levels for the parallelization, i.e., atom level, spin level, and diagonalization level as explained below: (i) Parallelization in the atom level. If the number of MPI processes is smaller than that of atoms, only the parallelization in the atom level is taken into account. The allocation of atoms to MPI processes is performed by a bisection method which is based on a projection of atoms onto a principal axis calculated from an inertia tensor and a modified binary tree of MPI processes to minimize memory usage and amount of MPI communications [36]. (ii) Parallelization in the spin level. If the number of MPI processes exceeds that of atoms, and the spin polarized calculation is performed, the parallelization in the spin level is introduced on top of the parallelization in the atom level, where a loop for the spin index is further parallelized. (iii) Parallelization in the diagonalization level. If the number of MPI processes is larger than the product of the number of atoms and the multiplicity of spin index, corresponding to 1, 2, and 1 for non-spin polarized, spin-polarized, and non-collinear calculations, respectively, a parallelization in the diagonalization level is further taken into account on top of both the parallelizations in the atom and spin levels. The parallelization in the diagonalization level is made by employing a parallel eigenvalue solver ELPA [73]. It is noted that the parallelization in the diagonalization level requires a considerable amount of MPI communications, while the parallelizations in the atom and spin levels have less MPI communications. So, one would expect a high parallel efficiency in the atom and spin levels, while the parallelization in the diagonalization level might be limited up to several tens of MPI processes. To achieve a better scaling for the parallelization in the diagonalization level, it is important to allocate CPU cores in the same computer node as MPI processes to avoid the inter-node communication as much as possible. We have implemented the multilevel parallelization so that amount of the inter-node communication can be minimized especially for the parallelization in the diagonalization level. In Fig. 8 the speed-up ratio in the MPI parallelization of the DC-LNO method is shown for non-spin polarized calculations of a diamond supercell containing 64 atoms. Since the multiplicity of spin index is 1, we see a nearly ideal behavior up to 64 MPI processes. Beyond 64 MPI processes the parallelization in the diagonalization level is taken into account on top of the parallelization in the atom level. A superlinear speed-up is observed at 128 and 256 MPI processes, which might be due to an effective use of cache by the reduction of memory usage, and a good scaling is achieved up to 1280 MPI processes at which the parallel efficiency is calculated to be 70% using the elapsed time at 1 MPI process as reference. Since each computer node has 20 CPU cores in this case, it would be reasonable to observe the good scaling up to 1280 (=64$\times$20) MPI processes. Thus, we see that the multilevel parallelization is very effective to minimize the computational time in accordance with a recent development of massively parallel computers.

Figure 8: Speed-up ratio in the MPI parallelization of the DC-LNO method for a diamond supercell containing 64 atoms, where the cutoff radius $r_{\mathrm{L}}$ of $8.0$ Å was used, leading to the numbers of atoms of 239 and 142 in the short and long range regions, and the dimension of matrices of 3675 for the truncated cluster problem. The calculations were performed using the same cluster machine as for Fig. 7.

4.5 Molecular dynamics simulations

To verify the accuracy of forces on atoms calculated by the DC-LNO method, results of NVE molecular dynamics (MD) simulations are shown in Fig. 9. We see that the sum of the kinetic energy of atomic nuclei (Kinetic) and the internal total energy ($E_{\mathrm{DFT}}$), being a conserved quantity, is reasonably conserved as a function of time, and the fluctuation is about one tenth of the kinetic energy or the internal total energy. It should be noted that the approximate conservation of the sum is achieved for not only Si being a semi-conductor, but also Al being a metal. Thus, it can be concluded that the accuracy of forces on atoms calculated by the DC-LNO method is sufficient for practical purposes, while it was remarked in the Sec. II that the forces on atoms in the DC-LNO method are not calculated variationally. The sufficient accuracy of the calculated forces is achieved by the use of large cutoff radii in constructing the truncation clusters, which is realized by both the introduction of LNOs and the massive parallelization with the multilevel parallelism.

Figure 9: The kinetic energy of atomic nuclei (Kinetic), the internal total energy ($E_{\mathrm{DFT}}$), and the sum of them as a function of time in NVE molecular dynamics simulations by the DC-LNO method for (a) Si and (b) Al, where the time step of 2 fs was used, and random velocities corresponding to 400 K were given for atomic nuclei at the first MD step. The simulation cells with experimental lattice constants for Si and Al contain 64 and 108 atoms, respectively. For the DC-LNO method the cutoff radii $r_{\mathrm{L}}$ of 11.3 and 10.1 Å were used for Si and Al, respectively, resulting in the truncated cluster containing 293 and 249 atoms in the ideal bulk structure. Each energy curve was shifted by adding a constant.

To further demonstrate the applicability of the DC-LNO method for MD simulations, we show radial distribution functions (RDFs) in liquid phases of silicon, aluminum, lithium, and SiO${}_2$ in Fig. 10. Since the electronic structures exhibit metallic features in the liquid phases of silicon, aluminum, and lithium, the MD simulations can be considered as a severe benchmark to validate the applicability of the DC-LNO method to metals. The cutoff radii $r_{\mathrm{L}}$ we used correspond to truncated clusters consisting of about 300 atoms in the ideal bulk structures. It turns out that in all the cases the DC-LNO method reproduces well the results by the conventional $O(N^3)$ diagonalization method, and that the obtained RDFs are well compared to other computational results [75, 76, 77, 78]. The considerable agreement between the DC-LNO and conventional methods strongly implies that a sufficient accuracy in reproducing at least RDF for MD simulations can be attainable with a cutoff radius $r_{\mathrm{L}}$ resulting in truncated clusters consisting of about 300 atoms for not only insulators but also metals. Thus, adjusting the cutoff radius $r_{\mathrm{L}}$ so that the number of atoms in a truncated cluster can be $\sim$ 300 atoms would be a compromise to balance the computational accuracy and efficiency, while the difference between the DC and DC-LNO methods in terms of the computational efficiency may not be significant for truncated clusters of this size. It is crucial to minimize the elapsed time for realization of long time MD simulations. With the computational condition the elapsed time per SCF step for silicon is 1.5 (sec.) on average using 1280 MPI processes on the same machine used for the calculations shown in Fig. 8.

Figure 10: Total radial distribution function (RDF) of (a) silicon at 3500 K, (b) aluminum at 2500 K, (c) lithium at 800 K, and (d) SiO${}_2$ at 3000 K, calculated by the conventional $O(N^3)$ diagonalization and the DC-LNO methods. The MD simulations were performed for cubic supercells containing 64, 108, 128 and 192 atoms with a fixed lattice constant of 10.86, 12.15, 14.04, 14.25 Å for silicon, aluminum, lithium, and SiO${}_2$, respectively, for 10 ps with the time step of 2 fs. The temperature was controlled by a velocity scaling scheme by Woodkock [74]. The coordinates for the first 1 ps were excluded to calculate RDF. In the DC-LNO method the cutoff radii $r_{\mathrm{L}}$ of 11.3, 10.1, 12.5, and 11.0 were used for silicon, aluminum, lithium, and SiO${}_2$, respectively, resulting in truncated clusters consisting of 293, 249, 339, and 344 (on average) atoms in the ideal bulk structures. In the conventional $O(N^3)$ diagonalization method, k-points of $7\times 7\times 7$, $8\times 8\times 8$, $7\times 7\times 7$, and $5\times 5\times 5$ were used for the Brillouin zone sampling in silicon, aluminum, lithium, and SiO${}_2$, respectively,

5 CONCLUSIONS

We have presented an efficient O($N$) method based on the DC approach and a coarse graining of basis functions by localized single-particle natural orbitals (LNOs) for large-scale DFT calculations. A straightforward way to attain sufficient accuracy in the DC method is to employ a relatively large cutoff radius for the truncation of a system, which is the most fundamental parameter in most of O($N$) methods to control the computational accuracy and efficiency. We have adopted the rather brute force approach, and attempted to decrease the computational cost by introducing LNOs as basis functions in the long range region of the truncated cluster, and to minimize the elapsed time in the computation with the help of a multilevel parallelization. The method of generating LNOs is based on a low-rank approximation to the projection operator for the occupied space by a local eigendecomposition at each atomic site, and the band structure calculations with PAOs and LNOs clearly show that the resultant LNOs span well the occupied space of not only gapped systems but also metals. It is also worth mentioning that the computational cost of generating LNOs is almost negligible thank to the independent calculation at each atomic site. By replacing PAOs with LNOs in the long range region of the truncated cluster in the DC method, the computational cost of the DC method can be reduced without largely sacrificing the accuracy. Nothing that the DC-LNO method holds the simple algorithm of the original DC method suited to the parallel calculation, we have implemented a multilevel parallelization using MPI by taking account of the atom level, spin level, and diagonalization level. It was demonstrated that the speed-up of the DC-LNO method by the multilevel parallelization can be expected up to a specific number of MPI processes which corresponds to the product of the number of atoms, the multiplicity of spin index, and the number of CPU cores in a single computer node. For example, if a spin-polarized calculation is performed for a system consisting of 1000 atoms on a parallel computer with 20 CPU cores per node, a high parallel efficiency might be expected up to 40000 MPI processes. As a validation of the applicability of the DC-LNO method, we have performed MD simulations for liquid phases of an insulator, semi-conductor, and metals, and confirmed that the RDFs calculated by the DC-LNOs are in good agreement with those by the conventional $O(N^3)$ diagonalization method, which may lead to its various applications to structural determinations of amorphous and liquid structures of complicated materials [79, 80, 81, 82, 83]. Considering the simplicity and robustness of the algorithm, we conclude that the DC-LNO method is an efficient and accurate approach to large-scale DFT calculations for a wide variety of materials including metals.

ACKNOWLEDGMENT

This work was supported by Priority Issue (Creation of new functional devices and high-performance materials to support next-generation industries) to be tackled by using Post ’K’ Computer, MEXT, Japan. Part of the computation in the study was performed using the computational facility of the Japan Advanced Institute of Science and Technology.

Appendix A: Asymptotic behaviors of the off-diagonal Green functions

As an example we show asymptotic behaviors of the off-diagonal Green functions for an one-dimensional (1D) tight-binding (TB) model with a single $s$-orbital on each site, and relate the asymptotic behaviors to electronic structures in gapped and metallic systems. The analysis interprets evidently the oscillating behavior of the error in the total energy of the metallic systems as shown in Fig. 5, and the rapid convergence in a high electronic temperature [33, 24].

Let us consider an orthogonal chain model with the nearest neighbor interaction $t$ and the on-site energy $\epsilon$ as defined by

$\widehat{H}=\epsilon {\displaystyle \sum _i}{\widehat{c}}_i^{†}{\widehat{c}}_i+t{\displaystyle \sum _i}({\widehat{c}}_i^{†}{\widehat{c}}_{i+1}+\mathrm{h}.\mathrm{c}.),$ (29)

where $t$ is assumed to be positive. By tri-diagonalizing the Hamiltonian with a Lanczos algorithm starting from a site ($i=0$), and calculating the diagonal Green function via a continued fraction using the recursion method [84, 85], one obtains a well known result for the diagonal Green function $G_{00}$ as follows:

$G_{00}(Z)={\displaystyle \frac{1}{\sqrt{{(Z-\epsilon )}^2-4t^2}}}.$ (30)

The off-diagonal Green functions can be obtained by using a recurrence relation [86] derived from $G^{(L)}(Z)(Z-H^{(L)})=I$, where $G^{(L)}$ and $H^{(L)}$ are the Green function and Hamiltonian matrices represented by the Lanczos vectors, and by performing a back unitary transformation as

$G_{01}(Z)$ $=$ $G_{00}(Z){\displaystyle \frac{\gamma }{2}}-{\displaystyle \frac{1}{2t}},$ (31)

$G_{02}(Z)$ $=$ $G_{00}(Z)\left({\displaystyle \frac{{\gamma }^2}{2}}-1\right)-{\displaystyle \frac{\gamma }{2t}},$ (32)

$G_{03}(Z)$ $=$ $G_{00}(Z)\left({\displaystyle \frac{{\gamma }^3}{2}}-{\displaystyle \frac{3}{2}}\gamma \right)-{\displaystyle \frac{{\gamma }^2-1}{2t}},$ (33)

$G_{04}(Z)$ $=$ $G_{00}(Z)\left({\displaystyle \frac{{\gamma }^4}{2}}-2{\gamma }^2+1\right)-{\displaystyle \frac{{\gamma }^3-2\gamma }{2t}},$ (34)

where $\gamma =(Z-\epsilon )/t$ and $G_{0j}$ is the off-diagonal element of Green function between the sites $0$ and $j$. It turns out that the off-diagonal Green functions can be expressed by $G_{00}$ and $\gamma$. To see asymptotic behaviors of the off-diagonal Green functions, by employing the following formula [87]:

${\displaystyle \frac{1}{\sqrt{a^2-x^2}}}={\displaystyle \frac{1}{a}}{\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{}_{2n}C_n{4}^n{a}^{2n}{x}^{2n}$ (35)

with the radius of convergence $|a|$, we Taylor exand $G_{00}$ at ${\gamma }^{-1}=0$ as

$G_{00}(Z)$ $=$ ${\displaystyle \frac{1}{t}}{\displaystyle \sum _{n=0}^{\mathrm{\infty }}}{\displaystyle \frac{{}_{2n}C_n}{4^n{2}^{-2n}}}{\gamma }^{-(2n+1)},$ (36)

$=$ ${\displaystyle \frac{1}{t}}\left({\displaystyle \frac{1}{\gamma }}+{\displaystyle \frac{2}{{\gamma }^3}}+{\displaystyle \frac{6}{{\gamma }^5}}+{\displaystyle \frac{20}{{\gamma }^7}}+\mathrm{\cdots }\right),$

where the convergence is guaranteed for $|\gamma |>2$. By inserting Eq. (36) into Eqs. (31)-(34), and taking the leading terms, we obtain the following relation:

$G_{0j}(Z)\propto {\displaystyle \frac{1}{t{\gamma }^{j+1}}}.$ (37)

Thus, we see that $G_{0j}$ approaches to zero asymptotically for $|\gamma |>2$ as $j\to \mathrm{\infty }$. On the other hand the Green functions at $Z=\epsilon$ corresponding to $\gamma =0$ are given by

$G_{0(2k-1)}(\epsilon )$ $=$ ${\displaystyle \frac{{(-1)}^k}{2t}},$ (38)

$G_{0(2k)}(\epsilon )$ $=$ ${(-1)}^k{G}_{00}(\epsilon ),$ (39)

where $G_{00}(\epsilon )=\frac{-i}{2t}$. It is found that $G_{0j}$ at $\gamma =0$ exhibits an oscillating behavior as a function of $j$, and never decays.

Figure A.1: Relation between the position of Matsubara poles (black filled circels), the spectrum range on the real axis (blue filled rectangles), and the convergent region of the off-diagonal Green functions whose boundaries are shown by red circles in the complex plane for (a) metal and (b) insulator. In the simple TB model the band width $w$ is given by $4t$, and the off-diagonal Green functions $G_{0j}$ decays asymptotically as $j$ increases in the exterior region of the red circles.

We now relate the asymptotic behaviors of Green functions to the calculation of density matrix which is defined by Eq. (24). Introducing the Matsubara expansion of the Fermi-Dirac function, and changing the integration path with the Cauchy theorem, one has [88]

${\rho }_{0j}={\displaystyle \frac{1}{2}}{\delta }_{0j}+\mathrm{Im}\left[{\displaystyle \frac{2i}{\beta }}{\displaystyle \sum _{p=1}^{\mathrm{\infty }}}{G}_{0j}({\alpha }_p)\right],$ (40)

where ${\alpha }_p$ are Matsubara poles located at $\mu +i\frac{(2p-1)\pi }{\beta }$ with a chemical potential of $\mu$ and $\beta =\frac{1}{k_{\mathrm{B}}T}$. The expression allows us to figure out a relation between the Matsubara poles, where the Green functions are evaluated, and the convergent region of the off-diagonal Green functions as illustrated in Fig. A.1. Remembering that in the simple TB model the band width $w$ is given by $4t$, and assuming that the single band is half filled in the metallic case, we may have Matsubara poles in the red circle which is the non-convergent region of the off-diagonal Green functions as shown in Fig. A.1(a). Since the off-diagonal Green functions evaluated at the Matsubara poles in the red circle do not simply decay in real space as $j$ increases, the truncation scheme commonly adopted in most of O($N$) methods should suffer from the long range characteristics of the off-diagonal Green functions, while the effect can appear in a different way depending on underlying principles of each O($N$) method. In the DC-LNO method the truncated eigenvalue problem $H^{(i)}{c}_{\mu }^{(i)}={\epsilon }_{\mu }^{(i)}{S}^{(i)}{c}_{\mu }^{(i)}$ is solved for each atom $i$, and the integrations of Eqs. (24) and (25) can be easily performed on the real axis since we have the approximate spectrum representation of Eq. (27). The way of evaluating the density matrix is numerically equivalent to the computational method via a generalized formula of Eq. (40) to the non-orthogonal basis set, where the Green functions for the truncated problem are computed at each Matsubara pole by the inverse calculation, since the Green function computed through the spectrum representation is exactly the same as one computed by the inverse calculation. Therefore, the oscillating behavior of error in the total energy calculation observed in Figs. 5(d), (e), and (f) should be attributed to the long range characteristics of the off-diagonal Green functions. It can also be understood that the use of a higher electronic temperature suppresses the deficiency since all the Matsubara poles can be placed in the exterior region of the red circles beyond a critical temperature [33, 24]. On the other hand we model an insulator by considering two bands as shown in Fig. 1(b), where each of them is expressed by the 1D TB model and the bands are separated by a finite gap. Unlike the metallic case, all the Matsubara poles are located in the exterior region of the red circles. The feature guarantees that ${\rho }_{0j}$ decays as $j$ increases since all the Green functions in the summation of Eq. (40) decay as $j$ increases, theoretically justifying that the truncation scheme is valid for gapped systems, although our benchmark calculations imply that the use of a large cutoff radius diminishes the effect of the long range characteristics of the off-diagonal Green functions even to metals at least for the calculations of density matrix and energy density matrix in a practical sense.