A hybrid method for the parallel computation of Green’s functions

Abstract

Quantum transport models for nanodevices using the non-equilibrium Green’s function method require the repeated calculation of the block tridiagonal part of the Green’s and lesser Green’s function matrices. This problem is related to the calculation of the inverse of a sparse matrix. Because of the large number of times this calculation needs to be performed, this is computationally very expensive even on supercomputers. The classical approach is based on recurrence formulas which cannot be efficiently parallelized. This practically prevents the solution of large problems with hundreds of thousands of atoms. We propose new recurrences for a general class of sparse matrices to calculate Green’s and lesser Green’s function matrices which extend formulas derived by Takahashi and others. We show that these recurrences may lead to a dramatically reduced computational cost because they only require computing a small number of entries of the inverse matrix. Then, we propose a parallelization strategy for block tridiagonal matrices which involves a combination of Schur complement calculations and cyclic reduction. It achieves good scalability even on problems of modest size.

Introduction

Commercial electronic devices are rapidly approaching the scale where quantum mechanical effects are affecting the device characteristics [1]. For example, nano-transistors in the 10–30 nm range are strongly affected by quantum mechanical effects such as tunneling and leakage. To allow for continued scaling of silicon technology, new materials are introduced [2] and new nanomaterials such as nanowires [3] or carbon nanotubes [4] are currently being considered as wires or as active parts in future generations of transistors. Also, more exotic devices using DNA [5] or other small molecules [6] are currently under investigation for future use in electronic devices.

To model such devices, it is necessary to apply a quantum mechanical framework [2]. In a system where current flows from a large reservoir (source) to a large drain, the non-equilibrium Green’s function approach is applicable [7], [8], [9]. This is a very general framework which can be applied to all the systems mentioned above. This method involves the coupled solution of the Schrödinger equation for the quantum wave function and the Poisson equation for the electrostatic potential. The Schrödinger equation can be solved in multiple ways, including using a Kohn–Sham Hamiltonian. The main computational difficulty of this problem is the self-consistent solution of the Kohn–Sham Hamiltonian and the Poisson equation, which usually requires an iterative scheme with many iterations. The solution of this system requires the calculation of a Green’s function which, once discretized, is often written in the matrix form:

$$\mathbf{G}\stackrel{\mathrm{def}}{=}{\mathbf{A}}^{-1}\text{,}\mathrm{where}\mathbf{A}\stackrel{\mathrm{def}}{=}E{\mathbf{S}}^{\mathrm{o}}-\mathbf{H}-{\mathbf{\Sigma }}^{\mathrm{L}}-{\mathbf{\Sigma }}^{\mathrm{R}}$$

where E is an energy point (a scalar), ${\mathbf{S}}^{\mathrm{o}}$ is an overlap matrix, H is the Hamiltonian of the system, and ${\mathbf{\Sigma }}^{\mathrm{L}}$ and ${\mathbf{\Sigma }}^{\mathrm{R}}$ are the self-energies. The Green’s function matrix is computed at many energy points and this is computationally quite expensive even for small matrices. New nanodevices usually require a description at the atomic level; however, their size is large enough that we may have thousands or even millions of atoms in the active part of the device and thus the dimension of A becomes very large. In such cases, direct inversion of A becomes impractical, and the application of the non-equilibrium Green’s function framework requires new efficient parallel algorithms that exploit the sparsity of A. The focus of this paper is the derivation and benchmarking of a new parallel algorithm that efficiently calculates select entries in G and the lesser Green’s function matrix

$${\mathbf{G}}^{<}\stackrel{\mathrm{def}}{=}\mathbf{G}\mathbf{\Gamma }{\mathbf{G}}^{†}$$

where $\mathbf{\Gamma }$ is at this point an arbitrary matrix. Assumptions regarding the non-zero entries of $\mathbf{\Gamma }$ will be formulated later on. The lesser Green’s function matrix is needed to account for non-equilibrium and scattering effects. The notations in this paper are indicated in Table 1.

In most cases not all entries in G and ${\mathbf{G}}^{<}$ are needed. Usually only the diagonal entries of the matrix are required in the iterative process. This leads to a significant reduction in computational cost, which can be realized using a number of different methods [10]. Most of them are related to a technique developed by Takahashi et al. [11], [12], [13]. In this approach, a block LDU factorization of A is computed. Simple algebra shows that:

$$\mathbf{G}={\mathbf{D}}^{-1}{\mathbf{L}}^{-1}+(\mathbf{I}-\mathbf{U})\mathbf{G}$$

$$={\mathbf{U}}^{-1}{\mathbf{D}}^{-1}+\mathbf{G}(\mathbf{I}-\mathbf{L})$$

where L, U, and D correspond, respectively, to the lower block unit triangular, upper block unit triangular, and diagonal block LDU factorization of A. The dense inverse G is treated conceptually as having a block structure based on that of A.

For example, consider the first equation and $j>i$. The block entry ${\mathbf{g}}_{\mathit{ij}}$ resides above the block diagonal of the matrix and therefore $[{\mathbf{D}}^{-1}{\mathbf{L}}^{-1}{]}_{\mathit{ij}}=0$. The first equation can then be written as:

$${\mathbf{g}}_{\mathit{ij}}=-\sum _{k>i}{\mathbf{u}}_{\mathit{ik}}{\mathbf{g}}_{\mathit{kj}}$$

This allows computing the entry ${\mathbf{g}}_{\mathit{ij}}$ if the entries ${\mathbf{g}}_{\mathit{kj}}(k>i)$ below it are known. Two similar equations can be derived for the case $j<i$ and $j=i$:

$${\mathbf{g}}_{\mathit{ij}}=-\sum _{k>i}{\mathbf{g}}_{\mathit{ik}}{\mathbf{\ell }}_{\mathit{kj}}$$

$${\mathbf{g}}_{\mathit{ii}}={\mathbf{d}}_{\mathit{ii}}^{-1}-\sum _{k>i}{\mathbf{u}}_{\mathit{ik}}{\mathbf{g}}_{\mathit{ki}}$$

This approach leads to fast backward recurrences, as shown in [12].

The generalization of Takahashi’s method to computing ${\mathbf{G}}^{<}$ is not straightforward. In particular Eqs. (3), (4) do not extend to this case. In this paper, we will provide a new way of deriving Eqs. (5), (6), (7). This derivation will then be extended to ${\mathbf{G}}^{<}$ for which similar relations will be derived. These recurrences are most efficient when the sparsity graph of $\mathbf{\Gamma }$ (Eq. (2)) is a subset of the graph of A, i.e., ${\mathit{\gamma }}_{\mathit{ij}}\ne 0\Rightarrow {\mathbf{a}}_{\mathit{ij}}\ne 0$.

For the purpose of computer implementation, an important distinction must be made to distinguish what we term 1D, 2D and 3D problems. In principle, all problems are 3 dimensional. However, if the mesh or device is elongated in one direction, say z, then the mesh can be split into “slices” along the z direction. This gives rise to a matrix A with block tridiagonal structure. The problem is then termed 1 dimensional. Similarly, if the mesh is elongated in two directions, the matrix A assumes a block penta-diagonal form and the problem is called 2 dimensional.

This paper is organized as follows. We first describe general relations to compute Green’s functions (Sections 2 Recurrence formulas for the inverse matrix, 3 Sequential algorithm for 1D problems). These are applicable to meshes with arbitrary connectivity. In Section 4, we calculate the computational cost of these approaches for 1D and 2D Cartesian meshes and for discretizations involving nearest neighbor nodes only, e.g., a 3 point stencil in 1D or a 5 point stencil in 2D. We also compare the computational cost of this approach with a recently published method by Li and Darve (the FIND algorithm [14]).

In Section 5, a parallel implementation of the recurrences for 1 dimensional problems is proposed. The original algorithms by Takahashi [11] and Svizhenko [10] are not parallel since they are based on intrinsically sequential recurrences. However, we show that an appropriate reordering of the nodes and definition of the blocks in A lead to a large amount of parallelism. In practical cases for nano-transistor simulations, we found that the communication time was small and that the scalability was very good, even for small benchmark cases. This scheme is based on a combination of domain decomposition and cyclic reduction techniques¹ (see, for example, Varga and Hageman [16]). Section 6 has some numerical results.