Guide 5: Workflow Recipes

Once you have a converged SCF calculation (see 1D Chain Contacts in gauNEGF or Bethe Lattice Contacts), these recipes show how to compute IV curves, study temperature effects, and handle long energy scans efficiently.

Overview

This guide covers four advanced workflows for transport calculations:

  • IV Curve Sweep - Loop over bias voltages to measure current-voltage characteristics

  • Multi-Temperature Sweep - Compare transport at different temperatures (e.g., 0K vs 300K)

  • Transport Checkpointing - Save intermediate results during long energy scans and resume gracefully

  • Warm-Start SCF - Hot-start a new SCF from a prior converged calculation

Each recipe includes working code adapted from the gauNEGF example suite.

Recipe 1: IV Curve Sweep

Measuring current as a function of applied bias voltage (IV curves) is one of the most common transport calculations. This recipe shows how to loop over a voltage list, perform SCF at each voltage, and save results to MAT files for easy resume.

The key insight: if you interrupt a run, you can resume from saved checkpoint files without recalculating previous voltages.

Code Example

# adapted from AuStudies/AuBetheFerrocene.py
from gauNEGF.scfE import NEGFE
from gauNEGF.transport import calculate_current, SigmaCalculator
from scipy import io
import os
import numpy as np

fn = 'AuBetheFerrocene'
extra = 'ESCF_BetheSOC_contFermi'
spin = 'g'
Vhi = 0.5
dV = 0.1
conv = 1e-3

# Initialize NEGF object with Bethe contact
negf = NEGFE(fn=fn, func='b3lyp', basis='chkbasis', spin=spin,
             fullSCF=False, route='integral=grid=superfine')
negf.setContactBethe([[1, 2, 3, 4, 5, 6],
                      [36, 37, 38, 39, 40, 41]], 'AuSOC')

# Build voltage list: up, then down, then back to zero
Vlist = np.concatenate((
    np.arange(dV, Vhi, dV),
    np.arange(Vhi, -1*Vhi, -1*dV),
    np.arange(-1*Vhi, 0, dV)
))

# Create output directory and log file
os.makedirs(f"{fn}_{extra}", exist_ok=True)
f = open(f'{fn}_{extra}/IV_{Vhi}_{dV}.log', 'w')
f.write("V (volts),I (Amps),Iuu,Iud,Idu,Idd\n")
f.write("0.0,0.0,0.0,0.0,0.0,0.0\n")
f.flush()

# Loop over voltage points
IList = np.zeros_like(Vlist)
IsList = np.zeros((len(Vlist), 4))

for i, V in enumerate(Vlist):
    matName = f'{fn}_{extra}/num{i}_{V:.2f}V.mat'
    # Resume from checkpoint if it exists
    if os.path.exists(matName):
        A = io.loadmat(matName)
        negf.setDen(A['den'])
        negf.setVoltage(V, A['fermi'][0][0])
    else:
        negf.setVoltage(V, negf.fermi)

    # SCF at this voltage
    negf.SCF(conv, 0.01, 1000, checkpoint=False)
    # Save density and Fock for resume
    negf.saveMAT(matName)

    # Compute current (returns total I and spin-resolved Is)
    har_to_eV = 27.211386
    I, Is = calculate_current(
        negf.F * har_to_eV, negf.S,
        SigmaCalculator(negf.g),
        negf.fermi, V, spin=spin
    )
    print(f"V = {V:.2f} V, I = {I:.3E} A")
    f.write(f"{V:.2f},{I:.3E},{Is[0]:.3E},{Is[1]:.3E},{Is[2]:.3E},{Is[3]:.3E}\n")
    f.flush()

    IList[i] = I
    IsList[i, :] = Is

# Save results to MAT file
io.savemat(f'{fn}_{extra}/IV_{Vhi}_{dV}.mat',
           {'Vlist': Vlist, 'IList': IList, 'IsList': IsList})
f.close()

Key Points

  • Voltage list construction: Use np.arange to build the up, down, and back portions. This explores hysteresis if present.

  • Resume logic: Before SCF, check if the MAT file exists using os.path.exists(). If so, load the prior density and Fermi level to hot-start.

  • Current returns: calculate_current() returns a tuple (I, Is):

    • I is the total current (scalar).

    • Is is spin-resolved current with shape (4,) for spin=’g’: [Iuu, Iud, Idu, Idd]. The indices represent the four possible spin configurations in non-collinear spin.

  • Logging: Write results to a plain-text log file for real-time monitoring.

Cross-references:

Recipe 2: Multi-Temperature Sweep

Temperature affects electron transport through the Fermi-Dirac distribution in the contacts. This recipe shows how to compute transmission (or current) at two different temperatures and compare.

Critical insight: You must run a full SCF at each temperature, not just change a parameter. The contact self-energy is temperature-dependent, so the density matrix must re-equilibrate.

Code Example

# adapted from NEGFTests/CNanowire.py
from gauNEGF.scfE import NEGFE
from gauNEGF.transport import cohTransE
from scipy import io
import numpy as np

fn = 'C2'
har_to_eV = 27.211386
Elist = np.linspace(-5, 5, 1000)

# First temperature: 0K (default)
print("SCF at T=0K...")
negf = NEGFE(fn=fn, func='b3lyp', basis='lanl2dz', fullSCF=False)
negf.setContact1D([[1], [2]], eta=1e-4)  # Default T=0
negf.setVoltage(0.0, fermiMethod='poly')
negf.SCF(1e-2, 0.02, 100)
negf.SCF(1e-4, 0.02, 1000, pulay=False)
negf.saveMAT('CNanowire_ESCF.mat')

T_0K = cohTransE(Elist + negf.fermi, negf.F * har_to_eV, negf.S, negf.g)
io.savemat('CNanowire_TESCF_0K.mat',
           {'Elist': Elist, 'fermi': negf.fermi, 'T': T_0K})

# Second temperature: 300K
print("SCF at T=300K...")
negf = NEGFE(fn=fn, func='b3lyp', basis='lanl2dz', fullSCF=False)
negf.setContact1D([[1], [2]], eta=1e-4, T=300)  # Set T=300K
negf.setVoltage(0.0, fermiMethod='poly')
negf.SCF(1e-2, 0.02, 100)
negf.SCF(1e-4, 0.02, 1000, pulay=False)
negf.saveMAT('CNanowire_ESCF_300K.mat')

T_300K = cohTransE(Elist + negf.fermi, negf.F * har_to_eV, negf.S, negf.g)
io.savemat('CNanowire_TESCF_300K.mat',
           {'Elist': Elist, 'fermi': negf.fermi, 'T': T_300K})

Key Points

  • Two independent SCF runs: Each temperature gets its own NEGFE object and full SCF. You cannot reuse one SCF across multiple temperatures.

  • Temperature in setContact1D(): Pass T=<value> (in Kelvin) to gauNEGF.scfE.NEGFE.setContact1D(). The default is 0K (T=0). This sets the Fermi-Dirac distribution in the contact self-energy.

  • Fermi level: Each temperature has its own Fermi level (stored in negf.fermi). This is computed during SCF and reflects the electron population at that temperature.

  • Transmission comparison: Plot T_0K vs T_300K to see how thermal broadening affects transport. Thermal effects typically broaden the transmission near the Fermi level.

Cross-references:

Recipe 3: Transport Checkpointing (Long Energy Scans)

Long transmission or DOS scans (e.g., 5000+ energy points) can take hours or days. This recipe shows how to use checkpointing to save intermediate results and resume after an interruption without losing progress.

Code Example

# adapted from tests/test_transport_checkpointing.py
from gauNEGF.transport import (
    SigmaCalculator, calculate_transmission, calculate_dos
)
import numpy as np
import tempfile
import os

# Assume you have F, S (Fock and overlap matrices)
# and sigma_calc (SigmaCalculator object)

energy_list = np.linspace(-5, 5, 5000)  # Long scan
checkpoint_file = 'transmission_checkpoint.npz'
checkpoint_interval = 100  # Save every 100 energies

# First call: compute transmission with checkpointing
print("Computing transmission with checkpointing...")
transmission = calculate_transmission(
    F, S, sigma_calc, energy_list,
    spin='r',
    checkpoint_file=checkpoint_file,
    checkpoint_interval=checkpoint_interval
)

# If interrupted, just call again with the same arguments.
# The function will:
# 1. Load the checkpoint file
# 2. Identify which energies are already computed (-1 indicates not done)
# 3. Resume from the next missing energy
# 4. Continue saving at the same interval

print(f"Transmission computed: {len(transmission)} energies")

# For DOS with checkpointing:
dos_total, dos_per_site = calculate_dos(
    F, S, sigma_calc, energy_list,
    spin='r',
    checkpoint_file='dos_checkpoint.npz',
    checkpoint_interval=100
)

Key Points

  • checkpoint_file: Path to an HDF5-like (.npz) file. The file is created on first call and reused on resume.

  • checkpoint_interval: How many energy points to compute before saving. Smaller values (e.g., 10) save frequently but add I/O overhead. Larger values (e.g., 500) are faster but lose more progress if interrupted. A value of 50-100 is typical.

  • Resume behavior: If you kill the script and restart with identical arguments, calculate_transmission() and calculate_dos() will:

    1. Load the saved checkpoint.

    2. Detect uncalculated energies (marked -1).

    3. Compute only those missing energies.

    4. Return the complete array.

  • No manual checkpoint management: You do not need to manage file I/O yourself. The function handles loading, checking for completeness, and saving automatically.

Cross-references:

Recipe 4: Warm-Start SCF from Prior Run

When running similar geometries or exploring slight parameter variations, starting SCF from a prior converged density matrix (warm-start) can halve convergence time compared to starting from scratch.

Code Example

# adapted from AuStudies/AuBetheFerrocene.py
from gauNEGF.scfE import NEGFE
from scipy import io

# Load a prior converged result
A = io.loadmat('AuBetheFerrocene_ESCF_BetheSOC_contFermi.mat')
prior_density = A['den']
prior_fermi = A['fermi'][0][0]

# Create a new NEGF object (same or similar geometry)
negf = NEGFE(fn='AuBetheFerrocene', func='b3lyp', basis='chkbasis',
             spin='g', fullSCF=False)
negf.setContactBethe([[1, 2, 3, 4, 5, 6],
                      [36, 37, 38, 39, 40, 41]], 'AuSOC')

# Hot-start with prior density and Fermi
negf.setDen(prior_density)  # Load density matrix
negf.setVoltage(0.0, prior_fermi)  # Set initial Fermi level
# Add `negf.setVoltage(0.0)` here if you want to re-enable Fermi search
# during the warm-started SCF (see contacts_bethe.rst warm-start pattern).

# SCF from warm-start (converges faster)
negf.SCF(1e-3, 0.02, 200, checkpoint=False)

Key Points

  • setDen(): Load a prior density matrix using gauNEGF.scf.NEGF.setDen(). This initializes the density in the SCF loop.

  • setVoltage() second argument: The second positional argument to gauNEGF.scf.NEGF.setVoltage() is the Fermi level (or chemical potential). Pass a scalar float (not a tuple).

  • MAT file keys: Standard keys in a gauNEGF MAT file are:

    • 'den' - Density matrix (ndarray)

    • 'fermi' - Fermi level as a 1x1 array (access with [0][0])

    • 'F' - Fock matrix

    • 'S' - Overlap matrix

  • checkpoint=False: Ensure that the density matrix will not be overwritten by existing f'{fn}_P.mat' file. Conversely, copy your warm-start file to that filename to skip all of the previous steps!

  • Speed improvement: Warm-start typically converges in 10-30% of the iterations needed for a cold start, especially if the geometry is similar.

Cross-references:

Integration Methods Reference

The recipes above use the real-axis integration method for transport (via calculate_current() and cohTransE()). For a comparison with complex contour integration and when to use each, see examples/IntegralDemo.ipynb – run it in Jupyter on a compute node.

The demo shows:

  • Accuracy vs. speed trade-offs

  • How to set energy windows via setIntegralLimits()

  • Convergence criteria and numerical precision

Summary

These four recipes cover the most common advanced workflows in gauNEGF:

  1. IV sweeps - Fast current-voltage measurements with resume capability.

  2. Temperature studies - Independent SCF runs at different T to study thermal effects.

  3. Checkpointing - Long scans with automatic save/resume for robustness.

  4. Warm-start - Fast re-convergence when starting from a prior result.

Combine these recipes as needed. For example:

  • IV sweep at T=300K using checkpointing for transmission calculations.

  • Multi-temperature comparison of warm-started SCF runs.

  • Checkpoint very long DOS scans using energy-dependent sigma from contacts.