Adapt to new version of SciMLOperators.jl

src/qobj/superoperators.jl

@@ -39,42 +39,29 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(B_T, Id)
 end
 
-## intrinsic liouvillian
-_liouvillian(
+## intrinsic liouvillian 
+function _liouvillian(
     H::MT,
     Id::AbstractMatrix,
-    assume_hermitian::Val{true},
-) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_Herm(H, Id)
-_liouvillian(
-    H::MT,
-    Id::AbstractMatrix,
-    assume_hermitian::Val{false},
-) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} = _liouvillian_NonHerm(H, Id)
-
-## intrinsic liouvillian (assume H is Hermitian)
-_liouvillian_Herm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
-    -1im * (_spre(H, Id) - _spost(H, Id))
-_liouvillian_Herm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_Herm(H.A, Id))
-_liouvillian_Herm(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian_Herm(H.L, Id))
-_liouvillian_Herm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_Herm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
-
-## intrinsic liouvillian (assume H is Non-Hermitian)
-function _liouvillian_NonHerm(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
-    CFType = _complex_float_type(H)
-    return CFType(-1.0im) * _spre(H, Id) + CFType(1.0im) * _spost(H', Id)
+    ::Val{assume_hermitian},
+) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator},assume_hermitian}
+    H_spre = _spre(H, Id)
+    H_spost = assume_hermitian ? _spost(H, Id) : _spost(H', Id)
+    return -1im * (H_spre - H_spost)
 end
-_liouvillian_NonHerm(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian_NonHerm(H.A, Id))
-function _liouvillian_NonHerm(H::ScaledOperator, Id::AbstractMatrix)
-    CFType = _complex_float_type(H)
-    return CFType(-1.0im) * ScaledOperator(H.λ, _spre(H.L, Id)) +
-           CFType(1.0im) * ScaledOperator(conj(H.λ), _spost(H.L', Id))
+function _liouvillian(H::MatrixOperator, Id::AbstractMatrix, assume_hermitian::Val)
+    isconstant(H) || error("Liouvillian can only be constructed from constant MatrixOperator.")
+    return MatrixOperator(_liouvillian(H.A, Id, assume_hermitian))
 end
-_liouvillian_NonHerm(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian_NonHerm(op, Id), +, H.ops) # use mapreduce to sum (+) all ops
+_liouvillian(H::ScaledOperator, Id::AbstractMatrix, assume_hermitian::Val{true}) =
+    ScaledOperator(H.λ, _liouvillian(H.L, Id, assume_hermitian))
+_liouvillian(H::AddedOperator, Id::AbstractMatrix, assume_hermitian::Val) =
+    AddedOperator(map(op -> _liouvillian(op, Id, assume_hermitian), H.ops))
 
 # intrinsic lindblad_dissipator
 function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
     Od_O = O' * O
-    return _sprepost(O, O') + _complex_float_type(O)(-0.5 + 0.0im) * (_spre(Od_O, Id) + _spost(Od_O, Id))
+    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
 end
 function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
     _O = O.A
@@ -211,25 +198,17 @@ function liouvillian(
     return L
 end
 
-liouvillian(
-    H::Nothing,
-    c_ops::Union{AbstractVector,Tuple},
-    Id_cache::Diagonal = I(prod(c_ops[1].dims));
-    assume_hermitian::Union{Bool,Val} = Val(true),
-) = _sum_lindblad_dissipators(c_ops, Id_cache)
+liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal = I(prod(c_ops[1].dims)); kwargs...) =
+    _sum_lindblad_dissipators(c_ops, Id_cache)
 
-liouvillian(H::Nothing, c_ops::Nothing; assume_hermitian::Union{Bool,Val} = Val(true)) = 0
+liouvillian(H::Nothing, c_ops::Nothing; kwargs...) = 0
 
 liouvillian(
     H::AbstractQuantumObject{Operator},
     Id_cache::Diagonal = I(prod(H.dimensions));
     assume_hermitian::Union{Bool,Val} = Val(true),
 ) = get_typename_wrapper(H)(_liouvillian(H.data, Id_cache, makeVal(assume_hermitian)), SuperOperator(), H.dimensions)
 
-liouvillian(
-    H::AbstractQuantumObject{SuperOperator},
-    Id_cache::Diagonal;
-    assume_hermitian::Union{Bool,Val} = Val(true),
-) = H
+liouvillian(H::AbstractQuantumObject{SuperOperator}, Id_cache::Diagonal; kwargs...) = H
 
 _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal) = sum(op -> lindblad_dissipator(op, Id_cache), c_ops)

src/time_evolution/time_evolution_dynamical.jl

@@ -571,8 +571,9 @@ function _DSF_mcsolve_Affect!(integrator)
     c_ops0_herm .= map(op -> op' * op, c_ops0)
 
     H_nh = convert(eltype(ψt), 0.5im) * sum(c_ops0_herm)
-    # By doing this, we are assuming that the system is time-independent and f is a MatrixOperator
-    copyto!(integrator.f.f.A, lmul!(-1im, H(op_l2, dsf_params).data - H_nh))
+    # By doing this, we are assuming that the system is time-independent and f is a ScaledOperator
+    # of the form -1im * (H - H_nh)
+    copyto!(integrator.f.f.L, H(op_l2, dsf_params).data - H_nh)
     return u_modified!(integrator, true)
 end
 

src/utilities.jl

@@ -175,7 +175,6 @@ end
 
 # functions for getting Float or Complex element type
 _float_type(::AbstractArray{T}) where {T<:Number} = _float_type(T)
-_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _float_type(::Type{Int32}) = Float32
 _float_type(::Type{Int64}) = Float64
 _float_type(::Type{Float32}) = Float32
@@ -187,7 +186,6 @@ _float_type(::Type{Complex{Float64}}) = Float64
 _float_type(::Type{Complex{T}}) where {T<:Real} = T
 _float_type(T::Type{<:Real}) = T # Allow other untracked Real types, like ForwardDiff.Dual
 _complex_float_type(::AbstractArray{T}) where {T<:Number} = _complex_float_type(T)
-_complex_float_type(::AbstractSciMLOperator{T}) where {T<:Number} = _complex_float_type(T)
 _complex_float_type(::Type{Int32}) = ComplexF32
 _complex_float_type(::Type{Int64}) = ComplexF64
 _complex_float_type(::Type{Float32}) = ComplexF32

test/core-test/time_evolution.jl

@@ -62,6 +62,7 @@ end
 
 @testitem "sesolve" setup=[TESetup] begin
     using SciMLOperators
+    import SciMLOperators: ScaledOperator
 
     # Get parameters from TESetup to simplify the code
     H = TESetup.H
@@ -82,7 +83,7 @@ end
     amp_rabi = TESetup.g^2 / Ω_rabi^2
     ##
 
-    @test prob.prob.f.f isa MatrixOperator
+    @test prob.prob.f.f isa ScaledOperator
     @test sum(abs.(sol.expect[1, :] .- amp_rabi .* sin.(Ω_rabi * tlist) .^ 2)) / length(tlist) < 0.1
     @test length(sol.times) == length(tlist)
     @test length(sol.times_states) == 1
@@ -390,6 +391,7 @@ end
 
 @testitem "mcsolve" setup=[TESetup] begin
     using SciMLOperators
+    import SciMLOperators: ScaledOperator
     using Statistics
 
     # Get parameters from TESetup to simplify the code
@@ -432,7 +434,7 @@ end
     expect_mc_states_mean = expect.(Ref(e_ops[1]), average_states(sol_mc_states))
     expect_mc_states_mean2 = expect.(Ref(e_ops[1]), average_states(sol_mc_states2))
 
-    @test prob_mc.prob.f.f isa MatrixOperator
+    @test prob_mc.prob.f.f isa ScaledOperator
     @test sum(abs, sol_mc.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, sol_mc2.expect .- sol_me.expect) / length(tlist) < 0.1
     @test sum(abs, average_expect(sol_mc3) .- sol_me.expect) / length(tlist) < 0.1