API documentation

compute_rates(eigvals, t_minus, t_plus, biases, temperature)

Compute lead-resolved transition rates for a bias sweep.

Parameters

eigvals : np.ndarray Eigenvalues E_n of shape (N,) aligned with the eigenvectors in t_*. t_minus : np.ndarray Lead tunneling matrices (L, N, N) for removal processes. t_plus : np.ndarray Lead tunneling matrices (L, N, N) for addition processes. biases : np.ndarray Bias grid of shape (S, L) (broadcast from (S,) or (L,)) where S is sweep points and L is number of leads. temperature : float Positive temperature (same units as eigvals and biases).

Returns

tuple[np.ndarray, np.ndarray] (w_plus, w_minus) of shape (S, L, N, N) with add/remove rates per sweep point and lead.

Source code in qrate/solvers.py

def compute_rates(
    eigvals: np.ndarray,
    t_minus: np.ndarray,
    t_plus: np.ndarray,
    biases: np.ndarray,
    temperature: float,
) -> Tuple[np.ndarray, np.ndarray]:
    """Compute lead-resolved transition rates for a bias sweep.

    Parameters
    ----------
    eigvals : np.ndarray
        Eigenvalues ``E_n`` of shape ``(N,)`` aligned with the eigenvectors in ``t_*``.
    t_minus : np.ndarray
        Lead tunneling matrices ``(L, N, N)`` for removal processes.
    t_plus : np.ndarray
        Lead tunneling matrices ``(L, N, N)`` for addition processes.
    biases : np.ndarray
        Bias grid of shape ``(S, L)`` (broadcast from ``(S,)`` or ``(L,)``) where
        ``S`` is sweep points and ``L`` is number of leads.
    temperature : float
        Positive temperature (same units as ``eigvals`` and ``biases``).

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        ``(w_plus, w_minus)`` of shape ``(S, L, N, N)`` with add/remove rates per
        sweep point and lead.
    """
    if len(t_minus) == 0:
        return np.array([]), np.array([])

    n_leads_from_t = t_minus.shape[0]
    biases = np.asarray(biases)

    # Handle 1D biases
    if biases.ndim == 1:
        if n_leads_from_t == 1:
            # Ambiguity: Is it 1 lead with S sweep points? Or L leads with 1 sweep point?
            # Since L=1, it must be S sweep points.
            biases = biases[:, None]
        elif biases.shape[0] == n_leads_from_t:
            # Assume L leads, 1 sweep point
            biases = biases[None, :]
        else:
            # Ambiguous. E.g. L=2, biases len=5. Assume 5 sweep points, broadcast to 2 leads?
            # Let's default to S sweeps, 1 bias (broadcast later) if length matches S?
            # Safer to treat as S sweeps for 1 virtual lead -> broadcast to L leads?
            # We conform to (S, 1)
            biases = biases[:, None]

    biases = np.atleast_2d(biases)
    n_sweeps, n_leads_bias = biases.shape

    # Delta E matrix: delta_E[n', n] = E_n' - E_n
    # Note: Using row=final, col=initial convention for Rate Matrix W.
    delta_E = eigvals[:, None] - eigvals[None, :]

    # biases: (S, L) -> (S, L, 1, 1) to broadcast against (N, N)
    mu = biases[:, :, None, None]

    # dE: (1, 1, N, N)
    dE = delta_E[None, None, :, :]

    # Fermi functions
    # f_plus = f(E_final - E_initial - mu) = f(dE - mu)
    f_plus = n_fermi(dE - mu, temperature)

    # f_minus_term = 1 - f(E_initial - E_final - mu) = 1 - f(-dE - mu)
    # This corresponds to hole occupation or empty state prob?
    arg_minus = -dE - mu
    f_minus_c = 1.0 - n_fermi(arg_minus, temperature)

    # Prepare T squared
    # t_plus: (L, N, N). Indices (lead, final, initial)
    abs_t_plus_sq = np.abs(t_plus) ** 2  # (L, N, N)
    abs_t_minus_sq = np.abs(t_minus) ** 2  # (L, N, N)

    # Broadcast to (S, L, N, N)
    # abs_t_plus_sq[None, ...] broadcasts to (1, L, N, N) which matches f_plus (S, L, N, N)
    w_plus = abs_t_plus_sq[None, :, :, :] * f_plus
    w_minus = abs_t_minus_sq[None, :, :, :] * f_minus_c

    return w_plus, w_minus

compute_t_matrices(eigvecs, lead_couplings)

Transform lead couplings into the eigenbasis.

Parameters

eigvecs : np.ndarray Unitary eigenvector matrix U of shape (N, N) (columns = eigenstates). lead_couplings : list[np.ndarray] Coupling matrices C_alpha of shape (N, N) for each lead alpha. If empty, returns empty arrays.

Returns

tuple[np.ndarray, np.ndarray] t_minus and t_plus with shape (L, N, N) where t_minus[alpha] = U^dagger C_alpha U and t_plus[alpha] = t_minus[alpha]^dagger.

Source code in qrate/solvers.py

def compute_t_matrices(
    eigvecs: np.ndarray, lead_couplings: List[np.ndarray]
) -> Tuple[np.ndarray, np.ndarray]:
    """Transform lead couplings into the eigenbasis.

    Parameters
    ----------
    eigvecs : np.ndarray
        Unitary eigenvector matrix ``U`` of shape ``(N, N)`` (columns = eigenstates).
    lead_couplings : list[np.ndarray]
        Coupling matrices ``C_alpha`` of shape ``(N, N)`` for each lead ``alpha``.
        If empty, returns empty arrays.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        ``t_minus`` and ``t_plus`` with shape ``(L, N, N)`` where
        ``t_minus[alpha] = U^dagger C_alpha U`` and
        ``t_plus[alpha] = t_minus[alpha]^dagger``.
    """
    U = eigvecs
    U_dagger = U.conj().T

    n_leads = len(lead_couplings)
    # Check if we have leads, if not handle gracefullly?
    if n_leads == 0:
        return np.array([]), np.array([])

    n_states = U.shape[1]

    t_minus = np.zeros((n_leads, n_states, n_states), dtype=complex)

    for i, C in enumerate(lead_couplings):
        # T^- = U^dagger @ C @ U
        t_minus[i] = U_dagger @ C @ U

    # T^+ = (T^-)^dagger = conjugate transpose of each matrix
    t_plus = np.transpose(t_minus.conj(), axes=(0, 2, 1))

    return t_minus, t_plus

g_matrix(solver, biases, step=1e-05)

Compute conductance via central finite differences of a current solver.

For each sweep point s and lead index j we evaluate

.. math:: G_{ij}(s) pprox rac{I_i(V_s + h e_j) - I_i(V_s - h e_j)}{2h}

using the provided solver to recompute currents at perturbed biases.

Parameters

solver : callable Function mapping biases (S, L) to currents (S, L). biases : np.ndarray Bias grid of shape (S, L) (or (S,)). step : float, optional Perturbation size h for the central difference.

Returns

np.ndarray Conductance tensor (S, L, L) where G[s, i, j] = dI_i/dV_j.

Source code in qrate/solvers.py

def g_matrix(
    solver: Callable[[np.ndarray], np.ndarray],
    biases: np.ndarray,
    step: float = 1e-5,
) -> np.ndarray:
    """Compute conductance via central finite differences of a current solver.

    For each sweep point ``s`` and lead index ``j`` we evaluate

    .. math::
       G_{ij}(s) \approx \frac{I_i(V_s + h e_j) - I_i(V_s - h e_j)}{2h}

    using the provided ``solver`` to recompute currents at perturbed biases.

    Parameters
    ----------
    solver : callable
        Function mapping biases ``(S, L)`` to currents ``(S, L)``.
    biases : np.ndarray
        Bias grid of shape ``(S, L)`` (or ``(S,)``).
    step : float, optional
        Perturbation size ``h`` for the central difference.

    Returns
    -------
    np.ndarray
        Conductance tensor ``(S, L, L)`` where ``G[s, i, j] = dI_i/dV_j``.
    """

    biases = np.asarray(biases)
    if biases.ndim == 1:
        biases = biases[:, None]

    if biases.ndim != 2:
        raise ValueError("biases must be of shape (S, L) or (S,)")

    S, L = biases.shape

    # Stack all + and - perturbations to call the solver only twice.
    deltas = np.zeros((L, 1, L))
    deltas[np.arange(L), 0, np.arange(L)] = step

    biases_plus = (biases[None, :, :] + deltas).reshape(L * S, L)
    biases_minus = (biases[None, :, :] - deltas).reshape(L * S, L)

    I_plus = solver(biases_plus).reshape(L, S, L)
    I_minus = solver(biases_minus).reshape(L, S, L)

    diff = I_plus - I_minus  # (L, S, L)
    G = np.transpose(diff, (1, 2, 0)) / (2 * step)  # (S, L, L)

    return G

n_fermi(energy_diff, temperature)

Fermi-Dirac occupations clipped to avoid overflow.

Parameters

energy_diff : np.ndarray Energy offsets E_f - E_i - mu; any broadcastable shape. temperature : float Positive temperature (same energy units). Must be > 0.

Returns

np.ndarray Occupations in [0, 1] with energy_diff.shape.

Source code in qrate/solvers.py

def n_fermi(energy_diff: np.ndarray, temperature: float) -> np.ndarray:
    """Fermi-Dirac occupations clipped to avoid overflow.

    Parameters
    ----------
    energy_diff : np.ndarray
        Energy offsets ``E_f - E_i - mu``; any broadcastable shape.
    temperature : float
        Positive temperature (same energy units). Must be ``> 0``.

    Returns
    -------
    np.ndarray
        Occupations in ``[0, 1]`` with ``energy_diff.shape``.
    """
    arg = energy_diff / temperature
    arg = np.clip(arg, -100.0, 100.0)
    return 1.0 / (np.exp(arg) + 1.0)

solve_master_equation(w_plus, w_minus)

Solve the steady-state master equation for each bias point.

Parameters

w_plus : np.ndarray Addition rates (S, L, N, N) from compute_rates. w_minus : np.ndarray Removal rates (S, L, N, N) from compute_rates.

Returns

tuple[np.ndarray, np.ndarray, np.ndarray] occupations (S, N state probabilities), currents (S, L net current into the system per lead), and W (S, N, N total transition matrix summed over leads).

Source code in qrate/solvers.py

def solve_master_equation(
    w_plus: np.ndarray, w_minus: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Solve the steady-state master equation for each bias point.

    Parameters
    ----------
    w_plus : np.ndarray
        Addition rates ``(S, L, N, N)`` from ``compute_rates``.
    w_minus : np.ndarray
        Removal rates ``(S, L, N, N)`` from ``compute_rates``.

    Returns
    -------
    tuple[np.ndarray, np.ndarray, np.ndarray]
        ``occupations`` (``S, N`` state probabilities), ``currents`` (``S, L`` net
        current into the system per lead), and ``W`` (``S, N, N`` total transition
        matrix summed over leads).
    """
    if w_plus.size == 0:
        return np.array([]), np.array([]), np.array([])

    n_sweeps, _, n_states, _ = w_plus.shape

    # Total rates per lead: Gamma_{alpha, n->m} (final, initial)
    Gamma = w_plus + w_minus  # (S, L, N, N)

    # Sum over leads to get total transition rates between states
    W = np.sum(Gamma, axis=1)  # (S, N, N)

    # Rate matrix A for linear system A p = 0
    # A_{nm} = W_{nm} (for n != m)
    # A_{nn} = - sum_{k} W_{kn} (escape rate)

    # Column sums (sum over rows) give total escape rates from state n
    # sum_k W_{kn} is sum of column n (fixed initial state n, sum over final states k)
    escape_rates = np.sum(W, axis=1)  # (S, N)

    A = W.copy()
    diag_idx = np.arange(n_states)
    A[:, diag_idx, diag_idx] -= escape_rates[:, diag_idx]

    # Enforce normalization sum p_n = 1 by appending a row of ones
    # A = np.concatenate((A, np.ones((n_sweeps, 1, n_states))), axis=1)
    # rhs = np.zeros((n_sweeps, n_states + 1))
    # rhs[:, -1] = 1.0

    A[:, -1, :] = 1.0  # Replace last row with ones
    rhs = np.zeros((n_sweeps, n_states))
    rhs[:, -1] = 1.0

    occupations = np.empty((n_sweeps, n_states), dtype=A.dtype)
    for s in range(n_sweeps):
        occupations[s] = np.linalg.solve(A[s], rhs[s])

    # Calculate currents.
    # Current I_alpha (into system) = sum_nm (Gamma^+_{alpha, nm} p_m - Gamma^-_{alpha, nm} p_m)
    # Sum over n(final), m(initial)

    # w_plus: (S, L, N, N)
    # p: (S, N)

    current_in = np.einsum("sanm, sm -> sa", w_plus, occupations)
    current_out = np.einsum("sanm, sm -> sa", w_minus, occupations)

    currents = current_in - current_out

    return occupations, currents, W

solve_system(hamiltonian, leads, temperature, parity_op=None, *, system_parameters=None, lead_parameters=None)

Build a callable that returns conductance for given biases.

Parameters

hamiltonian : array-like or callable Static Hamiltonian matrix (N, N) or callable returning one; called with system_parameters filtered to matching keywords. leads : list[array-like or callable] Lead coupling matrices C_alpha (N, N) or callables returning them; each is called with lead_parameters filtered to matching keywords. temperature : float Positive temperature used in Fermi factors (same energy units as biases). parity_op : np.ndarray, optional Parity operator commuting with H; required (no automatic fallback). system_parameters : dict, optional Keyword arguments forwarded to the Hamiltonian callable. lead_parameters : dict, optional Keyword arguments forwarded to each lead callable.

Returns

callable Function solver(biases) -> G that maps a bias grid (S, L) to the conductance tensor (S, L, L) using g_matrix.

Source code in qrate/solvers.py

def solve_system(
    hamiltonian: Any,
    leads: List[Any],
    temperature: float,
    parity_op: Optional[np.ndarray] = None,
    *,
    system_parameters: Optional[Dict[str, Any]] = None,
    lead_parameters: Optional[Dict[str, Any]] = None,
) -> Callable[[np.ndarray], np.ndarray]:
    """Build a callable that returns conductance for given biases.

    Parameters
    ----------
    hamiltonian : array-like or callable
        Static Hamiltonian matrix ``(N, N)`` or callable returning one; called
        with ``system_parameters`` filtered to matching keywords.
    leads : list[array-like or callable]
        Lead coupling matrices ``C_alpha`` ``(N, N)`` or callables returning them;
        each is called with ``lead_parameters`` filtered to matching keywords.
    temperature : float
        Positive temperature used in Fermi factors (same energy units as biases).
    parity_op : np.ndarray, optional
        Parity operator commuting with ``H``; required (no automatic fallback).
    system_parameters : dict, optional
        Keyword arguments forwarded to the Hamiltonian callable.
    lead_parameters : dict, optional
        Keyword arguments forwarded to each lead callable.

    Returns
    -------
    callable
        Function ``solver(biases) -> G`` that maps a bias grid ``(S, L)`` to the
        conductance tensor ``(S, L, L)`` using ``g_matrix``.
    """

    if parity_op is None:
        raise ValueError("parity_op is required; no default eigh fallback is used.")

    H = _call_if_needed(hamiltonian, system_parameters, "hamiltonian")
    lead_mats = [
        _call_if_needed(ld, lead_parameters, f"lead[{i}]") for i, ld in enumerate(leads)
    ]

    evals, _, evecs = parity_diag(H, parity_op)
    t_minus, t_plus = compute_t_matrices(evecs, lead_mats)

    def current_solver(biases: np.ndarray) -> np.ndarray:
        w_plus, w_minus = compute_rates(evals, t_minus, t_plus, biases, temperature)
        _, currents, _ = solve_master_equation(w_plus, w_minus)
        return currents

    def solver(biases: np.ndarray) -> np.ndarray:
        return g_matrix(current_solver, biases)

    return solver