Automatically convert TimeDependentOperatorSum to liouvillian (#235)

src/qobj/operator_sum.jl

@@ -3,7 +3,7 @@ export OperatorSum
 @doc raw"""
     struct OperatorSum
 
-A structure to represent a sum of operators ``\sum_i c_i \hat{O}_i`` with a list of coefficients ``c_i`` and a list of operators ``\hat{O}_i``.
+A constructor to represent a sum of operators ``\sum_i c_i \hat{O}_i`` with a list of coefficients ``c_i`` and a list of operators ``\hat{O}_i``.
 
 This is very useful when we have to update only the coefficients, without allocating memory by performing the sum of the operators.
 """
@@ -47,3 +47,7 @@ end
     end
     return y
 end
+
+function liouvillian(A::OperatorSum, Id_cache = I(prod(A.operators[1].dims)))
+    return OperatorSum(A.coefficients, liouvillian.(A.operators, Ref(Id_cache)))
+end

src/time_evolution/mesolve.jl

@@ -139,6 +139,10 @@ function mesolveProblem(
         is_empty_e_ops = isempty(e_ops)
     end
 
+    if H_t isa TimeDependentOperatorSum
+        H_t = liouvillian(H_t)
+    end
+
     p = (
         L = L,
         progr = progr,

src/time_evolution/time_evolution.jl

@@ -140,6 +140,10 @@ end
     return mul!(y, A.operator_sum, x, α, β)
 end
 
+function liouvillian(A::TimeDependentOperatorSum, Id_cache = I(prod(A.operator_sum.operators[1].dims)))
+    return TimeDependentOperatorSum(A.coefficient_functions, liouvillian(A.operator_sum, Id_cache))
+end
+
 #######################################
 
 ### LIOUVILLIAN ###

test/core-test/steady_state.jl

@@ -63,7 +63,7 @@
     e_ops = [a_d * a]
     psi0 = fock(N, 3)
     t_l = LinRange(0, 100 * 2π, 1000)
-    H_t_f = TimeDependentOperatorSum([(t, p) -> sin(t)], [liouvillian(H_t)])
+    H_t_f = TimeDependentOperatorSum((((t, p) -> sin(t)),), [H_t]) # It will be converted to liouvillian internally
     sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, H_t = H_t_f, progress_bar = Val(false))
     ρ_ss1 = steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetLinearSystem())[1]
     ρ_ss2 =