From d07350741d808554e8dadf0ea7f7bf2825e12023 Mon Sep 17 00:00:00 2001 From: Alberto Mercurio Date: Tue, 4 Nov 2025 20:11:45 +0100 Subject: [PATCH] Adapt to new version of SciMLOperators.jl --- src/qobj/superoperators.jl | 59 ++++++------------- .../time_evolution_dynamical.jl | 5 +- src/utilities.jl | 2 - test/core-test/time_evolution.jl | 6 +- 4 files changed, 26 insertions(+), 46 deletions(-) diff --git a/src/qobj/superoperators.jl b/src/qobj/superoperators.jl index a649be7cf..5b8765d03 100644 --- a/src/qobj/superoperators.jl +++ b/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,14 +198,10 @@ 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}, @@ -226,10 +209,6 @@ liouvillian( 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) diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl index 305822611..79e392f72 100644 --- a/src/time_evolution/time_evolution_dynamical.jl +++ b/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 diff --git a/src/utilities.jl b/src/utilities.jl index 5bbf41064..6d647c08b 100644 --- a/src/utilities.jl +++ b/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 diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl index 4feba8b2d..6d67495e6 100644 --- a/test/core-test/time_evolution.jl +++ b/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