Error running master equation in QuantumOptics.jl

When running the following code:

using QuantumOptics
using Plots
b = SpinBasis(1//2)
γ = 1
Ω₀ = .2
δ = .2
σ₋ = sigmam(b)#DenseOperator(b,[ 0 1; 0 0])
σ₊ =sigmap(b)# DenseOperator(b, [ 0 0; 1 0])
σz = sigmaz(b)#DenseOperator(b, [1 0 ; 0 -1])
σx = sigmax(b)#DenseOperator(b,[0 1;1 0])
H = -.5*δ*σz + .5*Ω₀*σx
J = γ*σ₋
Jd = γ*σ₊
ψ₀ = spindown(b)
ρ₀ = DenseOperator(b, [1 0; 0 1])
ts = range(0,10,100)
tout, ρ = timeevolution.master(ts, ρ₀, H, [σ₋],Jdagger=[σ₊],rates = [γ])

I get the error:

MethodError: no method matching recast!(::QuantumOpticsBase.Operator{QuantumInterface.SpinBasis{1//2, Int64}, QuantumInterface.SpinBasis{1//2, Int64}, Matrix{Int64}}, ::Matrix{Float64})

The function `recast!` exists, but no method is defined for this combination of argument types.


Closest candidates are:

  recast!(::QuantumOpticsBase.Operator{B, B, T}, ::Union{SubArray, Vector}) where {B, T}

   @ QuantumOptics ~/.julia/packages/QuantumOptics/ItxYh/src/stochastic_master.jl:185

  recast!(::QuantumOpticsBase.Operator{B, B, T}, ::T) where {B, T}

   @ QuantumOptics ~/.julia/packages/QuantumOptics/ItxYh/src/master.jl:263

  recast!(::T, ::QuantumOpticsBase.Operator{B, B, T}) where {B, T}

   @ QuantumOptics ~/.julia/packages/QuantumOptics/ItxYh/src/master.jl:266

  ...
Stack trace

Here is what happened, the most recent locations are first:

    df_(dx::Matrix{…}, x::Matrix{…}, p::SciMLBase.NullParameters, t::Float64) ...show types...
    from 
    QuantumOptics → timeevolution_base.jl:22
     

Void(::Matrix{…}, ::Vararg{…}) ...show types...
from 
SciMLBase → utils.jl:486
CallWrapper(f::SciMLBase.Void{…}, arg1::Matrix{…}, arg2::Matrix{…}, arg3::SciMLBase.NullParameters, arg4::Float64) ...show types...
from 
FunctionWrappers → FunctionWrappers.jl:65
macro expansion
from 
FunctionWrappers.jl:137
do_ccall
from 
FunctionWrappers.jl:125
FunctionWrapper
from 
FunctionWrappers.jl:144
_call
from 
FunctionWrappersWrappers.jl:12
FunctionWrappersWrapper
from 
FunctionWrappersWrappers.jl:10
ODEFunction
from 
scimlfunctions.jl:2355
initialize!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqLowOrderRK.DP5Cache{…}) ...show types...
from 
OrdinaryDiffEqLowOrderRK → low_order_rk_perform_step.jl:630
#__init#76(prob::SciMLBase.ODEProblem{…}, alg::OrdinaryDiffEqLowOrderRK.DP5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::SciMLBase.CallbackSet{…}, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Float64, reltol::Float64, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{}) ...show types...
from 
OrdinaryDiffEqCore → solve.jl:528
__init
from 
solve.jl:11
#__solve#75
from 
solve.jl:6
__solve
from 
solve.jl:1
#solve_call#44(_prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEqLowOrderRK.DP5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…}) ...show types...
from 
DiffEqBase → solve.jl:612
#solve_up#53(prob::SciMLBase.ODEProblem{…}, sensealg::Nothing, u0::Matrix{…}, p::SciMLBase.NullParameters, args::OrdinaryDiffEqLowOrderRK.DP5{…}; kwargs::@Kwargs{…}) ...show types...
from 
DiffEqBase → solve.jl:1092
solve_up
from 
solve.jl:1078
#solve#51(prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEqLowOrderRK.DP5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…}) ...show types...
from 
DiffEqBase → solve.jl:1015
#integrate#1(tspan::StepRangeLen{…}, df::QuantumOptics.timeevolution.var"#dmaster_nh_#24"{…}, x0::Matrix{…}, state::QuantumOpticsBase.Operator{…}, dstate::QuantumOpticsBase.Operator{…}, fout::QuantumOptics.timeevolution.var"#fout#7"; alg::OrdinaryDiffEqLowOrderRK.DP5{…}, steady_state::Bool, tol::Float64, save_everystep::Bool, saveat::StepRangeLen{…}, callback::Nothing, kwargs::@Kwargs{}) ...show types...
from 
QuantumOptics → timeevolution_base.jl:65
 
integrate
from 
timeevolution_base.jl:14
#integrate#6
from 
timeevolution_base.jl:81
integrate
from 
timeevolution_base.jl:76
#integrate_master#42
from 
master.jl:272
integrate_master
from 
master.jl:268
#master#22(tspan::StepRangeLen{…}, rho0::QuantumOpticsBase.Operator{…}, H::QuantumOpticsBase.Operator{…}, J::Vector{…}; rates::Vector{…}, Jdagger::Vector{…}, fout::Nothing, kwargs::@Kwargs{}) ...show types...
from 
QuantumOptics → master.jl:102
 
from 
This cell: line 1

tout, ρ = timeevolution.master(ts, ρ₀, H, [σ₋],Jdagger=[σ₊],rates = [γ])

As far as I can tell this code is correct, but I haven’t used the Lindblad master equation in Julia before so there could be a mistake somewhere that I’m missing. I’d greatly appreciate any help debugging this. I already tried using DenseOperators instead of Pauli matrices as you can see from my comments, and I also tried switching from a spin 1/2 to a 2 level system and all those approaches gave me the same errors.

Hi @DominicDams, it looks like the issue is that ρ₀ is a DenseOperator wrapped around an array of type Matrix{Int64} rather than Matrix{ComplexF64}. Try setting

ρ₀ = DenseOperator(b, ComplexF64[1 0; 0 1])

and see if that works.

On a related and slightly pedantic note: the comments here are incorrect. If you run any of these lines of code in the REPL, you’ll see that you are creating an object of type Operator wrapped around a matrix with complex entries.

julia> σ₋ = sigmam(b)
Operator(dim=2x2)
  basis: Spin(1/2)sparse([2], [1], ComplexF64[1.0 + 0.0im], 2, 2)

Thank you so much! I wonder if it’s worth trying to convert denseoperators automatically to matricies over complex numbers.

I’m sure the QuantumOptics maintainers have thoughts on this, but IMO it’s probably not worth it. As discussed in the documentation (Operators · QuantumOptics.jl), the data field of Operator supports any type that follows the Julia array interface. This enables composability with other packages that have their own array types, such as StaticArrays.jl and CUDA.jl. Automatically converting the data field of Operator to a matrix over complex numbers would pretty much remove this feature.