# ELEQTRONeX: A GPU-Accelerated Exascale Framework for Non-Equilibrium Quantum
Transport in Nanomaterials

Saurabh S. Sawant<EMAIL_ADDRESS>François Léonard<EMAIL_ADDRESS>Zhi Yao<EMAIL_ADDRESS>Andrew Nonaka<EMAIL_ADDRESS>
###### Abstract

Non-equilibrium electronic quantum transport is crucial for the operation of
existing and envisioned electronic, optoelectronic, and spintronic devices.
The ultimate goal of encompassing atomistic to mesoscopic length scales in the
same nonequilibrium device simulation approach has traditionally been
challenging due to the computational cost of high-fidelity coupled
multiphysics and multiscale requirements. In this work, we present ELEQTRONeX
(ELEctrostatic Quantum TRansport modeling Of Nanomaterials at eXascale), a
massively-parallel GPU-accelerated framework for self-consistently solving the
nonequilibrium Green’s function formalism and electrostatics in complex device
geometries. By customizing algorithms for GPU multithreading, we achieve
orders of magnitude improvement in computational time, and excellent scaling
on up to 512 GPUs and billions of spatial grid cells. We validate our code by
computing band structures, current-voltage characteristics, conductance, and
drain-induced barrier lowering for various 3D configurations of carbon
nanotube field-effect transistors. We also demonstrate that ELEQTRONeX is
suitable for complex device/material geometries where periodic approaches are
not feasible, such as modeling of arrays of misaligned carbon nanotubes
requiring fully 3D simulations.

###### keywords:

NEGF , Broyden’s algorithm , CNTFET , GPU , AMReX

††journal: arXiv

[figure]style=plain,subcapbesideposition=top

[label1]organization=The Center for Computational Sciences and Engineering,
Lawrence Berkeley National Laboratory,addressline=1 Cyclotron Road,
city=Berkeley, postcode=94720, state=California, country=United States

[label2]organization=Sandia National Laboratories,addressline=7011 East Ave,
city=Livermore, postcode=94551, state=California, country=United States

##  1 Introduction

Non-equilibrium electronic transport determines the properties of many modern
electronic devices. As device dimensions are reduced and new nanomaterials
introduced, novel quantum phenomena are being explored for improved
performance. In parallel, the inherent three-dimensional (3D) nature of nano-
and meso-scale devices, including their constituent nanomaterials, requires a
device modeling approach that can not only capture quantum phenomena but do so
in a complex 3D geometry. In particular, such devices may involve multiple
densely packed materials/channels, necessitating modeling the cross-talk
caused by neighboring channels to ensure accurate performance evaluation.
Thus, there is a need for an efficient approach capable of simulating
multiple-channel devices on the order of hundreds of nanometer in length,
while incorporating all of the atoms of the constituent materials, and
capturing refined electrostatics over three orders of magnitude in size scale.

The Non-Equilibrium Green’s Function (NEGF) formalism has emerged as a
promising approach for device simulations. It has been applied to the study of
devices based on carbon nanotubes[1, 2], graphene [3], two-dimensional
materials [4], nanowires [5], silicon nanosheets [6], and single molecules
[7]. In addition it can describe time-dependent transport [8] as well as
response under external stimuli [9, 10]. However, despite its versatility,
NEGF simulations incur significant computational costs, primarily due to the
need for self-consistently solving NEGF equations with electrostatics in
complex 3D geometries. While traditionally implemented through parallel
deployment on CPUs, these simulations may still require hours to converge for
a single device operating point [6]. Furthermore, this computational demand is
further amplified when modeling multiple channels, as it necessitates
simultaneous self-consistent computations of charge density across all
channels. Although numerous parallel NEGF implementations exist for CPUs,
their GPU counterparts are notably scarce, with none currently capable of
modeling multiple channels.

Sophisticated NEGF implementations include Transiesta [11, 12], QuantumATK
[13], NEMO5 [14, 15], the latter one is extended to GPUs using libraries like
MAGMA and cuSPARSE [16], while other solvers use message passing interface
(MPI) for parallelizations. All of these implementations use independent
computational units for parallel distribution such as energy points used for
contour integration, discrete wavevectors, and voltage bias points. However,
this strategy suffers from scalability and memory limitations since each rank
must compute the entire matrices and require storage of the entire non-
distributed matrices on the order of tens of gigabytes in the processor
memory. These requirements can be alleviated by at most an order of magnitude
using OpenMP threads per MPI rank [11, 16]. However, for a more substantial
improvement and extension to modeling of multiple channels, a low-level
optimization approach is needed that allows for multiple MPI ranks and GPUs to
take part in computation of these matrices, say, at each energy point.

Initial attempts at low-level optimization for NEGF were made in Ref. [17]
with a GPU-based implementation of the recursive Green’s function (RGF)
algorithm. However, this strategy is limited to decomposing the computations
into two MPI ranks, each with its own GPU. Others [18] are leading efforts to
port the libNEGF solver to GPUs using OpenACC. Currently, the implementation
supports parallelization across energy and wavevector points, while ongoing
investigations focus on domain decomposition. These efforts have not yet
addressed a parallelization strategy for self-consistently coupling with
electrostatics, let alone extension to multiple channels.

We introduce ELEQTRONeX (ELEctrostatic Quantum TRansport modeling Of
Nanomaterials at eXascale), an open-source self-consistent NEGF implementation
built on the DOE Exascale Computing Project AMReX library [19, 20],
emphasizing low-level MPI/GPU parallelization strategies across key
components: electrostatics, NEGF, and the self-consistency algorithm. The
electrostatics module utilizes AMReX’s GPU-accelerated multigrid capabilities
and can handle intricate shapes of terminal leads represented as embedded
boundaries. For NEGF parallelization, large matrices such as Green’s and
spectral functions are decomposed on distributed across MPI ranks, with GPU
acceleration applied to kernel computations. We employ Broyden’s modified
second method for self-consistency and present an MPI/GPU parallelization
approach suitable for modeling multiple channel materials, simultaneously.
Overall we demonstrate excellent scaling up to 512 GPUs.

In this work, we demonstrate the capability of ELEQTRONeX by modeling
nanodevices with arrays of up to 20 carbon nanotubes of 10 and 100 nm channel
lengths with varied average spacing. Our focus is on investigating the effect
of angular misalignment and non-uniform spacing. These non-idealities closely
reflect experimental device configurations [21, 22] and result in non-periodic
configurations, requiring fully 3D simulations. By rigorously comparing
several key device performance metrics, such as ON and OFF currents and
subthreshold swing, against perfectly parallel configurations, we observe for
the first time that carbon nanotube field effect transistors (CNTFETs) are
robust against these non-idealities. Additionally, we highlight that the
modular design of the code allows for straightforward extension to model a
wide range of materials beyond nanotubes.

The rest of this paper is organized as follows: §2 describes our computational
approach, highlighting MPI/GPU optimization strategies for each sub-module
including electrostatics (§2.1), NEGF (§2.2), and self-consistency (§2.3).
§3.1 provides a rigorous validation of ELEQTRONeX for gate-all-around and
planar CNTFET configurations. §3.2 demonstrates the parallel performance of
the code for self-consistent simulations of gate-all-around CNTFET
configurations using up to 512 GPUs and 400 nm channel length. §3.3 discusses
the current-voltage characteristics from the nonequilibrium simulations run
for the parallel performance cases and highlights the challenges in simulating
$O(100$ nm) nanotubes, self-consistently. In §3.4, we study the effect of
modeling multiple CNTs in a fully 3D planar CNTFET configuration, focusing on
the impact of CNT misalignment compared to perfectly parallel alignment,
simulated with and without periodic boundary conditions.

##  2 Computational Approach

[] []

Figure 1:  (a) System setup for the NEGF method with multiple material
channels, denoted by superscripts $a$ and $b$. The structure of the main
matrix for computations is also shown. (b) Workflow for computing self-
consistent quantum transport through nanomaterials.

In the application of the NEGF method for modeling quantum transport through
nanodevices, the system is decomposed into material channels connected to
semi-infinite external leads such as source and drain contacts, as described
in Figure 1. Here, Hamiltonian $\boldsymbol{H}$ describes the electronic
structure of each material channel, self-energies
$\boldsymbol{\Sigma}_{l}^{R}$ for each lead $l$, represent its coupling to
external leads, and $V_{gs}$ and $V_{ds}$ represent the gate-source and drain-
source voltages, respectively.

To model the complete system, the ELEQTRONeX framework comprises three major
components, as shown in Figure 1: the electrostatic module, the quantum
transport module, and the self-consistent coupling between the two modules.
The electrostatic module, described in §2.1, computes the electrostatic
potential due to induced charges in the nanomaterial, as well as by device
terminals such as source, drain, and gate, which can be modeled as embedded
boundaries with intricate shapes. Using the electrostatic potential on each
site as input, the quantum transport module uses the NEGF method to compute
retarded Green’s function $\boldsymbol{G}^{R}(\boldsymbol{r},E)$ for each
material channel, as described in §2.2. This function can then be used to
compute induced charges at the $N$ total sites in the material, also referred
to as the system size. In this work, we model multiple carbon nanotubes using
the mode-space approximation [23], where for each nanotube, $N$ corresponds to
the total number of carbon rings along the length of the nanotube. The induced
charge computed by the NEGF method, in turn, affects the electrostatic
potential, and therefore, the two must be solved iteratively until self-
consistency is achieved. This is achieved using Broyden’s modified second
algorithm, as described in §2.3.

In what follows, we provide a description of the major distinctive features of
these modules, focusing on MPI and GPU parallelization. Each component
description also includes a discussion of the algorithmic complexity and
amenability to parallelization, with performance tests detailed in §3.2. In
our scaling tests we increase the nanotube length, making the number of carbon
rings $N$ an appropriate parameterization of system size, as the number of
grid cells and total atoms are proportional to $N$. We describe our
algorithmic components in terms of serial time complexity, as well as our
parallel implementation time complexity assuming the number of computational
resources increases in proportion to $N$.

All simulations in this work were performed using AMReX [24] (git hash $\leq$
ae3af4339) and ELEQTRONeX [25] (git hash $\leq$ 850e464).

###  2.1 Electrostastics, Potential Interpolation, and Charge Deposition

The underlying spatial discretization is a finite-volume mesh with cuboid grid
cells that store the electrostatic potential and permittivity. Coupling with
NEGF involves the use of static Lagrangian particles to model atoms/sites,
which interact with the underlying spatial mesh through interpolation and
deposition kernels.

The electrostatic module solves Poisson’s equation,

$\nabla\cdot\left[\epsilon(\vec{r})\nabla\phi(\vec{r})\right]=-e[\rho(\vec{r})-\rho_{0}(\vec{r})]$
(1)

where $\epsilon(\vec{r})$ is the permittivity, $\phi(\vec{r})$ is the
spatially varying electrostatic potential, $e$ is the electronic charge,
$\rho(\vec{r})$ is the charge density, and $\rho_{0}(\vec{r})$ is the
background charge density.

To solve equation (1) robustly, the electrostatic module uses the geometric
multigrid-based linear solvers implemented in AMReX [20]. The algorithm
employs a series of V- or W-cycles, which are iterative procedures involving
relaxation techniques, grid transfer operations, and coarse grid computations
[26]. Specifically, at each iteration, the algorithm transfers residual
information between grid levels, applies relaxation methods to reduce errors,
and performs computations at coarser grid levels [27]. Notably, at the
coarsest level of the grid hierarchy, the discretized equations are directly
solved using a strategy termed as the ‘bottom solver’. While the default
choice for this task is BiCGSTAB (Bi-Conjugate Gradient STABilized) [28],
users have the flexibility to select from a range of methods or integrate with
external libraries such as hypre (High-Performance Preconditioners) for
enhanced robustness [29]. To handle complex, 3D device configurations, the
solver capabilities include support for diverse domain boundary conditions as
well as embedded boundaries (EBs, also referred to as cut-cell
representations) to represent electric potential boundary conditions on the
surfaces of intricate device terminal shapes. Specifying a spatially-varying
function of voltages on domain or EBs allows the computation of complete
current-voltage characteristics of a device when coupled with NEGF.
Additionally, the solver supports heterogeneous permittivity distributions and
allows precise specification of atom locations within computational cells,
facilitated by electrostatic potential interpolation and charge deposition
schemes based on the cloud-in-cell algorithm [30]. In this algorithm, the
potential at each atomic location is obtained by tri-linear interpolation of
the potential at eight neighboring cell-centered grid points. Similarly, the
algorithm employs a tri-linear deposition scheme to distribute charge density
from an atom to the eight neighboring grid points. Since we use the mode-space
approximation for modeling carbon nanotubes, we average the potential computed
over all atoms around the ring before using it in the NEGF. Similarly, we
distribute the charge computed at each ring equally over all atoms around the
ring. Furthermore, this framework offers flexible domain decomposition
strategies for parallelization, allowing the computational domain to be
decomposed irrespective of EBs, regions of varying permittivity, and atomic
structure of a material.

The serial time complexity of the linear solvers, as well as the charge
deposition and potential interpolation algorithms, are $O(N)$. These
algorithms are highly parallelizable and scaling is bounded by inter-processor
communication. In an ideal parallel implementation, the parallel time
complexity is $O(1)$.

###  2.2 Quantum transport using the NEGF method

####  2.2.1 Overview of the NEGF method

In this work we represent Hamiltonian $\boldsymbol{H}$ of size $(N\times N)$
using a tight-binding model, but the scheme is amenable to any representation.
The Hamiltonian has block tri-diagonal form, with the diagonal entries of
$\boldsymbol{H}$ corresponding to the electrostatic potential energy $-e\phi$
at each site within the material; in the present case, at each carbon ring.

Central to computing the DC transport properties of the device is the
calculation of the retarded Green’s function, expressed as

$\begin{split}\boldsymbol{G}^{R}(\vec{r},E)&=\left[(E+i\eta)\boldsymbol{I}-\boldsymbol{H}(\vec{r})-\sum_{l}^{L}{\boldsymbol{\Sigma}_{l}^{R}}\right]^{-1}\\\
\boldsymbol{\Sigma}_{l}^{R}(E)&=\boldsymbol{\tau}_{l}^{\dagger}\boldsymbol{g}_{l}^{R}\boldsymbol{\tau}_{l}\\\
\end{split}$ (2)

where $E$ denotes the electron energy, $\eta$ is an infinitesimal constant,
$\boldsymbol{I}$ is the identity matrix, $L$ denotes the number of leads in
the system, and $\boldsymbol{\tau}$ and $\boldsymbol{g}^{R}$ are matrices for
the device-lead coupling and the surface Green’s function, respectively. From
equation (2) and Figure 1, we note that
$\alpha_{j}=(E+i\eta)+H_{(j,j)}-\sum_{l}^{L}\Sigma_{l,(j,j)}$ for $j=0$ to
$N-1$, and $\beta_{j}=H_{(j,j+1)}$, $\gamma_{j}=H_{(j+1,j)}$ for $j=0$ to
$N-2$. For a system with 2 semi-infinite leads, only $\Sigma_{1,(0,0)}$ and
$\Sigma_{2,(N-1,N-1)}$ are non-zero.

Using the retarded Green’s function, we can compute charge density matrix
$\boldsymbol{\rho}$ as

$\begin{split}\boldsymbol{\rho}(\vec{r})&=-\frac{g_{s}g_{b}}{2\pi
i}\int_{-\infty}^{\infty}\boldsymbol{G}^{<}(\vec{r},E)dE\\\ \end{split}$ (3)

where $g_{s}=2$ and $g_{b}$ are spin and band degeneracies,
$\boldsymbol{G}^{<}=i\sum_{l}^{L}\boldsymbol{A}_{l}F_{l}$ represents the
lesser Green’s function,
$\boldsymbol{A}_{l}=\boldsymbol{G}^{R}\boldsymbol{\Gamma}_{l}\boldsymbol{G}^{R\dagger}$
denotes the lead-specific spectral function matrix, and
$\boldsymbol{\Gamma}_{l}=i(\boldsymbol{\Sigma}_{l}^{R}-\boldsymbol{\Sigma}_{l}^{R\dagger})$
stands for the broadening matrix. $F_{l}\equiv F_{l}(E-\mu_{l},k_{b}T_{l})$
represents the lead-specific Fermi function, depending on the electrochemical
potential $\mu_{l}$ and temperature $T_{l}$, where $k_{b}$ is the Boltzmann
constant. Using the above identities and following the approach described in
Refs. [31, 32, 11], We calculate $\boldsymbol{\rho}$ in two parts as

$\begin{split}\boldsymbol{\rho}(\vec{r})&=\frac{1}{\pi}\Im\left[\int_{-\infty}^{E_{min}}\boldsymbol{G}^{R}(\vec{r},E)F_{min}\;dE\right]-\frac{1}{2\pi}\int_{E_{min}}^{E_{max}}\sum_{l}\boldsymbol{A}_{l}(\vec{r},E)F_{l}\;dE\\\
\end{split}$ (4)

where $F_{min}$ is the lowest Fermi level,
$E_{min}=\mu_{min}-f_{1}k_{b}T_{min}$, $E_{max}=\mu_{max}+f_{2}k_{b}T_{max}$,
($\mu_{min}$, $T_{min}$) and ($\mu_{max}$, $T_{max}$) are minimum and maximum
electrochemical potentials and temperatures, respectively. The factors $f_{1}$
and $f_{2}$ are problem-dependent and typically set to $14$ to capture the
tail of the Fermi levels in the integration. We efficiently calculate the
first part of the integral using contour integration with the residue theorem,
and use Gauss-Legendre mapping for accuracy with fewer integration points
[11]. At equilibrium, i.e. when all leads have equal $\mu$ and $T$, the charge
density can be calculated solely using the residue theorem by setting
$E_{l}=E_{u}=\mu+fk_{b}T$. Similarly, the neutral charge density matrix
$\boldsymbol{\rho}_{0}(\vec{r})$ is computed with $E_{l}=E_{u}=0$. The
diagonal entries of the density matrix, $\boldsymbol{\rho}_{mm}$ and
$\boldsymbol{\rho}_{0,mm}$, represent the charge at each material site, which
are spread on the computational grid using the cloud-in-cell algorithm before
recomputing potential using equation (1).

Equations (1), (2), and (4) are evaluated until self-consistency is achieved,
as described in §2.3. Subsequently, essential device properties such as
current and transmission can be determined as

$\begin{split}I_{l}&=\frac{g_{s}e}{h}\int_{E_{min}}^{E_{max}}\text{Tr}\left[F_{l}\boldsymbol{\Gamma}_{l}\boldsymbol{A}+i\boldsymbol{\Gamma}_{l}\boldsymbol{G}^{<}\right]dE\\\
T_{pq}(E)&=\text{Tr}\left[\boldsymbol{\Gamma}_{p}\boldsymbol{G}^{R}\boldsymbol{\Gamma}_{q}\boldsymbol{G}^{R\dagger}\right]\\\
\end{split}$ (5)

where $I_{l}$ is the steady-state current through lead $l$, $h$ is Plank’s
constant, $\boldsymbol{A}=\sum_{l}^{L}\boldsymbol{A}_{l}$ is the total
spectral function, and $T_{pq}$ is the transmission from lead $p$ to $q$. At
zero source-drain bias, we may compute conductance as

$\begin{split}G(E)&=G_{0}\int_{-\infty}^{\infty}T(E)\left(-\frac{\partial
F}{\partial E}\right)dE\\\ \end{split}$ (6)

where $T$ is the transmission at equilibrium and $G_{0}=g_{s}e^{2}/h$ is the
quantum of conductance.

####  2.2.2 Flexible data structure to accommodate different materials

The composition of the matrices described above depends on the materials under
consideration, the device structure, and the chosen representation for the
Hamiltonian. Each block or element of these matrices can be a single number, a
submatrix, or an array storing entries of a diagonal submatrix. For example,
when modeling a CNT under the mode-space approximation with a single subband,
each block is represented as a single number. However, when considering
multiple subbands, it is represented as an array, with each entry
corresponding to a subband.

To ensure flexibility in data structure and facilitate future extensions to
other materials, we employed various C++ object-oriented programming
techniques, including templates, virtual functions, and operator overloading.
The general procedure of the NEGF method is implemented as a templated base
class, with specialized classes derived from it for specific materials and
their unique data structures. The general code in the base class, such as the
one used to compute equation (2), works regardless of the underlying data
structure of the materials-specific matrices. This is achieved by overloading
mathematical operators, such as matrix multiplication, for the specific data
structures of a material. Furthermore, during the execution of specific steps
in the NEGF algorithm, the code searches for material-specific specializations
of that step. If none are found, it defaults to the general procedure. For
instance, the surface Green’s function $g^{R}$ can be computed using a
decimation technique [33, 34], but may be overridden, for example when an
analytical expression is available (e.g. for CNTs [23]). This approach
facilitates easy extension to new materials, as only the aspects specific to
the new material need to be coded.

####  2.2.3 Overview of parallelization strategy

In the NEGF approach, the two core computations are the matrix inversion to
compute the retarded Green’s function $\boldsymbol{G}^{R}$ and the calculation
of the spectral function $\boldsymbol{A}$. These computations are performed
multiple times, for each integration point on all contour paths used to
compute density matrix (see equation (4)). Typically this means
$\mathcal{O}(100-10k)$ of these computations per iteration of the self-
consistency algorithm, depending on the problem. Existing parallelization
approaches map each energy point to a MPI process, leading to an efficient
parallel scenario with no communication overhead for these tasks. However,
this approach has limitations. Firstly, it restricts parallelism to the number
of energy points, limiting scalability. Secondly, it imposes substantial
memory demands, as each MPI process must store entire matrices for
$\boldsymbol{G}^{R}$ and $\boldsymbol{A}$, alongside additional memory for
other quantities. This constraint can significantly restrict the utilization
of available MPI processes, even on supercomputers. For instance, on the
Perlmutter supercomputer [35], CPU-only nodes feature 64 cores with 512 GB
DDR4 DRAM memory, allowing for the storage of only two complex floating-point
matrices of size (16k × 16k). For complex 3D configurations, we require a
scalable solver capable of managing matrices substantially larger in size—by
one or two orders of magnitude—while also allocating sufficient memory for
computational grids related to electrostatics on each node. Additionally, any
future extension of the code to model time-dependent equations will further
increase the requirements for the number of matrices stored [8].

To address these challenges and enable efficient computations on MPI and GPUs,
we adopt a novel parallelization strategy. In our approach, each MPI process
is assigned the responsibility of computing a nearly equal predetermined
number of columns of the Green’s or spectral function matrix. This
distribution ensures that all MPI processes contribute to computations at each
energy point, while significantly reducing memory requirements per process.
Consequently, each MPI process only needs to store Hamiltonian entries and
compute charges for the assigned range of columns. Moreover, this strategy
facilitates a high degree of parallelization for large systems, when the
number of columns ($N$) in the Hamiltonian matrix far exceeds the number of
energy points.

In extending this strategy to GPUs, each MPI process is bound to a GPU device,
with GPU threads invoked to handle the assigned columns, independently. The
GPU-accelerated portion of the code is designed with portability in mind,
ensuring compatibility across different GPU architectures and vendors. This
approach maximizes computational efficiency while accommodating diverse
hardware configurations.

####  2.2.4 Parallel computations of Green’s and spectral functions

To compute $\boldsymbol{G}^{R}$ and $\boldsymbol{A}$, we employ a block tri-
diagonal matrix inversion algorithm [36, 37], typically implemented as a
serial algorithm in existing NEGF numerical methods [11]. Here, we adopt this
algorithm for hybrid MPI/GPU architectures with one GPU per MPI rank, such
that individual columns of $\boldsymbol{G}^{R}$ and $\boldsymbol{A}$ are
computed independently by parallelizing across all invoked MPI ranks and GPU
threads, as illustrated in Algorithm 1 and Figure 2. For the computation of
$\boldsymbol{G}^{R}$, each GPU-thread invoked on each MPI-rank begins by
calculating the diagonal element of a column and subsequently computing the
off-diagonal elements of that column using a recursive approach. The spectral
function $\boldsymbol{A}$ is computed by aggregating the lead-specific
spectral functions $\boldsymbol{A}_{l}$, for each lead $l$, where each column
of $\boldsymbol{A}_{l}$ is computed recursively, similarly to the Green’s
function.

Input:
$\boldsymbol{\alpha},\boldsymbol{x},\boldsymbol{\tilde{x}},\boldsymbol{y},\boldsymbol{\tilde{y}}$
as computed from Algorithm 2.

Output: $\boldsymbol{G}^{R},\boldsymbol{A}$.

Define: Let $N_{p}$ be the portion of system size evaluated by MPI rank $p$.

Set: $\boldsymbol{G}^{R},\boldsymbol{A}$ to zero.

foreach  _MPI rank $p$_ do

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

1 $T\leftarrow t+pN_{p}$;

// global index, assuming $N_{p}$ is the same for all $p$

2
$\boldsymbol{G}_{(T,T)}^{R}\leftarrow[\boldsymbol{\alpha}_{(T)}-\boldsymbol{x}_{(T)}-\boldsymbol{y}_{(T)}]^{-1}$;

// diagonal element

3
$\boldsymbol{G}_{(n+1,T)}^{R}\leftarrow-\boldsymbol{\tilde{x}}_{(n)}\boldsymbol{G}_{(n,T)}^{R}\quad\text{for
}n\geq T$;

// elements below diagonal

4
$\boldsymbol{G}_{(n-1,T)}^{R}\leftarrow-\boldsymbol{\tilde{y}}_{(n)}\boldsymbol{G}_{(n,T)}^{R}\quad\text{for
}n\leq T$;

// elements above diagonal

foreach  _lead $l$ | $l\in[1,L]$_ do

5 set:  $k\leftarrow\text{block index for lead }l$;

6
$\boldsymbol{A}_{l,(n,T)}\leftarrow\boldsymbol{G}_{(n,k)}^{R}\Gamma_{l,(k,k)}\boldsymbol{G}_{(T,k)}^{R\dagger}\quad\text{for
}n=k$;

7
$\boldsymbol{A}_{l,(n+1,T)}\leftarrow-\boldsymbol{\tilde{x}}_{(n)}\boldsymbol{A}_{l(n,T)}\quad\text{for
}n\geq T$;

8
$\boldsymbol{A}_{l,(n-1,T)}\leftarrow-\boldsymbol{\tilde{y}}_{(n)}\boldsymbol{A}_{l,(n,T)}\quad\text{for
}n\leq T$;

9 $\boldsymbol{A}\leftarrow\boldsymbol{A}+\boldsymbol{A}_{l}$;

// aggregate $\boldsymbol{A}_{l}$ into $\boldsymbol{A}$.

Algorithm 1 MPI/GPU parallelization of $\boldsymbol{G}^{R}$ and
$\boldsymbol{A}$ computations.

To initiate the recursion for $\boldsymbol{A}_{l}$, each GPU thread computes
$\boldsymbol{A}_{l,(n,T)}$, setting $n=k$, subsequently, computing the
remaining entries in a given column using a suitable recursion method, as
shown in Algorithm 1. Here, $T$ represents the global index depending on MPI
rank $p$ and GPU thread index $t$, and $k$ corresponds to the block index for
lead $l$. To begin the recursion, each GPU thread needs to store some lead-
specific quantities and perform some additional computations. For example, in
a system with two semi-infinite leads having block indices $k=0$ and $N-1$,
respectively, the second GPU thread ($t=1$) on MPI rank $p=0$ initially
computes entries
$\boldsymbol{A}_{1,(0,1)}=\boldsymbol{G}_{(0,0)}^{R}\boldsymbol{\Gamma}_{1,(0,0)}\boldsymbol{G}_{(1,0)}^{R\dagger}$
and
$\boldsymbol{A}_{2,(N-1,1)}=\boldsymbol{G}_{(N-1,N-1)}^{R}\boldsymbol{\Gamma}_{2,(N-1,N-1)}\boldsymbol{G}_{(1,N-1)}^{R\dagger}$.
To achieve this, the GPU thread needs lead specific quantities,
$\Gamma_{1,(0,0)}$ and $\Gamma_{2,(N-1,N-1)}$, and entries
$\boldsymbol{G}_{(0,0)}$, $\boldsymbol{G}_{(N-1,N-1)}$, which can in turn be
used to compute, $\boldsymbol{G}_{(1,0)}$ and $\boldsymbol{G}_{(1,N-1)}$ using
a suitable recursion.

In practice, for NEGF method with tight-binding approximation for Hamiltonian,
we only need the diagonal entries of $\boldsymbol{G}^{R}$ and
$\boldsymbol{A}$, and therefore, we offer a way to turn off storage of off-
diagonal elements of these matrices. However, computing the diagonal elements
of $\boldsymbol{A}_{l}$ still requires elements $\boldsymbol{G}_{(T,k)}^{R}$,
necessitating $O(N)$ recursive computations per thread.

[] [] []

Figure 2: (a) Illustration of the MPI/GPU parallelized part of the algorithm
for computation of $\boldsymbol{G}^{R}$ and $\boldsymbol{A}$. (b) Optimization
of recursive part of the algorithm involving overlap of CPU computations and
CPU-to-GPU asynchronous copy. (c) Illustration of MPI/GPU parallelization
strategy for steps 14 to 16 in Algorithm 3 to compute
$\boldsymbol{s}_{p}=\boldsymbol{W}_{p}\boldsymbol{V}_{p}^{T}\boldsymbol{f}_{p}$,
assuming iterator $m$=3, total MPI ranks 2, and $N_{0}=N_{1}=2$. Solid and
dashed edges of matrices and vectors distinguish elements stored and processed
by MPI ranks 0 and 1, respectively, while the elements read or written by GPU
threads $0$, $1$, or both $0$ and $1$, of MPI rank 0, are represented by a
circle (o), cross (x), or by a circle-cross overlap, respectively.

From Figure 2, we see that to perform aforementioned recursions in parallel,
each MPI-rank and GPU-thread needs access to vectors $\boldsymbol{\tilde{x}}$
and $\boldsymbol{\tilde{y}}$. These vectors can be computed using a recursion
shown in Algorithm 2 using vectors $\boldsymbol{\alpha}$,
$\boldsymbol{\beta}$, $\boldsymbol{\gamma}$, defined earlier in § 2.2.1 and
Figure 1, representing the diagonal, upper and lower adjacent diagonals of the
block-tridiagonal matrix $(G^{R})^{-1}$ that is being inverted.

Input: Block vectors
$\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}$ as defined in
Figure 1.

Output: Block vectors
$\boldsymbol{x},\boldsymbol{\tilde{x}},\boldsymbol{y},\boldsymbol{\tilde{y}}$
of size $N$.

Set: $\boldsymbol{y}_{(0)}=\boldsymbol{x}_{(N-1)}=0$.

foreach  _vector index $n$ | $n\in[0,N)$_ do

$\boldsymbol{y}_{(n)}\leftarrow\boldsymbol{\beta}_{(n-1)}\boldsymbol{\tilde{y}}_{(n)}$;

$\boldsymbol{\tilde{y}}_{(n)}\leftarrow[\boldsymbol{\alpha}_{(n-1)}-\boldsymbol{y}_{(n-1)}]^{-1}\boldsymbol{\gamma}_{(n-1)}$;

$\boldsymbol{x}_{(n)}\leftarrow\boldsymbol{\gamma}_{(n)}\boldsymbol{\tilde{x}}_{(n)}$;

$\boldsymbol{\tilde{x}}_{(n)}\leftarrow[\boldsymbol{\alpha}_{(n+1)}-\boldsymbol{x}_{(n+1)}]^{-1}\boldsymbol{\beta}_{(n)}$;

Algorithm 2 Recursive computations of vectors needed as input for Algorithm 1.

The recursive computation of $\boldsymbol{\tilde{x}}$ and
$\boldsymbol{\tilde{y}}$ poses a bottleneck for this algorithm, as explained
shortly, limiting its time complexity to be proportional to $O(N)$, regardless
of number of MPI ranks used. This limitation can be partially mitigated by
implementing long recursion in smaller segments and overlapping parts of
recursion computations with CPU-to-GPU asynchronous copies, as illustrated in
Figure 2. For instance, assuming a system size $N$ simulated with $P$ total
MPI ranks, doubling the system size to $2N$ would typically result in double
of the algorithm time, even if computational resources are doubled to $2P$ MPI
ranks. Yet, by overlapping the computation of the latter half of the array
with the CPU-to-GPU copy of the former half, the algorithm time can be reduced
to a factor of approximately $1.5$. Scaling this approach to simulate system
size $ZN$ leads to a time-increase factor of $(Z+1)/2$, instead of $Z$. In
practice, the exact value of this factor depends on the size of the recursive
array segment chosen for overlap optimization. In general, if the segment size
is such that the time for computation equals the CPU-to-GPU copy time, we may
approach the ideal limit.

The time complexity of the serial block-tridiagonal algorithm is $O(N^{2})$,
whereas the parallel time complexity can be expressed as,

$T(N)=O\left(\frac{N}{s_{1}}\right)+O\left(\frac{N^{2}}{s_{2}P}\right).$ (7)

Here, the first and second terms represent time complexities of the serial
CPU-based recursion and GPU-parallelizable parts of the algorithm,
respectively, $s_{1}\in[1,2)$ represents the speedup obtained from the overlap
optimization discussed above, $s_{2}$ represents the speedup of GPUs over CPUs
in evaluating a column of length $N$ (typically, $s_{2}\sim O(100)$), and $P$
indicates the total number of parallel processing units. In our GPU
implementation we invoke as many GPU threads as the number of columns per
processor, thus $P$ is typically equal to $N$ and the resulting parallel time
complexity is $T(N)=O(N/s_{1})+O(N/s_{2})$. Since $s_{1}<s_{2}$, we expect the
algorithm to scale as $O(N/s_{1})$, and $s_{1}$ to approach 2 in the best case
scenario.

In the future, we plan to pre-compute and store the vectors
$\boldsymbol{\tilde{x}}$ and $\boldsymbol{\tilde{y}}$ for each energy point
over which $\boldsymbol{G}^{R}$ and $\boldsymbol{A}$ are computed. This way
the linear cost of recursion can be further amortized over the number of
independent energy points, leading to potentially significant savings in time,
albeit at the cost of additional memory requirement on the order of $O(NE_{\rm
pts})$, where $E_{\rm pts}$ is the number of energy points.

###  2.3 Parallel self-consistency module using Broyden’s second algorithm

The self-consistent calculation involves finding the nearest zero of the
function $\boldsymbol{f}=\boldsymbol{\rho}^{in}-\boldsymbol{\rho}^{out}$ where
the charge density reaches a self-consistent value. The superscripts ‘$in$’
and ‘$out$’ denote the input charge densities used for computing the
electrostatic potential and those generated from the NEGF algorithm,
respectively. Note that for this algorithm we only need the diagonal entries
of the density matrix. To achieve faster convergence, we employ Broyden’s
modified second algorithm, which does not require the storage of the inverse
Jacobian matrix, thus significantly reducing memory usage. Since the algorithm
details are provided in various references [38, 39, 40, 41], we focus solely
on presenting the MPI/GPU parallelization of this algorithm here.

The pseudo-code for the MPI/GPU parallelization of this algorithm is outlined
in Algorithm 3. Adhering to a strategy consistent with the parallelization of
the NEGF algorithm, each MPI rank is assigned to compute a designated subset
of the charge density vector, denoted as $\boldsymbol{\rho}_{p}$, with a size
of $N_{p}$, where the subscript $p$ denotes quantities specific to MPI rank
$p$. Step 1 of the algorithm involves using $\boldsymbol{\rho}_{p}^{in}$ and
$\boldsymbol{\rho}_{p}^{out}$ as input to electrostatics and output from the
NEGF algorithm for MPI rank $1$, respectively. Steps 2 to 4 entail calculating
local vectors, which can be easily parallelized across GPU threads. Steps 5
and 6 involve reduction operations, where all threads attempt to modify a
single variable, resulting in a read-modify-write race condition [42]. To
address this, GPU-atomic operations are utilized to ensure that the read-
modify-write operation sequence on a memory location by one thread does not
overlap with that by the other thread. Additionally, more advanced reduction
algorithms such as warp-shuffle [43] may offer improved efficiency, although
they were not implemented in this work. Steps 7 and 8 involve reduction
operations across all MPI ranks, which can be achieved using the MPI_Allreduce
primitive.

Steps 9 to 16 are intermediate steps to compute the vector
$\boldsymbol{s}_{p}$, which is used to update the input charge density
$\boldsymbol{\rho}_{p}^{in}$ for the next iteration. Among these, steps 12,
13, 17, and 18 can be easily parallelized across GPU threads without requiring
any synchronization or communication across GPU threads or MPI ranks. The most
computationally intensive steps are steps 9 to 11 and steps 14 to 16, which
are very similar in nature. Here, we provide elaboration on steps 14 to 16,
which entail computation of
$\boldsymbol{s}_{p}=\boldsymbol{W}_{p}\boldsymbol{V}_{p}^{T}\boldsymbol{f}_{p}$,
using the illustration provided in Figure 2.

We begin with computing the matrix-vector product,
$\boldsymbol{k}_{p}=\boldsymbol{V}_{p}^{T}\boldsymbol{f}_{p}$, where each
thread is assigned to process one column of $\boldsymbol{V}_{p}^{T}$. However,
this step introduces a read-modify-write race condition, when multiple threads
(in this case, 2) attempt to modify the same elements of $\boldsymbol{k}_{0}$
concurrently, similar to steps 5 and 6. To address this issue, we employ GPU-
atomic operations to ensure proper synchronization. The subsequent step
involves aggregating all local copies of $\boldsymbol{k}_{p}$, into a global
copy $\boldsymbol{k}$ on each MPI rank using the MPI_Allreduce primitive.
Finally, in the third step, we evaluate
$\boldsymbol{s}_{p}=\boldsymbol{W}_{p}\boldsymbol{k}$, where each thread is
responsible for computing one row of $\boldsymbol{W}_{p}$ multiplied by
$\boldsymbol{k}$. The computational time complexity of the serial algorithm is
$O(m_{\rm avg}N)$, where $m_{\rm avg}$ represents the average number of
iterations required for convergence for a system size of $N$. In an ideal
implementation the parallel time complexity is $O(m_{\rm avg})$, noting that
as the system size increases $m_{\rm avg}$ increases sub-linearly in our
demonstration cases.

This algorithm can be easily extended to accommodate simulations involving
multiple channel materials without requiring modifications to the core
algorithm. In this scenario, $N_{p}$ represents the sum of portions of all
system sizes evaluated by MPI rank $p$, where $N_{p}$ is defined as
$\sum_{c=1}^{c=C}N_{p}^{c}$, where $C$ represents number of material channels,
and $N_{p}^{c}$ represents portion of the system size of material channel $c$
evaluated by MPI rank $p$. For example, $\boldsymbol{\rho}_{p}^{in}$,
represents a vector composed of all charge density vectors
$\boldsymbol{\rho}_{p}^{c}$ for $c\in[1,C]$, denoted as
$\boldsymbol{\rho}_{p}^{in}\equiv(\boldsymbol{\rho}_{p}^{in,1},\boldsymbol{\rho}_{p}^{in,2},...,\boldsymbol{\rho}_{p}^{in,C})$.
This array can be formed using the MPI_Allgatherv primitive before proceeding
to step $2$. The subsequent steps ($2$ to $18$) are executed as usual. Before
repeating step $1$ for the next iteration, we decompose the vector
$\boldsymbol{\rho}_{p}^{in}$ into channel-specific charge densities
$\boldsymbol{\rho}_{p}^{in,c}$ for all $c\in[1,C]$. This strategy ensures that
during minimization of the function
$\boldsymbol{f}=\boldsymbol{\rho}^{in}-\boldsymbol{\rho}^{out}$ we account for
contribution from all material channels.

Input: A guess for charge density $\boldsymbol{\rho}^{in}$, Broyden fraction
$\zeta$, maximum number of iterations, $M$.

Output: Self-consistent charge
$\boldsymbol{\rho}^{out}\approx\boldsymbol{\rho}^{in}$.

Define: $\boldsymbol{V}_{p}^{T}$ and $\boldsymbol{W}_{p}$ as matrices of size
$M\times N_{p}$ and $N_{p}\times M$, respectively, $\vec{k}_{p}$ as vector of
size $M$, and all other quantities denoted by lowercase bold letters as
vectors of size $N_{p}$, where $N_{p}$ is the portion of system size evaluated
by MPI rank $p$.

Set: iterator, $m\leftarrow 0$, maximum relative norm,
$\epsilon_{max}\leftarrow\infty$, and all matrices, vectors to $0$.

while  _$\epsilon_{max} <10^{-5}$ and $m<M$ _ do

foreach  _MPI rank $p$_ do

1
$\boldsymbol{\rho}_{p}^{in}\rightarrow\text{Electrostatics}\rightarrow\text{NEGF}\rightarrow\boldsymbol{\rho}_{p}^{out}$;

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

2
$\boldsymbol{f}_{p,(t)}\leftarrow\boldsymbol{\rho}_{p,(t)}^{in}-\boldsymbol{\rho}_{p,(t)}^{out}$;

// compute minimization function

3 $\boldsymbol{\Delta
f}_{p,(t)}\leftarrow\boldsymbol{f}_{p,(t)}-\boldsymbol{f}_{p,(t)}^{\text{old}}$;

// compute change in minimization function

4
$\boldsymbol{\epsilon}_{p,(t)}\leftarrow\left|\boldsymbol{f}_{p,(t)}/(\boldsymbol{\rho}_{p,(t)}^{in}+\boldsymbol{\rho}_{p,(t)}^{out})\right|$;

// compute relative norm

5
$\epsilon_{max,p}\overset{\text{atomic}}{\leftarrow}\max(\boldsymbol{\epsilon}_{p,(t)})$;

// assign maximum element of $\boldsymbol{\epsilon}_{p}$ to $\epsilon_{max,p}$

6 $d_{p}\overset{\text{atomic}}{\leftarrow}||\boldsymbol{\Delta
f}_{p,(t)}||^{2}$;

// aggregate square of each element of $\boldsymbol{\Delta f}_{p}$ into
$d_{p}$

7MPI_Allreduce($\epsilon_{max}\leftarrow\epsilon_{max,p}$, max);

// find maximum $\epsilon_{max,p}$ over all processes

8 MPI_Allreduce($d\leftarrow d_{p}$, sum);

// aggregate sum of $d_{p}$ over all processes

if  _m > 0_ then

if  _m > 1_ then

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

foreach  _index $j$ | $j\in[1,m-1]$_ do

9
$\boldsymbol{k}_{p,(j)}\overset{\text{atomic}}{\leftarrow}\boldsymbol{V}_{p,(j,t)}^{T}\boldsymbol{\Delta
f}_{p,(t)}$;

// read-modify-write race condition

10MPI_Allreduce($\boldsymbol{k}\leftarrow\boldsymbol{k}_{p}$, sum);

// aggregate sum of $\boldsymbol{k}_{p}$ over all processes

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

foreach  _index $j$ | $j\in[1,m-1]$_ do

11
$\boldsymbol{s}_{p,(t)}\leftarrow\boldsymbol{W}_{p,(t,j)}\boldsymbol{k}_{(j)}$;

// each thread reads $\boldsymbol{k}$

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

12 $\boldsymbol{V}_{p,(m,t)}^{T}\leftarrow d^{-1}\boldsymbol{\Delta
f}_{p,(t)}$;

// Evaluate $m^{th}$ row of $\boldsymbol{V}_{p}^{T}$

13 $\boldsymbol{W}_{p,(t,m)}\leftarrow-\eta\boldsymbol{\Delta
f}_{p,(t)}+\boldsymbol{\rho}_{p,(t)}^{in}-\boldsymbol{\rho}_{p,(t)}^{in,old}-\boldsymbol{s}_{p,(t)}$;

// Evaluate $m^{th}$ column

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

foreach  _index $j$ | $j\in[1,m]$_ do

14
$\boldsymbol{k}_{p,(j)}\overset{\text{atomic}}{\leftarrow}\boldsymbol{V}_{p,(j,t)}^{T}\boldsymbol{f}_{p,(t)}$;

// read-modify-write race condition

15MPI_Allreduce($\boldsymbol{k}\leftarrow\boldsymbol{k}_{p}$, sum);

// aggregate sum of $\boldsymbol{k}_{p}$ over all processes

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

foreach  _index $j$ | $j\in[1,m]$_ do

16
$\boldsymbol{s}_{p,(t)}\leftarrow\boldsymbol{W}_{p,(t,j)}\boldsymbol{k}_{(j)}$;

// each thread reads $\boldsymbol{k}$

foreach  _GPU thread $t$ | $t\in[0,N_{p})$_ do

17
$\boldsymbol{\rho}_{p,(t)}^{in,old}\leftarrow\boldsymbol{\rho}_{p,(t)}^{in}$;

// update $\boldsymbol{n}_{p,(t)}^{in,old}$ for next iteration

18
$\boldsymbol{\rho}_{p,(t)}^{in}\leftarrow\boldsymbol{\rho}_{p,(t)}^{in,old}-\zeta\boldsymbol{f}_{p,(t)}-\boldsymbol{s}_{p,(t)}$;

// predict $\boldsymbol{n}_{p,(t)}^{in}$ for next iteration

19$m\leftarrow m+1$;

// increment iterator $m$

Algorithm 3 MPI/GPU Parallelization of Broyden’s modified second algorithm.

##  3 Results

Here we provide validation of ELEQTRONeX for gate-all-around and planar CNTFET
configurations (§3.1). Then we demonstrate the parallel performance of the
code for self-consistent simulations of gate-all-around CNTFET configurations
(§3.2). Next, we discuss the current-voltage characteristics from these
nonequilibrium simulations and discuss challenges in modeling longer nanotubes
(§3.3). Finally, we study the effect of modeling multiple CNTs, both aligned
and non-aligned, in fully 3D planar CNTFET configurations compared to
corresponding periodic configurations (§3.4).

###  3.1 Validation Studies

We validate ELEQTRONeX in obtaining DC transport properties by comparing our
results with those reported in Refs. [1] and [2] for gate-all-around and
planar CNTFET configurations, respectively. Further details specific to
modeling CNTFETs can be found in the references. In brief, we consider a
(17,0) zigzag CNT of diameter 1.3 nm and nearest-neighbor tight-binding
coupling of 2.5 eV, giving a bandgap of 0.55 eV. The Hamiltonian is written in
terms of nearest-neighbor coupling between parallel carbon rings. The Fermi
level of the source and drain contacts is 1 eV below the CNT midgap before
self-consistency, as would be expected from Pd. The gate oxide has a
dielectric constant of 3.9 to simulate SiO2.

Figures 3, 3, and 3 demonstrate good agreement for the drain-source current
$I_{ds}$ versus drain-source voltage $V_{ds}$, $I_{ds}$ versus gate-source
voltage $V_{gs}$, and drain-induced barrier lowering (DIBL) versus channel
length $L$, respectively, for a gate-all-around CNTFET configuration. The
simulations agree despite methodological differences compared to Ref. [1],
including differences in solving Poisson’s equation (finite-difference,
successive overrelaxation method versus finite-volume, multigrid), charge
deposition to mesh (Gaussian spreading versus cloud-in-cell), calculation of
surface Green’s function (iterative layer doubling technique [33] versus
analytical [23]), and the type of Broyden’s algorithm employed (first versus
parallelized modified second method). Furthermore, in Figure 3, DIBL remains
flat for channel lengths exceeding 20 nm, consistent with the semi-empirical
fit reported by Ref. [1]. These additional runs are discussed in later in §3.3
on performance studies.

[] [] [] []

Figure 3:  Current-voltage characteristics of CNTFETs. (a) Current versus
source-drain voltage for a gate-all-around CNTFET with channel length of
$L=20$ nm. (b) Current versus gate-source voltage for the gate-all-around
CNTFET. (c) Drain-induced barrier lowering (DIBL) as a function of $L$. (d)
Conductance $G$ versus $V_{gs}$ for planar CNTFET configurations with $d=3.2$
nm. In all panels insets show valence band edge and schematics with dimensions
in nanometers. Comparison is made with the data in Refs [1] and [2].

Next, we validate the code’s capability to model multiple channels in a full
3D geometry by simulating an array of nanotubes and comparing it with the
results of [2] for a planar CNTFET configuration, as shown in Figure 3.
Initially, we model a single nanotube with periodic boundary conditions in the
lateral direction (domain width $d=3.2$ nm) to emulate an infinite array of
nanotubes. The obtained conductance ($G$) from this simulation compares well
with the results of [2]. Subsequently, we simulate five nanotubes (domain
width $5d$), once again imposing periodic boundary conditions, to ensure
consistent behavior of conductance for each nanotube. This validates the
capability of our model to accurately simulate multiple nanotubes. Later in
§3.4, we will discuss the results from modeling only five CNTs to illustrate
how finite-size effects can impact device performance.

###  3.2 Parallel Performance

In this section, we conduct performance studies of ELEQTRONeX in parallel
using up to 512 MPI ranks with one GPU per MPI rank, focusing on how the code
runtime is impacted as computational resources increase in proportion to
system size $N$. We define computational load by the length of the nanotube
simulated; increase in the load involves increasing the channel length as well
as computational cells used for electrostatic calculation, proportionately.

We first simulate a gate-all-around CNTFET [see inset of Figure 3] with a
carbon nanotube of 54.53 nm overall length ($L$=23.9 nm) using 64 GPUs,
followed by simulations with 2, 4, and 8 times longer nanotubes, employing a
proportional increase in the number of GPUs. We analyze changes in ’time per
iteration’ across three key components of the code—electrostatics, NEGF, and
self-consistency. These times and other parameters are summarized in Table 1.
For this study, we used a number of integration points ($\sim$12k) that was
determined based on the need to obtain numerically converged solution for the
largest channel length, as explained in detail in §3.3. These simulations were
carried out using GPU resources at the National Energy Research Scientific
Computing Center’s (NERSC’s) Perlmutter supercomputer [35], where each node
consists of four NVIDIA A100 GPUs.

Table 1: Parameters used and times obtained from the weak scaling study.

∗Channel length of nanotube modeled / (nm) | 23.9 | 78.4 | 187.4 | 405.6
---|---|---|---|---
Number of carbon rings (system size), $N$ | 512 | 1024 | 2048 | 4096
Number of MPI ranks and GPUs used | 64 | 128 | 256 | 512
Computational cells in the length-wise direction | 1536 | 3072 | 6144 | 12288
Average number of Broyden’s iterations for convergence, $m_{\rm avg}$ | 21 | 22 | 28 | 46
Time for electrostatics / (s) | 1.45 | 1.54 | 1.59 | 1.63
Time for interpolation of potential from mesh to atoms / (ms) | 3.0 | 4.0 | 3.4 | 5.8
Time for charge deposition to mesh / (ms) | 0.63 | 0.67 | 0.79 | 1.20
Time for NEGF / (s) | 2.04 | 3.12 | 5.28 | 8.76
Time for self-consistency / (ms) | 0.50 | 0.59 | 0.74 | 1.00

  * ∗

Contact length was kept constant to 15.336 nm, computational cells in the
transverse direction were $192\times 192$, and approximate number of
integration points used were 12k.

Table 1 illustrates that the time for electrostatics remains nearly constant
indicating near-perfect parallelization of the electrostatic components. The
small 12% increase in the largest case could be attributed to increase in the
overall communication time. The time for the NEGF module increases
proportionally to $N/2$, consistent with parallel time complexity discussed in
§2.2.4, where a factor of $1/2$ is achieved due to overlapping a portion of
the recursive array computations with simultaneous CPU-to-GPU asynchronous
copying. In terms of absolute times, Table 1 reveals that the times per
iteration for Broyden’s algorithm, interpolating the electrostatic potential
from mesh to atomic sites, and depositing charge from atomic locations to
mesh, are orders of magnitude smaller compared to that for electrostatics and
NEGF. Among these, the self-consistency module with the Broyden’s algorithm
exhibits time increase roughly proportional to the average number of
iterations required for convergence, $m_{\rm avg}$, as expected.

In the above calculation, we used only one doubly-degenerate transport mode.
Therefore, next, we demonstrate the effect of increasing the number of modes
on computational time using a simpler simulation involving a CNT embedded in
an all-around circular lead, same as the CNT-lead portion of the gate-all-
around configuration shown in Figure 3. We vary the potential on the lead
metal from to 0 to 1 V, while at each condition, the Broyden iterations are
carried out to achieve self-consistency. Since this is an equilibrium Green’s
function calculation, we simulate a very long CNT (55.84 $\mu$m), and vary the
number of modes from 1 to 8 (all doubly-degenerate). For modes greater than
one, each block of the block-tridiagonal matrix is represented in the form of
an array with each array element corresponding to a mode. Table 2 shows that
with increase in number of modes from 1 to 8, the time for NEGF increases by
only 2.84 times. The sub-linear increase in time may be caused by two major
reasons that benefit calculations on both CPUs and GPUs: better cache and
memory bandwidth utilization, since the data for each mode is stored
contiguously in an array, and better core-pipeline latency hiding due to
increase in the load per thread.

Table 2: Parameters and times for the CNT-lead equilibrium calculations with multiple transport modes. Number of carbon rings (system size), $N$ | 524,288
---|---
Computational cells | $72\times 1,572,864\times 72$
Contour integration points | $90$
Number of modes | 1 | 2 | 4 | 8
Time for NEGF / (s) | 1.98 | 2.75 | 3.24 | 5.63
Time for electrostatics / (s) | 2.2 - 2.4

###  3.3 Long-Channel Gate-All-Around CNTFET Configurations

Here, we discuss the results obtained from the nonequilibrium simulations run
for the parallel performance cases with gate-all-around configurations,
discussed in Table 1, as well as some considerations in choosing integration
paths and points to obtain converged self-consistent solution for longer
nanotubes.

[] [] [] []

Figure 4: (a, b) Current-voltage characteristics for four channel lengths
analyzed for parallel performance studies. (c) Comparison of band bending
across the four cases. (d) Integrand at the channel center for nonequilibrium
charge density calculation. Legends for all figures are shown in (a),
representing channel lengths.

Figures 4 and 4 depict the current-voltage characteristics for the four
simulated cases, while Figure 4 illustrates the band bending at
$V_{gs}=0.8~{}$V and $V_{ds}=-0.1~{}$V. Initially, we compute the $I-V_{ds}$
characteristics at $V_{gs}=0.8~{}$V when the transistor is nearly in the OFF
state. Next, utilizing the charge density profile at $V_{ds}=-0.1~{}$V, we
begin simulations to compute the $I-V_{gs}$ characteristics, sweeping $V_{gs}$
from $0.8~{}$V to $-0.3~{}$V while using the converged solution at the
previous $V_{gs}$ to initialize the computation at the new $V_{gs}$.

In the subthreshold region the $I-V_{ds}$ characteristics reveal a p-type FET
behavior with saturation of the current at -0.1 V. The increase in the drain-
source saturation current magnitude with decreasing channel length is a
consequence of the DIBL discussed earlier. Of note is that the impact of DIBL
on the source-drain current is observed even when the channel is much longer
than the oxide thickness of 10 nm. This impact of the channel length is also
obversed in the subthreshold region of the $I-V_{gs}$ characteristics, while
the ON state current remains comparable to that of the shortest channel
length.

Turning to the question of integration points, we note that the band structure
of an infinitely long carbon nanotube exhibits van hove singularities at the
band edges. For a finite CNT, as the length gets longer, the singularity gets
sharper, requiring more integration points for accurate calculation of the
second part of the integral to compute charge density, shown in equation (4).
To illustrate, the integrand $A_{1}F_{1}+A_{2}F_{2}$ at the center of the
channel is shown in Figure 4 as a function of energy, $E$. We observe that as
the channel length increases, not only does the singularity get sharper, but
it also exhibits oscillating structures on the lower energy side.

Therefore, to reduce the number of integration points, we use a semi-adaptive
method and breakdown the integration path from $E_{min}$ to $E_{max}$ into
three subparts, $E_{min}$ to $E_{l1}$, $E_{l1}$ to $E_{l2}$, and $E_{l2}$ to
$E_{max}$. In each of these subregions we can specify the density of
integration points, defined by the number of integration points per $kT$.
Before starting the simulation, we do not know the exact location of the
singularity, and therefore, we start the simulation with a guess for $E_{l1}$
and $E_{l2}$, compute the integrand at the center of the channel as a function
of energy, obtain the peak of the singularity $E_{s}$, and refine limits as
$E_{l1}=E_{l1}-akT$, $E_{l2}=E_{s}+bkT$, where constants $a$ and $b$ can be
specified by the user, typically set to $3$ and $1$, respectively. The number
of integration points change according to the density of integration points
set by the user in each of these subregions. Note that during the self-
consistency iterations the location and peak of singularity may change, and
therefore, we update these limits at the first iteration for a given set of
conditions and update it at a set period of 50 iterations. For the evaluation
of integration in each subregion, we use Gauss-Legendre quadrature to further
improve the accuracy. For instance, to simulate the largest channel length we
used densities of integration points of (500, 1600, 50), respectively, which
lead to approximately 12,000 integration points in total. For the performance
analyses described above, we kept this total number of integration points
constant for consistent scaling of the workload across various cases. However,
in reality, to obtain a numerically accurate solution at 23.9, 78.4, and 187.4
channel lengths we only need approximately 400, 1000, and 6000 total
integration points, respectively. In fact, the smallest channel length can be
simulated with just one integration path from $E_{min}$ to $E_{max}$. We
suspect that a fully adaptive scheme might be necessary for simulating micron-
scale CNTs self-consistently.

###  3.4 3D Nonparallel Planar CNTFET Configurations

In §3.1, Figure 3 showed the conductance from a planar CNTFET configuration
with periodic boundary conditions in the lateral direction, effectively
emulating an infinite array of CNTs spaced apart by 3.2 nm. Although this
setup is instructive for theoretical analysis, experimental devices typically
contain a finite number of CNTs in the channel, necessitating 3D non-periodic
simulations to reflect more realistic conditions. To showcase the capabilities
of our approach, we simulated five nanotubes in the channel by applying
Neumann boundary conditions instead of periodic ones. An additional spacing of
$d$ was introduced between the lateral boundaries and the outermost nanotubes.

The simulation results, compared in Figure 3, show that the conductance
profiles of the middle three nanotubes closely resemble those obtained with
periodic boundary conditions. However, the conductance of the outermost
nanotubes, labeled as ‘second adjacent,’ is significantly higher. These
results are shown again in Figure 5, where we present the conductance obtained
from the non-periodic configuration as a shaded region around the mean
conductance to provide a clearer visual representation of variability. The
finite CNT case carries more current at a given gate voltage, $11.5\%$ more
current in the ON state (2.9 mS/$\mu$m vs 2.6 mS/$\mu$m), and a reduced
substhrehsold swing (from 2986 mV/dec to 2011 mV/dec) due to the different
gate capacitive coupling.

[] [] [] []

Figure 5: Comparison of zero-bias conductance $G$ versus gate-source voltage
$V_{gs}$ for periodic, non-periodic and parallel, and non-parallel
configurations with different CNT spacing $d$, and channel length $L$. (a)
$d=3.2$ nm, $L=10$ nm, (b) $d=12.8$ nm, $L=10$ nm. For nonperiodic
simulations, we show the mean conductance with the colored spread around it
marking the maximum and minimum conductance for that configuration. (c) and
(d) show electrostatic potential in non-parallel configurations at $V_{gs}=-2$
V.

Next, in the same Figure 5, we also compare the impact of CNT misalignment–a
common occurrence in CNT-array devices as noted in various studies [21, 22].
In this configuration, 5 CNTs are rotated around their midpoint by angles
sampled randomly in the range [$-5^{\circ}$, $5^{\circ}$], such that CNTs do
not touch or overlap other for this 10 nm channel case, as seen from the
contours of electrostatic potential in Figure 5. Such configurations lack
periodicity, requiring fully 3D simulations.

It is seen that the conductance for the non-parallel configuration only
differs slightly from the parallel configuration when the FET is in the ON
state (2.8 mS/$\mu$m vs 2.9 mS/$\mu$m). This small difference is attributed to
the small range of possible angles due to the short channel and close
proximity between CNTs. In addition, the CNT misalignment has little impact on
subthreshold swing (2010 mV/dec vs 2011 mV/dec) or threshold voltage compared
to the parallel 5 CNT case, while the variation in the variation of the
subthreshold swing from each CNT is only about 70 mV/decs.

[] [] [] [] [] []

Figure 6: Planar configurations of 20 CNTs with an average spacing of $d=12.8$
nm and a channel length of $L=100~{}$nm, illustrating the effects of
variations in non-parallel orientation angles and non-uniform spacing on key
characteristics compared to a periodic arrangement. For angular variations:
(a) comparison of $G$ versus $V_{gs}$, with inset showing the same plot on a
$y-log$ scale, (c) electrostatic potential at $V_{gs}=-4$ V on a plane passing
through CNTs, (e) comparison of band structure. For spacing variations:
similar characteristics are shown in (b), (d), and (f), respectively.

Next, in Figure 5, we study the effect of modeling a finite number of parallel
and non-parallel nanotubes when nanotubes are spaced farther apart ($d$=12.8
nm) such that $d$ is larger than the gate-oxide thickness of 10 nm. In this
case, contrary to the simulation with the smaller CNT spacing, the finite
array of 5 parallel CNTs does not show much deviation from the corresponding
infinite array because the CNT-CNT cross-talk is screened by the gate. (ON
stage conductance 3.6 mS/$\mu$m in both cases, subthreshold swing 1778 mV/dec
vs 1785 mV/dec). However, for the case with misalignment with angles in the
range [$-34^{\circ}$, $34^{\circ}$], we observe a decrease of the average ON
state conductance by about 8.3% to 3.4 mS/$\mu$m This decrease is caused by
cross-talk between the CNTs when the spacing between them decreases for part
of the channel length, as seen from the contours of electrostatic potential in
Figure 5.

To understand how the amount of variation in conductance changes with channel
length, we increased the channel length to 100 nm and modeled 20 CNTs with an
average spacing of $d$=12.8 nm and angular rotation within the range
[$-6.5^{\circ}$, $6.5^{\circ}$]. Figures 6, 6, and 6 show conductance,
electrostatic potential, and band structure for this simulation. From the
conductance plot, we again note that the misalignment results in only about 3%
lower conductance with respect to a periodic arrangement when the channel is
ON (3.4 mS/$\mu$m vs 3.5 mS/$\mu$m), and has no impact on the subthreshold
swing (76.7 mV/dec in both cases).

The insert shows the conductance on a y-log scale, revealing a minimum in the
conductance and an increase in the current for higher gate-source voltages.
This increase is caused by band-to-band tunneling as a result of the
conduction band-edge approaching the Fermi level at -1 eV. We further note
that the band-structure obtained from non-parallel simulation differs from a
periodic simulation only when the channel is ON; when there is more charge in
the channel, we see more cross-talk resulting in slightly lower conductance.

The above results suggest that CNT array devices have good resiliency to
misalignment over a broad parameter range. This can be understood primarily
based on screening arguments: CNT-CNT interactions due to misalignment occur
primarily near the source and drain electrodes, but are effectively screened
by the contact metal and the gate.

We also explored the effect of non-uniform spacing between nanotubes while
keeping them parallel. The average spacing is kept constant at $d=12.8$ nm,
and the spacing between each nanotube is varied in the range [0.12d, 0.88d].
Figures 6, 6, and 6 show conductance, electrostatic potential, and band
structure for this case. This simulation showed a broader variation in the
conductance of individual CNTs compared to a periodic simulation, starting
from the linearly increasing region of conductance ($V_{gs}=1~{}V$ to
$-1~{}V$). Yet, the total conductance of the device is close to that of the
periodic system; the maximum variation seen is still only about 3% when the
channel is ON. The band-structure shows a maximum variation with respect to
periodic configuration, when the channel is ON. Variations in the subthreshold
swing or the threshold voltage are also minimal in this case.

Another factor of importance for digital electronics is the CNTFET OFF
current, which has been shown to be dominated by phonon-assisted band-to-band
tunneling [44]. While we plan to add electron-phonon scattering in the future,
here we provide a simple estimate based on the direct tunneling current. In
CNTs, the photon-assisted tunneling current is dominated by scattering with
optical phonons of energy 0.2 eV [45]. Thus, we found the gate-source voltage
at which the conduction band edge is 0.2 eV above the lead Fermi level and
estimated the OFF state conductance from the conductance calculated at this
value without electron-phonon scattering. The results can be found in Table 3.
While the values are approximate, they show that the OFF state current is not
strongly affected by the finite size of the array and misorientation.

Table 3: Impact of non-idealities on planar CNTFETs∗.

L | d | Geometry | $G_{ON}$ | SS | $G_{OFF}$
---|---|---|---|---|---
(nm) | (nm) |  | ($m$S/$\mu$m) | (mV/dec) | ($\mu$S/$\mu$m)
10 | 3.2 | parallel, infinite | 2.6 | 2986 | 6.7
10 | 3.2 | parallel, 5 CNTs | 2.9 $\pm$ 0.21 | 2011 $\pm$ 66 | 6.9 $\pm$ 0.4
10 | 3.2 | misaligned, 5 CNTs | 2.8 $\pm$ 0.17 | 2010 $\pm$ 70 | 6.8 $\pm$ 0.4
10 | 12.8 | parallel, infinite | 3.6 | 1785 | 1.8
10 | 12.8 | parallel, 5 CNTs | 3.6 $\pm$ 3.2e-3 | 1778 $\pm$11 | 1.8 $\pm$ 3.5e-3
10 | 12.8 | misaligned, 20 CNTs | 3.3 $\pm$ 0.2 | 1668 $\pm$56 | 1.4 $\pm$ 0.34
100 | 12.8 | parallel, infinite | 3.5 | 76.7 | 3.9e-4
100 | 12.8 | misaligned, 20 CNTs | 3.4 $\pm$ 0.11 | 76.7 $\pm$ 2.0 | 3.7e-4 $\pm$ 4.6e-5
100 | 12.8 | varying pitch, 20 CNTs | 3.4 $\pm$ 0.15 | 77.2 $\pm$ 1.0 | 4.2e-4 $\pm$ 4.4e-5

  * ∗

Oxide thickness was kept constant to 10 nm. For cases with multiple nanotubes,
one standard deviation variation is shown.

In terms of computational cost, the simulations with 20 CNTs employed 10560,
768, and 768 cells in the X, Y, and Z directions, respectively, with a cell
size of 0.026625 nm in each direction. The computation utilized 512 MPI ranks
each associated with one GPU, averaging approximately 3.6 s per iteration for
electrostatics and 0.65 seconds per iteration with 90 integration points for
NEGF, while the time required for other steps in the iteration was an order of
magnitude lower. The code required 42 iterations on average, i.e.
approximately 3 minutes per gate-source voltage condition.

##  4 Conclusions

In summary, we present a robust, portable, and open-source implementation of
the NEGF method self-consistently coupled with electrostatics, employing the
AMReX library and MPI/GPU parallelization strategies to leverage modern
supercomputing architectures. We demonstrate rigorous validation of our code
against established literature for computing DC transport properties across
various CNTFET configurations, including multi-nanotube simulations.

The performance studies underscore ELEQTRONeX’s computational efficiency,
particularly in modeling larger systems such as long device channels and
larger computational domains for electrostatics. The time complexity of the
multigrid electrostatics is $O(N)$ for a system of size $N$, while the
parallel implementation achieves a near-ideal $O(1)$ time complexity, bounded
only by inter-processor communication. The time complexity of the serial NEGF
algorithm is $O(N^{2})$, yet we demonstrate that it can achieve $O(N)$
performance with heterogeneous MPI and GPU implementation, while optimizations
such as asynchronous CPU-to-GPU copying can further halve the scaling factor.
Additionally, we introduce an MPI/GPU parallel implementation of Broyden’s
modified second algorithm for self-consistency, enhancing parallel time
complexity to $O(m_{\rm avg})$ compared to $(Nm_{\rm avg})$ serial
implementation, where $m_{\rm avg}$ represents average number iterations
required for convergence.

Furthermore, we demonstrate our capability to model fully 3D non-periodic
configurations of CNTFETs and investigate the accuracy of simple periodic
calculations compared to realistic experimental setups. Such setups may
involve a finite number of nanotubes that are not perfectly parallel or
equally spaced. Across a broad range of configurations and device geometries
we find that CNTFETs are robust against misalignment and varying pitch for
several key device performance metrics.

In the future we plan to add several additional capabilities to ELEQTRONeX
including electron-phonon interactions, trap charges, top-gate geometries,
SRAM cell configurations, bulk materials, and time-dependent NEGF to address
computationally challenging state-of-the-art problems.

## Acknowledgements

Work by S. S. Sawant and Z. Yao was supported by the U.S. Department of
Energy, Office of Science, Office of Applied Scientific Computing Research,
the Microelectronics Co-Design Research Program, under Contract DE-
AC02-05-CH11231 (Co-Design and Integration of Nano-sensors on CMOS). Work by
F. Léonard was supported by the U.S. Department of Energy, Office of Science,
the Microelectronics Co-Design Research Program, under Contract DE-NA-0003525
(Co-Design and Integration of Nano-sensors on CMOS). Sandia National
Laboratories is a multimission laboratory managed and operated by National
Technology and Engineering Solutions of Sandia, LLC., a wholly owned
subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-
NA-0003525. Work by A. Nonaka was supported by the Exascale Computing Project
(17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office
of Science and the National Nuclear Security Administration; the U.S.
Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research, Exascale Computing Project under contract DE-
AC02-05CH11231. This research used resources of the National Energy Research
Scientific Computing Center (NERSC), a Department of Energy Office of Science
User Facility using NERSC award ASCR-ERCAP0026882.

Authors would also like to thank Dr. Dmitri Nikonov for valuable discussions
during the development of the solver.

## References

  * [1] F. Léonard, D. A. Stewart, Properties of short channel ballistic carbon nanotube transistors with ohmic contacts, Nanotechnology 17 (18) (2006) 4699.  doi:10.1088/0957-4484/17/18/029.
  * [2] F. Léonard, Crosstalk between nanotube devices: contact and channel effects, Nanotechnology 17 (9) (2006) 2381.  doi:10.1088/0957-4484/17/9/051.
  * [3] Y. M. Banadaki, A. Srivastava, Investigation of the width-dependent static characteristics of graphene nanoribbon field effect transistors using non-parabolic quantum-based model, Solid-State Electronics 111 (2015) 80–90.  doi:https://doi.org/10.1016/j.sse.2015.05.003.
  * [4] K. Ganapathi, Y. Yoon, S. Salahuddin, Monolayer MoS2 transistors - ballistic performance limit analysis, in: 69th Device Research Conference, 2011, pp. 79–80.  doi:10.1109/DRC.2011.5994421.
  * [5] B. Nag Chowdhury, S. Chattopadhyay, Dual-gate GaAs-nanowire FET for room temperature charge-qubit operation: A NEGF approach, Advanced Quantum Technologies 6 (4) (2023) 2200072.  doi:https://doi.org/10.1002/qute.202200072.
  * [6] H.-H. Park, W. Choi, M. A. Pourghaderi, J. Kim, U. Kwon, D. S. Kim, NEGF simulations of stacked silicon nanosheet FETs for performance optimization, in: 2019 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 2019, pp. 1–3.  doi:10.1109/SISPAD.2019.8870365.
  * [7] G. Cohen, M. Galperin, Green’s function methods for single molecule junctions, The Journal of Chemical Physics 152 (9) (2020) 090901.  arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5145210/15572359/090901\\_1\\_online.pdf, doi:10.1063/1.5145210.
  * [8] D. Kienle, F. Léonard, Terahertz response of carbon nanotube transistors, Physical review letters 103 (2) (2009) 026601.  doi:https://doi.org/10.1103/PhysRevLett.103.026601.
  * [9] F. Léonard, C. D. Spataru, M. Goldflam, D. W. Peters, T. E. Beechem, Dynamic wavelength-tunable photodetector using subwavelength graphene field-effect transistors, Scientific Reports 7 (APR 4 2017).  doi:10.1038/srep45873.
  * [10] D. A. Stewart, F. Léonard, Photocurrents in nanotube junctions, Phys. Rev. Lett. 93 (2004) 107401.  doi:10.1103/PhysRevLett.93.107401.
  * [11] N. Papior, N. Lorente, T. Frederiksen, A. García, M. Brandbyge, Improvements on non-equilibrium and transport green function techniques: The next-generation TRANSIESTA, Computer Physics Communications 212 (2017) 8–24.  doi:10.1016/j.cpc.2016.09.022.
  * [12] N. Papior, T. Gunst, D. Stradi, M. Brandbyge, Manipulating the voltage drop in graphene nanojunctions using a gate potential, Physical Chemistry Chemical Physics 18 (2) (2016) 1025–1031.  doi:10.1039/C5CP04613K.
  * [13] S. Smidstrup, T. Markussen, P. Vancraeyveld, J. Wellendorff, J. Schneider, T. Gunst, B. Verstichel, D. Stradi, P. A. Khomyakov, U. G. Vej-Hansen, et al., QuantumATK: an integrated platform of electronic and atomic-scale modelling tools, Journal of Physics: Condensed Matter 32 (1) (2019) 015901.  doi:10.1088/1361-648X/ab4007.
  * [14] G. Klimeck, M. Luisier, Atomistic modeling of realistically extended semiconductor devices with NEMO and OMEN, Computing in Science & Engineering 12 (2) (2010) 28–35.  doi:10.1109/MCSE.2010.32.
  * [15] S. Steiger, M. Povolotskyi, H.-H. Park, T. Kubis, G. Klimeck, Nemo5: A parallel multiscale nanoelectronics modeling tool, IEEE Transactions on Nanotechnology 10 (6) (2011) 1464–1474.
  * [16] H. Sahasrabudhe, J. Fonseca, G. Klimeck, Accelerating nano-scale transistor innovation with NEMO5 on Blue Waters, Tech. rep., University of Illinois at Urbana-Champaign, https://bluewaters.ncsa.illinois.edu/liferay-content/document-library/NEIS-P2-C1-Reports/Klimeck_NEISP2_Final_Report.pdf (2014).
  * [17] Y. Jeong, H. Ryu, High performance simulations of quantum transport using manycore computing, in: The International Conference on High Performance Computing in Asia-Pacific Region Companion, HPCAsia’21 Companion, Association for Computing Machinery, New York, NY, USA, 2021, p. 21–28.  doi:10.1145/3440722.3440879.
  * [18] Alessandro Pecchia et al., libNEGF: A library for Non-Equilibrium Green’s Function (NEGF) simulations, https://github.com/libnegf/libnegf (Accessed: April 2024).
  * [19] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, et al., AMReX: a framework for block-structured adaptive mesh refinement, Journal of Open Source Software 4 (37) (2019) 1370–1370.  doi:10.21105/joss.01370.
  * [20] W. Zhang, A. Myers, K. Gott, A. Almgren, J. Bell, AMReX: Block-structured adaptive mesh refinement for multiphysics applications, The International Journal of High Performance Computing Applications 35 (6) (2021) 508–526.  doi:10.1177/10943420211022811.
  * [21] K. R. Jinkins, S. M. Foradori, V. Saraswat, R. M. Jacobberger, J. H. Dwyer, P. Gopalan, A. Berson, M. S. Arnold, Aligned 2D carbon nanotube liquid crystals for wafer-scale electronics, Science Advances 7 (37) (2021) eabh0640.  doi:10.1126/sciadv.abh0640.
  * [22] L. Liu, J. Han, L. Xu, J. Zhou, C. Zhao, S. Ding, H. Shi, M. Xiao, L. Ding, Z. Ma, C. Jin, Z. Zhang, L.-M. Peng, Aligned, high-density semiconducting carbon nanotube arrays for high-performance electronics, Science 368 (6493) (2020) 850–856.  doi:10.1126/science.aba5980.
  * [23] J. Guo, S. Datta, M. Lundstrom, M. Anantam, Toward multiscale modeling of carbon nanotube transistors (2004).  doi:10.1615/IntJMultCompEng.v2.i2.60.
  * [24] AMReX team, AMReX (2024).
URL https://github.com/AMReX-Codes/amrex

  * [25] ELEQTRONeX team, ELEQTRONeX (2024).
URL https://github.com/AMReX-Microelectronics/ELEQTRONeX

  * [26] AMReX Team, AMReX Documentation, Release 24.03-dev, Lawrence Berkeley National Laboratory, https://amrex-codes.github.io/amrex/docs_html/ (2024).
  * [27] A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, M. L. Welcome, A conservative adaptive projection method for the variable density incompressible Navier–Stokes equations, Journal of computational Physics 142 (1) (1998) 1–46.  doi:10.1006/jcph.1998.5890.
  * [28] H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM Journal on Scientific and Statistical Computing 13 (2) (1992) 631–644.  doi:10.1137/0913035.
  * [29] R. D. Falgout, U. M. Yang, hypre: A library of high performance preconditioners, in: P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, J. J. Dongarra (Eds.), Computational Science — ICCS 2002, Springer Berlin Heidelberg, Berlin, Heidelberg, 2002, pp. 632–641.  doi:10.1007/3-540-47789-6_66.
  * [30] C. K. Birdsall, D. Fuss, Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation, Journal of Computational Physics 3 (4) (1969) 494–511.  doi:10.1016/0021-9991(69)90058-8.
  * [31] J. Taylor, H. Guo, J. Wang, Ab initio modeling of quantum transport properties of molecular electronic devices, Physical Review B 63 (24) (2001) 245407\.  doi:10.1103/PhysRevB.63.245407.
  * [32] M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, K. Stokbro, Density-functional method for nonequilibrium electron transport, Physical Review B 65 (16) (2002) 165401.  doi:10.1103/PhysRevB.65.165401.
  * [33] M. P. L. Sancho, J. M. L. Sancho, J. M. L. Sancho, J. Rubio, Highly convergent schemes for the calculation of bulk and surface green functions, Journal of Physics F: Metal Physics 15 (4) (1985) 851.  doi:10.1088/0305-4608/15/4/009.
  * [34] J.-S. Wang, N. Zeng, J. Wang, C. K. Gan, Nonequilibrium Green’s function method for thermal transport in junctions, Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 75 (6) (2007) 061128.  doi:10.1103/PhysRevE.75.061128.
  * [35] System details, https://docs.nersc.gov/systems/perlmutter/system_details/ (2021).
  * [36] E. M. Godfrin, A method to compute the inverse of an n-block tridiagonal quasi-Hermitian matrix, Journal of Physics: Condensed Matter 3 (40) (1991) 7843\.  doi:10.1088/0953-8984/3/40/005.
  * [37] O. Hod, J. E. Peralta, G. E. Scuseria, First-principles electronic transport calculations in finite elongated systems: A divide and conquer approach, The Journal of chemical physics 125 (11) (2006).  doi:10.1063/1.2349482.
  * [38] G. P. Srivastava, Broyden’s method for self-consistent field convergence acceleration, Journal of Physics A: Mathematical and General 17 (6) (1984) L317.  doi:10.1088/0305-4470/17/6/002.
  * [39] D. Singh, H. Krakauer, C. Wang, Accelerating the convergence of self-consistent linearized augmented-plane-wave calculations, Physical Review B 34 (12) (1986) 8391.  doi:10.1103/PhysRevB.34.8391.
  * [40] S. Ihnatsenka, I. Zozoulenko, M. Willander, Electron-electron interaction effects in transport through open dots: Pinning of resonant levels, Physical Review B 75 (23) (2007) 235307.  doi:10.1103/PhysRevB.75.235307.
  * [41] D. A. Areshkin, B. K. Nikolić, Electron density and transport in top-gated graphene nanoribbon devices: First-principles Green function algorithms for systems containing a large number of atoms, Physical Review B 81 (15) (2010) 155450\.  doi:10.1103/PhysRevB.81.155450.
  * [42] W. H. Wen-Mei, D. B. Kirk, I. El Hajj, Programming Massively Parallel Processors: A Hands-on Approach, 4th Edition, Morgan Kaufmann, 2022.
  * [43] NVIDIA, Using CUDA Warp-level primitives, https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/, accessed: 2023-05-18 (2023).
  * [44] Q. Lin, C. Gilardi, S.-K. Su, Z. Zhang, E. Chen, P. Bandaru, A. Kummel, I. Radu, S. Mitra, G. Pitner, H.-S. P. Wong, Band-to-band tunneling leakage current characterization and projection in carbon nanotube transistors, ACS Nano 17 (21) (2023) 21083–21092.  doi:10.1021/acsnano.3c04346.
  * [45] S. O. Koswatta, M. S. Lundstrom, M. P. Anantram, D. E. Nikonov, Simulation of phonon-assisted band-to-band tunneling in carbon nanotube field-effect transistors, Applied Physics Letters 87 (25) (2005) 253107.  doi:10.1063/1.2146065.