Efficient Integration of Green's functions: Ver. 1.1

The original has been published in Phys. Rev. B 75, 035123 (2007)

Taisuke Ozaki, RCIS, JAIST

Abstract:

An efficient and accurate contour integration method is presented for large-scale electronic structure calculations based on the Green function. By introducing a continued fraction representation of the Fermi-Dirac function derived from a hypergeometric function, the Matsubara summation is generalized with respect to distribution of poles so that the integration of the Green function can converge rapidly. Numerical illustrations, evaluation of the density matrix for a simple model Green function and a total energy calculation for aluminum bulk within density functional theory, clearly show that the method provides remarkable convergence with a small number of poles, indicating that the method can be applied to not only the electronic structure calculations, but also a wide variety of problems.

INTRODUCTION

The integration of the Green function associated with the Fermi-Dirac function is one of the most important ingredients for accurate and efficient implementation of electronic structure calculation methods regardless of ab initio or semi-empirical schemes. One can realize that the integration arises in lots of electronic structure calculation methods such as the Korringa-Kohn-Rostoker (KKR) Green function method,[1,2,3] Green function based O($N$) methods such as a recursion method,[4,5] a surface Green function method[6] such as an embedded method,[7] and Green function methods based on a many body perturbation theory such as the GW method.[8,9] Therefore, it is obvious that a highly accurate and efficient method must be developed for the integration to extend applicability of these electronic structure calculation methods. Since the well-known Matsubara summation [15,16] is insufficient for this purpose due to the slow convergence, considerable efforts have been devoted for development of efficient integration methods.[10,11,12,13,14] A way of developing an efficient method beyond the Matsubara summation is to seek another expression for the Fermi-Dirac function. Along this line, an expression $(1+(1+x/n)^n)^{-1}$ is used by Nicholson et al.,[10] and an integral representation is proposed by Goedecker.[11] In this paper I seek a different expression for the Fermi-Dirac function, and present an accurate and efficient method for the integration, based on a continued fraction representation of the Fermi-Dirac function which is derived from a hypergeometric function. Because of a simple distribution of poles in the continued fraction similar to the Matsubara poles, it is anticipated that the proposed method can be easily incorporated into many applications with the same procedure as for the Matsubara summation but with remarkable efficiency. It will be shown that the rapid convergence in the integration can be attributed to an interesting non-uniform distribution of the poles on the imaginary axis, and that even a small number of poles in the contour integration enables us to achieve convergent results with precision of more than ten digits. This paper is organized as follows: In Sec. II the theory of the proposed method is presented with comparison to other methods. In Sec. III in order to illustrate the rapid convergence properties of the proposed method I give two numerical examples: a calculation of the density matrix for a simple model Green function and a total energy calculation of aluminum bulk within a density functional theory. In Sec. IV the theory and applicability of the proposed method are summarized.

THEORY

Continued fraction representation of the Fermi-Dirac function

Let us introduce a continued fraction representation of the Fermi-Dirac function $f(x)$ which is derived from a hypergeometric function:

$$\frac{1}{1+\exp(x)} = \frac{1}{2} - \frac{x}{4} ( \frac{\strut 1}... ... (\frac{x}{2})^2} {\displaystyle \frac{\strut \cdots} {(2M-1)+}_{\ddots}, }}}})$$ (1)

where $x=\beta(z-\mu)$ with $\beta=\frac{1}{k_{\rm B}T}$, $\mu$ is chemical potential, $T$ electronic temperature, $z$ and $x$ are complex variables. The derivation of Eq. (1) is given in the Appendix A. Figure 1 shows the approximate Fermi-Dirac function terminated at 60 levels of the continued fraction given by Eq. (1) as a function of $x$ together with the exact Fermi-Dirac function, $(1+(1+x/n)^n)^{-1}$ with $n=120$, and the Matsubara expansion. Note that all the approximants possess the same number of poles, 120 poles, on the whole complex plane for substantial comparison. It is found that the approximate function by Eq. (1) can reproduce accurately the Fermi-Dirac function for a wide range of energy compared to the function $(1+(1+x/n)^n)^{-1}$ which has been previously used to approximate the Fermi-Dirac function in a contour integration method,[10] while the Matsubara expansion poorly approximates the Fermi-Dirac function when curtailed at a finite number of poles. [17]

Now I consider the termination of the continued fraction Eq. (1) to a finite level $M$, and assume $M$ is an even number, i.e., $M=2N (N=1,2,3,...)$. It can be pointed out that Eq. (1) with the finite continued fraction may correspond to the Padé approximant $P_{M}(x)/Q_{M}(x)$ of the Fermi-Dirac function itself rather than the exponential function, where $P_{M}(x)$ and $Q_{M}(x)$ are relevant polynomials of M-th order. Thus, $P_{M}(x)-Q_{M}(x)f(x)$ is zero up to $2M$-th order of $x$, where $f(x)$ is the Taylor series of the Fermi-Dirac function.[19] In this sense the proposed method is similar to the integral representation of the Fermi-Dirac function proposed by Goedecker,[11] while the overall feature of the derivation looks very different. In case of $N=1$, for instance, Eq. (1) with the finite continued fraction is decomposed into partial fractions via the Padé approximant as follows:

$$\frac{1}{2} - \frac{x}{4}( \frac{\strut 1} {\displaystyle 1 + \frac{\strut (\frac{x}{2})^2} {\displaystyle 3} } ) = \frac{x^2-6x+12}{2x^2+24},$$

$$= \frac{1}{2} + \frac{-\frac{3}{2}}{x+2\sqrt{3}i} + \frac{-\frac{3}{2}}{x-2\sqrt{3}i}.$$ (2)

In general, for the finite continued fraction terminated at $M=2N$, the Fermi-Dirac function can be approximated by the sum of partial fractions with the poles $z_p$ and residues $R_p$ as follows:

$$\frac{1}{1+\exp(x)} \simeq \frac{1}{2} + \sum_{p=1}^{N}\frac{R_{p}}{x-iz_p} + \sum_{p=1}^{N}\frac{R_{p}}{x+iz_p}.$$ (3)

The expression becomes exact in case of the limit: $n\to \infty$. Noting that $1/(1+\exp(x))=1/2-1/2\tanh(x/2)$ and $\tanh(x)$ is an odd function, it is found that the poles are symmetrically located on the complex plane with respect to the real axis. Furthermore, it will be shown in later discussion that the poles are located on only the imaginary axis and that the residues are real numbers. Thus, note that $z_p$ is a real number in Eq. (3). In case of an odd number of $M=2N-1$, an additional term which is proportional to $x$ is added to Eq. (3), and I will not discuss the odd case in this paper.

Although $z_p$ and $R_p$ may be analytically expressed for any termination to arbitrary order of $N$, however, the formulae tend to be overly complicated to handle as $N$ increases. Thus, I present an alternative method of finding these values by introducing an eigenvalue problem. Taking account of a fact that the (1,1) element $(C^{-1})_{11}$ of the inverse of a tridiagonal matrix $C$ defined by

$$C = \left( \begin{array}{ccccccc} c_{11} & c_{12} & & & & \smash{... ...dots \\ \smash{\hbox{\bg 0}}& & & & c_{M(M-1)} & c_{MM} \\ \end{array}\right)$$ (4)

is given by

$$(C^{-1})_{11} = \frac{\strut 1} {\displaystyle c_{11} - \frac{\st... ... - \frac{\strut c_{34}c_{43}} {\displaystyle \frac{\strut \cdots} {c_{MM}} }}}}$$ (5)

and comparing the continued fraction in the parenthesis of the right hand side of Eq. (1) with Eq. (5), the matrix elements of $C$ can be determined as follows:

$$c_{pp} = 2p - 1,$$ (6)

$$c_{p(p+1)} = c_{(p+1)p} = i\frac{x}{2},$$ (7)

where $p$ runs from 1 to $M$ for $c_{pp}$, and from 1 to $M-1$ for $c_{p(p+1)}$ and $c_{(p+1)p}$. Thus, one obtain an interesting relation that the Fermi-Dirac function can be expressed by the (1,1) element of the resolvent of $(ixB-A)$:

$$\frac{1}{1+\exp(x)} \simeq \frac{1}{2} -\frac{x}{4}\{(ixB-A)^{-1}\}_{11}$$ (8)

with

$$A = \left( \begin{array}{ccccccc} -1 & 0 & & & & \smash{\lower1.7... ...dots & \ldots & \\ \smash{\hbox{\bg 0}}& & & & &-(2M-1) \\ \end{array}\right)$$ (9)

and

$$B = \left( \begin{array}{ccccccc} 0 & \frac{1}{2} & & & & \smash{... ... \ldots \\ \smash{\hbox{\bg 0}}& & & & \frac{1}{2} & 0 \\ \end{array}\right),$$ (10)

where $A$ and $B$ are matrices of $M \times M$ in size, and Eq. (8) becomes exact as $M$ increases. The relation Eq. (8) tells us that finding the singular points of Eq. (1) is equivalent to solving the generalized eigenvalue problem $Ab=ixBb$, where $b$ is the eigenvector. Moreover, replacing $ix$ with $y$ yields $Ab=yBb$. Then, it can be proven that the eigenvalue $x$ is a pure imaginary number, giving that $z_p$ is a real number in Eq. (3), since $y$ is a real number due to the symmetric real matrices $A$ and $B$.[20] The imaginary eigenvalues go away from the real axis and the distribution becomes sparse as the distance between the pole and the real axis increases as shown in the inset of Fig. 2, and this is a key feature for the fast convergence in the integrations of the Green function as shown later, since the Green function becomes structureless, while going away from the real axis. In addition, noting that the residue $R_p'$ calculated from the eigenvectors $b$ for $Ab=yBb$ is a real number, since the eigenvectors $b$ consist of real numbers, and that

$$-\frac{x}{4}\{(ixB-A)^{-1}\}_{11} =$$

$$-\frac{x}{4} \sum_{p=1}^{N} \left( \frac{R_p'}{y-z_p} + \frac{-R_p'}{y+z_p} \right),$$

$$= \sum_{p=1}^{N}\frac{\frac{z_pR_p'}{4}}{x-iz_p} + \sum_{p=1}^{N}\frac{\frac{z_pR_p'}{4}}{x+iz_p},$$

$$= \sum_{p=1}^{N}\frac{R_p}{x-iz_p} + \sum_{p=1}^{N}\frac{R_p}{x+iz_p},$$

one can find that the residue $R_p$ in Eq. (3) is a real number, where $R_p'$ and $-R_p'$ appearing in the first line of Eq. (11) are due to the symmetry of $\tanh(x)$. Thus, the poles and residues in Eq. (1) can be easily obtained by solving the eigenvalue problem $Ab=yBb$ with the real symmetric matrices $A$ and $B$ without any numerical difficulty.[21]

The characteristic feature of Eq. (1) that the poles are located on the imaginary axis and the residues are real is similar to that of the Matsubara summation. It is noted that the Fermi-Dirac function in the Matsubara summation is expressed by a development:

$$\frac{1}{1+\exp(x)} = \frac{1}{2} + \sum_{n=1}^{\infty} \frac{-2x}{x^2+\pi^2(2n-1)^2}$$ (11)

with the poles $\pm i\pi(2n-1)$ and the residues $-1$. This relation Eq. (12) can be derived by taking account of Eq. (A10) and the Mittag-Leffler expansion [22] for $\tanh(x)$. Therefore, one may think that their convergence properties should be comparable to each other. However, it will be shown later that the continued fraction representation by Eq. (1) provides a remarkable convergence speed compared to the Matsubara summation. The roots of the rapid convergence lies on an interesting non-uniform distribution of the poles on the imaginary axis. It is observed from Fig. 3 that in case of $N=100$ the interval between neighboring poles of Eq. (1) are uniformly located up to around 61th poles with the same interval $2\pi$ as that of the Matsubara poles, and from then onward it increases very rapidly as the distance between the pole and the real axis increases, while the Matsubara poles distribute uniformly on the imaginary axis everywhere. The distribution of the poles seems to be almost independent of the number of poles $N$, provided that the pole index is normalized by the number of poles $N$ on the upper half complex plane, i.e., the interval starts to increase rapidly at about 61 % of the total number of poles $N$. The sparse distribution of poles of Eq. (1) in the faraway region from the real axis well matches the fact that the Green function becomes a smooth function, while going away from the real axis, and by covering a wide range of the imaginary axis with a small number of poles, the continued fraction representation Eq. (1) can reproduce very well the Fermi-Dirac function on the real axis as shown in Fig. 1. The property is the reason why Eq. (1) is much superior to the Matsubara summation for the integration of the Green function.

Calculation of the density matrix

We are now ready to make full use of Eq. (3) to evaluate integrated quantities associated with a Green function with the Fermi-Dirac function. Here, as an example, I consider evaluation of the one-particle density matrix $\rho$, which is one of the most important quantities in electronic structure calculations, calculated by the Green function $G(z)$ being an analytic function apart from the real axis $E$ with the Fermi-Dirac function. Considering the expression of Eq. (3) and the contour integral with an integration path on the upper half complex plane as shown in Fig. 2, the density matrix $\rho$ can be evaluated by

$$\rho = \Im\left[-\frac{2}{\pi}\lim_{R\to \infty} \int_{-R+\mu}^{R+\mu} dE G(E+i0^+)f(E)\right],$$

$$= \Im\left[ \frac{1}{\pi}\int_{c} dz G(z) +\frac{2}{\beta\pi} \sum_{p=1}^{N} \oint_{{\bar c}_p} dz G(z) \frac{R_{p}}{z-\alpha_p} \right],$$

$$= \mu^{(0)} + \Im\left[ -\frac{4i}{\beta} \sum_{p=1}^{N} G(\alpha_p)R_p \right],$$ (12)

where the factor 2 in the first line is for spin multiplicity, $\alpha_{p} = \mu + i\frac{z_p}{\beta}$, $\mu^{(0)}$ the zeroth order moment of the Green function, and $0^+$ a positive infinitesimal. This expression Eq. (13) is valid for Green functions such as $G(z)=(zS-H)^{-1}$ and $G(z)=(zS-H-\Sigma(z))^{-1}$, where $H$ and $S$ are real symmetric Hamiltonian and overlap matrices, and $\Sigma(z)$ is an energy-dependent matrix such as a self-energy matrix which is real symmetric on the real axis. A more general case with energy dependent complex matrices is discussed in the Appendix B. The first term in the second line of Eq. (13) is transformed into the zeroth order moment $\mu^{(0)}$ of the Green function by the contour integration as follows:

$$\Im\left[ \frac{1}{\pi}\int_{c} dz G(z) \right] =$$

$$\Im\left[ -\frac{1}{\pi}\int_{-\infty}^{\infty}dE G(E+i0^+) \right],$$

$$= \mu^{(0)}.$$ (13)

Although the evaluation of the second term in the final line of Eq. (13) is straightforward, the calculation of the first term $\mu^{(0)}$ is nontrivial. However, it is shown that the zeroth order moment can be easily evaluated by utilizing the moment representation of the Green function as shown below. Starting from the Lehmann representation of the Green function $G(z)=\int_{-\infty}^{\infty} dE \frac{g(E)}{z-E}$ and taking account of $(z-E)^{-1}=z^{-1}(1-E/z)^{-1}=\sum_{p=0}^{\infty}\frac{E^p}{z^{p+1}}$ for $\vert E/z \vert<1$, one can obtain the moment representation of the Green function:

$$G(z)= \sum_{p=0}^{\infty} \frac{\mu^{(p)}}{z^{p+1}},$$ (14)

where $\mu^{(p)}$ is the $p$th moment of the Green function defined by

$$\mu^{(p)} = \int dE E^{p} g(E).$$ (15)

On the other hand, multiplying the both side of Eq. (15) by $z^q$ and integrating from $-\infty$ to $\infty$ along the real axis with an infinitesimal imaginary part $0^+$, one has an alternative expression for the moment:

$$\mu^{(q)} = \Im\left[ -\frac{1}{\pi}\int_{-\infty}^{\infty} dE E^{q} G(E+i0^+) \right].$$ (16)

It is found that the case of $q=0$ in Eq. (17) corresponds to Eq. (14). Writing $z=R\exp(i\theta)$ in Eq. (15), assuming that $R$ is large enough so that the terms associated with higher orders of $R$ can be negligible, and taking $\theta=\frac{\pi}{2}$, one can obtain an approximate expression for $\mu^{(0)}$:

$$\mu^{(0)} \simeq iR~G(iR),$$ (17)

where $R$ is a large real number and 1.0e+10 is used in this study. Therefore, one can see that just one evaluation of the Green function yields the zeroth order moment $\mu^{(0)}$. In Eq. (13) the integration on the real axis is transformed into a quadrature on an axis perpendicular to the real axis with the characteristic distribution of the poles except for the evaluation of the zeroth moment $\mu^{(0)}$. Thus, it is reasonable to anticipate that the density matrix evaluated by Eq. (13) quickly converges as a function of the number of poles, since the Green function becomes structureless by going away from own poles located on the real axis.

In case of a general density matrix

In this subsection an expression is given for the evaluation of a general density matrix $\rho^{(n)}$ defined by

$$\rho^{(n)} = \Im\left[-\frac{2}{\pi}\lim_{R\to \infty} \int_{-R+\mu}^{R+\mu} dE E^{n} G(E+i0^+)f(E)\right],$$ (18)

where it is assumed that the Green function $G$ possesses the same property as discussed for Eq. (13). It is worth mentioning that $\rho^{(0)}$ is a conventional density matrix given by Eq. (13), and $\rho^{(1)}$ is the so called energy density matrix. Noting that $z^n=(z-\alpha_p)(\sum_{m=0}^{n-1}z^{n-1-m}\alpha_p^m) +\alpha_p^n$, and taking account of Eq. (3), one has the following expression:

$${ \frac{z^n}{1+\exp(\beta(z-\mu))} \simeq}$$

$$\frac{1}{2}z^n \nonumber + \sum_{p=1}^{N} \frac{R_p}{\beta} \sum_{m=0}^{n-1} (\alpha_p^{n-1-m}+(\alpha_p^*)^{n-1-m}) z^m$$

$$\qquad + \sum_{p=1}^{N} \frac{\alpha_p^n \frac{R_{p}}{\beta}}{z-\alpha_p} + \sum_{p=1}^{N} \frac{(\alpha_p^*)^n \frac{R_{p}}{\beta}}{z-\alpha_p^*}.$$ (19)

This expression becomes exact in case of $N\to \infty$. Therefore, putting Eq. (20) into Eq. (19) and considering a contour integration with the same path given by Fig. 2 as discussed for Eq. (13), we see that $\rho^{(n)}$ can be evaluated by

$$\rho^{(n)} = \mu^{(n)} + 2\sum_{m=0}^{n-1} \gamma_m \mu^{(m)} + \Im\left[ -\frac{4i}{\beta} \sum_{p=1}^{N} G(\alpha_p)R_p \alpha_p^n \right]$$

with

$$\gamma_m = \sum_{p=1}^{N} \frac{R_p}{\beta} (\alpha_p^{n-1-m}+(\alpha_p^*)^{n-1-m}).$$ (20)

The moments $\mu^{(m)}$ up to $n$th orders in Eq. (21) can be evaluated by generalizing the way of calculating $\mu^{(0)}$ discussed in the previous subsection. If the summation in the moment representation Eq. (15) of the Green function is terminated at the $n$th level, the moments $\mu^{(m)}$ up to $n$th orders can be determined by solving a linear equation for each element (i,j):

$$\left( \begin{array}{ccccc} 1 & z_0^{-1} & \cdots & z_0^{-n} \\ ... ...ij}(z_0)\\ z_1 G_{ij}(z_1)\\ \cdots\\ z_n G_{ij}(z_n)\\ \end{array}\right),$$ (21)

where $z_m(m=0,1..,n)$ are large complex numbers, and can be given by $z_m = R\exp(i\frac{\pi}{2}(1+\frac{m}{n}))$ with a large $R$, while they can be chosen arbitrarily only if $\vert z_m \vert$ is large enough so that the higher order terms can be negligible. Thus, one can see that the general density matrix $\rho^{(n)}$ is evaluated by the proposed method as well as the conventional density matrix. In the Appendix C the integrations of the Green function with the derivative of the Fermi-Dirac function with respect to energy are discussed.

NUMERICAL RESULTS

To demonstrate the rapid convergence of the proposed method I give two numerical examples, (i) calculation of the density matrix for a simple model Green function and (ii) a total energy calculation of aluminum bulk in the face centered cubic (FCC) structure within the Harris density functional.[23]

As a simple model, let us consider a model Green function given by

$$G(z) = \frac{1}{z+10} + \frac{1}{z+5} + \frac{1}{z+2} + \frac{1}{z-5},$$ (22)

where the unit of energy is in eV. It is assumed that the chemical potential is 0 eV and that the electronic temperature is 300 K. Table I show the density matrix $\rho$ calculated by three method: the proposed method through Eq. (13), the method discussed in Ref. [10], and the Matsubara summation.[15] We see that the density matrix calculated by the proposed method rapidly converges at the analytic value (3.0) using only 40 poles, while the convergence speed by the second method is about three times slower than that of the proposed method. Also, the stepwise improvement of the density matrix in the second method is attributed to the plunge step of the approximate function $(1+(1+x/n)^n)^{-1}$ used in the method in the lower energy region as shown in Fig. 1. It is also confirmed that the convergence of the Matsubara summation is quite slow.

As well as the simple model, I show in Table II the convergence of the total energy calculated by the recursion method [5] with the Harris functional[23] for FCC aluminum as a function of the number of poles. In this calculation I use pseudo-atomic double valence numerical functions as basis set,[24] a norm-conserving pseudopotential,[25] and a generalized gradient approximation (GGA) to the exchange-correlation potential.[26] The Green function is evaluated by the recursion method with ten recursion levels.[5] Also, an electronic temperature of 600 K is used. The DFT calculation was done using the OpenMX code.[27] As shown in Table II, it is obvious that the proposed method provides remarkable rapid convergence. The use of only 80 poles yields the convergent result of 14 digits which corresponds to the limit of double precision. On the other hand, it turns out that the second and third methods are much slower than the proposed method, and that the full convergence is not achieved even in case of ten thousands poles. It is also worth mentioning that the proposed method can be applicable to insulators and semi-conductors as well as metallic systems without depending the magnitude of band gap in a system if the Fermi energy is determined so that the total number of electrons in the whole system can be conserved.

CONCLUSIONS

I have developed an efficient contour integration method for electronic structure calculation methods based on the Green function. By introducing a continued fraction representation of the Fermi-Dirac function, the integration of Green function with the Fermi-Dirac function on the real axis is transformed into a quadrature on an axis perpendicular to the real axis except for the evaluation of the moments. Although the feature is similar to that of the Matsubara summation, the most important distinction between those lies on the distribution of poles. The poles of the continued fraction distribute uniformly on an axis perpendicular to the real axis with the same interval as the Matsubara poles up to about 61 % of the total number of poles, and from then onward the interval between neighboring poles increases rapidly as the poles go away from the real axis, while the Matsubara poles are uniformly located on an axis perpendicular to the real axis. The characteristic distribution of poles in the continued fraction well matches the smoothness of the Green function in the distant region from the real axis, and therefore enables us to achieve the remarkable convergence using a small number of poles. The numerical examples clearly demonstrate the computational efficiency and accuracy of the proposed method. Since the contour integration by the continued fraction representation can be regarded as a natural generalization of the Matsubara summation with respect to the distribution of poles, it is anticipated that the proposed method can be applied to not only the electronic structure calculation methods, but also many other applications.

Acknowledgments

The author was partly supported by CREST-JST and the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan.

Derivation of the continued fraction representation of the Fermi-Dirac function

In this Appendix the continued fraction representation Eq. (1) of the Fermi-Dirac function is derived using a hypergeometric function. The generalized hypergeometric function is defined as a power series expansion by

$${_p\hspace{-0.5mm}F_q}(a_1,\cdots,a_p;b_1,\cdots,b_q;x) = \sum_{n=0}^{\infty} \frac{\alpha_n x^n} {n!}$$ (A.1)

with

$$\alpha_n = \frac{(a_1)_n(a_2)_n\cdots(a_p)_n} {(b_1)_n(b_2)_n\cdots(b_q)_n},$$ (A.2)

where $(a_1)_n$ is the Pochhammer symbol defined by

$$(a_1)_n = a_1(a_1+1)(a_1+2)\times \cdots \times(a_1+n-1)$$ (A.3)

with an exceptional definition of $(a_1)_0=1$ in case of $n=0$. The special cases of Eq. (A1) yield the following relations:

$$\sinh(x) = x~{_0\hspace{-0.5mm}F_1}(\frac{3}{2},\frac{x^2}{4}),$$ (A.4)

$$\cosh(x) = {_0\hspace{-0.5mm}F_1}(\frac{1}{2},\frac{x^2}{4}).$$ (A.5)

For later discussion, I will drop the suffixes of the hypergeometric function for simplicity of notation, and define to be $F\equiv {_0\hspace{-0.5mm}F_1}$. Furthermore, let us introduce a recurrence relation for $F$:

$$F(a-1,x) = F(a,x) + \frac{x}{(a-1)a} F(a+1,x).$$

This relation can be proven by comparing the term in the power series piece by piece. Dividing both the sides of Eq. (A6) by $F(a,x)$ and considering the reciprocal, one has

$$\frac{F(a,x)}{F(a-1,x)} = \frac{1}{1+\displaystyle\frac{x}{(a-1)a}\displaystyle\frac{F(a+1,x)}{F(a,x)}}.$$ (A.6)

Replacing the ratio of $F$ in the right hand side of Eq. (A7) by Eq. (A7) with increment of $a$ by one recursively, one can obtain the following continued fraction:

$$\frac{F(a,x)}{F(a-1,x)} = \frac{1} { 1+ \displaystyle \frac{ \fra... ...}{a(a+1)} } { 1+ \displaystyle \frac{ \frac{x}{(a+1)(a+2)} } {}_{\ddots}. } } }$$ (A.7)

Now $\tanh(x)$ can be expressed by a continued fraction using Eqs. (A4), (A5), and (A8) as follows:

$$\tanh(x) = \frac{\sinh(x)}{\cosh(x)},$$

$$= x\left( \frac{F(\displaystyle\frac{3}{2},\displaystyle\frac{x^2}{4})} {F(\displaystyle\frac{3}{2}-1,\displaystyle\frac{x^2}{4})} \right),$$

$$= \frac{x} { 1+ \displaystyle \frac{ x^2 } { 3+ \displaystyle \frac{ x^2 } { 5+ \displaystyle \frac{ x^2 } {}_{\ddots}. } } }$$ (A.8)

Noting that $\exp(x)=\frac{1+\tanh(\frac{x}{2})}{1-\tanh(\frac{x}{2})}$, one has

$$\frac{1}{1+\exp(x)} = \frac{1} { 1+\displaystyle\frac{1+\tanh(\frac{x}{2})}{1-\tanh(\frac{x}{2})} },$$

$$= \frac{1}{2} - \frac{1}{2}\tanh(\frac{x}{2}).$$ (A.9)

Inserting Eq. (A9) into Eq. (A10) yields the continued fraction representation Eq. (1) of the Fermi-Dirac function.

In case of the Green function consisting of complex matrices

If the Green function consists of energy dependent complex matrices, a more general expression, different from Eq. (13), is used to calculate the density matrix. As such a case, here let us assume that the Green function consists of matrices depending on energy $z$ and k-point $\bf k$, i.e.,

$$G^{(\bf k)}(z) = (zS^{(\bf k)}-H^{(\bf k)}-\Sigma^{(\bf k)}(z))^{-1},$$ (B.1)

$$= \int dE \frac{g^{(\bf k)}(E)}{z-E},$$ (B.2)

$$= \sum_{n=0}^{\infty} \frac{\mu^{(n,\bf k)}}{z^{n+1}},$$ (B.3)

where $z$ and $E$ are complex and real variables, respectively, and the second and third lines are the Lehmann and moment representations of the Green function. It is assumed that poles of the Green function Eq. (B1) are located on the real axis, which is verified at least in case the self-energy matrix $\Sigma$ is calculated using the surface Green function.[6] Considering that off-diagonal elements of $g^{(\bf k)}(E)$ are complex in general, the density of states $n^{(\bf k)}(E)$ depending on ${\bf k}$ is given using both the retarded and advanced Green functions, $G^{(\bf k)}(E+i0^+)$ and $G^{(\bf k)}(E-i0^+)$, by

$$n^{(\bf k)}(E) = \Im\left[ \frac{-1}{\pi} \left\{ G^{(\bf k)}(E+i0^+) - G^{(\bf k)}(E-i0^+) \right\}\right],$$

where the factor 2 due to spin multiplicity was considered. Thus, the density matrix $\rho$ can be evaluated by

$$\rho = \frac{1}{V_{\rm c} } \int_{\rm BZ} dk^3 \int_{-\infty}^{\infty} dE n^{(\bf k)}(E) f(E),$$

$$= \frac{1}{V_{\rm c}} \int_{\rm BZ} dk^3 \left\{ \rho_r^{(\bf k)} - \rho_a^{(\bf k)} \right\},$$ (B.4)

where $V_{\rm c}$ is the volume of the unit cell, $\int_{\rm BZ}$ represents the integration over the first Brillouin zone, and $\rho_r^{(\bf k)}$ and $\rho_a^{(\bf k)}$ are defined by

$$\rho_r^{(\bf k)} = \Im\left[ \frac{-1}{\pi} \int_{-\infty}^{\infty} dE G^{(\bf k)}(E+i0^+) f(E) \right],$$ (B.5)

$$\rho_a^{(\bf k)} = \Im\left[ \frac{-1}{\pi} \int_{-\infty}^{\infty} dE G^{(\bf k)}(E-i0^+) f(E) \right].$$ (B.6)

Thus, we see that from the Green function defined by Eq. (B1) the density matrix $\rho$ can be easily evaluated using the proposed method, since $\rho_r^{(\bf k)}$ and $\rho_a^{(\bf k)}$ are computed in the same way as for Eq. (13) by considering the contour integrals in the upper and lower complex planes, respectively.

In case of the derivative of the Fermi-Dirac function

The calculations of physical quantities such as electric and thermal conductances require integrals with the derivative of the Fermi-Dirac function with respect to energy. In the Appendix I give the outline for applying the proposed method to such cases. By making use of the continued fraction representation Eq. (1) of the Fermi-Dirac function, the derivative of the Fermi-Dirac function can be expressed as follows:

$$\frac{df(z)}{dz} = \beta \left( -f(z) + \left(f(z)\right)^2 \right),$$

$$\simeq -\frac{\beta}{4} + \frac{1}{\beta} \sum_{p,p'=1}^{2N} \frac{R_pR_{p'}}{(z-\alpha_p)(z-\alpha_{p'})},$$

$$= -\frac{\beta}{4} + 2\sum_{p=1}^{2N}\frac{R_pw_p}{z_p-\alpha_p} + \frac{1}{\beta}\sum_{p=1}^{2N}\frac{R_p^2}{(z_p-\alpha_p)^2},$$

where $\alpha_{N+i}=(\alpha_{i})^*$, $R_{N+i}=R_{i}$, and $w_p$ is given by

$$w_p = \sum_{p'\neq p}^{2N}\frac{iR_p'}{z_p'-z_p}.$$ (C.1)

It is found from Eq. (C1) that the derivative of the Fermi-Dirac function approximated by the finite continued fraction possesses the second order poles of the continued fraction as well as the first order poles. Therefore, if the Green function is analytic around the pole, by Taylor expanding it at the pole $\alpha_p$ as

$$G(z) = G(\alpha_p) + G'(\alpha_p)(z-\alpha_p) + \cdots,$$ (C.2)

one has to include contributions of the derivative of the Green function $G'$ with respect to energy which can be evaluated through a response function calculated from the Green functions.

Bibliography

1

J. Korringa, Physica 13, 392 (1947).

2

W. Kohn and N. Rostoker, Phys. Rev. 94, 1111 (1954).

3

N. Papanikolaou, R. Zeller, and P. H. Dederichs, J. Phys.: Condens. Matter 14, 2799 (2002).

4

R. Haydock, V. Heine, and M.J. Kelly, J. Phys. C 5, 2845 (1972); 8, 2591 (1975); R. Haydock, Solid State Phys. 35, 216 (1980).

5

T. Ozaki and K. Terakura, Phys. Rev. B 64, 195126 (2001).

6

M.P. Lopez Sancho, J.M. Lopez Sancho, J.M.L. Sancho, and J. Rubio, J. Phys. F 15, 851 (1985).

7

J. Inglesfield, J. Phys. C 14, 3795 (1981).

8

L. Hedin, Phys. Rev. A 139, 796 (1965).

9

F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61, 237 (1998).

10

D.M.C. Nicholson, G.M. Stocks, Y. Wang, W.A. Shelton, Z. Szotek, and W.M. Temmerman, Phys. Rev. B 50, 14686 (1994).

11

S. Goedecker, Phys. Rev. B 48, 17573 (1993).

12

K. Wildberger, P. Lang, R. Zeller, and P.H. Dederichs, Phys. Rev. B 52, 11502 (1995).

13

V. Drchal, J. Kudrnovský, P. Bruno, I. Turek, P.H. Dederichs, and P. Weinberger, Phys. Rev. B 60, 9588 (1999).

14

F. Gagel, J. Comput. Phys. 139, 399 (1998).

15

T. Matsubara, Prog. Theor. Phys. 14, 351 (1955).

16

G. D. Mahan, Many Particle Physics, (Plenum, New York, 1981).

17

In addition to Eq. (1), I investigated the convergence properties of other continued fraction representations of the Fermi-Dirac function, which are given in Ref. [17], and it turned out that Eq. (1) is the best convergent approximant among these continued fraction representations. However, there is a possibility of further improvement for the convergence of the continued fraction representation, e.g., a continued fraction derived from a multipoints Padé approximant.

18

L. Lorentzen and H. Waadeland, Continued Fractions with Applications, Amsterdam, Netherlands: North-Holland, 1992.

19

Although the rigorous proof remains unavailable, I confirmed validity of the statement for the cases up to $M=10$ by MATHEMATICA.

20

Consider an eigenvalue problem $A'b=yB'b$ having the same eigenvalues and eigenvectors as for $Ab=yBb$, where $A'=-A$ and $B'=-B$. Noting that $A'$ is positive definite, the eigenvalue problem can be transformed into $[A'^{-1/2}B'A'^{-1/2}]^{-1}d=yd$, where $A'^{-1/2}$ is a diagonal matrix of which diagonal elements are the inverse of square root of the corresponding diagonal element of $A'$, and $d=A'^{1/2}b$. If the dimension of $B'$ is an even number $2N$, the determinant of $B'$ is $(-1/4)^N$, while it is zero in case of an odd number, $[A'^{-1/2}B'A'^{-1/2}]^{-1}$ in this case is a well defined real symmetric matrix. Thus, it is found that the eigenvalues $y$ and eigenvectors $b$ of $Ab=yBb$ are purely real, provided that the dimension of $B'$ is an even number $2N$.

21

Tables for the poles and residues, and a program code of generating those are freely available on a web site (http://staff.aist.go.jp/t-ozaki/) under the term of GNU-GPL.

22

I.S. Gradshtein, I.M. Ryzhik, Table of integrals, series, and products, A. Jeffrey, Ed. (Academic Press, New York, 1980).

23

J. Harris, Phys. Rev. B 31, 1770 (1985).

24

T. Ozaki, Phys. Rev. B. 67, 155108 (2003); T. Ozaki and H. Kino, ibid. 69, 195113 (2004).

25

N. Troullier and J.L. Martins, Phys. Rev. B 43, 1993 (1991).

26

J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).

27

The code, OpenMX, pseudo-atomic basis functions, and pseudopotentials are available on a web site (http://www.openmx-square.org/).

Table: The density matrix of a model Green function given by Eq. (24) as a function of the number of poles in the upper half plane used in three contour integration methods. The analytic value is 3.0.

Poles	Proposed	$\frac{1}{1+(1+\frac{x}{n})^n}$	Matsubara
10	2.897457365704	0.000000000704	2.268430836092
20	2.999785910601	0.938582044963	2.424349652146
30	2.999999992975	1.000000000000	2.520372160464
40	3.000000000000	1.000000000001	2.588358024187
60	-	2.000000000000	2.681479492145
80	-	2.000000000000	2.742612466824
100	-	2.999998800183	2.785347036205
120	-	3.000000000000	2.816567540132
200	-	-	2.885375367071
500	-	-	2.953166094829
5000	-	-	2.995297020881

Table: The total energy (Hartree) of FCC aluminum bulk calculated by the recursion method with the Harris functional as a function of the number of poles in the upper half plane used in three contour integration methods, where the electronic temperature is 600 (K) corresponding to a smearing of 0.052 (eV). It is noted that the width of the valence band consisting of $3s$ and $3p$ orbitals up to the Fermi level is about 10 (eV).

Poles	Proposed	$\frac{1}{1+(1+\frac{x}{n})^n}$	Matsubara
10	-42.933903047211	-33.734015919550	-39.612354360046
20	-47.224346653790	-33.623477214678	-39.849746603905
40	-48.323790725570	-33.346245616679	-40.216055898502
60	-48.324441992259	-33.143128624551	-39.676965494522
80	-48.324441994952	-32.870752577236	-43.523770052176
150	-48.324441994952	-33.837428496424	-41.836938942518
200	-	-33.418012271726	-42.543354202255
250	-	-34.003411636691	-43.024756221080
300	-	-34.003236479262	-43.466729654170
350	-	-48.324440028792	-43.834528739677
400	-	-48.324440274509	-44.185100655185
600	-	-48.324440847749	-45.233651519749
1000	-	-48.324441306517	-46.331692884149
2000	-	-48.324441650693	-47.202779497545
5000	-	-48.324441857239	-47.921384128418
10000	-	-48.324441926094	-48.122496320516

Figure: (Color online) The Fermi-Dirac function and its approximants as a function of the real axis $x$. For substantial comparison all the approximants possess the same number of poles, 120 poles.

Figure: The path of the contour integration associated with Eq. (1). The inset shows the distribution of zero points of Eq. (1).

Figure: (Color online) The interval between neighboring poles of Eq. (1) terminated at $N=100$ as a function of an index of poles, where the index of poles is determined in ascending order with respect to the absolute value of poles. The interval for the Matsubara poles is $2\pi$.