Skip to content

Open Systems Module

Overview

In this section, we outline the main functionality of the open_systems module.

First, we will provide some conceptual explanations that provide the user with a necessary background to understand the code. Then we provide some illustrative examples that demonstrate how the code can be used. Finally, we provide the source code as an API reference to the source code.

Concepts

In this section, we can add some theoretical background/explanation of relevant concepts.

Examples

In this section, we can add some illustrative examples.

Source Code

qflux.open_systems

DynamicsOS

Class for open-system dynamics (Lindblad equation).

This class provides methods to simulate open-system dynamics described by the Lindblad equation.

Attributes:

Name Type Description
Nsys int

System Hilbert Space Dimension.
Hsys ndarray

Hamiltonian of the system (shape (N, N)).
rho0 ndarray

Initial density matrix (shape (N, N)).
c_ops List[ndarray]

List of collapse operators (each of shape (N, N)).
Source code in src/qflux/open_systems/numerical_methods.py 
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 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
 72
 73
 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
110
111
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
class DynamicsOS:
    """Class for open-system dynamics (Lindblad equation).

    This class provides methods to simulate open-system dynamics described by the Lindblad equation.

    Attributes:
        Nsys (int): System Hilbert Space Dimension.
        Hsys (np.ndarray): Hamiltonian of the system (shape (N, N)).
        rho0 (np.ndarray): Initial density matrix (shape (N, N)).
        c_ops (List[np.ndarray]): List of collapse operators (each of shape (N, N)).
    """

    def __init__(
        self,
        Nsys: int,
        Hsys: np.ndarray,
        rho0: np.ndarray,
        c_ops: Optional[List[np.ndarray]] = None
    ) -> None:
        """
        Initialize the DynamicsOS instance.

        Args:
            Nsys (int): System Hilbert Space Dimension.
            Hsys (np.ndarray): Hamiltonian of the system.
            rho0 (np.ndarray): Initial density matrix.
            c_ops (Optional[List[np.ndarray]]): List of collapse operators. Defaults to an empty list.
        """
        if c_ops is None:
            c_ops = []
        self.Nsys: int = Nsys
        self.Hsys: np.ndarray = Hsys
        self.rho0: np.ndarray = rho0
        self.c_ops: List[np.ndarray] = c_ops

    def Gt_matrix_expo(self, time_arr: List[float], Is_show_step: bool = False) -> List[np.ndarray]:
        """
        Compute the propagator of the Lindblad equation using the matrix exponential.

        The propagator is computed by exponentiating the Liouvillian operator defined by the
        system Hamiltonian and collapse operators.

        Args:
            time_arr (List[float]): Array of time values for the simulation.
            Is_show_step (bool, optional): If True, prints the current simulation step. Defaults to False.

        Returns:
            List[np.ndarray]: List of propagators corresponding to each time in `time_arr`.
        """
        ident_h: np.ndarray = np.eye(self.Nsys, dtype=np.complex128)

        # Build the A matrix for time-derivation of the vectorized density matrix.
        Amat: np.ndarray = -1j * (np.kron(self.Hsys, ident_h) - np.kron(ident_h, self.Hsys.T))
        for i in range(len(self.c_ops)):
            Amat += 0.5 * (
                2.0 * np.kron(self.c_ops[i], self.c_ops[i].conj())
                - np.kron(ident_h, (self.c_ops[i].T @ self.c_ops[i].conj()))
                - np.kron((self.c_ops[i].T.conj() @ self.c_ops[i]), ident_h)
            )

        G_prop: List[np.ndarray] = []
        for i, t in enumerate(time_arr):
            if Is_show_step:
                print("step", i, "time", t)
            Gt: np.ndarray = LA.expm(Amat * t)
            G_prop.append(Gt)
        return G_prop

    def propagate_matrix_exp(
        self,
        time_arr: List[float],
        observable: np.ndarray,
        Is_store_state: bool = False,
        Is_show_step: bool = False,
        Is_Gt: bool = False,
    ) -> Any:
        """
        Solve the Lindblad equation using matrix exponential.

        This method computes the propagator, evolves the initial density matrix, and calculates
        the expectation value of the observable over time. Optionally, it stores the evolved density matrices.

        Args:
            time_arr (List[float]): Time array for dynamic simulation.
            observable (np.ndarray): Observable matrix for which the expectation value is computed.
            Is_store_state (bool, optional): If True, stores the density matrices at each time step.
                                             Defaults to False.
            Is_show_step (bool, optional): If True, prints the current simulation step. Defaults to False.
            Is_Gt (bool, optional): If True, includes the propagators in the result. Defaults to False.

        Returns:
            Result: An object with the following attributes:
                - expect (List[float]): List of expectation values over time.
                - density_matrix (List[np.ndarray], optional): List of density matrices (if `Is_store_state` is True).
                - Gprop (List[np.ndarray], optional): List of propagators (if `Is_Gt` is True).
        """

        class Result:
            """Class for storing propagation results."""
            def __init__(self, store_state: bool, include_Gt: bool) -> None:
                self.expect: List[float] = []
                if store_state:
                    self.density_matrix: List[np.ndarray] = []
                if include_Gt:
                    self.Gprop: Optional[List[np.ndarray]] = None

        result = Result(Is_store_state, Is_Gt)

        # Compute the propagator of the Lindblad equation.
        G_prop: List[np.ndarray] = self.Gt_matrix_expo(time_arr, Is_show_step)
        if Is_Gt:
            result.Gprop = G_prop

        # Initialize the vectorized density matrix.
        vec_rho0: np.ndarray = self.rho0.reshape(self.Nsys**2)

        for i, _ in enumerate(time_arr):
            vec_rhot: np.ndarray = G_prop[i] @ vec_rho0
            # Reshape back to density matrix form.
            rhot: np.ndarray = vec_rhot.reshape(self.Nsys, self.Nsys)

            if Is_store_state:
                result.density_matrix.append(rhot)
            result.expect.append(np.trace(rhot @ observable).real)

        return result

    def propagate_qt(self, time_arr: List[float], observable: Any, **kwargs: Any) -> List[float]:
        """
        Propagate the system using QuTiP's `mesolve` function.

        This method solves the Lindblad master equation using QuTiP's `mesolve` to compute the expectation
        values of the observable over time.

        Args:
            time_arr (List[float]): Time array for dynamic simulation.
            observable (Any): Observable operator(s) for which the expectation value is computed.
                              Can be a single operator or a list of operators.
            **kwargs: Additional keyword arguments to pass to `mesolve`.

        Returns:
            List[float]: List of expectation values of the observable over time.
        """
        c_ops: List[Qobj] = [Qobj(c_op) for c_op in self.c_ops]

        if isinstance(observable, list):
            obs = [Qobj(op) for op in observable]
        else:
            obs = Qobj(observable)

        result = mesolve(Qobj(self.Hsys), Qobj(self.rho0), time_arr, c_ops, obs, **kwargs)
        return result.expect

Gt_matrix_expo(time_arr, Is_show_step=False)

Compute the propagator of the Lindblad equation using the matrix exponential.

The propagator is computed by exponentiating the Liouvillian operator defined by the system Hamiltonian and collapse operators.

Parameters:

Name Type Description Default
time_arr List[float]

Array of time values for the simulation.

 required
Is_show_step bool

If True, prints the current simulation step. Defaults to False.

 False

Returns:

Type Description
List[ndarray]

List[np.ndarray]: List of propagators corresponding to each time in time_arr.
Source code in src/qflux/open_systems/numerical_methods.py 
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
72
73
74
75
76
77
def Gt_matrix_expo(self, time_arr: List[float], Is_show_step: bool = False) -> List[np.ndarray]:
    """
    Compute the propagator of the Lindblad equation using the matrix exponential.

    The propagator is computed by exponentiating the Liouvillian operator defined by the
    system Hamiltonian and collapse operators.

    Args:
        time_arr (List[float]): Array of time values for the simulation.
        Is_show_step (bool, optional): If True, prints the current simulation step. Defaults to False.

    Returns:
        List[np.ndarray]: List of propagators corresponding to each time in `time_arr`.
    """
    ident_h: np.ndarray = np.eye(self.Nsys, dtype=np.complex128)

    # Build the A matrix for time-derivation of the vectorized density matrix.
    Amat: np.ndarray = -1j * (np.kron(self.Hsys, ident_h) - np.kron(ident_h, self.Hsys.T))
    for i in range(len(self.c_ops)):
        Amat += 0.5 * (
            2.0 * np.kron(self.c_ops[i], self.c_ops[i].conj())
            - np.kron(ident_h, (self.c_ops[i].T @ self.c_ops[i].conj()))
            - np.kron((self.c_ops[i].T.conj() @ self.c_ops[i]), ident_h)
        )

    G_prop: List[np.ndarray] = []
    for i, t in enumerate(time_arr):
        if Is_show_step:
            print("step", i, "time", t)
        Gt: np.ndarray = LA.expm(Amat * t)
        G_prop.append(Gt)
    return G_prop

__init__(Nsys, Hsys, rho0, c_ops=None)

Initialize the DynamicsOS instance.

Parameters:

Name Type Description Default
Nsys int

System Hilbert Space Dimension.

 required
Hsys ndarray

Hamiltonian of the system.

 required
rho0 ndarray

Initial density matrix.

 required
c_ops Optional[List[ndarray]]

List of collapse operators. Defaults to an empty list.

 None
Source code in src/qflux/open_systems/numerical_methods.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def __init__(
    self,
    Nsys: int,
    Hsys: np.ndarray,
    rho0: np.ndarray,
    c_ops: Optional[List[np.ndarray]] = None
) -> None:
    """
    Initialize the DynamicsOS instance.

    Args:
        Nsys (int): System Hilbert Space Dimension.
        Hsys (np.ndarray): Hamiltonian of the system.
        rho0 (np.ndarray): Initial density matrix.
        c_ops (Optional[List[np.ndarray]]): List of collapse operators. Defaults to an empty list.
    """
    if c_ops is None:
        c_ops = []
    self.Nsys: int = Nsys
    self.Hsys: np.ndarray = Hsys
    self.rho0: np.ndarray = rho0
    self.c_ops: List[np.ndarray] = c_ops

propagate_matrix_exp(time_arr, observable, Is_store_state=False, Is_show_step=False, Is_Gt=False)

Solve the Lindblad equation using matrix exponential.

This method computes the propagator, evolves the initial density matrix, and calculates the expectation value of the observable over time. Optionally, it stores the evolved density matrices.

Parameters:

Name Type Description Default
time_arr List[float]

Time array for dynamic simulation.

 required
observable ndarray

Observable matrix for which the expectation value is computed.

 required
Is_store_state bool

If True, stores the density matrices at each time step. Defaults to False.

 False
Is_show_step bool

If True, prints the current simulation step. Defaults to False.

 False
Is_Gt bool

If True, includes the propagators in the result. Defaults to False.

 False

Returns:

Name Type Description
Result Any

An object with the following attributes: - expect (List[float]): List of expectation values over time. - density_matrix (List[np.ndarray], optional): List of density matrices (if Is_store_state is True). - Gprop (List[np.ndarray], optional): List of propagators (if Is_Gt is True).
Source code in src/qflux/open_systems/numerical_methods.py 
 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
110
111
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
def propagate_matrix_exp(
    self,
    time_arr: List[float],
    observable: np.ndarray,
    Is_store_state: bool = False,
    Is_show_step: bool = False,
    Is_Gt: bool = False,
) -> Any:
    """
    Solve the Lindblad equation using matrix exponential.

    This method computes the propagator, evolves the initial density matrix, and calculates
    the expectation value of the observable over time. Optionally, it stores the evolved density matrices.

    Args:
        time_arr (List[float]): Time array for dynamic simulation.
        observable (np.ndarray): Observable matrix for which the expectation value is computed.
        Is_store_state (bool, optional): If True, stores the density matrices at each time step.
                                         Defaults to False.
        Is_show_step (bool, optional): If True, prints the current simulation step. Defaults to False.
        Is_Gt (bool, optional): If True, includes the propagators in the result. Defaults to False.

    Returns:
        Result: An object with the following attributes:
            - expect (List[float]): List of expectation values over time.
            - density_matrix (List[np.ndarray], optional): List of density matrices (if `Is_store_state` is True).
            - Gprop (List[np.ndarray], optional): List of propagators (if `Is_Gt` is True).
    """

    class Result:
        """Class for storing propagation results."""
        def __init__(self, store_state: bool, include_Gt: bool) -> None:
            self.expect: List[float] = []
            if store_state:
                self.density_matrix: List[np.ndarray] = []
            if include_Gt:
                self.Gprop: Optional[List[np.ndarray]] = None

    result = Result(Is_store_state, Is_Gt)

    # Compute the propagator of the Lindblad equation.
    G_prop: List[np.ndarray] = self.Gt_matrix_expo(time_arr, Is_show_step)
    if Is_Gt:
        result.Gprop = G_prop

    # Initialize the vectorized density matrix.
    vec_rho0: np.ndarray = self.rho0.reshape(self.Nsys**2)

    for i, _ in enumerate(time_arr):
        vec_rhot: np.ndarray = G_prop[i] @ vec_rho0
        # Reshape back to density matrix form.
        rhot: np.ndarray = vec_rhot.reshape(self.Nsys, self.Nsys)

        if Is_store_state:
            result.density_matrix.append(rhot)
        result.expect.append(np.trace(rhot @ observable).real)

    return result

propagate_qt(time_arr, observable, **kwargs)

Propagate the system using QuTiP's mesolve function.

This method solves the Lindblad master equation using QuTiP's mesolve to compute the expectation values of the observable over time.

Parameters:

Name Type Description Default
time_arr List[float]

Time array for dynamic simulation.

 required
observable Any

Observable operator(s) for which the expectation value is computed. Can be a single operator or a list of operators.

 required
**kwargs Any

Additional keyword arguments to pass to mesolve.

 {}

Returns:

Type Description
List[float]

List[float]: List of expectation values of the observable over time.
Source code in src/qflux/open_systems/numerical_methods.py 
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def propagate_qt(self, time_arr: List[float], observable: Any, **kwargs: Any) -> List[float]:
    """
    Propagate the system using QuTiP's `mesolve` function.

    This method solves the Lindblad master equation using QuTiP's `mesolve` to compute the expectation
    values of the observable over time.

    Args:
        time_arr (List[float]): Time array for dynamic simulation.
        observable (Any): Observable operator(s) for which the expectation value is computed.
                          Can be a single operator or a list of operators.
        **kwargs: Additional keyword arguments to pass to `mesolve`.

    Returns:
        List[float]: List of expectation values of the observable over time.
    """
    c_ops: List[Qobj] = [Qobj(c_op) for c_op in self.c_ops]

    if isinstance(observable, list):
        obs = [Qobj(op) for op in observable]
    else:
        obs = Qobj(observable)

    result = mesolve(Qobj(self.Hsys), Qobj(self.rho0), time_arr, c_ops, obs, **kwargs)
    return result.expect

QubitDynamicsOS

Bases: DynamicsOS

Class for simulating quantum dynamics using either vectorized density matrix or Kraus operator representations.

This class provides methods to initialize state vectors, construct quantum circuits, and perform quantum simulations using Qiskit backends.

Source code in src/qflux/open_systems/quantum_simulation.py 
 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
110
111
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
class QubitDynamicsOS(DynamicsOS):
    """
    Class for simulating quantum dynamics using either vectorized density matrix or Kraus operator representations.

    This class provides methods to initialize state vectors, construct quantum circuits,
    and perform quantum simulations using Qiskit backends.
    """

    def __init__(self, rep: str = 'Density', **kwargs: Any) -> None:
        """
        Initialize a QubitDynamicsOS instance.

        Depending on the representation, either "Density" or "Kraus", the number of qubits is computed.
        Additional keyword arguments are passed to the base DynamicsOS class.

        Args:
            rep (str, optional): Representation type, either 'Density' or 'Kraus'. Defaults to 'Density'.
            **kwargs: Additional keyword arguments for the DynamicsOS initializer.
        """
        super().__init__(**kwargs)

        if rep == 'Density':
            # Vectorized density matrix representation
            self.rep: str = 'Density'
            self.Nqb: int = int(np.log2(self.Nsys**2))
        elif rep == 'Kraus':
            # Kraus operator representation
            self.rep = 'Kraus'
            self.Nqb = int(np.log2(self.Nsys))

        # The counting qubits bit string and observable matrix are initialized to None.
        self.count_str: Optional[List[str]] = None
        self.observable: Optional[np.ndarray] = None

        # Default dilation method for quantum simulation.
        self.dilation_method: str = 'Sz-Nagy'

    def set_dilation_method(self, method: str) -> None:
        """
        Set the dilation method for quantum simulation.

        Args:
            method (str): The dilation method, e.g., 'Sz-Nagy', 'SVD', or 'SVD-Walsh'.
        """
        self.dilation_method = method

    def set_count_str(self, count_str: List[str]) -> None:
        """
        Set the counting bit string for measurement.

        Args:
            count_str (List[str]): The counting bit string.
        """
        self.count_str = count_str

    def set_observable(self, observable: np.ndarray) -> None:
        """
        Set the observable for the quantum simulation.

        Args:
            observable (np.ndarray): The observable matrix.
        """
        self.observable = observable

    def init_statevec_vecdens(self) -> Tuple[np.ndarray, float]:
        """
        Initialize the state vector from the initial density operator using vectorized representation.

        The initial density matrix is reshaped into a vector and normalized.

        Returns:
            Tuple[np.ndarray, float]: A tuple containing the normalized state vector and the norm
            of the original vectorized density matrix.
        """
        vec_rho0 = self.rho0.reshape(self.Nsys**2)
        norm0 = LA.norm(vec_rho0, 2)
        statevec = vec_rho0 / norm0
        return statevec, norm0

    def init_statevec_Kraus(self, tol: float = 1e-6) -> Tuple[List[np.ndarray], List[float]]:
        """
        Initialize state vectors from the initial density operator using the Kraus operator representation.

        The density matrix is decomposed using eigenvalue decomposition, and eigenstates with eigenvalues
        below the specified tolerance are ignored.

        Args:
            tol (float, optional): Tolerance for ignoring eigenstates with small eigenvalues. Defaults to 1e-6.

        Returns:
            Tuple[List[np.ndarray], List[float]]: A tuple containing a list of state vectors and a list of
            corresponding probabilities.
        """
        eigenvalues, eigenvectors = LA.eigh(self.rho0)
        statevec: List[np.ndarray] = []
        prob: List[float] = []
        for i in range(len(eigenvalues) - 1, -1, -1):
            if abs(eigenvalues[i]) < tol:
                break
            prob.append(eigenvalues[i])
            statevec.append(eigenvectors[:, i])
        return statevec, prob

    def _get_qiskit_observable(self, Isdilate: bool = False, tol: float = 5e-3) -> SparsePauliOp:
        """
        Prepare and return the Qiskit observable operator.

        Converts the observable matrix to its Pauli representation and returns a SparsePauliOp.

        Args:
            Isdilate (bool, optional): Flag indicating whether to use the dilated observable.
                Defaults to False.
            tol (float, optional): Tolerance for the Pauli decomposition. Defaults to 5e-3.

        Returns:
            SparsePauliOp: The Qiskit representation of the observable.
        """
        if self.observable is None:
            print('Error: observable is None')

        if Isdilate:
            num_qubits = self.Nqb + 1
            Obs_mat = np.zeros((2 * self.Nsys, 2 * self.Nsys), dtype=np.complex128)
            Obs_mat[:self.Nsys, :self.Nsys] = self.observable[:self.Nsys, :self.Nsys]
        else:
            num_qubits = self.Nqb
            Obs_mat = self.observable

        Obs_paulis_dic = tb.ham_to_pauli(Obs_mat, num_qubits, tol=tol)

        # Prepare the Qiskit observable from the Pauli strings of the observable matrix.
        data: List[str] = []
        coef: List[float] = []
        for key in Obs_paulis_dic:
            data.append(key)
            coef.append(Obs_paulis_dic[key])
        obs_q = SparsePauliOp(data, coef)
        return obs_q

    def qc_simulation_kraus(
        self,
        time_arr: List[float],
        shots: int = 1024,
        Kraus: Optional[Dict[int, List[np.ndarray]]] = None,
        Gprop: Optional[List[np.ndarray]] = None,
        tolk: float = 1e-5,
        tolo: float = 1e-5,
        **kwargs: Any
    ) -> np.ndarray:
        """
        Perform quantum simulation using the Kraus operator representation.

        This method simulates the quantum system dynamics over a series of time steps using a Kraus operator-based approach.
        It constructs quantum circuits for each Kraus operator and initial state, runs the simulation using Qiskit's Estimator,
        and accumulates the measurement results.

        Args:
            time_arr (List[float]): List of time steps for simulation.
            shots (int, optional): Number of shots for each measurement. Defaults to 1024.
            Kraus (Optional[Dict[int, List[np.ndarray]]], optional): Dictionary mapping time step index to a list of Kraus operators.
                If None, Kraus operators are generated from the propagator. Defaults to None.
            Gprop (Optional[List[np.ndarray]], optional): Propagator matrix (or list of matrices) for simulation.
                If None, it will be calculated. Defaults to None.
            tolk (float, optional): Tolerance for generating Kraus operators. Defaults to 1e-5.
            tolo (float, optional): Tolerance for observable decomposition. Defaults to 1e-5.
            **kwargs: Additional keyword arguments for propagator calculation.

        Returns:
            np.ndarray: Array containing the quantum simulation results.
        """
        nsteps = len(time_arr)

        # Generate Kraus operators if not provided.
        if Kraus is None:
            Kraus = {}
            if Gprop is None:
                print('Calculating the propagator')
                Gprop = self.Gt_matrix_expo(time_arr, **kwargs)
            print('Generating the Kraus operators')
            for i in range(nsteps):
                print('At step', i, 'of', nsteps)
                Kraus[i] = gen_Kraus_list(Gprop[i], self.Nsys, tol=tolk)
        print('Kraus operator generation complete')

        # Perform Qiskit simulation using the Estimator.
        estimator = Estimator()

        statevec, prob = self.init_statevec_Kraus()
        n_inistate = len(statevec)
        print('Number of initial states in the density matrix:', n_inistate)
        print('Probabilities:', prob)

        # Obtain the Qiskit observable.
        obs_q = self._get_qiskit_observable(Isdilate=True, tol=tolo)

        print('Starting quantum simulation')
        result_simulation = np.zeros(nsteps, dtype=np.float64)

        for i in range(nsteps):
            print('Simulation step', i, 'of', nsteps)
            current_kraus_list = Kraus[i]
            print('Number of Kraus operators:', len(current_kraus_list))
            for kraus_op in current_kraus_list:
                for istate in range(n_inistate):
                    qc = self._create_circuit(kraus_op, statevec[istate], Isscale=False)
                    result = estimator.run(qc, obs_q, shots=shots).result()
                    result_simulation[i] += result.values[0] * prob[istate]

        return result_simulation

    def qc_simulation_vecdens(
        self,
        time_arr: List[float],
        shots: int = 1024,
        backend: Any = AerSimulator(),
        Gprop: Optional[List[np.ndarray]] = None,
        **kwargs: Any
    ) -> np.ndarray:
        """
        Perform quantum simulation using the vectorized density matrix representation.

        This method simulates the quantum system dynamics over a series of time steps by constructing circuits
        based on the vectorized density matrix representation, performing measurements, and processing the results.

        Args:
            time_arr (List[float]): List of time steps for simulation.
            shots (int, optional): Number of measurement shots. Defaults to 1024.
            backend (Any, optional): Quantum simulation backend. Defaults to AerSimulator().
            Gprop (Optional[List[np.ndarray]], optional): Propagator matrix (or list of matrices) for simulation.
                If None, it will be calculated. Defaults to None.
            **kwargs: Additional keyword arguments for propagator calculation.

        Returns:
            np.ndarray: Array containing the quantum simulation results.
        """
        if Gprop is None:
            Gprop = self.Gt_matrix_expo(time_arr, **kwargs)

        nsteps = len(time_arr)

        if self.count_str is None:
            print("Error: count_str is not assigned")

        n_bitstr = len(self.count_str)
        statevec, norm0 = self.init_statevec_vecdens()
        result = np.zeros((nsteps, n_bitstr), dtype=np.float64)

        for i in range(nsteps):
            if i % 100 == 0:
                print('Quantum simulation step', i)
            Gt = Gprop[i]
            circuit, norm = self._create_circuit(Gt, statevec, Isscale=True)
            circuit.measure(range(self.Nqb + 1), range(self.Nqb + 1))
            if self.dilation_method == 'SVD-Walsh':
                circuit = transpile(circuit, backend)
            counts = backend.run(circuit, shots=shots).result().get_counts()
            for j in range(n_bitstr):
                bitstr = self.count_str[j]
                if bitstr in counts:
                    result[i, j] = np.sqrt(counts[bitstr] / shots) * norm * norm0
                else:
                    print('At time', i, 'with shots =', shots, "no counts for", bitstr)

        return result

    def _create_circuit(
        self,
        array: np.ndarray,
        statevec: Union[np.ndarray, List[np.ndarray]],
        Isscale: bool = True
    ) -> QuantumCircuit:
        """
        Construct and return the quantum circuit.

        This method wraps the call to the dilation circuit construction function.

        Args:
            array (np.ndarray): Array used for circuit construction (e.g., propagator or Kraus operator).
            statevec (Union[np.ndarray, List[np.ndarray]]): State vector(s) to be used in the circuit.
            Isscale (bool, optional): Flag indicating whether scaling should be applied. Defaults to True.

        Returns:
            QuantumCircuit: The constructed quantum circuit.
        """
        return dc.construct_circuit(self.Nqb, array, statevec, method=self.dilation_method, Isscale=Isscale)

__init__(rep='Density', **kwargs)

Initialize a QubitDynamicsOS instance.

Depending on the representation, either "Density" or "Kraus", the number of qubits is computed. Additional keyword arguments are passed to the base DynamicsOS class.

Parameters:

Name Type Description Default
rep str

Representation type, either 'Density' or 'Kraus'. Defaults to 'Density'.

 'Density'
**kwargs Any

Additional keyword arguments for the DynamicsOS initializer.

 {}
Source code in src/qflux/open_systems/quantum_simulation.py 
 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
110
def __init__(self, rep: str = 'Density', **kwargs: Any) -> None:
    """
    Initialize a QubitDynamicsOS instance.

    Depending on the representation, either "Density" or "Kraus", the number of qubits is computed.
    Additional keyword arguments are passed to the base DynamicsOS class.

    Args:
        rep (str, optional): Representation type, either 'Density' or 'Kraus'. Defaults to 'Density'.
        **kwargs: Additional keyword arguments for the DynamicsOS initializer.
    """
    super().__init__(**kwargs)

    if rep == 'Density':
        # Vectorized density matrix representation
        self.rep: str = 'Density'
        self.Nqb: int = int(np.log2(self.Nsys**2))
    elif rep == 'Kraus':
        # Kraus operator representation
        self.rep = 'Kraus'
        self.Nqb = int(np.log2(self.Nsys))

    # The counting qubits bit string and observable matrix are initialized to None.
    self.count_str: Optional[List[str]] = None
    self.observable: Optional[np.ndarray] = None

    # Default dilation method for quantum simulation.
    self.dilation_method: str = 'Sz-Nagy'

_create_circuit(array, statevec, Isscale=True)

Construct and return the quantum circuit.

This method wraps the call to the dilation circuit construction function.

Parameters:

Name Type Description Default
array ndarray

Array used for circuit construction (e.g., propagator or Kraus operator).

 required
statevec Union[ndarray, List[ndarray]]

State vector(s) to be used in the circuit.

 required
Isscale bool

Flag indicating whether scaling should be applied. Defaults to True.

 True

Returns:

Name Type Description
QuantumCircuit QuantumCircuit

The constructed quantum circuit.
Source code in src/qflux/open_systems/quantum_simulation.py 
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
def _create_circuit(
    self,
    array: np.ndarray,
    statevec: Union[np.ndarray, List[np.ndarray]],
    Isscale: bool = True
) -> QuantumCircuit:
    """
    Construct and return the quantum circuit.

    This method wraps the call to the dilation circuit construction function.

    Args:
        array (np.ndarray): Array used for circuit construction (e.g., propagator or Kraus operator).
        statevec (Union[np.ndarray, List[np.ndarray]]): State vector(s) to be used in the circuit.
        Isscale (bool, optional): Flag indicating whether scaling should be applied. Defaults to True.

    Returns:
        QuantumCircuit: The constructed quantum circuit.
    """
    return dc.construct_circuit(self.Nqb, array, statevec, method=self.dilation_method, Isscale=Isscale)

_get_qiskit_observable(Isdilate=False, tol=0.005)

Prepare and return the Qiskit observable operator.

Converts the observable matrix to its Pauli representation and returns a SparsePauliOp.

Parameters:

Name Type Description Default
Isdilate bool

Flag indicating whether to use the dilated observable. Defaults to False.

 False
tol float

Tolerance for the Pauli decomposition. Defaults to 5e-3.

 0.005

Returns:

Name Type Description
SparsePauliOp SparsePauliOp

The Qiskit representation of the observable.
Source code in src/qflux/open_systems/quantum_simulation.py 
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def _get_qiskit_observable(self, Isdilate: bool = False, tol: float = 5e-3) -> SparsePauliOp:
    """
    Prepare and return the Qiskit observable operator.

    Converts the observable matrix to its Pauli representation and returns a SparsePauliOp.

    Args:
        Isdilate (bool, optional): Flag indicating whether to use the dilated observable.
            Defaults to False.
        tol (float, optional): Tolerance for the Pauli decomposition. Defaults to 5e-3.

    Returns:
        SparsePauliOp: The Qiskit representation of the observable.
    """
    if self.observable is None:
        print('Error: observable is None')

    if Isdilate:
        num_qubits = self.Nqb + 1
        Obs_mat = np.zeros((2 * self.Nsys, 2 * self.Nsys), dtype=np.complex128)
        Obs_mat[:self.Nsys, :self.Nsys] = self.observable[:self.Nsys, :self.Nsys]
    else:
        num_qubits = self.Nqb
        Obs_mat = self.observable

    Obs_paulis_dic = tb.ham_to_pauli(Obs_mat, num_qubits, tol=tol)

    # Prepare the Qiskit observable from the Pauli strings of the observable matrix.
    data: List[str] = []
    coef: List[float] = []
    for key in Obs_paulis_dic:
        data.append(key)
        coef.append(Obs_paulis_dic[key])
    obs_q = SparsePauliOp(data, coef)
    return obs_q

init_statevec_Kraus(tol=1e-06)

Initialize state vectors from the initial density operator using the Kraus operator representation.

The density matrix is decomposed using eigenvalue decomposition, and eigenstates with eigenvalues below the specified tolerance are ignored.

Parameters:

Name Type Description Default
tol float

Tolerance for ignoring eigenstates with small eigenvalues. Defaults to 1e-6.

 1e-06

Returns:

Type Description
List[ndarray]

Tuple[List[np.ndarray], List[float]]: A tuple containing a list of state vectors and a list of
List[float]

corresponding probabilities.
Source code in src/qflux/open_systems/quantum_simulation.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def init_statevec_Kraus(self, tol: float = 1e-6) -> Tuple[List[np.ndarray], List[float]]:
    """
    Initialize state vectors from the initial density operator using the Kraus operator representation.

    The density matrix is decomposed using eigenvalue decomposition, and eigenstates with eigenvalues
    below the specified tolerance are ignored.

    Args:
        tol (float, optional): Tolerance for ignoring eigenstates with small eigenvalues. Defaults to 1e-6.

    Returns:
        Tuple[List[np.ndarray], List[float]]: A tuple containing a list of state vectors and a list of
        corresponding probabilities.
    """
    eigenvalues, eigenvectors = LA.eigh(self.rho0)
    statevec: List[np.ndarray] = []
    prob: List[float] = []
    for i in range(len(eigenvalues) - 1, -1, -1):
        if abs(eigenvalues[i]) < tol:
            break
        prob.append(eigenvalues[i])
        statevec.append(eigenvectors[:, i])
    return statevec, prob

init_statevec_vecdens()

Initialize the state vector from the initial density operator using vectorized representation.

The initial density matrix is reshaped into a vector and normalized.

Returns:

Type Description
ndarray

Tuple[np.ndarray, float]: A tuple containing the normalized state vector and the norm
float

of the original vectorized density matrix.
Source code in src/qflux/open_systems/quantum_simulation.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def init_statevec_vecdens(self) -> Tuple[np.ndarray, float]:
    """
    Initialize the state vector from the initial density operator using vectorized representation.

    The initial density matrix is reshaped into a vector and normalized.

    Returns:
        Tuple[np.ndarray, float]: A tuple containing the normalized state vector and the norm
        of the original vectorized density matrix.
    """
    vec_rho0 = self.rho0.reshape(self.Nsys**2)
    norm0 = LA.norm(vec_rho0, 2)
    statevec = vec_rho0 / norm0
    return statevec, norm0

qc_simulation_kraus(time_arr, shots=1024, Kraus=None, Gprop=None, tolk=1e-05, tolo=1e-05, **kwargs)

Perform quantum simulation using the Kraus operator representation.

This method simulates the quantum system dynamics over a series of time steps using a Kraus operator-based approach. It constructs quantum circuits for each Kraus operator and initial state, runs the simulation using Qiskit's Estimator, and accumulates the measurement results.

Parameters:

Name Type Description Default
time_arr List[float]

List of time steps for simulation.

 required
shots int

Number of shots for each measurement. Defaults to 1024.

 1024
Kraus Optional[Dict[int, List[ndarray]]]

Dictionary mapping time step index to a list of Kraus operators. If None, Kraus operators are generated from the propagator. Defaults to None.

 None
Gprop Optional[List[ndarray]]

Propagator matrix (or list of matrices) for simulation. If None, it will be calculated. Defaults to None.

 None
tolk float

Tolerance for generating Kraus operators. Defaults to 1e-5.

 1e-05
tolo float

Tolerance for observable decomposition. Defaults to 1e-5.

 1e-05
**kwargs Any

Additional keyword arguments for propagator calculation.

 {}

Returns:

Type Description
ndarray

np.ndarray: Array containing the quantum simulation results.
Source code in src/qflux/open_systems/quantum_simulation.py 
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
def qc_simulation_kraus(
    self,
    time_arr: List[float],
    shots: int = 1024,
    Kraus: Optional[Dict[int, List[np.ndarray]]] = None,
    Gprop: Optional[List[np.ndarray]] = None,
    tolk: float = 1e-5,
    tolo: float = 1e-5,
    **kwargs: Any
) -> np.ndarray:
    """
    Perform quantum simulation using the Kraus operator representation.

    This method simulates the quantum system dynamics over a series of time steps using a Kraus operator-based approach.
    It constructs quantum circuits for each Kraus operator and initial state, runs the simulation using Qiskit's Estimator,
    and accumulates the measurement results.

    Args:
        time_arr (List[float]): List of time steps for simulation.
        shots (int, optional): Number of shots for each measurement. Defaults to 1024.
        Kraus (Optional[Dict[int, List[np.ndarray]]], optional): Dictionary mapping time step index to a list of Kraus operators.
            If None, Kraus operators are generated from the propagator. Defaults to None.
        Gprop (Optional[List[np.ndarray]], optional): Propagator matrix (or list of matrices) for simulation.
            If None, it will be calculated. Defaults to None.
        tolk (float, optional): Tolerance for generating Kraus operators. Defaults to 1e-5.
        tolo (float, optional): Tolerance for observable decomposition. Defaults to 1e-5.
        **kwargs: Additional keyword arguments for propagator calculation.

    Returns:
        np.ndarray: Array containing the quantum simulation results.
    """
    nsteps = len(time_arr)

    # Generate Kraus operators if not provided.
    if Kraus is None:
        Kraus = {}
        if Gprop is None:
            print('Calculating the propagator')
            Gprop = self.Gt_matrix_expo(time_arr, **kwargs)
        print('Generating the Kraus operators')
        for i in range(nsteps):
            print('At step', i, 'of', nsteps)
            Kraus[i] = gen_Kraus_list(Gprop[i], self.Nsys, tol=tolk)
    print('Kraus operator generation complete')

    # Perform Qiskit simulation using the Estimator.
    estimator = Estimator()

    statevec, prob = self.init_statevec_Kraus()
    n_inistate = len(statevec)
    print('Number of initial states in the density matrix:', n_inistate)
    print('Probabilities:', prob)

    # Obtain the Qiskit observable.
    obs_q = self._get_qiskit_observable(Isdilate=True, tol=tolo)

    print('Starting quantum simulation')
    result_simulation = np.zeros(nsteps, dtype=np.float64)

    for i in range(nsteps):
        print('Simulation step', i, 'of', nsteps)
        current_kraus_list = Kraus[i]
        print('Number of Kraus operators:', len(current_kraus_list))
        for kraus_op in current_kraus_list:
            for istate in range(n_inistate):
                qc = self._create_circuit(kraus_op, statevec[istate], Isscale=False)
                result = estimator.run(qc, obs_q, shots=shots).result()
                result_simulation[i] += result.values[0] * prob[istate]

    return result_simulation

qc_simulation_vecdens(time_arr, shots=1024, backend=AerSimulator(), Gprop=None, **kwargs)

Perform quantum simulation using the vectorized density matrix representation.

This method simulates the quantum system dynamics over a series of time steps by constructing circuits based on the vectorized density matrix representation, performing measurements, and processing the results.

Parameters:

Name Type Description Default
time_arr List[float]

List of time steps for simulation.

 required
shots int

Number of measurement shots. Defaults to 1024.

 1024
backend Any

Quantum simulation backend. Defaults to AerSimulator().

 AerSimulator()
Gprop Optional[List[ndarray]]

Propagator matrix (or list of matrices) for simulation. If None, it will be calculated. Defaults to None.

 None
**kwargs Any

Additional keyword arguments for propagator calculation.

 {}

Returns:

Type Description
ndarray

np.ndarray: Array containing the quantum simulation results.
Source code in src/qflux/open_systems/quantum_simulation.py 
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
def qc_simulation_vecdens(
    self,
    time_arr: List[float],
    shots: int = 1024,
    backend: Any = AerSimulator(),
    Gprop: Optional[List[np.ndarray]] = None,
    **kwargs: Any
) -> np.ndarray:
    """
    Perform quantum simulation using the vectorized density matrix representation.

    This method simulates the quantum system dynamics over a series of time steps by constructing circuits
    based on the vectorized density matrix representation, performing measurements, and processing the results.

    Args:
        time_arr (List[float]): List of time steps for simulation.
        shots (int, optional): Number of measurement shots. Defaults to 1024.
        backend (Any, optional): Quantum simulation backend. Defaults to AerSimulator().
        Gprop (Optional[List[np.ndarray]], optional): Propagator matrix (or list of matrices) for simulation.
            If None, it will be calculated. Defaults to None.
        **kwargs: Additional keyword arguments for propagator calculation.

    Returns:
        np.ndarray: Array containing the quantum simulation results.
    """
    if Gprop is None:
        Gprop = self.Gt_matrix_expo(time_arr, **kwargs)

    nsteps = len(time_arr)

    if self.count_str is None:
        print("Error: count_str is not assigned")

    n_bitstr = len(self.count_str)
    statevec, norm0 = self.init_statevec_vecdens()
    result = np.zeros((nsteps, n_bitstr), dtype=np.float64)

    for i in range(nsteps):
        if i % 100 == 0:
            print('Quantum simulation step', i)
        Gt = Gprop[i]
        circuit, norm = self._create_circuit(Gt, statevec, Isscale=True)
        circuit.measure(range(self.Nqb + 1), range(self.Nqb + 1))
        if self.dilation_method == 'SVD-Walsh':
            circuit = transpile(circuit, backend)
        counts = backend.run(circuit, shots=shots).result().get_counts()
        for j in range(n_bitstr):
            bitstr = self.count_str[j]
            if bitstr in counts:
                result[i, j] = np.sqrt(counts[bitstr] / shots) * norm * norm0
            else:
                print('At time', i, 'with shots =', shots, "no counts for", bitstr)

    return result

set_count_str(count_str)

Set the counting bit string for measurement.

Parameters:

Name Type Description Default
count_str List[str]

The counting bit string.

 required
Source code in src/qflux/open_systems/quantum_simulation.py 
121
122
123
124
125
126
127
128
def set_count_str(self, count_str: List[str]) -> None:
    """
    Set the counting bit string for measurement.

    Args:
        count_str (List[str]): The counting bit string.
    """
    self.count_str = count_str

set_dilation_method(method)

Set the dilation method for quantum simulation.

Parameters:

Name Type Description Default
method str

The dilation method, e.g., 'Sz-Nagy', 'SVD', or 'SVD-Walsh'.

 required
Source code in src/qflux/open_systems/quantum_simulation.py 
112
113
114
115
116
117
118
119
def set_dilation_method(self, method: str) -> None:
    """
    Set the dilation method for quantum simulation.

    Args:
        method (str): The dilation method, e.g., 'Sz-Nagy', 'SVD', or 'SVD-Walsh'.
    """
    self.dilation_method = method

set_observable(observable)

Set the observable for the quantum simulation.

Parameters:

Name Type Description Default
observable ndarray

The observable matrix.

 required
Source code in src/qflux/open_systems/quantum_simulation.py 
130
131
132
133
134
135
136
137
def set_observable(self, observable: np.ndarray) -> None:
    """
    Set the observable for the quantum simulation.

    Args:
        observable (np.ndarray): The observable matrix.
    """
    self.observable = observable

DVR_grid

Class for Discrete Variable Representation (DVR) grid methods.

This class handles grid-based representations for systems where the potential is expressed on grid points.

Attributes:

Name Type Description
Ngrid int

Number of grid points.
xmin float

Minimum value of the grid.
xmax float

Maximum value of the grid.
mass float

Mass of the particle.
xgrid ndarray

Array of grid points in position space.
dx float

Spacing between grid points.
dk float

Spacing in momentum space.
kgrid ndarray

Array of grid points in momentum space.
ak2 ndarray

Kinetic energy array in momentum space.
hamk ndarray

Kinetic Hamiltonian matrix in position space.
potential Optional[ndarray]

Potential energy array on the grid.
Source code in src/qflux/open_systems/numerical_methods.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
class DVR_grid:
    """Class for Discrete Variable Representation (DVR) grid methods.

    This class handles grid-based representations for systems where the potential is expressed on grid points.

    Attributes:
        Ngrid (int): Number of grid points.
        xmin (float): Minimum value of the grid.
        xmax (float): Maximum value of the grid.
        mass (float): Mass of the particle.
        xgrid (np.ndarray): Array of grid points in position space.
        dx (float): Spacing between grid points.
        dk (float): Spacing in momentum space.
        kgrid (np.ndarray): Array of grid points in momentum space.
        ak2 (np.ndarray): Kinetic energy array in momentum space.
        hamk (np.ndarray): Kinetic Hamiltonian matrix in position space.
        potential (Optional[np.ndarray]): Potential energy array on the grid.
    """

    def __init__(self, xmin: float, xmax: float, Ngrid: int, mass: float) -> None:
        """
        Initialize the DVR_grid instance.

        Args:
            xmin (float): Minimum x-value.
            xmax (float): Maximum x-value.
            Ngrid (int): Number of grid points.
            mass (float): Mass of the particle.
        """
        self.Ngrid: int = Ngrid
        self.xmin: float = xmin
        self.xmax: float = xmax
        self.mass: float = mass

        # Set up the position grid.
        self.xgrid: np.ndarray = np.array([])
        self._set_xgrid()
        self.dx: float = self.xgrid[1] - self.xgrid[0]

        # Set up the momentum grid.
        self.dk: float = 2.0 * np.pi / (self.Ngrid * self.dx)
        self.kgrid: np.ndarray = np.array([])
        self.ak2: np.ndarray = np.array([])  # Kinetic energy array.
        self.hamk: np.ndarray = np.array([])  # Kinetic Hamiltonian matrix.
        self._set_kinet_ham()

        # Potential energy array (to be set later).
        self.potential: Optional[np.ndarray] = None

    def _set_xgrid(self) -> None:
        """
        Set up the position space grid.

        Initializes the `xgrid` attribute using a linear space between `xmin` and `xmax`.
        """
        self.xgrid = np.linspace(self.xmin, self.xmax, self.Ngrid)

    def set_potential(self, func_pot: Callable[[float], float]) -> None:
        """
        Set up the potential energy array on the grid.

        Args:
            func_pot (Callable[[float], float]): Function that returns the potential value at a given x.
        """
        self.potential = np.zeros_like(self.xgrid)
        for i in range(self.Ngrid):
            self.potential[i] = func_pot(self.xgrid[i])

    def _set_kinet_ham(self) -> None:
        """
        Set up the kinetic Hamiltonian matrix in position space.

        This method computes the momentum grid and the corresponding kinetic energy values,
        and constructs the kinetic Hamiltonian matrix in position space using a Fourier transform.
        """
        self.kgrid = np.zeros(self.Ngrid, dtype=np.float64)
        self.ak2 = np.zeros(self.Ngrid, dtype=np.float64)

        coef_k: float = pa.hbar**2 / (2.0 * self.mass)

        for i in range(self.Ngrid):
            if i < self.Ngrid // 2:
                self.kgrid[i] = i * self.dk
            else:
                self.kgrid[i] = -(self.Ngrid - i) * self.dk
            self.ak2[i] = coef_k * self.kgrid[i]**2

        akx0: np.ndarray = sfft.ifft(self.ak2)
        self.hamk = np.zeros((self.Ngrid, self.Ngrid), dtype=np.complex128)

        for i in range(self.Ngrid):
            for j in range(self.Ngrid):
                if i < j:
                    self.hamk[i, j] = akx0[i - j].conj()
                else:
                    self.hamk[i, j] = akx0[i - j]

    def get_eig_state(self, Nstate: int) -> Tuple[np.ndarray, np.ndarray]:
        """
        Get the eigenstates for the potential in x-space.

        Args:
            Nstate (int): Number of eigenstates to output.

        Returns:
            Tuple[np.ndarray, np.ndarray]: A tuple containing:
                - Eigenvalues (np.ndarray) for the first `Nstate` states.
                - Eigenvectors (np.ndarray) for the first `Nstate` states, normalized by sqrt(dx).
        """
        Mata: np.ndarray = self.hamk.copy()
        for i in range(self.Ngrid):
            Mata[i, i] += self.potential[i]

        val, arr = LA.eigh(Mata)
        return val[:Nstate], arr[:, :Nstate] / np.sqrt(self.dx)

    def x2k_wave(self, psi: np.ndarray) -> np.ndarray:
        """
        Transform the wavefunction from position space to momentum space.

        Args:
            psi (np.ndarray): Wavefunction in position space.

        Returns:
            np.ndarray: Wavefunction in momentum space.
        """
        return tb.x2k_wave(self.dx, psi)

    def k2x_wave(self, psik: np.ndarray) -> np.ndarray:
        """
        Transform the wavefunction from momentum space to position space.

        Args:
            psik (np.ndarray): Wavefunction in momentum space.

        Returns:
            np.ndarray: Wavefunction in position space.
        """
        return tb.k2x_wave(self.dx, psik)

__init__(xmin, xmax, Ngrid, mass)

Initialize the DVR_grid instance.

Parameters:

Name Type Description Default
xmin float

Minimum x-value.

 required
xmax float

Maximum x-value.

 required
Ngrid int

Number of grid points.

 required
mass float

Mass of the particle.

 required
Source code in src/qflux/open_systems/numerical_methods.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def __init__(self, xmin: float, xmax: float, Ngrid: int, mass: float) -> None:
    """
    Initialize the DVR_grid instance.

    Args:
        xmin (float): Minimum x-value.
        xmax (float): Maximum x-value.
        Ngrid (int): Number of grid points.
        mass (float): Mass of the particle.
    """
    self.Ngrid: int = Ngrid
    self.xmin: float = xmin
    self.xmax: float = xmax
    self.mass: float = mass

    # Set up the position grid.
    self.xgrid: np.ndarray = np.array([])
    self._set_xgrid()
    self.dx: float = self.xgrid[1] - self.xgrid[0]

    # Set up the momentum grid.
    self.dk: float = 2.0 * np.pi / (self.Ngrid * self.dx)
    self.kgrid: np.ndarray = np.array([])
    self.ak2: np.ndarray = np.array([])  # Kinetic energy array.
    self.hamk: np.ndarray = np.array([])  # Kinetic Hamiltonian matrix.
    self._set_kinet_ham()

    # Potential energy array (to be set later).
    self.potential: Optional[np.ndarray] = None

_set_kinet_ham()

Set up the kinetic Hamiltonian matrix in position space.

This method computes the momentum grid and the corresponding kinetic energy values, and constructs the kinetic Hamiltonian matrix in position space using a Fourier transform.

Source code in src/qflux/open_systems/numerical_methods.py 
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def _set_kinet_ham(self) -> None:
    """
    Set up the kinetic Hamiltonian matrix in position space.

    This method computes the momentum grid and the corresponding kinetic energy values,
    and constructs the kinetic Hamiltonian matrix in position space using a Fourier transform.
    """
    self.kgrid = np.zeros(self.Ngrid, dtype=np.float64)
    self.ak2 = np.zeros(self.Ngrid, dtype=np.float64)

    coef_k: float = pa.hbar**2 / (2.0 * self.mass)

    for i in range(self.Ngrid):
        if i < self.Ngrid // 2:
            self.kgrid[i] = i * self.dk
        else:
            self.kgrid[i] = -(self.Ngrid - i) * self.dk
        self.ak2[i] = coef_k * self.kgrid[i]**2

    akx0: np.ndarray = sfft.ifft(self.ak2)
    self.hamk = np.zeros((self.Ngrid, self.Ngrid), dtype=np.complex128)

    for i in range(self.Ngrid):
        for j in range(self.Ngrid):
            if i < j:
                self.hamk[i, j] = akx0[i - j].conj()
            else:
                self.hamk[i, j] = akx0[i - j]

_set_xgrid()

Set up the position space grid.

Initializes the xgrid attribute using a linear space between xmin and xmax.

Source code in src/qflux/open_systems/numerical_methods.py 
214
215
216
217
218
219
220
def _set_xgrid(self) -> None:
    """
    Set up the position space grid.

    Initializes the `xgrid` attribute using a linear space between `xmin` and `xmax`.
    """
    self.xgrid = np.linspace(self.xmin, self.xmax, self.Ngrid)

get_eig_state(Nstate)

Get the eigenstates for the potential in x-space.

Parameters:

Name Type Description Default
Nstate int

Number of eigenstates to output.

 required

Returns:

Type Description
Tuple[ndarray, ndarray]

Tuple[np.ndarray, np.ndarray]: A tuple containing: - Eigenvalues (np.ndarray) for the first Nstate states. - Eigenvectors (np.ndarray) for the first Nstate states, normalized by sqrt(dx).
Source code in src/qflux/open_systems/numerical_methods.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
def get_eig_state(self, Nstate: int) -> Tuple[np.ndarray, np.ndarray]:
    """
    Get the eigenstates for the potential in x-space.

    Args:
        Nstate (int): Number of eigenstates to output.

    Returns:
        Tuple[np.ndarray, np.ndarray]: A tuple containing:
            - Eigenvalues (np.ndarray) for the first `Nstate` states.
            - Eigenvectors (np.ndarray) for the first `Nstate` states, normalized by sqrt(dx).
    """
    Mata: np.ndarray = self.hamk.copy()
    for i in range(self.Ngrid):
        Mata[i, i] += self.potential[i]

    val, arr = LA.eigh(Mata)
    return val[:Nstate], arr[:, :Nstate] / np.sqrt(self.dx)

k2x_wave(psik)

Transform the wavefunction from momentum space to position space.

Parameters:

Name Type Description Default
psik ndarray

Wavefunction in momentum space.

 required

Returns:

Type Description
ndarray

np.ndarray: Wavefunction in position space.
Source code in src/qflux/open_systems/numerical_methods.py 
293
294
295
296
297
298
299
300
301
302
303
def k2x_wave(self, psik: np.ndarray) -> np.ndarray:
    """
    Transform the wavefunction from momentum space to position space.

    Args:
        psik (np.ndarray): Wavefunction in momentum space.

    Returns:
        np.ndarray: Wavefunction in position space.
    """
    return tb.k2x_wave(self.dx, psik)

set_potential(func_pot)

Set up the potential energy array on the grid.

Parameters:

Name Type Description Default
func_pot Callable[[float], float]

Function that returns the potential value at a given x.

 required
Source code in src/qflux/open_systems/numerical_methods.py 
222
223
224
225
226
227
228
229
230
231
def set_potential(self, func_pot: Callable[[float], float]) -> None:
    """
    Set up the potential energy array on the grid.

    Args:
        func_pot (Callable[[float], float]): Function that returns the potential value at a given x.
    """
    self.potential = np.zeros_like(self.xgrid)
    for i in range(self.Ngrid):
        self.potential[i] = func_pot(self.xgrid[i])

x2k_wave(psi)

Transform the wavefunction from position space to momentum space.

Parameters:

Name Type Description Default
psi ndarray

Wavefunction in position space.

 required

Returns:

Type Description
ndarray

np.ndarray: Wavefunction in momentum space.
Source code in src/qflux/open_systems/numerical_methods.py 
281
282
283
284
285
286
287
288
289
290
291
def x2k_wave(self, psi: np.ndarray) -> np.ndarray:
    """
    Transform the wavefunction from position space to momentum space.

    Args:
        psi (np.ndarray): Wavefunction in position space.

    Returns:
        np.ndarray: Wavefunction in momentum space.
    """
    return tb.x2k_wave(self.dx, psi)