Pushing Back the Limit of Ab-initio Quantum Transport Simulations on Hybrid Supercomputers

TL;DR

This paper demonstrates a significant advancement in ab-initio quantum transport simulations by integrating CP2K and OMEN on hybrid supercomputers, achieving large-scale, high-efficiency calculations with novel algorithms.

Contribution

The development of the SplitSolve algorithm enables efficient hybrid CPU-GPU computations, significantly reducing simulation time and scaling to large nanostructures.

Findings

Over one order of magnitude reduction in time-to-solution

Simulation of structures with over 50,000 atoms

Parallel efficiency of 97% on thousands of nodes

Abstract

The capabilities of CP2K, a density-functional theory package and OMEN, a nano-device simulator, are combined to study transport phenomena from first-principles in unprecedentedly large nanostructures. Based on the Hamiltonian and overlap matrices generated by CP2K for a given system, OMEN solves the Schroedinger equation with open boundary conditions (OBCs) for all possible electron momenta and energies. To accelerate this core operation a robust algorithm called SplitSolve has been developed. It allows to simultaneously treat the OBCs on CPUs and the Schroedinger equation on GPUs, taking advantage of hybrid nodes. Our key achievements on the Cray-XK7 Titan are (i) a reduction in time-to-solution by more than one order of magnitude as compared to standard methods, enabling the simulation of structures with more than 50000 atoms, (ii) a parallel efficiency of 97% when scaling from 756 up to 18564 nodes, and (iii) a sustained performance of 15 DP-PFlop/s.

Tables

Table 1. TABLE I: Technical specifications of Piz Daint and Titan.

	Piz Daint [49]	Titan [50]
Hybrid nodes	5272	18688
GPUs	5272	18688
GPU model	Tesla K20X	Tesla K20X
Cores	42176	299008
CPU model	Intel Xeon E5-2670	Opteron 6274
Node perf. (GFlop/s)	166.4+1311	134.4+1311

Table 2. TABLE II: Weak scaling data corresponding to Fig.  11 (a).

Hybrid Nodes	Time (s)	Avg. $E$/node	Avg. Time/$E$ (s)
588	1277	14.1	90.8
1176	1197	13.4	89
2352	1281	13.8	92.7
4704	1213	13.8	87.7
9408	1204	13.3	90.3
18564	1130	12.9	87.5

Table 3. TABLE III: Strong scaling data corresponding to Fig.  11 (b) and to the run that reached 15 PFlop/s (last line).

Hybrid Nodes	Time (s)	$\left|\right|$ Eff. (%)	PFlop/s
756	26975	100	0.54
1512	13593	99.2	1.06
3024	6806	99.1	2.12
6048	3415	98.7	4.23
12096	1711	98.5	8.45
18564	1130	97.3	12.8
18564	912.5	-	15.01

Equations

$$\left(-\frac{{\hslash }^2}{2m_0}{\nabla }^2+V(r)+V_H(r)+V_{xc}(r)\right)\psi (r)=E\psi (r)$$

$$\psi (r)$$

$$H\cdot c$$

$$\left(E\cdot S-H-{\Sigma }^{RB}(E)\right)\cdot G^R(E)$$

$$\left(E\cdot S-H-{\Sigma }^{RB}\right)\cdot c(E)$$

$$\sum _{l=-N_{\mathrm{BW}}}^{+N_{\mathrm{BW}}}{\text{e}}^{i\cdot l\cdot k_B}(H_{q,q+l}-E\cdot S_{q,q+l})u_B=0,$$

$$\left[Q_F^{†}{A}_F{Q}_F\right]u_B$$

$$A_F$$

$$B_F$$

$$Q_F$$

$$(A+BC{)}^{-1}=A^{-1}-A^{-1}B(1+C{A}^{-1}B{)}^{-1}C{A}^{-1}.$$

$$x$$

$$=$$

Full text

Pushing Back the Limit of Ab-initio Quantum Transport

Simulations on Hybrid Supercomputers

Mauro Calderara*∗, Sascha Brück∗, Andreas Pedersen∗, Mohammad H. Bani-Hashemian†*,

Joost VandeVondele*†, and Mathieu Luisier∗*

*∗*Integrated Systems Laboratory, ETH Zürich, 8092 Zürich, Switzerland

*†*Nanoscale Simulations, ETH Zürich, 8093 Zürich, Switzerland

ABSTRACT

The capabilities of CP2K, a density-functional theory package and OMEN, a nano-device simulator, are combined to study transport phenomena from first-principles in unprecedentedly large nanostructures. Based on the Hamiltonian and overlap matrices generated by CP2K for a given system, OMEN solves the Schrödinger equation with open boundary conditions (OBCs) for all possible electron momenta and energies. To accelerate this core operation a robust algorithm called SplitSolve has been developed. It allows to simultaneously treat the OBCs on CPUs and the Schrödinger equation on GPUs, taking advantage of hybrid nodes. Our key achievements on the Cray-XK7 Titan are (i) a reduction in time-to-solution by more than one order of magnitude as compared to standard methods, enabling the simulation of structures with more than 50000 atoms, (ii) a parallel efficiency of 97% when scaling from 756 up to 18564 nodes, and (iii) a sustained performance of 15 DP-PFlop/s.

1. INTRODUCTION

The fabrication of nanostructures has considerably improved over the last couple of years and is rapidly approaching the point where individual atoms are reliably manipulated and assembled according to desired patterns. There is still a long way to go before such processes enter mass production, but exciting applications have already been demonstrated: a low-temperature single-atom transistor [1], atomically precise graphene nanoribbons [2], or van der Waals heterostructures based on metal-dichalcogenides [3] were recently synthesized.

Despite these promising advances the realization of nano-devices remains a very tedious task, more complicated than it was at the micrometer scale. Researchers can no longer rely on their sole intuition and past experience to conceive properly working nanostructures. Doing so could lead them to wrong assumptions or make them miss relevant physical effects. Advanced technology computer aided design (TCAD) platforms are needed to support the experimental work and accelerate the emergence of novel device concepts. This requires the development of accurate simulation approaches, whose key ingredient is the bandstructure model they rely on. The latter accounts for the material properties of the considered systems, which might completely determine the device functionality.

The effective mass approximation (EMA), k$\cdot$p models [4], tight-binding [5], or pseudopotential methods [6] offer a satisfactory level of accuracy in many applications, but generally they suffer from several deficiencies such as their necessary parameterization (all are empirical models), the transferability of the parameters from bulk to nanostructures, the treatment of heterostructures, or the absence of atomic resolution for EMA and k$\cdot$p. Ab-initio methods, e.g. density-functional theory (DFT) based on the Kohn-Sham equations [7] appear to be more promising solutions since they address the shortcomings mentioned above and their known band gap underestimations can be corrected with hybrid functionals [8].

First-principles codes such as VASP [9], ABINIT [10], Quantum ESPRESSO [11], SIESTA [12], or CP2K [13] extensively use DFT and therefore lend themselves ideally to the calculation of electronic and crystal structures, phase diagrams, charge densities, or vibrational frequencies in solids. However, they are not well-suited to deal with out-of-equilibrium situations, where, for example, an external voltage or temperature gradient is applied to a nanostructure, inducing an electron current between its contacts. Being able to rapidly and efficiently engineer the magnitude and the direction of this current is a goal of utmost importance in transistor [14], molecular switch [15], quantum well solar cell [16], quantum dot light-emitting diode [17], nanowire thermoelectric generator [18], switching resistive memory [19], or lithium-ion battery [20] research. All these devices rely on optimized current flows that can be directly compared to measurements, contrary to electronic structures, charge distributions, or potential profiles.

To study the transport rather than the static properties of matter, standard DFT packages must be augmented with quantum transport (QT) simulation capabilities, which implies replacing the computation of large eigenvalue problems by the solution of linear systems of equations, as produced by the versatile Non-equilibrium Green’s Function (NEGF) formalism [21]. While ab-initio tools combining DFT and NEGF have been formally demonstrated many years ago [22, 23, 24, 25], they remain computationally very intensive, thus limiting their application to small systems composed of up to 1000 atoms like molecules, nanotubes, or nanoribbons, all simulated in the ballistic limit of transport (no scattering) [26, 27, 28, 29, 30]. There are two notable exceptions where the authors claim to have reached 20000 atoms with a first-principles NEGF scheme [31, 32]. In both studies a linearized muffin-tin-orbital (LMTO) basis is employed, which exhibits a tight-binding-like sparsity pattern that is mostly suitable for close-packed metallic crystals, not for arbitrary materials, as inspected here.

To allow for the simulation of realistic nano-devices made of any elements we propose in this paper a new ab-initio quantum transport solver that goes beyond the existing solutions. The selected approach links two state-of-the-art, massively parallel codes, CP2K [13] and OMEN [33], and leverages their DFT and transport capabilities, respectively, in order to handle systems composed of tens of thousands of atoms. The main target applications are the design of nanoscale transistors and the enhancement of the electronic conductivity in lithium ion battery electrodes, as shown in Fig. 1. These are two research areas with a tremendous need for novel simulation tools. The implementation of the code is general enough to treat other types of nanostructures too.

OMEN already possesses advanced numerical algorithms capable of breaking the petascale barrier while performing quantum transport calculations, but they are optimized for tight-binding bases, not for DFT ones, and they are restricted to CPUs only [33]. Hence, a more powerful sparse linear solver called SplitSolve has recently been developed and integrated into OMEN. It supports the usage of CPUs and GPUs at the same time, significantly reducing the time-to-solution of DFT+QT problems and opening the door for the exploration of larger design spaces. With SplitSolve, the simulation of a Si nanowire transistor composed of 55488 atoms has been successfully achieved, which is, to the best of the authors’ knowledge, at least 10 times larger than what others have reported so far in the literature for similar 3-D semiconducting structures. On the Cray-XK7 Titan at Oak Ridge National Laboratory (ORNL) we have also been able to demonstrate a sustained double-precision performance of 15 PFlop/s in production mode.

The paper is organized as follows: in Section 2, an overview of the ab-initio quantum transport approach is given, followed by the algorithmic innovation in Section 3. A short description of the CP2K and OMEN applications is proposed in Section 4. The time-to-solution, scalability, and peak performance of the code are presented and analyzed in Section 5 before conclusion.

2. SIMULATION APPROACH

A. Density-Functional Theory with a Localized Basis

The coupling between OMEN and CP2K is schematized in Fig. 2. For a given nanostructure, CP2K starts by solving the Kohn-Sham DFT equation [7]

[TABLE]

where the first term refers to the electron kinetic operator, the second one, $V(r)$, to the electron-nuclei interactions, the third one, $V_{H}(r)$, to the Hartree (Coulomb) potential, and the last one, $V_{xc}(r)$, to the exchange-correlation energy. The wave function $\psi(r)$ is expanded in a localized basis made of contracted Gaussian orbitals

[TABLE]

In Eq. (2) $\phi_{\mu}(r)$ is a contracted Gaussian function of type $\mu$ and $c_{\mu}$ the corresponding expansion coefficient. Inserting this expression into the Kohn-Sham equation Eq. (1) gives rise to a generalized eigenvalue problem

[TABLE]

with the Hamiltonian $H$ and overlap $S$ matrices as well as the unknown expansion coefficients $c$ (now a vector) and eigenvalues $E$. All results presented later in this paper rely on a 3SP basis within the local density approximation (LDA) to model the exchange-correlation energy [34]. However, the SplitSolve algorithm works with any basis set and functional, as shown in Fig. 1: the current flowing through a lithiated SnO battery anode was calculated with a double-$\zeta$ basis and the PBE general gradient approximation (GGA) [35], while the transmission probability through a Si nanowire was obtained with the hybrid HSE06 functional [8].

B. Quantum Transport

For efficient simulations of transport through nanostructures, it is very convenient to work within a localized basis, as the one provided by CP2K. A sparse matrix, usually block tri-diagonal, then describes the on-site and inter-atomic interactions. The main difference to Eq. (3) is that open boundary conditions (OBCs) must be introduced to inject electrons at the predefined contacts. These so-called reservoirs or leads might experience different chemical potentials, depending on the externally applied bias. This split induces a current flow. The resulting system of equations has the following form in the Non-equilibrium Green’s Function formalism

[TABLE]

where $S$ and $H$ are directly imported from CP2K, the OBCs are cast into the boundary self-energy $\Sigma^{RB}(E)$, $\mathbb{1}$ is the identity matrix, and the unknowns are the retarded Green’s Functions $G^{R}(E)$. In the ballistic limit of transport, as here, it is computationally more efficient to transform Eq. (4) into the Wave Function formalism, which takes the form of a linear system of equations [38]

[TABLE]

The vector $Inj(E)$ denotes the injection mechanism into the out-of-equilibrium devices, whereas $c(E)$ has the same meaning as in Eq. (3). The structure of Eq. (5) is plotted in Fig. 4, highlighting its sparse linear pattern. The charge and current densities can be directly derived from $G^{R}(E)$ or $c(E)$ after Eq. (4) or (5) have been solved for all potential electron energies $E$ and wave vectors $k$ in cases of periodicity along at least one direction [39].

The matrices $H$ and $S$ from CP2K contain about 100 times more non-zero entries than their tight-binding or LMTO counterparts, as shown in Fig. 3. This is why the standard algorithms of OMEN [33] do not perform well in large ab-initio quantum transport calculations. We note that CP2K currently does not provide any momentum- or $k$-dependence for $H$ and $S$ in periodic systems. This issue is resolved by first cutting all the needed blocks from 3-D simulations and then generating $H(k)$ and $S(k)$ in OMEN.

3. ALGORITHMIC INNOVATIONS

A. Open Boundary Conditions: the FEAST Algorithm

The self-energy $\Sigma^{RB}$ and the injection vector $Inj$ in Eq. (5) are calculated from the wave vectors $k_{B}$ and eigenmodes $u_{B}$ of the leads/reservoirs of the considered systems [39]. The following polynomial eigenvalue problem must be solved to obtain them

[TABLE]

where $H_{q,q+l}$ and $S_{q,q+l}$ are the parts of the Hamiltonian and overlap matrix that describe the interaction of unit cell $q$ and $q+l$ within the leads and $N_{\mathrm{BW}}$ indicates the range of the inter-cell interactions, typically $N_{\mathrm{BW}}\geq$2. It has been demonstrated that in a tight-binding basis, $\Sigma^{RB}$ can be more rapidly evaluated with Eq. (6) and a shift-and-invert approach [38] than with the standard iterative decimation technique used in NEGF [40]. However, in a DFT basis, the open boundary conditions start to represent a serious computational bottleneck due to the size increase of the involved matrices and the difficulty to parallelize the shift-and-invert method.

In practice it is not necessary to determine the $N_{BC}$ possible phase factors $\lambda=\textrm{e}^{i\cdot k_{B}}$ as the contribution from fast decaying modes is negligible. It suffices to find the $m$ eigenvalues inside an annulus around $|\lambda|=1$ in the complex plane centered at the origin. Hence, a contour integration method can be employed to find a subspace projector $Q_{F}$ that spans the same space as the eigenvectors corresponding to the eigenvalues enclosed by the shaded region in Fig. 5. The FEAST algorithm [41] and its extension to generalized eigenvalue problems with non-Hermitian matrices [42] has been chosen to produce $Q_{F}$. When applying the Rayleigh-Ritz method, the eigenvalue problem reduces to a $m\times m$ system

[TABLE]

with the following matrix definitions

[TABLE]

[TABLE]

[TABLE]

Here, $\tilde{H}_{j}=H_{q,q+j}-E\cdot S_{q,q+j}$, the not-shown entries are equal to 0, the size of $A_{F}$ and $B_{F}$ is $N_{BC}\times N_{BC}$, that of $Q_{F}$ and $Y_{F}$ (matrix of random numbers) $N_{BC}\times m$ with $N_{BC}\gg m$, and the $z_{p}$ are the integration points in the trapezoid rule. The linear systems of equations in Eq. (18) determine the execution time of FEAST. Through an analytical block LU decomposition, their size can be decreased to $N_{BC}/(2N_{\mathrm{BW}})$. Furthermore, they can be solved in parallel since the solution at each of the $N_{p}$ integration points is independent from the others. FEAST can be modified according to Ref. [43] to further reduce the calculation time.

B. Schrödinger Equation: the SplitSolve Algorithm

A parallel direct sparse linear solver such as MUMPS [46] or a custom-made block cyclic reduction (BCR) [33] are typically needed to solve the Schrödinger equation with OBCs. Iterative solvers cannot be efficiently used due to the presence of multiple right-hand-sides. Since our BCR method relies on the sparsity provided by a tight-binding basis, it does not work with DFT, whereas MUMPS becomes slow when the number of non-zero entries in the Hamiltonian and overlap matrices increases drastically. To address these issues we have developed the SplitSolve algorithm that can solve Eq. (5) on accelerators (GPUs or others) and leverage its particular structure displayed in Fig. 4. The matrix $T=(E\cdot S-H-\Sigma^{RB})$ is of size $N_{SS}$ (total number of atoms times orbital per atom) with a block tridiagonal shape. The right-hand-side $Inj$ has multiple columns and non-zero elements only in the top and bottom block rows. Consequently, solving such a system is equivalent to obtaining the first and last block columns of $T^{-1}$ and multiplying them with the non-zero rows of $Inj$.

The goal of SplitSolve is twofold: (i) efficiently computing only the required parts of $T^{-1}$ and (ii) decoupling the calculation of the open boundary conditions $\Sigma^{RB}$ from the solution of $T^{-1}$. This can be achieved by choosing a suitable low rank product representation of the boundary self-energy $\Sigma^{RB}\vcentcolon=BC$, as in Fig. 4, and by recalling the Sherman-Morrison-Woodbury formula

[TABLE]

By doing so the computation of $\Sigma^{RB}$ can be interleaved with the solution of the full problem. A similar approach has been successfully tested in Refs. [44, 45] for tight-binding. Defining $x\vcentcolon=c$, $b\vcentcolon=Inj$ and introducing $Q\vcentcolon=A^{-1}B$ and $y\vcentcolon=A^{-1}b$, it follows that

[TABLE]

where $R\vcentcolon=(\mathbb{1}-CQ)$ and $z\vcentcolon=R^{-1}Cy$. The final solution $x$ is then obtained in four steps

• •

Step 1 Solve $AQ=B$ for $Q$.

• •

Step 2 Solve $Ay=b$ for $y$.

• •

Step 3 Solve $Rz=(\mathbb{1}-CQ)z=Cy$ for $z$.

• •

Step 4 Compute the full solution as $x=y+Qz$.

The choice of the $B$ and $C$ matrices is free to a large extent, but critical to minimize the computational burden. In SplitSolve $B$ is a $N_{SS}\times 2s$ zero matrix, except the $s\times s$ top left and bottom right subblocks which are equal to $\mathbb{1}$, while $C$ is a $2s\times N_{SS}$ matrix with the top left block set to $\Sigma^{RB,L}$ and the bottom right one to $\Sigma^{RB,R}$, the boundary self-energies of the left ($L$) and right ($R$) contacts, both of size $s\times s$. Given this selection of $B$ and $C$ the following observations are made about SplitSolve:

• •

The algorithm can be decomposed into a pre- and post-processing phase. The former consists of Step 1 and can be computed in parallel with the evaluation of $C$ and $b$. The latter comprises Steps 2, 3, and 4 and must be performed after the completion of the $Q$, $C$ and $b$ computation.

• •

The structure of $A$ prior to the addition of the boundary conditions is preserved and can be leveraged in the preprocessing phase: $A$ is usually real symmetric in 3-D structures and complex Hermitian in 1-D and 2-D.

• •

In terms of numerics, $Q$ consists of the first and last $s$ columns of $A^{-1}$, $y$ can be obtained as $Q\cdot b$, $R$ is a system of comparably small size $2s\times 2s$, and $x$ can be computed as $x=Q\cdot(b^{\prime}+z)$ where $b^{\prime}$ denotes the nonzero rows of $b$.

**SplitSolve Preprocessing: ** To calculate the first and last $s$ columns of $A^{-1}$ in Step 1, a modified version of the recursive Green’s Function (RGF) algorithm [47] has been implemented. If $A_{i,j}$ denotes the $i^{\textrm{th}}$ block row and $j^{\textrm{th}}$ block column of the block tridiagonal matrix $A$ with $n_{B}$ diagonal blocks, then the preprocessing part of SplitSolve can be summarized as in Algorithm 1.

This gives $Q_{i,1:s}$=$A^{-1}_{i,1}$, the first block column of $A^{-1}$. The last one, $Q_{i,s+1:2s}$=$A^{-1}_{i,n_{B}}$, can be derived similarly and in parallel with $A^{-1}_{i,1}$ since the two tasks are independent from each other and naturally scale to two accelerators. Alternatively, the matrix $A$ can be partitioned horizontally into two parts of equal size and the preprocessing in SplitSolve includes four phases, $P_{1}$ to $P_{4}$, as illustrated in Fig. 6.

To study transport through realistic nanostructures it is crucial to parallelize SplitSolve beyond two accelerators so that large system of equations can be efficiently solved. In our implementation, the number of partitions must be a power of 2. Each pair of accelerators computes the first and last block columns of the inverse corresponding to its local partition using Algorithm 1. In a second step adjacent partitions are recursively merged together based on a modified and optimized variant of the SPIKE algorithm [48]. With the proposed scheme the merging steps have a constant cost, their computation is evenly distributed over all the accelerators, and their number grows logarithmically with the number of partitions. The parallelization incurs very little memory overhead so that the extra memory gained by the presence of additional accelerators allows to treat bigger devices.

SplitSolve Postprocessing: Upon availability of $Q$, $\Sigma^{RB}$, and $b$, the matrices $R$ and $Cy$ are constructed on the two accelerators storing the first and last partition and calculating $z$. This quantity is then distributed to all the accelerators to solve Eq. (21) with one single matrix-vector multiplication per block, i.e. $x=Q\cdot(b^{\prime}+z)$.

C. Node-Level FEAST and SplitSolve Performance

The performance of the newly implemented FEAST+SplitSolve approach has been tested on the Cray-XC30 Piz Daint at the Swiss National Supercomputing Centre (CSCS) [49] and on the Cray-XK7 Titan at Oak Ridge National Laboratory (ORNL) [50]. More details about both machines are given in Section 5A. Relevant is that they both offer hybrid nodes with CPUs and GPUs. The main features of FEAST+SplitSolve are then the following:

• •

The calculation of the OBCs and the solution of Eq. (5) are interleaved, FEAST being executed on the CPUs, while SplitSolve employs the GPUs;

• •

All the operations in FEAST, e.g. the solution of Eq. (6) rely on BLAS [51] and LAPACK [52], mostly zgemm, zggev, and zgesv;

• •

In SplitSolve, the matrix $A=(E\cdot S-H)$ is distributed over all the available GPUs and stored in their memory. Half of the matrix $Q$ is stored on the GPUs and half on the CPUs to reduce the memory consumption. The induced CPU$\leftrightarrow$GPU data transfer overlaps with computation (no cost);

• •

The calculation of each $Q_{i}$ block in Alg. 1 requires two matrix-matrix multiplications, one LU factorization, and one backward substitution. All the involved matrices are either dense or treated as such to leverage the cuBLAS (d/zgemm) [53] and MAGMA (zgesv_nopiv_gpu) [54] libraries.

The weak and strong scaling of FEAST+SplitSolve on Piz Daint are reported in Fig. 7 for a 2-D ultra-thin-body field-effect transistor (UTBFET), as in Fig. 1(c). All CPUs and GPUs per hybrid node are used. The calculation of the OBCs with FEAST is completely hidden by the solution of Eq. (5). The relatively low weak scaling efficiency originates from the extra calculation of spikes, which are required to parallelize the workload beyond 2 GPUs. These spikes are depicted in Fig. 6. Their generation takes 10 sec per recursive step so that the total simulation time increases from 30 sec on 2 GPUs (1 partition) up to 70 sec on 32 GPUs (16 partitions, 4 recursive steps). In the strong scaling case, the limitation comes from the impossibility to have a device structure that is large enough to create enough work on 16 GPUs and small enough to fit onto 2 GPUs.

Our strategy is to choose the minimum number of GPUs that can accommodate the desired nanostructure. This is what has been done in Fig. 8 on Titan where the time-to-solution of FEAST+SplitSolve is compared to FEAST+MUMPS [46] and shift-and-invert+MUMPS [38] for a 3-D Si nanowire field-effect transistor (NWFET) composed of 55488 atoms (16 GPUs) and a Si UTBFET with 23040 atoms (4 GPUs). MUMPS_5.0 has been selected because it is faster than SuperLU_dist [55] for these examples. The speedup between shift-and-invert+MUMPS and FEAST+SplitSolve is larger than 50 in both cases, i.e. our algorithms optimized for ab-initio transport calculations significantly outperform those designed for tight-binding problems. Also, SplitSolve alone is between 6 and 16 times faster than MUMPS on the same number of hybrid nodes. These speedups demonstrate that the scaling in Fig. 7 is not a limiting factor.

4. APPLICATION DESCRIPTION

The work flow of the OMEN+CP2K DFT-based quantum transport simulator is summarized in Fig. 2. It is basically a Schrödinger-Poisson solver with open boundary conditions to account for electron flows between contacts. CP2K, which is a freely available DFT package [56], is used here to construct the structure of the investigated nano-devices, relax their atom positions, and generate the corresponding Hamiltonian and overlap matrices. The latter are then transferred to OMEN, which performs quantum transport calculations based on them. This part consumes 99% of the total simulation time and is the focus of this work.

OMEN is a massively parallel, one-, two-, and three-dimensional quantum transport simulator that self-consistently solves the Schrödinger and Poisson equations in nanostructures. The tool was originally based on different flavors of the nearest-neighbor tight-binding model, but it has been extended towards ab-initio capabilities through its link with CP2K and the incorporation of the FEAST+SplitSolve approach presented here. It is important to realize that the OMEN version that produced all the results discussed in the next Section is the same as the one that is used to simulate nano-devices on a daily basis, either in the ballistic limit of transport or in the presence of scattering [57]. OMEN is a real production code and not a software that was only tuned to reach the highest possible performance on a specific instance of a general problem.

Although four natural levels of parallelism are available in OMEN, only three have been utilized in this paper, as shown in Fig. 9: the momentum $k$ and energy $E$ points are almost embarrassingly parallel, while FEAST+SplitSolve provides a 1-D spatial domain decomposition. The distribution of the workload is controlled with MPI [58] and a hierarchical organization of communicators. To avoid any work imbalance between sub-communicators corresponding to different $k$ points, a dynamical allocation of the number of nodes per momentum has been developed and verified before [45].

On the implementation side, OMEN is written in C++ and CP2K in Fortran 2003. The coupling between the two packages currently occurs through a transfer of binary files. Not all the nodes running OMEN load the Hamiltonian and overlap matrices, but only those necessary to gather all the unique parts of $H$ and $S$. The resulting data are then distributed to all the available MPI ranks with MPI_Bcast and copied to the GPUs. The parts of CP2K and OMEN ported to GPUs are implemented in the CUDA language from NVIDIA.

5. PERFORMANCE RESULTS

A. System Description

All the OMEN+CP2K simulations have been run on the Cray-XC30 Piz Daint at CSCS and Cray-XK7 Titan at ORNL. The technical specifications of both machines are summarized in Table I. The codes are compiled with GNU 4.8.2 and CUDA 5.5/6.5 and use the Cray MPICH 7.0.4, libsci 4.9, and MAGMA 1.6.2 libraries with customizations. On Piz Daint, all the CPUs per node are active, while at least half of them remain idle on Titan. The MAGMA function zgesv_nopiv_gpu is responsible for that: it performs a hybrid LU factorization on a CPU core and a GPU. On Titan it has been observed that the factorization time deteriorates if all the CPUs work because the competition between them negatively affects the execution of zgesv_nopiv_gpu. In spite of that SplitSolve is still about 10% slower per node on Titan than on Piz Daint.

B. Measurement Methodology

As can be seen in Algorithm 1 the number of floating point operations (FLOPs) involved in SplitSolve is deterministic and can be accurately estimated. Still, we decided to employ PAPI [59] to measure the CPU FLOPs and the CUPTI library [60] for the GPUs. In the CPU case, the function PAPI_start_counters with Events=PAPI_DP_OPS is placed right after the loading of the Hamiltonian and overlap matrices imported from CP2K, i.e. the input reading time is not accounted for in the performance measurement. Our limited access to the computational resources of Titan is the reason for this neglection: in a full-scale run, loading $H$ and $S$ and setting up the simulation environment takes about 4 minutes. Then each self-consistent Schrödinger-Poisson iteration lasts between 15 and 20 minutes. Knowing that an entire simulation involves roughly 40-50 iterations for 10 bias points, the 4 initial minutes turn negligible. However, due to our limited allocation, we could only simulate 1 or 2 iterations so that accounting for the initialization phase is not representative of the actual code performance. The function PAPI_stop_counters is set at the end of the simulation and includes the output writing time (a couple of seconds).

On the GPUs the FLOPs are obtained by sampling the device counters every second with the function cuptiEventGroupReadEvent. The accuracy of the method has been verified with examples where the number of FLOPs is known, e.g. matrix-matrix multiplications. This approach cannot be used “on-the-fly” since it considerably slows down the execution time. Hence, what is measured is the number of FLOPs that arise from the solution of the OBCs on CPUs and of Eq. (5) on GPUs for one single representative energy point. To find the total FLOPs this result is then multiplied by the total number of energy points that have been calculated.

C. Time-to-Solution

As explained earlier, ab-initio quantum transport simulations are usually limited to 1000 atoms due to their heavy computational burden. The FEAST+SplitSolve combination makes possible the investigation of much larger nanostructures such as a Si gate-all-around nanowire transistor with a diameter $d$=3.2 nm, a source and drain extension $L_{s}$=$L_{d}$=20 nm, a gate length $L_{g}$=64.3 nm, and a total of 55488 atoms. A schematic view of this device is given in Fig. 1(a), while its atomically resolved charge and current distributions as well as its spectral current are shown in Fig. 10 for one bias point. The structure dimensions are at least one order of magnitude larger than what is reported in the literature and comparable to what is fabricated in laboratories [61]. As indicated in Fig. 8 the computational time per energy point for this nanowire reduces to 102 sec with FEAST+SplitSolve using 16 hybrid nodes of Titan. Hence, it has been measured that a self-consistent Schrödinger-Poisson iteration including 2000 energy points takes less than 10 minutes on 8192 nodes. Knowing that the time per energy point with FEAST+MUMPS is in the order of 30 minutes on 16 nodes, a CPU machine with four times as many nodes would still be 3$\times$ slower than our refined approach.

D. Scalability

The weak and strong scalability of OMEN on Titan, after importing the required Hamiltonian and overlap matrices from CP2K, are presented in Fig. 11 and the data summarized in Tables II and III. A Si double-gate ultra-thin-body transistor, as in Fig. 1(c), with a body thickness $t_{body}$=5 nm, $L_{s}$=$L_{d}$=20 nm, $L_{g}$=38.2 nm, and a total of 23040 atoms has been chosen as test structure of realistic dimensions [62]. In each run, OMEN computed one Schrödinger-Poisson iteration for one single bias point with 21 $k$-points and FEAST+SplitSolve on 4 hybrid nodes. Increasing the number of bias points or self-consistent iterations does not change the scalability of the code since each point/iteration is processed sequentially, one after the other, and the workload is dynamically redistributed after each step [45].

In the weak scaling experiment the average number of energy points that each node deals with should ideally remain constant, but it varies here between 12.9 and 14.1, as detailed in Table II. These differences stem from the fact that the number of calculated energy points is not a direct input parameter, but depends on others, as explained in the caption of Fig. 11. If the extracted simulation times are normalized with the average number of energy point per node, as in the fourth column of Table II, it can be seen that the weak scaling efficiency is good with $\sim$5% variation across all the nodes.

Figure 11(b) shows the strong scaling results of OMEN. The total number of energy points is the same in all the runs (59908), but it varies with the momentum, i.e. $E$ depends on $k$. The simulation time decreases almost linearly when scaling from 756 up to 18564 nodes, reaching a parallel efficiency of 97%. This is not a surprise since the loop over the momentum and energy points is embarrassingly parallel, but it demonstrates that the code performs as expected.

E. Peak Performance

OMEN running FEAST+SplitSolve for the same UTB transistor as before initially reached a sustained performance of 12.8 PFlop/s in double-precision, as indicated in Table III. The solution of the OBCs and Eq. (5) for each of the 59908 computed energy points consumes 241 TFLOPs, 11 for the OBCs on the CPUs (5%) and 230 for SplitSolve on the GPUs (95%). The bottleneck is the LU factorization and backward substitution in Alg. 1. They are done with the MAGMA function zgesv_nopiv_gpu that does not perform as well on Titan as on Piz Daint. Replacing zgesv_nopiv_gpu with zhesv_nopiv_gpu and profiting from the property that the matrix $A=(E\cdot S-H)$ is Hermitian in 2-D structures (UTBFET) helps decrease SplitSolve’s execution time per energy point [63]. By combining this trick with further profiling and tuning of the code as well as algorithm adaptations to Titan, the double-precision sustained performance increased to 15.01 PFlop/s in production mode. Exactly the same structure as the one producing 12.8 PFlop/s was used for that purpose, but the simulation time decreased from 1130 down to 912.5 sec, while the number of operations per energy point went down from 241 to 228 TFLOPS.

The machine- and GPU-level power profiles of the simulation that reached 15 PFlop/s are plotted in Fig. 12(a). The 13 energy points that each group of 4 GPUs treats can be identified at both levels, but the GPU-level matches better with the different algorithm phases. The peak power consumption of Titan during this run is equal to 8.8 MW, its average to 7.6 MW, which corresponds to 1975 MFLOPS/W. At the GPU level, the computational efficiency increases to 5396 MFLOPS/W. Finally, in Fig. 12(b), the typical GPU utilization (measured with nvprof) is reported for the calculation of one energy point. Note the high utilization and well-overlapped computation-communication pattern.

6. CONCLUSION

By linking two existing applications, CP2K and OMEN, developing SplitSolve, an innovative, GPU-based, algorithm, and combining it with a parallel eigenvalue solver, FEAST, we have been able to redefine the limit of ab-initio quantum transport simulations. Significant improvements have been achieved in terms of structure dimensions (10$\times$ larger), sustained performance (15 PFlop/s), and time-to-solution ($>$50$\times$ shorter). While standard approaches are typically limited to 1000 atoms, we have reached more than 50000 for a 3-D nanowire structure and 23040 for a 2-D ultra-thin-body.

SplitSolve heavily relies on the structure of the matrices encountered in quantum transport calculations (block tri-diagonal + sparse right-hand-side) to deliver its best efficiency, but these properties can be found in other research fields such as computational fluid dynamics or in the solution of the Poisson equation. Hence, our multi-GPU sparse linear solver is not limited to one single problem, rather it may be applied to others.

From 2011 to 2015 the sustained performance of OMEN increased by one order of magnitude, going from 1.28 PFlop/s on the former Cray-XT5 Jaguar up to 15 PFlop/s on Titan, mainly through algorithmic and hardware improvements. A roofline analysis of SplitSolve and FEAST shows that both algorithms have high arithmetic intensity and are clearly compute bound. It can thus be expected that OMEN will run efficiently on future supercomputing systems offering lower relative memory bandwidth, but higher computational power such as the next-generation CORAL machines.

ACKNOWLEDGEMENT

This work was supported by SNF Grant No. PP00P2_133591, by the Hartmann Müller-Fonds on ETH-Research Grant ETH-34 12-1, by the Platform for Advanced Scientific Computing in Switzerland (ANSWERS), by the European Research Council under Grant Agreement No 335684-E-MOBILE, by the EU FP7 DEEPEN project, and by a grant from the Swiss National Supercomputing Centre under Project No. s579. This research also used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. The authors would like to thank Jack Wells, James Rogers, Don Maxwell, Scott Atchley, Devesh Tiwari, and colleagues at ORNL for giving them access to Titan and for greatly supporting their work on this machine, Ichitaro Yamazaki from UT Knoxville for his help with the MAGMA library, and Peter Messmer from NVIDIA for his GPU expertise.