RSDFT-NEGF transport simulations in realistic nanoscale transistors

The paper presents a device simulator for computing transport characteristics from first principles. The developed computer program effectively performs large-scale parallel calculation of quasi-one-dimensional quantum transport in realistic nanoscale devices with thousands of atoms in the cross section area of the device channel. Our simulator is based on the real-space Kohn–Sham Hamiltonian in the density functional theory and improved numerical algorithms for reducing computational burden in non-equilibrium Green’s function (NEGF) method. Several computational improvements have been introduced in constructing a reduced quantum transport model from the original Kohn-Sham Hamiltonian and implementing the R-matrix computational scheme in the NEGF simulations.


Introduction
Constant progress in the integrated circuit technology has been supported by the miniaturization of semiconductor silicon devices.Recent technological innovations make it possible to create semiconductor structures with the channel lengths of metal-oxide-semiconductor field-effect-transistors (MOSFETs) in the order of 10 nm or even smaller.However, the microelectronic scaling has been showing signs of reaching the end point due to production cost of chip fabrication, limitations in processing chemistry and physical properties of the materials.
Continued progress is no longer a matter of further reduction but may require novel approaches in science behind materials, device architectures and MOSFET technology for future semiconductor industry.The increase of the integration level has compelled the development of a variety of novel devices such as double-gate [1,2] and gate-all-around (GAA) transistors [3][4][5][6], FinFET [7][8][9], heterojunction tunneling devices [10][11][12], carbon nanotubes and nanoribbons.Fundamental physical limitations of silicon has also stimulated the search and development of innovative materials with higher carrier mobility and reduced sensitivity to high temperatures.New types of gate electrodes, various modifications of silicon dioxide and high-k dielectrics are being studied in connecting with applications in MOSFET technology.The modern growth mechanisms [13] and processing technologies enable one to synthesize composite materials with advanced material characteristics, low densities of defects and good control in uniformity.
In spite of this progress, continued downscaling meets with serious technological challenges as well as increased costs of fabricating and testing possible candidates for future electronic devices.In order to make an educated guess among many available options, it is desirable to build a proper device simulation environment that does not require any empirical data or parameters optimization.Detailed theoretical studies can help in better understanding the correlation of material properties and transport characteristics and addressing other practical issues which are needed to accelerate the development of MOSFETs with ultra short channels.
At the subnanometer scale, quantum effects play a crucial role in the device transport characteristics.Classical and semiclassical approaches fail and must be replaced by computationally much more intensive quantum-mechanical methods.Moreover, accurate transport modeling must capture details of the semiconductor band structure and dynamic properties of mobile carries in the presence of impurities, crystal imperfections or lattice mismatch at the channel interface.This stimulates the development of an accurate atomistic quantum device simulator based on the first-principle density functional theory (DFT) [14,15] for electronic-structure calculations.As of now, although various first-principle device simulators has been reported [16][17][18][19][20], there remain many practical issues regarding high computational cost and applicability to heterogeneous nanostructures which is essential for developing predictive tool and assisting the next-generation device development.
Significant efforts has been focused on developing effective methods for computing quantum transport in nanoscale devices [21][22][23][24].The quantum transmission boundary method [25][26][27][28] and the nonequilibrium Green's function (NEGF) formalism [29][30][31] are commonly used in quantum device simulations.At the device dimensions less than the dephasing length, ballistic transport is expected to dominate the device behavior.In this regime, the non-elastic processes can be neglected and the two methods are equivalent since the nonequilibrium state of the device can be described in terms of appropriate one-particle scattering wave functions.In quasi-one-dimensional transport, a device can be partitioned into small blocks in the transport direction and the Schrödinger or Dyson equations can be solved recursively without huge operations of the order of the entire device.However, the recursive algorithms requires complex-valued matrix inversion operations which may become numerically costly [32].On the contrary, the R-matrix method [33,34] only operates on real-valued matrices and offers more options for numerical optimization.
In this paper we present a device simulator which efficiently performs large-scale first-principle transport simulations in realistic systems with tens of thousands of atoms in the device channel.Our computational scheme is based on the real-space representation of the density functional theory which is suitable for parallel computations.In this scheme, singular ionic potentials are replaces by smooth non-local pseudopotentials V NL and the Kohn-Sham equations are treated as finite-dif- ference equations on N grid mesh points in the computational domain.Fast Fourier transformation in this scheme is unnecessary and due to the sparse nature of the Hamiltonian matrix the computational time for the main matrix operation H KS Ψ DFT scales as O(N grid ∕N CPU ) which makes it suitable for parallel simulations in realistic nanostructures.The real-space DFT (RSDFT) simulator [35] can perform atomic configuration optimization and electronic structure calculation in a supercell of up to 10,000 atoms with periodic boundary conditions.The self-consistent Bloch Hamiltonian H KS (q) contains the nonlocal terms ∼ q ∕ r and e iqr V NL e −iqr which in the real grid representation can be redefined as ∑ n W n exp(iqa n ) where q is the wave vector and a n is a translation vector in the Bravais lattice.A set of N grid ⊕N grid matrices W n can be used to define a device Hamiltonian and compute non-equilibrium electronic states with arbitrary boundary conditions.However, direct application of this approach to the NEGF simulations in realistic nanostructures are prohibitively difficult because of too heavy computational burden.For example, in a 10 nmdiameter Si nanowire channel with cutoff energy 21 E h ( E h is the Hartree energy), the size of the supercell is of the order of ∼ 10 7 .In order to improve the numerical accuracy of the finite difference approximation one has to deal with a long range kinetic energy matrix [36] which complicates the application of commonly used recursive Green's function algorithm.The R-matrix approach is much more suitable in such systems, since all the numerical operations in this method are real-valued and the recursion sequence can be chosen in arbitrary way thus giving full control over the size of the required inversion matrix operations.The propagation algorithm in the R-matrix method is constructed by adding an arbitrary small group of mesh points to a computation domain where the R-matrix is assumed to has already been calculated.The size of the R-matrix is determined by the number of mesh points with non-vanishing interaction with the rest of the device area.For a long range kinetic energy matrix, it is of the same order of magnitude as the whole supercell and the transport simulations may be prohibitively difficult in terms of both processing time and memory requirements.The heavy computational burden can be reduced by introducing a low rank basis representation for the Green's function.In the effective mass approximation, the well known mode space approach makes use of few lowest eigenstates (subband modes) computed at fixed value of coordinate along the current.The scattering problem in this representation is thus reduced to solving one-dimensional matrix equations in the transport direction.In atomistic simulations, the physically relevant transport modes (conducting electrons or holes) are located far from the edge of the total spectrum of atomistic Hamiltonian.Hence, even though one can always extract enough atomistic Bloch eigenstates to reproduce all relevant modes, there is no energy minimization principle to guarantee correct spectrum and energy band gap in the basistransformed transport Hamiltonian.This problem has been addressed in Ref. [37] by constructing a low-dimensional equivalent model (EM) for atomistic transport Hamiltonians.The method has been recently implemented and extended by several groups [19,20,[38][39][40] and the results have confirmed good accuracy of the reduced representation and applicability of the method to realistic nanodevices.The method has been often referred to as atomistic mode space approach.Here we prefer the original method name in order to indicate the difference with the ordinary selected eigenstates representation for bounded spectrum.The abbreviation EM should not be confused with the effective mass approximation, although the method plays a similar role, providing a simple approximation for the low-energy electrons in nanostructures.
The paper is divided into five sections.In Sect.2, we introduce necessary modifications to the EM scheme in order to incorporate the method into the real-space first-principle density functional program and obtain a reduced transport model suitable for large-scale parallel calculations.In Sect.3, we briefly discuss general properties of the R-matrix approach and demonstrate compatibility of the EM representation.In Sect.4, we discuss practical application of the R-matrix method in the EM representation and present test calculations in realistic nanostructures.In Sect.5, we discuss an alternative application of the EM approach which may offer an interesting option to facilitate numerical studies of quantum transport in more realistic inhomogeneous nanostructures.

Equivalent model for large-scale real-space DFT calculations
In this section, we discuss construction of the equivalent transport model for the RSDFT Hamiltonian using a Si nanowire as an example.The RSDFT simulator code [35] performs structural optimization and computes self-consistently an one-particle DFT Hamiltonian in a supercell with periodic boundary conditions in all three dimensions.A suitable device Hamiltonian can be defined by lifting the periodicity restriction along the wire and imposing zero boundary conditions at the boundary of the channel cross section.It is expected that these changes in the supercell geometry and boundary conditions do not affect drastically the transport and dielectric properties of the nanostructure.In the RSDFT simulator, the interaction with the nucleus and inert core electrons is described by pseudopotential [41] using the Fourier-filtering [42] method in order to obtain the nonlocal part of the pseudopotential in the form of short range projectors [43].Thus, one can redefine the unit cell by including a narrow vacuum region of ∼ 0.1 − 0.2 nm width outside the passivation hydrogen atoms and ignore the rest of the computational domain as long as the corresponding projector operators are fully present in the modified system.The energy band structure can now be found by solving the eigen-value problem for the Bloch Hamiltonian where q is the normalized wave vector in the transport direction, H 0 is the part of the Kohn-Sham Hamiltonian for an isolated supercell and the remaining part of the nonlocal interaction forms the coupling terms W(W T ) .In the (1) local-density approximation (LDA), only the non-local part of the pseudopotential for the boundary ions and the kinetic energy contributes to the coupling terms leading to the nearest neighbor interaction in the transport Hamiltonian.
The size of the RSDFT eigen-value problem can be rather large depending on the cutoff energy and the number of ions.In the device simulations, one is normally interested in a rather narrow energy spectrum of mobile carriers in the vicinity of the band gap which suggests using a projection type methods [44,45] for computing the relevant Bloch states Ψ q where numerates the energy states at given wave vector.The main numerical challenge in such methods comes from the complex-valued resolvent ∼ [z − H(q)] −1 at a set of z-points in complex-valued energy space and the computational time for M eigen states scales as ∼ N grid M N inv N iter , where the number of steps N iter does not strongly depend on the system size N grid .On the other hand, numerical tests show that the required number of iterations in the numerical inversion N inv grows faster than the size of the system and the cal- culations become rather time consuming in systems with hundreds of ions.Another way of solving the Kohn-Sham eigenvalue problem in electronic-structure calculations is a subspace iteration method [36,46].The algorithm is based on the recursive construction of a subspace spanned by a set of M ≪ N DFT low energy eigenstates.Given approxi- mate set of orbitals Ψ q , one applies the conjugate gradi- ent (CG) method to each member of this set in order to reduce the corresponding expectation energies (Rayleigh quotient) q ≡ ⟨Ψ q �H(q)�Ψ q ⟩∕⟨Ψ q �Ψ q ⟩ .Thus obtained set of states is used to construct a new ortho-normalized basis for the "improved" low-energy subspace.Projecting the original Hamiltonian onto this subspace and solving the corresponding low-dimensional eigenvalue problem generates a new set of orbitals which are used to start the next iteration.The most time-consuming part in this approach is due to subspace projection ( t ∼ N grid M 2 ).The computational cost can be reduced by regrouping the larger portion of the required operations into a matrixby-matrix product [36] and making use of a highly tuned linear algebra library (BLAS) [47].The subspace diagonalization (SD) can be done by using standard LAPACK (SCALAPACK) library [48] and the computational time for this step ∼ M 3 is relatively small in supercells up to few thousands atoms.The projection subspace is much larger in the CGSD method but the number of CG interactions should be small (typically ∼ 2 − 3 ) in order to preserve lin- ear independence of the orbitals which makes the method suitable for large systems.At the same time, in order to use the Bloch eigenstates in building the EM representation it is essential to ensure high numerical accuracy.Hence, in our simulations we combine these two methods.The CGSD method is used to generate all the states in the valence band as well as a low energy part of the conduction band.The numerical accuracy in the conduction band is further improved by applying the FEAST algorithm [45] for a narrow energy range relevant for transport.
Figures 1 and 2 presents a comparison between the band structures in two nanowires with diameters 1 nm and 8 nm for the Bloch Hamiltonian in Eq. ( 1) with zero boundary condition and the original periodic supercell.Our calculations confirm that different choice of the supercell geometry does not lead to any significant changes in the valence band spectra.A small positive shift of the conduction band spectrum can be seen in both cases but the energy difference in the band gap does not exceed 0.1 percent which is quite acceptable for our present purpose.
The new supercell geometry naturally defines the boundary of the conducting channel.In the device simulations, the Poisson equation is solved by assuming a continuous dielectric material in the gate region outside the RSDFT computational domain.More realistic description of the channel boundary requires incorporation of the area of dielectric into the RSDFT domain as discussed briefly in the following section.
The first step in building the EM representation for the RSDFT Hamiltonian is to construct an orthonormal basis Φ from a representative set of Bloch eigenstates within physically relevant energy interval [E 1 , E 2 ] .The number of the orbitals and the size of the preliminary basis N b is kept as small as possible but it has to be large enough to ensure that any Bloch state Ψ q with the energy (q) ∈ E 1 , E 2 can be represented by this basis for all q with enough accuracy.The coefficients in the basis representation are found by solving the corresponding eigenvalue problem for the N b × N b model Hamiltonian Hereafter, we employ matrix notations and omits the indices for grid points, basis states etc.To make a distinction with the original real-space grid representation all the quantities (Hamiltonian, self-energy, wave function etc.) in the basis representation are denoted using lowercase letters.
The eigenstates approximation provides a good representation for a bottom part of the spectrum of bounded operators.However, for an arbitrary energy interval the energy variation principle does not apply and the method fails due to appearance of spurious states with close energies.The false states can be identified by large ( ≥ 1 ) values of the residual factor which for the physical states is typically of the order of ∼ 10 −2 -10 −3 depending on the numerical accuracy in the CGSD scheme and the number of the primary representative RSDFT orbitals.In principle, an accurate low rank model can be constructed by simply adding extra electronic states from the valence band below E 1 .Such a direct construction can easily be done in the case of two-dimensional materials, (2) Ψ q = Φ q (3) h(q) q = (q) q , (4) h(q) = Φ T H(q)Φ.
(5) but it is impractical for realistic 3D nanostructures.The EM method [37]  In this case h (q) = (q) and X(q) represents the matrix elements Here we introduced a new vector where P = 1 − ΦΦ T is the projector operator to the orthogo- nal complement V ⟂ Φ of the subspace spanned by Φ.It has been shown [37] that the change of the band structure in Eq. ( 6) compared to the previous representation h(q) can be "measured" by the variational functional where the first term is a sum over a set of n q representative wavenumbers q i and ΔN(…) is defined as where and ⟨… ⟩ stands of for the average value over a set of points in the complex z-plane along the contour with center c = (E 1 + E 2 )∕2 and radius = (E 2 − E 1 )∕2 .ΔN(q, E 1 , E 2 ) represents a change in the number of energy levels at the wavenumber q within the ( 6) . Finding the minimum of the variational functional Eq. ( 9) is equivalent to constructing the model Hamiltonian Eq. ( 6) with less density of states.
The number of spurious states at each q within arbitrary energy range can change at most by one and the best solution is obtained when one unphysical state eliminated (shifted deep into the valence band) at all q i which corresponds to Dealing with the RSDFT Hamiltonian in realistic nanostructures requires massively parallel computations.The numerical performance strongly depends on computational details including the definition of orthogonal variational subspace and algorithms for choosing parameters and starting guess for effective minimization.In Ref. [37], we constructed Φ as a linear combination of the vectors (q) in which case the variational freedom is limited within a rather small ≤ 3N b dimensional orthogonal subspace spanned by all the functions with maximum coupling X (q) over the Brillouin zone.In this case, Eq. ( 9) becomes a simple analytical function of the expansion coefficients which is suitable for numerical applications.The reasoning behind this definition is that shifting the energy of a particular spurious state 0 at the boundaries of the Brillouin zone requires large value of the corresponding coupling terms X 0 in Eq. ( 6) which leads to a natural choice for starting (real-valued) guess Φ ∼ 0 and the variational subspace which incorporate all such functions.However, the RSDFT first principle calculations in systems with hundreds or thousands of atoms reveal strong dispersion of the spurious states especially in nanostructures with asymmetrical cross section (e.g.nanosheets).In such cases the number of spurious states grows and the previous simple choice often become insufficient for eliminating unphysical states at all q i .In this work we generalize the previous method by allowing for more variational freedom and formulate a reliable strategy for adjusting parameters in the variational calculations.In many cases this significantly reduces the size of the final EM basis and facilitates the transport simulations.
Let us consider a particular unphysical state 0 (q i ) in Eq. ( 3).In large systems, with rare possible exceptions, the diagonal term √ N grid which does not make any noticeable energy change Δ 0 .In this case ΔN(q, …) reduces to a small positive contribution from a sin- gle level located outside the targeted energy interval and the minimization of F[ Φ] fails to produce the desired solution.
In order to improve the numerical efficiency of the variational calculations we consider two types of states.The simplest candidates for Φ are the vectors with the maximum matrix elements X 0 at all q i .Introducing the real and imagi- nary components 0 = 1 + i 2 we find (13) where x 1 x 2 is the largest eigenvector of the 2 × 2 matrix: At the boundary of the Brillouin zone Eq. ( 13) reduces to the real-valued vectors (0) and ( ) used previously to construct an orthogonal subspace for variational calculations [37].One more trial vector of similar kind can be obtained by maximizing all the coupling terms at once.In a similar manner the solution is given by a linear combination of 2n q vectors 1 , 2 at all q i with the coefficients computed as the largest eigenvector for the corresponding 2n q × 2n q matrix.The second type of trial vectors can be obtained by estimating the maximum shift of the unphysical level Δ 0 .The condition det[ − h] = 0 gives and in the simplest approximation we set 0 ≈ 0 .We introduce the new vector 0 and find instead of Eq. ( 13) where x 1 x 2 is the largest eigenvector of the 2 × 2 matrix: Again, one might consider an additional vector of similar type by maximizing the energy shifts at all q i at once.In computing the auxiliary vector , we do not require high accuracy and use a simple estimate obtained by a moderate number of iterations in a standard CG scheme for real-valued symmetric matrix.Equations (13,16) give an initial set of trial states which are likely to have strong effect on the band structure within the targeted energy interval.Similar to Ref. [37] one can now introduce an auxiliary orthonormal basis for the Krylov subspace spanned by these vectors and obtain from Eq. ( 10) a rational variational function of the expansion coefficients.However, the corresponding basis is rather large and the variational freedom restriction is no longer justified.Instead, in the RSDFT simulations we use the original definition Eq. ( 9) and construct the new basis state Φ ∈ V ⟂ Φ without additional limitations.We also introduce a small modification of the variational functional by allowing for q-dependent energy target intervals in the definition of ⟨… ⟩ in Eq. (10).The calculations are most effective when each term ΔN(q i , …) in Eq. ( 9) returns a negative ( 14) value ≈ −1 which corresponds to shifting the unphysical energy at given q i outside the interval In a typical case of large H Φ Φ the energy shift is negative and much better results are obtained from the variational functional where at each q the targeting interval is placed above the spurious level 0 (q) which needs to be eliminated.Thus, we define 1 (q i ) = 0 (q i ) , 2 (q i ) = 0 (q i ) where E 1 <   0 (q i ) < E 2 and 0 (q i ) is the next spurious level in the reduce model Hamiltonian.If there are no unphysical states Minimization of the variational functional F[ Φ] is per- formed by the standard conjugate gradient algorithm [49,50].At each step the new conjugate directions is found as a a linear combination of the previous direction and the residual vector F∕ Φ calculated at the local directional minimum at the previous step.The residual vector is calculated in the form where the coefficients s and s are found from Eqs. (9)(10)(11).At each step one only need to perform simple RSDFT operations ∼ H Φ since all the other terms are linear com- binations of the previously computed RSDFT vectors with known coefficients.Most of the large matrix operations in the expansion coefficients and/or the line minimization are performed only once and the most time-consuming part of the numerical simulations is optimization of the initial trial states in Eq. ( 16).
Let us now summarize the main steps in the variational calculation.The first step is to compute the exact Bloch orbitals of the RSDFT Hamiltonian Eq. ( 1) at n q equi- distant wavenumbers 0 ≤ q i ≤ .We use the real and imaginary part of all the eigenstates with the energies within the targeted interval E 1 , E 2 to form a primary N b dimensional real-valued basis Φ [37].Compute the new set of the RSDFT vectors which defines a direct sum Φ 1 ≡ PH 0 Φ ⊕ PWΦ ⊕ PW T Φ which will be used further to form all the vectors (q) in Eq. ( 8).The variational cal- culations of the additional basis vectors in the EM model proceed as follows.
(2) From the obtained set of vectors select the trial state Φ which returns the minimum value of the original variational functional F[ Φ] , Eq. ( 9).Compute the new set of three RSDFT vectors � (3) For each q i solve the reduced eigenvalue problem in the basis representation Φ ⊕ � Φ and redefine 0 (q i ) as well as the next unphysical level 0 (q i ) .Define the vari- ational functional in Eq. ( 18) by setting 1 (q i ) = 0 (q i ) and 2 (q i ) = 0 (q i ) .If there is no unphysical level (4) Minimize the variational functional F[ Φ] by the CG algorithm as outlined above.
Repeat the steps ( 2) and ( 3) to reach the minimum value of 1) -( 4) until there are no more spurious states at all q i .As an auxiliary useful options one may add the band gap verification at step (2).This is done by solving the eigenvalue problem for computing the self-energy (see below) in the model Φ ⊕ � Φ below the conduction band edge and counting the number of open channels.Then the selection of the trial state in (2) is to be made among the candidates with the smallest number of open channels.This option helps to ensure that steep unphysical branched are not missed by the simulator even if the q-grid is not dense enough.
Figures 3, 4 show two examples of the EM construction in Si wires with diameters 1 nm and 8 nm.The black points in these figures represent the exact eigenvalues of the Bloch Hamiltonians at n q = 13 representative wavenumbers.The left (right) panel show the band structure in the reduced representation before (after) the variational calculations.In these two examples, the correct transport models are obtained after 7 and 64 iteration steps reducing the size of the scattering problem by the factor of N b ∕N grid ≈ 5 × 10 −3 and 4 × 10 −4 respectively.

R-matrix method for ballistic quantum transport
The atomistic first-principle transport simulations require the full quantum description of non-equilibrium electronic state of the device.In the two-probe configuration the system under consideration consists of a central area (device) and two semi-infinite electrodes (leads) which are periodic in the transport direction.The part of the device at the contacts with the leads is assumed to be a seamless extension of the atomistic structure of the corresponding electrodes.In the NEGF formalism, the electric current and density are calculated from the Green's functions which need to be found from the following system of equations where H D represents the Kohn-Sham Hamiltonian for the electrons in the device area and Σ R,< (E) are the self-energy terms which describe possible scattering mechanisms.H D also includes a potential term which depends on the applied bias and mobile density distribution which need to be computed self-consistently.Here we only consider ballistic transport in which case the retarded self-energy Σ R comes from the outgoing boundary conditions at the contacts and the lesser self-energy term Σ < describes the flow of mobile carriers from the corresponding reservoirs.In the ballistic regime, the Green's functions are computed independently for each energy value and Σ R,< are spatially localized in the contact region which greatly reduces computational cost.A current-carrying device state can be computed in terms of the boundary part of the retarded Green's function which is equivalent to solving an one-particle scattering problem.In the real space representation, the electronic states are calculated on a three-dimensional spatial grid and the electronic Hamiltonian is represented by a sparse (or block tridiagonal) matrix.The transport simulations are often conducted by the recursive Green's function algorithm [51,52] or quantum transmitting boundary method [25].The RSDFT Hamiltonian generally has a long range nondiagonal part and the recursive algorithms require large complex-valued matrix inversion which yield prohibitively large computational loads.One the other hand, the R-matrix approach [34] only uses real-valued arrays and enables one to optimize the size of the inversion operations.
The retarded Green's function in Eq. ( 20) can be calculated in terms of the corresponding Green's function G 0 in the close system without leads.In the ballistic simulations, one only needs a small boundary part of the Green's function and all the results can be obtained from the realvalued symmetric R-matrix where P is the boundary projector to the contact region of the device D with non-vanishing Hamiltonian matrix elements H DL where L indicates the grid points outside the device area.The symbol P should not be confused with the projection operator in the previous section where we studied trial functions orthogonal to the EM space.
For example, the drain current is given by the Landauer formula [29] (in a.u.)where ( 22) the scattering rate in the p-th probe and f p is the correspond- ing Fermi factor.
From the boundary projection of the Dyson's equation one finds and, since the self-energy is spatially localized Σ R = Σ R P = PΣ R , the transmission coefficients Eq. ( 24) can be computed in terms of the R-matrix.Although the real-valued Hamiltonian matrix H D in Eq. ( 22) is often very large, the R-matrix can be calculated recursively, like in other popular computational schemes [34].For completeness, we outline here the R-matrix propagation scheme in the most general form.
One may consider an arbitrary fragment C ⊂ D of the device area and define the corresponding local R-matrix R cc by the equation similar to Eq. ( 22) with P, H D being substituted by the corresponding operators P c , H c .One next adds a small area a ≪ C from the rest of the device and calculate the new R-matrix as the boundary projection R = (P c + P a ) G(P c + P a ) of the Green's function for the closed system C = C + a The calculations are straightforward.Two blocks Gcc and Gac in Eq. ( 26) are calculated in terms of matrix products which contain the Green's function G cc = [E − H c ] −1 .Tak- ing the boundary projection turns G cc into R cc and gives the final result Equation ( 27) is just the R-matrix propagation recipe Eq. ( 19) in Ref. [34] written in a more compact form.In the above equations, the main Hamiltonian matrix operation X ac ≡ H ac R cc needs to be performed only once.In practical application one can keep the size of a small and use columnlike distribution of the R-matrix among the processors.Computing X ac scales as ∼ N a N c ∕N CPU with a proportionality fac- tor depending on the RSDFT parameters.Computing the total correction X ac H ca in Eq. (27a) and Rca in Eq. (27b) also requires inter-mode communications which scale in a similar way and does not add computational burden in practice.The most time consuming part of the simulations is the realvalued rectangular matrix multiplication in Eq. (27c) Rca X ac (25 which in parallel computations scales as N a N 2 c ∕N CPU .This can be effectively performed even for a large R-matrix by making use of a highly tuned linear algebra library (BLAS) [47].Finally, one redefines the boundary P c + P a → P c and reduce R by eliminating the rows and columns for the redun- dant points ∉ P c .The propagation start from the R-matrix R aa = E − H a −1 for an arbitrary initial fragment.
After the full R-matrix is computed one can add the contact self-energies and obtain the boundary part of the retarded Green's function Eq. ( 25).For the p-th probe of nanowire the corresponding self-energy term is calculated as [53] where {�� ⃗  p , � ⃗ Z p } is the outgoing/decaying solution of the generalized scattering eigenvalue problem [24] in the p-th probe.Even in the simplest cases, such as the one in Fig. 1, solving this problem in the original real-space grid representation is extremely difficult due to the large size of the supercell [54].On the contrary, there is no difficulty in computing the 2 N b possible scattering states in the EM representation.The choice of the physical characteristics of the probes outside the conducting channel is a matter of convenience.The only essential condition is that the difference in the materials does not cause noticeable reflection from the contact area.The probes with the low rank EM Hamiltonian satisfies this condition by the construction.As a test, we compute the transmission function in a 1 nm-diameter Si nanowire (SiNW) in Fig. 1.Adding the EM contacts can be understood as an extra step of the R-matrix propagation.On can consider a final fragment as a composite of the unconnected unit cells at two contacts and Eq.(27a) gives the 2 × 2 block R -matrix in the EM representation where h i and p i stand for the unit cell EM Hamiltonian and the corresponding projector at the EM contact 1 (left) or 2 (right).The transfer part of the Hamiltonian at the contacts with the leads is determined by the one-sided EM basis transformation H DL p 2 = WΦ , H DL p 1 = W T Φ and their transpose.The transmission function can now be calculated without referring to the original real space grid representation.Figure 5 illustrates the outlines procedure.The computed transmission function exhibits a step-function like behavior which clearly indicates that the contacts between the original first principle model and the probes in the EM representation cause no unphysical backscattering.
It has also been shown that the R-matrix provides all the necessary information for computing the diagonal part of (28) the lesser Green's function G < .The boundary projection of the retarded Green's function and the sequence of the local R-matrix elements Rac , Raa in Eq. ( 27) can be used to con- struct any scattering solution in the reverse order.Thus, one can calculate the mobile charge distribution and obtain self-consistently the electric potential in the device area.Application of the R-matrix method to the self-consistent device simulations is discussed briefly in the next section.
The above matching procedure is equivalent to computing the retarded Green's function in a RSDFT + EM composite system with boundary unit cells being treated in the EM representation.A unit cell at the contact corresponds to the diagonal block h + R in the composite Hamiltonian matrix coupled to the rest of the system by the off-diagonal blocks WΦ and Φ T W T .Calculating the inverse in Eq. ( 20), for the part of the system in the original real-space representation one obtains the Green's function in the same form as Eq. ( 20) with the RSDFT self-energy correction Thus, one can use the EM representation only as a tool to compute the self-energy terms and treat the entire simulation domain using the original DFT basis.Note that Eq. ( 30) differs from the direct transformation ΣR = Φ R Φ T suggested in Ref. [19].Using the exact identity R = w(E − h − R ) −1 w T , which is valid for any self-energy in the form of Eq. ( 28), we obtain

R-matrix method in quasi-one-dimensional systems
In realistic nanostructures, computing the R-matrix in the original real-space grid representation may become prohibitively time consuming.Although the size of the local R-matrix remains nearly constant in the course of the propagation, the boundary of arbitrary domain is determined by the range of the nonlocal interaction in the RSDFT Hamiltonian.Within the local density approximation (LDA), the decisive factor is the highest order in the finite-difference expression for the kinetic energy operator [36] and for a typical choice of parameters in the RSDFT geometry optimization simulations the kinetic energy has a long range up to 6 a.u.As a result, the size of the boundary region is of the same order of magnitude as the supercell of the wire.For example, the supercell in a 8 nm-diameter ⟨100⟩ SiNW in Fig. 2 with the cutoff energy ∼ 20 E h contains more than N grid ∼ 1.1 × 10 6 grid points.For a 10 nm channel length, the total number of the multiplication operations in the R-matrix propagation is huge ∼ 2 × 10 19 which makes it next to impossible to perform actual self-consistent device simulations.
In homogeneous nanostructures, one considers a device Hamiltonian in the form where n numerate unit structures (supercells) in homogeneous wire and V n is an extra diagonal potential term which needs to be computed self-consistently with respect to the mobile charge distribution in the device area.In this case the EM basis representation can be used through the whole device structure.The corresponding low-dimensional device Hamiltonian is obtained by the EM basis transformation Eq. ( 4) which reduces the simulation time by the factor of ∼ (N grid ∕N b ) 3 and enable device research.Thus, the EM model in right panel in Fig. 4 contains N b = 367 real-valued basis functions which reproduce the band structure and the corresponding microscopic solution with 2 -3 significant digits of accuracy at the bottom of the conduction band.The effective energy range in this representation ∼ 0.4 eV is more than enough for accurate ballistic transport simulations [37].
We consider a n-SiNW GAA MOSFET shown on the inset of Fig. 6 with a gate oxide thickness t ox = 2 nm .In these simulations, the oxide layer is considered to be a continuous media with SiO 2 = 3.8 which only affects the bound- ary conditions in the nonlinear Poisson equations.The validity of this approximation is briefly discussed below.The whole simulation domain is 30 nm in the transport direction which corresponds to M = 55 supercells.Other parameters are: V SD = 0.1 V , T = 300 K , Si = 11.9 , dopant concentra- tion in the source/drain regions 10 20 cm −3 .The first (source) (32) where �� ⃗  1,2 , ⃗ z 1,2 are the matrices of the outgoing/decaying Bloch states and the Bloch factors in the corresponding (left or right) leads [53].
The NEGF equations in the EM representation do not involve any large quantities and can be solved easily.M supercells in the device are used as independent blocks of the R-matrix propagation.As become clear below, the propagation sequence must be chosen such that the 1-st and M-th contact blocks are added last.For example, one can start with the supercell n = 2 and add the rest of the blocks from the right one by one.Adding the last block n = 1 at the source contact completes the propagation sequence.An intermediate cluster in the propagation algorithms is a sequence of n − 1 unit cells 2, 3, … , n .The boundary of the cluster consists of two unit cells and the R-matrix is defined as a 2 × 2 N b -dimensional block matrix r 11 r 12 r 21 r 22 , analogous to the EM-transformed R-matrix in Eq. ( 29).After adding the next (n + 1)-th block from the right, the R-matrix is redefined as which is a particular case of the previous general result Eq. ( 27) written in a form of recurrence relations.The initial conditions at the first step are given by . The final step of the calculations is to add the contact unit cell n = 1 .The R-matrix transformation is given by Eq. ( 34) with h n+1 being replaced by h 1 and opposite propagation direction which corresponds to the transformation 1 ↔ 2 and w ↔ w T .The final R-matrix and the contact self energies Eq. ( 33) determines the boundary projection of the retarded Green's function in the form of Eq. ( 25).The drain current is calculated directly from the low-dimensional model without referring to the original real-space representation.On the contrary, the electronic density distribution on the original three dimensional grid in the computation domain is essential for calculating the electric potential in the device.In realistic systems, as the density of states grows, computing the electronic density becomes the most time-consuming part of the device simulations which is deserving of more consideration.

Optimization of the carrier density calculations
The electronic density in the NEGF formalism is determined by the diagonal part of the lesser Green's function and can be evaluated as [29] where p is the EM scattering rate at the p-th contact.In practical implementations, the required number of operations grows as ∼ N grid N b which becomes the major challenge for large systems.The R-matrix formalism enables one to significantly reduce the computational cost.We consider partial scattering eigenstates of the EM contact scattering operator with positive Λ p .Apart from rare exceptions of energy values located close to a scattering threshold in the band structure, the spectrum of Λ s has a finite gap and we use a simple cutoff criteria Λ > Λ 0 ≈ 10 −3 a.u. to identify the scattering states.The number of these open channels N ch is typically small compared to the EM dimensions N ch ≪ N b .Equation ( 35) can now be rewritten in the form [34] which is similar to the quantum transmission-boundary method [25].Here the wave functions Ψ p are defined as which is different from ordinary scattering modes but it is particularly useful in the backward wave function calculations.Equation ( 38) satisfies the Schrödinger equation in the EM representation everywhere except at the contact unit cells which explains the choice of the R-matrix propagation scheme.The sequence of the local R-matrices in Eq. ( 34) contains all the necessary information to construct all scattering solutions at given energy.To see that, we calculate a boundary projection of the Lippmann-Schwinger equation in a sequence of blocks 2, 3, … , n to obtain for the wave function in the last added block.Equation (39) gives a numerically stable double-end recurrence relation for the backward wave function propagation.A set of the connection matrices r 21 w T and r 22 w must be stored during the forward recursion Eq. (34).In turn, the boundary projection of the EM partial solution p in Eq. (38) gives the boundary condition and the wave function in all other blocks is calculated from Eq. ( 39) in the reverse order.Since the R-matrix refers to a close system without probes, the backward connection matrices do not depend on the boundary conditions which is important for practical implementation of the method.In realistic systems, one needs to compute contribution to the mobile charge at huge number of grid points from many partial scattering solution which becomes computationally costly.The use of the EM representation in the channelindependent backward propagation can significantly reduce the computational load.At each propagation step n → (n − 1) the contribution to the electronic charge from all the partial states is computed by performing three matrix operations = T ; Ψ = Φ ; = ΨΨ † diag , where Ψ is a complex- valued (N grid ∕N CPU ) × N ch array of RSDFT scattering solu- tions, is the corresponding array in the EM picture and T represents the real-valued connection matrices in Eq. (39).Transforming the density calculation into matrix-by-matrix products and making use of a highly tuned linear algebra library (BLAS) [47] greatly reduces the computational time.As an illustration, we show in Fig. 7 the CPU time for ( 40) computing the integrant Eq. ( 37) as a function of total energy.The reference time t 0 corresponds to the forward R-matrix propagation.For comparison, we also show the simulation time in performing the EM-RSDFT basis transformation for each partial scattering solution before computing its contribution to the electronic density (black curve).The increase in the computational time is due to the larger number of partial contributions in Eq. ( 37) which is up to ∼ 140 t 0 at the highest energy values.Although the total number of floating point operations is the same in both cases, regrouping Eq. ( 37) into the matrix form suitable for the BLAS matrix operations makes the simulations almost two orders of magnitude faster.An important feature in the original RSDFT simulator is an efficient real-space parallelization [36].Because of the modified supercell geometry the message passing interface (MPI)-communication scheme in the device simulator differs from the original code.The supercell is split into a set of clusters keeping the number of mesh points within a cluster small and nearly constant.The same fragmentation scheme is used both in the exact R-matrix propagation scheme of the previous section and for the global numeration of the mesh points though the supercell.In general, the maximum number of points in a cluster is a numerical parameter which should be optimized based on the CPU time for the matrix multiplication and the inversion operation in Eq. ( 27).In the large-scale NEGF EM device simulations we choose this parameter based on the condition that the number of clusters in the supercell is of the order of the number of nodes.The numeration of the mesh points through the clusters is chosen is such a way as to ensure geometrical proximity of the points with consecutive numbers in order to reduce the number of relevant nodes required for the MPI-communication in the non-local part of the Hamiltonian operator.The optimization of the MPI-communication scheme in the modified geometry has not been fully addressed yet and the non-local operations on the real grid are presently about twice as slow as in the original code.
The Poisson equation is solved on the same three-dimensional equidistant grid used in the RSDFT structure calculations.A dielectric layer in the gate area is introduced by extending the original grid beyond the boundary of the modified supercell.We thus obtain a layer of mesh points with a desired thickness t ox which is used to represent a die- lectric layer in the device structure.The non-linear Poisson equation is solved in the integral form by using the Newton-Raphson method.panel also shows the one-dimensional profile of the electric potential averaged over the cross section of the silicon channel.The number of atoms in the simulation domain in this case N at = 83, 875 and, although the total number of mesh points exceeds 6.5 × 10 7 , the self-consistent device simulations can be performed in a reasonable timeframe due to the EM basis representation and the properties of the R-matrix method.For a given potential distribution, the energy-dependent R-matrix propagation is performed independently at energy quadrature points evenly distributed among available nodes.The obtained local R-matrix data is sent back to all the nodes and independent cumulative contributions to the electronic density from all the scattering states are computed as a single matrix-by-matrix multiplication operation ∼ (N grid ∕N CPU ) N b N ch .The device simulations in this section are performed using N CPU = 1, 536 CPUs of the Fugaku supercomputer.Computing the I-V characteristics in Fig. 6 take ∼ 20 hours in total.This does not includes the preliminary RSDFT first-principle structure optimization of the silicon crystal structure and construction of the EM basis representation.

Numerical illustration: RSDFT-EM transport simulations to silicon FETs
The developed simulator is applicable to any geometry of the conducting channel.The only difference may come from the boundary conditions in the Poisson solver as well as the choice of the supercell which is defined by the original RSDFT input file containing the atomic coordinates.As an illustration, we outline the main steps in the device simulations for a n-Ge nanosheet FET with a cross-section of 4 nm by 8 nm.The original supercell of the Ge diamond structure in [001] direction contains 999 atoms including 172 hydrogen atoms necessary to fully passivate the surface of the sample.The equidistant grid of 968,000 mesh points in the supercell has been used to perform the RSDFT structural optimization of the atomistic structure which corresponds to a cutoff energy ∼ 20 E h .These simulations takes ∼ 60 hours.One thus obtains the optimized atomic configuration and the self-consistent local potential which fully specifies the one-particle Kohn-Sham Hamiltonian and the corresponding energy band structure.As a test, we repeat the self-consistent calculation of the electronic density/potential at higher cutoff energy ∼ 30 E h using the previously obtained atomic struc- ture.Apart from an immaterial energy shift, the conduction band structure is found nearly identical within ∼ 0.01 − 0.02 eV level of accuracy.The device simulations are performed using the smaller RSDFT data set.The geometry of the unit cell is modified in order to impose zero boundary conditions outside the channel cross section which includes a narrow 0.2 nm region outside the hydrogen atoms.This gives a smaller supercell with N grid ∼ 8 × 10 5 mesh points.Figure 8 shows a part of the energy band structure at the bottom of the conduction band.The black scattered points represents the Block states at a set of 13 equidistant wave numbers which are used to form an initial scattering eigenstate basis.Compared to the previous example of the silicon wire, this figure indicates a stronger q-dependence on the Bloch state across the Brillouin zone.The eigenstate basis transformation produces a large number of spurious branches with strong dispersion.In this case the second type of the trial basis states Eq. ( 16) play an important role in reducing the number of steps in the variational simulations.The red solid lines in Fig. 8 shows the band structure in the final EM model with N b = 389 basis states.The variational cal- culations take ∼ 14 hours using 1,536 CPUs of the Fugaku supercomputer which includes computing the exact Bloch states and 56 supplementary basis functions The device Hamiltonian is obtained by performing the EM basis transformation of the chain Hamiltonian Eq. ( 32).The electric potential on the three dimensional grid is found from the Poisson equation assuming extra homogeneous dopant concentration N + = 10 20 cm −3 in source and drain regions both consisting of 16 supercells.The most time consuming part of the simulations is the density calculation in Eq. (35).In order to reduce the number of energy points we use Gauss-Jacobi quadrature [0, −1∕2] in a sequence of narrow intervals between the ordered threshold energies in two probes.The lowest threshold corresponds to the bottom of the conduction band where N ch changes from zero to a positive number.The next thresholds are determined as the energy values where N ch changes its value.We do not Fig. 8 (Color online) Conduction band in a nanosheet n-Ge FET with a cross-section of 4 nm by 8 nm.The black squares shows the results of RSDFT simulation with 9 × 10 5 grid points in a supercell.The red lines represent EM with 389 basis functions distinguish close threshold energies ΔE thr < 10 −3 eV which are treated as (nearly) degenerate branches.The number of quadrature points in each interval can be kept small ∼ 20 without losing numerical accuracy.Equidistant grid is used for the extra energy interval −0.3 eV above the highest threshold energy.
Optimizations of the density calculation has been discussed in the previous example.Another practically useful scheme is based on separation of the constant basis in Eq. ( 35) from the energy integration.In this case, the backward propagation is used to compute partial contributions to the N b × N b diagonal blocks of the density matrix in the EM representation D ∼ ∑ Ep f p (E) p (E)⊕ T p (E) which only needs an additional inter-node MPI_ALLREDUCE communications.The real space density is then obtained by a single EM transformation ∼ [Φ T DΦ] diag which scales as (N grid ∕N CPU )N 2 b .Figure 9 shows I-V characteristic in two n-Ge nanosheeet FETs with different gate length.As expected, the device with larger gate area L g exhibits far better gate control but shows similar value of the drain current in ON state.This can also be clearly seen from the electric potential shown in the right panel of Fig. 9.The total simulation time for the I-V characteristics is ∼ 18 -23 hours depending on the numerical parameters and integration scheme which is quite suitable for a practical device research.

Si-SiO 2 interface in the RSDFT-EM transport simulations
In the previous examples, the material properties of the dielectric layer in the gate area have not been considered and the DFT calculations are performed in the semiconductor area only.This area however cannot be defined uniquely as there is certain freedom in how one specifies the boundary of the semiconductor domain.The previous choice has been made in order to keep the computational domain as small as possible without losing any pseudopotential terms.
The question remains of how different choice of the boundary may affect the electric potential within the channel or whether actual mobile density distribution at the boundary semiconductor/dielectric may have a noticeable effect in the self-consistent device simulations.
As the last example we consider an oxidized n-SiNW channel by including the silicon dioxide in the first-principle RSDFT structure simulations.Even for a thin oxidized wire the DFT simulations are much more demanding compared to the previous examples.Due to a shorter O-H bond in silicon dioxide the RSDFT structural optimization generally requires a much denser three dimensional grid.The energy band structure also contains many dioxide bands at low energies separated from the valence band by a large gap ∼ 7 eV which make the subspace diagonalization in the CGSD calculations more demanding.Yet another obstacle comes from a trial atomistic configuration which needs to be specified in order to start the geometry optimization.The simplest trial configuration can be obtained by inserting oxygen atoms in the mid positions in the silicon lattice outside a specified radius.However, numerical tests indicate that such trial configurations are totally useless.Large lattice mismatch at the Si-SiO 2 interface in the wire geometry causes strong local strain in the boundary region and produces multiple spurious energy bands at the Fermi level leading to a very poor convergence in the structural optimization.In order to reduce the strain and insure the finite band gap in the band structure, we perform multiple steps of oxidation-optimization process.Figure 10 illustrates a few steps in such calculations.The simulations start by selecting a 2 nm-diameter cylindrical core in a larger silicon sample with perfect diamond structure.In order to reduce the local lattice mismatch one also has to increase the size of the supercell in the wire direction.For a [001] nanowire, we include eight atomic layers along z-direction in the definition of the initial silicon supercell.At each following step we introduce a small number of oxygen atoms into vacant midpoint positions between silicon atoms outside the core region and perform the RSDFT optimization to obtain a stable sample with converged total energy.Even for such a thin sample, the size of the supercell is as large as 3.9 × 10 6 which corresponds to a cutoff energy ∼ 69 E h .The simulation for one oxidation step with 50 extra oxygen atoms take about one day.The simulations proceed until there are no more vacant positions.The final nanostructure consists of 226 Si atoms, 174 O atoms and 104 H atoms to fully passivate the surface.The structure corresponds to a 2 nm-diameter Si wire with a 1 nm-thick silicon dioxide (SiO 2 ) layer. Figure 10b shows the supercell in the final oxidized Si-SiO 2 nanowire as well as the core region of strained silicon which can be used for comparison.The core area is defined based on the atomic positions of the nearest oxygen atoms and the corresponding bonds are passivated by attaching a new set of hydrogen atoms.We further reduce the size of the supercell by extracting four atomic layers in [001] direction and thus obtain a supercell of 69 Si atoms which can be used in simplified transport simulations in a nanowire with "equivalent" central silicon atomic structure.
Figure 11(a) shows the energy band structure in the Si-SiO 2 supercell (called "Si-SiO_2" hereafter) and in the equivalent silicon wire.In the latter case we use 16 CPUs of Inter 2.6 GHz workstation to compute the local RSDFT potential and perform all the necessary steps in device simulations.Figure 11 also shows the band structure in three Si supercells with the same number of atoms which we use for comparison.Case (b) is obtained by reducing the size of the supercell by extracting four atomic layers in [001] direction ("Str Si").Case (c) is obtained by relaxing the strain in the extracted central core region ("Opt Si").Case (d) corresponds to relocating the silicon atoms into the nearest position in an ideal diamond lattice ("Si").The main feature of the EM method is that the size of the reduced device Hamiltonian depends only on complexity of the energy band structure and not on the original exact representation.Although the supercell of the oxidized sample is a few times larger compared to the previous examples, the low energy part of the conduction band structure is as simple as the one in thin Si wires which make the EM construction rather trivial.The primary basis of the Bloch eigenstates already generate a reduced model with correct band gap and there are just three spurious branches within the transport energy range which are easily eliminated by three extra steps.We thus obtains a 61-dimensional EM representation (N b ∕N grid ∼ 2 × 10 −5 ) suitable for device simulations.We have computed the I-V characteristics of a 2 nm-diameter n-SiNW FET using the constructed RSDFT-EM model for the oxidized silicon with 1 nm-thick SiO 2 layer.Other parameters are the same as in the previous examples.The results are compared with simple models using continuous dielectric media in the gate region.Our simulations show that the drain current in the ON states have similar values but the subthreshold slope in the continuous dioxide models is significantly underestimated.The likely reason for this discrepancy is the difference in the charge distribution and the boundary condition in the area close to the gate contact.The potential profile for two simulation models in Fig. 12 clearly indicates different width of the effective potential barrier.In order to reduce the discrepancy, we have changed the gate geometry in continuous modeling by extending the gate area by 1 nm.The right panel in Fig. 12 shows the potential profile in the new configuration and Fig. 13 present our final results.The calculated drain current in the strained Si channel with equivalent atomistic configuration is in good agreement with the original Si-SO 2 modeling.Introducing extra strain relaxation in the atomic structure leads to the increased drain current in ON state, but show similar subthreshold slope.Despite generally fairly good agreement between these models, both the band structure and the electronic charge distribution in the channel are quite different.Figure 14 shows the electronic density at two representative values of gate voltage in different models for the same Si channel.The strain relaxation leads to a more oriented atomic structure with reduced density fluctuations across the supercell.In idealistic homogeneous nanostructures this does not seem to lead to significant changes in the simulation results.However, in a more realistic devices with a spatially inhomogeneous Si/SiO 2 interface the full ab initio description may be necessary in order to properly reproduce the device characteristics.Unfortunately, the present EM method is not well suited for inhomogeneous materials and the computational strategy needs to be modified.In the next section we briefly discuss a possible way to overcome the limitations of the EM method.

RSDFT-EM representation in scattering-matrix method
The EM method originally assumes that a nonequilibrium state of homogeneous device is formed as a result of various scattering processed in a reference periodic nanostructure.Thus, the same set of basis functions can be used to represent transport properties of all part of the device channel.The method is suitable for inelastic transport simulations since reducing the size of the device Hamiltonian enables one to compute the full Green's function in the basis representation [37,38].Scattering by lattice imperfections or impurities, treating heterogeneous devices of varying cross section and/or materials requires different approach.In principle, physically relevant modes can be extracted locally from a set of scattering problems in a sequence of larger supercells along the channels.Therefore it may be possible to represent defects or impurities by a local basis.However, such a local basis representation may produces extra backscattering at the boundaries between different representations thus leading to enhanced localization and underestimated transmission coefficients [55].On the other hand, atomistic simulations in a tight binding model show that physical solutions in the entire channel can be well reproduced even if the basis representation does not completely eliminate the unphysical backscattering [56].This suggests that the EM representation can be used as an approximate solution to facilitate more time consuming atomistic simulations.Since the regular perturbation expansion may not be suitable for capturing the behavior of a nonequilibrium open system, we consider a more straightforward approach and use the EM approximation as a preconditioner in recursive solution of inhomogeneous scattering equation.The partial contribution to the electric charge in Eq. ( 37) is defined in terms of the eigenstates of the scattering operator which ensures correct normalization of the corresponding wave function.In this section we prefer the standard definition of scattering wave functions which correspond to the incoming wave in particular channels.With proper normalization, one can find these solutions approximately and compute the carrier density and drain current.We consider a structure in Fig. 5 with the device Hamiltonian H D and two probes in the EM representation.The scattering mode Ψ is found from the equation where ⃖���   is the incoming wave in open channel with the corresponding phase factor As usual, we assume semi-infinite probes of nanowire geometry and square matrices of Bloch eigenstates �� ⃗  and ⃖��  [24].To keep the nota- tions simple, we use a single index to numerate possible modes in both probes.The second term in the brackets in the right hand sides represents outdoing/decaying waves in the probes after scattering.The unknown expansion coefficients a are related to the elements of the S-matrix and therefore to the transmission coefficient.The unitarity of the S-matrix requires unit-current normalization of the Bloch eigenstates where v is the group velocity and � ⃗ Z  is the corresponding Bloch factor [24].Under this condition, the mobile charge density is calculated from Eq. ( 37) and the transmission function is given by where the prime indicates that and run over open channels in different leads.The corresponding elements of the scattering matrix are obtained by projecting the wave function in the opposite EM contact to the outgoing states S  = ��� ⃗   −1 P  where ��� ⃗   −1 is the -th row in the inverse matrix �� ⃗  −1 .At the contacts, PΨ is given by the expression in the right hand side of Eq. ( 41) except that all the Bloch states s are multiplied by the inverse of the corresponding Bloch eigenvalue z −1 .Hence, one finds the expansion coefficients and obtains the closed equation For a sparse banded Hamiltonian matrix H D Eq. ( 45) can be solved recursively.However, similar to other recursive numerical schemes, the exact solution may become prohibitively difficult for the RSDFT Hamiltonian in realistic devices.On the other hand, solving this equation in the EM representation presents no difficulty even for fairly large systems.Let P EM be a projector to the subspace of the EM basis functions.If the accuracy of the EM basis is good enough, the coupling term P EM H D (1 − P EM ) is supposed to be a small correction and the Green's function G EM = ΦgΦ T in the EM representation nearly diagonalizes the left hand size of the equation.Thus one can apply an iterative algorithm from the family of complex conjugate gradient methods [57] and use the preconditioner in the form G EM + (1 − P EM ).
As a test, we consider a thin 1 nm-diameter nanowire in Fig. 1 and take g = g 0 as a Green's function in the close sys- tem.Applying the preconditioning is equivalent to solving the linear equation where x is the P EM projection of the residual in the original equation.The equation can be again solved by the "doublesided" R-matrix recursion technique.Assuming the x n , n be the components within the n-th EM unit cell one calculates n by the recursion where we introduced an auxiliary vector y which needs to be evaluated in the course of the R-matrix propagation Eq. (34).The corresponding recursive formulas are easily found from Eq. ( 46): At the propagation step n → n + 1 after computing the new r 22 in Eq. (34a) two additional operations need to be performed before completing the R-matrix propagation step Eq.(34b-d).The propagation starts at n = 2 with the ini- tial condition r ij = (E − h 2 ) −1 and y 2 = r 22 x 2 .After the last propagation step the solution in the boundary unit cells are found as and all other s are computed from Eq. ( 47) in the reverse order.The EM transformation completes the preconditioning step.The calculations for all the open channels should be performed in as a single matrix operation in order to improve the computer performance in much the same way as for the previously discussed electric density calculations.Similar approach can used for the preconditioner g R .The only dif- ference is that the complex-valued self-energies must be taken into account at the last propagation step.Most time consuming part of solving Eq. ( 45) is the preconditioning and the simulation time for one iteration is analogous to computing the contribution to the carrier density in the EM representation.Figure 15 presents an example of the I-V characteristics in a n-SiNW FET.The right panels show the potential and one-dimensional profile of the mobile density at low gate voltage.As a demonstration, we show the convergence of partial density contribution to this OFF state from incoming scattering modes at the source probe at low energy E = −1.1 eV. Figure 16 illustrates the behavior or the residual error in the bicongugate iterative scheme for solving Eq. (45).The solution of the scattering problem Eq. ( 45) is well localized in the source region and the contribution to the electronic density in the drain contact from these states can be neglected.However, as may be seen from the left panel, even in this deep tunneling regime, the transmission probability Eq. (43) shows steady convergence to the exponentially small tunneling value ( T ∼ 3 × 10 −3 ). Figure 17 shows the partial contribution   45) after 10 and 20 iterations.In this particular case N iter = 30 iterations give a few percent level of accuracy.The simulations have been performed using 16 CPUs of the Intel workstation and computing the contribution from one energy point takes ≤ 20 sec which is ∼ 10 times faster compared to the exact R-matrix propagation.
The present example may seem not to be informative enough as the EM representation in this case has been constructed to fit the quantum states of the periodic nanostructure.One should expect the simulations in inhomogeneous or amorphous nanostructures to be much more time consuming.However, since the time dependence on the total size N grid is linear, the approximate iterative schemes may offer a promising alternative to the exact recursive schemes.Convergence of the method and simulation time depends on the ab initio model for describing the crystal imperfections, the method for constructing the local basis representation and the choice of the iterative algorithm.These issues need to be addressed in future studies.

Summary
We have presented a first principle device simulator which is suitable for massively-parallel computers.The developed computer code makes it possible to calculate transport characteristics within a reasonable simulation time in large quasi-one-dimensional nanostructures with thousands of atoms in the channel cross-section.It is based on the previously reported real-space first-principle density functional program which efficiently performs large-scale parallel calculations in realistic systems.Necessary modifications are introduced to the geometry of the simulation domain as well as the boundary condition for the electronic structure calculations to be combined with the appropriate Poisson solver.An effective algorithm for constructing equivalent model (EM) representation has been developed in order to transform the original RSDFT data into a low rank quantum transport model suitable for practical parallel simulations.Coupled with the real-valued R-matrix propagation algorithm, the method greatly reduces the computational cost and enabled one to perform self-consistent device simulations within a reasonable time frame.

Fig. 1 (Fig. 2 (
Fig. 1 (Color online) The valence band (left) and conduction band (right) in a 1 nm-diameter Si wire.Solid lines represent the band structure in the original RSDFT simulator and squares correspond to the device Hamiltonian for the modified supercell geometry

Fig. 3 (Fig. 4 (
Fig. 3 (Color online) Constructing the EM basis in a 1 nm-diameter Si wire.The solid red lines represents the band structure in the reduced model before (left) and after (right) minimizing the number of states

Fig. 5 (
Fig. 5 (Color online) a Mixing RSDFT and the contact self-energy in the EM representations.b Transition probability and (c) the band structure in a 1 nm-diameter ideal Si wire

Fig. 6 (
Fig. 6 (Color online) I-V characteristics of a 8 nm-diameter n-SiNW FET whose schematic is shown in the inset of the left panel.The right panel shows the electric potential along the wire at four values of the gate voltage as indicated by the colored dots in the left panel

Figure 6 Fig. 7 (
Figure 6 presents an example of calculated I-V characteristics in a n-SiWN MOSFET shown in the inset.The right

Fig. 9 (
Fig. 9 (Color online) I-V characteristics of nanosheet n-Ge FETs with the conduction band in Fig. 8 with gate lengths of 10 nm and 14 nm.The right panel shows the corresponding electric potentials at selected values of the gate voltage

Fig. 10 (
Fig. 10 (Color online) a Oxidation of a Si wire in the RSDFT structure optimization simulations.At each step ≤ 50 oxygen atoms are introduced at random vacant positions outside the core region.b Extracting the core region for the reduced first principle transport simulations

Fig. 12 (
Fig.12 (Color online) Electric potential at various gate voltage in a SiNW FET using models (a) and (b) in Fig.11shown respectively by black and red lines.The results for a strained silicon in the right panel correspond to extending dielectric layer in the gate area Fig.13(Color online) The I-V characteristics in an n-SiNW MOS-FET using four types of RSDFT data in Fig.11

Fig. 14
Fig. 14 Electric charge in OFF (upper panels) and ON (low panels) states in a Si channel before (left) and (right) strain relaxation

2Fig. 15 (Fig. 16 (Fig. 17
Fig. 15 (Color online) The I-V characteristics in a 1 nm-diameter n-SiNW FET.The right panels shows the electric potential and mobile density in the OFF state used in the demonstration of iterative method Fig.17 (Color online) The partial contribution to the mobile charge in Fig. 15 from the low energy electrons coming from the left contact (red solid line) after 10 and 20 iterations.The black solid line shows the exact result−i[G R Σ < S G R † ] diag allows one to eliminate the spurious states from the targeted energy range [E 1 , E 2 ] and obtain the cor-