6. Modules

6.1. hamitonians

This module holds the class that holds the Hamiltonian.

class hamiltonians.hamiltonian(H0, V, A=None, g=None, gradg=None)[source]

A class to hold and manipulate the Hamiltonian and dissipation terms.

The Hamiltonian should contain a static part, a list of possible perturbations, and also a list of dissipative terms for the Lindblad equation (which is not, of course, a proper term of the Hamiltonian)

There are two options for the user to pass the Hamiltonian to qocttools: as a list of Qobj objects, or as a Qobj-valued user-defined function.

In the first case, the Hamiltonian has the form:

\[H(t) = H_0 + \sum_{i=1}^N f_i(t) V_i\]

The pulses \(f_i(t)\) are not described in the hamiltonian object. The user must then pass H0 and a list of perturbations V (it may also be just one). In addition, for open systems, the hamiltonian object should also contain the Lindblad operators, which should be supplied by the user as a list of Qobj objects.
In the second case, the user should supply Qobj-valued functions. In this case, the H0 function contains the full Hamiltonian, and depends on time and on the set of “pulses” \(f_i(t)\):

\[H_0 = H_0(t, f_1(t), \dots, f_N(t))\]

The user must then also supply, in the argument V, a list with the derivatives of \(H0\) with respect to each \(f_i\) argument:

\[V_i(t, f_1(t), \dots, f_N(t)) = \frac{\partial H_o}{\partial f_i}(t, f_1(t), \dots, f_N(t))\]

The function H0 (and the functions in the V list) should have the following inteface:
```
def H0(t, args):
    f = args["f"]
    # f[0], f[1], f[2], ..., f[N-1] are the time-dependent functions
```
The arguments are a float t that means time, and a dictionary args with just one element with name “f”, and whose value is a list of functions of time.

Parameters:

H0 (Qobj or Qobj-valued function) – The qutip Qobj instance holding the static part of the Hamiltonian.
V (Qobj or list of Qobj or Qobj-valued function or list of Qobj-valued functions.) – One or various Qobj instances holding external perturbations.
A (Qobj or list of Qobj, default = None) – One or various Qobj instances holding the jump operators defining the Lindbladian.
g (still undocumnted, default = None) –
gradg (still undocumented, default = None) –

H0

The qutip Qobj instance holding the static part of the Hamiltonian.

Type:: Qobj

V

One or various Qobj instances holding external perturbations.

Type:: list of Qobj

A

One or various Qobj instances holding the jump operators defining the Lindblad operator.

Type:: list of Qobj

has_dissipative_terms()[source]

Returns True if the dissipative terms are not None

Returns:: True if the attribute holding the dissipative terms is not None, False otherwise.
Return type:: bool

hamiltonians.toliouville(h)[source]

This function transforms a Hamiltonian, defined in Hilbert space, together with the set of Lindblad operators, to the corresponding Liouvillian in Liouville space. It is also returned as a hamiltonian object

Parameters:: h (hamiltonian) – The input Hamiltonian, to be transformed.
Returns:: A hamiltonian object with the transformed Hamiltonian.
Return type:: hamiltonian

6.2. qoct

class qoct.Qoct(H, T, ntsteps, tg, f, y0, interaction_picture=False, solve_method='cfmagnus4', equality_constraints=None, sparse=False, output_file=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>, new_parametrization=False, mpi_comm=None, floquet_mode='qoct')[source]

The class that is used to setup the optimization.

It stores the Hamiltonian, the definition of the target, the pulses, the initial state, and some parameters that determine how the optimization is done.

Parameters:

H (hamiltonian or list of hamiltonian) – A hamiltonian class instance, or list of instances for multitarget optimizations.
T (float) – total propagation time
ntsteps (int) – number of time steps
tg (Target) – The object that defines the target functional.
f (pulse or list of pulse) – The pulse object, or list of pulse objects (in the case the Hamiltonian(s) have more than one perturbations)
y0 (Qobj or list of Qobj) – The initial state(s) for the optimization. It should be a list for multitarget optimizations. It can be a Hilbert space state, a density matrix, or a propagator.
interaction_picture (bool, default = False) – The calculations can be done in the interaction picture, or in Schrödinger’s picture.
solve_method (str, default = 'cfmagnus4') – The method to use for the propagations.
equality_constraints (list, default = None) – The optimization can be done with constraints, that should be specified as a list of [func, tol] objects: func is a function that defines the constraint, and tol is a float that sets the tolerance accepted for the fulfillment of that constraint.
sparse (bool, default = False) – It forces the use of sparse algebra only (otherwise, some parts of the calculation are done with dense matrices).
output_file (file, default = sys.stdout) – One may ask the internal procedures of the qoct object to output messages to a file different from standard output

y0

Type:: list of Qobj

H

Type:: list of hamiltonian

almost_converged(fraction=0.99)[source]

Returns the iteration number at which the convergence was almost achieved

It analyzes the convergence history of an optimization, and looks for the iteration at which the value of the functional was already a fraction (0.99 by default) of the maximum.

Parameters:

fraction (float, default = 0.99) – The fraction of the maximum value at which the process is considered to be almost converged. If, at iteration 76, the value of the functional is already fraction * maximum, then this function will return 76 (along with the number of propagations necessary to reach iteration 76.

Returns:

float – The iteration at which the process was almost converged
float – The number of propagations done until that iteration.

check_grad(u, m=None, delta=0.01)[source]

A check on the accuracy of the gradient calculation

It computes one component of the gradient of the target functional using both the QOCT expression, and a finite-difference algorithm (the Ridders algorithm).

Parameters:

u (ndarray) – A numpy array holding all the control parameters
m (int, default = None) – Which gradient component to compute. If None, it will use the one with the largest absolute value
delta (float, default = 0.01) – A parameter determining the starting finite difference in the Ridders algorithm.

Returns:

float – The gradient component as computed with the QOCT formula.
float – The gradient component as computed with the finite-difference formula.
float – An estimation of the gradient accuracy computed with the finite difference formula.

dgfunc(u)[source]

Returns the value and gradient of the target functional for a set of parameters u.

It just calls the dgfunction function, using all the information contained in the class object.

Parameters:

u (ndarray) – A numpy array holding all the control parameters

Returns:

float – The value of the target parameter.
ndarray – The gradient of the target functional.

gfunc(u)[source]

Returns the value of the target functional for a set of parameters u.

It just calls the gfunction function, using all the information contained in the class object.

Parameters:: u (ndarray) – A numpy array holding all the control parameters
Returns:: The value of the target parameter.
Return type:: float

maximize(maxeval=100, stopval=None, verbose=False, algorithm=11, local_algorithm=None, upper_bounds=None, lower_bounds=None, tolerance=1e-06)[source]

Performs the QOCT maximization

This is the main procedure of the class, as it is the one that launches the optimization process. There are a number of parameters controlling the way in which this optimization is done, that can be controlled here.

The starting guess parameters for the optimizatio are read from the values of the pulses, that were associated to the class object when it was created.

Parameters:

maxeval (int, default = 100) – The maximum number of function evaluations that the optimization algorithm is allowed to do. Note that, depending on the algorithm that is used, this may not be strictly enforced.
stopval (float, default = None) – The optimization will stop if the value of the target function reaches a value above stopval.
verbose (bool, default = False) – The code will output more or less data, depending on the value of this parameter
algorithm (int, default = nlopt.LD_LBFGS) –
The function maximization algorith to be used. qocttools relies on the nlopt optimization library:

https://nlopt.readthedocs.io/en/latest/

and therefore one can use any of the algorithms that are defined there. See

https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/

for a list. If you include the nlopt module in your driver Python code, you can use the parameters that are defined there. Although note that in the case of Python, the nomenclature would be nlopt.LD_MMA, nlopt.LN_COBYLA, etc., and not NLOPT_LD_MMA, NLOPT_LN_COBYLA, etc. as it is done when using C.

In addition, if you set algorithm = -1, the code will not use any of the nlopt algorithms, but a fairly basic implementation of Krotov’s algorithm.
local_algorithm (int, default = None) – Some nlopt algorithms use a global algorithm and a local algorithm. If this is the case, you need to specify here the local one (whereas the argument of “algorithm” would be the global one).
upper_bounds (ndarray, default = None) – You may specify upper bound for the parameters – although this is only possible with some of the nlopt algorithms; the code will stop if bounds are requested and the algorithm is not capable. See the documetation of nlopt for guidance.
lower_bounds (ndarray, default = None) – The lower bounds for the parameters.
tolerance (float, default = 1E-06) – The algorithm will consider that a successful optimization has been achieved, if a certain tolerance is obtained, characterized by this parameter. This is a tolerance imposed on the function vale: when the algorithm is not capable of modifying the function value beyond this tolerance, the code will stop and, in principle, a local maximum has been found. Note, however, that this local maximum may not be a good one. If that is the case, one should change the initial guess, the algorithm, or any of the other choices that is done.

Returns:

ndarray – A numpy array with the optimized control parameters.
float – The value of the local maximum that has been achieved
int – An error code; if it is not zero, something went wrong.

6.3. pulses

This module includes the ‘pulse’ class, that is used to contain the time-dependent real function that is used to define the perturbation part of the Hamiltonian.

Also, it contains definitions of ‘typical’ pulse forms, such as the pi-pulse.

pulses.pi_pulse(t, u)[source]

A pi-pulse value at time t, parameters specifiey by the ndarray u

Returns the pi-pulse value at a specific time

Input: A: pi-pulse amplitude w: pi-pulse frecuency t0: time in which start the propagation of the pulse t: time in which the pulse is calculated t_duration: pi-pulse duration
Output:: if the t value in which the pulse calculation is requested is between t0 and t0 + t_duration, return the value of a pi-pulse with frecuency w and amplitude A, otherwise return 0.

pulses.pi_pulse_chain(t, u)[source]: Returns the value of a chain of pi pulse at a specific time

class pulses.pulse(type, T, u, bound=None)[source]

Definition of a class to hold control functions, a.k.a. pulses.

The “pulses” are real time-dependent functions defined in the interval \([0, T]\), defined as functions \(f(u_1, u_2, ..., u_P, T)\), where \(u = u_1, ..., u_P\) are parameters (the control parameters). These pulses are the control functions used by qocttools.

Depending on the type of the pulse, the parametrized form of these pulses is:

fourier

In this case, the pulse is a simple Fourier expansion:

\[f(u, t) = \frac{u_0}{\sqrt{T}} + \frac{2}{\sqrt{T}} \sum_{i=1}^M \left[ u_{2i}\cos(\omega_i t) + u_{2i-1}\sin(\omega_i t)\right]\]

The pulse therfore requires an odd number of parameters, \(P = 2M + 1\). The frequencies are \(\omega_i = i\omega_0 = i\frac{2\pi}{T}\) for \(i=1, \dots, M\). The cutoff M will be decided by the number of parameters that are passed when initializing the object. Therefore, in order to initialize a Fourier pulse with \(M=5\):
```
M = 5
u = np.zeros(2*M+1)
f = pulses.pulse('fourier', T, u)
```
bound_fourier

This is a modification of the normal Fourier expansion, that ensures that the absolute value of the pulse is never larger than a certain bound. The definition is:

\[f(u, t) = \Phi(g(u, t))\]

where \(g(u, t)\) is a normal Fourier expansion as the one given above, and:

\[\Phi(x) = \frac{\kappa x}{\kappa+\vert x\vert}\]

The value of \(\kappa\) must be supplied by the bound argument:
```
M = 5
kappa = 3.0
u = np.zeros(2*M+1)
f = pulses.pulse('fourier', T, u, bound = kappa)
```
realtime

In this case, the pulse is given by a piecewise-constant function. If the time discretization is given by \(t_0 = 0, t_1 = \Delta t, t_2 = 2\Delta t,\dots,t_{N-1}=(N-1)\Delta t = T\), then the function is given by:

\[f(u, t) = \sum_{i=0}^{N-1} u_i {\bf 1}_{[t_i,t_{i+1}]}\]

where \({\bf 1}_{[t_i,t_{i+1}]}\) is the indicator function in the interval \([t_i,t_{i+1}]\)

The number \(N\) is given by the number of elements in the numpy array u used when creating th epulse, and therefore \(\Delta t\) is computed accordingly, given \(T\).

In fact, although still experimental and therefore undocumented, one can also use some interpolating function to obtain the values in the interior of the intervals, so that the function is smooth
enveloped_fourier

This is a normal Fourier expansion, but multiplied by a function \(S(t)\):

\[f(u, t) = S(t) g_{\rm Fourier}(u, t)\]

The function \(S(t)\) must be passed to the function assign_envelope_function(), and is simply a function of time:
```
def S(t):
    return np.sin(t)
f = pulses.pulse('enveloped_fourier', T, u, bound = kappa)
f.assign_envelope_function(S)
```

user_defined

This is the most general possibility: the user defines the pulse as a a python function. The function must be communicated to the pulse object using the assign_user_defined_function(). One should also define a function with the derivatives of the function with respect to the parameters. Here is an example using only three of parameters:

\[f(u, t) = u_0 \cos(u_1 t+ u_2)\]

def fu(t, u):
    return u[0] * np.cos(u[1]*t + u[2])
def dfdu(t, u, m):
    # m is an integer that signals what is the parameter with respect to which the
    # derivative is taken
    if m == 0:
        return np.cos(u[1]*t + u[2])
    elif m == 1:
        return -u[0] * t * np.sin(u[1]*t + u[2])
    elif m == 2:
        return -u[0] * np.sin(u[1]*t + u[2])
u = np.array([1.0, 0.5, 0.0])
f = pulses.pulse('user_defined', T, u)
f.assign_user_defined_function(fu, dfdu)

If one is going to use a gradient-less QOCT algorithm, there is no need to define the derivative function, and one can just pass the None value.

Parameters:

type (string) – The type of function: ‘fourier’, ‘bound_fourier’, ‘realtime’, ‘user_defined’, or ‘enveloped’ fourier
T (float) – The duration of the pulse
u (ndarray) – The parameters
bound (float, default = None) – For the ‘bound_fourier’ type, the bound

type

Type:: string

T

Type:: float

u

Type:: ndarray

nu

Type:: int

bound

Type:: float

assign_envelope_function(efunc)[source]: Sets the envelope function that may be multiplied by the pulse itself.

assign_user_defined_function(fu, dfu)[source]: For user-defined functions, set the function definition, and the definition of the gradient of the function.

dfu(t, m, u=None)[source]: Derivative of the pulse with respect to the m-th parameter, at time t

envelope(t)[source]: Returns the value of the envelope function at time t

f(t, args)[source]: The value of the pulse at time t

fitparams(fref, times, u0)[source]: Fit the parameters of a function to the best mpossible much wrt a reference function

fu(t, u=None)[source]

Returns the value of pulse at time t

Optionally, one can pass the control parameters u and in that way those are updated before the computation.

fw(t, u, w=None, explicit_formula=False)[source]: Fourier transform of the pulse

gradf(t)[source]

Grdient of the function at time t

“Gradient of the function” means the gradient with respect to all the control parameters.

print(filename)[source]

Prints a pulse to a file called ‘filename’.

It does not print the possible constraint functions.

Parameters:: filename (str) – Name of the file where the info about the pulse is written.

set_constraint(which, val=0.0)[source]: Sets a constraint on the shape of the pulse

set_parameters(u)[source]: Sets the parameters \(u = u_1, \dots, u_P\) on the pulse.

pulses.pulse_collection_get_parameters(f)[source]

Sets the parameters of a collection of pulses

It receives a numpy array of parameters u, whose dimension should be equal to the summ of parameters of the pulses in the collection f.

pulses.pulse_collection_parameter_range(f, j)[source]

Returns the ‘parameter range’ of a given pulse in a collection

Given a list of pulses ‘f’, the full control parameter list will be an array joining all the control parameters of each of of them. This function returns the indexes that would correspond, in that array, to one of the pulses.

Thus, for example, if we have two pulses with nu1 and nu2 parameters each, the range of the first one would be (0, nu1), whereas the range of the second would be (nu1, nu1+nu2).

Parameters:

f (list of pulse) – A list with pulse objects
j (int) – The pulse for which we want to get the parameter range.

Returns:

a list of two integer numbers, with the starting index of the parameters corresponding to the first pulse, and the final index.

Return type:

list of int

pulses.pulse_collection_set_parameters(f, u)[source]

Sets the parameters of a collection of pulses

It receives a numpy array of parameters u, whose dimension should be equal to the summ of parameters of the pulses in the collection f.

pulses.pulse_constraint_functions(f, param_range)[source]

Given a pulse, returns a list with all its constraint functions.

The constraint functions contained in the pulse object are functions of the nu control parameters on which the pulse depends. In contrast, the functions returned by this function are functions of all the control parameters in a list of pulses. Therefore, we need as an input a ‘parameter range’, i.e. the starting and finishin indexes of the parameters of this particular pulse in the list.

pulses.read_pulse(filename)[source]

Reads a pulse from the info contained in file.

Parameters:: filename (str) – Name of the file where the info about the pulse is written.
Returns:: The pulse whose information was previously written to filename.
Return type:: pulse

pulses.rotation(theta, n, dim=2, i=None, j=None)[source]

Returns a rotation operator between two levels, given an angle an axis.

The rotation affects level i and j (or the only two levels, if dim = 2, which is the default). The angle of rotation is given by theta, and the axis is given by the unitary vector n.

Parameters:

theta (float) – The rotation angle.
n (ndarray) – A three-dimensional float array containing a unit vector.
dim (int, default = 2) – The dimension of the rotation operator that will be created.
i (int, default = None) – One of the states affected by the rotation.
j (int, default = None) – One of the states affected by the rotation.

Returns:

A Qobj operator with the rotation operator.

Return type:

Qobj

pulses.rotationpulse(axis, A, omega, mu0, theta, t0=None)[source]

Returns a pulse that implements a rotation betwen two states.

Given an amplitude for a pulse A, and two levels characterized by an energy difference omega and a coupling mu0, creates a pulse object that implements a rotation between those two levels in direction given by axis (0 => x, 1 => y, 2 => z), and angle theta.

Parameters:

axis (int) – The rotation index (0 => x, 1 => y, 2 => z).
A (float) – The amplitude of the pulse.
omega (float) – The frequency difference between the levels.
mu0 (float) – The coupling matrix element
theta (float) – The rotation angle

Returns:

A pulse object with the pulse that produces the rotation.

Return type:

pulse

6.4. solvers

Solvers for the Schrödinger/Lindblad equation

The code may use the qutip solvers, or use some internal solvers that may be faster in some circumstances.

solvers.cfmagnus2solver(H, f, y0, time, returnQoutput=True, interaction_picture=False, cops=None)[source]

Implementation of the exponential midpoint rule.

It is also the second-order commutator-free Magnus method.

solvers.cfmagnus4solver(H, f, psi0, time, returnQoutput=True, interaction_picture=False, cops=None)[source]: Implementation of the fourth order commutator-free Magnus method.

solvers.solve(solve_method, H, f, u, y0, time)[source]

This function propagate the object y0 using either the Schrödinger equation with H = H0 + f(t)*V as Hamiltonian, or Lindblad’s equation. If an inverse time array is given, the solver take y0 as the state at final time, making a backward propagation.

The problem may be solved in the interaction picture if interaction_picture = True. Note that in that case H0 should be diagonal on entry (the system should be represented in the eigenbasis of H0). In this case, the output will also be in the interaction picture.

Parameters:

solve_method (string indicating the propagation method) –
H0 (hamiltonian) – Hamiltonian’s time independet part
V (list of functions) – Hamiltonian’s perturbation component
f (list of pulse) – pulse class defined in typical_pulse.py based on Fourier expansion parametrization
y0 – initial state
time (ndarray) – array that contain each time step.
u (ndarray) – array that contain the control parameters of the pulse, needed to build the args diccionary defined in pulse class
Returns –
....... –
qutip.result – qutiip.result type object with the propagation data

solvers.solve_h_is_function(solve_method, H, f, y0, time, returnQoutput=True, options=None, interaction_picture=False)[source]: Solver for the case in which the hamiltonian is given as a function

6.5. target

This module contains the Target class, and associated procedures.

class target.Target(targettype, Fyu=None, dFdy=None, dFdu=None, operator=None, Utarget=None, alpha=None, S=None, targeteps=None, T=None, fepsilon=None, dfdepsilon=None)[source]

The class that holds the definition of the target.

The user can choose between (1) writing a Python function that defines the target functional (in general, it must be accompanied of a Python function that computes the functional derivative with respect to the wavefunction, evolution operator or density matrix), and (2) using one of the predefined target types defined by the code. The former is chosen by setting the targettype to “generic”, whereas for the latter one must use either “expectationvalue”, “evolutionoperator”, or “floquet”.

Parameters:

targettype (string) –
The code admits the following types of targets:
1. generic
  
  The user may input any function as target, defined by the driver program, using the Fyu, dFdy and dFdu parameters. See below for information about the proper definition of these functions.
2. expectationvalue
  
  The target is the expectation value of some operator:
  
  \[F(\Psi, u) = \langle \Psi(T)\vert O \vert \Psi(T)\rangle\]
  
  or, if the multitarget option is used,
  
  \[F(\Psi, u) = \frac{1}{Q} \sum_i \langle \Psi_i(T)\vert O_i \vert \Psi_i(T)\rangle\]
  
  Here, \(Q\) is the number of targets. The operator(s) should be provided using the operator argument (below). If one is using density matrices instead of wave functions to describe open systems, then:
  
  \[F(\rho, u) = {\rm Tr}{\rho(T)O}\]
  
  or, if the multitarget option is used,
  
  \[F(\rho_i, u) = \frac{1}{Q} \sum_i {\rm Tr}{\rho_i(T)O_i}\]
  
  This option is simple to use, one just needs to instantiate a Target class object as:
```
tg = target.Target('expectationvalue', operator = O)
```
  where “O” is a qutip object containing the operator whose expectation value is to be maximized. It must be a list of operators, if optimizing in multitarget mode.
3. evolutionoperator
  
  The goal now is to find a set of control parameters that leads to an evolution of the system characterized by a fixed target evolution operator. This target evolution operator (or operators, if a multitarget setup is used), is provided by the argument Utarget. The target definition is:
  
  \[F(U, u) = \frac{1}{Q} \sum_i \vert U(T) \cdot U^{(i)}_{\rm target} \vert^2\]
  
  In this case, the creation of the target object would be:
```
tg = target.Target('evolutionoperator', Utarget = U)
```
  where U would be the (list of) unitary operators.
4. floquet
  
  (undocumented)
Fyu (function) –
If targettype == “generic”, the user should pass the function used to define the target. The interface should be:
```
def Fyu(y: qutip.Qobj or list of qutip.Qobj, u: ndarray) -> float:
```
The function Fyu should receive as argument a system state (be it a Hilbert state, evolution operator, or density matrix), that is supposed to be the state at the end of the propagation, and a control parameters set. It should output the target functional value. It may also take as argument a list of states, if one is running the optimizaion in multitarget mode.
dFdy (function) –
If targettype == “generic”, the user should pass the function used to define the target (Fyu), and the derivative of that function with respect to the state, which would be given by this dFdy parameter. The interface should be:
```
def dFdy(y: qutip.Qobj or list of qutip.Qobj, u: ndarray) -> list of float:
```
dFdu (function) –
If targettype == “generic”, the user should pass the function used to define the target (Fyu), and, if it is not zero, the derivative of that function with respect to the control parameters, which would be given by this dFdu parameter. The interface should be:
```
def dFdu(u: ndarray, m: int) -> float:
```
operator (qutip.Qobj or list of qutip.Qobj) – If targettype == “expectationvalue”, this argument should be the operator whose expectation value is to be maximized. In multitarget mode, this would be a list of operators, one for each system.
Utarget (qutip.Qobj or list of qutip.Qobj) – If targettype == “evolutionoperator”, this argument should be the unitary operator that is set as target. In multitarget mode, this would be a list of unitaries, one for each system.

Examples

Let us show a simple example of the use of the “generic” target type. The following would be equivalent to using the “expectationvalue” targe type with operator O:

def Fpsi(psi, u):
    return qt.expect(O, psi)
def dFdpsi(psi, u):
    return O * psi
tg = target.Target("generic", Fyu = Fpsi, dFdy = dFdpsi)

Fyu(y, u)[source]

The functional F of the trayectory y that is to be maximized

It may also be an explicit function of the control parameters u.

Parameters:

y (qutip.result) – The trajectory of the quantum system
u (ndarray) – The control parameters that were used to generate the trajectory.

Returns:

The value of the target functional.

Return type:

float

dFdy(y, u)[source]

Derivative of the target functional wrt the quantum trajectory.

Right now, it in fact assumes that the functional only depends on the final state of the quantum trajectory. This functional derivative is used to determine the boundary condition used in the definition of the costate.

Parameters:

y (qutip.Qobj or list of qutip.Qobj) – The state of the system at the final time of the propagation
u (ndarray) – The control parameters

Returns:

The derivative of F with respect to the quantum state.

Return type:

qutip.Qobj or list of qutip.Qobj

6.6. floquet

This module contains all the procedures needed to do Floquet optimization.

class floquet.Nessopt(target_operator, T, nts, L, glist)[source]: class used for optimization of NESSs

floquet.depsi(UT, T, delta=0.0001)[source]

Given a unitary UT, returns the derivatives of the Floquet quasienergies with respect to it. Specifically, it returns the derivative of each Floquet quasienergy with respect to each matrix element \(U_{ij}^*\) (technically, a Wirtinger derivative, with respect to the complex conjugate of the parameter).

Parameters:

UT (ndarray) – A unitary matrix, representing an evolution operator.
T (float) – The Floquet period
delta (float, default = 1.0e-4) – The small displacement used to compute the derivatives through a finite-difference formula.

Returns:

The numpy array containing the derivatives; the shape is (dim, dim, dim): depsi[alpha, :, :] is the matrix of derivatives of the alpha-th quasienergy with respect to all the elements of the unitary.

Return type:

ndarray

floquet.dmtovector(rho)[source]

Transforms a Qobj [[d], [d]] shaped operator object into a [[d^2], [1]] shaped ket.

Parameters:: rho (Qobj) – Density matrix
Returns:: A Qobj vector containing the vectorized version of the input density matrix
Return type:: Qobj

floquet.epsi(UT, T)[source]

Given a matrix UT, returns the Floquet quasienergies associated to it.

In principle, the matrix should be unitary (or otherwise, there is no guarantee that the matrix will have unit eigenvalues). However, the routine accepts non-unitary matrix, as it computes the Schur decomposition, and it uses the diagonal elements of the Schur form. If the matrix is indeed unitary, the output will be proper Floquet quasienergies.

Parameters:

UT (ndarray) – A unitary matrix, representing an evolution operator.
T (float) – The Floquet period

Returns:

The Floquet quasienergies (ordered, and in the Floquet Brillouin zone)

Return type:

ndarray

floquet.epsilon(H, T, args=None)[source]

Computes Floquet quasienergies corresponding to the Hamiltonian H and the period T.

For this to make sense, the Hamiltonian time-dependence must be periodic. The Hamiltonian must be given as a “qutip Hamiotonian”, i.e. anything that could be accepted by the qutip propagation functions as valid Hamiltonians.

Parameters:

H (qutip Hamiltonian) – The (periodic) Hamiltonian
T (float) – The period
args (dict, default = None) – An “args” argument that can b passed to the qutip sesolve function that is used to compute the unitary evolution operator at time T.

Returns:

The ordered set of quasienergies.

Return type:

ndarray

floquet.epsilon3(H, f, u, T)[source]

Given an Hamiltonian H (hamiltonians.hamiltonian), a set of driving periodic perturbations f, returns the Floquet quasienergies

Parameters:

H (hamiltonians.hamiltonian) – The Hamiltonian of the system, given as a hamiltonian class object.
f (list of pulses.pulse) – A list of periodic drivings – the number should match the number of perturbations of the Hamiltonian
u (ndarray) – The control parameters
T (float) – The period

Returns:

An array with the Floquet quasienergies

Return type:

ndarray

floquet.gradepsilon(H, f, u, T)[source]

Returns the gradient of the Floquet quasienergies with respect to the control parameters of the set of periodic drivings.

Parameters:

H (hamiltonians.hamiltonian) – The Hamiltonian of the system, given as a hamiltonian class object.
f (list of pulses.pulse) – A list of periodic drivings – the number should match the number of perturbations of the Hamiltonian
u (ndarray) – The control parameters
T (float) – The period

Returns:

An array shaped (dim, nu) where dim is the Hilert space dimension, and nu is the number of control parameters, containing the derivatives of the Floquet quasienergies with respect to those parameters.

Return type:

ndarray

floquet.ness(T, ntsteps, L, f, compute_gradient=False)[source]

Computes the NESS and, if required, the gradient of the NESS with respect to the control parameters of the periodic drivings.

See FLoquet optimizations: non-equilibrium steady states for details about how the gradient is computed.

floquet.vectortodm(rhovec)[source]

Transforms a numpy array of dimension (d^2) into an “operator” Qobj object; this object has dimension [[d], [d]]. The input 1D array is supposed to store the operator with stacked columns.

Parameters:: rhovec (ndarray) – A one-dimensional numpy complex array with d^2 elements
Return type:: A [[d], [d]] shaped qutip Qobj

6.7. krotov

This module implements the Krotov’s algorithm [RNK12].

The current implementation is rather restrictive, and it can only be applied when:

The state is a wavefunction (no density matrices or evolution operators).
The Hamiltonian is linear in the control field, with only one perturbation, i.e. it has the form:

\[H(t) = H_0 + f(t) V\]
The pulse is parametrized in real time.
The only allowed propagators are RK4 and cfmagnus4.
The target has the form:

\[F(\psi(T)) = \langle \psi(T)\vert O \vert\psi(T)\rangle - \int_0^T\!\!{\rm d}t\; \frac{\alpha}{S(t)}f^2(t)\]

i.e. the objective is to maximize the expectation value of some operator \(O\), and we have to add a penalty function, depending on the parameter \(\alpha > 0\) and on the weigth function \(S(t)\), that can be supplied by the user.

The algorithm that is implemented here is the following:

\(k = 0\); \(f^{(0)}\) is the initial guess pulse. Solve:

\[i\frac{\partial}{\partial t} \psi^{(0)}(t) = [H_0 + f^{(0)}V(t)]\psi^{0}(t)\]

\[\psi^{(0)}(0) = \psi_0\]

\[i\frac{\partial}{\partial t} \chi^{(0)}(t) = [H_0 + f^{(0)}V(t)]\chi^{0}(t)\]

\[\chi^{(0)}(T) = O\psi^{(0)}(T)\]
For \(k=1,\dots\) until convergence:

\[i\frac{\partial}{\partial t} \psi^{(k)}(t) = [H_0 + \frac{S(t)}{\alpha} {\rm Im} [\langle \chi^{(k-1)}\vert V \vert\psi^{(k)}(t)\rangle ]\psi^{(k)}(t)\]

\[\psi^{(k)}(0) = \psi_0\]

\[f^{(k)}(t) = \frac{S(t)}{\alpha} {\rm Im} [\langle \chi^{(k-1)}\vert V \vert\psi^{(k)}(t)\rangle ]\]

\[i\frac{\partial}{\partial t} \chi^{(k)}(t) = [H_0 + f^{(k)}V(t)]\chi^{k}(t)\]

\[\chi^{(k)}(T) = O\psi^{(k)}(T)\]

Convergence is checked by comparing \(f^{(k+1)}\) with \(f^{(k)}\).

krotov.cfmagnus4_nonlinear(H0, V, f, chi, psi0, time, S, alpha=2, interaction_picture=False, cops=None)[source]: This function solves the nonlinear schordinger equation using cfmagnus4 as the integrator algorithm.

krotov.check_pulse(f, times)[source]

This function checks if the pulse parametrization is realtime.

If it is not, it convert the pulse into realtime parametrization.

krotov.coestate_eq(solve_method, H, f, psiF, times, O, obj, interaction_picture=False)[source]

This function calculate the coestate equation with the final condition described by

\[\vert\chi(T)\rangle = O\vert \Psi(t=T)\rangle\]

In the case of having a state propagating, and

\[\vert B(T)\rangle = Tr [U(T)+*O]*O/d^2\]

in the case of having an operator propagating.

The main idea is to solve it as an linear schrödinger equation but backwards, so we have to take care of the time and the pulse in order to make this correctly (with realtime parametrized pulses only is needed to invert the time). In principle this should work with any solve_method: rk4, cfm2 or cmf4.

krotov.coestate_interpolator(chi, times, coestate_flag, solve_method)[source]

This function interpolates the value of the coestate.

In this way, we can solve properly the non-linear Schrodinger equation with the rk4 and cfmagnus4 integrators. This function returns the value of the coestate in the time that´s required

krotov.have_converged(old, new, times, epsilon=1e-06)[source]

This function checks if the algorithm have converged enough.

The idea is to check every time of the pulse.

Parameters:

old – The pulse before solving the non-linear schrodinger
new – The pulse after solving the non-linear schrodinger
epsilon – The tolerance we have so we can say the pulse have converged enough

Returns:

A boolean which is ‘True’ if every amplitude have converged, and is ‘False’ if any
of the amplitudes have not converged enough

krotov.krotov(solve_method, H, f_initial, psi0, times, O, S, alpha=2.0, tolerance=0.001, verbose=False, maxeval=100, interaction_picture=False)[source]

This is the main function on the module.

It performs the optimization and returns the optimal value of the functional, the optimal parameters, and a code signaling the success or failure of the process.

Parameters:

solve_method (string) – string indicating the propagation method
H (hamiltonian) – Hamiltonian
f_initial (pulse) – pulse to be used as initial guess (on output, it will contain the optimized pulse.
psi0 (Qobj) – initial state.
time (ndarray) – array that contains each time step.
O (Qobj) – The target operator that we want to optimize.
S (function) – The penalty function \(S(t)\)
alpha (float, default = 2.0) – Penalty factor \(alpha\)
tolerance (float, default = 1.0e-3) – Precision for the algorithm
verbose (bool, default = True) – Whether or not the algorithm will print info to stdout
maxeval (int, default = 100) – The maximum number of iterations that will be done before the algorithm stops.
interaction_picture (bool, default = False) – The code will transform the equations to the interaction picture if this is set to True

Returns:

ndarray – The optimal parameters
float – The optimal target function value
list – A list with data about the convergence history of the algorith
int – An integer flag about the outcome: zero if successful, non-zero otherwise

krotov.rk4_nonlinear(H0, V, f, chi, psi0, times, S, alpha=2, interaction_picture=False)[source]

This function solves the nonlinear schordinger equation using rk4 as the integrator algorithm.

It is a modification of the function rk4solver in the solvers module of qocttools, so it can solve a slightly different equation and it can return the new state and the new pulse.

krotov.state_eq(solve_method, H0, V, f, chi, psi0, times, S, alpha=2, interaction_picture=False)[source]: This function solves the non-linear Schrödinger equation that is part of Krotov’s algorithm.

6.8. math_extra

math_extra.diff_ridders(function, origin, stepsize, con=1.4, safe=2.0, maxiter=15)[source]

Estimate first-order derivative with Ridders’ finite difference method.

This implementation is based on the one from the book Numerical Recipes. The code is pythonized and no longer using fixed-size arrays. Also, the output of the function can be an array.

Parameters:

function (function) – The function to be differentiated.
origin (float) – The point at which must be differentiated.
stepsize (float) – The initial step size.
con (float) – The rate at which the step size is decreased (contracted). Must be larger than one.
safe (float) – The safety check used to terminate the algorithm. If Errors between successive orders become larger than safe times the error on the best estimate, the algorithm stop. This happens due to round-off errors.
maxiter (int) – The maximum number of iterations, equals the maximum number of function calls and also the highest polynomial order in the Neville method.

Returns:

estimate (float) – The best estimate of the first-order derivative.
error (float) – The (optimistic) estimate of the error on the derivative.

math_extra.frobenius_product(A, B)[source]: return the Frobenius product between the operators A and B

math_extra.infidelity(U, dmi, dmo)[source]

Returns a measure of the infidelity of a quantum process.

(…)

math_extra.timegrid(H0, T, delta)[source]: Generates a time grid in the interval [0, T]