When i use the package `QuantumOptics.jl`

to simulate the evolution of Liouvillian superoperator, i find it cannot run by the command `timeevolution.master_dynamic`

. Because after obtaining the superoperator, i ready to solve the propagator of this superoperator. I believe it cannot use the command `liouvillian`

, and i cannot find out which type of program my program belongs to. Thus, i need to ask if it can help me solve my difficulties.

My program is

```
using QuantumOptics
using PyPlot
using SpecialFunctions
using BenchmarkTools
using LinearAlgebra
elapsed_time_1 = @elapsed begin
b = NLevelBasis(4)
s0 = nlevelstate(b,1)
s1 = nlevelstate(b,2)
sr = nlevelstate(b,3)
se = nlevelstate(b,4)
plus = (s0 + s1)/sqrt(2)
minus = (s0 - s1)/sqrt(2)
eye4 = identityoperator(b)
eye64 = identityoperator(b⊗b⊗b)
sigma_x = s0⊗dagger(s1) + s1⊗dagger(s0)
sigma_y = -1im*s0⊗dagger(s1) + 1im*s1⊗dagger(s0)
sigma_z = s0⊗dagger(s0) - s1⊗dagger(s1)
GHZ_00 = ( tensor(plus, s0, plus) + tensor(minus, s1, minus) )/sqrt(2)
GHZ_01 = ( tensor(plus, s0, plus) - tensor(minus, s1, minus) )/sqrt(2)
GHZ_10 = ( tensor(plus, s0, minus) + tensor(minus, s1, plus) )/sqrt(2)
GHZ_11 = ( tensor(plus, s0, minus) - tensor(minus, s1, plus) )/sqrt(2)
GHZ_20 = ( tensor(plus, s1, plus) + tensor(minus, s0, minus) )/sqrt(2)
GHZ_21 = ( tensor(plus, s1, plus) - tensor(minus, s0, minus) )/sqrt(2)
GHZ_30 = ( tensor(plus, s1, minus) + tensor(minus, s0, plus) )/sqrt(2)
GHZ_31 = ( tensor(plus, s1, minus) - tensor(minus, s0, plus) )/sqrt(2)
# parameters
Omega_r = 0.2*gamma_P1; Omega_a = 0.2*gamma_P2;
Omega = 2*pi*1*10^6; Urr = 2*pi*150*10^6;
alpha = 43.184769223265;
omega_0 = Urr/14;
Delta_0 = 0*omega_0;
delta = alpha*omega_0;
m0 = Delta_0/omega_0;
m1 = (Delta_0-Urr)/omega_0;
m2 = (Delta_0-2*Urr)/omega_0;
phase = m0*pi/2-alpha-pi/2 + pi*sign(besselj(m0,alpha));
t1 = abs(pi/(Omega*besselj(m0,alpha)))
maxtime = t1
# initial and final states
psi_0 = (GHZ_00 + GHZ_01 + GHZ_10 + GHZ_11 + GHZ_20 + GHZ_21 + GHZ_30 + GHZ_31)/sqrt(8);
psi_f = GHZ_00;
psi0_l = tensor(psi_0,psi_0)
psif_l = tensor(psi_f,psi_f)
propag_0 = identityoperator(b⊗b⊗b⊗b⊗b⊗b)
tlist = range(0, stop=maxtime, length=1000)
# Hamiltonian
Hv = Urr*tensor(sr⊗dagger(sr), sr⊗dagger(sr), eye4) + Urr*tensor(sr⊗dagger(sr), eye4, sr⊗dagger(sr)) + Urr*tensor(eye4, sr⊗dagger(sr), sr⊗dagger(sr))
function Liouvillian(t,rho)
detuning = Delta_0*t + (2*delta*sin((omega_0*t)/2)^2)/omega_0;
HP = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(plus) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(plus));
H0 = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s0) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s0));
H1 = Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s1) + dagger(Omega/2*exp(-1im*phase)*exp(-1im*detuning)*sr⊗dagger(s1));
SEP_123 = tensor(H0, eye4, eye4) + tensor(eye4, HP, eye4) + tensor(eye4, eye4, H1) + Hv;
Hami = SEP_123;
Liouvillian = -1im*(tensor(Hami, eye64) - tensor(eye64, transpose(Hami)))
return Liouvillian
end
tout, rho = timeevolution.master_dynamic(tlist, propag_0, Liouvillian)
exp_psi0_l = real(expect(psif_l ⊗ dagger(psif_l), rho*(psi0_l ⊗ dagger(psi0_l))))
figure(figsize=(9,3))
# subplot(1,2,1)
# ylim([0, 2])
plot(tlist, exp_psi0_l);
xlabel(L"T")
ylabel(L"Fidelity")
tight_layout()
end
println("run time：", elapsed_time_1)
gcf()
```

And it gives the Error:

```
ERROR: MethodError: no method matching OrdinaryDiffEq.DP5Cache(::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::Diagonal{…}, ::typeof(OrdinaryDiffEq.trivial_limiter!), ::typeof(OrdinaryDiffEq.trivial_limiter!), ::Static.False)
Closest candidates are:
OrdinaryDiffEq.DP5Cache(::var"#2018#uType", ::var"#2018#uType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2019#rateType", ::var"#2018#uType", ::var"#2018#uType", ::var"#2020#uNoUnitsType", ::var"#2021#StageLimiter", ::var"#2022#StepLimiter", ::var"#2023#Thread") where {var"#2018#uType", var"#2019#rateType", var"#2020#uNoUnitsType", var"#2021#StageLimiter", var"#2022#StepLimiter", var"#2023#Thread"}
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\caches\low_order_rk_caches.jl:552
Stacktrace:
[1] alg_cache(alg::OrdinaryDiffEq.DP5{…}, u::Diagonal{…}, rate_prototype::Diagonal{…}, ::Type{…}, ::Type{…}, ::Type{…}, uprev::Diagonal{…}, uprev2::Diagonal{…}, f::SciMLBase.ODEFunction{…}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, ::Val{…})
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\caches\low_order_rk_caches.jl:600
[2] __init(prob::SciMLBase.ODEProblem{…}, alg::OrdinaryDiffEq.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(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::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:342
[3] __solve(::SciMLBase.ODEProblem{…}, ::OrdinaryDiffEq.DP5{…}; kwargs::@Kwargs{…})
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:5
[4] solve_call(_prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEq.DP5{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:609
[5] solve_up(prob::SciMLBase.ODEProblem{…}, sensealg::Nothing, u0::Diagonal{…}, p::SciMLBase.NullParameters, args::OrdinaryDiffEq.DP5{…}; kwargs::@Kwargs{…})
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:1058
[6] solve(prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEq.DP5{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:981
[7] integrate(tspan::StepRangeLen{…}, df::QuantumOptics.timeevolution.var"#dmaster_#41"{…}, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::QuantumOptics.timeevolution.var"#fout#7"; alg::OrdinaryDiffEq.DP5{…}, steady_state::Bool, tol::Float64, save_everystep::Bool, saveat::StepRangeLen{…}, callback::Nothing, kwargs::@Kwargs{})
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:59
[8] integrate(tspan::StepRangeLen{…}, df::Function, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:14
[9] #integrate#6
@ C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:75 [inlined]
[10] integrate(tspan::StepRangeLen{…}, df::Function, x0::Diagonal{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:70 [inlined]
[11] integrate_master(tspan::StepRangeLen{…}, df::Function, rho0::Operator{…}, fout::Nothing; kwargs::@Kwargs{})
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:264
[12] integrate_master
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:260 [inlined]
[13] #master_dynamic#40
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:217 [inlined]
[14] master_dynamic(tspan::StepRangeLen{…}, rho0::Operator{…}, f::Function)
@ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:211
[15] macro expansion
@ d:\OneDrive(重重重重重重重重)\OneDrive\博士期间------工作\第二篇\程序\Liouvillian方法\GHZ\Liouville 传播子\stabilzer\request_help.jl:78 [inlined]
[16] top-level scope
@ .\timing.jl:395
Some type information was truncated. Use `show(err)` to see complete types.
```

Thanks for your help!