⇦Chinese Physics Letters, 2024, Vol. 41, No. 5, Article code 053102Express Letter Random Green's Function Method for Large-Scale Electronic Structure Calculation Mingfa Tang (汤明发)1, Chang Liu (刘畅)2, Aixia Zhang (张爱霞)1, Qingyun Zhang (张青云)1, Jiayu Zhai (翟佳羽)3, Shengjun Yuan (袁声军)4,5, and Youqi Ke (柯友启)1* Affiliations 1School of Physical Science and Technology, ShanghaiTech University, Shanghai 201210, China 2Xiaogan Sichuang Information Technology Co., LTD, Xiaogan 432000, China 3Institute of Mathematical Sciences, ShanghaiTech University, Shanghai 201210, China 4Key Laboratory of Artificial Micro- and Nano-structures of Ministry of Education, and School of Physics and Technology, Wuhan University, Wuhan 430072, China 5Wuhan Institute of Quantum Technology, Wuhan 430206, China Received 17 April 2024; accepted manuscript online 23 April 2024; published online 28 April 2024 *Corresponding author. Email: keyq@shanghaitech.edu.cn Citation Text: Tang M F, Liu C, Zhang A X et al. 2024 Chin. Phys. Lett. 41 053102    Abstract We report a linear-scaling random Green's function (rGF) method for large-scale electronic structure calculation. In this method, the rGF is defined on a set of random states and is efficiently calculated by projecting onto Krylov subspace. With the rGF method, the Fermi–Dirac operator can be obtained directly, avoiding the polynomial expansion to Fermi–Dirac function. To demonstrate the applicability, we implement the rGF method with the density-functional tight-binding method. It is shown that the Krylov subspace can maintain at small size for materials with different gaps at zero temperature, including H$_{2}$O and Si clusters. We find with a simple deflation technique that the rGF self-consistent calculation of H$_{2}$O clusters at $T=0$ K can reach an error of $\sim$ $1$ meV per H$_{2}$O molecule in total energy, compared to deterministic calculations. The rGF method provides an effective stochastic method for large-scale electronic structure simulation.
cpl-41-5-053102-fig1.png
cpl-41-5-053102-fig2.png
 DOI:10.1088/0256-307X/41/5/053102 © 2024 Chinese Physics Society Article Text An accurate and linear-scaling electronic structure (ES) method is essential for large-scale applications, and has being highly pursued in the materials simulation community.[1,2] The standard implementation of the Kohn–Sham density functional theory (DFT)[3,4] is difficult to deal with systems containing more than hundreds or thousands of atoms, due to the cubic-scaling computational complexity of diagonalization. To circumvent the diagonalization, many classic linear-scaling methods have been developed by employing various approximations based on Kohn's nearsightedness principle,[5] including divide-and-conquer[6] and other approximate density matrix (DM) techniques.[5,7-9] Nonetheless, as cost for achieving linear scaling, these O(N) ES methods suffer from the large prefactor and the limited accuracy, and thus present limitations on the dimensionality and the character of systems.[1,2] Recently, the random-state method[10] has presented a new avenue for developing linear-scaling (or even sub-linear) simulations of materials' ES, including the reports of stochastic DFT,[11] and time-dependent random-state propagation based DFT.[12] The random-state first-principles simulations of over $10^4$ atoms have been reported.[12] It has been proven that the random-state method bears only the unbiased distributed statistical error of the order $1/\sqrt{N_{\rm s}}$ ($N_{\rm s}$ is the number of random states), and the error in the averaged global property, e.g., energy per electron, scales as $1/\sqrt{N_{\rm e}N_{\rm s}}$ ($N_{\rm e}$ is the electron number) due to the self-average.[13,14] Various techniques, e.g., the deflation[15] and probing[16-20] approaches, exist for effectively reducing the statistical errors, providing the important basis for applying random-state method for ES simulations. It has been shown that the statistical error in random-state based DFT can be effectively reduced by about an order of magnitude by combining with the techniques of embedded fragment,[21-23] energy windows,[24] a universal mixed stochastic-deterministic approach,[25] to achieve important accuracy in charge density, energy, and forces. The present implementation of random-state based DFT utilizes the Chebyshev polynomial expansion to represent the Fermi–Dirac (FD) operator,[11-13] and it is thus suitable for high temperature simulation[25] and large band-gap systems (where the smearing of FD filter with high value of temperature can be applied to approximate the zero-temperature results). However, for narrow-gap systems (including the non-uniform mixtures of large and small gap materials), a large number of Chebyshev polynomials are required to converge the low-temperature results, because of the sharp step-like behavior of Fermi–Dirac function.[11,12] It is noted that a recent report of the rational approximation for implementing random-state methods presents the computational cost similar to the polynomial approximation for gold metal clusters at finite temperature.[19] Besides the applications to DFT, the random-state technique has also been applied to accelerate the important calculations of the GW (G refers to Green's function and W denotes the dynamically screened Coulomb interaction),[26] Bethe–Salpeter equation[27] and time-dependent DFT,[28] where the Chebyshev expansion is generally applied to the Heaviside step functions. In this Letter, we present a linear-scaling random Green's function (rGF) method for simulation of large-scale systems by combining with the Krylov subspace approach. For the application to DFT, we show the Fermi–Dirac operator can be directly obtained with the rGF method, avoiding the polynomial expansion to the Fermi–Dirac function. We demonstrate the rGF method by implementing it with density functional tight-binding (DFTB) approach. To start, a general rGF $\tilde{G}$ can be introduced as \begin{equation} \tilde{G}=G\cdot\tilde{I}, \tag {1} \end{equation} where $G$ is the exact GF which models the system's spatial and temporal response possessing wide applications in quantum physics. For example, Keldysh's GF $G^{\mathcal{K}}$ and the retarded GF $G^{\mathcal{R}}$ are the central quantities for simulating the respective nonequilibrium and equilibrium systems.[29,30] Here, $\tilde{I}$ is an approximate random resolution of the projection operator or identity matrix,[11,31] given by \begin{equation} \tilde{I}=\frac{1}{N_{\rm s}}\sum_{i=1}^{N_{\rm s}}| X_i \rangle \langle X_i |, \tag {2} \end{equation} where $|X_i \rangle$ is the random state with the random elements 1 or $-1$ of equal probability (other choice of $\tilde{I}$ is possible), $N_{\rm s}$ is the number of random states. It is clear that $\tilde{I}=I+\delta$ with $\delta_{ii}=0$. The random $\delta_{ij} (i\neq j)$ satisfies the Gaussian distribution (due to the central limit theorem) with the average $\langle \delta_{ij} \rangle=0$, and a standard deviation $\langle \delta^2_{ij} \rangle = 1/N_{\rm s}$. For an application to the equilibrium noninteracting case, we consider a physical system with the overlap $O$ (for generality) and Hamiltonian $H$ defined with $N$ localized bases. Then, the random retarded GF, namely $\tilde{G}^{\mathcal{R}}(z)=G^{\mathcal{R}}(z)\tilde{I}$, where $G^{\mathcal{R}}(z)=[zO-H]^{-1}$ ($z=\epsilon+i\eta$, $\eta$ is infinitesimal), can be written as \begin{equation} \tilde{G}^{\mathcal{R}}(z) = \frac{1}{N_{\rm s}}\sum_{i=1}^{N_{\rm s}}|R_i(z) \rangle \langle X_i|, \tag {3} \end{equation} where the vector $|R_i(z)\rangle$ satisfies the linear equation \begin{equation} [zO-H]|R_i(z)\rangle=|X_i\rangle, \tag {4} \end{equation} which features a shifted linear system for different $z$, and can be solved by Krylov subspace methods.[32-34] For applications for many $z$ points and the rGF integration, we can adopt a common Krylov subspace,[33,34] which is constructed with the Arnoldi process,[33] denoted by $K^{(i),\,v} = K_{v}(H',\,|\mathcal{X}_i\rangle)$, where $v$ is the size of subspace, $H'=O^{-1}H$ and $|\mathcal{X}_i\rangle=O^{-1}|X_i\rangle$. Here, $K^{(i),\,v}$ contains $\{|\varPhi^{(i),\,v}_m\rangle,\,m=1,\,\ldots ,\,v\}$ in which the involved $O^{-1}$-vector product is solved efficiently with the iterative conjugate gradient (CG) linear solver. Then, we can obtain a projected $v \times v$ subspace overlap $O^{(i),\,v}_{K,\,mn}=\langle \varPhi^{(i),\,v}_m|O|\varPhi^{(i),\,v}_n\rangle$ and Hamiltonian $H^{(i),\,v}_{K,\,mn}=\langle \varPhi^{(i),\,v}_m|H|\varPhi^{(i),\,v}_n\rangle$. Then, in $K^{(i),\,v}$, the subspace GF is obtained as $\tilde{G}^{(i),\,v}_K(z)=(zO^{(i),\,v}_K-H^{(i),\,v}_K)^{-1}$. We have \begin{equation} |R_i(z)\rangle = \sum^{v}_{m,n}G^{(i),v}_{K,mn}(z) |\varPhi^{(i),v}_m\rangle \langle\varPhi^{(i),v}_n|X_i\rangle, \tag {5} \end{equation} for giving rGF in Eq. (3). For a large scale system, $N_{\rm s}$ in Eq. (3) for a good accuracy can be much smaller than the dimension of $H$, providing an efficient linear-scaling approach for calculating GF of the system. For an application to electronic structure self-consistent calculation, the randomized FD operator, namely $\tilde{F}(O^{-1}H)$ and $F(X)=1/(\exp[\beta(X-\mu)]+1)$ ($\beta=1/k_{\scriptscriptstyle{\rm B}}T$, $\mu$ is chemical potential), can be connected to the rGF by \begin{equation} \tilde{F}=-\frac{1}{\pi}{\rm Im}\Big\{O\int^{+\infty}_{-\infty} F(\epsilon) \tilde{G}^{\mathcal{R}}(z)d\epsilon\Big\}=F\tilde{I}. \tag {6} \end{equation} By introducing the polynomial expansion to $F({\epsilon})$ in Eq. (6), one can obtain the polynomial approximation of $\tilde{F}(O^{-1}H)$, e.g., the Chebyshev expansion for the DM,[11,12] avoiding the complexity for calculating rGF and the integration. However, Polynomial expansion of Fermi–Dirac operator requires high orders for low-temperature properties of narrow-gap materials. In this work, by combining Eqs. (3) and (5), Eq. (6) can be projected onto a Krylov subspace integration, namely \begin{equation} \tilde{F}^{(i),v}_{K}= -\frac{1}{\pi}{\rm Im}\Big[\int dz F(z) \tilde{G}_K^{(i),v}(z)\Big], \tag {7} \end{equation} which can be obtained analytically with subspace diagonalization (see the Supplementary Material[35]). Finally, we can obtain \begin{equation} O^{-1}\tilde{F}=\frac{1}{N_{\rm s}}\sum^{N_{\rm s}}_{i}\sum^{v}_{m,n}\tilde{F}_{K,mn}^{(i),v}|\varPhi^{(i),v}_m\rangle \langle\varPhi^{(i),v}_n|X_i\rangle\langle X_i |. \tag {8} \end{equation} It should be noted that Eq. (8) essentially provides a Krylov expansion with $v$ polynomials of $H'$ that is free from the temperature dependence (distinct from the Chebyshev expansion), featuring a linear scaling computational cost. As an application of rGF, the generalized non-interacting DM can be obtained as $\tilde{\mathcal{D}}=O^{-1}\tilde{F}=\mathcal{D}\tilde{I}$ ($\mathcal{D}=O^{-1}F$ is the exact DM). The Krylov subspace method provides a controllable approach to calculate $\tilde{D}$ with high accuracy. The convergence of the Krylov subspace for $\tilde{D}$ can be related to the number of distinct clusters of eigenvalues of the material system and their width.[36] Moreover, for the linear-algebraic feature of the rGF method, a large volume of techniques in linear algebra, e.g., various preconditioning techniques, can be combined to further improve the efficiency. For the converged $v$, the DM $\mathcal{D}$ only bears a statistical error, namely $\mathcal{D}\cdot\delta$, which scales as $1/\sqrt{N_{\rm s}}$. For an effective error reduction, the general deflation and probing approach can be applied.[15-20] To demonstrate the applicability of the rGF method, we implement it with DFTB[37-41] to calculate the ES and total energy and compare with the deterministic method. We use the non-orthogonal basis which usually possesses the higher accuracy and transferability than orthogonal basis. In each DFTB self-consistent iteration, the Mulliken charge $q_{\scriptscriptstyle{R}}$ for site $R$, required for updating the Hamiltonian $H$, is calculated as \begin{equation} q_{\scriptscriptstyle{R}}= {\rm Tr}(\{\tilde{\mathcal{D}}\cdot O\}_{\scriptscriptstyle{RR}}). \tag {9} \end{equation} In DFTB, the total orbital energy is obtained as \begin{equation} E_{\rm orbit}={\rm Tr}(\tilde{ \mathcal{D}}\cdot H), \tag {10} \end{equation} and then the DFTB total energy is given by \begin{equation} E_{\rm tot}=E_{\rm orbit}+E_{\rm rep}-E_{\rm coul}, \tag {11} \end{equation} where $E_{\rm rep}$ is pairwise energy in which the parameters are fitted to the high-level theoretical calculation (e.g., the Kohn–Sham DFT), and $E_{\rm coul}$ is coulomb energy that corrects the double counting in $E_{\rm orbit}$.[37-41] In the rGF-DFTB calculations, to ensure the convergence of the self-consistent procedure, the same set of random states are used in all the iterations, avoiding the stochastic noise. For applications, all the codes are compiled with the intel Math Kernel Library, and we use the CG solver package PETSc[42-44] for $O^{-1}$-vector product and matrix-vector multiplication in constructing the Krylov subspace. The present rGF-DFTB solver is implemented based on the open-source package LATTE. All the calculations are performed on the Intel Xeon Gold 6226 CPU (2.70 GHz). We adopt exact diagonalization for the Krylov subspace integration in Eq. (7) for different $T$.[35]
cpl-41-5-053102-fig1.png
Fig. 1. Deviation of $E_{\rm orbit}$ versus $v$ for (H$_{2}$O)$_{526}$, Si$_{1024}$H$_{384}$ and $\overline{\rm Si}_{1024}$H$_{384}$ clusters at $T=0$. Inset: total CPU time per iteration for the full-DFTB and rGF ($N_{\rm s}=1000$ and the fixed $v=35$ for safely converging all larger systems) calculations performed on a single CPU core.
In Fig. 1, we investigate the convergence of Krylov subspace, namely deviation of $E_{\rm orbit}$ versus $v$ ($N_{\rm s}=1000$), for (H$_2$O)$_{526}$, Si$_{1024}$H$_{384}$, and $\overline{\rm Si}_{1024}$H$_{384}$ clusters at $k_{\scriptscriptstyle{\rm B}}T=0$ eV (based on the self-consistent charges of full DFTB in LATTE).[45] We use the minimal 6 and 4 atomic bases for the H$_2$O and Si, respectively, with the associated DFTB parameters as in Refs. [46-48] H$_2$O clusters are generated with the solvation box module in VMD,[49] and Si- and $\overline{\rm Si}$-clusters are cut from the diamond bulk with the lattice constants 3.86 Å and 4.06 Å, respectively. The HOMO-LUMO energy gaps are 12.63 eV, 2.0 eV, and 0.89 eV for the H$_2$O, Si, and $\overline{\rm Si}$ clusters, respectively. It is found that $E_{\rm orbit}$ can be converged with an error $\le$ 5 meV with $v=20$, 50, and 100 for the H$_2$O, Si, and $\overline{\rm Si}$ clusters, respectively. However, as shown in the Supplementary Material (where we compare the convergence of $E_{\rm orbit}$ versus Krylov and Chebyshev polynomials at different temperatures for FD operator),[35] to safely converge the energy within a few meV at $k_{\scriptscriptstyle{\rm B}}T=0.03$ eV, Chebyshev expansion requires 2000 polynomials for H$_2$O and 1400 polynomials for both Si$_{1024}$H$_{384}$ and $\overline{\rm Si}_{1024}$H$_{384}$, consistent with Ref. [11]. It should be noted that investigations of the convergence at different low temperatures for H$_2$O and Si clusters are meaningful for their mixture with small-gap materials, which determine the least bound for the polynomials, e.g., Chebyshev expansion. As shown in the Supplementary Material,[35] to safely produce the $E_{\rm orbit}$ at $T=0$ (by positioning $\mu$ at the center of gap), we find the maximal $k_{\scriptscriptstyle{\rm B}}T_{\max}\approx E_{\rm gap}/20$, for example, $k_{\scriptscriptstyle{\rm B}}T_{\max}$ values for H$_2$O, Si$_{1024}$H$_{384}$, and $\overline{\rm Si}_{1024}$H$_{384}$ are estimated to be about 0.5, 0.1, and 0.05 eV, respectively. It is found that increasing the temperature can significantly reduce the required number of Chebyshev polynomials, while the Krylov polynomials present slight decreases for FD operator of gapped systems (since the Krylov subspace method explores the intrinsic structure of the materials' ES). For the present (H$_2$O)$_{526}$, Si$_{1024}$H$_{384}$, and $\overline{\rm Si}_{1024}$H$_{384}$ clusters at their $k_{\scriptscriptstyle{\rm B}}T_{\max}$, Chebyshev expansion requires the numbers of 140, 300, and 800 polynomials, respectively, which are still larger than the corresponding Krylov polynomials $v$ at $T=0$. The small size of Krylov subspace $v$ provides the important efficiency for large-scale ES simulations with the rGF method.
cpl-41-5-053102-fig2.png
Fig. 2. Stochastic charge deviation from full DFTB calculation for (H$_2$O)$_{99}$ cluster: [(a), (c)] for $N_{\rm s}=1000$, [(b), (d)] for $N_{\rm s}=4000$. The results with (circles) and without (squares) deflation correction are compared for H [(c), (d)] and O [(a), (b)] atoms.
In the following, we investigate the accuracy of self-consistent rGF ES calculations for H$_2$O clusters with different sizes, in comparison with the deterministic full DFTB calculation. To effectively reduce the statistical error in the rGF method, we adopt the deflation technique (denoted by rGF-deflation), namely $\tilde{\mathcal{D}}= \mathcal{D}_0 + (\mathcal{D}-\mathcal{D}_0)\tilde{I}$, where $\mathcal{D}_0$ denotes a good sparse approximation to $\mathcal{D}$. We use the $D_{{\rm H}_2{\rm O}}$ of each single H$_2$O in the cluster to construct a sparse $\mathcal{D}_0$, and keep it fixed in the self-consistency. All self-consistent calculations are converged to the tolerance $10^{-5}$ in $q_{\scriptscriptstyle{R}}$. Figure 2 presents the Mulliken charge deviation $\Delta q_{\scriptscriptstyle{R}}$ from the full DFTB results for both rGF and rGF-deflation calculations of (H$_2$O)$_{99}$ cluster, and the results for $N_{\rm s}=1000$ and 4000 are compared. As is shown, for both O and H atoms, rGF-deflation results present almost one order of magnitude smaller than the rGF calculations, presenting the important correction to rGF ES calculations. In particular, for $N_{\rm s}=1000$, the rGF calculation presents the values of averaged charge deviation, namely $\langle |\Delta q| \rangle=\sum_{R=1}^{N_{\rm Atom}}|\Delta q_{\scriptscriptstyle{R}}|/N_{\rm Atom}$, 0.02 and 0.038 for the H and O atoms, respectively, compared to the corresponding values 0.002 and 0.006 of rGF-deflation results. By increasing $N_{\rm s}$ to 4000, the charge deviation $\langle |\Delta q| \rangle$ is reduced to 0.01 and 0.02 with the rGF method, and 0.001 and 0.0025 with rGF-deflation for the respective H and O atoms. By comparing the $N_{\rm s}=1000$ and $N_{\rm s}=4000$ results of both rGF and rGF-deflation calculations, we can see that the statistical error presents a $1/\sqrt{N_{\rm s}}$ dependence. In addition, the statistical error in both rGF and rGF-deflation results shows unbiased distribution throughout the water cluster. Nevertheless, in percentage, $\langle |\Delta q| \rangle/q_{\scriptscriptstyle{R}}$ is rather small, for example, for O atom, only 0.6$\%$ and 0.09$\%$ with the respective rGF and rGF-deflation of $N_{\rm s}=1000$ (even smaller with $N_{\rm s}=4000$).
Table 1. Total energy per H$_2$O for the deterministic full DFTB, rGF, and rGF-deflation calculations (with $N_{\rm s}=1000$ and 4000) for different H$_2$O clusters.
System $E$/H$_{2}$O full-DFTB rGF rGF rGF-deflation rGF-deflation
$N_{\rm s}$ 1000 4000 1000 4000
(H$_2$O)$_{99}$ $-110.4922$ $-110.3071$ $-110.3428$ $-110.4904$ $-110.4926$
(H$_2$O)$_{233}$ $-110.5078$ $-110.3473$ $-110.4652$ $-110.4976$ $-110.5065$
(H$_2$O)$_{526}$ $-110.5069$ $-110.4630$ $-110.4927$ $-110.5070$ $-110.5042$
(H$_2$O)$_{995}$ $-110.5103$ $-110.4202$ $-110.4967$ $-110.5130$ $-110.5114$
(H$_2$O)$_{1981}$ $-110.5153$ $-110.4860$ $-110.5162$ $-110.5173$ $-110.5165$
(H$_2$O)$_{4087}$ – – – $-110.5206$ $-110.5203$
(H$_2$O)$_{9984}$ – – – $-110.5229$ $-110.5231$
In Table 1, we present the total energy per H$_2$O for different clusters self-consistently calculated by rGF, rGF-deflation with $N_{\rm s}=1000$ and 4000, together with the benchmark full DFTB results for comparison. Compared to full DFTB, the rGF ($N_{\rm s}=1000$) results present large error in the total energy, for example, 185 meV and 160 meV/H$_2$O in (H$_2$O)$_{99}$ and (H$_2$O)$_{233}$, respectively. One can also easily find that the total energy error for rGF ($N_{\rm s}=1000$) presents the important self-average by the cluster size, namely the error decreases with increasing the cluster size (for example, the error is reduced to 29 meV/H$_2$O in (H$_2$O)$_{1981}$), consistent with Refs. [11,12]. The large total energy error in the rGF ($N_{\rm s}=1000$) results is consistent with the large charge deviation as shown in Fig. 2. However, compared to rGF ($N_{\rm s}=1000$) calculations, the rGF ($N_{\rm s}=4000$) results present significant reduction in the total energy error, presenting better agreement with the full DFTB. For example, the energy errors in rGF ($N_{\rm s}=4000$) are 42.6 meV and 14.2 meV/H$_2$O in (H$_2$O)$_{233}$ and (H$_2$O)$_{526}$, respectively, much smaller than the values 160 meV and 44 meV/H$_2$O meV of rGF ($N_{\rm s}=1000$). With increasing cluster size to (H$_2$O)$_{1981}$, the energy error in rGF ($N_{\rm s}=4000$) is reduced to 5.9 meV/H$_2$O, presenting the important self-average effect. When it comes to the rGF-deflation ($N_{\rm s}=1000$) calculations, we find that the error in total energy is a few meV, for example, only 1.8 meV/H$_2$O for the cluster (H$_2$O)$_{99}$, which is dramatically improved over the rGF calculations. By increasing $N_{\rm s}$ to 4000, the total energy deviation from full DFTB calculation falls to $\sim$ $1$ meV or even smaller, demonstrating the important effectiveness of deflation. In addition, it is found that, with increasing the cluster size to (H$_2$O)$_{4086}$ and (H$_2$O)$_{9984}$, rGF-deflation calculations with $N_{\rm s}=1000$ and 4000 can present a deviation less than 1 meV. With better choices of $\mathcal{D}_0$[5,6,8,9] or applying the probing technique,[16-20] the accuracy of rGF is expected to be further improved. Finally, the inset of Fig. 1 shows the total CPU time per iteration versus system size of H$_2$O cluster for the rGF-deflation ($N_{\rm s}=1000$, and $v=35$, which is enough to converge all clusters) and full DFTB calculations. We can find that the CPU time of the rGF method scales linearly with the size of H$_2$O clusters, while full DFTB presents almost cubic scaling. The CPU time of rGF-deflation intersects with that of full DFTB at about (H$_2$O)$_{1500}$. For calculating the (H$_2$O)$_{9984}$ cluster, the rGF takes about 1200 S on a single CPU core. For the rGF calculations, significant acceleration is expected by introducing effective preconditioners[50,51] and developing the multiple Krylov subspaces as reported in Ref. [34] to avoid the CG loop O$^{-1}$H-vector product in constructing Krylov subspace for non-orthogonal basis. In summary, we have developed a linear scaling rGF method for simulating large-scale systems. The rGF is defined on a set of random states and is calculated by utilizing the linear-scaling Krylov subspace method. By combining with the DFTB, we have demonstrated that the Krylov subspace can maintain at small size for materials with different gaps, providing the high efficiency of rGF calculations. We find that the rGF self-consistent electronic structure calculation can reach an error of $\sim$ $1$ meV per H$_{2}$O in total energy, in comparison with the deterministic calculation. Due to the wide applicability of GF, the present random GF method is generally applicable with many other electronic structure methods for accelerating calculations of large-scale material properties, including density functional theory, tight-binding model and nonequilibrium Green's function method for nanoelectronics, quantum many-body perturbation methods, and possesses potential applications in other fields. Acknowledgments. Y. Ke acknowledges financial support from the National Natural Science Foundation of China (Grant No. 12227901), and S. Yuan thanks the financial support from the National Natural Science Foundation of China (Grant Nos. 11974263 and 12174291). The authors thank the HPC platform of ShanghaiTech University for providing the computational facility. References Linear scaling electronic structure methods\mathcalO(N) methods in electronic structure calculationsInhomogeneous Electron GasSelf-Consistent Equations Including Exchange and Correlation EffectsDensity-matrix electronic-structure method with linear system-size scalingDirect calculation of electron density in density-functional theoryGeneralization of the density-matrix method to a nonorthogonal basisExpansion algorithm for the density matrixTight-binding electronic-structure calculations and tight-binding molecular dynamics with localized orbitalsExact results for a three-dimensional alloy with site diagonal disorder: comparison with the coherent potential approximationSelf-Averaging Stochastic Kohn-Sham Density-Functional TheoryA Time-Dependent Random State Approach for Large-Scale Density Functional CalculationsStochastic density functional theoryFast algorithm for finding the eigenvalue distribution of very large matricesDeflation as a Method of Variance Reduction for Estimating the Trace of a Matrix InverseA probing method for computing the diagonal of a matrix inverseHierarchical Probing for Estimating the Trace of the Matrix Inverse on Toroidal LatticesGradient-based stochastic estimation of the density matrixAssessment of localized and randomized algorithms for electronic structureExtending Hierarchical Probing for Computing the Trace of Matrix InversesCommunication: Embedded fragment stochastic density functional theoryEquilibrium configurations of large nanostructures using the embedded saturated-fragments stochastic density functional theoryOverlapped embedded fragment stochastic density functional theory for covalently-bonded materialsStochastic density functional theory: Real- and energy-space fragmentation for noise reductionFast and Universal Kohn-Sham Density Functional Theory Algorithm for Warm Dense Matter to Hot Dense PlasmaBreaking the Theoretical Scaling Limit for Predicting Quasiparticle Energies: The Stochastic G W ApproachTime-dependent stochastic Bethe-Salpeter approachSublinear scaling for time-dependent stochastic density functional theoryA stochastic estimator of the trace of the influence matrix for laplacian smoothing splinesEfficient and accurate linear algebraic methods for large-scale electronic structure calculations with nonorthogonal atomic orbitalsAn order- N electronic structure theory with generalized eigenvalue equations and its application to a ten-million-atom systemThe Idea Behind Krylov MethodsGMRES and the minimal polynomialConstruction of tight-binding-like potentials on the basis of density-functional theory: Application to carbonCalculations of molecules, clusters, and solids with a simplified LCAO-DFT-LDA schemeSelf-consistent-charge density-functional tight-binding method for simulations of complex materials propertiesDensity-functional tight-binding for beginnersDensity functional tight bindingModern Software Tools for Scientific ComputingNumerical Optimization of Density Functional Tight Binding Models: Application to Molecules Containing Carbon, Hydrogen, Nitrogen, and OxygenStoichiometric and non-stoichiometric (101̄0) and (112̄0) surfaces in 2H–SiC: a theoretical studyVMD: Visual molecular dynamicsIterative refinement method for the approximate factorization of a matrix inverseRecursive Factorization of the Inverse Overlap Matrix in Linear-Scaling Quantum Molecular Dynamics Simulations
[1] Goedecker S 1999 Rev. Mod. Phys. 71 1085
[2] Bowler D R and Miyazaki T 2012 Rep. Prog. Phys. 75 036503
[3] Hohenberg P and Kohn W 1964 Phys. Rev. 136 B864
[4] Kohn W and Sham L J 1965 Phys. Rev. 140 A1133
[5] Li X P, Nunes R W, and Vanderbilt D 1993 Phys. Rev. B 47 10891
[6] Yang W 1991 Phys. Rev. Lett. 66 1438
[7] Nunes R W and Vanderbilt D 1994 Phys. Rev. B 50 17611
[8] Niklasson A M N 2002 Phys. Rev. B 66 155115
[9] Goedecker S and Teter M 1995 Phys. Rev. B 51 9455
[10] Alben R, Blume M, Krakauer H, and Schwartz L 1975 Phys. Rev. B 12 4090
[11] Baer R, Neuhauser D, and Rabani E 2013 Phys. Rev. Lett. 111 106402
[12] Zhou W and Yuan S 2023 Chin. Phys. Lett. 40 027101
[13] Fabian M D, Shpiro B, Rabani E, Neuhauser D, and Baer R 2019 WWIREs: Comput. Mol. Sci. 9 e1412
[14] Hams A and De Raedt H 2000 Phys. Rev. E 62 4365
[15] Gambhir A S, Stathopoulos A, and Orginos K 2017 SIAM J. Sci. Comput. 39 A532
[16] Tang J M and Saad Y 2012 Numer. Linear Alg. Appl. 19 485
[17] Stathopoulos A, Laeuchli J, and Orginos K 2013 SIAM J. Sci. Comput. 35 S299
[18] Wang Z T, Chern G W, Batista C D, and Barros K 2018 J. Chem. Phys. 148 094107
[19] Moussa J E and Baczewski A D 2019 Electron. Struct. 1 033001
[20] Laeuchli J and Stathopoulos A 2020 SIAM J. Sci. Comput. 42 A1459
[21] Neuhauser D, Baer R, and Rabani E 2014 J. Chem. Phys. 141 041102
[22] Arnon E, Rabani E, Neuhauser D, and Baer R 2017 J. Chem. Phys. 146 224111
[23] Chen M, Baer R, Neuhauser D, and Rabani E 2019 J. Chem. Phys. 150 034106
[24] Chen M, Baer R, Neuhauser D, and Rabani E 2021 J. Chem. Phys. 154 204108
[25] White A J and Collins L A 2020 Phys. Rev. Lett. 125 055002
[26] Neuhauser D, Gao Y, Arntsen C, Karshenas C, Rabani E, and Baer R 2014 Phys. Rev. Lett. 113 076402
[27] Rabani E, Baer R, and Neuhauser D 2015 Phys. Rev. B 91 235302
[28] Gao Y, Neuhauser D, Baer R, and Rabani E 2015 J. Chem. Phys. 142 034106
[29]Mahan G D 2013 Many-Particle Physics (New York: Springer)
[30]Kadanoff L P and Baym G A 1962 Quantum Statistical Mechanics (New York: Benjamin W A, Inc.)
[31] Hutchinson M F 1990 Commun. Stat. -Simul. Comput. 19 433
[32] Teng H, Fujiwara T, Hoshi T, Sogabe T, Zhang S L, and Yamamoto S 2011 Phys. Rev. B 83 165103
[33]Sogabe T 2022 Krylov Subspace Methods for Linear Systems (Singapore: Springer)
[34] Hoshi T, Yamamoto S, Fujiwara T, Sogabe T, and Zhang S L 2012 J. Phys.: Condens. Matter 24 165502
[35]Supplementary material includes more information about the convergence of Krylov subspace and Chebyshev polynomial expansion for H$_{2}$O- and Si-clusters.
[36] Ipsen I C F and Meyer C D 1998 Am. Math. Mon. 105 889
Campbell S L, Ipsen I C F, Kelley C T, and Meyer C D 1996 BIT Numer. Math. 36 664
[37] Porezag D, Frauenheim T, Köhler T, Seifert G, and Kaschner R 1995 Phys. Rev. B 51 12947
[38] Seifert G, Porezag D, and Frauenheim T 1996 Int. J. Quantum Chem. 58 185
[39] Elstner M, Porezag D, Jungnickel G, Elsner J, Haugk M, Frauenheim T, Suhai S, and Seifert G 1998 Phys. Rev. B 58 7260
[40] Koskinen P and Mäkinen V 2009 Comput. Mater. Sci. 47 237
[41] Elstner M and Seifert G 2014 Philos. Trans. R. Soc. A 372 20120483
[42]Balay S, Abhyankar S, Mark F et al. 2023 [Online] https://petsc.org/
[43]Balay S, Abhyankar S, Mark F et al. 2023 PETSc/TAO Users Manual.
[44] Balay S, Gropp W D, McInnes L C, and Smith B F 1997 Efficient Management of Parallelism in Object Oriented Numerical Software Libraries. In: Modern Software Tools for Scientific Computing, Arge E, Bruaset A M, and Langtangen H P (eds) (Birkhäuser Press) pp 163–202
[45]Bock N, Cawkwell M J, Coe J D, Krishnapriyan A, Kroonblawd M P, Lang A, Liu C, Saez E M, Mniszewski S M, Negre C F A, Niklasson A M N, Sanville E, Wood M A, and Yang P 2008 “Latte” [Online] https://github.com/lanl/LATTE
[46] Krishnapriyan A, Yang P, Niklasson A M N, and Cawkwell M J 2017 J. Chem. Theory Comput. 13 6191
[47] Rauls E, Elsner J, Gutierrez R, and Frauenheim T 1999 Solid State Commun. 111 459
[48]Seifert G, Eschrig H, and Bieger W 1986 Z. Phys. Chem. 267 529
[49] Humphrey W, Dalke A, and Schulten K 1996 J. Mol. Graphics 14 33
[50] Niklasson A M N 2004 Phys. Rev. B 70 193102
[51] Negre C F A, Mniszewski S M, Cawkwell M J, Bock N, Wall M E, and Niklasson A M N 2016 J. Chem. Theory Comput. 12 3063