FLiMESolve

VERSION

@@ -1 +1 @@
-5.0.0.dev
+5.0.0a1

doc/changes/2186.feature

@@ -0,0 +1 @@
+Added Floquet-Lindblad Master Equation for sinuisoidally driven time-dependent open system dynamics

doc/guide/dynamics/dynamics-floquet.rst

@@ -93,7 +93,7 @@ Consider for example the case of a strongly driven two-level atom, described by
 In QuTiP we can define this Hamiltonian as follows:
 
 .. plot::
-   :context: reset
+   :context: close-figs
 
    >>> delta = 0.2 * 2*np.pi
    >>> eps0 = 1.0 * 2*np.pi
@@ -107,7 +107,7 @@ In QuTiP we can define this Hamiltonian as follows:
 The :math:`t=0` Floquet modes corresponding to the Hamiltonian :eq:`eq_driven_qubit` can then be calculated using the :class:`qutip.FloquetBasis` class, which encapsulates the Floquet modes and the quasienergies:
 
 .. plot::
-   :context: close-figs
+   :context:
 
    >>> T = 2*np.pi / omega
    >>> floquet_basis = FloquetBasis(H, T, args)
@@ -131,7 +131,7 @@ For certain driving amplitudes the quasienergy levels cross.
 Since the quasienergies can be associated with the time-scale of the long-term dynamics due that the driving, degenerate quasienergies indicates a "freezing" of the dynamics (sometimes known as coherent destruction of tunneling).
 
 .. plot::
-   :context: close-figs
+   :context:
 
    >>> delta = 0.2 * 2 * np.pi
    >>> eps0  = 0.0 * 2 * np.pi
@@ -175,7 +175,7 @@ The purpose of calculating the Floquet modes is to find the wavefunction solutio
 To do that, we first need to decompose the initial state in the Floquet states, using the function :meth:`FloquetBasis.to_floquet_basis`
 
 .. plot::
-   :context: close-figs
+   :context:
 
    >>> psi0 = rand_ket(2)
    >>> f_coeff = floquet_basis.to_floquet_basis(psi0)
@@ -186,7 +186,7 @@ To do that, we first need to decompose the initial state in the Floquet states,
 and given this decomposition of the initial state in the Floquet states we can easily evaluate the wavefunction that is the solution to :eq:`eq_driven_qubit` at an arbitrary time :math:`t` using the function :meth:`FloquetBasis.from_floquet_basis`:
 
 .. plot::
-   :context: close-figs
+   :context:
 
    >>> t = 10 * np.random.rand()
    >>> psi_t = floquet_basis.from_floquet_basis(f_coeff, t)
@@ -284,7 +284,24 @@ Finally, :func:`qutip.solver.floquet.fmmesolve`  always expects the ``e_ops`` to
     output = fmmesolve(H, psi0, tlist, [sigmax()], e_ops=[num(2)], spectra_cb=[noise_spectrum], T=T, args=args)
     p_ex = output.expect[0]
 
-.. plot::
-    :context: reset
-    :include-source: false
-    :nofigs:
+The Floquet-Lindblad master equation in QuTiP
+-------------------------------------------
+The QuTiP function :func:`qutip.floquet.flimesolve` implements the Floquet-Lindblad master equation. It calculates the dynamics of a system given its initial state, a time-dependent Hamiltonian, a list of operators through which the system couples to its environment and a list of corresponding time-independent decay rates for the associated operators. Flimesolve is much like the :func:`qutip.mesolve` and :func:`qutip.mcsolve`, in how the system couples to the environment.
+
+Flimesolve expects the collapse operators and rates to be provided as a list of [c_op, rate] pairs. For example:
+
+.. code-block:: python
+
+    gamma1 = 0.1
+    c_ops_and_rates = [[sigmax(), gamma1]]
+
+The other parameters are similar to the :func:`qutip.mesolve` and :func:`qutip.mcsolve`, and the same format for the return value is used :class:`qutip.solve.solver.Result`. The following example extends the example studied above, and uses :func:`qutip.floquet.flimesolve` to introduce dissipation into the calculation
+
+.. plot:: guide/scripts/floquet_ex4.py
+   :width: 4.0in
+   :include-source:
+
+Finally, :func:`qutip.solver.floquet.fmmesolve`, similar to :func:`qutip.solver.floquet.flimesolve`, always expects the ``e_ops`` to be specified in the laboratory basis:
+
+    output = flimesolve(H, psi0, tlist, [[sigmax(),gamma1]], e_ops=[num(2)], T=T, args=args)
+    p_ex = output.expect[0]

doc/guide/scripts/floquet_ex3.py

@@ -35,7 +35,6 @@ def noise_spectrum(omega):
     # Alternatively
     psi_t = output.states[idx]
     p_ex[idx] = qutip.expect(qutip.num(2), psi_t)
-
 # For reference: calculate the same thing with mesolve
 output = qutip.mesolve(H, psi0, tlist,
                        [np.sqrt(gamma1) * qutip.sigmax()], [qutip.num(2)],

doc/guide/scripts/floquet_ex4.py

@@ -0,0 +1,53 @@
+import numpy as np
+from matplotlib import pyplot
+import qutip
+
+delta = 0.0 * 2 * np.pi
+eps0 = 1.0 * 2 * np.pi
+A = 0.25 * 2 * np.pi
+omega = 1.0 * 2 * np.pi
+T = 2 * np.pi / omega
+tlist = np.linspace(0.0, 20 * T, 301)
+psi0 = qutip.basis(2, 0)
+
+H0 = - delta/2.0 * qutip.sigmax() - eps0/2.0 * qutip.sigmaz()
+H1 = A/2.0 * qutip.sigmax()
+args = {'w': omega}
+H = [H0, [H1, lambda t, w: np.sin(w * t)]]
+
+# collapse [operator, rate] pair
+gamma1 = 0.1
+c_op_and_rate = [[qutip.sigmax(), gamma1]]
+
+# solve the floquet-markov master equation
+output = qutip.flimesolve(
+    H, psi0, tlist,
+    c_ops_and_rates=[[qutip.sigmax(), gamma1]],
+    T=T,
+    args=args, options={"store_floquet_states": True}
+)
+
+# calculate expectation values in the computational basis
+p_ex = np.zeros(tlist.shape, dtype=np.complex128)
+for idx, t in enumerate(tlist):
+    f_coeff_t = output.floquet_states[idx]
+    psi_t = output.floquet_basis.from_floquet_basis(f_coeff_t, t)
+    # Alternatively
+    psi_t = output.states[idx]
+    p_ex[idx] = qutip.expect(qutip.num(2), psi_t)
+
+# For reference: calculate the same thing with mesolve
+output = qutip.mesolve(H, psi0, tlist,
+                       [np.sqrt(gamma1) * qutip.sigmax()], [qutip.num(2)],
+                       args)
+p_ex_ref = output.expect[0]
+
+
+# plot the results
+pyplot.plot(tlist, np.real(p_ex), 'r--', tlist, 1-np.real(p_ex), 'b--')
+pyplot.plot(tlist, np.real(p_ex_ref), 'r', tlist, 1-np.real(p_ex_ref), 'b')
+pyplot.xlabel('Time')
+pyplot.ylabel('Occupation probability')
+pyplot.legend(("Floquet $P_1$", "Floquet $P_0$",
+              "Lindblad $P_1$", "Lindblad $P_0$"))
+pyplot.show()

qutip/solver/correlation.py

@@ -11,11 +11,12 @@
     Qobj, QobjEvo, liouvillian, spre, unstack_columns, stack_columns,
     tensor, expect, qeye_like, isket
 )
+from .floquet import FloquetBasis
+from .flimesolve import FLiMESolver
 from .mesolve import MESolver
 from .mcsolve import MCSolver
 from .brmesolve import BRSolver
 from .heom.bofin_solvers import HEOMSolver
-
 from .steadystate import steadystate
 from ..ui.progressbar import progress_bars
 
@@ -346,7 +347,8 @@ def coherence_function_g1(
         n = solver.run(state0, taulist, e_ops=[a_op.dag() * a_op]).expect[0]
 
     # calculate the correlation function G1 and normalize with n to obtain g1
-    G1 = correlation_3op(solver, state0, [0], taulist, None, a_op.dag(), a_op)[0]
+    G1 = correlation_3op(
+        solver, state0, [0], taulist, None, a_op.dag(), a_op)[0]
 
     g1 = G1 / np.sqrt(n[0] * np.array(n))[0]
     return g1, G1
@@ -416,8 +418,23 @@ def coherence_function_g2(H, state0, taulist, c_ops, a_op, solver="me",
 
 
 def _make_solver(H, c_ops, args, options, solver):
-    H = QobjEvo(H, args=args)
-    c_ops = [QobjEvo(c_op, args=args) for c_op in c_ops]
+    if solver == "fme":
+        if isinstance(H, FloquetBasis):
+            floquet_basis = H
+        else:
+            T = options['T']
+            # are for the open system evolution.
+            floquet_basis = FloquetBasis(H, T, args, precompute=None)
+        time_sense = options.get('time sense', 0.)
+        solver_instance = FLiMESolver(
+            floquet_basis,
+            c_ops,
+            args,
+            time_sense=time_sense
+        )
+    else:
+        H = QobjEvo(H, args=args)
+        c_ops = [QobjEvo(c_op, args=args) for c_op in c_ops]
     if solver == "me":
         solver_instance = MESolver(H, c_ops, options=options)
     elif solver == "es":