Skip to content

API documentation

current_solver(f_hamiltonian, f_leads, system_parameters, lead_parameters, kbT)

Constructs a solver function for quantum transport problems. The returned function computes the current for a given bias, using the provided Hamiltonian and lead callables, system and lead parameters, and temperature.

Parameters:

Name Type Description Default
f_hamiltonian callable

Callable that returns the Hamiltonian matrix for given system parameters.

 required
f_leads list[callable]

List of callables, each returning a lead matrix for given lead parameters.

 required
system_parameters dict

Dictionary of parameters for the system Hamiltonian.

 required
lead_parameters dict

Dictionary of parameters for the leads.

 required
kbT float

Temperature in energy units.

 required

Returns:

Name Type Description
callable

A function f(biases) that computes the current for the specified biases.
Source code in qrate/solvers.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
def current_solver(f_hamiltonian, f_leads, system_parameters, lead_parameters, kbT):
    """
    Constructs a solver function for quantum transport problems. The returned function computes the current for a given bias, using the provided Hamiltonian and lead callables, system and lead parameters, and temperature.

    Args:
        f_hamiltonian (callable): Callable that returns the Hamiltonian matrix for given system parameters.
        f_leads (list[callable]): List of callables, each returning a lead matrix for given lead parameters.
        system_parameters (dict): Dictionary of parameters for the system Hamiltonian.
        lead_parameters (dict): Dictionary of parameters for the leads.
        kbT (float): Temperature in energy units.

    Returns:
        callable: A function f(biases) that computes the current for the specified biases.
    """
    # Compute the Hamiltonian matrix and lead matrices
    hamiltonian_matrix = f_hamiltonian(**system_parameters)
    lead_matrices = []

    for f_lead in f_leads:
        sig = inspect.signature(f_lead)
        valid_params = sig.parameters
        filtered_kwargs = {k: v for k, v in lead_parameters.items() if k in valid_params}
        result = f_lead(**filtered_kwargs)
        lead_matrices.append(result)

    # Setup a callable for the W matrices
    f_w_matrices = setup_system(hamiltonian_matrix, lead_matrices, kbT)

    # Since the same eigenstates are reused for all biases
    # we can create a single function of biases
    def f(biases):
        w_plus, w_minus = f_w_matrices(biases)
        current = get_current(w_plus, w_minus)
        return current

    return f

get_current(w_plus, w_minus)

Solves the rate equations and computes the current for the system.

Parameters:

Name Type Description Default
w_plus ndarray

Transition rates for electron addition.

 required
w_minus ndarray

Transition rates for electron removal.

 required

Returns:

Type Description

np.ndarray: Current for each lead and bias configuration.
Source code in qrate/solvers.py 
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def get_current(w_plus, w_minus):
    """
    Solves the rate equations and computes the current for the system.

    Args:
        w_plus (np.ndarray): Transition rates for electron addition.
        w_minus (np.ndarray): Transition rates for electron removal.

    Returns:
        np.ndarray: Current for each lead and bias configuration.
    """
    A = w_plus + w_minus
    n_l, _, n_e, _ = A.shape

    U = np.sum(A, axis=-2)
    for i in range(A.shape[-1]):
        A[:, :, i, i] -= U[:, :, i]

    p_lhs = np.concatenate(
        (A, np.ones(shape=(n_l, n_l, 1, n_e))), axis=2
    )  # shape: n_l, n_l, n_e+1, n_e
    p_lhs = p_lhs.reshape(n_l * n_l, n_e + 1, n_e)  # shape: n_l*n_l,n_e+1,n_e
    p_lhs_inv = list(
        map(np.linalg.pinv, [p_lhs[i, :, :] for i in range(p_lhs.shape[0])])
    )  # Here we can also use multiprocessing
    p_rhs = np.hstack(
        (np.zeros(shape=(n_l * n_l, n_e)), np.ones(shape=(n_l * n_l, 1)))
    ).flatten()  # shape: n_l*n_l*(n_e+1)

    p_vec = scipy.sparse.block_diag(p_lhs_inv) @ p_rhs
    p_vec = p_vec.reshape(n_l, n_l, n_e)  # shape: n_l, n_l, n_e

    current = np.einsum("abcd,ebd -> aec", (w_plus - w_minus), p_vec)
    current = np.sum(current, axis=-1)

    return current

n_fermi(eigs, mu, kbT)

Computes the Fermi-Dirac distribution for a set of eigenvalues and chemical potentials.

Parameters:

Name Type Description Default
eigs ndarray

Array of eigenvalues.

 required
mu ndarray

Array of chemical potentials (biases).

 required
kbT float

Temperature in energy units.

 required

Returns:

Type Description

np.ndarray: Fermi-Dirac occupation probabilities with shape (n_l, n_l, n_e, n_e).
Source code in qrate/solvers.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
def n_fermi(eigs: np.ndarray, mu: np.ndarray, kbT: float):
    """
    Computes the Fermi-Dirac distribution for a set of eigenvalues and chemical potentials.

    Args:
        eigs (np.ndarray): Array of eigenvalues.
        mu (np.ndarray): Array of chemical potentials (biases).
        kbT (float): Temperature in energy units.

    Returns:
        np.ndarray: Fermi-Dirac occupation probabilities with shape (n_l, n_l, n_e, n_e).
    """
    delta_eigs = eigs[:, None] - eigs[None, :]
    delta_mu = mu[None, :] - mu[:, None]
    arg = (delta_eigs[None, None, ...] - delta_mu[:, :, None, None]) / kbT
    return 1 / (np.exp(arg) + 1)

setup_system(hamiltonian_matrix, lead_matrices, kbT)

Diagonalizes the Hamiltonian and prepares the system for transport calculations.

Parameters:

Name Type Description Default
hamiltonian_matrix ndarray

The system Hamiltonian matrix.

 required
lead_matrices list[ndarray]

List of lead coupling matrices.

 required
kbT float

Temperature in energy units.

 required

Returns:

Name Type Description
callable

Function that computes W matrices for given biases.
Source code in qrate/solvers.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def setup_system(hamiltonian_matrix, lead_matrices, kbT):
    """
    Diagonalizes the Hamiltonian and prepares the system for transport calculations.

    Args:
        hamiltonian_matrix (np.ndarray): The system Hamiltonian matrix.
        lead_matrices (list[np.ndarray]): List of lead coupling matrices.
        kbT (float): Temperature in energy units.

    Returns:
        callable: Function that computes W matrices for given biases.
    """
    # Diagonalize the Hamiltonian matrix
    eigvals, eigvecs = np.linalg.eigh(hamiltonian_matrix)

    # Compute the transition matrices
    t_minus, t_plus = t_matrices(eigvecs, np.array(lead_matrices))

    def f(biases):
        w_plus = np.abs(t_plus[np.newaxis, ...]) ** 2 * n_fermi(eigvals, biases, kbT)
        # shape: n_l,n_l,n_e,n_e Hewre is where sum over leads needs to be taken
        w_minus = np.abs(t_minus[np.newaxis, ...]) ** 2 * (
            1 - n_fermi(-eigvals, biases, kbT)
        )
        return w_plus, w_minus

    return f

t_matrices(vecs, lead_coupling)

Computes the transition matrices t_minus and t_plus for the system.

Parameters:

Name Type Description Default
vecs ndarray

Eigenvectors of the Hamiltonian.

 required
lead_coupling ndarray

Lead coupling matrices.

 required

Returns:

Name Type Description
tuple

(t_minus, t_plus) transition matrices.
Source code in qrate/solvers.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def t_matrices(vecs: np.ndarray, lead_coupling: callable):
    """
    Computes the transition matrices t_minus and t_plus for the system.

    Args:
        vecs (np.ndarray): Eigenvectors of the Hamiltonian.
        lead_coupling (np.ndarray): Lead coupling matrices.

    Returns:
        tuple: (t_minus, t_plus) transition matrices.
    """
    t_minus = np.array(
        [
            vecs.conj().T @ lead_coupling[i, :, :] @ vecs
            for i in range(lead_coupling.shape[0])
        ]
    )
    t_plus = np.transpose(t_minus.conj(), axes=(0, 2, 1))  # n_l,n_e,n_e
    return t_minus, t_plus