Merge branch 'simpler-bloch-redfield-tensor' of https://github.com/sd…

…109/QuantumOptics.jl into simpler-bloch-redfield-tensor

Project.toml

@@ -1,6 +1,6 @@
 name = "QuantumOptics"
 uuid = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
-version = "v0.8.7"
+version = "v0.8.8"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"

src/master.jl

@@ -105,6 +105,33 @@ function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
     end
 end
 
+function master(tspan, rho0::Operator, L::SuperOperator; fout=nothing, kwargs...)
+    # Rewrite rho as Ket and L as Operator
+    dim = length(rho0.basis_l)*length(rho0.basis_r)
+    b = GenericBasis(dim)
+    rho_ = Ket(b,reshape(rho0.data, dim))
+    L_ = Operator(b,b,L.data)
+    dmaster_(t,rho,drho) = dmaster_liouville(rho,drho,L_)
+
+    # Rewrite into density matrix when saving
+    tmp = copy(rho0)
+    if fout===nothing
+        fout_ = function(t,rho)
+            tmp.data[:] = rho.data
+            return copy(tmp)
+        end
+    else
+        fout_ = function(t,rho)
+            tmp.data[:] = rho.data
+            return fout(t,tmp)
+        end
+    end
+
+    # Solve
+    return integrate_master(tspan, dmaster_, rho_, fout_; kwargs...)
+end
+
+
 """
     timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)
 
@@ -163,7 +190,7 @@ end
 
 
 # Automatically convert Ket states to density operators
-master(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master(tspan, dm(psi0), H, J; kwargs...)
+master(tspan, psi0::Ket{B}, args...; kwargs...) where B<:Basis = master(tspan, dm(psi0), args...; kwargs...)
 master_h(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_h(tspan, dm(psi0), H, J; kwargs...)
 master_nh(tspan, psi0::Ket{B}, Hnh::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_nh(tspan, dm(psi0), Hnh, J; kwargs...)
 master_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_dynamic(tspan, dm(psi0), f; kwargs...)
@@ -176,12 +203,12 @@ function recast!(x::T, rho::Operator{B,B,T}) where {B<:Basis,T}
 end
 recast!(rho::Operator{B,B,T}, x::T) where {B<:Basis,T} = nothing
 
-function integrate_master(tspan, df::Function, rho0::T,
-                        fout::Union{Nothing, Function}; kwargs...) where {B<:Basis,T<:Operator{B,B}}
+function integrate_master(tspan, df::Function, rho0,
+                        fout::Union{Nothing, Function}; kwargs...)
     tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
     x0 = rho0.data
-    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
-    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
+    state = deepcopy(rho0)
+    dstate = deepcopy(rho0)
     integrate(tspan_, df, x0, state, dstate, fout; kwargs...)
 end
 
@@ -282,6 +309,11 @@ function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOper
     return drho
 end
 
+function dmaster_liouville(rho,drho,L)
+    mul!(drho,L,rho)
+    return drho
+end
+
 function dmaster_h_dynamic(t, rho::T, f::Function,
                     rates::DecayRates,
                     drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}

test/test_timeevolution_master.jl

@@ -177,4 +177,9 @@ tout, Ψket_t = timeevolution.schroedinger(T, Ψ₀, Hdense; reltol=1.e-7)
 tout, Ψbra_t = timeevolution.schroedinger(T, dagger(Ψ₀), Hdense; reltol=1.e-7)
 @test tracedistance(Ψket_t[end]⊗Ψbra_t[end], ρ) < 1e-5
 
+L = liouvillian(Hdense, J)
+tout, rho_t = timeevolution.master(T, Ψ₀, L)
+tout, h_exp = timeevolution.master(T, dm(Ψ₀), L; fout=(t,rho)->expect(Hdense, rho)::ComplexF64)
+@test isequal(expect(Hdense, rho_t), h_exp)
+
 end # testset