`timeevolution.master_dynamic` cannot run correctly

fqguo · June 20, 2024, 10:44am

I have another question on timeevolution.master_dynamic and timeevolution.schroedinger_dynamic. I run the same code on time-evolution system by timeevolution.master_dynamic with density operator and timeevolution.schroedinger_dynamic and ket state. However, i can obtain correct answer with timeevolution.schroedinger_dynamic, while show an error with timeevolution.master_dynamic.
Below is my program:

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms


b = NLevelBasis(dim) # Basis for a single atom

s0 = nlevelstate(b,1) # |0> state
s1 = nlevelstate(b,2) # |1> state
sr = nlevelstate(b,3) # |r> state
plus = s0 # |+> state
minus = s1 # |-> state
# plus = (s0 + s1)/sqrt(2) # |+> state
# minus = (s0 - s1)/sqrt(2) # |-> state
Identity4 = identityoperator(b)
Hzeros = 0*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)

XZ = tensor(sigma_x, sigma_z, Identity4)
ZXZ = tensor(sigma_z, sigma_x, sigma_z)
ZX = tensor(Identity4, sigma_z, sigma_x)

Omega_1 = 2*pi*3*10^6;   Omega_2 = 2*pi*3*10^6;
Delta = 20*Omega_1;      U_rr = Delta;
NN = 6;
MM = 16;
t2 = pi/NN/Omega_2;
t1 = 7*t2;
t_pump = t1 + t2;
t_decay = 3.02*10^-6;

t_cycle = 3*NN*t_pump;

gamma_r = 2*pi*0*10^3;
gammar_P = 1/(26.2*10^-9);
Omega_pr = 0.2*gammar_P;
gammar_eff = Omega_pr^2/gammar_P;
gamma_0r = 2/3*gammar_eff;
gamma_1r = 1/3*gammar_eff;

Hv = U_rr*( tensor(sr⊗dagger(sr), sr⊗dagger(sr), Identity4) + tensor(sr⊗dagger(sr), Identity4, sr⊗dagger(sr)) + tensor(Identity4, sr⊗dagger(sr), sr⊗dagger(sr)) );

function URP123(t,args=nothing)
        
    H0_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(s0) + Omega_1/2*exp(1im*Delta*t)*s0⊗dagger(sr) - (Omega_1/2)^2/Delta*(s0⊗dagger(s0) - sr⊗dagger(sr));
    H0_ome2 = 1im*Omega_2/2*sr⊗dagger(s0) - 1im*Omega_2/2*s0⊗dagger(sr);

    HP_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(plus) + Omega_1/2*exp(-1im*Delta*t)*plus⊗dagger(sr) - (Omega_1/2)^2/Delta*(plus⊗dagger(plus) - sr⊗dagger(sr));
    HP_ome2 = 1im*Omega_2/2*sr⊗dagger(plus) - 1im*Omega_2/2*plus⊗dagger(sr);

    H1_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(s1) + Omega_1/2*exp(-1im*Delta*t)*s1⊗dagger(sr) - (Omega_1/2)^2/Delta*(s1⊗dagger(s1) - sr⊗dagger(sr));
    H1_ome2 = 1im*Omega_2/2*sr⊗dagger(s1) - 1im*Omega_2/2*s1⊗dagger(sr);

    HM_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(minus) + Omega_1/2*exp(-1im*Delta*t)*minus⊗dagger(sr) - (Omega_1/2)^2/Delta*(minus⊗dagger(minus) - sr⊗dagger(sr));
    HM_ome2 = 1im*Omega_2/2*sr⊗dagger(minus) - 1im*Omega_2/2*minus⊗dagger(sr);

    URP123_ome1 = tensor(H0_ome1, Identity4, Identity4) + tensor(Identity4, HP_ome1, Identity4) + tensor(Identity4, Identity4, H1_ome1) + Hv;
    URP123_ome2 = tensor(H0_ome2, Identity4, Identity4) + tensor(Identity4, HP_ome2, Identity4) + tensor(Identity4, Identity4, H1_ome2) + Hv;

    return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1)
end

psi_0 = tensor(s0,minus,s0);
psi_f = tensor(sr,minus,s0);

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

tlist_pump = range(0, stop=t_pump, length=1000);

exp_psi0 = [];

for jj = 1:NN
    global tout, rho_XZ, dm0, psi_0

    tout, rho_XZ = timeevolution.master_dynamic(tlist_pump, dm0, URP123, []);
    append!(exp_psi0, real(expect(rho_final, rho_XZ)))

    dm0 = rho_XZ[end];
end
plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(exp_psi0)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()

which show the error

ERROR: MethodError: no method matching master_dynamic(::StepRangeLen{…}, ::Operator{…}, ::typeof(URP123), ::Vector{…})

Closest candidates are:
  master_dynamic(::Any, ::Ket, ::Any...; kwargs...)
   @ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:228
  master_dynamic(::Any, ::Operator, ::AbstractTimeDependentOperator, ::Any; kwargs...)
   @ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:220
  master_dynamic(::Any, ::Operator, ::Any; rates, fout, kwargs...)
   @ QuantumOptics C:\Users\Administrator\.julia\packages\QuantumOptics\9adUY\src\master.jl:211

Stacktrace:
 [1] top-level scope
   @ d:\fqguo\software\Onedrive_education\OneDrive - nenu.edu.cn\A级重要---第二篇论文程序\验证URP-N次演化\URP_unitary.jl:100
Some type information was truncated. Use `show(err)` to see complete types.

and another program

using QuantumOptics
using PyPlot
# elapsed_time_1 = @elapsed begin

dim = 3; # Dimension of the single atom Hilbert space
Natom = 2; # Number of atoms


b = NLevelBasis(dim) # Basis for a single atom

s0 = nlevelstate(b,1) # |0> state
s1 = nlevelstate(b,2) # |1> state
sr = nlevelstate(b,3) # |r> state
plus = s0 # |+> state
minus = s1 # |-> state
# plus = (s0 + s1)/sqrt(2) # |+> state
# minus = (s0 - s1)/sqrt(2) # |-> state
Identity4 = identityoperator(b)
Hzeros = 0*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)

XZ = tensor(sigma_x, sigma_z, Identity4)
ZXZ = tensor(sigma_z, sigma_x, sigma_z)
ZX = tensor(Identity4, sigma_z, sigma_x)

Omega_1 = 2*pi*3*10^6;   Omega_2 = 2*pi*3*10^6;
Delta = 20*Omega_1;      U_rr = Delta;
NN = 6;
MM = 16;
t2 = pi/NN/Omega_2;
t1 = 7*t2;
t_pump = t1 + t2;
t_decay = 3.02*10^-6;

t_cycle = 3*NN*t_pump;

gamma_r = 2*pi*0*10^3;
gammar_P = 1/(26.2*10^-9);
Omega_pr = 0.2*gammar_P;
gammar_eff = Omega_pr^2/gammar_P;
gamma_0r = 2/3*gammar_eff;
gamma_1r = 1/3*gammar_eff;

Hv = U_rr*( tensor(sr⊗dagger(sr), sr⊗dagger(sr), Identity4) + tensor(sr⊗dagger(sr), Identity4, sr⊗dagger(sr)) + tensor(Identity4, sr⊗dagger(sr), sr⊗dagger(sr)) );

function URP123(t,args=nothing)
        
    H0_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(s0) + Omega_1/2*exp(1im*Delta*t)*s0⊗dagger(sr) - (Omega_1/2)^2/Delta*(s0⊗dagger(s0) - sr⊗dagger(sr));
    H0_ome2 = 1im*Omega_2/2*sr⊗dagger(s0) - 1im*Omega_2/2*s0⊗dagger(sr);

    HP_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(plus) + Omega_1/2*exp(-1im*Delta*t)*plus⊗dagger(sr) - (Omega_1/2)^2/Delta*(plus⊗dagger(plus) - sr⊗dagger(sr));
    HP_ome2 = 1im*Omega_2/2*sr⊗dagger(plus) - 1im*Omega_2/2*plus⊗dagger(sr);

    H1_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(s1) + Omega_1/2*exp(-1im*Delta*t)*s1⊗dagger(sr) - (Omega_1/2)^2/Delta*(s1⊗dagger(s1) - sr⊗dagger(sr));
    H1_ome2 = 1im*Omega_2/2*sr⊗dagger(s1) - 1im*Omega_2/2*s1⊗dagger(sr);

    HM_ome1 = Omega_1/2*exp(-1im*Delta*t)*sr⊗dagger(minus) + Omega_1/2*exp(-1im*Delta*t)*minus⊗dagger(sr) - (Omega_1/2)^2/Delta*(minus⊗dagger(minus) - sr⊗dagger(sr));
    HM_ome2 = 1im*Omega_2/2*sr⊗dagger(minus) - 1im*Omega_2/2*minus⊗dagger(sr);

    URP123_ome1 = tensor(H0_ome1, Identity4, Identity4) + tensor(Identity4, HP_ome1, Identity4) + tensor(Identity4, Identity4, H1_ome1) + Hv;
    URP123_ome2 = tensor(H0_ome2, Identity4, Identity4) + tensor(Identity4, HP_ome2, Identity4) + tensor(Identity4, Identity4, H1_ome2) + Hv;

    return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1)
end

psi_0 = tensor(s0,minus,s0);
psi_f = tensor(sr,minus,s0);

dm0 = psi_0⊗dagger(psi_0);
rho_final = psi_f⊗dagger(psi_f);

tlist_pump = range(0, stop=t_pump, length=1000);

exp_psi0 = [];

for jj = 1:NN
    global tout, rho_XZ, dm0, psi_0
    tout, rho_XZ = timeevolution.schroedinger_dynamic(tlist_pump, psi_0, URP123);
    append!(exp_psi0, real(expect(rho_final, rho_XZ)))
    psi_0 = rho_XZ[end];
    
end
plt.close()
figure(figsize=(9,3))
subplot(1,2,1)
plot(exp_psi0)
xlabel(L"T")
ylabel(L"Fidelity")

tight_layout()
gcf()

Krastanov · June 20, 2024, 7:25pm

It is helpful to start by interpreting what the error says. For master_dynamic it says there is no method implemented for the set of arguments you have and then it lists a bunch of methods that are indeed implemented.

The first one has a Ket for a second argument so it is probably not what you want.

The second one has something called AbstractTimeDependentOperator for a third argument, while you seem to have your own function, so it is also probably not what you want.

The third one seems to match your signature most closely, but it is only a three-argument call.

So just from looking at the error message, the most probable issue is that you should not have that empty list at the end of your function call, rather you should have only three arguments.

The documentation seems to confirm that: API · QuantumOptics.jl – this lists a third argument that is supposed to be a function and no more arguments after that (there are optional keyword arguments, but they need to be explicitly set with a keyword).

fqguo · June 21, 2024, 1:05am

Thanks for your help. You may be right. I find " This version takes the Hamiltonian H and jump operators J as time-dependent operators. The jump operators may be <: AbstractTimeDependentOperator or other types of operator." in API · QuantumOptics.jl. However, In fact we may encounter only Hamiltonian H is time dependent while jump operators J is time-independent and even empty (corresponding to Von Neumann equation). Obviously, my goal is the von Neumann equation. According to your answer, I realized my mistake that I added extra empty list. However, as I have asked you before, in reality, the timeevolution.master can accept empty list. I am aware of the difference between these two, but I am still not clear enough on how to express time-dependent Von Neumann equation.

Krastanov · June 21, 2024, 1:30am

I am looking at this part of the documentation API · QuantumOptics.jl

You already have a correct tspan and rho0.

You need an f function for your third argument that returns a tuple (H, J, Jdagger) where H is the hamiltonian at time t and J is the list of collapse operators at that time. I guess the lists J and Jdagger will be empty in your case.

Maybe you have been looking at the other part of the documentation:

This is referring to a particular way to represent time-dependent operators much more efficiently and “lazily”. E.g. like this one API · QuantumOptics.jl

These are quite a bit more advanced. Happy to discuss them as they are a pretty interesting piece of functionality, but I think it would be better to first have the simpler implementation from above done, so that we have a good starting point.

fqguo · June 21, 2024, 2:02am

Thanks for your help very much. I may understand your answer that I confused its three uses. According your suggestion, I revise the code return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1) to return URP123_ome2 * (t >= 0 && t <= t2) + URP123_ome1 * (t > t2 && t <= t2 + t1), [], []. It shows the correct result. I appreciate your assistance！

Krastanov · June 21, 2024, 12:09pm

Happy to help! Feel free to create another post if you want to discuss performance optimization. I would suggest using @benchmark from BenchmarkTools to measure how fast your code is. It is a good idea to try to test just the URP123 function on its own first, without testing the entire solver.

fqguo · June 21, 2024, 12:35pm

Thank you for your suggestion. After I have gained some understanding of BenchmarkTools, I will try to seek your help in another post!