PaperlayerView exact revision on arXiv

Source-linked AI summary

Improvements on non-equilibrium and transport Green function techniques: the next-generation transiesta

Nick Papior, Nicolás Lorente, Thomas Frederiksen, Alberto García, Mads Brandbyge

arXiv:1607.04464v1cond-mat.mes-hall

TL;DR

The paper addresses the challenge of efficient and flexible atomistic transport calculations for equilibrium and multi-electrode non-equilibrium systems. It rewrites transiesta and develops complementary tbtrans/phtrans methods for generalized electrodes, contour integration, electrostatic effects, and scalable transport computation. The resulting codes support customized electrode conditions and calculations exceeding 10,000 orbitals for DFT-NEGF and 1,000,000 orbitals for tight-binding transport.

  • Problem

    Atomistic non-equilibrium transport requires efficient methods for open systems with multiple electrodes, differing chemical potentials and temperatures, and demanding matrix calculations.

  • Method

    The paper completely rewrites transiesta and implements generalized multi-electrode Green-function methods alongside optimized contours, electrostatic gates, matrix inversion, pivoting, and tbtrans/phtrans post-processing.

  • Results

    The reimplemented codes provide scalable N_e ≥1 electrode techniques and enable calculations exceeding 10,000 orbitals with DFT-NEGF and 1,000,000 orbitals with tight-binding methods.

  • Takeaways & Limitations

    The codes support equilibrium surface calculations, non-equilibrium thermoelectric studies, gated devices, electron and phonon transport, and flexible user-defined tight-binding interfaces.

Abstract

from arXiv · show

We present novel methods implemented within the non-equilibrium Green function code (NEGF) transiesta based on density functional theory (DFT). Our flexible, next-generation DFT-NEGF code handles devices with one or multiple electrodes ($N_e\ge1$) with individual chemical potentials and electronic temperatures. We describe its novel methods for electrostatic gating, contour opti- mizations, and assertion of charge conservation, as well as the newly implemented algorithms for optimized and scalable matrix inversion, performance-critical pivoting, and hybrid parallellization. Additionally, a generic NEGF post-processing code (tbtrans/phtrans) for electron and phonon transport is presented with several novelties such as Hamiltonian interpolations, $N_e\ge1$ electrode capability, bond-currents, generalized interface for user-defined tight-binding transport, transmission projection using eigenstates of a projected Hamiltonian, and fast inversion algorithms for large-scale simulations easily exceeding $10^6$ atoms on workstation computers. The new features of both codes are demonstrated and bench-marked for relevant test systems.

I. INTRODUCTION

The paper motivates a DFT-NEGF rewrite for atomistic transport, extending calculations to multi-electrode non-equilibrium systems while improving numerical integration and implementation efficiency.

  • Motivation: Atomistic transport remains difficult because realistic devices require non-equilibrium many-body treatment beyond static one-electron scattering.The paper notes that DFT is fundamentally a ground-state theory and faces limitations beyond linear conductance.
  • Contributions: The complete transiesta rewrite targets performance, multi-terminal equations, temperature gradients, electrostatic gating, contour convergence, and flexible tight-binding support.These improvements are presented as the main contributions of the new code.
  • Framework: DFT-NEGF seeks a self-consistent electron density and effective Kohn-Sham Hamiltonian for open systems coupled to electrodes in local equilibrium.When electrodes differ in distributions, the device develops a non-equilibrium electron density.
  • Framework: Bloch k-point sampling is applied to periodic directions in both self-consistent DFT-NEGF and subsequent transport calculations.The electrode k mesh must be correspondingly denser when a smaller cell is repeated in the device unit cell.
  • Framework: The density matrix uses electrode spectral functions weighted by Fermi distributions, while the energy density matrix enables non-equilibrium force calculations.The supplied equations define the density and energy density matrices within the Green-function framework.
  • Complex contour integration: Equilibrium integrations are evaluated through complex-plane contours because the Green function is smoother away from real-axis poles.Circle and square contours are mathematically equivalent, with retarded and advanced integrations taken in opposite half-planes.

B. Non-equilibrium (NEGF)

The multi-electrode non-equilibrium formulation constructs the density matrix from equilibrium contributions plus corrections associated with differing electrode distributions. This confines real-axis integrations to the bias window but increases matrix-operation demands.

  • Non-equilibrium formulation: Non-equilibrium transport arises when electrodes differ in chemical potential, electronic temperature, or both.The bias window is bounded by the minimum and maximum electrode chemical potentials, with appropriate Fermi-function tails.
  • Non-equilibrium formulation: The density matrix can be expressed as an equilibrium electrode contribution plus correction terms from electrodes with different Fermi distributions.These corrections confine the real-axis integral to the bias window.
  • Non-equilibrium formulation: Equivalent Fermi distributions eliminate corresponding correction terms and reduce the number of distinct density-matrix expressions.For N_e > 2 electrodes with two different Fermi distributions, only two mathematically equivalent expressions are needed.
  • Non-equilibrium formulation: The same electrode-distribution decomposition applies to the non-equilibrium energy density matrix.The paper states that the energy-density treatment follows the density-matrix construction similarly.
  • Computational cost: Compared with equilibrium, non-equilibrium density calculations require additional triple matrix products for electrode spectral functions.These operations are performed in addition to Green-function inversions.

III. IMPLEMENTATION DETAILS IN TRANSIESTA

Transiesta introduces optimized contour integration and a weighted treatment of non-equilibrium density matrices to improve convergence and reduce electrode-choice arbitrariness.

  • A. Complex contour optimization: The right-side quadrature scheme slightly reduces equilibrium contour points without loss of accuracy.It integrates over the connected retarded and advanced contour using the right half of the abscissa.
  • A. Complex contour optimization: In the gold-chain test, the circle contour benefits from the right-side scheme, whereas the square contour does not.The circle right-side scheme performs best in this setup.
  • A. Complex contour optimization: The Simpson quadrature method converges slowly, while the continued-fraction scheme converges fast and is simple.Changing the initial L+ points or pole count changes convergence properties.
  • B. Weighing ρ and bound states: Bound states make electrode-based density expressions unequal because they do not couple to electrodes through spectral functions.Their filling depends on the equilibrium electrode used, motivating equal weighting of equilibrium contributions when bound states are present.
  • B. Weighing ρ and bound states: The non-equilibrium density matrix is averaged over numerically unequal electrode-based expressions to avoid arbitrary equilibrium-electrode selection.The weighting is extended to multi-terminal systems and chosen to minimize the final density variance.

C. Inversion algorithms and performance

Transiesta implements direct, sparse, and block-tri-diagonal inversion strategies, with recursive propagation methods that reduce memory and computation for NEGF quantities.

  • C. Inversion algorithms and performance: Three inversion algorithms—LAPACK direct, MUMPS sparse, and block-tri-diagonal—are implemented in Γ-only and arbitrary-k variants.The variants limit memory usage where applicable.
  • C. Inversion algorithms and performance: Equation (16) is implemented for all three inversion algorithms, substantially reducing memory requirements, particularly for BTD.The approach targets the Green-function block products needed for non-equilibrium spectral functions.
  • C. Inversion algorithms and performance: The block-tri-diagonal algorithm is transiesta’s default and represents the recursive Green function method for quasi-1D block-tridiagonal matrices.It propagates self-energies through the matrix to obtain Green-function blocks.
  • C. Inversion algorithms and performance: The BTD matrix is expressed through diagonal and off-diagonal blocks Ai, Bi, and Ci representing the non-zero Hamiltonian and overlap structure.Downfolded self-energies connect neighboring sequences and may be split among submatrices.
  • C. Inversion algorithms and performance: The BTD formulation computes only selected diagonal and neighboring Green-function blocks required by the density matrix.Spectral-function propagation reuses intermediate matrices and avoids recalculating them for each electrode.

2. Orbital pivoting for minimizing bandwidth

Transiesta reduces block-tridiagonal bandwidth through graph-based and ordering-based pivoting, improving inversion performance while accommodating multi-electrode structures.

  • 2. Orbital pivoting for minimizing bandwidth: BTD performance is determined solely by Hamiltonian bandwidth, represented by the size of the Ai blocks.Bandwidth reflects the quasi-1D size of the system.
  • 2. Orbital pivoting for minimizing bandwidth: Five pivoting methods minimize matrix bandwidth, including Hamiltonian connectivity, peripheral connectivity, and Cuthill-McKee ordering.The implemented methods address sparsity patterns that rarely have minimum bandwidth under the atomic input sequence.
  • 2. Orbital pivoting for minimizing bandwidth: Atomic-graph pivoting can greatly improve performance because atomic and orbital connectivity graphs differ little in these systems.The method uses the atomic graph rather than the orbital graph.
  • 2. Orbital pivoting for minimizing bandwidth: Sorting each Ai block can split the broadening matrix Γe into two consecutive memory blocks and reduce tri-diagonal bandwidth.Splitting Γe is more beneficial than retaining a single block when possible.

3. Performance of inversion algorithms

Benchmarking compares MUMPS, NEGF, and BTD inversion methods on graphene systems with different lengths and bias conditions, showing strong BTD scaling.

  • 3. Performance of inversion algorithms: The benchmark uses pristine graphene with square electrodes containing 48 atoms and lengths up to approximately 6,000 orbitals.A high-bias case uses 0.75 V and equal non-equilibrium and equilibrium contour contributions.
  • 3. Performance of inversion algorithms: For quasi-1D systems with two electrodes, naive partitioning is already the best partitioning shown, while other electrode counts benefit from improved pivoting.Self-energy ordering can be swapped freely in the illustrated systems.
  • 3. Performance of inversion algorithms: BTD reaches around 100 times better performance than the reference methods at 5,000 orbitals.Larger systems are expected to increase the speedup further.
  • 3. Performance of inversion algorithms: MUMPS performs very well for EGF calculations and improves for systems larger than approximately 5,000 orbitals as sparsity increases.NEGF is relatively inefficient in the tested setup because of column clustering.

4. Parallelization

transiesta combines energy-point parallelism with hybrid OpenMP/MPI threading to reduce inversion cost and memory pressure. The BTD approach delivers major speedups, while threading efficiency depends on block size.

  • Hybrid parallelization: Hybrid parallelization reduces memory pressure by combining MPI energy-point distribution with OpenMP threading across matrix operations.Pure energy-point parallelization requires each MPI process to hold complete matrices, which can reach gigabyte scale.
  • Threading performance: 95% parallel fraction was obtained across the investigated thread counts for the 9,130-orbital graphene test system.The system used 11 BTD blocks with average block size NB = 830 and showed extremely good scaling up to 4 threads.
  • Benchmark: More than 40× speedup over direct LAPACK was achieved by BTD for the largest pristine-graphene system.MUMPS began gaining speed after 5,000 orbitals.
  • Self-energy calculations: Smaller electrode unit cells reduce self-energy matrix size, making cubic scaling more favorable despite requiring calculations at n k-points.The method uses Bloch expansion when transverse periodicity is present.

E. Electrostatics in NEGF

The electrostatics implementation addresses open-boundary conditions, gating, and charge neutrality in multi-electrode NEGF systems. Custom Hartree guesses, multigrid solutions, and potential shifts improve the treatment of difficult electrostatic configurations.

  • Boundary conditions: Custom Hartree guesses are recommended for Ne > 2 because they can improve charge conservation and SCF convergence.The initial Hartree guess is linear in electrode chemical-potential differences, so only one guess calculation is needed.
  • Boundary conditions: Multigrid Poisson solutions provide initial Hartree potentials for complex six-electrode geometries with chemical potentials [−V/2, 0, V/2].The Poisson solution varies linearly with V, and the setup includes buffer atoms behind each electrode.
  • Electrostatic gating: Electrostatic gates are implemented as additional non-interacting electrodes, extending gating methods to DFT and transiesta calculations.The implementation addresses electrostatic terms associated with charge in the gate material.
  • Electrostatic gating: λ = 2.4 Å was obtained for gate-induced electronic screening in a strongly gated 10-layer graphene stack.The decay length depends on electric field and doping level and agrees with reported values for similar gate levels.
  • Charge conservation: A device-region potential shift enforces qD = 0 to mitigate excess charge and difficult SCF convergence.The shift is estimated from the device DOS at a reference energy and updated during the SCF loop.

G. Thermoelectric effects under NEGF

The NEGF implementation allows each electrode to have an independent chemical potential and electronic temperature, enabling self-consistent thermoelectric calculations. The post-processing tools extend transport analysis to multi-electrode electronic and phononic systems.

  • Thermoelectric effects: Independent electrode temperatures enable self-consistent NEGF calculations with temperature differences between electrodes.The example fixes TL = 3000 K and TR = 300 K for a weakly coupled 1D carbon junction.
  • Thermoelectric effects: Uniform-temperature self-consistency cannot determine the fully self-consistent IV curve when electrode temperatures differ.In Fig. 10b, the fully self-consistent result is the black curve, whereas uniform-temperature calculations produce red or blue curves.
  • Post-processing: tbtrans and phtrans calculate electronic, thermal-energy, and phonon transport for general multi-electrode setups.They also provide density-of-states, projected-transmission, transmission-eigenvalue, and bond-current analyses.
  • Post-processing: Transmission analysis includes scattering-matrix or generalized Fisher–Lee formulations, eigenvalues, currents, and thermal energy transfer.Transmission eigenvalues support quantities such as Fano factors, while current and heat transfer are computed from transmission functions.
  • Transport relations: Time-reversal symmetry yields relations between opposite-direction transmissions and connects net work with net heat supplied.The stated relation is Qee′ + Qe′e = (µ′e − µe)/|e| Iee′ ≡ W.

B. Inversion algorithm — again

The BTD inversion strategy down-folds electrode self-energies into a smaller region and uses separate pivoted matrices for multi-terminal transport. Hamiltonian interpolation and large-scale tight-binding construction further reduce repeated transport costs.

  • BTD inversion: Down-folding lets tbtrans calculate transmission and current in a smaller subspace rather than the full device space.For multi-terminal systems, the selected region cannot directly couple electrode branches or their self-energy origins become entangled.
  • BTD inversion: The smallest suitable block maximizes throughput when only transmission or current is required.Larger retained regions are needed for analyses such as Mulliken charges.
  • BTD inversion: A six-terminal calculation uses Ne + 1 BTD matrices, with separate electrode matrices and one matrix for the central region.Each matrix has its own pivoting scheme to reduce bandwidth and improve performance.
  • Hamiltonian interpolation: Spline Hamiltonian interpolation agrees better than linear interpolation with the self-consistent transmission for the Cu-tip–C60/Cu(111) junction.The comparison is made at V = −1.5 V using converged Hamiltonians from −2 V through 2 V.

D. tbtrans as transport back-end and feature generalization

tbtrans is presented as a flexible, stand-alone NEGF transport back-end supporting user-defined Hamiltonian corrections, subspace transport analysis, and molecular-state transmission projections.

  • Feature generalization – δH: tbtrans accepts Hamiltonians through NetCDF-4 files, enabling generic tight-binding transport without Fortran or SIESTA binary-format expertise.The sisl Python package supports creation and manipulation of these Hamiltonians.
  • Feature generalization – δH: Users can define transport features by supplying δHk(ϵ) Hamiltonian corrections, including four combinations of energy and k-point dependence.The variants are δH, δHk, δH(ϵ), and δHk(ϵ).
  • Molecular state projection transmission: tbtrans decomposes transport into selected device subspaces and stores MPSH-state coupling quantities for flexible transmission analysis.The method uses orthogonalized molecular projected self-consistent Hamiltonian states.
  • Molecular state projection transmission: For a two-molecule bridge, projected transmission can resolve electrons injected through one molecular state and extracted through specified states of the other molecule.The schematic contains two states for molecule A and three for molecule B.
  • Molecular state projection transmission: Incoming, outgoing, and simultaneous projectors generally differ because of asymmetric electrode coupling or electrode hybridization.The projection states may also be k-resolved or Γ-point-only.
  • Molecular state projection transmission: In C60/Cu(111), one LUMO orbital carries most transmission, while the other two nearly degenerate LUMOs contribute negligibly.The comparison uses 13×13 Monkhorst-Pack sampling and includes Γ- and k-point projector choices.

F. Phonon transport — phtrans

phtrans extends the Green-function transport framework to phonons by replacing electronic matrices and energies with dynamical matrices and squared frequencies.

  • Phonon transport — phtrans: Phonon transport in phtrans uses the electronic transmission formalism with the Hamiltonian replaced by the dynamical matrix, unity overlap, and ε + iη replaced by ω^2 + iη^2.The implementation retains tbtrans functionality, including Ne ≥1 terminals.
  • Phonon transport — phtrans: sisl extracts dynamical matrices from gulp and writes phtrans-compatible files, enabling empirical-potential phonon transport for very large systems.This provides an interface to third-party tools.
  • Phonon transport — phtrans: At 600 K, the GB558 grain boundary reduces graphene’s heat transport to approximately 60% of pristine graphene’s value.The comparison uses the Brenner potential and q-averaged transmissions.
  • Phonon transport — phtrans: GB558 hosts localized out-of-plane modes near ℏω = 115 meV and in-plane modes near ℏω = 200 meV.These modes appear as peaks in the atom-resolved projected phonon DOS.

V. CONCLUSIONS

The reimplemented transiesta, tbtrans, and phtrans codes generalize DFT-NEGF transport across electrode configurations, gating, transport projections, phonons, and scalable computation.

  • CONCLUSIONS: The framework covers equilibrium single-electrode and non-equilibrium multi-electrode calculations with independently specified electrode chemical potentials and temperatures.This supports equilibrium surface calculations and non-equilibrium thermoelectric studies.
  • CONCLUSIONS: Contour-integration schemes improve density-matrix convergence while reducing equilibrium integration abscissas, and generalized weighting supports Ne > 2.Flexible charge/Hartree gate geometries enable gate-effect studies in the device region.
  • CONCLUSIONS: A BTD-based inversion method with bandwidth-reducing pivoting and OpenMP threading produced an approximately 100× speed-up for one test system.The implementation also includes an efficient spectral-function propagation algorithm.
  • CONCLUSIONS: tbtrans supports subspace transport quantities, bond-currents, and molecular-orbital projected transmissions, while phtrans applies the framework to phonons.Projection supports both Γ- and k-point projectors.
  • CONCLUSIONS: Finite-bias Hamiltonian interpolation reduces the computational burden of full I−V curves, with spline interpolation reported as accurate for complex systems.sisl provides a generic Python interface and can pass dynamical matrices from gulp to phtrans.
  • CONCLUSIONS: The codes support Ne ≥1 electrodes and everyday calculations exceeding 10,000 orbitals in DFT-NEGF and 1,000,000 orbitals in tight-binding transport.The conclusion attributes this capability to scalable and efficient algorithms.
Loading 1607.04464v1…