Merge pull request #1952 from christian512/floquet_transform

Add missing transform to Floquet-Markov Solver in v5

doc/changes/1952.bugfix

@@ -0,0 +1 @@
+Add missing state transformation to floquet_markov_mesolve

qutip/solve/floquet.py

@@ -906,7 +906,7 @@ def transform(rho, t):
         if not r.successful():
             break
 
-        rho = Qobj(unstack_columns(r.y), rho0.dims)
+        rho = transform(Qobj(unstack_columns(r.y), rho0.dims), t)
 
         if expt_callback:
             e_ops(t, rho)

qutip/tests/solve/test_floquet.py

@@ -191,6 +191,70 @@ def noise_spectrum(omega):
 
         np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
 
+    def testFloquetMasterEquation3(self):
+        """
+        Test Floquet-Markov Master Equation for a two-level system
+        subject to dissipation with internal transform of fmmesolve
+        """
+
+        delta = 1.0 * 2 * np.pi
+        eps0 = 1.0 * 2 * np.pi
+        A = 0.5 * 2 * np.pi
+        omega = np.sqrt(delta ** 2 + eps0 ** 2)
+        T = (2 * np.pi) / omega
+        tlist = np.linspace(0.0, 2 * T, 101)
+        psi0 = rand_ket(2)
+        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
+        e_ops = [num(2)]
+        gamma1 = 1
+
+        A = 0. * 2 * np.pi
+        psi0 = rand_ket(2)
+        H1 = A / 2.0 * sigmax()
+        args = {'w': omega}
+        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
+
+        # Collapse operator for Floquet-Markov Master Equation
+        c_op_fmmesolve = sigmax()
+
+        # Collapse operator for Lindblad Master Equation
+        def noise_spectrum(omega):
+            if omega > 0:
+                return 0.5 * gamma1 * omega/(2*np.pi)
+            else:
+                return 0
+
+        ep, vp = H0.eigenstates()
+        op0 = vp[0]*vp[0].dag()
+        op1 = vp[1]*vp[1].dag()
+
+        c_op_mesolve = []
+        gamma = np.zeros([2, 2], dtype=complex)
+        for i in range(2):
+            for j in range(2):
+                if i != j:
+                    gamma[i][j] = 2*np.pi*c_op_fmmesolve.matrix_element(
+                        vp[j], vp[i])*c_op_fmmesolve.matrix_element(
+                        vp[i], vp[j])*noise_spectrum(ep[j]-ep[i])
+
+        for i in range(2):
+            for j in range(2):
+                c_op_mesolve.append(np.sqrt(gamma[i][j])*(vp[i]*vp[j].dag()))
+
+        # Solve the floquet-markov master equation
+        output1 = fmmesolve(
+                            H, psi0, tlist, [c_op_fmmesolve], [num(2)],
+                            [noise_spectrum], T, args, floquet_basis=False)
+        p_ex = output1.expect[0]
+        # Compare with mesolve
+        output2 = mesolve(H, psi0, tlist, c_op_mesolve, [num(2)], args)
+        p_ex_ref = output2.expect[0]
+
+        np.testing.assert_allclose(np.real(p_ex), np.real(p_ex_ref), atol=1e-4)
+
     def testFloquetRates(self):
         """
         Compare rate transition and frequency transitions to analytical