# QuantumOptics.jl time-dependent Hamiltonian cannot run

I cannot modefiy my programma, in which the Hamiltonian is time-dependent, then the “timeevolution.master” show some bug. I cannot run. Can someone help me?

``````using QuantumOptics

using PyPlot

using SpecialFunctions

b = NLevelBasis(4)

s0 = nlevelstate(b,1)

s1 = nlevelstate(b,2)

sr = nlevelstate(b,3)

se = nlevelstate(b,4)

eye4 = identityoperator(b)

plus = (s0 + s1)/sqrt(2)

minus = (s0 - s1)/sqrt(2)

Omega = 2*pi*1*10^6

Urr = 2*pi*100*10^6

beta = 3

delta_0 = Urr

omega_0 = Urr/beta

alpha = 2.404825

delta = alpha*omega_0

time = pi/(Omega*besselj(2*beta,alpha))

psi_0 = tensor(s0,s1)

psi_f = tensor(s0,s1)

tlist = range(0, stop=time, length=1000)

function Hami(t)

detuning = delta_0*t + (2*delta*sin((omega_0*t)/2)^2)/omega_0

H1 = Omega/2*exp(1i*(-beta*pi+alpha+pi/2))*exp(-1i*detuning)*sr*dagger(s1) + dagger((Omega/2*exp(1i*(-beta*pi+alpha+pi/2))*exp(-1i*detuning)*sr*dagger(s1)))

H2 = Omega/2*exp(1i*(-beta*pi+alpha+pi/2))*exp(-1i*detuning)*sr*dagger(s1) + dagger((Omega/2*exp(1i*(-beta*pi+alpha+pi/2))*exp(-1i*detuning)*sr*dagger(s1)))

return tensor(H1, eye4) + tensor(eye4, H2) + Urr*tensor(sr*dagger(sr), sr*dagger(sr))

end

function calc_pops(psi, t, psi0=psi_0)

fidelity = abs(dagger(psi) * psi_f)^2

return fidelity

end

tout, pops = timeevolution.master(tlist, psi_0, Hami; fout=calc_pops)

p1 = [p[1] for p=pops]

p2 = [p[2] for p=pops]

p3 = [p[3] for p=pops]

figure(figsize=(6, 3))

plot(tout, p1, label="Initial ground state")

plot(tout, p2, "--", label="Excited state")

plot(tout, p3, label="Other ground state")

axis([0, tmax, 0, 1])

legend()

gcf()
``````

the software VS code shows the error:

``````ERROR: MethodError: no method matching master(::StepRangeLen{…}, ::Operator{…}, ::typeof(Hami); fout::typeof(calc_pops))

Closest candidates are:
master(::Any, ::Ket, ::Any...; kwargs...)
master(::Any, ::Operator, ::SuperOperator; fout, kwargs...)
master(::Any, ::Operator, ::AbstractOperator, ::Any; rates, Jdagger, fout, kwargs...)

Stacktrace:
[1] master(tspan::StepRangeLen{…}, psi0::Ket{…}, args::Function; kwargs::@Kwargs{…})
[2] top-level scope
@ f:\optimze\using QuantumOptics.jl:69
Some type information was truncated. Use `show(err)` to see complete types
``````

Hi! Could you reformat your post to make it easier to follow your code. In particular, surrounding your code with backticks formats it in a way much easier to read. The backtick symbol is

```````
``````

Could you also post the trace you are getting (the error message). Usually these are extremely valuable in understanding what is happening.

can you help me to solve this problem?

Yes, happy to help! Could you please make the two changes (code formatting, using backticks, so that I can read the code without having some of the code garbled, and posting the error message you are getting so that I have a reference point on what the issue is (it would be best if the error message is in backticks as well so that it does not format weirdly)).

Thanks for your guidance. I am a beginner and not familiar with Julia and the rules of this community, so I used the wrong format. I have now revised the format to make it look clearer. This is the first program I wrote using Julia, and I really hope you can help me solve this problem.

No worries, we were all novices in the language at some point! The changes I asked for are just habits that make it much easier and quicker for other folks to help.

I actually get some errors earlier than you. The line that creates the variable `time` does not work for me because the variable `time` already exists (it is the function `Base.time`). I am not sure why we see this difference, maybe because I am running this in the REPL and not as a script, maybe something else. Anyway, not important, I just renamed `time` to `maxtime`. With that change, I now observe the same behavior as you.

Let’s read the error message you are getting:

`MethodError` means that there is no method `master` that can work on the type of arguments you have given to it. Let’s see what argument types you have given (reading what is reported in the error message):

• a range (great, that is the times over which to compute state)
• an operator (also good, presumably that is the initial density matrix)
• a function of singleton type `typeof(Hami)` - each function in julia has its own type, this also seems fine. Presumably this function calculates your time dependent Hamiltonian. The way the function is written will probably make things quite slow, but it should work. We can improve performance after we fix the current error though.

So all of these seem fine, but then in the list of suggested methods do not seem to contain a method with the signature you want. And looking at the documentation for `timeevolution.master` also does not seem very helpful in suggesting alternatives.

The issue is that you should use `timeevolution.master_dynamic` for time-dependent Hamiltonians. But the fact that the error message did not make that obvious and that the documentation of `master` did not have a link to `master_dynamic` are documentation issues in QuantumOptics.jl that should be fixed.

Anyway, after switching to `master_dynamic` there is another problem. Your function `Hami` does not have a method of the type required by `master_dynamic`. The error message is another `MethodError` saying `no method matching Hami(::Float64, ::Operator{CompositeBasis{…}, CompositeBasis{…}, Matrix{…}})`. So it seems `Hami` needs to take time and state. This can be fixed by creating an additional method that uses the already existing method:

`Hami(t,rho) = Hami(t)`

Then there is another error, complaining that `i` is not defined. This is because the imaginary unit in julia is not written as `i`, but as `im`. I will let you fix that yourself as it is in a lot of places.

In summary, just switch to `timeevolution.master_dynamic`. There are the following related topics:

• it is very helpful to understand what `MethodError` is telling you
• the documentation of `timeevolution.master` can be improved – I created an issue on github to track that (if you are interested in fixing this so future users have more useful error messages, let me know on the issue tracker – I would be happy to help you with potentially your first open source contribution)
• to read the documentation for a function type `?` in the REPL to switch it to help mode and then type the name of the object you want to learn about
• there are things we can do to make your `Hami` much faster
1 Like

Also, consider using the specific domain - quantum category for future posts – that way `QuantumOptics.jl` users will notice your post much faster.

Thank you very much! However, it seems I’m not very smart. After modifying the code based on your advice, I encountered new problems again. I still haven’t been able to solve it successfully. If you could help me once again, I would be truly grateful.

The cuurent program is fixed as

``````using QuantumOptics

using PyPlot

using SpecialFunctions

b = NLevelBasis(4)

s0 = nlevelstate(b,1)

s1 = nlevelstate(b,2)

sr = nlevelstate(b,3)

se = nlevelstate(b,4)

eye4 = identityoperator(b)

plus = (s0 + s1)/sqrt(2)

minus = (s0 - s1)/sqrt(2)

Omega = 2*pi*1*10^6

Urr = 2*pi*10*10^6

beta = 3

delta_0 = Urr

omega_0 = Urr/beta

alpha = 2.404825

delta = alpha*omega_0

maxtime = pi/abs(Omega*besselj(beta,alpha))

psi_0 = tensor(s0,s1)

psi_f = tensor(s0,s1)

tlist = range(0, stop=maxtime, length=1000)

function Hami(t)

detuning = delta_0*t + (2*delta*sin((omega_0*t)/2)^2)/omega_0

H1 = Omega/2*exp(1im*(-beta*pi+alpha+pi/2))*exp(-1im*detuning)*sr⊗dagger(s1) + (Omega/2*exp(-1im*(-beta*pi+alpha+pi/2))*exp(1im*detuning)*s1⊗dagger(sr))

H2 = Omega/2*exp(1im*(-beta*pi+alpha+pi/2))*exp(-1im*detuning)*sr⊗dagger(s1) + (Omega/2*exp(-1im*(-beta*pi+alpha+pi/2))*exp(1im*detuning)*s1⊗dagger(sr))

return tensor(H1, eye4) + tensor(eye4, H2) + Urr*tensor(sr⊗dagger(sr), sr⊗dagger(sr))

end

Hami(t,rho) = Hami(t)

function calc_pops(psi, t, psi0=psi_0)

fidelity = abs(dagger(psi) * psi_f)^2

return fidelity

end

tout, pops = timeevolution.master_dynamic(tlist, psi_0, Hami; fout=calc_pops)

p1 = [p[1] for p=pops]

p2 = [p[2] for p=pops]

p3 = [p[3] for p=pops]

figure(figsize=(6, 3))

plot(tout, p1, label="Initial ground state")

plot(tout, p2, "--", label="Excited state")

plot(tout, p3, label="Other ground state")

axis([0, tmax, 0, 1])

legend()

gcf()
``````

And the error shows as

``````ERROR: AssertionError: 3 <= length(result) <= 4
Stacktrace:
[1] dmaster_h_dynamic!(drho::Operator{…}, f::typeof(Hami), rates::Nothing, rho::Operator{…}, drho_cache::Operator{…}, t::Float64)
[2] dmaster_
[3] df_
[4] Void
@ SciMLBase C:\Users\Simon\.julia\packages\SciMLBase\aft1j\src\utils.jl:481 [inlined]
[5] (::FunctionWrappers.CallWrapper{…})(f::SciMLBase.Void{…}, arg1::Matrix{…}, arg2::Matrix{…}, arg3::SciMLBase.NullParameters, arg4::Float64)
@ FunctionWrappers C:\Users\Simon\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65
[6] macro expansion
@ C:\Users\Simon\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined]
[7] do_ccall
@ C:\Users\Simon\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined]
[8] FunctionWrapper
@ C:\Users\Simon\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined]
[9] _call
@ C:\Users\Simon\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:12 [inlined]
[10] FunctionWrappersWrapper
@ C:\Users\Simon\.julia\packages\FunctionWrappersWrappers\9XR0m\src\FunctionWrappersWrappers.jl:10 [inlined]
[11] ODEFunction
@ C:\Users\Simon\.julia\packages\SciMLBase\aft1j\src\scimlfunctions.jl:2180 [inlined]
[12] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.DP5Cache{…})
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\perform_step\low_order_rk_perform_step.jl:914
[13] __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(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::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
@ OrdinaryDiffEq C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:512
[14] __init (repeats 5 times)
@ C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:10 [inlined]
[15] #__solve#749
@ C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:5 [inlined]
[16] __solve
@ C:\Users\Simon\.julia\packages\OrdinaryDiffEq\pWCEq\src\solve.jl:1 [inlined]
[17] 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
[18] solve_call
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:567 [inlined]
[19] #solve_up#42
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:1058 [inlined]
[20] solve_up
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:1044 [inlined]
[21] #solve#40
@ DiffEqBase C:\Users\Simon\.julia\packages\DiffEqBase\eLhx9\src\solve.jl:981 [inlined]
[22] integrate(tspan::StepRangeLen{…}, df::QuantumOptics.timeevolution.var"#dmaster_#41"{…}, x0::Matrix{…}, state::Operator{…}, dstate::Operator{…}, fout::typeof(calc_pops); alg::OrdinaryDiffEq.DP5{…}, steady_state::Bool, tol::Float64, save_everystep::Bool, saveat::StepRangeLen{…}, callback::Nothing, kwargs::@Kwargs{})
[23] integrate(tspan::StepRangeLen{…}, df::Function, x0::Matrix{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
[24] integrate_master(tspan::StepRangeLen{…}, df::Function, rho0::Operator{…}, fout::Function; kwargs::@Kwargs{})
[25] integrate_master
[26] #master_dynamic#40
[27] master_dynamic
[28] master_dynamic(tspan::StepRangeLen{…}, psi0::Ket{…}, args::Function; kwargs::@Kwargs{…})
[29] top-level scope
@ Untitled-1:76
Some type information was truncated. Use `show(err)` to see complete types.
``````

You are not running any Lindbladians, only a Hamiltonian, so you need `schroedinger_dynamic` not `master_dynamic`. That too could use a better error message (reported here better error message when the Hamiltonian or Lindblad function are not of the necessary output type (currently we just get an confusing assert error) · Issue #387 · qojulia/QuantumOptics.jl · GitHub).

Then there is another problem, your `calc_pops` is backwards. The first argument should be `t`, the second argument should be the state.

`tout, pops = timeevolution.schroedinger_dynamic(tlist, psi_0, Hami)` works for me.

1 Like

Thank you sincerely for your assistance. You have guided me through understanding how to compose posts using Markdown, and most importantly, helped me complete my first Julia program. Your diligent attitude is worth emulating. I am more than happy to provide you with suggestions for improving the documentation of `timeevolution.master`, which I will be offering to you shortly. Once again, heartfelt thanks for your assistance.

1 Like

Thanks for saying this Happy to help!