Skip to content

Source Code

qflux.closed_systems

DynamicsCS

Class for closed-system dynamics. All input parameters must be in atomic units to ensure consistency. Please be sure to convert your parameters to atomic units prior to instantiation.

Source code in src/qflux/closed_systems/classical_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
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
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
class DynamicsCS:
    """
        Class for closed-system dynamics. **All input parameters must be in
        atomic units to ensure consistency. Please be sure to convert your
        parameters to atomic units prior to instantiation.**

    """

    def __init__(self, n_basis: int = 128, xo: float = 1.0, po: float = 0.0, mass: float = 1.0,
                 omega: float = 1.0) -> None:
        """
        Args:
            n_basis (int): Number of states to include in the chosen representation. If basis
                = 'ladder', this is the Fock cutoff and defines the number of states
                used for representing the ladder operators. If basis = 'coordinate',
                this defines the number of points for the position and momenta.

            xo (float, optional): Defines the displacement of the initial state in the position
                coordinate. Default is 1.0 Bohr.

            po (float, optional): Defines the displacement of the initial state in the position
                coordinate. Default is 1.0 au.

            mass (float, optional): Defines the mass of the particle/system of interest.
                Default is 1.0 au.

            omega (float, optional): Frequency of harmonic oscillator.
                Default is 1.0 au.

        """
        #--------- Required Attributes Populated During Execution ----------#
        self.n_basis                   = n_basis
        self.xo                        = xo
        self.po                        = po
        self.mass                      = mass
        self.hbar                      = 1.0
        self.omega                     = omega
        #--------- Below are Attributes Populated During Execution ---------#
        self.total_time                = 0.
        self.n_tsteps                  = 0.
        self._KE_op                    = None
        self._PE_op                    = None
        self.H_op                      = None
        self.prop_KE_op                = None
        self.prop_PE_op                = None
        self.prop_H_op                 = None
        # Grid operators
        self._KE_grid                  = None
        self._PE_grid                  = None
        self.H_grid                    = None
        self.PE_prop_grid              = None
        self.KE_prop_grid              = None


    def _get_xgrid(self, x_min: float, x_max: float) -> None:
        """
        Populate the `self.x_grid` and `self.dx` attributes. This function
        generates an array of `self.n_basis` evenly spaced values between
        `x_min` and `x_max`.

        Args:
            x_min (float): Minimum value of x-coordinates
            x_max (float): Maximum value of x-coordinates

        Returns:
            self.dx (float): Spacing between points in the x-coordinate grid.
            self.xgrid (array_like): Array of grid points from x_min to x_max with spacing of dx
        """
        dx = (x_max - x_min) / self.n_basis
        x_grid = np.arange(-self.n_basis / 2, self.n_basis / 2) * dx
        self.dx = dx
        self.x_grid = x_grid
        return


    def _get_pgrid(self, x_min: float, x_max: float, reorder: bool = True) -> None:
        """
        Populate the `self.p_grid` and `self.dp` attributes. This function
        generates an array of `self.n_basis` evenly spaced values.

        Args:
            x_min (float): Minimum value of x-coordinates
            x_max (float): Maximum value of x-coordinates
            reorder (bool): Boolean flag to determine whether points should be reordered to be
                compatible with the FFT routine or not.

        Returns:
            self.dp (float): Spacing between points in the p-coordinate grid.
            self.pgrid (array_like): Array of momentum grid points
        """
        dp = 2 * np.pi / (x_max - x_min)
        pmin = -dp * self.n_basis / 2
        pmax = dp * self.n_basis / 2
        plus_pgrid = np.linspace(0, pmax, self.n_basis // 2 + 1)
        minus_pgrid = - np.flip(np.copy(plus_pgrid))
        if reorder:
            p_grid = np.concatenate((plus_pgrid[:-1], minus_pgrid[:-1]))
        else:
            p_grid = np.concatenate((minus_pgrid, plus_pgrid))
        self.p_grid = p_grid
        self.dp = dp
        return


    def set_coordinate_operators(self, x_min: float = -7., x_max: float = 7., reorder_p: bool = True) -> None:
        """
        Populate the `self.x_grid`, `self.p_grid`, `self.dx`, and `self.dp`
        attributes. This functions generates an array of `self.n_basis`
        evenly spaced values.

        Args:
            x_min : float
                Minimum value of x-coordinates
            x_max : float
                Maximum value of x-coordinates
            reorder_p : bool
                Boolean flag to determine whether momentum values should be
                reordered to be compatible with the FFT routine or not.

        Returns:
            self.dx : float
                Spacing between points in the x-coordinate grid.
            self.xgrid : array_like
                Array of x-values
            self.dp : float
                Spacing between points in the p-coordinate grid.
            self.pgrid : array_like
                Array of p-values
        """
        self._get_xgrid(x_min, x_max)
        self._get_pgrid(x_min, x_max, reorder=reorder_p)
        return


    def initialize_operators(self):
        """
            Function to initialize core operators in the chosen basis.

        """

        self.a_op = qt.destroy(self.n_basis)
        self.x_op = qt.position(self.n_basis)
        self.p_op = qt.momentum(self.n_basis)
        return


    def _set_hamiltonian_grid(self, potential_type: str = 'harmonic', **kwargs):
        if potential_type == 'harmonic':

            # Set attributes for the coordinate basis
            self._PE_grid = self.mass * self.omega ** 2 * self.x_grid ** 2 / 2.
            self._KE_grid = self.p_grid ** 2 / (2. * self.mass)
            self.H_grid = self._PE_grid + self._KE_grid
        elif potential_type == 'quartic':
            if kwargs:
                if 'a0' in kwargs:
                    a0 = kwargs['a0']
                if 'a1' in kwargs:
                    a1 = kwargs['a1']
                if 'a2' in kwargs:
                    a2 = kwargs['a2']
                if 'a3' in kwargs:
                    a3 = kwargs['a3']
                if 'a4' in kwargs:
                    a4 = kwargs['a4']
                if 'x0' in kwargs:
                    x0 = kwargs['x0']
                # Assume that all inputs have the proper atomic units:
                cf = 1.0
                xi = self.x_op

            else:
                # Define relevant parameters
                cf = convert_eV_to_au(1.)
                x0 = 1.9592
                a0 = 0.0
                a1 = 0.429
                a2 = -1.126
                a3 = -0.143
                a4 = 0.563
                # Do calculation for ladder basis
                xi = self.x_grid / x0
            self._PE_grid = cf * (a0 + a1 * xi + a2 * xi ** 2 + a3 * xi ** 3 + a4 * xi ** 4)
            self._KE_grid = self.p_grid ** 2 / (2. * self.mass)
            self.H_grid = self._PE_grid + self._KE_grid
        return


    def _set_hamiltonian_qt(self, potential_type: str = 'harmonic', **kwargs):
        if potential_type == 'harmonic':
            # Set attributes for the ladder basis
            self.H_op = self.hbar * self.omega * (self.a_op.dag() * self.a_op + 0.5)
            self._KE_op = self.p_op ** 2 / (2. * self.mass)
            self._PE_op = self.mass * self.omega ** 2 * self.x_op ** 2 / 0.5
            self.H_xp_op = self._PE_op + self._KE_op
        elif potential_type == 'quartic':
            if kwargs:
                if 'a0' in kwargs:
                    a0 = kwargs['a0']
                if 'a1' in kwargs:
                    a1 = kwargs['a1']
                if 'a2' in kwargs:
                    a2 = kwargs['a2']
                if 'a3' in kwargs:
                    a3 = kwargs['a3']
                if 'a4' in kwargs:
                    a4 = kwargs['a4']
                if 'x0' in kwargs:
                    x0 = kwargs['x0']
                # Assume that all inputs have the proper atomic units:
                cf = 1.0
                xi = self.x_op

            else:
                # Define relevant parameters
                cf = convert_eV_to_au(1.)
                x0 = 1.9592
                a0 = 0.0
                a1 = 0.429
                a2 = -1.126
                a3 = -0.143
                a4 = 0.563
                # Do calculation for ladder basis
                xi = self.x_op / x0
            self.x0 = x0
            self._PE_op = cf * (a0 + a1 * xi + a2 * xi ** 2 + a3 * xi ** 3 + a4 * xi ** 4)
            self._KE_op = self.p_op ** 2 / (2. * self.mass)
            self.H_op = self._PE_op + self._KE_op
            return

    def set_hamiltonian(self, potential_type: str = 'harmonic', **kwargs):
        """
        Function to define Hamiltonian.

        Args:
            potential_type : str
                String defining the type of potential energy surface.
                Available options are: ('harmonic', 'quartic', ...)

                Note: You can manually define your potential energy using the functions:
                    - set_H_grid_with_custom_potential
                    - set_H_op_with_custom_potential

        """

        if potential_type == 'harmonic':
            self._set_hamiltonian_grid(potential_type=potential_type, **kwargs)
            self._set_hamiltonian_qt(potential_type=potential_type, **kwargs)
        elif potential_type == 'quartic':
            self._set_hamiltonian_grid(potential_type=potential_type, **kwargs)
            self._set_hamiltonian_qt(potential_type=potential_type, **kwargs)
        else:
            print('Error, this potential type has not yet been implemented!')
            print('Set your parameters with the custom functions!')
        return


    def set_H_grid_with_custom_potential(self, custom_function: Callable, **kwargs):
        """
        Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

        Args:
            custom_function (Callable): Function that defines the custom potential
                energy. Must return an array

        """
        potential = custom_function(**kwargs)
        self._PE_grid = potential
        self._KE_grid = self.p_grid ** 2 / (2. * self.mass)
        self.H_grid = self._PE_grid + self._KE_grid
        return


    def set_H_op_with_custom_potential(self, custom_function: Callable, **kwargs):
        """
        Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

        Args:
            custom_function (Callable): Function that defines the potential
                energy in terms of qutip QObj operators. Must return a qutip.Qobj
        """
        potential = custom_function(**kwargs)
        self._PE_op = potential
        self._KE_op = self.p_op ** 2 / (2. * self.mass)
        self.H_op = self._PE_op + self._KE_op
        return


    def set_initial_state(self, wfn_omega: float = 1.0):
        """
        Function to define the initial state. By default, a coherent state is
        used as the initial state defined in the basis chosen upon instantiation

        Args:
            wfn_omega (float, optional): Defines the frequency/width of the initial state.
                Default is 1.0 au.
        """

        alpha_val = (self.xo + 1j * self.po) / np.sqrt(2)
        psio = qt.coherent(self.n_basis, alpha=alpha_val)
        # Now populate the initial state in the grid basis
        normalization = (self.mass * wfn_omega / np.pi / self.hbar) ** (0.25)
        exponential = np.exp(-1 * (self.mass * wfn_omega / self.hbar / 2) *
                             ((self.x_grid - self.xo) ** 2)
                             + 1j * self.po * self.x_grid / self.hbar
                             )

        coherent_state = normalization * exponential
        # Set the attributes
        self.psio_grid = coherent_state
        self.psio_op = psio
        return



    def custom_grid_state_initialization(self, function_name: Callable, **kwargs):
        """
        Function to allow for customized grid state initialization.

        Args:
            function_name (Callable): name of user-defined function that returns
                the initial state. Must return an array
        """

        self.psio_grid = function_name(**kwargs)
        return

    def custom_ladder_state_initialization(self, function_name: Callable, **kwargs):
        """
        Function to allow for customized ladder state initialization.

        Args:
            function_name (Callable): name of user-defined function that returns
                the initial state. Must return a qutip.Qobj.
        """

        self.psio_op = function_name(**kwargs)
        return

    def set_propagation_time(self, total_time: float, n_tsteps: int):
        """
        Function to define the propagation time, an array of times from
        t=0 to total_time, with n_tsteps equally-spaced steps.

        Args:
        total_time : float
            The total time for which we wish to compute the dynamics.
        n_tsteps : int
            The number of equally-spaced time steps used to compute the dynamics

        Returns:
        self.tlist : array-like

        """

        self.tlist = np.linspace(0., total_time, n_tsteps+1)
        self.dt = self.tlist[1] - self.tlist[0]
        self.n_tsteps = n_tsteps
        return


    def propagate_qt(self, solver_options : dict = None):
        """
        Function used to propagate with qutip.

        Args:
            solver_options (dict): A dictionary of arguments to pass to the qutip.sesolve function

        Returns:
            dynamics_results (array-like): array containing the propagated state

        """

        options = {'nsteps': len(self.tlist),
                    'progress_bar': True}

        if solver_options:
            for key in solver_options:
                options[key] = solver_options[key]

        results = qt.sesolve(self.H_op, self.psio_op, self.tlist,
                             options=options)

        self.dynamics_results_op = results.states
        return


    def propagate_SOFT(self):
        """
        Function used to propagate with the 2nd-Order Trotter Expansion.

        $$
        e^{- \frac{i}{\\hbar} H t} \approx e^{- \frac{i}{\\hbar} V t/2} e^{- \frac{i}{\\hbar} T t} e^{- \frac{i}{\\hbar} V t/2} + \\mathcal{O}^{3}
        $$

        Returns:
            dynamics_results_grid (array-like): array containing the propagated state
                                                shape (n_tsteps x self.n_basis)

        """
        self.tau = self.tlist[1] - self.tlist[0]
        PE_prop = np.exp(-1.0j * self._PE_grid / 2 * self.tau / self.hbar)
        KE_prop = np.exp(-1.0j * self._KE_grid * self.tau / self.hbar)

        self.PE_prop_grid = PE_prop
        self.KE_prop_grid = KE_prop

        propagated_states = [self.psio_grid]
        psi_t = self.psio_grid
        for ii in range(1, len(self.tlist)):
            psi_t_position_grid = PE_prop * psi_t
            psi_t_momentum_grid = KE_prop * np.fft.fft(psi_t_position_grid, norm="ortho")
            psi_t = PE_prop * np.fft.ifft(psi_t_momentum_grid, norm="ortho")
            propagated_states.append(psi_t)

        self.dynamics_results_grid = np.asarray(propagated_states)
        return

__init__(n_basis=128, xo=1.0, po=0.0, mass=1.0, omega=1.0)

Parameters:

Name Type Description Default
n_basis int

Number of states to include in the chosen representation. If basis = 'ladder', this is the Fock cutoff and defines the number of states used for representing the ladder operators. If basis = 'coordinate', this defines the number of points for the position and momenta.

 128
xo float

Defines the displacement of the initial state in the position coordinate. Default is 1.0 Bohr.

 1.0
po float

Defines the displacement of the initial state in the position coordinate. Default is 1.0 au.

 0.0
mass float

Defines the mass of the particle/system of interest. Default is 1.0 au.

 1.0
omega float

Frequency of harmonic oscillator. Default is 1.0 au.

 1.0
Source code in src/qflux/closed_systems/classical_methods.py 
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
def __init__(self, n_basis: int = 128, xo: float = 1.0, po: float = 0.0, mass: float = 1.0,
             omega: float = 1.0) -> None:
    """
    Args:
        n_basis (int): Number of states to include in the chosen representation. If basis
            = 'ladder', this is the Fock cutoff and defines the number of states
            used for representing the ladder operators. If basis = 'coordinate',
            this defines the number of points for the position and momenta.

        xo (float, optional): Defines the displacement of the initial state in the position
            coordinate. Default is 1.0 Bohr.

        po (float, optional): Defines the displacement of the initial state in the position
            coordinate. Default is 1.0 au.

        mass (float, optional): Defines the mass of the particle/system of interest.
            Default is 1.0 au.

        omega (float, optional): Frequency of harmonic oscillator.
            Default is 1.0 au.

    """
    #--------- Required Attributes Populated During Execution ----------#
    self.n_basis                   = n_basis
    self.xo                        = xo
    self.po                        = po
    self.mass                      = mass
    self.hbar                      = 1.0
    self.omega                     = omega
    #--------- Below are Attributes Populated During Execution ---------#
    self.total_time                = 0.
    self.n_tsteps                  = 0.
    self._KE_op                    = None
    self._PE_op                    = None
    self.H_op                      = None
    self.prop_KE_op                = None
    self.prop_PE_op                = None
    self.prop_H_op                 = None
    # Grid operators
    self._KE_grid                  = None
    self._PE_grid                  = None
    self.H_grid                    = None
    self.PE_prop_grid              = None
    self.KE_prop_grid              = None

_get_pgrid(x_min, x_max, reorder=True)

Populate the self.p_grid and self.dp attributes. This function generates an array of self.n_basis evenly spaced values.

Parameters:

Name Type Description Default
x_min float

Minimum value of x-coordinates

 required
x_max float

Maximum value of x-coordinates

 required
reorder bool

Boolean flag to determine whether points should be reordered to be compatible with the FFT routine or not.

 True

Returns:

Type Description
None

self.dp (float): Spacing between points in the p-coordinate grid.
None

self.pgrid (array_like): Array of momentum grid points
Source code in src/qflux/closed_systems/classical_methods.py 
 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
def _get_pgrid(self, x_min: float, x_max: float, reorder: bool = True) -> None:
    """
    Populate the `self.p_grid` and `self.dp` attributes. This function
    generates an array of `self.n_basis` evenly spaced values.

    Args:
        x_min (float): Minimum value of x-coordinates
        x_max (float): Maximum value of x-coordinates
        reorder (bool): Boolean flag to determine whether points should be reordered to be
            compatible with the FFT routine or not.

    Returns:
        self.dp (float): Spacing between points in the p-coordinate grid.
        self.pgrid (array_like): Array of momentum grid points
    """
    dp = 2 * np.pi / (x_max - x_min)
    pmin = -dp * self.n_basis / 2
    pmax = dp * self.n_basis / 2
    plus_pgrid = np.linspace(0, pmax, self.n_basis // 2 + 1)
    minus_pgrid = - np.flip(np.copy(plus_pgrid))
    if reorder:
        p_grid = np.concatenate((plus_pgrid[:-1], minus_pgrid[:-1]))
    else:
        p_grid = np.concatenate((minus_pgrid, plus_pgrid))
    self.p_grid = p_grid
    self.dp = dp
    return

_get_xgrid(x_min, x_max)

Populate the self.x_grid and self.dx attributes. This function generates an array of self.n_basis evenly spaced values between x_min and x_max.

Parameters:

Name Type Description Default
x_min float

Minimum value of x-coordinates

 required
x_max float

Maximum value of x-coordinates

 required

Returns:

Type Description
None

self.dx (float): Spacing between points in the x-coordinate grid.
None

self.xgrid (array_like): Array of grid points from x_min to x_max with spacing of dx
Source code in src/qflux/closed_systems/classical_methods.py 
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def _get_xgrid(self, x_min: float, x_max: float) -> None:
    """
    Populate the `self.x_grid` and `self.dx` attributes. This function
    generates an array of `self.n_basis` evenly spaced values between
    `x_min` and `x_max`.

    Args:
        x_min (float): Minimum value of x-coordinates
        x_max (float): Maximum value of x-coordinates

    Returns:
        self.dx (float): Spacing between points in the x-coordinate grid.
        self.xgrid (array_like): Array of grid points from x_min to x_max with spacing of dx
    """
    dx = (x_max - x_min) / self.n_basis
    x_grid = np.arange(-self.n_basis / 2, self.n_basis / 2) * dx
    self.dx = dx
    self.x_grid = x_grid
    return

custom_grid_state_initialization(function_name, **kwargs)

Function to allow for customized grid state initialization.

Parameters:

Name Type Description Default
function_name Callable

name of user-defined function that returns the initial state. Must return an array

 required
Source code in src/qflux/closed_systems/classical_methods.py 
326
327
328
329
330
331
332
333
334
335
336
def custom_grid_state_initialization(self, function_name: Callable, **kwargs):
    """
    Function to allow for customized grid state initialization.

    Args:
        function_name (Callable): name of user-defined function that returns
            the initial state. Must return an array
    """

    self.psio_grid = function_name(**kwargs)
    return

custom_ladder_state_initialization(function_name, **kwargs)

Function to allow for customized ladder state initialization.

Parameters:

Name Type Description Default
function_name Callable

name of user-defined function that returns the initial state. Must return a qutip.Qobj.

 required
Source code in src/qflux/closed_systems/classical_methods.py 
338
339
340
341
342
343
344
345
346
347
348
def custom_ladder_state_initialization(self, function_name: Callable, **kwargs):
    """
    Function to allow for customized ladder state initialization.

    Args:
        function_name (Callable): name of user-defined function that returns
            the initial state. Must return a qutip.Qobj.
    """

    self.psio_op = function_name(**kwargs)
    return

initialize_operators()

Function to initialize core operators in the chosen basis.

Source code in src/qflux/closed_systems/classical_methods.py 
145
146
147
148
149
150
151
152
153
154
def initialize_operators(self):
    """
        Function to initialize core operators in the chosen basis.

    """

    self.a_op = qt.destroy(self.n_basis)
    self.x_op = qt.position(self.n_basis)
    self.p_op = qt.momentum(self.n_basis)
    return

propagate_SOFT()

Function used to propagate with the 2nd-Order Trotter Expansion.

\[ e^{- rac{i}{\hbar} H t} pprox e^{- rac{i}{\hbar} V t/2} e^{- rac{i}{\hbar} T t} e^{- rac{i}{\hbar} V t/2} + \mathcal{O}^{3} \]

Returns:

Name Type Description
dynamics_results_grid array - like

array containing the propagated state shape (n_tsteps x self.n_basis)
Source code in src/qflux/closed_systems/classical_methods.py 
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
def propagate_SOFT(self):
    """
    Function used to propagate with the 2nd-Order Trotter Expansion.

    $$
    e^{- \frac{i}{\\hbar} H t} \approx e^{- \frac{i}{\\hbar} V t/2} e^{- \frac{i}{\\hbar} T t} e^{- \frac{i}{\\hbar} V t/2} + \\mathcal{O}^{3}
    $$

    Returns:
        dynamics_results_grid (array-like): array containing the propagated state
                                            shape (n_tsteps x self.n_basis)

    """
    self.tau = self.tlist[1] - self.tlist[0]
    PE_prop = np.exp(-1.0j * self._PE_grid / 2 * self.tau / self.hbar)
    KE_prop = np.exp(-1.0j * self._KE_grid * self.tau / self.hbar)

    self.PE_prop_grid = PE_prop
    self.KE_prop_grid = KE_prop

    propagated_states = [self.psio_grid]
    psi_t = self.psio_grid
    for ii in range(1, len(self.tlist)):
        psi_t_position_grid = PE_prop * psi_t
        psi_t_momentum_grid = KE_prop * np.fft.fft(psi_t_position_grid, norm="ortho")
        psi_t = PE_prop * np.fft.ifft(psi_t_momentum_grid, norm="ortho")
        propagated_states.append(psi_t)

    self.dynamics_results_grid = np.asarray(propagated_states)
    return

propagate_qt(solver_options=None)

Function used to propagate with qutip.

Parameters:

Name Type Description Default
solver_options dict

A dictionary of arguments to pass to the qutip.sesolve function

 None

Returns:

Name Type Description
dynamics_results array - like

array containing the propagated state
Source code in src/qflux/closed_systems/classical_methods.py 
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
def propagate_qt(self, solver_options : dict = None):
    """
    Function used to propagate with qutip.

    Args:
        solver_options (dict): A dictionary of arguments to pass to the qutip.sesolve function

    Returns:
        dynamics_results (array-like): array containing the propagated state

    """

    options = {'nsteps': len(self.tlist),
                'progress_bar': True}

    if solver_options:
        for key in solver_options:
            options[key] = solver_options[key]

    results = qt.sesolve(self.H_op, self.psio_op, self.tlist,
                         options=options)

    self.dynamics_results_op = results.states
    return

set_H_grid_with_custom_potential(custom_function, **kwargs)

Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

Parameters:

Name Type Description Default
custom_function Callable

Function that defines the custom potential energy. Must return an array

 required
Source code in src/qflux/closed_systems/classical_methods.py 
268
269
270
271
272
273
274
275
276
277
278
279
280
281
def set_H_grid_with_custom_potential(self, custom_function: Callable, **kwargs):
    """
    Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

    Args:
        custom_function (Callable): Function that defines the custom potential
            energy. Must return an array

    """
    potential = custom_function(**kwargs)
    self._PE_grid = potential
    self._KE_grid = self.p_grid ** 2 / (2. * self.mass)
    self.H_grid = self._PE_grid + self._KE_grid
    return

set_H_op_with_custom_potential(custom_function, **kwargs)

Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

Parameters:

Name Type Description Default
custom_function Callable

Function that defines the potential energy in terms of qutip QObj operators. Must return a qutip.Qobj

 required
Source code in src/qflux/closed_systems/classical_methods.py 
284
285
286
287
288
289
290
291
292
293
294
295
296
def set_H_op_with_custom_potential(self, custom_function: Callable, **kwargs):
    """
    Function to allow for user-defined potential defined by custom_function. Must be a function of qutip operators.

    Args:
        custom_function (Callable): Function that defines the potential
            energy in terms of qutip QObj operators. Must return a qutip.Qobj
    """
    potential = custom_function(**kwargs)
    self._PE_op = potential
    self._KE_op = self.p_op ** 2 / (2. * self.mass)
    self.H_op = self._PE_op + self._KE_op
    return

set_coordinate_operators(x_min=-7.0, x_max=7.0, reorder_p=True)

Populate the self.x_grid, self.p_grid, self.dx, and self.dp attributes. This functions generates an array of self.n_basis evenly spaced values.

Parameters:

Name Type Description Default
x_min

float Minimum value of x-coordinates

 -7.0
x_max

float Maximum value of x-coordinates

 7.0
reorder_p

bool Boolean flag to determine whether momentum values should be reordered to be compatible with the FFT routine or not.

 True

Returns:

Type Description
None

self.dx : float Spacing between points in the x-coordinate grid.
None

self.xgrid : array_like Array of x-values
None

self.dp : float Spacing between points in the p-coordinate grid.
None

self.pgrid : array_like Array of p-values
Source code in src/qflux/closed_systems/classical_methods.py 
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
def set_coordinate_operators(self, x_min: float = -7., x_max: float = 7., reorder_p: bool = True) -> None:
    """
    Populate the `self.x_grid`, `self.p_grid`, `self.dx`, and `self.dp`
    attributes. This functions generates an array of `self.n_basis`
    evenly spaced values.

    Args:
        x_min : float
            Minimum value of x-coordinates
        x_max : float
            Maximum value of x-coordinates
        reorder_p : bool
            Boolean flag to determine whether momentum values should be
            reordered to be compatible with the FFT routine or not.

    Returns:
        self.dx : float
            Spacing between points in the x-coordinate grid.
        self.xgrid : array_like
            Array of x-values
        self.dp : float
            Spacing between points in the p-coordinate grid.
        self.pgrid : array_like
            Array of p-values
    """
    self._get_xgrid(x_min, x_max)
    self._get_pgrid(x_min, x_max, reorder=reorder_p)
    return

set_hamiltonian(potential_type='harmonic', **kwargs)

Function to define Hamiltonian.

Parameters:

Name Type Description Default
potential_type

str String defining the type of potential energy surface. Available options are: ('harmonic', 'quartic', ...)

Note: You can manually define your potential energy using the functions: - set_H_grid_with_custom_potential - set_H_op_with_custom_potential

 'harmonic'
Source code in src/qflux/closed_systems/classical_methods.py 
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
def set_hamiltonian(self, potential_type: str = 'harmonic', **kwargs):
    """
    Function to define Hamiltonian.

    Args:
        potential_type : str
            String defining the type of potential energy surface.
            Available options are: ('harmonic', 'quartic', ...)

            Note: You can manually define your potential energy using the functions:
                - set_H_grid_with_custom_potential
                - set_H_op_with_custom_potential

    """

    if potential_type == 'harmonic':
        self._set_hamiltonian_grid(potential_type=potential_type, **kwargs)
        self._set_hamiltonian_qt(potential_type=potential_type, **kwargs)
    elif potential_type == 'quartic':
        self._set_hamiltonian_grid(potential_type=potential_type, **kwargs)
        self._set_hamiltonian_qt(potential_type=potential_type, **kwargs)
    else:
        print('Error, this potential type has not yet been implemented!')
        print('Set your parameters with the custom functions!')
    return

set_initial_state(wfn_omega=1.0)

Function to define the initial state. By default, a coherent state is used as the initial state defined in the basis chosen upon instantiation

Parameters:

Name Type Description Default
wfn_omega float

Defines the frequency/width of the initial state. Default is 1.0 au.

 1.0
Source code in src/qflux/closed_systems/classical_methods.py 
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
def set_initial_state(self, wfn_omega: float = 1.0):
    """
    Function to define the initial state. By default, a coherent state is
    used as the initial state defined in the basis chosen upon instantiation

    Args:
        wfn_omega (float, optional): Defines the frequency/width of the initial state.
            Default is 1.0 au.
    """

    alpha_val = (self.xo + 1j * self.po) / np.sqrt(2)
    psio = qt.coherent(self.n_basis, alpha=alpha_val)
    # Now populate the initial state in the grid basis
    normalization = (self.mass * wfn_omega / np.pi / self.hbar) ** (0.25)
    exponential = np.exp(-1 * (self.mass * wfn_omega / self.hbar / 2) *
                         ((self.x_grid - self.xo) ** 2)
                         + 1j * self.po * self.x_grid / self.hbar
                         )

    coherent_state = normalization * exponential
    # Set the attributes
    self.psio_grid = coherent_state
    self.psio_op = psio
    return

set_propagation_time(total_time, n_tsteps)

Function to define the propagation time, an array of times from t=0 to total_time, with n_tsteps equally-spaced steps.

total_time : float The total time for which we wish to compute the dynamics. n_tsteps : int The number of equally-spaced time steps used to compute the dynamics

Returns: self.tlist : array-like

Source code in src/qflux/closed_systems/classical_methods.py 
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
def set_propagation_time(self, total_time: float, n_tsteps: int):
    """
    Function to define the propagation time, an array of times from
    t=0 to total_time, with n_tsteps equally-spaced steps.

    Args:
    total_time : float
        The total time for which we wish to compute the dynamics.
    n_tsteps : int
        The number of equally-spaced time steps used to compute the dynamics

    Returns:
    self.tlist : array-like

    """

    self.tlist = np.linspace(0., total_time, n_tsteps+1)
    self.dt = self.tlist[1] - self.tlist[0]
    self.n_tsteps = n_tsteps
    return

QubitDynamicsCS

Bases: DynamicsCS

Class to extend DynamicsCS by adding qubit-based methods for dynamics.

Source code in src/qflux/closed_systems/qubit_methods.py 
 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
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
class QubitDynamicsCS(DynamicsCS):
    """
    Class to extend `DynamicsCS` by adding qubit-based methods for dynamics.

    """
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.n_qubits        = int(np.log2(self.n_basis))
        self.quantum_circuit = None


    def _create_QSOFT_Circuit(self, psio: npt.ArrayLike=None):
        """
        Function to construct the QSOFT Circuit.

        Args:
            psio (npt.ArrayLike): initial state that we wish to propagate
        """
        tgrid = self.tlist
        time_step = self.dt
        n_qubits = self.n_qubits
        # Qubit-Basis Propagators
        self.prop_PE_qubit = np.diag(np.exp(-1j*self._PE_grid/2*time_step))
        self.prop_KE_qubit = np.diag(np.exp(-1j*self._KE_grid*time_step))

        q_reg = QuantumRegister(n_qubits)
        c_reg = ClassicalRegister(n_qubits)
        qc = QuantumCircuit(q_reg)
        if type(psio) == type(None):
            qc.initialize(self._psio_grid, q_reg[:], normalize=True)
        else:
            qc.initialize(psio, q_reg[:], normalize=True)
        # Define our PE and KE propagators in Qiskit-friendly manner
        PE_cirq_op = Operator(self.prop_PE_qubit)
        KE_cirq_op = Operator(self.prop_KE_qubit)
        qc.append(PE_cirq_op, q_reg)
        qc.append(QFT(self.n_qubits, do_swaps=True, inverse=False), q_reg)
        qc.append(KE_cirq_op, q_reg)
        qc.append(QFT(self.n_qubits, do_swaps=True, inverse=True), q_reg)
        qc.append(PE_cirq_op, q_reg)
        self.quantum_circuit = qc
        return(qc)


    def _execute_circuit(self, QCircuit: QuantumCircuit, backend=None, shots: int = None, real_backend: bool = False):
        """
            Function to replace the now-deprecated Qiskit
            `QuantumCircuit.execute()` method.

            Args:
                QCircuit (qiskit.QuantumCircuit): qiskit.QuantumCircuit object
                backend (qiskit.Backend): qiskit backend instance
                shots (int): the number of shots to use for circuit sampling

            Returns:
                job: an executed quantum circuit job
        """
        if shots:
            n_shots = shots
        else:
            n_shots = 1024 # Use the qiskit default if not specified
        backend_type = type(backend)
        sv_type = qiskit_aer.backends.statevector_simulator.StatevectorSimulator
        if backend_type == sv_type:
            real_backend = False
        else:
            real_backend = True

        if real_backend:
            QCircuit.measure_all()
            qc = transpile(QCircuit, backend=backend)
            sampler = Sampler(backend)
            job = sampler.run([qc], shots=n_shots)
        else:
            # Transpile circuit with statevector backend
            tmp_circuit = transpile(QCircuit, backend)
            # Run the transpiled circuit
            job = backend.run(tmp_circuit, n_shots=shots)
        return(job)


    def propagate_qSOFT(self, backend=None, n_shots: int = 1024):
        """Function to propagate dynamics object with the qubit SOFT method.

            Args:
                backend (qiskit.Backend): qiskit backend object
                n_shots (int): specifies the number of shots to use when
                    executing the circuit

            Example for using the Statevector Simulator backend:
                >>> from qiskit_aer import Aer
                >>> backend = Aer.get_backend('statevector_simulator')
                >>> self.propagate_qSOFT(backend=backend)
        """
        if backend is None:
            print('A valid backend must be provided ')
        backend_type = type(backend)
        sv_type = qiskit_aer.backends.statevector_simulator.StatevectorSimulator
        if backend_type != sv_type:
            self._propagate_qSOFT_real(backend=backend, n_shots=n_shots)
            return
        else:

            psi_in = self.psio_grid
            # Get initial state from qiskit routine
            q_reg = QuantumRegister(self.n_qubits)
            c_reg = ClassicalRegister(self.n_qubits)
            qc = QuantumCircuit(q_reg, c_reg)
            qc.initialize(self.psio_grid, q_reg[:], normalize=True)
            qc_result = self._execute_circuit(qc, backend=backend, shots=n_shots)
            psio_cirq = qc_result.result().get_statevector().data
            psi_in = psio_cirq
            # Now do propagation loop
            qubit_dynamics_results = [psio_cirq]
            for ii in trange(1, len(self.tlist)):
                circuit = self._create_QSOFT_Circuit(psio=psi_in)
                executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
                psi_out = executed_circuit.result().get_statevector().data
                qubit_dynamics_results.append(psi_out)
                psi_in = psi_out

            self.dynamics_results_qSOFT = np.asarray(qubit_dynamics_results)
            return


    def get_statevector_from_counts(self, counts, n_shots):
        new_statevector = np.zeros_like(self.psio_grid)

        for key in counts:
            little_endian_int = int(key, 2)
            new_statevector[little_endian_int] = counts[key]/n_shots
        return(new_statevector)


    def _propagate_qSOFT_real(self, backend='statevector_simulator', n_shots=1024):
        """
            Function to propagate dynamics object with the qubit SOFT method.

            Args:
                backend (qiskit.Backend): qiskit backend object
                n_shots (int): specifies the number of shots to use when
                    executing the circuit

            Example for using the Statevector Simulator backend:
                >>> from qiskit_aer import Aer
                >>> backend = Aer.get_backend('statevector_simulator')
                >>> self.propagate_qSOFT(backend=backend)
        """


        psi_in = self.psio_grid
        # Get initial state from qiskit routine
        q_reg = QuantumRegister(self.n_qubits)
        c_reg = ClassicalRegister(self.n_qubits, name='c')
        qc = QuantumCircuit(q_reg, c_reg)
        qc.initialize(self.psio_grid, q_reg[:], normalize=True)
        # Now do propagation loop
        qubit_dynamics_results = []
        for ii in trange(len(self.tlist)):
            circuit = self._create_QSOFT_Circuit(psio=psi_in)
            executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
            circuit_result = executed_circuit.result()
            measured_psi = circuit_result[0].data['meas'].get_counts()
            self._last_measurement = measured_psi
            psi_out = self.get_statevector_from_counts(measured_psi, n_shots)
            psi_in = psi_out
            qubit_dynamics_results.append(psi_out)
            psi_in = psi_out

        self.dynamics_results_qSOFT = np.asarray(qubit_dynamics_results)
        return


    def _construct_pauli_gate(self, hamiltonian_matrix=None):
        """
            Function to construct a pauli evolution gate from Hamiltonian

            Args:
                hamiltonian_matrix (npt.ArrayLike): array-like matrix representing the hamiltonian of interest
                    If not provided, use the operator representation of the Hamiltonian by default.

        """

        if type(hamiltonian_matrix) == type(None):
            decomposed_H = decompose(self.H_op.full())
        else:
            decomposed_H = decompose(hamiltonian_matrix)
        H_pauli_sum  = pauli_strings_2_pauli_sum(decomposed_H)
        synthesizer  = LieTrotter(reps=2)
        prop_pauli_H = PauliEvolutionGate(operator=H_pauli_sum, time=self.dt, synthesis=synthesizer)
        self.pauli_prop = prop_pauli_H
        self._pauli_hamiltonian = decomposed_H
        return


    def _construct_pauli_cirq(self, psio=None):
        q_reg = QuantumRegister(self.n_qubits)
        c_reg = ClassicalRegister(self.n_qubits)
        qc = QuantumCircuit(q_reg, c_reg)
        qc.initialize(psio, q_reg[:], normalize=True)
        qc.append(self.pauli_prop, q_reg)
        self.quantum_circuit = qc
        return(qc)


    def propagate_qmatvec(self, backend=None, n_shots: int = 1024, hamiltonian_matrix=None, initial_state=None):
        r"""
            Function to propagate dynamics object with the qubit matvec method.

            Args:
                backend (qiskit.Backend): qiskit backend object
                n_shots (int): specifies the number of shots to use when
                    executing the circuit
                hamiltonian_matrix (npt.ArrayLike): array-like matrix representing the Hamiltonian
                    Used to construct the propagator:

                    $$ U(t) = e^{- i H t / \hbar} $$ 

                    By default, the operator representation of the hamiltonian `self.H_op` is used.
                initial_state (npt.ArrayLike): array-like vector representing the initial state
            Example for using the Statevector Simulator backend:
                >>> from qiskit_aer import Aer
                >>> backend = Aer.get_backend('statevector_simulator')
                >>> self.propagate_qSOFT(backend=backend)
        """

        # Create the Pauli propagator:
        self._construct_pauli_gate(hamiltonian_matrix=hamiltonian_matrix)

        q_reg = QuantumRegister(self.n_qubits)
        c_reg = ClassicalRegister(self.n_qubits)
        qc = QuantumCircuit(q_reg, c_reg)
        # Initialize State
        if type(initial_state) == type(None):
            qc.initialize(self.psio_op.full().flatten(), q_reg[:], normalize=True)
        else:
            qc.initialize(initial_state, q_reg[:], normalize=True)
        qc_result = self._execute_circuit(qc, backend=backend, shots=n_shots)
        psio_cirq = qc_result.result().get_statevector().data
        psi_in = psio_cirq
        new_qubit_dynamics_result = [psio_cirq]
        for ii in trange(1, len(self.tlist)):
            circuit = self._construct_pauli_cirq(psio=psi_in)
            executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
            psi_out = executed_circuit.result().get_statevector().data
            new_qubit_dynamics_result.append(psi_out)
            psi_in = psi_out
        self.dynamics_results_qubit = np.asarray(new_qubit_dynamics_result)
        return

_construct_pauli_gate(hamiltonian_matrix=None)

Function to construct a pauli evolution gate from Hamiltonian

Parameters:

Name Type Description Default
hamiltonian_matrix ArrayLike

array-like matrix representing the hamiltonian of interest If not provided, use the operator representation of the Hamiltonian by default.

 None
Source code in src/qflux/closed_systems/qubit_methods.py 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
def _construct_pauli_gate(self, hamiltonian_matrix=None):
    """
        Function to construct a pauli evolution gate from Hamiltonian

        Args:
            hamiltonian_matrix (npt.ArrayLike): array-like matrix representing the hamiltonian of interest
                If not provided, use the operator representation of the Hamiltonian by default.

    """

    if type(hamiltonian_matrix) == type(None):
        decomposed_H = decompose(self.H_op.full())
    else:
        decomposed_H = decompose(hamiltonian_matrix)
    H_pauli_sum  = pauli_strings_2_pauli_sum(decomposed_H)
    synthesizer  = LieTrotter(reps=2)
    prop_pauli_H = PauliEvolutionGate(operator=H_pauli_sum, time=self.dt, synthesis=synthesizer)
    self.pauli_prop = prop_pauli_H
    self._pauli_hamiltonian = decomposed_H
    return

_create_QSOFT_Circuit(psio=None)

Function to construct the QSOFT Circuit.

Parameters:

Name Type Description Default
psio ArrayLike

initial state that we wish to propagate

 None
Source code in src/qflux/closed_systems/qubit_methods.py 
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
def _create_QSOFT_Circuit(self, psio: npt.ArrayLike=None):
    """
    Function to construct the QSOFT Circuit.

    Args:
        psio (npt.ArrayLike): initial state that we wish to propagate
    """
    tgrid = self.tlist
    time_step = self.dt
    n_qubits = self.n_qubits
    # Qubit-Basis Propagators
    self.prop_PE_qubit = np.diag(np.exp(-1j*self._PE_grid/2*time_step))
    self.prop_KE_qubit = np.diag(np.exp(-1j*self._KE_grid*time_step))

    q_reg = QuantumRegister(n_qubits)
    c_reg = ClassicalRegister(n_qubits)
    qc = QuantumCircuit(q_reg)
    if type(psio) == type(None):
        qc.initialize(self._psio_grid, q_reg[:], normalize=True)
    else:
        qc.initialize(psio, q_reg[:], normalize=True)
    # Define our PE and KE propagators in Qiskit-friendly manner
    PE_cirq_op = Operator(self.prop_PE_qubit)
    KE_cirq_op = Operator(self.prop_KE_qubit)
    qc.append(PE_cirq_op, q_reg)
    qc.append(QFT(self.n_qubits, do_swaps=True, inverse=False), q_reg)
    qc.append(KE_cirq_op, q_reg)
    qc.append(QFT(self.n_qubits, do_swaps=True, inverse=True), q_reg)
    qc.append(PE_cirq_op, q_reg)
    self.quantum_circuit = qc
    return(qc)

_execute_circuit(QCircuit, backend=None, shots=None, real_backend=False)

Function to replace the now-deprecated Qiskit QuantumCircuit.execute() method.

Parameters:

Name Type Description Default
QCircuit QuantumCircuit

qiskit.QuantumCircuit object

 required
backend Backend

qiskit backend instance

 None
shots int

the number of shots to use for circuit sampling

 None

Returns:

Name Type Description
job

an executed quantum circuit job
Source code in src/qflux/closed_systems/qubit_methods.py 
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
def _execute_circuit(self, QCircuit: QuantumCircuit, backend=None, shots: int = None, real_backend: bool = False):
    """
        Function to replace the now-deprecated Qiskit
        `QuantumCircuit.execute()` method.

        Args:
            QCircuit (qiskit.QuantumCircuit): qiskit.QuantumCircuit object
            backend (qiskit.Backend): qiskit backend instance
            shots (int): the number of shots to use for circuit sampling

        Returns:
            job: an executed quantum circuit job
    """
    if shots:
        n_shots = shots
    else:
        n_shots = 1024 # Use the qiskit default if not specified
    backend_type = type(backend)
    sv_type = qiskit_aer.backends.statevector_simulator.StatevectorSimulator
    if backend_type == sv_type:
        real_backend = False
    else:
        real_backend = True

    if real_backend:
        QCircuit.measure_all()
        qc = transpile(QCircuit, backend=backend)
        sampler = Sampler(backend)
        job = sampler.run([qc], shots=n_shots)
    else:
        # Transpile circuit with statevector backend
        tmp_circuit = transpile(QCircuit, backend)
        # Run the transpiled circuit
        job = backend.run(tmp_circuit, n_shots=shots)
    return(job)

_propagate_qSOFT_real(backend='statevector_simulator', n_shots=1024)

Function to propagate dynamics object with the qubit SOFT method.

Parameters:

Name Type Description Default
backend Backend

qiskit backend object

 'statevector_simulator'
n_shots int

specifies the number of shots to use when executing the circuit

 1024
Example for using the Statevector Simulator backend

from qiskit_aer import Aer backend = Aer.get_backend('statevector_simulator') self.propagate_qSOFT(backend=backend)

Source code in src/qflux/closed_systems/qubit_methods.py 
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
def _propagate_qSOFT_real(self, backend='statevector_simulator', n_shots=1024):
    """
        Function to propagate dynamics object with the qubit SOFT method.

        Args:
            backend (qiskit.Backend): qiskit backend object
            n_shots (int): specifies the number of shots to use when
                executing the circuit

        Example for using the Statevector Simulator backend:
            >>> from qiskit_aer import Aer
            >>> backend = Aer.get_backend('statevector_simulator')
            >>> self.propagate_qSOFT(backend=backend)
    """


    psi_in = self.psio_grid
    # Get initial state from qiskit routine
    q_reg = QuantumRegister(self.n_qubits)
    c_reg = ClassicalRegister(self.n_qubits, name='c')
    qc = QuantumCircuit(q_reg, c_reg)
    qc.initialize(self.psio_grid, q_reg[:], normalize=True)
    # Now do propagation loop
    qubit_dynamics_results = []
    for ii in trange(len(self.tlist)):
        circuit = self._create_QSOFT_Circuit(psio=psi_in)
        executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
        circuit_result = executed_circuit.result()
        measured_psi = circuit_result[0].data['meas'].get_counts()
        self._last_measurement = measured_psi
        psi_out = self.get_statevector_from_counts(measured_psi, n_shots)
        psi_in = psi_out
        qubit_dynamics_results.append(psi_out)
        psi_in = psi_out

    self.dynamics_results_qSOFT = np.asarray(qubit_dynamics_results)
    return

propagate_qSOFT(backend=None, n_shots=1024)

Function to propagate dynamics object with the qubit SOFT method.

Parameters:

Name Type Description Default
backend Backend

qiskit backend object

 None
n_shots int

specifies the number of shots to use when executing the circuit

 1024
Example for using the Statevector Simulator backend

from qiskit_aer import Aer backend = Aer.get_backend('statevector_simulator') self.propagate_qSOFT(backend=backend)

Source code in src/qflux/closed_systems/qubit_methods.py 
 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
def propagate_qSOFT(self, backend=None, n_shots: int = 1024):
    """Function to propagate dynamics object with the qubit SOFT method.

        Args:
            backend (qiskit.Backend): qiskit backend object
            n_shots (int): specifies the number of shots to use when
                executing the circuit

        Example for using the Statevector Simulator backend:
            >>> from qiskit_aer import Aer
            >>> backend = Aer.get_backend('statevector_simulator')
            >>> self.propagate_qSOFT(backend=backend)
    """
    if backend is None:
        print('A valid backend must be provided ')
    backend_type = type(backend)
    sv_type = qiskit_aer.backends.statevector_simulator.StatevectorSimulator
    if backend_type != sv_type:
        self._propagate_qSOFT_real(backend=backend, n_shots=n_shots)
        return
    else:

        psi_in = self.psio_grid
        # Get initial state from qiskit routine
        q_reg = QuantumRegister(self.n_qubits)
        c_reg = ClassicalRegister(self.n_qubits)
        qc = QuantumCircuit(q_reg, c_reg)
        qc.initialize(self.psio_grid, q_reg[:], normalize=True)
        qc_result = self._execute_circuit(qc, backend=backend, shots=n_shots)
        psio_cirq = qc_result.result().get_statevector().data
        psi_in = psio_cirq
        # Now do propagation loop
        qubit_dynamics_results = [psio_cirq]
        for ii in trange(1, len(self.tlist)):
            circuit = self._create_QSOFT_Circuit(psio=psi_in)
            executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
            psi_out = executed_circuit.result().get_statevector().data
            qubit_dynamics_results.append(psi_out)
            psi_in = psi_out

        self.dynamics_results_qSOFT = np.asarray(qubit_dynamics_results)
        return

propagate_qmatvec(backend=None, n_shots=1024, hamiltonian_matrix=None, initial_state=None)

Function to propagate dynamics object with the qubit matvec method.

Parameters:

Name Type Description Default
backend Backend

qiskit backend object

 None
n_shots int

specifies the number of shots to use when executing the circuit

 1024
hamiltonian_matrix ArrayLike

array-like matrix representing the Hamiltonian Used to construct the propagator:

\[ U(t) = e^{- i H t / \hbar} \]

By default, the operator representation of the hamiltonian self.H_op is used.

 None
initial_state ArrayLike

array-like vector representing the initial state

 None

Example for using the Statevector Simulator backend: >>> from qiskit_aer import Aer >>> backend = Aer.get_backend('statevector_simulator') >>> self.propagate_qSOFT(backend=backend)

Source code in src/qflux/closed_systems/qubit_methods.py 
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
def propagate_qmatvec(self, backend=None, n_shots: int = 1024, hamiltonian_matrix=None, initial_state=None):
    r"""
        Function to propagate dynamics object with the qubit matvec method.

        Args:
            backend (qiskit.Backend): qiskit backend object
            n_shots (int): specifies the number of shots to use when
                executing the circuit
            hamiltonian_matrix (npt.ArrayLike): array-like matrix representing the Hamiltonian
                Used to construct the propagator:

                $$ U(t) = e^{- i H t / \hbar} $$ 

                By default, the operator representation of the hamiltonian `self.H_op` is used.
            initial_state (npt.ArrayLike): array-like vector representing the initial state
        Example for using the Statevector Simulator backend:
            >>> from qiskit_aer import Aer
            >>> backend = Aer.get_backend('statevector_simulator')
            >>> self.propagate_qSOFT(backend=backend)
    """

    # Create the Pauli propagator:
    self._construct_pauli_gate(hamiltonian_matrix=hamiltonian_matrix)

    q_reg = QuantumRegister(self.n_qubits)
    c_reg = ClassicalRegister(self.n_qubits)
    qc = QuantumCircuit(q_reg, c_reg)
    # Initialize State
    if type(initial_state) == type(None):
        qc.initialize(self.psio_op.full().flatten(), q_reg[:], normalize=True)
    else:
        qc.initialize(initial_state, q_reg[:], normalize=True)
    qc_result = self._execute_circuit(qc, backend=backend, shots=n_shots)
    psio_cirq = qc_result.result().get_statevector().data
    psi_in = psio_cirq
    new_qubit_dynamics_result = [psio_cirq]
    for ii in trange(1, len(self.tlist)):
        circuit = self._construct_pauli_cirq(psio=psi_in)
        executed_circuit = self._execute_circuit(circuit, backend=backend, shots=n_shots)
        psi_out = executed_circuit.result().get_statevector().data
        new_qubit_dynamics_result.append(psi_out)
        psi_in = psi_out
    self.dynamics_results_qubit = np.asarray(new_qubit_dynamics_result)
    return