Nested Dissection Solver for Transport in 3D Nano-Electronic Devices

TL;DR

This paper extends the Hierarchical Schur Complement method to efficiently solve non-equilibrium Green's functions in three-dimensional nanoscale devices, achieving over 1000x speed-up compared to traditional algorithms.

Contribution

The work applies the HSC-extension to 3D NEGF problems, reducing computational complexity from N^7 to N^6 and enabling practical simulations of large nanoscale systems.

Findings

HSC-extension reduces operation count from N^7 to N^6.

Runtime speed-ups exceed three orders of magnitude.

Enables simulation of large realistic 3D nanoscale devices.

Abstract

The Hierarchical Schur Complement method (HSC), and the HSC-extension, have significantly accelerated the evaluation of the retarded Green's function, particularly the lesser Green's function, for two-dimensional nanoscale devices. In this work, the HSC-extension is applied to determine the solution of non-equilibrium Green's functions (NEGF) on three-dimensional nanoscale devices. The operation count for the HSC-extension is analyzed for a cuboid device. When a cubic device is discretized with $N \times N \times N$ grid points, the state-of-the-art Recursive Green Function (RGF) algorithm takes $\mathcal{O}(N^7)$ operations, whereas the HSC-extension only requires $\mathcal{O}(N^6)$ operations. %Realistic operation counts also depend on the system dimensions in $xyz$-directions and the form of contact self-energy matrix. Operation counts and runtimes are also studied for three-dimensional nanoscale devices of practical interest: a graphene-boron\- nitride-graphene multilayer system, a silicon nanowire, and a DNA molecule. The numerical experiments indicate that the cost for the HSC-extension is proportional to the solution of one linear system (or one LU-factorization) and that the runtime speed-ups over RGF exceed three orders of magnitude when simulating realistic devices, such as a graphene-boron nitride-graphene multilayer system with 40,000 atoms.

Tables

Table 1. Table 1: Operation counts of HSC-extension and RGF for various configurations of cuboid mesh.

Configuration	HSC-extension	RGF
Cubic mesh ($N_x=N_z=N_y=N$)	$𝒪\left(N^6\right)$	$𝒪\left(N^7\right)$
Elongated mesh ($N_x=N_z=N\ll N_y$)	$𝒪\left(N^5{N}_y\right)$	$𝒪\left(N^6{N}_y\right)$
Flattened mesh ($N_z\ll N_x=N_y=N$)	$𝒪\left(N_z^3{N}^3\right)$	$𝒪\left(N_z^3{N}^4\right)$

Table 2. Table 2: Extrapolated CPU timings of transmission calculation for SiNW devices at one energy point with L = 20 𝐿 20 L=20 nm for various shapes.

Shapes	HSC-extension (s)	RGF (s)	speed-up
$L=10D$	62.9	85.3	1.4
$L=6D$	1,064	3,030	2.8
$L=3D$	19,920	147,706	7.4

Table 3. Table 3: CPU timings of transmission calculation for DNA molecules at one energy point.

Number of Base Pairs	7	9	11	13	15
HSC-extension (s)	5.0	8.7	11.3	12.5	14.6
RGF (s)	5.5	7.8	9.7	11.6	13.7

Table 4. Table 4: CPU timings of transmission calculation for DNA molecules at one energy point.

Number of Base Pairs	9
HSC-extension (s)	8.7
RGF (s)	7.8
HSC-extension with customized partition (s)	7.8
Customized RGF with two 2-pairs layers (s)	11.8

Equations

$${\mathbf{G}}^r(E)=\mathbf{A}(E{)}^{-1}\text{mbox}with\mathbf{A}(E)=E\mathbf{I}-\mathbf{H}-{\mathbf{\Sigma }}^r$$

$${\mathbf{G}}^{<}(E)={\mathbf{G}}^r(E){\mathbf{\Sigma }}^{<}(E){\left({\mathbf{G}}^r(E)\right)}^{†}$$

$${\mathbf{L}}^{-1}\mathbf{A}{\mathbf{L}}^{-T}=\mathbf{D},$$

$${\mathbf{G}}^r=\left(\mathbf{I}-{\mathbf{L}}^T\right){\mathbf{G}}^r+{\mathbf{D}}^{-1}{\mathbf{L}}^{-1}$$

$${\mathbf{D}}^{-1}{\mathbf{L}}^{-1}{\mathbf{\Sigma }}^{<}{\mathbf{L}}^{-†}{\mathbf{D}}^{-†}$$

$${\mathbf{G}}^{<}=\left(\mathbf{I}-{\mathbf{L}}^T\right){\mathbf{G}}^{<}+{\mathbf{D}}^{-1}{\mathbf{L}}^{-1}{\mathbf{\Sigma }}^{<}{\mathbf{L}}^{-†}{\mathbf{D}}^{-†}{\left({\mathbf{L}}^T\right)}^{-†}$$

$${\mathbf{D}}^{-1}{\mathbf{L}}^{-1}{\mathbf{\Sigma }}^{<}{\left({\mathbf{G}}^r\right)}^{†}$$

$${\mathbf{G}}^{<}=\left(\mathbf{I}-{\mathbf{L}}^T\right){\mathbf{G}}^{<}+{\mathbf{D}}^{-1}{\mathbf{L}}^{-1}{\mathbf{\Sigma }}^{<}{\left({\mathbf{G}}^r\right)}^{†}$$

$$T=\{\begin{array}{ll}\mathcal{O}(N_x^{1.5}{N}_y^{1.5}{N}_z^2)& \text{when}{N}_z\ll N_x\simeq N_y\\ \mathcal{O}(N_x^2{N}_y{N}_z^2)& \text{when}{N}_z\ll N_x\ll N_y\end{array}$$

Full text

1 Introduction

With the downscaling of nanoscale electronic devices, the non-equilibrium Green’s function (NEGF) method has become the most effective and accurate modeling approach in predicting the electronic transport performance. The self-consistent solution of NEGF and Poisson’s equations is capable of accurately modeling coherent as well as decoherent transport by integrating atomic/molecular orbital level physics datta2002non . An aspect of NEGF, that currently prevents a broader use, is the large computational cost associated with the method.

Researchers and analysts have used several approximations to alleviate this cost. Examples include the mode-space approximation ren2003nanomos , which couples 1D NEGF transport simulation to transverse states from solving 2D Schrödinger equation on the cross-section, and the restriction to levels of lower accuracy (such as tight-binding or effective mass levels). Although these approximations of NEGF have successfully predicted the transport characteristics of nanoscale devices barker1989theory ; luisier2006quantum , these alternatives remain unable to capture atomic-scale inhomogeneities, such as surface roughness, unintentional doping, and trapped charges (see, for example, asenov2003simulation for an illustration of the resulting inaccuracy). The modeling of such inhomogeneities requires a full 3D real-space NEGF simulation. Recent studies martinez2007study ; martinez2006development ; martinez2007self ; martinez2009comparison ; martinez2010detailed have performed 3D NEGF simulations to handle these inhomogeneities but with “coarse discretization” (i.e. small number of atomic orbitals or small number of grid points) to control the computational cost. To enhance the predictability of these 3D NEGF simulations, finer discretizations have to be considered. Therefore the computational cost of NEGF needs to be addressed. The goal of this paper is to present a numerical method that significantly reduces the cost of 3D NEGF simulations.

During a NEGF simulation, the main computational bottleneck is the repeated evaluations of the retarded Green’s function, $G^{r}$, and the lesser Green’s function, $G^{<}$, over a wide range of energy values Anantram2008aa . Physical quantities, such as density of states, carrier density, and transmission coefficients, are obtained from the evaluations of these Green’s functions. After discretization of the nanoscale device, these repeated evaluations amount to computing entries in the matrices $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$ that approximate, respectively, the retarded and lesser Green’s functions. Recall that the matrix $\mathbf{G}^{r}$ is defined as

[TABLE]

where $\mathbf{H}$ is the system Hamiltonian matrix and $\mathbf{I}$ is the identity matrix. The self-energy matrix $\boldsymbol{\Sigma}^{r}$ consists of two parts: the contact self-energy $\boldsymbol{\Sigma}^{r}_{C}$ and a self-energy to represent scattering mechanisms $\boldsymbol{\Sigma}^{r}_{scant}$. The matrix $\mathbf{G}^{<}$ is defined as

[TABLE]

where $\boldsymbol{\Sigma}^{<}$ corresponds to the lesser self-energy matrix Anantram2008aa . To compute diagonal and desired off-diagonal entries of $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$, the recursive Green’s function algorithm (RGF) Svizhenko2002aa has been the algorithm of choice for many years. Recently, the hierarchical Schur complement (HSC) method Lin2009ab ; Lin2011ac and the fast inverse using nested dissection (FIND) method Li2008ab ; Li2011ab have exhibited significant speed-ups over RGF to evaluate the diagonal entries of $\mathbf{G}^{r}$. Both of these methods have a smaller operation count than the RGF method. Our HSC-extension to compute diagonal and desired off-diagonal entries of $\mathbf{G}^{<}$ hetmaniuk2013nested has also demonstrated significant speed-ups over RGF simulations for 2D nanoscale devices. In this paper, the efficiency of the HSC-extension hetmaniuk2013nested over the RGF method is studied for 3D nanoscale devices.

The rest of the paper is organized as follows. Section 2 gives a brief description of the RGF and HSC-extension approaches and discusses the asymptotic operation counts for both approaches when simulating 3D devices. Section 3 presents the runtime analysis for a cuboid structure with a regular mesh and three state-of-the-art nanoscale devices: a graphene-boron nitride device with a multilayer hexagonal mesh; a silicon nanowire structure with 10 localized atomic orbitals; and a DNA molecule sequence with an unstructured atomic connectivity. Section 4 provides a brief summary of the numerical experiments.

2 Algorithms overview and cost analysis

This section briefly describes the RGF method Svizhenko2002aa and the HSC-extension approach hetmaniuk2013nested . Then the operation counts for both methods are discussed when simulating cuboid or brick-like nanoscale devices.

2.1 Overview of the RGF algorithm

The RGF algorithm Svizhenko2002aa evaluates the specific entries of the matrices $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$ when the system matrix $\mathbf{A}$ is structured as a block-tridiagonal matrix. The off-diagonal and diagonal blocks of $\mathbf{A}$ are uniquely associated with layers (of atoms or of grid points) orthogonal to the transport direction, usually the $y$-direction. The RGF algorithm is composed of two passes to compute diagonal entries of $\mathbf{G}^{r}$ and two passes to compute diagonal and desired off-diagonal entries of $\mathbf{G}^{<}$. The two passes to compute entries of $\mathbf{G}^{r}$ (or of $\mathbf{G}^{<}$) are interpreted as follows:

• •

the first pass marches one layer at a time from left to right along the $y$-direction and, recursively, folds the effect of left layers into the current layer;

• •

the second pass marches one layer at a time from right to left along the $y$-direction and, recursively, extracts the diagonal blocks and the nearest neighbor off-diagonal blocks for the final result.

Mathematically, when computing diagonal entries of $\mathbf{G}^{r}$, the first pass amounts to transforming $\mathbf{A}$ into a block diagonal matrix $\mathbf{D}$,

[TABLE]

where $\mathbf{L}$ is a lower block-triangular matrix with identity blocks on the diagonal. Recall that $\mathbf{A}$ being a complex symmetric matrix, the block $\mathrm{LDL}^{T}$-factorization is a particular form of a block LU-factorization. The second pass in the evaluation of entries for $\mathbf{G}^{r}$ is expressed by the formula

[TABLE]

(see Takahashi1973 ; Erisman1975 for further details). Starting from the knowledge of the block on the diagonal of $\mathbf{G}^{r}$ for the rightmost layer111The block on the diagonal of $\mathbf{G}^{r}$ for the rightmost layer is the corresponding block in the matrix $\mathbf{D}^{-1}$., the formula (4) provides, one layer at a time, the diagonal block and the nearest neighbor off-diagonal blocks when marching from right to left along the $y$-direction.

When computing $\mathbf{G}^{<}$, the first pass amounts to evaluating the matrix-matrix products,

[TABLE]

where the matrices $\mathbf{L}$ and $\mathbf{D}$ form the block $\mathrm{LDL}^{T}$-factorization of $\mathbf{A}$. The second pass in the evaluation of entries for $\mathbf{G}^{<}$ is expressed by the formula

[TABLE]

(see Takahashi1973 for further details). Starting from the knowledge of the block on the diagonal of $\mathbf{G}^{<}$ for the rightmost layer, the formula (6) provides, one layer at a time, the diagonal block and the nearest neighbor off-diagonal blocks when marching from right to left along the $y$-direction.

2.2 Overview of the HSC algorithm and its extension for $\mathbf{G}^{<}$

The HSC method and its extension Lin2009ab ; Lin2011ac ; hetmaniuk2013nested evaluate specific entries of the matrices $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$ when the matrix $\mathbf{A}$ is symmetric and sparse. In contrast to RGF, the matrix $\mathbf{A}$ is not required to be block-tridiagonal.

The first step of the HSC method and of its extension gathers the grid points (or atoms) into arbitrarily-shaped clusters. These clusters are then organized in a binary tree. The ordering of rows of $\mathbf{A}$ according to this hierarchy of clusters arranges the non-entries of $\mathbf{A}$ into a block arrow-shaped structure. Such choice of ordering allows to fold and to extract in any physical direction when following the vertical hierarchy of the binary tree.

Fig. 1 illustrates a partition for a cuboid device of size $(2a+1)\times(2a+1)\times(2a+1)$ and the corresponding binary tree222In practice, the binary tree is likely to be balanced.. Three levels of separators are depicted in Fig. 1 and the colors in the binary tree match the colors in the partition. Note that the three separators are orthogonal to three distinct physical directions. The fold process represents the computation from the lowest level to the top-most level, while the extract process marches one level at a time from the top-most level to the lowest one.

The mathematical interpretation of the two passes for $\mathbf{G}^{r}$ remains unchanged with the stipulation that the rows of $\mathbf{A}$, $\mathbf{G}^{r}$, and $\mathbf{G}^{<}$ are ordered according to the hierarchy of clusters in the binary tree. The first pass for $\mathbf{G}^{<}$ is interpreted as the evaluation of the matrix-matrix products,

[TABLE]

and the second pass is expressed by the formula

[TABLE]

where the matrix rows remain ordered according to the hierarchy of clusters in the binary tree.

Next we provide two pseudo-codes describing, respectively, the HSC method and its extension. On the basis of the hierarchical structure of a binary tree with $L$ levels, let $P_{i}=\{{\rm clusters}\ j\}$ denote the set of all cluster indices $j$ such that cluster $j$ is an ancestor of cluster $i$. Similarly, let $C_{i}=\{{\rm clusters}\ j\}$ denote the set of all cluster indices $j$ such that cluster $j$ is an descendant of cluster $i$. Algorithm 1 describes the computation of specific entries of $\mathbf{G}^{r}$ (HSC method).

Algorithm 1 (HSC: computation of $\mathbf{G}^{r}$)

for $l=1$ to $L-1$ **do

$\mathbf{A}^{(l)}=\mathbf{A}^{(l-1)}$** $\mathbf{A}^{(0)}=\mathbf{A}$

for all the clusters $i$ on level $l$ **do

$\boldsymbol{\Psi}_{i,j}=-\left(\mathbf{A}_{i,i}^{(l)}\right)^{-1}\mathbf{A}_{i,j}^{(l)}$** for all $j$ in $P_{i}$

$\mathbf{A}_{j,k}^{(l)}=\mathbf{A}_{j,k}^{(l)}+\boldsymbol{\Psi}_{i,j}^{T}\mathbf{A}_{i,k}^{(l)}$ for all $j$ and $k$ in $P_{i}$

$\mathbf{A}_{k,j}^{(l)}=\left(\mathbf{A}_{j,k}^{(l)}\right)^{T}$ for all $j$ and $k$ in $P_{i}$

$\mathbf{A}_{i,j}^{(l)}=\mathbf{0}$ and $\mathbf{A}_{j,i}^{(l)}=\mathbf{0}$ for all $j$ in $P_{i}$

**end for

end for

$\mathbf{G}^{(L-1)}=\left(\mathbf{A}^{(L-1)}\right)^{-1}$** $\mathbf{A}^{(L-1)}$ is block diagonal

for $l=L-2$ to [math] **do

$\mathbf{G}^{(l)}=\mathbf{G}^{(l+1)}$

for** all the clusters $i$ on level $l$ **do

$\mathbf{G}_{i,j}^{(l)}=\mathbf{G}_{i,j}^{(l)}+\sum_{k\in P_{i}}\boldsymbol{\Psi}_{i,k}\mathbf{G}_{k,j}^{(l)}$** for all cluster indices $j$ in $P_{i}$

$\mathbf{G}_{j,i}^{(l)}=\left(\mathbf{G}_{i,j}^{(l)}\right)^{T}$ for all cluster indices $j$ in $P_{i}$

$\mathbf{G}_{i,i}^{(l)}=\mathbf{G}_{i,i}^{(l)}+\sum_{j\in P_{i}}\boldsymbol{\Psi}_{i,j}\mathbf{G}_{j,i}^{(l)}$

**end for

end for

**Computed entries in $\mathbf{G}^{(0)}$ match corresponding entries in $\mathbf{G}^{r}$

Algorithm 2 describes the computation of specific entries of $\mathbf{G}^{<}$ with the HSC-extension. For further details about the HSC method for $\mathbf{G}^{r}$ and about its extension for $\mathbf{G}^{<}$, the reader is referred to the references Lin2009ab ; Lin2011ac ; hetmaniuk2013nested .

Algorithm 2 (HSC-extension: computation of $\mathbf{G}^{<}$)

$\mathbf{N}=\boldsymbol{\Sigma}^{<}\left(\mathbf{G}^{(0)}\right)^{\dagger}$ $\boldsymbol{\Sigma}^{<}$ block diagonal

for $l=1$ to $L-1$ **do

$\mathbf{N}^{(l)}=\mathbf{N}^{(l-1)}$** $\mathbf{N}^{(0)}=\mathbf{N}$

for all the clusters $i$ on level $l$ **do

$\mathbf{N}_{j,k}^{(l)}=\mathbf{N}_{j,k}^{(l)}+\boldsymbol{\Psi}_{i,j}^{T}\mathbf{N}_{i,k}^{(l)}$** for all $j$ and $k$ in $P_{i}$

**end for

end for

$\mathbf{P}^{(L-1)}=\mathbf{G}^{(L-1)}\mathbf{N}^{(L-1)}$

for** $l=L-2$ to [math] **do

$\mathbf{P}^{(l)}=\mathbf{P}^{(l+1)}$

for** all the clusters $i$ on level $l$ **do

$\mathbf{P}_{i,j}^{(l)}=\mathbf{P}_{i,j}^{(l)}+\sum_{k\in P_{i}}\boldsymbol{\Psi}_{i,k}\mathbf{P}_{k,j}^{(l)}$** for all cluster indices $j$ in $P_{i}$

$\mathbf{P}_{j,i}^{(l)}=-\left(\mathbf{P}_{i,j}^{(l)}\right)^{\dagger}$ for all cluster indices $j$ in $P_{i}$

$\mathbf{P}_{i,i}^{(l)}=\mathbf{P}_{i,i}^{(l)}+\sum_{j\in P_{i}}\boldsymbol{\Psi}_{i,j}\mathbf{P}_{j,i}^{(l)}$

**end for

end for

**Computed entries in $\mathbf{P}^{(0)}$ match corresponding entries in $\mathbf{G}^{<}$

The key element controlling the efficiency of the HSC method and its extension is the definition of clusters in the binary tree. When the binary tree is replaced by a degenerate (or pathological) tree where each parent node has only one associated child node, the HSC method becomes equivalent to the RGF method (see, for example, (hetmaniuk2013nested, , p. 596)). For the HSC algorithm and its extension, the multilevel nested dissection defines the binary tree and the resulting clusters. The nested dissection George1973aa approach divides the domain of grid points or atoms into two uncoupled sub-domains (clusters of grid points/atoms) with a separator, such that the two disconnected sub-domains are of roughly the same size and the separator is minimized. Multilevel nested dissection repeats this division on each sub-domain, reducing the clusters down to a size small enough (manageable for an inversion of the corresponding block matrix), or until the clusters cannot be divided by nested dissection. In practice, the METIS graph partition library Karypis:1998:FHQ:305219.305248 provides the implementation of the multilevel nested dissection. This well-established partition library allows the HSC extension approach to handle any arbitrarily-discretized nanoscale electronic devices.

2.3 Operation count analysis when simulating 3D brick-like devices

Next we discuss the operation counts for the RGF method and for the HSC-extension. Consider a cuboid device, covered by a three-dimensional orthogonal mesh with $N_{x}$, $N_{y}$, and $N_{z}$ grid points per direction. The discretization of the Hamiltonian is obtained via a 7-point stencil. The self-energy matrix $\boldsymbol{\Sigma}^{r}$ is assumed to be represented via a similar 7-point stencil (for example, with a crude diagonal approximation or with a PML-like approximation nissen2011optimized ). The resulting matrix $\mathbf{A}$ is of dimension $N_{x}\times N_{z}\times N_{y}$, where the $y$-direction is the transport direction.

The RGF approach groups the grid points into $N_{y}$ disjoint layers, each layer holding $N_{x}\times N_{z}$ grid points. Fig. 2 illustrates these layers for $N_{x}=3$, $N_{z}=3$, and $N_{y}=5$.

By ordering the grid points one layer at a time, the matrix $\mathbf{A}$ exhibits the block-tridiagonal structure, required by the RGF approach. The operation count for the RGF method on this cuboid device is $\mathcal{O}\left(N_{x}^{3}N_{z}^{3}N_{y}\right)$.

The HSC-extension employs a multilevel nested dissection to gather the grid points into a hierarchy of clusters. Fig. 3 illustrates the separators obtained at each level and the eight subdomains. The resulting binary tree is also depicted in Fig. 3 with matching colors for the separators.

The operation count for the HSC-extension is derived in the appendix. For the sake of conciseness, the final value is summarized in Table 1.

Remark

When $N_{z}=1$, the operation counts for RGF and HSC reduce to their expression for 2D devices (with mesh $N\times N$). Namely, for RGF, the operation count becomes $\mathcal{O}(N^{4})$ and, for HSC, $\mathcal{O}(N^{3})$.

In practice the self-energy matrix $\boldsymbol{\Sigma}^{r}$ contains dense blocks for the grid points on open boundary conditions. The analysis in the appendix and the counts (Table 1) do not cover such cases. The next section will study numerically the operation counts for practical nanoscale devices with open boundary conditions.

3 Numerical experiments

Next we demonstrate and analyze numerically the performance of HSC-extension approach when evaluating entries of the matrices $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$. First a cuboid device is used to illustrate the cost analysis presented in Section 2.3. The impact of dense blocks in $\boldsymbol{\Sigma}^{r}$ when modeling open boundary conditions is also discussed. Then three nanoscale devices of practical importance are considered: a graphene-boron nitride-graphene multilayer system, a silicon nanowire (SiNW), and a DNA molecule. The discretizations for these three devices yield matrices with different sparsity pattern and provide different challenges for the HSC-extension.

Both algorithms (RGF and HSC-extension) are implemented as C codes. All the runtime data corresponds to the total CPU time for the evaluation of the diagonal and desired off-diagonal entries of $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$ at a single energy point. For reference, timings for the LU-factorization (lu routine in MATLAB 2011b matlab calling UMFPACK v5.0.4 umfpack ) of $\mathbf{A}$ are also included. The operation count for the LU-factorization represents the optimal cost complexity because every solution of a linear system with $\mathbf{A}$ requires, at least, the cost of one LU-factorization. All numerical experiments are performed with one thread on a machine with Intel i7-2600 3.40GHz CPU and 12GB memory.

3.1 Cuboid nanoscale device

As described in section 2.3, a cuboid nanoscale device is considered where the Hamiltonian is constructed by effective-mass approximation and with 7-point stencil finite difference. A three-dimensional orthogonal mesh is used with $N_{x}$, $N_{y}$, and $N_{z}$ grid points per direction. Two distinct treatments for the self-energy matrices are studied: a diagonal approximation, referred to as SPARSE, and a DENSE approximation for modeling open boundary conditions.

3.1.1 Results with SPARSE self-energy matrices

Here diagonal self-energy matrices are considered. The matrix $\mathbf{A}$ has the same sparsity as the Hamiltonian matrix $\mathbf{H}$. Fig. 4 illustrates the pattern of non-zero entries in the matrix $\mathbf{A}$, when the grid points are ordered one layer at a time (as in Fig. 2).

First, when $N_{x}=N_{y}=N_{z}=N$, the CPU times for evaluating $\mathbf{G}^{r}$ and $\mathbf{G}^{<}$ at one energy point are plotted as a function of $N$ in Fig. 5.

The slopes are consistent with the analysis of section 2.3, namely $\mathcal{O}(N^{7})$ for RGF and $\mathcal{O}(N^{6})$ for HSC-extension. The LU-factorization of $\mathbf{A}$ exhibits also a complexity $\mathcal{O}(N^{6})$. When $N=32$, the HSC-extension exhibits a speed-up of 10 times. Note that, on this 12GB machine, RGF can solve problems only up to $N=32$ (the resulting matrix $\mathbf{A}$ is of dimension 32,768), while the HSC-extension can solve these problems up to $N=40$ (dimension of matrix $\mathbf{A}$ is 64,000).

Next the case of an elongated device is considered, i.e. $N_{x}=N_{z}=N\ll N_{y}$. Fig. 6 illustrates timings for different elongated devices with square cross-section ($N_{x}=N_{z}=N$).

The asymptotic slopes (black dashed curves) match the analysis, namely $\mathcal{O}(N^{5}N_{y})$ for the HSC-extension and $\mathcal{O}(N^{6}N_{y})$ for RGF. Here again the HSC-extension and the LU-factorization have numerically the same complexity. When $N_{y}$ is fixed at 200, the CPU time of HSC-extension is initially higher than the one for RGF at small cross-sections and becomes smaller than the one for RGF when $N_{x}=N_{z}\geq 12$, eventually reaching a speed-up of 2 for the largest structure studied.

Finally, for the case of flattened devices, the CPU results are shown in Fig. 7.

When $N_{z}=4$ and $N_{x}=N_{y}=N\gg N_{z}$, the costs of $\mathcal{O}(N^{3})$ for HSC-extension and $\mathcal{O}(N^{4})$ for RGF are observed. These asymptotic behaviors are consistent with the analysis in section 2.3 and with the conclusions for 2D devices with $N_{x}\times N_{y}$ grid points hetmaniuk2013nested .

Our numerical experiments in Fig. 4-7 illustrate the asymptotic operation count of HSC-extension as a function of system dimensions for various cuboidal shapes. In all three cases, the HSC-extension and the LU-factorization have identical asymptotic operation counts. These numerical experiments strongly suggest that the HSC-extension reaches the ideal complexity for 3D nanoscale devices.

3.1.2 Effect of DENSE self-energy matrices

Next dense self-energy matrices are considered to model open boundary conditions. Fig. 8 illustrates the pattern of non-zero entries in the matrix $\mathbf{A}$, when the grid points are ordered one layer at a time (as in Fig. 2).

RGF does not exploit the sparsity present in most diagonal blocks of matrix $\mathbf{A}$. So using a dense self-energy matrix does not impact the performance of RGF. On the other hand, the HSC-extension aims to exploit as much as possible the sparsity of $\mathbf{A}$. So it is important to study the impact of a dense self-energy matrix on the HSC-extension. The analysis in section 2.3 does not handle this situation.

First consider the case where $N_{x}=N_{y}=N_{z}=N$. Fig. 9 plots the CPU timings as a function of $N$.

Here the RGF calculation stops at $N=32$ due to memory limitation, while the HSC-extension can solve problems up to $N=36$ with dense self-energy matrices. The presence of a dense self-energy matrix yields larger CPU times for the HSC-extension but the asymptotic operation count is not modified.

Fig. 10 plots the CPU timings when $N_{y}$ is modified and $N_{x}$ and $N_{z}$ remain constant. Timings for RGF, for HSC-extension with sparse self-energy matrix, for HSC-extension with dense self-energy matrix, and for the LU-factorization of $\mathbf{A}$ are reported.

Note that in this numerical experiment, the layered-structure decomposition employed in RGF is kept along $y$-direction (even when $N_{y}<N_{x}$). When $N_{y}$ is comparable to $N_{x}$ and $N_{z}$, the speed-up for HSC-extension over RGF is reduced when a dense self-energy matrix is considered. As $N_{y}$ gets larger, the timings of HSC-extension with the two forms of self-energy are closer and, asymptotically, approaching the complexity $\mathcal{O}(N_{y})$.

3.2 Graphene - boron nitride - graphene multilayer system

Graphene, stacked with boron nitride insulating material, is a promising material to build next generation transistors because of its extraordinary thermal and electronic properties gastro2009 . Here a multilayer heterostructure is considered, as shown in Fig. 11.

The device consists of two semi-infinitely long monolayer armchair-edged graphene nanoribbon (AGNR) electrodes sandwiching an ultra-thin hexagonal boron nitride (hBN) multilayer film, yielding a vertical tunneling heterostructure with hBN acting as a potential barrier britnell2013resonant . The hBN film is a few atomic layers thick and the central graphene-hBN-graphene (G-BN-G) overlapping heterostructure/multilayer region is stacked in AB-order (Bernal stacking). For this problem, the number of 2D vertical layers is denoted by $N_{z}$ in units of atomic layers. The system width is $N_{x}$ and the length of the multilayer stacking region is $N_{y}$, also in units of atomic layers. The semi-infinitely long AGNR monolayer electrodes at the top and bottom layers are treated as open boundary conditions, their effect is folded into dense self-energy blocks (extreme blocks) of dimension $N_{x}\times N_{x}$.

The system Hamiltonian is constructed using the nearest neighbor tight binding approximation with parameters from zhao2015negative . Only the low energy $p_{z}$ orbitals are considered here; thus the Hamiltonian has the same dimension as the total number of atoms simulated. The geometric lattice complexity of the multilayer system yields an average 5-point stencil Hamiltonian sparsity (multiple hexagonal-meshed layers stacked in AB order). The complexity of RGF remains $\mathcal{O}\left(N_{x}^{3}N_{z}^{3}N_{y}\right)$ because the sparsity inside each block is not exploited. The complexity of HSC-extension is studied numerically.

Fig. 12 plots the CPU timings when $N_{x}=N_{y}$ and $N_{z}=5$. The HSC-extension and the LU-factorization of $\mathbf{A}$ have the same asymptote, indicating an operation count $\mathcal{O}(N_{x}^{3})$.

Fig. 13 plots the CPU timings for different configurations of $N_{x}$, $N_{y}$, and $N_{z}$. The experiments illustrate that HSC-extension still exhibits a complexity similar to the LU-factorization of $\mathbf{A}$. For the largest devices simulated in Fig. 13(a), $N_{x}=256$, the HSC-extension method offers a speed-up of 3 orders of magnitude over RGF.

The numerical experiments indicate that the asymptotic cost of HSC-extension is $\mathcal{O}\left(N_{x}^{1.5}N_{z}^{2}N_{y}^{1.5}\right)$. This cost of HSC-extension can be compared to the cost (Table 1) for the flattened device. The term $\mathcal{O}(N_{x}^{1.5}N_{y}^{1.5})$ matches well with $\mathcal{O}(N^{3})$ by assuming $N_{x}=N_{y}=N$. The term $\mathcal{O}(N_{z}^{2})$ is due to the Bernal stacking order for the multilayer structure in the $z$-direction.

Finally, Fig. 14 illustrates CPU timings for an elongated device, i.e. $N_{x}\ll N_{y}$. The timings for the HSC-extension behave like $\mathcal{O}(N_{x}^{2}N_{y})$, demonstrating a lower order of complexity over the RGF method.

As a summary, the runtime cost of HSC-extension for the G-BN-G multilayer structure is

[TABLE]

while the runtime cost of RGF behaves like $\mathcal{O}(N_{x}^{3}N_{z}^{3}N_{y})$.

3.3 Silicon nanowire structure

Silicon nanowire devices have shown promises to become key components in the next generation computer chips cui2001functional . Solving efficiently the NEGF equations for such devices is therefore important.

In order to investigate the scaling of computational runtime as a function of SiNW lateral dimensions, specifically the number of atoms in each layer and the number of unit cells, we consider a SiNW device depicted in Fig. 15(a). The number of atomic layers in $y$-direction is denoted as $N_{y}$ and the number of silicon atoms within each atomic layer (cross-section) is denoted as $N_{cs}$. So a $N_{cs}\times N_{y}$ SiNW structure contains an array of $N_{y}/2$ unit cells with $2N_{cs}$ atoms per unit cell.

Next the $sp^{3}d^{5}s^{*}$ tight-binding formalism luisier2006atomistic is used to discretize the system. Each silicon atom is represented by a $10\times 10$ diagonal block, thereby interconnecting with up to 40 orbitals of the nearest-neighbor silicon atoms. The resulting Hamiltonian matrix exhibits a stencil involving more than 40 points, which results in a particular computational challenge. Dense self-energy matrices of dimension $10N_{cs}\times 10N_{cs}$ are employed. Fig. 15(b) illustrates the sparsity of $\mathbf{A}$ when $N_{cs}=108$ and $N_{y}=8$.

In this section, we consider nanowire whose length $L$ is proportional to the diameter $D$ of the cross-section, namely $L=\alpha D$. When ordering the atoms one layer at a time, the Hamiltonian matrix, as well as matrix $\mathbf{A}$, has a block-tridiagonal structure, where each block is of dimension $40N_{cs}\times 40N_{cs}$. The operation count for RGF becomes $\mathcal{O}\left(N_{cs}^{3}N_{y}\right)$. Since $N_{y}=\mathcal{O}(L)$, $N_{cs}=\mathcal{O}(D^{2})$, and $L=\alpha D$, the operation count for RGF is $\mathcal{O}(N_{y}^{7})$. The operation count for HSC-extension will be studied numerically.

Fig. 16 plots CPU timings as a function of $N_{y}$ for structures shaped by $L=10D$, $L=6D$ and $L=3D$. Our numerical experiments exhibit an asymptotic cost of $\mathcal{O}(N_{y}^{7})$ for RGF and $\mathcal{O}(N_{y}^{6})$ for HSC-extension. We would like to emphasize that the complexity $\mathcal{O}(N_{y}^{6})$ is valid for HSC-extension as long as $L=\alpha D$, independent of the value $\alpha$. The HSC-extension has the same asymptotic behavior as the LU-factorization of $\mathbf{A}$.

In practice, analysts may consider nanowires of 20nm length. Table 2 lists CPU timings when the length is 20nm. Because these simulations did not fit in the 12GB RAM of our machine, the CPU timings were extrapolated from the asymptotes in Fig. 16.

This extrapolation indicates that the speedup of HSC-extension over RGF improves as $\alpha$ decreases. The reduction of $\alpha$ enlarges the cross-section $N_{cs}$ while $L$ is fixed, yielding a higher speedup of HSC-extension, which is consistent with the other experiments.

3.4 DNA molecule

Finally we test the algorithm for DNA-based structure, which represents a complex organic system. It has been shown that DNA is one of the promising candidates in the molecule devices Gohler18022011 . The study of electronic structures can be used to develop new sequencing techniques postma2010rapid , acting as DNA fingerprints. Another application is disease detection tsutsui2011electrical . Many diseases are linked with the mutation in DNA bases, resulting in different electronic properties which can be used to distinguish the mutated DNA. In our numerical experiments, the DNA molecule is described by density functional theory (DFT) method. Although the number of atoms contained in single DNA molecule is not huge, the decomposed Hamiltonian matrices are relatively dense, thus impeding an effective decomposition using the multilevel nested dissection.

The DNA molecule in our simulation is a double-helix structure containing 7 – 15 base pairs in each strand sketched in Fig. 17(a).

The Hamiltonian matrices are generated by DFT package GAUSSIAN 09 g09 at HG/6-31G ($d$, $p$) level qi2013unified . The number of orbitals (matrix dimension) for each base is about 250. For example, for a 9-mer DNA molecule, the Hamiltonian is of dimension $4500\times 4500$ as shown in Fig. 17(b). The Hamiltonian can be decomposed by treating each base pair as one layer, yielding a block tri-diagonal shape with 9 layers. Different from the structures studied above, the diagonal and nearest neighbor off-diagonal blocks in the Hamiltonian are fully dense.

The CPU timing results for various DNA molecules are summarized in Table 3.

For these configurations, the HSC-extension seems to be less efficient than RGF.

To better understand these CPU times, it is important to look at the clusters defined by the multilevel nested dissection. Fig. 18 illustrates the clusters used for the HSC-extension when the multilevel nested dissection is applied blindly to the graph of $\mathbf{A}$, and the layers used for RGF.

RGF employs layers with one base pair and the resulting blocks are of dimension $500\times 500$. The multilevel nested dissection works at the level of base pairs because one base pair corresponds to one fully-dense diagonal block in the matrix $\mathbf{A}$. The resulting partition introduces clusters with one base pair except for two bottom level (level 3) clusters that can not be partitioned with nested dissection. These two clusters result in two block matrices of dimension $1000\times 1000$. The time discrepancy arises from these two distinct choices of row gathering.

To further illustrate the impact of the row numbering or partitioning, Table 4 lists timings for two additional approaches.

The four different partitionings used for these approaches are depicted on Fig. 18. The customized RGF method with two 2-pairs layers gathers, twice, 2-pairs into one layer. This particular choice of layers indicates the impact of two blocks of dimension $1000\times 1000$ on the overall CPU time, i.e. an increase in CPU time. The customized partition for the HSC-extension allows a degenerate sub-tree333Any cluster with 2-pairs is partitioned according to a degenerate tree, where each parent node has only one child. to avoid any 2-pairs cluster. This customized partition makes the HSC-extension operate on blocks of dimension, at most, $500\times 500$. This choice results in a lower CPU time, on par with the original RGF approach. A similar behavior was observed for other DNA molecules, where METIS gathers 2-pairs into one cluster.

These numerical experiments suggest that the HSC-extension, combined with a multilevel nested dissection, is an efficient approach even for smaller but denser matrices. When the graph partitioning is allowed to insert degenerate sub-trees, the performance is comparable with that of the performance of RGF.

4 Conclusion

In this work, we demonstrate the HSC-extension based NEGF solver as a working methodology for various 3D systems. The cost analysis for HSC-extension is performed on a cuboid structure. HSC-extension exhibits operation count of $\mathcal{O}(N^{6})$ when simulating cubic device with dimension $N\times N\times N$, whereas a $\mathcal{O}(N^{7})$ count is observed for RGF. We also illustrate various asymptotic costs of HSC-extension when the device has an elongated shape ($N_{x},\,N_{z}\ll N_{y}$), when the device is flattened ($N_{z}\ll N_{x},\,N_{y}$), and when a dense self-energy is used to model the open boundary conditions.

The runtime performance of HSC-extension is further investigated for nano-electronic devices of practical interest: graphene-hBN-graphene multilayer heterostructure, silicon nanowire and DNA molecule. These devices exhibit distinct atomistic sparsity, indicating different computational efficiency for HSC-extension. The numerical experiments suggest that the HSC-extension exhibits asymptotic runtimes and operation counts proportional to the runtime of the LU-factorization. For all the nano-electronic devices considered, the HSC-extension becomes faster than the RGF method as the device gets larger. A 1,000 speed-up is observed for a graphene-hBN-graphene multilayer device with 40,000 atoms. Since the HSC-extension requires less operations than RGF, these speed-ups will increase as the device gets larger. A MATLAB version of HSC-extension algorithm is available at \urlhttps://els.comotion.uw.edu/express_license_technologies/glessuw.

Acknowledgements.

The authors acknowledge the support by the National Science Foundation under Grant ECCS 1231927. S. R. Patil acknowledges the support by UGC-India Grant No. F.5-50/2014(IC). M. P. Anantram and J. Qi acknowledge the support from the National Science Foundation under Grant No. 102781 for the DNA part of the paper.

Appendix A Complexity Derivation of HSC-extension for 3D Cuboidal Structures

To analyze the runtime complexity of HSC-extension, we consider a cuboid device with $N_{x}\times N_{y}\times N_{z}$ grid points per direction.

First we discuss the case of a cubic mesh, i.e. $N_{x}=N_{y}=N_{z}=N$. The operation count for evaluating diagonal entries in $\mathbf{G}^{r}$ was discussed by Lin et al. (Lin2009ab, , section 2.5.3). According to their analysis, the operation count grows as $\mathcal{O}(N^{6})$. For the HSC-extension, the complexity for evaluating diagonal entries of $\mathbf{G}^{<}$ is identical to the complexity for $\mathbf{G}^{r}$ (as discussed in Hetmaniuk et al. hetmaniuk2013nested for two-dimensional devices). The overall operation count will grow as $\mathcal{O}(N^{6})$.

Next we consider the case of an elongated device, where the numbers of grid points per direction satisfy $N_{x}=N_{z}=N\ll N_{y}$. The multilevel nested-dissection will identify $N_{y}/N_{z}$ subdomains, each discretized with $N\times N\times N$ grid points. When evaluating diagonal entries in $\mathbf{G}^{r}$, the operation count for each cubic subdomains will grow as $\mathcal{O}(N^{6})$. The remaining operations will involve dense matrices for the separators of dimension $N_{x}\times N_{z}$. The algebraic operations for one separator will include dense matrix-matrix multiplications and dense matrix inversions, yielding an asymptotic cost $\mathcal{O}(N_{x}^{3}N_{z}^{3})=\mathcal{O}(N^{6})$. The number of separators is also $\mathcal{O}(N_{y}/N)$. The overall operation count will grow as $\mathcal{O}(N^{6})\mathcal{O}(N_{y}/N)=\mathcal{O}(N^{5}N_{y})$.

Finally we consider the case of a flattened device, where the number of grid points per direction satisfy $N_{z}\ll N_{x}=N_{y}=N$. As discussed by Lin et al. (Lin2009ab, , section 2.5.3), this configuration is similar to a two-dimensional problem with $N\times N$ grid points. The prefactor will depend on $N_{z}$. The HSC algorithm and our extension can proceed as if the device is two-dimensional by replacing scalar algebraic operations with block algebraic operations (each block being of dimension $N_{z}\times N_{z}$). These block operations will cost $\mathcal{O}(N_{z}^{3})$. So the overall operation count will grow as $\mathcal{O}(N_{z}^{3})\mathcal{O}(N^{3})=\mathcal{O}(N_{z}^{3}N^{3})$.