QuantumOptics.jl time-dependent Hamiltonian cannot run

fqguo · April 24, 2024, 1:49pm

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...)
   @ QuantumOptics C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:228
  master(::Any, ::Operator, ::SuperOperator; fout, kwargs...)
   @ QuantumOptics C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:119
  master(::Any, ::Operator, ::AbstractOperator, ::Any; rates, Jdagger, fout, kwargs...)
   @ QuantumOptics C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:81

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

Krastanov · April 24, 2024, 2:32pm

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

This is typical markdown formatting. Consult also Please read: make it easier to help you

Krastanov · April 24, 2024, 2:33pm

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

fqguo · April 24, 2024, 2:52pm

can you help me to solve this problem?

Krastanov · April 24, 2024, 6:34pm

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)).

fqguo · April 25, 2024, 3:58am

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.

Krastanov · April 25, 2024, 10:42am

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.

Onto your code:

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

Krastanov · April 25, 2024, 10:44am

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

fqguo · April 25, 2024, 2:10pm

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)
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:398
  [2] dmaster_
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:216 [inlined]
  [3] df_
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:23 [inlined]
  [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{})
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:59
 [23] integrate(tspan::StepRangeLen{…}, df::Function, x0::Matrix{…}, state::Operator{…}, dstate::Operator{…}, fout::Function)
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\timeevolution_base.jl:14
 [24] integrate_master(tspan::StepRangeLen{…}, df::Function, rho0::Operator{…}, fout::Function; kwargs::@Kwargs{})
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:264
 [25] integrate_master
    @ C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:260 [inlined]
 [26] #master_dynamic#40
    @ C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:217 [inlined]
 [27] master_dynamic
    @ C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:211 [inlined]
 [28] master_dynamic(tspan::StepRangeLen{…}, psi0::Ket{…}, args::Function; kwargs::@Kwargs{…})
    @ QuantumOptics.timeevolution C:\Users\Simon\.julia\packages\QuantumOptics\9adUY\src\master.jl:228
 [29] top-level scope
    @ Untitled-1:76
Some type information was truncated. Use `show(err)` to see complete types.

Krastanov · April 25, 2024, 2:26pm

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.

fqguo · April 25, 2024, 3:35pm

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.

Krastanov · April 25, 2024, 3:43pm

Thanks for saying this Happy to help!