API Reference

This section provides detailed API documentation for all gauNEGF modules.

Core Modules

NEGF Base Class

Self-consistent field (SCF) implementation for Non-Equilibrium Green’s Function calculations.

This module provides the base NEGF class for performing self-consistent DFT+NEGF calculations using Gaussian quantum chemistry package. It implements the energy-independent self-energy approach (Damle et al. [1]) for molecular transport calculations. The module handles:

• Integration with Gaussian for DFT calculations

• SCF convergence with Pulay DIIS mixing [2]

• Contact self-energy calculations

• Voltage bias and electric field effects

• Spin-polarized calculations

References

[1] Damle, P., Ghosh, A. W., & Datta, S. (2002). First-principles analysis of molecular

conduction using quantum chemistry software. Chemical Physics, 281(2-3), 171-187. DOI: 10.1016/S0301-0104(02)00496-2

[2] Pulay, P. (1980). Convergence acceleration of iterative sequences. The case of SCF

iteration. Chemical Physics Letters, 73(2), 393-398. DOI: 10.1016/0009-2614(80)80396-4

class gauNEGF.scf.NEGF(fn, basis='chkbasis', func='hf', spin='r', fullSCF=True, route=None, section=None, nPulay=4)

Bases: object

Non-Equilibrium Green’s Function calculator integrated with Gaussian DFT.

Manages the self-consistent field cycle between NEGF transport calculations and Gaussian DFT. Uses constant (energy-independent) self-energies; for energy-dependent calculations see the NEGFE subclass.

Parameters:

• fn (str) – Base filename for Gaussian input/output files (without extension)

• basis (str, optional) – Gaussian basis set name (default: ‘chkbasis’)

• func (str, optional) – DFT functional to use (default: ‘hf’)

• spin ({'r', 'u', 'ro', 'g'}, optional) – Spin configuration: - ‘r’: restricted - ‘u’: unrestricted - ‘ro’: restricted open-shell - ‘g’: generalized open-shell (non-collinear) (default: ‘r’)

• fullSCF (bool, optional) – Whether to run full SCF or use Harris approximation (default: True)

• route (str, optional) – Additional Gaussian route commands (default: ‘’)

• section (str, optional) – Gaussian input section specification (default: None)

• nPulay (int, optional) – Number of previous iterations to use in Pulay mixing (default: PULAY_MIXING_SIZE from gauNEGF.config)

F

Fock matrix in Hartree (multiply by 27.211386 for eV)

Type:

ndarray

P

Density matrix

Type:

ndarray

S

Overlap matrix

Type:

ndarray

fermi

Fermi energy in eV

Type:

float

nelec

Number of electrons

Type:

float

FockToP()

Calculate density matrix from Fock matrix using energy-independent approach.

Diagonalizes the orthogonalized Fock+Sigma matrix, integrates analytically for constant self-energies (Damle et al. 2002), and transforms back to the non-orthogonal basis. Updates self.P and optionally self.fermi.

Returns:

(eigenvalues, occupations) sorted by energy

Return type:

tuple

PMix(damping, Pulay=False)

Mix old and new density matrices using damping or Pulay DIIS method.

Applies simple linear damping by default. When Pulay=True, applies the DIIS update using the accumulated history buffer (length set in __init__).

Parameters:

• damping (float) – Mixing parameter between 0 and 1

• Pulay (bool, optional) – Whether to use Pulay mixing (default: False)

Returns:

(RMSDP, MaxDP) - RMS and maximum density matrix differences

Return type:

tuple

PToFock()

Calculate new Fock matrix from current density matrix using Gaussian.

Returns:

Energy difference from previous iteration

Return type:

float

SCF(conv=0.001, damping=0.02, maxcycles=100, checkpoint=True, pulay=True)

Run self-consistent field calculation until convergence.

Alternates FockToP / PMix / PToFock until all three criteria (dE, RMSDP, MaxDP) fall below conv, or maxcycles is reached. Pulay DIIS is applied every nPulay iterations when pulay=True.

Parameters:

• conv (float, optional) – Convergence criterion for energy and density (default: SCF_CONVERGENCE_TOL)

• damping (float, optional) – Mixing parameter between 0 and 1 (default: SCF_DAMPING)

• maxcycles (int, optional) – Maximum number of SCF cycles (default: SCF_MAX_CYCLES)

• checkpoint (bool, optional) – Save density matrix at each iteration and load if job interrupted (default: True)

• pulay (bool, optional) – Whether to use Pulay DIIS mixing (default: True)

Returns:

(count, PP, TotalE) - cycle number, number of electrons, and DFT energy at each cycle

Return type:

tuple

__init__(fn, basis='chkbasis', func='hf', spin='r', fullSCF=True, route=None, section=None, nPulay=4)

Initialize NEGF calculator and run initial DFT calculation.

Sets up the calculator with specified parameters and runs an initial DFT calculation using Gaussian to obtain the starting Fock and overlap matrices.

Parameters are documented in class docstring.

getHOMOLUMO()

Calculate HOMO and LUMO energies.

Diagonalizes the Fock matrix in orthogonalized basis to get orbital energies, then identifies HOMO and LUMO based on electron occupation and spin configuration.

Returns:

Array of [HOMO, LUMO] energies in eV

Return type:

ndarray

getSigma(E=0)

runDFT(fullSCF=True)

Run DFT calculation using Gaussian.

Performs either a full SCF calculation or generates an initial Harris guess using Gaussian. Updates the Fock matrix and orbital indices after completion.

Parameters:

fullSCF (bool, optional) – If True, runs full SCF to convergence. If False, uses Harris guess only. (default: True)

Notes

• Attempts to load from checkpoint file first

• Falls back to full calculation if checkpoint fails

• Updates self.F and self.locs with new Fock matrix

saveMAT(matfile='out.mat')

Save calculation results to MATLAB format file.

Parameters:

matfile (str, optional) – Output filename (default: “out.mat”)

Returns:

Fock matrix in orthogonalized basis

Return type:

ndarray

setContacts(lContact=None, rContact=None)

Set contact atoms and get corresponding orbital indices.

Identifies the orbital indices corresponding to the specified contact atoms. If no contacts specified, uses all orbitals.

Parameters:

• lContact (int or array-like, optional) – Atom number(s) for left contact. None means all atoms. (default: None)

• rContact (int or array-like, optional) – Atom number(s) for right contact. None means all atoms. (default: None)

Returns:

(left_orbital_indices, right_orbital_indices)

Return type:

tuple of ndarrays

setDen(P_, enableSpinLock=False, spinLockList=None)

Set the density matrix and update dependent quantities.

Updates the density matrix, stores it in Gaussian format, recalculates the number of electrons, and updates the Fock matrix. Optionally extracts and stores spin information for spin-locking.

Parameters:

• P (ndarray) – New density matrix

• enableSpinLock (bool, optional) – If True, extract and lock spin orientations from the density matrix. This is an advanced feature for specialized use cases. Default: False

• spinLockList (list, optional) – List of atoms to apply spin-locking to

setFock(F_)

Set the Fock matrix, converting from eV to Hartree units.

Parameters:

F (ndarray) – Fock matrix in eV units

setSigma(lContact=None, rContact=None, sig=-0 - 0.1j, sig2=None)

Set self-energies for left and right contacts.

Configures the contact self-energies, handling various input formats and spin configurations. Self-energies can be scalar, vector, or matrix, with automatic handling of spin degrees of freedom.

Parameters:

• lContact (array-like) – Atom numbers for left contact, all atoms if None (default: None)

• rContact (array-like) – Atom numbers for right contact, all atoms if None (default: None)

• sig (scalar or array-like, optional) – Self-energy for left contact. Can be: - scalar: same value for all orbitals - vector: one value per orbital - matrix: full self-energy matrix (default: -0.1j)

• sig2 (scalar, array-like, or None, optional) – Self-energy for right contact. If None, uses sig. Same format options as sig. (default: None)

Returns:

Left and right contact orbital indices

Return type:

tuple

Notes

• Handles spin configurations (‘r’, ‘u’, ‘ro’, ‘g’)

• Automatically expands scalar/vector inputs to full matrices

• Verifies matrix dimensions match Fock matrix

• Updates self.sigma1, self.sigma2, self.sigma12

Raises:

Exception – If matrix dimensions don’t match or invalid input format

setVoltage(qV, fermi=None, Emin=None, Eminf=None)

Set voltage bias and Fermi energy, updating electric field.

Applies a voltage bias between contacts and updates the chemical potentials and electric field. Can optionally set the Fermi energy and integration limits.

Parameters:

• qV (float) – Voltage bias in eV

• fermi (float, optional) – Fermi energy in eV. If not provided, will be calculated or use existing value (default: None)

• Emin (float, optional) – Minimum energy for integration in eV (default: None)

• Eminf (float, optional) – Lower bound energy in eV (default: None)

Notes

• Requires contacts to be set first

• Updates chemical potentials as fermi ± qV/2

• Calculates and applies electric field between contacts

• If fermi not provided, uses (HOMO+LUMO)/2 for initial guess

updateN()

Update the total number of electrons from the density matrix.

Calculates the number of electrons by tracing the product of the density and overlap matrices. For restricted calculations, multiplies by 2 to account for spin degeneracy.

Returns:

Total number of electrons

Return type:

float

updatePulay(nPulay)

Reset Pulay mixing buffers to use a new history length.

Allows changing the number of previous iterations used in Pulay DIIS mixing without reinitializing the full NEGF object. Buffers are re-seeded from the current density matrix, matching the initialization in __init__. Any accumulated mixing history is discarded.

Parameters:

nPulay (int) – New number of previous iterations to use in Pulay mixing.

writeChk()

Write current state to a Gaussian binary checkpoint (.chk) file.

Uses gauopen’s bar.writefile(), which routes through the BAF/unfchk path internally – no .fchk intermediate is needed.

Energy-Dependent NEGF

See also

1D Chain Contacts in gauNEGF – 1D chain contact setup via NEGFE.setContact1D

Bethe Lattice Contacts – Bethe lattice contact setup via NEGFE.setContactBethe

Choosing a Contact Type – Choosing between contact types

Energy-dependent extensions to the NEGF class for quantum transport calculations.

This module extends the base NEGF class to handle energy-dependent self-energies, temperature effects, and advanced Fermi energy search methods. It provides support for Bethe lattice and 1D chain contacts [1] with proper energy integration.

References

[1]

Jacob, D. & Palacios, J. J. Chem. Phys. 134, 044118 (2011)

class gauNEGF.scfE.NEGFE(fn, basis='chkbasis', func='hf', spin='r', fullSCF=True, route=None, section=None, nPulay=4)

Bases: NEGF

Extended NEGF class with energy-dependent self-energies and temperature effects.

This class extends the base NEGF implementation to handle: - Energy-dependent contact self-energies - Finite temperature effects - Advanced Fermi energy search methods - Integration methods for the density matrix

Inherits all attributes and methods from NEGF class.

FockToP()

Calculate density matrix using energy-dependent integration.

Performs complex contour integration for the equilibrium part and real axis integration for the non-equilibrium part. Updates Fermi energy using specified method if required.

Returns:

(energies, occupations) - Sorted jnp.linalg.eigenvalues and occupations

Return type:

tuple

PToFock()

Calculate new Fock matrix and update surfG object.

Returns:

Energy difference from previous iteration

Return type:

float

getSigma(E, E2=None)

Get contact self-energies at specified energy.

Parameters:

E (float) – Energy in eV

Returns:

(left_sigma, right_sigma) - Contact self-energies

Return type:

tuple

integralCheck(cycles=10, damp=0.02, pauseFermi=False)

Check and optimize integration parameters.

Parameters:

• cycles (int, optional) – Number of SCF cycles to run (default: 10)

• damp (float, optional) – Damping parameter (default: 0.02)

• pauseFermi (bool, optional) – Whether to pause Fermi updates (default: False)

setContact1D(contactList, tauList=None, stauList=None, alphas=None, aOverlaps=None, betas=None, bOverlaps=None, neList=None, muList=None, eta=1e-05, T=0.0, symmetrize_contacts=None)

Set energy-dependent 1D chain contacts.

Parameters:

• contactList (list) – List of atom indices for contacts

• tauList (list or None, optional) – Coupling matrices or atom indices (default: None)

• stauList (list or None, optional) – Overlap matrices (default: None)

• alphas (array_like or None, optional) – On-site energies (default: None)

• aOverlaps (array_like or None, optional) – On-site overlaps (default: None)

• betas (array_like or None, optional) – Hopping energies (default: None)

• bOverlaps (array_like or None, optional) – Hopping overlaps (default: None)

• neList (list or None, optional) – Number of electrons per unit cell (default: None)

• muList (list or None, optional) – fermiEnergy for each contact in eV (default: None)

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

• T (float, optional) – Temperature in Kelvin (default: 0)

• symmetrize_contacts (bool or None, optional) – Enforce identical on-site parameters for both contacts during SCF. True for same-material contacts, False for junctions. None (default) auto-detects based on setup pattern.

Returns:

Left and right contact orbital indices

Return type:

tuple

setContactBethe(contactList, latFile='Au', eta=1e-05, T=0.0)

Set energy-dependent Bethe lattice contacts.

Parameters:

• contactList (list) – List of atom indices for contacts

• latFile (str, optional) – Lattice parameter file (default: ‘Au’)

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

• T (float, optional) – Temperature in Kelvin (default: 0)

Returns:

Left and right contact orbital indices

Return type:

tuple

setIntegralLimits(N1=None, N2=None, Nnegf=None, tol=0.0001, Emin=None)

Set integration parameters for density calculation.

Two modes:

• Default (Emin is None, tol given): place Emin below the band via calcEmin (true-Sigma DOS loop) and pin the deep bound Eminf to the wide config value ENERGY_MIN (TSW unused, set None). calcTSW / calcPseudoPoleFloor are NOT on this path – the lower contour is handled analytically by damleLowerDensity, whose cross-term warning in FockToP replaces calcTSW’s dTSW<0 detection. Used at initial setup (called by setContact1D) so the first FockToP cycle has sensible bounds; FockToP refreshes these each SCF cycle thereafter.

• Fixed-grid (N1/N2 given): DEPRECATED for production. The adaptive integration in FockToP replaces fixed-grid Emin/Eminf for the default tol path; N1/N2 is retained only for code paths that explicitly opt into fixed-point quadrature.

Parameters:

• N1 (int, optional) – Number of points for complex contour (default: None)

• N2 (int, optional) – Number of points for real axis (default: None)

• Nnegf (int, optional) – Number of points for non-equilibrium integration (default: None)

• tol (float, optional) – Tolerance for integration (default: 1e-4)

• Emin (float or None, optional) – Minimum energy for integration (default: None)

setSigma(lContact=None, rContact=None, sig=-0 - 0.1j, sig2=None, T=0.0)

Set constant self-energy contacts with temperature.

Parameters:

• lContact (list) – Atom numbers for left contact, all atoms if None (default: None)

• rContact (array-like) – Atom numbers for right contact, all atoms if None (default: None)

• sig (complex, optional) – Left contact self-energy (default: -0.1j)

• sig2 (complex or None, optional) – Right contact self-energy (default: None)

• T (float, optional) – Temperature in Kelvin (default: 0)

Returns:

Left and right contact orbital indices

Return type:

tuple

setVoltage(qV, fermi=None, Emin=None, Eminf=None, fermiMethod=None)

Set voltage bias and Fermi search method.

Parameters:

• qV (float) – Applied voltage in eV

• fermi (float, optional) – Fermi energy in eV (default: None)

• Emin (float, optional) – Minimum energy for integration (default: None)

• Eminf (float, optional) – Minimum energy for Fermi search (default: None)

• fermiMethod (str, optional) – Method for Fermi search: ‘muller’, ‘secant’, ‘predict’ or ‘poly’ (default: ‘muller’)

spawnNEGF(mu1=None, mu2=None)

Density Module

Density matrix calculation methods for quantum transport simulations.

This module provides functions for calculating density matrices in quantum transport calculations using various integration methods:

• Analytical integration for energy-independent self-energies

• Complex contour integration for equilibrium calculations

• Real-axis integration for non-equilibrium calculations

• Parallel processing support for large systems

Notes

The module supports both serial and parallel computation modes, automatically selecting the most efficient method based on system size and available resources. Temperature effects are included through Fermi-Dirac statistics.

gauNEGF.density.bisectFermi(V, Vc, D, Gam, Nexp, conv=0.001, Eminf=-1000000.0)

Find Fermi energy using bisection method.

Uses bisection to find the Fermi energy that gives the expected number of electrons. The search is performed using the analytical density matrix calculation.

Parameters:

• V (ndarray) – Eigenvectors of Fock matrix in orthogonalized basis

• Vc (ndarray) – Inverse conjugate transpose of V

• D (ndarray) – Eigenvalues of Fock matrix

• Gam (ndarray) – Broadening matrix Γ = i[Σ - Σ†] in orthogonalized basis

• Nexp (float) – Expected number of electrons

• conv (float, optional) – Convergence criterion for electron number (default: 1e-3)

• Eminf (float, optional) – Lower bound for integration in eV (default: -1e6)

Returns:

Fermi energy in eV that gives the expected electron count

Return type:

float

Notes

The search is bounded by the minimum and maximum eigenvalues. Warns if maximum iterations reached without convergence.

gauNEGF.density.calcEmin(F, S, g, tol=0.001, maxN=100, Emin=None, X=None)

Calculate minimum energy bound for numerical integration.

Finds the lower energy bound for density matrix integration by iteratively lowering Emin until the density of states at that energy falls below the specified tolerance. This ensures the integration bounds capture all occupied states without including energies too far below the band edge.

Parameters:

• F (ndarray) – Fock matrix in eV.

• S (ndarray) – Overlap matrix.

• g (surfG object) – Surface Green’s function calculator providing sigmaTot() and crossTermQTot() methods for energy-dependent self-energies.

• tol (float, optional) – Convergence tolerance for density of states (default: FERMI_CALCULATION_TOL).

• maxN (int, optional) – Maximum number of expansion iterations (default: MAX_CYCLES).

• Emin (float or None, optional) – Initial lower bound in eV. If None, initialized from minimum eigenvalue of X @ F @ X (symmetric Lowdin orthogonalization with X = S^(-1/2)) minus EMIN_BUFFER (default: None).

• X (ndarray or None, optional) – Lowdin orthogonalizer S^(-1/2). If None (default), computed locally via fractional_matrix_power(S, -0.5). Passing the caller’s cached X avoids recomputing it. Used only when Emin is None.

Returns:

Lower energy bound in eV where DOS is below tolerance.

Return type:

float

Notes

Prints a warning if maximum iterations are reached without convergence. X @ F @ X (X = S^(-1/2)) is used for the initial Emin seed because it is Hermitian by construction and its eigenvalues are the true device MO energies; eigh(inv(S) @ F) would silently fail for non-trivial S.

gauNEGF.density.calcFermi(g, ne, Emin, Ef, lBound=None, uBound=None, tol=0.0001, conv=0.001, maxcycles=10, T=0.0)

Calculate Fermi energy using bisection method with adaptive integration.

Parameters:

• g (surfG object) – Surface Green’s function calculator

• ne (float) – Target number of electrons

• Emin (float) – Lower bound for complex contour integration in eV

• Ef (float) – Initial guess for Fermi energy in eV

• lBound (float, optional) – Lower bound for bisection search in eV. If None, determined dynamically.

• uBound (float, optional) – Upper bound for bisection search in eV. If None, determined dynamically.

• tol (float, optional) – Tolerance for adaptive integration (default: ADAPTIVE_INTEGRATION_TOL)

• conv (float, optional) – Convergence criterion for electron count (default: FERMI_CALCULATION_TOL)

• maxcycles (int, optional) – Maximum number of iterations (default: FERMI_SEARCH_CYCLES)

• T (float, optional) – Temperature in Kelvin (default: TEMPERATURE)

Returns:

Calculated Fermi energy in eV

Return type:

float

gauNEGF.density.calcFermiBisect(g, ne, Emin, Ef, N, tol=0.0001, conv=0.001, maxcycles=10, T=0.0, uBound=None, lBound=None)

Calculate Fermi energy of system using bisection

gauNEGF.density.calcFermiMuller(g, ne, Emin, Ef, N, tol=0.0001, conv=0.001, maxcycles=10, T=0.0)

Calculate Fermi energy using Muller’s method, starting with 3 initial points

gauNEGF.density.calcFermiPolyFit(g, ne, Emin, Ef, N, tol=0.0001, conv=0.001, maxcycles=10, T=0.0, order=3)

Calculate Fermi energy using polynomial regression fit with accumulating points.

Uses all accumulated points to fit a polynomial of specified order through (n, E) pairs using least squares regression (polyfit) and finds where n = ne (target electron count).

Parameters:

• g (surfG object) – Surface Green’s function calculator

• ne (float) – Target number of electrons

• Emin (float) – Lower bound for complex contour in eV

• Ef (float) – Initial guess for Fermi energy in eV

• N (int or None) – Number of integration points (None for adaptive)

• tol (float, optional) – Tolerance for integration (default: 1e-5)

• conv (float, optional) – Convergence criterion for electron count (default: 1e-3)

• maxcycles (int, optional) – Maximum number of iterations (default: 20)

• T (float, optional) – Temperature in Kelvin (default: 300)

• order (int, optional) – Order of polynomial fit (default: 3)

Returns:

(Ef, dE, P, error) - Calculated Fermi energy, last energy step, density matrix, and final error

Return type:

tuple

Notes

This method builds up a history of (n, E) points and uses polynomial regression to fit an inverse polynomial (E as a function of n) to find E at n = ne. Uses numpy.polyfit for regression and numpy.roots to find the solution. Polynomial order is min(counter, order) to avoid overfitting with few points.

gauNEGF.density.calcFermiSecant(g, ne, Emin, Ef, N, tol=0.0001, conv=0.001, maxcycles=10, T=0.0)

Calculate Fermi energy using Secant method, updating dE at each step

gauNEGF.density.calcPseudoPoleFloor(F, S, g, alpha=0.1, min_buffer=50.0, E1=-1000.0, E2=-10000.0)

Compute a safe Eminf floor by detecting pseudo-poles of the device GF.

Pseudo-poles are spurious real-axis poles of G_DD^R(E) = (E*S - F - Sigma)^{-1} that arise from the asymptotically E-linear piece of Sigma when the contact is non-orthogonal. If the contour integration captures them, P_DD acquires negative natural-orbital occupations and is no longer DFT-compatible.

All matrices are treated as complex Hermitian throughout, so this works for GHF / SOC where F, S, X, Sigma_0 have nontrivial imaginary parts as well as the real-symmetric closed-shell case.

See docs/pseudo_pole_handling.md for the derivation and worked examples.

Algorithm:

1. Two-point asymptotic probe of Sigma:

   X = (Sigma(E2) - Sigma(E1)) / (E2 - E1) Sigma_0 = Sigma(E1) - E1 * X

2. Form the effective generalized eigenproblem H_eff v = E S_eff v that governs deep-|E| poles of G_DD^R:

   H_eff = F + Sigma_0 (Hermitian) S_eff = S - X (Hermitian, possibly indefinite)

3. Reduce to the standard non-Hermitian eigproblem T = S_eff^{-1} H_eff and solve with the JAX general eig (via gauNEGF.utils.eig). Eigenvalues are real for Hermitian (H, S) even with indefinite S; take Re(…) to drop floating-point imaginary noise.

4. Classify each eigenvector by its Hermitian S_eff-norm:

   v^H S_eff v > 0 -> physical pole (not returned) v^H S_eff v < 0 -> pseudo-pole (returned)

   Number of pseudo-poles equals the number of negative eigenvalues of S_eff (Sylvester’s law of inertia).

5. Filter to negative-energy pseudo-poles only (positive ones cannot be excluded by Eminf, which is a LOWER bound). For each, compute floor_i = E_pp_i + max(alpha * |E_pp_i|, min_buffer); return the shallowest (largest) floor_i, or ENERGY_MIN if no negative pseudo-poles exist or if the shallowest buffered floor is >= 0.

Parameters:

• F (ndarray (N, N)) – Device Fock matrix in eV. Complex Hermitian for GHF/SOC, real-symmetric otherwise; both are handled.

• S (ndarray (N, N)) – Device overlap matrix. Complex Hermitian, positive-definite.

• g (surfG-like object) – Surface Green’s function calculator with sigmaTot(E) method.

• alpha (float, optional) – Fractional buffer coefficient (default 0.1, i.e. 10%). For each pseudo-pole at energy E_pp, the buffer is max(alpha * |E_pp|, min_buffer). The fractional component reflects the scale-invariance of the asymptotic Sigma ~ E*X regime: pole “width” near the resolvent singularity scales with |E|, so a constant-eV buffer is overkill at deep poles and inadequate at moderate ones.

• min_buffer (float, optional) – Absolute floor for the buffer in eV (default 50). Prevents the fractional formula from placing Eminf arbitrarily close to a moderate-depth pseudo-pole. Set both alpha=0 and min_buffer=0 to disable buffering entirely (useful for tests of raw pole detection).

• E1 (float, optional) – Deep probe energies (eV) for asymptotic extraction. Default -1e3, -1e4.

• E2 (float, optional) – Deep probe energies (eV) for asymptotic extraction. Default -1e3, -1e4.

Returns:

Eminf_floor – Safe lower bound for Eminf in eV. Returns ENERGY_MIN from config if no pseudo-poles are detected.

Return type:

float

Notes

Buffer formula: per pseudo-pole, the gap is max(alpha * |E_pp|, min_buffer). The fractional term handles scale-invariance (deep poles need proportionally larger gaps); min_buffer prevents a meaninglessly tight gap near shallow poles. Defaults of 0.1 and 50 eV give ~300 eV gap at -3000 eV poles.

Positive pseudo-poles are filtered out: they cannot be excluded by Eminf (a LOWER bound) and cancel in calcTSW’s dTSW comparison.

A returned floor is GUARANTEED strictly negative (or exactly ENERGY_MIN). If the shallowest buffered floor is >= 0 (pseudo-pole mixed with physical states), returns ENERGY_MIN and lets calcTSW’s dTSW<0 path handle the system.

F, S, X, Sigma_0 are kept fully complex: taking Re() of a Hermitian matrix drops its antisymmetric imaginary part, giving a symmetric matrix with different eigenvalues – this matters for GHF/SOC.

gauNEGF.density.calcTSW(F, S, g, tol=0.001, maxN=10, Eminf=None, TSW=None, Emin_floor=-1000000.0)

Calculate integration bounds by stabilizing total spectral weight.

Iteratively doubles |Eminf| (Eminf <- 2*Eminf) until Tr(rho @ S) at the truncated contour matches a reference TSW computed at the configured Emin_floor. Uses a one-sided convergence test (TSW_ref - TSW_new < tol) so the loop only widens while the new contour is undercounting relative to the reference. dTSW < 0 (reference is smaller than the truncated contour) is the signature of a deep negative-weight pseudo-pole inside the reference contour, induced by asymptotically E-linear Sigma(E) from a non-orthogonal contact (see docs/tsw_convergence_notes.md); in that case the truncated contour is the physically correct choice and convergence is accepted to leave the pseudo-pole excluded.

The upper contour bound is held at -ENERGY_MIN (the wide config bound) independently of Emin_floor so that tightening Emin_floor (e.g. via calcPseudoPoleFloor) does not also chop the upper spectrum.

Parameters:

• F (ndarray) – Fock matrix in eV.

• S (ndarray) – Overlap matrix.

• g (surfG object) – Surface Green’s function calculator.

• tol (float, optional) – Convergence tolerance for relative TSW change (default: FERMI_CALCULATION_TOL).

• maxN (int, optional) – Maximum doubling iterations (default: FERMI_SEARCH_CYCLES).

• Eminf (float or None, optional) – Warm-start lower bound in eV. If None, initialized from min eigenvalue.

• TSW (float or None, optional) – Warm-start spectral weight. If None, recomputed from scratch using Emin_floor as the reference contour lower bound.

• Emin_floor (float, optional) – Hard lower limit for Eminf in eV (default: ENERGY_MIN from config). Also used as the reference contour lower bound when TSW is recomputed. If the doubling loop would expand past this limit without converging, the function warns and returns the current Eminf. Tighten this via calcPseudoPoleFloor to keep the reference contour from straddling deep pseudo-poles.

Returns:

(Eminf, TSW) – converged lower bound and reference spectral weight.

Return type:

tuple (float, float)

gauNEGF.density.damleCrossTerm(V, D, Vc, Y_eff, Q0, Q1, lo, hi)

Analytic lower-contour cross-term delta_N for the Damle method.

Reuses the eig (V, D, Vc) of Fbar already computed for the density matrix (no second eigendecomposition). Q(E) ~ Q0 + Q1*E is linearized between two anchors. The unphysical linear-Q ‘flat’ term b1*(hi-lo) is OMITTED – only the physical in-window-pole term is kept (see docs/lower_contour_math_and_probes.md sec 5):

delta_N = -(1/pi) Im sum_i (b0_i + b1_i*D_i)

• [log(1 - hi/D_i) - log(1 - lo/D_i)]

with b(E) = Vc^dagger Y_eff Q(E) Y_eff V, b0 = diag(…Q0…), b1 = diag(…Q1…). Poles outside [lo, hi] give a real log-difference (no contribution); in-window poles pick up the i*pi (the physical count).

Parameters:

• V (ndarray) – Eigendecomposition of Fbar: Fbar = V diag(D) Vc^dagger, Vc = inv(V^dagger).

• D (ndarray) – Eigendecomposition of Fbar: Fbar = V diag(D) Vc^dagger, Vc = inv(V^dagger).

• Vc (ndarray) – Eigendecomposition of Fbar: Fbar = V diag(D) Vc^dagger, Vc = inv(V^dagger).

• Y_eff (ndarray) – S_eff^(-1/2).

• Q0 (ndarray (N, N)) – Linear model of crossTermQTot: Q(E) ~ Q0 + Q1*E.

• Q1 (ndarray (N, N)) – Linear model of crossTermQTot: Q(E) ~ Q0 + Q1*E.

• lo (float) – Lower-contour energy limits in eV (lo = ENERGY_MIN, hi = Emin).

• hi (float) – Lower-contour energy limits in eV (lo = ENERGY_MIN, hi = Emin).

Returns:

Cross-term electron-count correction delta_N.

Return type:

float

gauNEGF.density.damleLowerDensity(F_eV, Y_eff, Sigma_0, g, Emin, ENERGY_MIN_=-1000000.0)

Analytic Damle lower-contour density on [ENERGY_MIN_, Emin].

Builds Fbar = Y_eff (F + Sigma_0 + i*ETA) Y_eff, diagonalizes it ONCE, and returns both the lower-contour density matrix (in the AO basis) and the analytic cross-term delta_N (which reuses the same eig; no second decomposition).

The density matrix is negated so its sign matches the densityComplex convention. The cross-term uses a linear model of crossTermQTot(E) anchored near the band bottom at (Emin, 2*Emin) and DROPS the unphysical linear-Q flat term (see damleCrossTerm).

Parameters:

• F_eV (ndarray (N, N)) – Device Fock matrix in eV.

• Y_eff (ndarray (N, N)) – S_eff^(-1/2) (from NEGFE._initAsymptoticSigma).

• Sigma_0 (ndarray (N, N)) – Asymptotic constant contact self-energy.

• g (surfG object) – Provides crossTermQTot(E) for the cross-term. If it returns None (orthogonal system), delta_N = 0.

• Emin (float) – Upper bound of the lower contour in eV (set by calcEmin upstream).

• ENERGY_MIN (float, optional) – Lower bound in eV (default config.ENERGY_MIN).

Returns:

(P_lower, delta_N).

Return type:

tuple (ndarray, float)

gauNEGF.density.density(V, Vc, D, Gam, Emin, mu)

Calculate density matrix using analytical integration for energy-independent self-energies.

Implements the analytical integration method from Eq. 27 in PRB 65, 165401 (2002). The method assumes energy-independent self-energies and uses the spectral function representation of the density matrix.

Parameters:

• V (ndarray) – Eigenvectors of Fock matrix in orthogonalized basis

• Vc (ndarray) – Inverse conjugate transpose of V

• D (ndarray) – Eigenvalues of Fock matrix

• Gam (ndarray) – Broadening matrix Γ = i[Σ - Σ†] in orthogonalized basis

• Emin (float) – Lower bound for integration in eV

• mu (float) – Chemical potential in eV

Returns:

Density matrix in orthogonalized basis

Return type:

ndarray

Notes

The integration is performed analytically using the residue theorem. The result includes contributions from poles below the Fermi energy.

gauNEGF.density.densityComplex(F, S, g, Emin, mu, tol=0.0001, T=0.0, debug=False)

Calculate equilibrium density matrix using complex contour integration.

Performs numerical integration along a complex contour that encloses the poles of the Fermi function. More efficient than real-axis integration for equilibrium calculations. Uses adaptive integration.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• Emin (float) – Lower bound for integration in eV

• mu (float) – Chemical potential in eV

• tol (float, optional) – Convergence tolerance (default: 1e-3)

• T (float, optional) – Temperature in Kelvin (default: 300)

• debug (bool, optional) – Whether to print debug information (default: False)

Returns:

(P, delta_N) where P is the equilibrium density matrix and delta_N is the Mulliken cross-term correction from device-lead overlap.

Return type:

tuple (ndarray, float)

Notes

Uses adaptive integration with GrIntCross to co-accumulate both the density matrix and cross-term scalar in a single pass. The cross-term converges at the same adaptive grid as the density matrix.

gauNEGF.density.densityComplexN(F, S, g, Emin, mu, N=100, T=0.0, showText=True, method='ant')

Calculate equilibrium density matrix using complex contour integration.

Performs numerical integration along a complex contour that encloses the poles of the Fermi function. More efficient than real-axis integration for equilibrium calculations. Uses vectorized integration for efficiency.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• Emin (float) – Lower bound for integration in eV

• mu (float) – Chemical potential in eV

• N (int, optional) – Number of integration points (default: 100)

• T (float, optional) – Temperature in Kelvin (default: 300)

• showText (bool, optional) – Whether to print progress messages (default: True)

• method ({'ant', 'legendre', 'chebyshev'}, optional) – Integration method to use (default: ‘ant’)

Returns:

(P, delta_N) where P is the equilibrium density matrix and delta_N is the Mulliken cross-term correction from device-lead overlap.

Return type:

tuple (ndarray, float)

Notes

The ‘ant’ method uses a modified Gauss-Chebyshev quadrature optimized for transport calculations, matching the ANT.Gaussian implementation.

gauNEGF.density.densityGrid(F, S, g, mu1, mu2, ind=None, tol=0.0001, T=0.0, debug=False)

Calculate non-equilibrium density matrix using real-axis integration.

Performs numerical integration for the non-equilibrium part of the density matrix when a bias voltage is applied. Uses ANT-modified Gauss-Chebyshev quadrature.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• mu1 (float) – Left contact chemical potential in eV

• mu2 (float) – Right contact chemical potential in eV

• ind (int or None, optional) – Contact index (None for total) (default: None)

• tol (float, optional) – Convergence tolerance (default: 1e-3)

• T (float, optional) – Temperature in Kelvin (default: 300)

• debug (bool, optional) – Whether to print debug information (default: False)

Returns:

Non-equilibrium contribution to density matrix

Return type:

ndarray

gauNEGF.density.densityGridN(F, S, g, mu1, mu2, ind=None, N=100, T=0.0, showText=True)

Calculate non-equilibrium density matrix using real-axis integration.

Performs numerical integration for the non-equilibrium part of the density matrix when a bias voltage is applied. Uses vectorized integration for efficiency.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• mu1 (float) – Left contact chemical potential in eV

• mu2 (float) – Right contact chemical potential in eV

• ind (int or None, optional) – Contact index (None for total) (default: None)

• N (int, optional) – Number of integration points (default: 100)

• T (float, optional) – Temperature in Kelvin (default: 300)

• showText (bool, optional) – Whether to print progress messages (default: True)

Returns:

Non-equilibrium contribution to density matrix

Return type:

ndarray

gauNEGF.density.densityGridTrap(F, S, g, mu1, mu2, ind=None, N=100, T=0.0)

Calculate non-equilibrium density matrix using trapezoidal integration on a real energy grid.

Alternative implementation to densityGridN() using trapezoidal rule instead of Gauss-Legendre quadrature. Provides direct loop-based integration for comparison.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• mu1 (float) – Left contact chemical potential in eV

• mu2 (float) – Right contact chemical potential in eV

• ind (int or None, optional) – Contact index (None for total) (default: None)

• N (int, optional) – Number of integration points (default: 100)

• T (float, optional) – Temperature in Kelvin (default: 300)

Returns:

Non-equilibrium contribution to density matrix

Return type:

ndarray

gauNEGF.density.densityReal(F, S, g, Emin, mu, tol=0.0001, T=0.0, debug=False)

Calculate equilibrium density matrix using adaptive real-axis ANT integration.

Uses integratePointsAdaptiveANT with the same Gauss-Chebyshev scheme as densityComplex, but integrating along the real axis instead of a complex contour. Suitable for zero-DOS regions (below all occupied states) where Im(G^R)=0 and the integral converges trivially in a small number of points.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• Emin (float) – Lower bound for integration in eV

• mu (float) – Chemical potential in eV

• tol (float, optional) – Convergence tolerance (default: 1e-3)

• T (float, optional) – Temperature in Kelvin (default: TEMPERATURE from gauNEGF.config)

• debug (bool, optional) – Whether to print per-iteration diagnostics (default: False)

Returns:

(P, delta_N) where P is the density matrix and delta_N is the Mulliken cross-term correction from device-lead overlap.

Return type:

tuple (ndarray, float)

gauNEGF.density.densityRealN(F, S, g, Emin, mu, N=100, T=0.0, showText=True)

Calculate equilibrium density matrix using real-axis integration on a specified grid.

Performs numerical integration along the real energy axis using Gauss-Legendre quadrature. Suitable for equilibrium calculations with energy-dependent self-energies.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• Emin (float) – Lower bound for integration in eV

• mu (float) – Chemical potential in eV

• N (int, optional) – Number of integration points (default: 100)

• T (float, optional) – Temperature in Kelvin (default: 300)

• showText (bool, optional) – Whether to print progress messages (default: True)

Returns:

(P, delta_N) where P is the density matrix and delta_N is the Mulliken cross-term correction from device-lead overlap.

Return type:

tuple (ndarray, float)

gauNEGF.density.fermi(E, mu, T)

Calculate Fermi-Dirac distribution.

Parameters:

• E (float) – Energy in eV

• mu (float) – Chemical potential in eV

• T (float) – Temperature in Kelvin

Returns:

Fermi-Dirac occupation number

Return type:

float

gauNEGF.density.getANTPoints(N)

Generate integration points and weights matching ANT.Gaussian implementation.

Follows the IntCompPlane subroutine in device.F90 from ANT.Gaussian package. Uses a modified Gauss-Chebyshev quadrature scheme optimized for transport calculations. _Note: Always generates an even number of points._

Parameters:

N (int) – Number of integration points

Returns:

(points, weights) - Arrays of integration points and weights

Return type:

tuple

gauNEGF.density.getFermiContact(g, ne, Emin=None, lBound=None, uBound=None, tol=0.0001, conv=0.001, maxcycles=10, T=0.0)

Calculate Fermi energy for a contact using adaptive integration.

Determines the Fermi energy for a contact system (3D lattice or 1D chain) by matching the electron count using bisection with adaptive complex contour integration.

Parameters:

• g (surfG object) – Surface Green’s function calculator

• ne (float) – Target number of electrons

• Emin (float, optional) – Lower bound for complex contour integration in eV. If None, calculated from DOS.

• lBound (float, optional) – Lower bound for bisection search in eV. If None, estimated from eigenvalues.

• uBound (float, optional) – Upper bound for bisection search in eV. If None, estimated from eigenvalues.

• tol (float, optional) – Tolerance for adaptive integration (default: ADAPTIVE_INTEGRATION_TOL)

• conv (float, optional) – Convergence tolerance for electron count (default: FERMI_CALCULATION_TOL)

• maxcycles (int, optional) – Maximum number of iterations (default: FERMI_SEARCH_CYCLES)

• T (float, optional) – Temperature in Kelvin (default: TEMPERATURE)

Returns:

Fermi energy in eV

Return type:

float

gauNEGF.density.integralFit(F, S, g, mu, Eminf=-1000000.0, tol=0.001, T=0.0, maxN=100)

Optimize integration parameters for density calculations.

Determines optimal integration parameters by iteratively testing convergence of the density matrix calculation.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• mu (float) – Equilibrium contact fermi energy in eV

• Eminf (float) – Lower bound for integration in eV (default: -1e6)

• tol (float, optional) – Convergence tolerance (default: 1e-5)

• T (float) – Temperature in Kelvin for Fermi broadening (default: 0)

• maxN (int, optional) – Max grid points and Emin search iterations(default: 1000)

Returns:

(Emin, N1, N2) - Optimized integration parameters: - Emin: Lower bound for complex contour - N1: Number of complex contour points - N2: Number of real axis points

Return type:

tuple

Notes

The optimization process: 1. Finds Emin by checking DOS convergence 2. Optimizes N1 for complex contour integration 3. Optimizes N2 for real axis integration

gauNEGF.density.integralFitNEGF(F, S, g, fermi, qV, Eminf=-1000000.0, tol=0.001, T=0.0, maxGrid=500)

Determines number of grid points for non-equilibrium density calculations.

Same procedure as integralFit() but applied to densityGrid()

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• fermi (float) – Equilibrium Fermi energy in eV

• qV (float) – Applied bias voltage in eV

• Eminf (float) – Lower bound for integration in eV (default: -1e6)

• tol (float, optional) – Convergence tolerance (default: 1e-5)

• T (float) – Temperature in Kelvin for Fermi broadening (default: 0)

• maxGrid (int, optional) – Maximum number of gridpoints (default: 1000)

Returns:

Number of grid points

Return type:

int

gauNEGF.density.integratePoints(computePointFunc, numPoints, parallel=False, numWorkers=None, chunkSize=None, debug=False)

Perform parallel or serial integration for quantum transport calculations.

This function provides a flexible integration framework that automatically chooses between numpy’s built-in parallelization for matrix operations and process-level parallelization for large workloads.

Parameters:

• computePointFunc (callable) – Function that computes a single integration point. Should take an integer index i and return a matrix/array.

• numPoints (int) – Total number of points to integrate

• parallel (bool, optional) – Whether to force process-level parallel processing (default: False)

• numWorkers (int, optional) – Number of worker processes. If None, automatically determined based on system resources and workload.

• chunkSize (int, optional) – Size of chunks for parallel processing. If None, automatically optimized based on numPoints and numWorkers.

• debug (bool, optional) – Whether to print debug information (default: False)

Returns:

Sum of all computed points

Return type:

ndarray

Notes

• Automatically detects SLURM environment for HPC compatibility

• Falls back to serial processing if parallel execution fails

• Uses numpy’s built-in parallelization for small workloads

• Switches to process-level parallelization for:

  • Large number of points (≥100)

  • Many available cores (≥32)

  • When explicitly requested via parallel=True

gauNEGF.density.integratePointsAdaptiveANT(computePoint, tol=0.0001, maxN=500, debug=False)

Adaptive integration using ANT-modified Gauss-Chebyshev quadrature (IntCompPlane subroutine from ANT.Gaussian package)

Parameters:

• computePoint (callable) – Function that computes integral over a list of weights and points. Should return a matrix/array.

• tol (float, optional) – Tolerance for the adaptive integration.

• maxN (int, optional) – Maximum number of points to integrate.

• debug (bool, optional) – Whether to print debug information.

Returns:

Integral of the function.

Return type:

ndarray

Transport Module

See also

Guide 5: Workflow Recipes – IV curve sweep, checkpointing, multi-temperature workflows

Transport calculations for quantum systems using Non-Equilibrium Green’s Functions.

This module provides functions for calculating quantum transport properties: - Coherent transmission through molecular junctions - Spin-dependent transport calculations - Current calculations at finite bias - Density of states calculations

The module supports both energy-independent and energy-dependent self-energies, with implementations for both spin-restricted and spin-unrestricted calculations. Spin-dependent transport follows the formalism described in [1].

References

[1]

Herrmann, C., Solomon, G. C., & Ratner, M. A. J. Chem. Theory Comput. 6, 3078 (2010) DOI: 10.1021/ct1006949

gauNEGF.transport.DOS(Elist, F, S, sig1, sig2)

Calculate density of states with energy-independent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate DOS at

• F (ndarray) – Fock matrix, NxN

• S (ndarray) – Overlap matrix, NxN

• sig1 (ndarray) – Left contact self-energy, vector (1xN) or matrix (NxN)

• sig2 (ndarray) – Right contact self-energy, vector (1xN) or matrix (NxN)

Returns:

(DOS, DOSList) where: - DOS: Total density of states at each energy - DOSList: Site-resolved DOS at each energy

Return type:

tuple

gauNEGF.transport.DOSE(Elist, F, S, g)

Calculate density of states with energy-dependent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate DOS at

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

Returns:

(DOS, DOSList) where: - DOS: Total density of states at each energy - DOSList: Site-resolved DOS at each energy

Return type:

tuple

class gauNEGF.transport.SigmaCalculator(sig1, sig2=None, energy_dependent=None)

Bases: object

Unified interface for energy-dependent and energy-independent self-energies.

Automatically detects whether sigma inputs are static arrays or energy-dependent surface Green’s function objects, and provides a consistent interface for both.

__init__(sig1, sig2=None, energy_dependent=None)

Initialize sigma calculator.

Parameters:

• sig1 (ndarray or surfG object) – Left contact self-energy (static) or surface Green’s function calculator

• sig2 (ndarray, optional) – Right contact self-energy (static arrays only)

• energy_dependent (bool, optional) – Force energy dependence detection. If None, auto-detect from sig1 type

get_Q_tot(E, spin=None, matrix_size=None)

Get cross-term Q_tot at energy E, with spin expansion matching get_sigma_total.

Returns None for energy-independent sigma or when surfG has no overlap coupling. Returns a JAX array of zeros_like(S) shape if Q is explicitly zero.

get_gamma(E, contact_index, spin=None, matrix_size=None)

Get gamma matrix (coupling) at energy E for specified contact.

get_sigma(E, contact_index, spin=None, matrix_size=None)

Get contact-specific self-energy at energy E.

get_sigma_total(E, spin=None, matrix_size=None)

Get total self-energy at energy E.

gauNEGF.transport.calculate_current(F, S, sigma_calculator, fermi, qV, T=0.0, spin=None, dE=0.001, **kwargs)

Calculate current at applied voltage with checkpointing.

Generates energy grid internally based on fermi, qV, T, and dE, then calculates transmission over integration window and integrates. Uses transmission checkpointing internally. The checkpoint file stores transmission data which can be reused.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sigma_calculator (SigmaCalculator) – Sigma calculator object

• fermi (float) – Fermi energy in eV

• qV (float) – Applied bias voltage in eV

• T (float, optional) – Temperature in Kelvin (default: TEMPERATURE from config)

• spin (str, optional) – Spin configuration (‘r’, ‘u’, ‘ro’, ‘g’). Defaults to ‘r’

• dE (float, optional) – Energy step for integration in eV (default: ENERGY_STEP from config)

• **kwargs – Additional parameters (such as checkpointing parameters) passed to calculate_transmission

Returns:

For restricted: current value (float) For open shell: (current_total (float), current_spin (4,))

Return type:

float or tuple

gauNEGF.transport.calculate_dos(F, S, sigma_calculator, energy_list, spin=None, checkpoint_file=None, checkpoint_interval=10)

Calculate DOS over energy range with checkpointing.

Uses -1 as placeholder for uncalculated energies. Saves checkpoint every checkpoint_interval energies. If checkpoint_file exists, resumes from first uncalculated energy.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sigma_calculator (SigmaCalculator) – Sigma calculator object

• energy_list (array_like) – List of energies in eV

• spin (str, optional) – Spin configuration (‘r’, ‘u’, ‘ro’, ‘g’). Defaults to ‘r’

• checkpoint_file (str, optional) – Path to checkpoint file (.npz format). If None, no checkpointing

• checkpoint_interval (int, optional) – Save checkpoint every N energies (default: 10)

Returns:

For restricted: (dos_total (N,), dos_per_site (N,M)) For open shell: (dos_total (N,), dos_per_site (N,M), dos_spin (N,2)) where M = F.shape[0] (number of sites)

Return type:

tuple

gauNEGF.transport.calculate_transmission(F, S, sigma_calculator, energy_list, spin=None, checkpoint_file=None, checkpoint_interval=10)

Calculate transmission over energy range with checkpointing.

Uses -1 as placeholder for uncalculated energies. Saves checkpoint every checkpoint_interval energies. If checkpoint_file exists, resumes from first uncalculated energy.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sigma_calculator (SigmaCalculator) – Sigma calculator object

• energy_list (array_like) – List of energies in eV

• spin (str, optional) – Spin configuration (‘r’, ‘u’, ‘ro’, ‘g’). Defaults to ‘r’

• checkpoint_file (str, optional) – Path to checkpoint file (.npz format). If None, no checkpointing

• checkpoint_interval (int, optional) – Save checkpoint every N energies (default: 10)

Returns:

For restricted: transmission array (N,) For open shell: (transmission (N,), spin_transmission (N,4))

Return type:

ndarray or tuple

gauNEGF.transport.cohTrans(Elist, F, S, sig1, sig2)

Calculate coherent transmission with energy-independent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate transmission at

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sig1 (ndarray) – Left contact self-energy (vector or matrix)

• sig2 (ndarray) – Right contact self-energy (vector or matrix)

Returns:

Transmission values at each energy

Return type:

list

Notes

Supports both vector and matrix self-energies. For vector self-energies, diagonal matrices are constructed internally.

gauNEGF.transport.cohTransE(Elist, F, S, g)

Calculate coherent transmission with energy-dependent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate transmission at

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

Returns:

Transmission values at each energy

Return type:

list

Notes

Uses the surface Green’s function calculator to compute energy-dependent self-energies at each energy point.

gauNEGF.transport.cohTransSpin(Elist, F, S, sig1, sig2, spin='u')

Calculate spin-dependent coherent transmission with energy-independent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate transmission at

• F (ndarray) – Fock matrix (2N x 2N for spin-unrestricted)

• S (ndarray) – Overlap matrix (2N x 2N for spin-unrestricted)

• sig1 (ndarray) –

  Left contact self-energy (spin independent vector (1xN) or matrix (NxN),

  spin dependent matrix (2Nx2N))

• sig2 (ndarray) –

  Right contact self-energy (spin independent vector (1xN) or matrix (NxN),

  spin dependent matrix (2Nx2N))

• spin (str, optional) – Spin basis {‘r’, ‘u’, ‘ro’, or ‘g’} (default: ‘u’)

Returns:

(Tr, Tspin) where: - Tr: Total transmission at each energy - Tspin: Array of spin-resolved transmissions [T↑↑, T↑↓, T↓↑, T↓↓]

Return type:

tuple

Notes

For collinear spin calculations (‘u’ or ‘ro’), the matrices are arranged in blocks: [F↑↑ 0 ] [S↑↑ 0 ] [0 F↓↓], [0 S↓↓] For generalized spin basis (‘g’), each orbital contains a 2x2 spinor block: [F↑↑ F↑↓] [S↑↑ S↑↓] [F↓↑ F↓↓], [S↓↑ S↓↓] which are then combined into a 2Nx2N matrix.

gauNEGF.transport.cohTransSpinE(Elist, F, S, g, spin='u')

Calculate spin-dependent coherent transmission with energy-dependent self-energies.

Parameters:

• Elist (array_like) – List of energies in eV to calculate transmission at

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• spin (str, optional) – Spin basis {‘r’, ‘u’, ‘ro’, or ‘g’} (default: ‘u’)

Returns:

(Tr, Tspin) where: - Tr: Total transmission at each energy - Tspin: Array of spin-resolved transmissions [T↑↑, T↑↓, T↓↑, T↓↓]

Return type:

tuple

gauNEGF.transport.current(F, S, sig1, sig2, fermi, qV, T=0.0, spin='r', dE=0.001)

Calculate coherent current using NEGF with energy-independent self-energies.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sig1 (ndarray) – Left contact self-energy (vector or matrix)

• sig2 (ndarray) – Right contact self-energy (vector or matrix)

• fermi (float) – Fermi energy in eV

• qV (float) – Applied bias voltage in eV

• T (float) – Temperature in Kelvin (default: TEMPERATURE from gauNEGF.config)

• spin (str, optional) – Spin configuration (‘r’ for restricted) (default: ‘r’)

• dE (float, optional) – Energy step for integration in eV (default: 0.01)

Returns:

Current in Amperes

Return type:

float

gauNEGF.transport.currentE(F, S, g, fermi, qV, T=0.0, spin='r', dE=0.001)

Calculate coherent current using NEGF with energy-dependent self-energies.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• g (surfG object) – Surface Green’s function calculator

• fermi (float) – Fermi energy in eV

• qV (float) – Applied bias voltage in eV

• T (float) – Temperature in Kelvin (default: TEMPERATURE from gauNEGF.config)

• spin (str, optional) – Spin configuration (‘r’ for restricted) (default: ‘r’)

• dE (float, optional) – Energy step for integration in eV (default: 0.01)

Returns:

Current in Amperes

Return type:

float

gauNEGF.transport.currentF(fn, dE=0.001, T=0.0)

Calculate current from saved SCF matrix file.

Parameters:

• fn (str) – Filename of .mat file containing SCF data

• dE (float, optional) – Energy step for integration in eV (default: 0.01)

• T (float, optional) – Temperature in Kelvin (default: TEMPERATURE from gauNEGF.config)

Returns:

Current in Amperes

Return type:

float

Notes

The .mat file should contain: - F: Fock matrix - S: Overlap matrix - sig1, sig2: Contact self-energies - fermi: Fermi energy - qV: Applied voltage - spin: Spin configuration

gauNEGF.transport.currentSpin(F, S, sig1, sig2, fermi, qV, T=0.0, spin='r', dE=0.001)

Calculate coherent spin current using NEGF with energy-independent self-energies.

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• sig1 (ndarray) – Left contact self-energy (vector or matrix)

• sig2 (ndarray) – Right contact self-energy (vector or matrix)

• fermi (float) – Fermi energy in eV

• qV (float) – Applied bias voltage in eV

• T (float) – Temperature in Kelvin (default: TEMPERATURE from gauNEGF.config)

• spin (str, optional) – Spin configuration for spin-dependent calculations (default: ‘r’)

• dE (float, optional) – Energy step for integration in eV (default: 0.01)

Returns:

Spin-currents (in Amperes) [I↑↑, I↑↓, I↓↑, I↓↓]

Return type:

list

gauNEGF.transport.dos_single_energy(E, F_jax, S_jax, sigma_calc, spin=None)

Calculate density of states at a single energy.

Parameters:

• E (float) – Energy in eV

• F_jax (jax.numpy.ndarray) – Fock matrix (already JAX array)

• S_jax (jax.numpy.ndarray) – Overlap matrix (already JAX array)

• sigma_calc (SigmaCalculator) – Sigma calculator object

• spin (str, optional) – Spin configuration using Gaussian labels (‘r’, ‘u’, ‘ro’, ‘g’) If None, defaults to ‘r’

Returns:

For ‘r’: (total_dos, dos_per_site) For ‘u’/’ro’/’g’: (total_dos, dos_per_site, dos_spin_up, dos_spin_down) where dos_per_site includes both spins, dos_spin_up/down are separate

Return type:

tuple

gauNEGF.transport.transmission_single_energy(E, F_jax, S_jax, sigma_calc, spin=None)

Calculate transmission at a single energy.

Parameters:

• E (float) – Energy in eV

• F_jax (jax.numpy.ndarray) – Fock matrix (already JAX array)

• S_jax (jax.numpy.ndarray) – Overlap matrix (already JAX array)

• sigma_calc (SigmaCalculator) – Sigma calculator object

• spin (str, optional) – Spin configuration using Gaussian labels: ‘r’ - restricted (closed shell) ‘u’ - unrestricted (open shell, block diagonal) ‘ro’ - restricted open (treated same as ‘u’) ‘g’ - generalized (non-collinear, spinor) If None, defaults to ‘r’

Returns:

For ‘r’: transmission value For ‘u’/’ro’/’g’: (total_transmission, [T_up_up, T_up_down, T_down_up, T_down_down])

Return type:

float or tuple

Contact Models

Bethe Lattice

See also

Bethe Lattice Contacts – Full Bethe lattice contact guide (Au, AuSOC, SOC workflows)

Surface Green’s function implementation for Bethe lattice contacts.

This module provides a Bethe lattice implementation for modeling semi-infinite metallic contacts in quantum transport calculations. It supports: - FCC [111] surface geometry with proper orbital symmetries - Slater-Koster parameterization for s, p, and d orbitals - Temperature-dependent calculations - Spin-restricted and unrestricted calculations

The implementation follows the ANT.Gaussian approach [1], using a minimal basis set with s, p, and d orbitals for each contact atom. The Bethe lattice model provides an efficient way to describe bulk metallic electrodes while maintaining proper orbital symmetries and electronic structure. This approach allows for accurate modeling of metal-molecule interfaces without the computational cost of explicit periodic boundary conditions.

References

[1] Jacob, D., & Palacios, J. J. (2011). Critical comparison of electrode models

in density functional theory based quantum transport calculations. The Journal of Chemical Physics, 134(4), 044118. DOI: 10.1063/1.3526044

class gauNEGF.surfGBethe.surfGB(F, S, contacts, bar, latFile='Au', spin='r', eta=1e-05, T=0.0)

Bases: object

Surface Green’s function calculator for Bethe lattice contacts.

This class implements the Bethe lattice approximation for modeling semi-infinite metallic contacts. It handles: - Contact geometry detection and setup - Slater-Koster parameter management - Energy-dependent self-energy calculations - Temperature effects - Spin configurations

The Bethe lattice model represents the contacts as a semi-infinite tree-like structure with proper coordination number and orbital symmetries matching those of bulk FCC metals. This provides an efficient way to compute surface Green’s functions and self-energies for the electrodes, as demonstrated by Jacob & Palacios in their 2011 paper [1].

Parameters:

• F (ndarray) – Fock matrix

• S (ndarray) – Overlap matrix

• contacts (list of lists) – List of atom indices for each contact

• bar (object) – Gaussian interface object containing basis set information

• latFile (str, optional) – Filename for Bethe lattice parameters (default: ‘Au’)

• spin (str, optional) – Spin configuration (‘r’ for restricted) (default: ‘r’)

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

• T (float, optional) – Temperature in Kelvin (default: 0)

cVecs

Normal vectors for each contact surface

Type:

list

latVecs

Lattice vectors for each contact

Type:

list

indsLists

Orbital indices for each contact atom

Type:

list

dirLists

Direction vectors for nearest neighbors

Type:

list

nIndLists

Nearest neighbor indices for each atom

Type:

list

gList

surfGBAt objects for each contact

Type:

list

References

[1] Jacob, D., & Palacios, J. J. (2011). Critical comparison of electrode models

in density functional theory based quantum transport calculations. The Journal of Chemical Physics, 134(4), 044118. DOI: 10.1063/1.3526044

constructMat(Mdict, dirCosines, SOC=False)

Construct hopping/overlap matrix using Slater-Koster formalism.

Builds a 9x9 matrix for s, p, and d orbital interactions based on the Slater-Koster two-center approximation. The matrix is first constructed assuming a [0,0,1] bond direction, then rotated to the given direction using direction cosines. When SOC=True, expands to 18x18 via kron(M, I2).

Parameters:

• Mdict (dict) – Dictionary of Slater-Koster parameters (sss, sps, pps, etc.)

• dirCosines (ndarray) – Array [l,m,n] of direction cosines for the bond

Returns:

9x9 matrix containing orbital interactions in the rotated frame

Return type:

ndarray

Notes

Matrix blocks: - [0,0]: s-s interaction - [0:4,0:4]: s-p block - [0:4,4:9]: s-d and p-d blocks - [4:9,4:9]: d-d block

crossTermQ(E, i, conv=1e-05)

Cross-term Q_sym for contact i in full device basis.

Mirrors surfGB.sigma: iterates atoms in contact i, calls gList[i].crossTermQSurf with per-atom active directions, assembles result, then applies de-orthonormalization (same as sigma).

crossTermQTot(E, conv=1e-05)

Total cross-term Q_sym from all contacts.

genNeighbors(plane_normal, first_neighbor)

Generate 12 nearest neighbor unit vectors for an FCC [111] surface.

Creates a list of unit vectors representing the 12 nearest neighbors in an FCC lattice: - 6 in-plane vectors forming a hexagonal pattern (3 pairs of opposite vectors) - 6 out-of-plane vectors forming triangular patterns (3 pairs of opposite vectors)

Parameters:

• plane_normal (ndarray) – Vector normal to the crystal plane (will be normalized)

• first_neighbor (ndarray) – Vector to one nearest neighbor (will be projected onto plane)

Returns:

12 unit vectors representing nearest neighbor directions

Return type:

list

getSigma(Elist=[None, None], conv=1e-05)

Helper method for getting the left and right contact self-energies

Parameters:

• Elist (tuple, optional) – A list of contact energies for selecting sigma, (default: use contact Fermi energy)

• conv (float, optional) – Convergence criterion for the self-energy matrix

Returns:

A tuple of both self-energy matrices (ndarrays)

Return type:

tuple

readBetheParams(filename)

Read Slater-Koster parameters from a .bethe file.

Reads and validates parameters for minimal basis with single s, p, and d orbitals. Parameters are stored in dictionaries for onsite energies, hopping integrals, and overlap matrices.

Parameters:

filename (str) – Name of the .bethe file (without extension)

Raises:

AssertionError – If parameters are missing or invalid

Notes

Parameters are sorted into: - Edict: Onsite energies (converted from Hartrees to eV) - Vdict: Hopping parameters (converted from Hartrees to eV) - Sdict: Overlap parameters

runAllTests()

Run all validation tests for surfGB.

Executes all test methods to validate: - d orbital angular functions - d orbital symmetry - p-d interactions - d-d interactions - General hopping physics

setF(F, muL, muR)

Update Fock matrix and contact chemical potentials.

Sets the Fock matrix and updates the Fermi levels of the left and right contacts if they have changed.

Parameters:

• F (ndarray) – New Fock matrix

• muL (float) – Chemical potential for left contact in eV

• muR (float) – Chemical potential for right contact in eV

sigma(E, i, conv=1e-05)

Calculate self-energy matrix for a specific contact.

Computes the self-energy matrix for contact i by: 1. Calculating surface self-energies for all 9 directions 2. Summing contributions from directions not connected to the device 3. Applying de-orthonormalization if needed 4. Handling spin configurations

Parameters:

• E (float) – Energy point for self-energy calculation (in eV)

• i (int) – Index of the contact to calculate self-energy for

• conv (float, optional) – Convergence criterion for self-energy calculation (default: SURFACE_GREEN_CONVERGENCE)

Returns:

Self-energy matrix for the specified contact, with dimensions: - (N, N) for restricted calculations - (2N, 2N) for unrestricted or generalized spin calculations

Return type:

ndarray

References

[1] Jacob, D., & Palacios, J. J. (2011). Critical comparison of electrode models

in density functional theory based quantum transport calculations. The Journal of Chemical Physics, 134(4), 044118. DOI: 10.1063/1.3526044

sigmaTot(E, conv=1e-05)

Calculate total self-energy matrix from all contacts.

Sums sigma(E, i) over all contacts and returns the result in the full device basis.

Parameters:

• E (float) – Energy point for Green’s function calculation (in eV)

• conv (float, optional) – Convergence criterion for self-energy calculation (default: 1e-5)

Returns:

Total self-energy matrix in the full device basis

Return type:

ndarray

testDDInteraction()

Test d-d orbital interactions.

Validates d-d orbital interactions by checking: - dyz-dyz interaction along x-axis (should be pure delta) - dz2-dz2 interaction along z-axis (should be pure sigma)

testDOrbitalFunctions()

Test d orbital angular functions.

Validates the angular dependence of d orbital interactions by checking: - dxy interaction along x-axis (should be zero) - dx2-y2 interaction along x-axis (should be sqrt(3)/2 * sds) - dz2 interaction along x-axis (should be -1/2 * sds)

testDOrbitalSymmetry()

Test d orbital symmetry properties.

Validates that d orbital interactions respect inversion symmetry by comparing interactions along opposite directions.

testHoppingPhysics()

Test physical properties of hopping matrices.

Validates hopping matrix physics by checking: - s-p hopping antisymmetry - Conservation of total s-p hopping magnitude - Proper angular dependence along principal axes and 45-degree rotations

testPDInteraction()

Test p-d orbital interactions.

Validates p-d orbital interactions by checking: - px-dxy interaction along x-axis (should be zero) - pz-dz2 interaction along z-axis (should be pure sigma)

updateFermi(i, Ef)

Update Fermi energy for a specific contact.

Shifts the Hamiltonian of contact i to align its Fermi level with the specified energy.

Parameters:

• i (int) – Contact index

• Ef (float) – New Fermi energy in eV

class gauNEGF.surfGBethe.surfGBAt(H, Slist, Vlist, eta, T=0.0, SOC=False)

Bases: object

Atomic-level Bethe lattice Green’s function calculator.

This class implements the surface Green’s function calculation for a single atom in the Bethe lattice, handling: - Onsite and hopping matrix construction - Self-energy calculations for bulk and surface - Temperature effects - Fermi energy optimization

Parameters:

• H (ndarray) – Onsite Hamiltonian matrix (9x9 for minimal basis)

• Slist (list of ndarray) – List of 12 overlap matrices for nearest neighbors

• Vlist (list of ndarray) – List of 12 hopping matrices for nearest neighbors

• eta (float) – Broadening parameter in eV

• T (float, optional) – Temperature in Kelvin (default: 0)

NN

Number of nearest neighbors (fixed to 12 for FCC)

Type:

int

fermi

Current Fermi energy

Type:

float

F

Extended Fock matrix including neighbors

Type:

ndarray

S

Extended overlap matrix including neighbors

Type:

ndarray

DOS(E)

Calculate surface density of states of the Bethe lattice.

Parameters:

E (float) – Energy point for DOS calculation (in eV)

Returns:

Density of states at energy E

Return type:

float

__init__(H, Slist, Vlist, eta, T=0.0, SOC=False)

Initialize surfGBAt with Hamiltonian and neighbor matrices.

Parameters:

• H (ndarray) – Onsite Hamiltonian matrix (9x9 for minimal basis)

• Slist (list of ndarray) – List of 12 overlap matrices for nearest neighbors

• Vlist (list of ndarray) – List of 12 hopping matrices for nearest neighbors

• eta (float) – Broadening parameter in eV

• T (float, optional) – Temperature in Kelvin (default: 0)

• SOC (bool, optional) – Whether to include spin-orbit coupling (default: False)

Raises:

AssertionError – If matrix dimensions are incorrect or number of neighbors != 12

calcFermi(ne, tol=0.001)

Calculate Fermi energy using bisection method.

Uses getFermiContact from density.py to find the Fermi energy that gives the correct number of electrons.

Parameters:

• ne (float) – Target number of electrons

• tol (float, optional) – Convergence tolerance (default: 1e-5)

Returns:

Calculated Fermi energy in eV

Return type:

float

crossTermQ(E, i, conv=1e-05)

Cross-term Q_sym for contact i (delegates to crossTermQBulk).

crossTermQBulk(E, conv=1e-05, mix=0.5)

Bulk cross-term Q_sym over all 12 directions.

Uses g_k = inv(A - sigTot + sigK[pair_k]) for each direction k, i.e. the neighbor’s Green’s function excluding the coupling back to the center atom.

crossTermQSurf(E, sigInds=None, conv=1e-05, mix=0.5)

Symmetrized cross-term Q_sym = sum_k (B_k g_k S_k + S_k g_k B_k^bar) / 2.

sigInds: list of surface direction indices to include (default: all 9). Uses the neighbor Green’s function g_k = inv(A - sigTot + sigSurf[pair_k]) for each direction k, not the center atom’s own Green’s function. For out-of-plane UP directions (3,4,5) whose pairs (9,10,11) are not surface directions, g_k = inv(A - sigTot) since the pair self-energy is already absent from sigTot.

crossTermQTot(E, conv=1e-05)

Total cross-term Q_sym (single contact = crossTermQ).

setF(F, mu1, mu2)

Empty function for compatibility with density.py methods.

Bethe lattice bulk properties are intrinsic (dependent on TB parameters).

Parameters:

• F (ndarray) – Fock matrix (unused)

• mu1 (float) – First chemical potential (unused)

• mu2 (float) – Second chemical potential (unused)

sigma(E, i, conv=1e-05)

Self-energy for contact i (only i=0, single bulk contact).

sigmaK(E, conv=1e-05, mix=0.5)

Calculate bulk self-energies for all 12 lattice directions.

Computes self-energies for an FCC lattice with the following geometry:

[3x out of plane dir]

|/

[3x plane dir] - o - [3x plane dir]

/|

[3x out of plane dir]

Uses a self-consistent iteration scheme with mixing to solve the Dyson equation.

Parameters:

• E (float) – Energy point for Green’s function calculation (in eV)

• conv (float, optional) – Convergence criterion for Dyson equation (default: 1e-5)

• mix (float, optional) – Mixing factor for Dyson equation (default: 0.5)

Returns:

Array of 12 self-energy matrices (9x9 each) in order by lattice direction

Return type:

ndarray

Notes

Uses previous solution as initial guess when energy point is close to previous calculation to improve convergence.

sigmaSurf(E, conv=1e-05, mix=0.5)

Calculate surface self-energies for an FCC lattice.

Computes self-energies for atoms at the surface with the geometry: [3x plane dir] - o - [3x plane dir]

/|

[3x out of plane dir]

Uses a self-consistent iteration scheme with mixing to solve the Dyson equation. The implementation follows the Bethe lattice approach described in Jacob & Palacios (2011), where the self-energy is computed recursively for a semi-infinite tree-like structure that preserves the proper coordination number and orbital symmetries of bulk FCC metals.

Parameters:

• E (float) – Energy point for Green’s function calculation (in eV)

• conv (float, optional) – Convergence criterion for Dyson equation (default: 1e-5)

• mix (float, optional) – Mixing factor for Dyson equation (default: 0.5)

Returns:

List of self-energy matrices for the surface atom.

Return type:

list

Notes

First calculates bulk self-energies using sigmaK, then iterates to find surface self-energies for the 9 surface directions. The recursive method ensures proper treatment of the metal-molecule interface while maintaining computational efficiency.

References

[1] Jacob, D., & Palacios, J. J. (2011). Critical comparison of electrode models

in density functional theory based quantum transport calculations. The Journal of Chemical Physics, 134(4), 044118. DOI: 10.1063/1.3526044

sigmaTot(E, conv=1e-05)

updateH(fermi=None)

Update Hamiltonian and Fock matrix.

Sets F = H (9x9) and S = I for protocol compatibility with density.py.

Parameters:

fermi (float, optional) – New Fermi energy setpoint in eV (default: None)

1D Chain

See also

1D Chain Contacts in gauNEGF – 1D chain contact setup guide (three usage patterns)

Surface Green’s function implementation for 1D chain contacts.

This module provides a surface Green’s function calculator for quasi-1D chain contacts in quantum transport calculations. It supports three usage patterns for specifying contact parameters: fully automatic extraction from Fock/Overlap matrices, custom coupling matrices with automatic onsite parameters, or fully specified contact parameters. The implementation uses iterative solvers with JAX JIT compilation for efficient computation.

The surfG class handles: - Semi-infinite 1D chain Green’s function calculations - Multiple contact geometries and coupling specifications - Orthogonal and non-orthogonal contact overlap matrices - De-orthonormalization for proper self-energy contributions - Contact regularization via congruent eigenvalue clipping - Self-energy and cross-term calculations with convergence control

Key exports: - surfG: Main surface Green’s function calculator class

class gauNEGF.surfG1D.surfG(Fock, Overlap, indsList, taus=None, staus=None, alphas=None, aOverlaps=None, betas=None, bOverlaps=None, eta=1e-05, spin='r')

Bases: object

Surface Green’s function calculator for 1D chain contacts.

This class implements the surface Green’s function calculation for 1D chain contacts. It supports three usage patterns:

1. Fully automatic extraction from Fock matrix: - Provide contact indices and connection indices - All parameters extracted from F/S matrices Example: surfG1D(F, S, [[c1], [c2]], [[c1conn], [c2conn]])

2. Fock matrix with custom coupling: - Provide contact indices and coupling matrices - Onsite contact parameters from F/S, coupling specified manually - If staus=None (default), assumes orthonormal coupling Example: surfG1D(F, S, [[c1], [c2]], [tau1, tau2])

3. Fully specified contacts: - All contact parameters provided manually - If aOverlaps=None, onsite overlap defaults to identity (orthonormal) - If bOverlaps=None, hopping overlap defaults to zeros (orthonormal) Example: surfG1D(F, S, [[c1], [c2]], [tau1, tau2], [stau1, stau2],

   [alpha1, alpha2], [salpha1, salpha2], [beta1, beta2], [sbeta1, sbeta2])

Parameters:

• Fock (ndarray) – Fock matrix for the extended system

• Overlap (ndarray) – Overlap matrix for the extended system

• indsList (list of lists) – Lists of orbital indices for each contact region

• taus (list or None, optional) – Either coupling matrices or connection indices (default: None) - If indices: [[contact1connection], [contact2connection]] - If matrices: [tau1, tau2]

• staus (list or None, optional) – Overlap matrices for coupling (default: None = orthonormal coupling) None entries trigger de-orthonormalization in sigma()

• alphas (list of ndarray or None, optional) – On-site energies for contacts, required for pattern (c) (default: None)

• aOverlaps (list of ndarray or None, optional) – On-site overlap matrices; defaults to identity when None (orthonormal)

• betas (list of ndarray or None, optional) – Hopping matrices between contact unit cells, required for pattern (c) (default: None)

• bOverlaps (list of ndarray or None, optional) – Overlap matrices between contact unit cells; defaults to zeros when None

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

F

Fock matrix

Type:

ndarray

S

Overlap matrix

Type:

ndarray

X

Inverse square root of overlap matrix for orthogonalization (S^-0.5)

Type:

ndarray

Xi

Square root of overlap matrix for de-orthonormalization (S^+0.5 = inv(X))

Type:

ndarray

tauList

Contact coupling matrices

Type:

list

stauList

Contact coupling overlap matrices (None entries -> orthonormal coupling)

Type:

list

aList

On-site energy matrices for contacts

Type:

list

aSList

On-site overlap matrices for contacts

Type:

list

bList

Hopping matrices between contact unit cells

Type:

list

bSList

Overlap matrices between contact unit cells

Type:

list

gPrev

Previous surface Green’s functions for convergence

Type:

list

__init__(Fock, Overlap, indsList, taus=None, staus=None, alphas=None, aOverlaps=None, betas=None, bOverlaps=None, eta=1e-05, spin='r')

Initialize the surface Green’s function calculator.

The initialization follows one of three patterns: a) Fully automatic: Only provide Fock, Overlap, indsList, and connection indices in taus b) Custom coupling: Provide Fock, Overlap, indsList, coupling matrices in taus

• staus=None (default) means orthonormal coupling -> de-ortho applied in sigma()

3. Fully specified: Provide all parameters including alphas, aOverlaps, betas, bOverlaps - aOverlaps=None defaults to identity (orthonormal onsite) - bOverlaps=None defaults to zeros (orthonormal hopping)

Parameters:

• Fock (ndarray) – Fock matrix for the extended system

• Overlap (ndarray) – Overlap matrix for the extended system

• indsList (list of lists) – Lists of orbital indices for each contact region

• taus (list or None, optional) – Either coupling matrices or connection indices (default: None) - If indices: [[contact1connection], [contact2connection]] - If matrices: [tau1, tau2]

• staus (list or None, optional) – Overlap matrices for coupling (default: None = orthonormal)

• alphas (list of ndarray or None, optional) – On-site energies for contacts, required for pattern (c) (default: None)

• aOverlaps (list of ndarray or None, optional) – On-site overlap matrices for contacts (default: None = identity)

• betas (list of ndarray or None, optional) – Hopping matrices between contact unit cells, required for pattern (c) (default: None)

• bOverlaps (list of ndarray or None, optional) – Overlap matrices between contact unit cells (default: None = zeros)

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

• spin (str, optional) – Spin configuration (‘r’ for restricted) (default: ‘r’)

crossTermQ(E, i, conv=1e-05)

Symmetrized cross-term matrix Q_sym_i in full device basis.

Q_sym = (t_eff @ g_surf @ S_LD + S_DL @ g_surf @ bar_t_eff) / 2

where t_eff is the regularized tau_DL and bar_t_eff is the regularized tau_LD (EOM coupling with unconjugated E, NOT t_eff^dagger). Returns None if contact i has orthogonal coupling (stauList[i] is None).

crossTermQTot(E, conv=1e-05)

Sum of Q_sym over all contacts. Returns None if all contacts orthogonal.

g(E, i, conv=1e-05, relFactor=0.5)

Calculate surface Green’s function for a contact.

Uses an iterative scheme to calculate the surface Green’s function for contact i at energy E. The iteration continues until the change in the Green’s function is below the convergence criterion.

Parameters:

• E (float) – Energy point in eV

• i (int) – Contact index (static argument for JAX JIT compilation)

• conv (float, optional) – Convergence criterion for iteration (default: 1e-5)

• relFactor (float, optional) – Relaxation factor for iteration mixing (default: 0.5)

Returns:

Surface Green’s function matrix for contact i

Return type:

ndarray

setF(F, mu1=None, mu2=None)

Update the Fock matrix and contact chemical potentials.

This method updates the system’s Fock matrix and optionally shifts the contact chemical potentials. If the contacts are extracted from the Fock matrix, their parameters are automatically updated.

Parameters:

• F (ndarray) – New Fock matrix for the system

• mu1 (float or None, optional) – Chemical potential for first contact in eV (default: None)

• mu2 (float or None, optional) – Chemical potential for second contact in eV (default: None)

sigma(E, i, conv=1e-05)

Calculate self-energy matrix for a contact.

Computes the self-energy matrix for contact i at energy E using the surface Green’s function. The self-energy represents the effect of the semi-infinite contact on the device region.

When stauList[i] is None (orthonormal coupling), applies de-orthonormalization: sig -> Xi[inds,inds] @ sig @ Xi[inds,inds] where Xi = S^+0.5 = inv(X).

Parameters:

• E (float) – Energy point in eV

• i (int (static)) – Contact index

• conv (float, optional) – Convergence criterion for surface Green’s function (default: 1e-5)

Returns:

Self-energy matrix for contact i

Return type:

ndarray

sigmaTot(E, conv=1e-05)

Calculate total self-energy matrix from all contacts.

Computes the total self-energy matrix at energy E by summing contributions from all contacts. This represents the combined effect of all semi-infinite contacts on the device region.

Parameters:

• E (float) – Energy point in eV

• conv (float, optional) – Convergence criterion for surface Green’s functions (default: 1e-5)

Returns:

Total self-energy matrix from all contacts

Return type:

ndarray

3D Contacts

class gauNEGF.surfG3D.surfG3(F, S, contacts, bar, latFile='Au', spin='r', eta=1e-05, T=0.0)

Bases: object

Surface Green’s function calculator for 3D lattice with [111] surface.

Parameters:

• F (ndarray) – Fock matrix from DFT calculation

• S (ndarray) – Overlap matrix from DFT calculation

• contacts (list of lists) – Lists of atom indices for each contact region

• bar (QCBinAr) – Gaussian interface object containing geometry and orbital information

• latFile (str, optional) – Name of .bethe file containing Slater-Koster parameters (default: ‘Au’)

• spin ({'r', 'u', 'ro', 'g'}, optional) – Spin treatment: restricted, unrestricted, or generalized (default: ‘r’)

• eta (float, optional) – Broadening parameter in eV (default: 1e-9)

• T (float, optional) – Temperature in Kelvin (default: 0)

F

Fock matrix

Type:

ndarray

S

Overlap matrix

Type:

ndarray

gList

List of atomic surface Green’s function calculators for each contact

Type:

list of surfGAt3D

constructMat(Mdict, dirCosines)

Construct hopping/overlap matrix using Slater-Koster formalism.

Builds a 9x9 matrix for s, p, and d orbital interactions based on the Slater-Koster two-center approximation. The matrix is first constructed assuming a [0,0,1] bond direction, then rotated to the given direction using direction cosines.

Parameters:

• Mdict (dict) – Dictionary of Slater-Koster parameters (ssσ, spσ, ppσ, etc.)

• dirCosines (ndarray) – Array [l,m,n] of direction cosines for the bond

Returns:

9x9 matrix containing orbital interactions in the rotated frame

Return type:

ndarray

Notes

Matrix blocks: - [0,0]: s-s interaction - [0:4,0:4]: s-p block - [0:4,4:9]: s-d and p-d blocks - [4:9,4:9]: d-d block

crossTermQ(E, i, conv=0.0001)

Cross-term Q_sym for contact i in full device basis.

Mirrors surfG3.sigma: computes G_AB via gSurf, then for each atom in contact i assembles Q using active directions, applies Xi.

crossTermQTot(E, conv=0.0001)

Total cross-term Q_sym from all contacts.

genNeighbors(plane_normal, first_neighbor)

Generate 12 nearest neighbor unit vectors for an FCC [111] surface.

Creates a list of unit vectors representing the 12 nearest neighbors in an FCC lattice: - 6 in-plane vectors forming a hexagonal pattern (3 pairs of opposite vectors) - 6 out-of-plane vectors forming triangular patterns (3 pairs of opposite vectors)

Parameters:

• plane_normal (ndarray) – Vector normal to the crystal plane (will be normalized)

• first_neighbor (ndarray) – Vector to one nearest neighbor (will be projected onto plane)

Returns:

12 unit vectors representing nearest neighbor directions

Return type:

list

getSigma(Elist=[None, None], conv=0.0001)

Helper method for getting the left and right contact self-energies

Parameters:

• Elist (tuple, optional) – A list of contact energies for selecting sigma, (default: use contact Fermi energy)

• conv (float, optional) – Convergence criterion for the self-energy matrix

Returns:

A tuple of both self-energy matrices (ndarrays)

Return type:

tuple

readBetheParams(filename)

Read Slater-Koster parameters from a .bethe file.

Reads and validates parameters for minimal basis with single s, p, and d orbitals. Parameters are stored in dictionaries for onsite energies, hopping integrals, and overlap matrices.

Parameters:

filename (str) – Name of the .bethe file (without extension)

Raises:

AssertionError – If parameters are missing or invalid

Notes

Parameters are sorted into: - Edict: Onsite energies (converted from Hartrees to eV) - Vdict: Hopping parameters (converted from Hartrees to eV) - Sdict: Overlap parameters

runAllTests()

Run all validation tests for surfG.

Executes all test methods to validate: - d orbital angular functions - d orbital symmetry - p-d interactions - d-d interactions - General hopping physics

setF(F, muL, muR)

Update Fock matrix and contact chemical potentials.

Sets the Fock matrix and updates the Fermi levels of the left and right contacts if they have changed.

Parameters:

• F (ndarray) – New Fock matrix

• muL (float) – Chemical potential for left contact in eV

• muR (float) – Chemical potential for right contact in eV

sigma(E, i, conv=0.0001)

Calculate self-energy matrix for a specific contact.

Computes the self-energy matrix for contact i by: 1. Calculating surface self-energies for all 9 directions 2. Summing contributions from directions not connected to the device 3. Applying de-orthonormalization if needed 4. Handling spin configurations

Parameters:

• E (float) – Energy point for self-energy calculation (in eV)

• i (int) – Index of the contact to calculate self-energy for

• conv (float, optional) – Convergence criterion for self-energy calculation (default: 1e-4)

Returns:

Self-energy matrix for the specified contact, with dimensions: - (N, N) for restricted calculations - (2N, 2N) for unrestricted or generalized spin calculations

Return type:

ndarray

References

[1] Jacob, D., & Palacios, J. J. (2011). Critical comparison of electrode models

in density functional theory based quantum transport calculations. The Journal of Chemical Physics, 134(4), 044118. DOI: 10.1063/1.3526044

sigmaTot(E, conv=0.0001)

Calculate total self-energy matrix from all contacts.

Sums sigma(E, i) over all contacts and returns the result in the full device basis.

Parameters:

• E (float) – Energy point for Green’s function calculation (in eV)

• conv (float, optional) – Convergence criterion for self-energy calculation (default: 1e-5)

Returns:

Total self-energy matrix in the full device basis

Return type:

ndarray

testDDInteraction()

Test d-d orbital interactions.

Validates d-d orbital interactions by checking: - dyz-dyz interaction along x-axis (should be pure delta) - dz2-dz2 interaction along z-axis (should be pure sigma)

testDOrbitalFunctions()

Test d orbital angular functions.

Validates the angular dependence of d orbital interactions by checking: - dxy interaction along x-axis (should be zero) - dx2-y2 interaction along x-axis (should be sqrt(3)/2 * sds) - dz2 interaction along x-axis (should be -1/2 * sds)

testDOrbitalSymmetry()

Test d orbital symmetry properties.

Validates that d orbital interactions respect inversion symmetry by comparing interactions along opposite directions.

testHoppingPhysics()

Test physical properties of hopping matrices.

Validates hopping matrix physics by checking: - s-p hopping antisymmetry - Conservation of total s-p hopping magnitude - Proper angular dependence along principal axes and 45-degree rotations

testPDInteraction()

Test p-d orbital interactions.

Validates p-d orbital interactions by checking: - px-dxy interaction along x-axis (should be zero) - pz-dz2 interaction along z-axis (should be pure sigma)

updateFermi(i, Ef)

Update Fermi energy for a specific contact.

Shifts the Hamiltonian of contact i to align its Fermi level with the specified energy.

Parameters:

• i (int) – Contact index

• Ef (float) – New Fermi energy in eV

class gauNEGF.surfG3D.surfGAt3D(H, Slist, Vlist, vecs, eta, T=0.0, kPoints=11)

Bases: object

Atomic-level 3D Green’s function calculator for a single atom.

This class implements the surface Green’s function calculation for a single atom in the 3D lattice with [111] surface, handling: - Onsite and hopping matrix construction - Self-energy calculations for bulk and surface

Parameters:

• H (ndarray) – Onsite Hamiltonian matrix (9x9 for minimal basis)

• Slist (list of ndarray) – List of 12 overlap matrices for nearest neighbors

• Vlist (list of ndarray) – List of 12 hopping matrices for nearest neighbors

• eta (float) – Broadening parameter in eV

• T (float, optional) – Temperature in Kelvin (default: 0)

NN

Number of nearest neighbors (fixed to 12 for FCC)

Type:

int

fermi

Current Fermi energy

Type:

float

F

Extended Fock matrix including neighbors

Type:

ndarray

S

Extended overlap matrix including neighbors

Type:

ndarray

DOS(E, conv=0.0001, mix=0.1)

Use surface Green’s function to calculate density of states.

Parameters:

• E (float) – Energy point for DOS calculation (in eV)

• conv (float, optional) – Convergence criterion for self-energy calculation (default: 1e-5)

• mix (float, optional) – Mixing parameter for self-energy convergence (default: 0.5)

Returns:

Density of states at energy E

Return type:

float

__init__(H, Slist, Vlist, vecs, eta, T=0.0, kPoints=11)

Initialize surfGAt3D with Hamiltonian and neighbor matrices.

Parameters:

• H (ndarray) – Onsite Hamiltonian matrix (9x9 for minimal basis)

• Slist (list of ndarray) – List of 12 overlap matrices for nearest neighbors

• Vlist (list of ndarray) – List of 12 hopping matrices for nearest neighbors, first six in plane,

• eta (float) – Broadening parameter in eV (recommended: >= 1e-4 for numerical stability)

• T (float, optional) – Temperature in Kelvin (default: 0)

• kPoints (int, optional) – Number of k-points per direction for BZ integration (default: 11) Recommended values: - 11: Fast, reasonable accuracy (121 points total) - 15: Good balance of speed/accuracy (225 points) - 21: High accuracy for production (441 points) Note: Odd values use Monkhorst-Pack grid, even use Gamma-centered

Raises:

AssertionError – If matrix dimensions are incorrect or number of neighbors != 12

calcFermi(ne, tol=1e-05)

Calculate Fermi energy using bisection method.

Uses getFermiContact from density.py to find the Fermi energy that gives the correct number of electrons.

Parameters:

• ne (float) – Target number of electrons

• tol (float, optional) – Convergence tolerance (default: 1e-5)

Returns:

Calculated Fermi energy in eV

Return type:

float

crossTermQ(E, i)

Cross-term Q_sym for contact i (delegates to crossTermQBulk).

Parameters:

• E (float) – Energy point in eV

• i (int) – Contact index (must be 0)

Returns:

Symmetrized cross-term matrix (shape: dim x dim)

Return type:

ndarray

crossTermQBulk(E)

Bulk cross-term Q_sym using G_AB propagator (dim x dim).

Computes the symmetrized cross-term for Mulliken population correction:

Q_sym = (tau @ G_AB @ S_LD + S_DL @ G_AB @ bar_tau) / 2

where tau and bar_tau use all 12 neighbor directions.

Parameters:

E (float) – Energy point in eV, pre-shifted by the caller to remove dFermi.

Returns:

Symmetrized cross-term matrix (shape: dim x dim)

Return type:

ndarray

crossTermQSurf(E, active_dirs=None, conv=0.0001, mix=0.1, G_AB=None)

Symmetrized cross-term Q_sym using G_AB propagator.

Q_sym = (tau @ G_sub @ S_LD + S_DL @ G_sub @ tau^dagger) / 2

where tau is the same phase-free coupling used in sigma().

crossTermQTot(E)

Total cross-term Q_sym (single bulk contact).

Parameters:

E (float) – Energy point in eV

Returns:

Symmetrized cross-term matrix (shape: dim x dim)

Return type:

ndarray

gBulk(E)

Calculate bulk Green’s function with full 3D periodicity via direct inversion.

For bulk with all neighbors included, no Dyson equation needed: g(k) = [(E + iη)S(k) - H(k)]^-1

Uses 3D k-space mesh with full periodicity in all directions.

Parameters:

E (float) – Energy point for Green’s function calculation (in eV)

Returns:

Real-space G_AB propagator matrix (12*dim x 12*dim). Block (A,B) at G_AB[A*dim:(A+1)*dim, B*dim:(B+1)*dim] gives the propagator G(R_A - R_B) between directions A and B.

Return type:

ndarray

gSurf(E, conv=0.0001, mix=0.1, maxIter=5000)

generate_band_plot(E_fermi=0.0, n_points=40, plot=True, save_path=None)

Generate band structure plot along high-symmetry path for FCC.

Uses the pre-computed reciprocal lattice vectors (b1_3D, b2_3D, b3_3D) to convert fractional k-coordinates to Cartesian k-space.

Parameters:

• E_fermi (float, optional) – Fermi energy to shift bands relative to (default: 0.0)

• n_points (int, optional) – Number of points between each high-symmetry point (default: 40)

• plot (bool, optional) – Whether to generate matplotlib plot (default: True)

• save_path (str, optional) – Path to save plot (default: None, display only)

Returns:

Dictionary containing: - ‘distances’: 1D array of k-path distances - ‘bands’: 2D array (n_kpoints x n_bands) of eigenvalues - ‘k_labels’: List of high-symmetry point labels - ‘k_positions’: Positions of high-symmetry points along path

Return type:

dict

setF(F, mu1, mu2)

Empty function for compatibility with density.py methods.

Bethe lattice bulk properties are intrinsic (dependent on TB parameters).

Parameters:

• F (ndarray) – Fock matrix (unused)

• mu1 (float) – First chemical potential (unused)

• mu2 (float) – Second chemical potential (unused)

sigma(E, i)

Self-energy for contact i (only i=0, single bulk contact).

Parameters:

• E (float) – Energy point in eV

• i (int) – Contact index (must be 0)

Returns:

Self-energy matrix (shape: dim x dim)

Return type:

ndarray

sigmaBulk(E)

Total bulk self-energy for the center atom (dim x dim).

Uses G_AB propagator with all 12 neighbors:

Sigma = tau @ G_AB @ bar_tau

where tau = [tau_0 | tau_1 | … | tau_11] is the horizontal concatenation of all 12 phase-free coupling matrices.

Parameters:

E (float) – Energy point for self-energy calculation (in eV), pre-shifted by the caller to remove dFermi.

Returns:

Total self-energy matrix (shape: dim x dim)

Return type:

ndarray

sigmaSurf(E, active_dirs=None, conv=0.0001, mix=0.1, G_AB=None)

Calculate surface self-energy using G_AB propagator.

Computes Sigma = tau @ G_S @ tau’ where: - G_S is the sub-block of G_AB for active directions - tau = [tau_0 | tau_1 | …] is the horizontal concatenation

of phase-free coupling matrices tau_a = (E+i*eta)*S_a - V_a

Parameters:

• E (float) – Energy point for Green’s function calculation (in eV)

• active_dirs (list of int, optional) – Direction indices to include. Default None means all 9 surface directions [0..8].

• conv (float, optional) – Convergence criterion for Dyson equation (default: 1e-4)

• mix (float, optional) – Mixing factor for Dyson equation (default: 0.1)

• G_AB (ndarray, optional) – Pre-computed G_AB propagator from gSurf(). If None, gSurf() is called internally.

Returns:

Self-energy matrix of shape (dim, dim).

Return type:

ndarray

sigmaTot(E)

Total self-energy for the center atom (dim x dim).

Delegates to sigmaBulk after removing the Fermi shift so that H0/Vlist0 (immutable reference frame) are used.

Parameters:

E (float) – Energy point for self-energy calculation (in eV)

Returns:

Total self-energy matrix (shape: dim x dim)

Return type:

ndarray

updateH(fermi=None)

Update Hamiltonian and protocol matrices F/S.

Updates onsite and hopping matrices. F and S are the 9x9 single-atom Hamiltonian and identity overlap, compatible with density.py and the SurfGProtocol interface. Self-energies from bulk neighbors are provided via sigmaTot()/sigma().

Parameters:

fermi (float, optional) – New Fermi energy setpoint in eV (default: None)

Notes

When fermi is provided and different from current value: - Shifts onsite energies by the Fermi level difference - Updates hopping matrices with overlap contributions

Constant Self Energy

Constant self-energy module for NEGF calculations.

This module provides energy-independent self-energies for NEGF calculations. While useful for testing and validation, it also serves production calculations where constant self-energies are appropriate, particularly for non-zero temperature calculations where energy-dependent effects may be less critical or when computational efficiency is prioritized over full energy dependence.

class gauNEGF.surfGTester.surfGTest(Fock, Overlap, indsList, sig1=None, sig2=None, spin='r')

Bases: object

Energy-independent self-energy calculator.

This class provides constant (energy-independent) self-energies for NEGF calculations. While useful for testing and validation, it also enables production calculations with constant self-energies, which can be appropriate for non-zero temperature calculations or when computational efficiency is prioritized. It implements the same interface as other surface Green’s function calculators.

Parameters:

• Fock (ndarray) – Fock matrix for the system

• Overlap (ndarray) – Overlap matrix for the system

• indsList (list) – List containing [left_contact_indices, right_contact_indices]

• sig1 (complex or array-like, optional) – Self-energy for left contact (default: None)

• sig2 (complex or array-like, optional) – Self-energy for right contact (default: None)

F

Fock matrix

Type:

ndarray

S

Overlap matrix

Type:

ndarray

X

Inverse square root of overlap matrix

Type:

ndarray

N

Size of the system matrices

Type:

int

indsList

Contact orbital indices

Type:

list

sig

List of self-energy matrices [left, right]

Type:

list

Notes

If no self-energies are provided, defaults to -0.05j on contact orbitals. This class is compatible with scfE.py for temperature-dependent calculations where constant self-energies provide a good approximation. Many method parameters are unused but maintained for interface consistency with other surface Green’s function calculators.

__init__(Fock, Overlap, indsList, sig1=None, sig2=None, spin='r')

Initialize constant self-energy calculator.

Parameters:

• Fock (ndarray) – Fock matrix for the system

• Overlap (ndarray) – Overlap matrix for the system

• indsList (list) – List containing [left_contact_indices, right_contact_indices]

• sig1 (complex or array-like, optional) – Self-energy for left contact (default: None)

• sig2 (complex or array-like, optional) – Self-energy for right contact (default: None)

crossTermQ(E, i, conv=1e-05)

Cross-term Q_sym for surfGTest. Always None (orthogonal basis).

crossTermQTot(E, conv=1e-05)

Total cross-term Q_sym for surfGTest. Always None (orthogonal basis).

setF(F, mu1, mu2)

Update Fock matrix and chemical potentials.

Parameters:

• F (ndarray) – New Fock matrix

• mu1 (float) – Chemical potential for left contact (unused - kept for API compatibility)

• mu2 (float) – Chemical potential for right contact (unused - kept for API compatibility)

Notes

Chemical potentials are ignored since self-energies are constant. This method exists purely for API compatibility with other surface Green’s function calculators.

sigma(E, i, conv=1e-05)

Get self-energy matrix for a specific contact.

Parameters:

• E (float) – Energy point (unused - kept for API compatibility)

• i (int) – Contact index (0 for left, 1 for right)

• conv (float, optional) – Convergence parameter (unused - kept for API compatibility)

Returns:

Self-energy matrix for contact i

Return type:

ndarray

sigmaTot(E, conv=1e-05)

Calculate total self-energy matrix from all contacts.

Parameters:

• E (float) – Energy point (unused - kept for API compatibility)

• conv (float, optional) – Convergence parameter (unused - kept for API compatibility)

Returns:

Total self-energy matrix (sum of all contact self-energies)

Return type:

ndarray

Utilities

Matrix Tools

Matrix manipulation utilities for quantum transport calculations.

This module provides helper functions for handling matrices in quantum transport calculations, with a focus on: - Self-energy matrix construction - Density and Fock matrix manipulation - Spin treatment (restricted, unrestricted, and generalized) - Integration with Gaussian’s matrix formats

The functions handle three types of spin treatments: - ‘r’: Restricted (same orbitals for alpha and beta electrons) - ‘ro’/’u’: Unrestricted (separate alpha and beta orbitals) - ‘g’: Generalized (non-collinear spin treatment)

gauNEGF.matTools.formSigma(inds, V, nsto, S=0)

Form a self-energy matrix for specified orbitals.

Creates a self-energy matrix of size nsto x nsto with values V at the specified orbital indices. Can handle both scalar and matrix-valued self-energies.

Parameters:

• inds (list of int) – Orbital indices where self-energy should be applied

• V (complex or ndarray) – Self-energy value(s) to insert - If scalar: Same value used for all specified orbitals - If matrix: Must match size of indices

• nsto (int) – Total number of orbitals (size of resulting matrix)

• S (ndarray or int, optional) – Overlap matrix for broadening term. If 0, identity used (default: 0)

Returns:

Complex self-energy matrix of size nsto x nsto

Return type:

ndarray

gauNEGF.matTools.getDen(bar, spin)

Build density matrix from Gaussian checkpoint file.

Extracts the density matrix from a Gaussian checkpoint file based on the specified spin treatment. Handles restricted, unrestricted, and generalized cases.

Parameters:

• bar (QCBinAr) – Gaussian interface object

• spin (str) – Spin treatment to use: - ‘r’: Restricted (same orbitals for alpha/beta) - ‘ro’/’u’: Unrestricted (separate alpha/beta) - ‘g’: Generalized (non-collinear)

Returns:

Density matrix in appropriate format for spin treatment: - Restricted: Single block - Unrestricted: Two blocks (alpha/beta) - Generalized: Single block with complex entries

Return type:

ndarray

Raises:

ValueError – If spin treatment is not recognized

gauNEGF.matTools.getEnergies(bar, spin)

Get orbital energies from Gaussian checkpoint file.

Extracts orbital energies from a Gaussian checkpoint file based on the specified spin treatment. Converts energies from Hartrees to eV.

Parameters:

• bar (QCBinAr) – Gaussian interface object

• spin (str) – Spin treatment to use: - ‘r’: Restricted (same orbitals for alpha/beta) - ‘ro’/’u’: Unrestricted (separate alpha/beta) - ‘g’: Generalized (non-collinear)

Returns:

All orbital energies in eV, sorted in ascending order. For restricted/unrestricted spin, alpha and beta levels are interleaved before the final sort, so each energy level appears twice; for generalized spin, each level appears once.

Return type:

ndarray

Raises:

ValueError – If spin treatment is not recognized

gauNEGF.matTools.getFock(bar, spin)

Build Fock matrix from Gaussian checkpoint file.

Extracts the Fock matrix and orbital indices from a Gaussian checkpoint file based on the specified spin treatment. Handles restricted, unrestricted, and generalized cases.

Parameters:

• bar (QCBinAr) – Gaussian interface object

• spin (str) – Spin treatment to use: - ‘r’: Restricted (same orbitals for alpha/beta) - ‘ro’/’u’: Unrestricted (separate alpha/beta) - ‘g’: Generalized (non-collinear)

Returns:

(Fock, locs) where: - Fock: ndarray, Fock matrix in appropriate format for spin treatment - locs: ndarray, Orbital indices with positive for alpha/paired and

negative for beta orbitals

Return type:

tuple

Raises:

ValueError – If spin treatment is not recognized

gauNEGF.matTools.storeDen(bar, P, spin)

Store density matrix in Gaussian checkpoint format.

Converts the density matrix to the appropriate format based on spin treatment and stores it in the Gaussian checkpoint file. Handles compression and proper typing of matrices.

Parameters:

• bar (QCBinAr) – Gaussian interface object

• P (ndarray) – Density matrix to store

• spin (str) – Spin treatment to use: - ‘r’: Restricted (same orbitals for alpha/beta) - ‘ro’/’u’: Unrestricted (separate alpha/beta) - ‘g’: Generalized (non-collinear)

Notes

For restricted calculations, the density is divided by 2 to account for double occupation. For generalized calculations, complex matrices are used.

Raises:

ValueError – If spin treatment is not recognized

Integration Tools

JAX-powered Green’s Functions Integration

Supports both retarded (Gr) and lesser (G<) Green’s functions.

Author: William Livernois

gauNEGF.integrate.GrInt(F, S, g, Elist, weights)

Integrate retarded Green’s function over energy using JAX parallelization.

Parameters:

• F (ndarray) – Fock matrix (NxN)

• S (ndarray) – Overlap matrix (NxN)

• g (surfG object) – Surface Green’s function calculator with sigmaTot(E) method

• Elist (ndarray) – Array of energies in eV (Mx1)

• weights (ndarray) – Array of weights for each energy (Mx1)

Returns:

Integrated retarded Green’s function (NxN)

Return type:

ndarray

gauNEGF.integrate.GrIntCross(F, S, g, Elist, weights)

Integrate G^R with co-accumulation of cross-term scalar.

Returns (lineInt, cross_scalar) where: - lineInt = sum_k w_k * G^R(z_k) (NxN matrix, same as GrInt) - cross_scalar = sum_k w_k * Tr(G^R(z_k) @ Q_tot(z_k)) (complex scalar)

The cross-term delta_N = -(1/pi) * Im(cross_scalar).

For orthogonal systems (crossTermQTot returns None), uses the fast vmap/scan path via GrInt with no cross-term computation. For non-orthogonal systems, uses a single-pass vmap that computes sigmaTot, G^R, and Q_tot once per energy point.

gauNEGF.integrate.GrLessInt(F, S, g, Elist, weights, ind=None)

Integrate lesser Green’s function over energy using JAX parallelization.

Parameters:

• F (ndarray) – Fock matrix (NxN)

• S (ndarray) – Overlap matrix (NxN)

• g (surfG object) – Surface Green’s function calculator

• Elist (ndarray) – Array of energies in eV (Mx1)

• weights (ndarray) – Array of weights for each energy (Mx1)

• ind (int, optional) – Contact index for partial density calculation (default: None)

Returns:

Integrated lesser Green’s function (NxN)

Return type:

ndarray

Spin Tools

Spin manipulation tools for non-collinear NEGF calculations.

This module provides utilities for handling spin degrees of freedom in non-collinear (generalized) spin-polarized transport calculations. It includes:

• Pauli matrix definitions as module-level constants

• SU(2) rotation matrix generation via matrix exponential

• Orthogonal spin direction generation for spin-torque analysis

• Hirshfeld charge and magnetic moment extraction from Gaussian logs

Rotation matrices are constructed using the standard spin-1/2 representation: R(n, omega) = exp(-i * omega/2 * (n . sigma)), where n is the rotation axis unit vector and sigma are the Pauli matrices.

gauNEGF.spinTools.constructSOCterm(lambdas)

Construct spin-orbit coupling term for s, p and d orbitals.

Builds the L·S operator in the real orbital basis by transforming from the spherical harmonic basis, then assembles a block-diagonal 18x18 matrix scaled by the per-shell coupling parameters.

Parameters:

lambdas (list) – List of spin-orbit coupling parameters [lambda_s, lambda_p, lambda_d] for s, p and d orbitals respectively

Returns:

18x18 matrix containing spin-orbit coupling terms in the basis (s_up, s_dn, px_up, px_dn, py_up, py_dn, pz_up, pz_dn,

d3z2r2_up, d3z2r2_dn, dxz_up, dxz_dn, dyz_up, dyz_dn, dx2y2_up, dx2y2_dn, dxy_up, dxy_dn)

Return type:

ndarray

gauNEGF.spinTools.genOrthRotFile(filename)

Use genOrthRots() to generate a set of orthogonal spin rotations based on the direction of maximum magnetic moment in a Gaussian output file.

Parameters:

filename (str) – Path to the Gaussian output (.log) file.

Returns:

• rotations (list of ndarray, each of shape (2, 2)) – 7 SU(2) rotation matrices in order: [identity, +vec1, -vec1, +vec2, -vec2, -original (via vec1), -original (via vec2)].

• directions (ndarray of shape (7, 3)) – Unit vectors corresponding to each rotation target.

• max_moment_index (int) – Index of the atom with the maximum magnetic moment.

gauNEGF.spinTools.genOrthRotGrids(spinVec=None, filename=None, dphi=0.7853981633974483)

Generate a grid of orthogonal spin rotations using genOrthRots() from the spinVec or filename, depending on which is provided. A list of 2 grids are returned, one for each of the orthogonal directions (+vec1, +vec2).

Parameters:

• spinVec (array_like, shape (3,), optional) – Input spin quantization axis as a unit vector.

• filename (str, optional) – If given, the direction is determined from the file’s largest magnetization.

• dphi (float, optional) – Step size in radians. Default is np.pi/4.

Returns:

• rotations_grids (list of list of ndarray) – Each entry is a list of SU(2) rotation matrices (shape [n,2,2]). There are two entries, corresponding to the two orthogonal directions (+vec1, +vec2).

• angles_grids (list of ndarray) – Angles (in radians) used to generate the respective rotation grid.

gauNEGF.spinTools.genOrthRots(spinVec)

Generate a set of orthogonal spin rotations for spin-torque analysis.

Given an input spin quantization axis, constructs 7 rotation matrices mapping to: the original axis (identity), two pairs of orthogonal axes (+/-vec1, +/-vec2), and two independent 180 degree rotations to the anti-parallel directions (-original via vec1, -original via vec2). The orthogonal basis is built to avoid singularities when the input is aligned with any coordinate axis.

Returns:

• rotations (list of ndarray, each of shape (2, 2)) – 7 SU(2) rotation matrices in order: [identity, +vec1, -vec1, +vec2, -vec2, -original (via vec1), -original (via vec2)].

• directions (ndarray of shape (7, 3)) – Unit vectors corresponding to each rotation target.

Notes

Indices 5 and 6 both target -original but rotate around different axes (vec1 and vec2 respectively).

gauNEGF.spinTools.genRot(n, omega)

Generate an SU(2) rotation matrix for spin-1/2.

Constructs R = exp(-i * omega/2 * (n . sigma)) using the module-level Pauli matrices. This represents a rotation by angle omega about axis n in spin space.

Parameters:

• n (array_like of shape (3,)) – Unit vector defining the rotation axis in (x, y, z) coordinates.

• omega (float) – Rotation angle in radians.

Returns:

Complex-valued SU(2) rotation matrix.

Return type:

ndarray of shape (2, 2)

gauNEGF.spinTools.genRotsGrid(spinVec, dphi=0.7853981633974483)

Generate a grid of spin rotations around a given axis. Starts at 0 rads and goes to 2*pi-dphi rads in steps of dphi.

Parameters:

• spinVec (array_like of shape (3,)) – Unit vector defining the rotation axis in (x, y, z) coordinates.

• dphi (float, optional) – Step size in radians. Default is np.pi/4.

Returns:

• rotations (list of ndarray, each of shape (2, 2)) – SU(2) rotation matrices.

• angles (ndarray of shape (n,)) – Angles in radians corresponding to each rotation.

gauNEGF.spinTools.get_hirshfeld_data(filename)

Extract Hirshfeld charges and magnetic moments from a Gaussian log file.

Parses the last occurrence of the Hirshfeld partition block in the log file, returning per-atom net charges and 3D magnetic moment vectors.

Parameters:

filename (str) – Path to the Gaussian output (.log) file.

Returns:

• charges (ndarray of shape (N,)) – Hirshfeld net charges for each atom.

• mag_moments (ndarray of shape (N, 3)) – Hirshfeld magnetic moment vectors (Mx, My, Mz) for each atom.

• None – Returned if no Hirshfeld section is found in the file.

gauNEGF.spinTools.spinorTF(rho_initial, rho_target)

Find unitary transformation U such that U @ rho_initial @ U.conj().T = rho_target using Bloch sphere rotation.

Parameters:

• rho_initial (ndarray) – Initial density matrix

• rho_target (ndarray) – Target density matrix

Returns:

Unitary transformation matrix U

Return type:

ndarray

JIT / Linear Algebra Helpers

Utility functions compiled with JIT for gauNEGF.

Contains commonly used pure mathematical functions that are reused across multiple modules in the gauNEGF package.

gauNEGF.utils.eig(A)

Compute eigenvalues and eigenvectors of a square matrix.

Wrapper around JAX’s general eigenvalue decomposition for non-Hermitian matrices. For Hermitian matrices, prefer eigh() which uses a more stable algorithm.

Parameters:

A (jax array) – Square matrix (NxN), may be complex or non-Hermitian.

Returns:

• eigenvalues (jax array of shape (N,)) – Eigenvalues of A (may be complex).

• eigenvectors (jax array of shape (NxN)) – Eigenvectors as columns, with A @ eigenvectors = eigenvectors @ diag(eigenvalues).

gauNEGF.utils.eigh(A)

Compute eigenvalues and eigenvectors of a Hermitian matrix.

Wrapper around JAX’s Hermitian eigenvalue decomposition. Assumes the input matrix is Hermitian (A = A^H) and uses a more stable algorithm than eig().

Parameters:

A (jax array) – Hermitian matrix (NxN).

Returns:

• eigenvalues (jax array of shape (N,)) – Real-valued eigenvalues of A in ascending order.

• eigenvectors (jax array of shape (NxN)) – Eigenvectors as columns, with A @ eigenvectors = eigenvectors @ diag(eigenvalues).

gauNEGF.utils.fractional_matrix_power(S, power)

Calculate matrix power S^p using eigendecomposition. Supports fractional powers including negative values like -0.5.

Parameters:

• S (jax array) – Input matrix (should be Hermitian for numerical stability)

• power (float) – Power to raise matrix to (e.g., 0.5 for sqrt, -0.5 for inverse sqrt)

Returns:

Matrix power S^p

Return type:

jax array

Notes

This function is optimized for Hermitian matrices (like overlap matrices) and uses eigendecomposition: S^p = V @ D^p @ V^H where S = V @ D @ V^H.

Unlike JAX’s matrix_power, this function properly handles fractional powers including negative values.

gauNEGF.utils.inv(A)

Compute matrix inverse using JAX linalg solve.

Solves the linear system A @ X = I for X, which is equivalent to computing the inverse A^(-1). This method is more numerically stable than direct inversion for ill-conditioned matrices.

Parameters:

A (jax array) – Square invertible matrix (NxN).

Returns:

Inverse matrix A^(-1) (NxN).

Return type:

jax array

gauNEGF.utils.inv_sqrt_general(M)

Inverse square root M^(-1/2) for a general diagonalizable matrix.

Diagonalizes M = V @ diag(D) @ V^(-1) with the general (non-symmetric) eigendecomposition and returns Y = V @ diag(D^(-1/2)) @ V^(-1). This satisfies Y @ M @ Y = I to machine precision for ANY diagonalizable M – Hermitian, complex-symmetric, or neither – since D^(-1/2) D D^(-1/2) = 1 holds for every branch of the scalar square root.

Use this for the effective overlap S_eff = S - X_asymp. The retarded contact self-energy (Sigma = A @ g_surf @ A^dagger, g_surf complex-symmetric) carries broadening Gamma = i(Sigma - Sigma^dagger) and is therefore NON-Hermitian, so X_asymp and S_eff are non-Hermitian. eigh assumes a Hermitian matrix and uses only one triangle, which silently computes the wrong S_eff^(-1/2); use this routine instead. The physical lower-contour density is independent of the sqrt branch chosen here, because G(E) = Y @ (E*I - Fbar)^(-1) @ Y reduces to [E*S_eff - H_eff]^(-1) for any Y with Y @ Y = S_eff^(-1).

Parameters:

M (jax array (N, N)) – General diagonalizable matrix (may be complex / non-Hermitian).

Returns:

Y = M^(-1/2), satisfying Y @ M @ Y = I.

Return type:

jax array (N, N), complex

Notes

Uses JAX’s general eig (jnp.linalg.eig, CPU backend) – the same wrapper used elsewhere in this module. Accuracy degrades if the eigenvector matrix is ill-conditioned (near-defective M); a Schur-based matrix power is the robust fallback for that case.

Configuration Reference

See also

Configuration and Tuning – When to change ETA, SCF_DAMPING, fermiMethod, and TEMPERATURE

Global configuration settings for gauNEGF calculations.

This module provides centralized default parameters used throughout the gauNEGF package.

gauNEGF.config.shard_array(array, axis=0)

Shard an array across all devices along the specified axis. Use this for input arrays before passing to jitted functions.

Parameters:

• array – JAX array to shard

• axis – Which axis to shard (default: 0, typically batch dimension)

Returns:

Sharded array distributed across devices

Developer / Extensibility Reference

Protocol definitions for the surfG interface.

SurfGProtocol: the single interface that density.py and integrate.py interact with.

Implemented by: surfG, surfG3, surfGB, surfGTest, surfGAt3D, surfGBAt

These are structural (duck-typed) protocols – classes do NOT need to inherit from them. They exist for documentation and optional type-checking with isinstance() via @runtime_checkable.

Limitation: @runtime_checkable only verifies that methods and attributes exist on the instance. It does not check signatures, return types, or attribute types. For full correctness, see test_cross_term.py and test_jit_fermi_shift.py which test functional behavior.

class gauNEGF.protocols.SurfGProtocol(*args, **kwargs)

Bases: Protocol

Protocol for surface Green’s function calculators.

This is the interface that density.py and integrate.py interact with.

Implemented by: surfG, surfG3, surfGB, surfGTest, surfGAt3D, surfGBAt

Note: Parameter names for setF() vary across implementations (mu1/mu2 vs muL/muR). This does not affect structural typing.

F: ndarray

S: ndarray

crossTermQ(E: complex, i: int, conv: float = Ellipsis) → ndarray | None

Symmetrized cross-term overlap matrix Q_sym_i in full device basis.

Returns None if contact i has orthogonal coupling (S_DL = 0).

The cross-term electron count correction is:

delta_N_i = -(1/pi) * Im(sum_k w_k * Tr(G_DD^R(z_k) @ Q_sym_i(z_k)))

crossTermQTot(E: complex, conv: float = Ellipsis) → ndarray | None

Sum of Q_sym over all contacts, or None if all orthogonal.

eta: float

num_contacts: int

setF(F: ndarray, mu1: float, mu2: float) → None

Update Fock matrix and contact chemical potentials.

sigma(E: complex, i: int, conv: float = Ellipsis) → ndarray

Self-energy for contact i in full device basis (N x N).

sigmaTot(E: complex, conv: float = Ellipsis) → ndarray

Total self-energy from all contacts (N x N).