From b4f126338569338036a1c9d04214df5ffa396a48 Mon Sep 17 00:00:00 2001 From: apkille Date: Sat, 10 Aug 2024 19:20:17 -0400 Subject: [PATCH 1/4] add nocache derivative methods --- src/master.jl | 8 ++++++++ src/mcwf.jl | 4 ++++ src/semiclassical.jl | 6 ++++++ 3 files changed, 18 insertions(+) diff --git a/src/master.jl b/src/master.jl index 0cf63ff6..307235cc 100644 --- a/src/master.jl +++ b/src/master.jl @@ -332,6 +332,8 @@ function dmaster_h!(drho, H, J, Jdagger, rates::AbstractMatrix, rho, drho_cache) return drho end +dmaster_h!(drho, H, J, Jdagger, rates, rho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, copy(drho)) + """ dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, drho_cache) @@ -373,6 +375,8 @@ function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::AbstractMatrix, r return drho end +dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho) = dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, copy(drho)) + """ dmaster_liouville!(drho,L,rho) @@ -410,6 +414,8 @@ function dmaster_h_dynamic!(drho, f::F, rates, rho, drho_cache, t) where {F} dmaster_h!(drho, H, J, Jdagger, rates_, rho, drho_cache) end +dmaster_h_dynamic!(drho, f::F, rates, rho, t) where {F} = dmaster_h_dynamic!(drho, f, rates, rho, copy(drho), t) + """ dmaster_nh_dynamic!(drho, f, rates, rho, drho_cache, t) @@ -433,6 +439,8 @@ function dmaster_nh_dynamic!(drho, f::F, rates, rho, drho_cache, t) where {F} dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates_, rho, drho_cache) end +dmaster_nh_dynamic!(drho, f::F, rates, rho, t) where {F} = dmaster_nh_dynamic!(drho, f, rates, rho, copy(drho), t) + function check_master(rho0, H, J, Jdagger, rates) isreducible = true # test if all operators are sparse or dense diff --git a/src/mcwf.jl b/src/mcwf.jl index 30ece69e..8a8d4ead 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -270,6 +270,8 @@ function dmcwf_h_dynamic!(dpsi, f::F, rates, psi, dpsi_cache, t) where {F} dmcwf_h!(dpsi, H, J, Jdagger, rates_, psi, dpsi_cache) end +dmcwf_h_dynamic!(dpsi, f::F, rates, psi, t) where {F} = dmcwf_h_dynamic!(dpsi, f, rates, psi, copy(dpsi), t) + """ dmcwf_nh_dynamic!(dpsi, f, psi, t) @@ -557,6 +559,8 @@ function dmcwf_h!(dpsi, H, J, Jdagger, rates::AbstractVector, psi, dpsi_cache) return dpsi end +dmcwf_h!(dpsi, H, J, Jdagger, rates, psi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, copy(dpsi)) + """ check_mcwf(psi0, H, J, Jdagger, rates) diff --git a/src/semiclassical.jl b/src/semiclassical.jl index 55505df2..9887817b 100644 --- a/src/semiclassical.jl +++ b/src/semiclassical.jl @@ -217,6 +217,9 @@ function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t return dstate end +dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t) where {F,G} + = dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate), t) + """ dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t) @@ -232,6 +235,9 @@ function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, tmp, t) return dpsi end +dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, t) where {F,G} + = dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi), t) + function jump_dynamic(rng, t, psi, fquantum::F, fclassical!::G, fjump_classical!::H, psi_new, probs_tmp, rates) where {F,G,H} result = fquantum(t, psi.quantum, psi.classical) QO_CHECKS[] && @assert 3 <= length(result) <= 4 From 3527910b5d707e5f7c7a901afd15d55939c73d11 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 11 Aug 2024 03:58:25 -0400 Subject: [PATCH 2/4] add tests and docstrings --- src/master.jl | 8 ++++++-- src/mcwf.jl | 4 +++- src/semiclassical.jl | 14 +++++++++----- test/test_semiclassical.jl | 12 ++++++++++++ test/test_timeevolution_master.jl | 25 +++++++++++++++++++++++++ test/test_timeevolution_mcwf.jl | 14 ++++++++++++++ 6 files changed, 69 insertions(+), 8 deletions(-) diff --git a/src/master.jl b/src/master.jl index 307235cc..d69c9a04 100644 --- a/src/master.jl +++ b/src/master.jl @@ -280,6 +280,7 @@ end """ dmaster_h!(drho, H, J, Jdagger, rates, rho, drho_cache) + dmaster_h!(drho, H, J, rates, rho) Update `drho` according to a master equation given in standard Lindblad form. A cached copy `drho_cache` of `drho` is used as a temporary saving step. @@ -332,10 +333,11 @@ function dmaster_h!(drho, H, J, Jdagger, rates::AbstractMatrix, rho, drho_cache) return drho end -dmaster_h!(drho, H, J, Jdagger, rates, rho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, copy(drho)) +dmaster_h!(drho, H, J, rates, rho) = dmaster_h!(drho, H, J, dagger.(J), rates, rho, copy(drho)) """ dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, drho_cache) + dmaster_nh!(drho, Hnh, J, rates, rho) Updates `drho` according to a master equation given in standard Lindblad form. The part of the Liuovillian which can be written as a commutator should be @@ -375,7 +377,7 @@ function dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates::AbstractMatrix, r return drho end -dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho) = dmaster_nh!(drho, Hnh, Hnh_dagger, J, Jdagger, rates, rho, copy(drho)) +dmaster_nh!(drho, Hnh, J, rates, rho) = dmaster_nh!(drho, Hnh, dagger(Hnh), J, dagger.(J), rates, rho, copy(drho)) """ dmaster_liouville!(drho,L,rho) @@ -393,6 +395,7 @@ end """ dmaster_h_dynamic!(drho, f, rates, rho, drho_cache, t) + dmaster_h_dynamic!(drho, f, rates, rho, t) Computes the Hamiltonian and jump operators as `H,J,Jdagger=f(t,rho)` and update `drho` according to a master equation. Optionally, rates can also be @@ -418,6 +421,7 @@ dmaster_h_dynamic!(drho, f::F, rates, rho, t) where {F} = dmaster_h_dynamic!(drh """ dmaster_nh_dynamic!(drho, f, rates, rho, drho_cache, t) + dmaster_nh_dynamic!(drho, f, rates, rho, t) Computes the non-hermitian Hamiltonian and jump operators as `Hnh,Hnh_dagger,J,Jdagger=f(t,rho)` and update `drho` according to a master diff --git a/src/mcwf.jl b/src/mcwf.jl index 8a8d4ead..99b796ad 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -251,6 +251,7 @@ end """ dmcwf_h_dynamic!(dpsi, f, rates, psi, dpsi_cache, t) + dmcwf_h_dynamic!(dpsi, f, rates, psi, t) Compute the Hamiltonian and jump operators as `H,J,Jdagger=f(t,psi)` and update `dpsi` according to a non-Hermitian Schrödinger equation. @@ -534,6 +535,7 @@ end """ dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, dpsi_cache) + dmcwf_h!(dpsi, H, J, rates, psi) Update `dpsi` according to a non-hermitian Schrödinger equation. The non-hermitian Hamiltonian is given in two parts - the hermitian part H and @@ -559,7 +561,7 @@ function dmcwf_h!(dpsi, H, J, Jdagger, rates::AbstractVector, psi, dpsi_cache) return dpsi end -dmcwf_h!(dpsi, H, J, Jdagger, rates, psi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, copy(dpsi)) +dmcwf_h!(dpsi, H, J, rates, psi) = dmcwf_h!(dpsi, H, J, dagger.(J), rates, psi, copy(dpsi)) """ check_mcwf(psi0, H, J, Jdagger, rates) diff --git a/src/semiclassical.jl b/src/semiclassical.jl index 9887817b..cfc47ec2 100644 --- a/src/semiclassical.jl +++ b/src/semiclassical.jl @@ -204,9 +204,10 @@ end """ dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, tmp, t) + dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, t) Update the semiclassical state `dstate` according to a time-dependent, -semiclassical master eqaution. +semiclassical master equation. See also: [`semiclassical.master_dynamic`](@ref) """ @@ -217,11 +218,13 @@ function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t return dstate end -dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t) where {F,G} - = dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate), t) +function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t) where {F,G} + dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate), t) +end """ dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, tmp, t) + dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, t) Update the semiclassical state `dpsi` according to a time-dependent, semiclassical and non-Hermitian Schrödinger equation (MCWF). @@ -235,8 +238,9 @@ function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, tmp, t) return dpsi end -dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, t) where {F,G} - = dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi), t) +function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, t) where {F,G} + dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi), t) +end function jump_dynamic(rng, t, psi, fquantum::F, fclassical!::G, fjump_classical!::H, psi_new, probs_tmp, rates) where {F,G,H} result = fquantum(t, psi.quantum, psi.classical) diff --git a/test/test_semiclassical.jl b/test/test_semiclassical.jl index 72198612..75d5948f 100644 --- a/test/test_semiclassical.jl +++ b/test/test_semiclassical.jl @@ -162,6 +162,18 @@ tout5, ut = semiclassical.mcwf_dynamic(T,ψ0,fquantum,fclassical,fjump_classical @test ψt1[end].classical[2] == u0[2] + 1.0 +# Test no cache derivative methods + +t0, t1 = (0.0, 1.0) + +fmaster_h_dynamic!(dstate, state, p, t) = semiclassical.dmaster_h_dynamic!(dstate, fquantum, fclassical!, nothing, state, t) +prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, ψ0, (t0, t1)) +@test_nowarn sol_master_h_dynamic! = solve(prob_master_h_dynamic!, DP5(); save_everystep=false) + +fmcwf_h_dynamic!(dpsi, psi, p, t) = semiclassical.dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, nothing, psi, t) +prob_mcwf_h_dynamic! = ODEProblem(fmcwf_h_dynamic!, ψ0, (t0, t1)) +@test_nowarn sol_mcwf_h_dynamic! = solve(prob_mcwf_h_dynamic!, DP5(); save_everystep=false) + # Test classical jump behavior before_jump = findfirst(t -> !(t∈T), tout3) after_jump = findlast(t-> !(t∈T), tout4) diff --git a/test/test_timeevolution_master.jl b/test/test_timeevolution_master.jl index b58c4c23..cf0802f5 100644 --- a/test/test_timeevolution_master.jl +++ b/test/test_timeevolution_master.jl @@ -38,6 +38,13 @@ J = [Ja, Jc] Jlazy = [LazyTensor(basis, 1, sqrt(γ)*sm), Jc] Hnh = H - 0.5im*sum([dagger(J[i])*J[i] for i=1:length(J)]) +function Ht(t, psi) + return H*exp(-(5-t)^2), J, dagger.(J) +end +function Hnht(t, psi) + Hnhtime = H*exp(-(5-t)^2) - 0.5im * sum(i' * i for i in J) + return Hnhtime, dagger(Hnhtime), J, dagger.(J) +end Hdense = dense(H) Hlazy = LazySum(Ha, Hc, Hint) @@ -109,6 +116,24 @@ tout, ρt = timeevolution.master_nh(T, ρ₀, Hnh_dense, J; reltol=1e-6) tout, ρt = timeevolution.master_nh(T, Ψ₀, Hnh_dense, Jdense; reltol=1e-6) @test tracedistance(ρt[end], ρ) < 1e-5 +# Test no cache derivative methods + +t0, t1 = (0.0, 1.0) +fmaster_h!(drho, rho, p, t) = timeevolution.dmaster_h!(drho, H, J, nothing, rho) +prob_master_h! = ODEProblem(fmaster_h!, ρ₀, (t0, t1)) +@test_nowarn sol_master_h! = solve(prob_master_h!, DP5(); save_everystep=false) + +fmaster_nh!(drho, rho, p, t) = timeevolution.dmaster_nh!(drho, Hnh, J, nothing, rho) +prob_master_nh! = ODEProblem(fmaster_nh!, ρ₀, (t0, t1)) +@test_nowarn sol_master_nh! = solve(prob_master_nh!, DP5(); save_everystep=false) + +fmaster_h_dynamic!(drho, rho, p, t) = timeevolution.dmaster_h_dynamic!(drho, Ht, nothing, rho, t) +prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, ρ₀, (t0, t1)) +@test_nowarn sol_master_h_dynamic! = solve(prob_master_h_dynamic!, DP5(); save_everystep=false) + +fmaster_nh_dynamic!(drho, rho, p, t) = timeevolution.dmaster_nh_dynamic!(drho, Hnht, nothing, rho, t) +prob_master_nh_dynamic! = ODEProblem(fmaster_nh_dynamic!, ρ₀, (t0, t1)) +@test_nowarn sol_master_nh_dynamic! = solve(prob_master_nh_dynamic!, DP5(); save_everystep=false) # Test explicit gamma vector rates_vector = [γ, κ] diff --git a/test/test_timeevolution_mcwf.jl b/test/test_timeevolution_mcwf.jl index d0465903..cbba11bf 100644 --- a/test/test_timeevolution_mcwf.jl +++ b/test/test_timeevolution_mcwf.jl @@ -115,6 +115,20 @@ tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh_dense, J; seed=UInt(1), reltol=1 tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh, J; seed=UInt(2), reltol=1e-6) @test norm(Ψt[end]-Ψ) > 0.1 +# Test no cache derivative methods + +function Ht(t, psi) + H*exp(-(5-t)^2), J, dagger.(J) +end +t0, t1 = (0.0, 1.0) + +fmcwf_h_dynamic!(dpsi, psi, p, t) = timeevolution.dmcwf_h_dynamic!(dpsi, Ht, nothing, psi, t) +prob_mcwf_h_dynamic! = ODEProblem(fmcwf_h_dynamic!, Ψ₀, (t0, t1)) +@test_nowarn sol_mcwf_h_dynamic! = solve(prob_mcwf_h_dynamic!, DP5(); save_everystep=false) + +fmcwf_h!(dpsi, psi, p, t) = timeevolution.dmcwf_h!(dpsi, H, J, nothing, psi) +prob_mcwf_h! = ODEProblem(fmcwf_h!, Ψ₀, (t0, t1)) +@test_nowarn sol_mcwf_h! = solve(prob_mcwf_h!, DP5(); save_everystep=false) # Test convergence to master solution From c230bcd43c8a0f23f27bca37f87e88c4f86d6284 Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 11 Aug 2024 04:32:28 -0400 Subject: [PATCH 3/4] add semiclassical tests --- src/semiclassical.jl | 4 ++-- test/test_semiclassical.jl | 8 +++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/semiclassical.jl b/src/semiclassical.jl index 7fcff86c..9a29256a 100644 --- a/src/semiclassical.jl +++ b/src/semiclassical.jl @@ -299,7 +299,7 @@ function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t end function dmaster_h_dynamic!(dstate, fquantum::F, fclassical!::G, rates, state, t) where {F,G} - dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate), t) + dmaster_h_dynamic!(dstate, fquantum, fclassical!, rates, state, copy(dstate.quantum), t) end """ @@ -319,7 +319,7 @@ function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, tmp, t) end function dmcwf_h_dynamic!(dpsi, fquantum::F, fclassical!::G, rates, psi, t) where {F,G} - dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi), t) + dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, rates, psi, copy(dpsi.quantum), t) end function jump_dynamic(rng, t, psi, fquantum::F, fclassical!::G, fjump_classical!::H, psi_new, probs_tmp, rates) where {F,G,H} diff --git a/test/test_semiclassical.jl b/test/test_semiclassical.jl index 4ce4f974..24bbb344 100644 --- a/test/test_semiclassical.jl +++ b/test/test_semiclassical.jl @@ -2,6 +2,7 @@ using Test using QuantumOptics using LinearAlgebra using QuantumOpticsBase: IncompatibleBases +using OrdinaryDiffEq @testset "semiclassical" begin @@ -166,12 +167,13 @@ tout5, ut = semiclassical.mcwf_dynamic(T,ψ0,fquantum,fclassical,fjump_classical # Test no cache derivative methods t0, t1 = (0.0, 1.0) +rho0 = semiclassical.State(dm(spinup(ba)⊗fockstate(bf,0)),u0) -fmaster_h_dynamic!(dstate, state, p, t) = semiclassical.dmaster_h_dynamic!(dstate, fquantum, fclassical!, nothing, state, t) -prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, ψ0, (t0, t1)) +fmaster_h_dynamic!(dstate, state, p, t) = semiclassical.dmaster_h_dynamic!(dstate, fquantum_master, fclassical, nothing, state, t) +prob_master_h_dynamic! = ODEProblem(fmaster_h_dynamic!, rho0, (t0, t1)) @test_nowarn sol_master_h_dynamic! = solve(prob_master_h_dynamic!, DP5(); save_everystep=false) -fmcwf_h_dynamic!(dpsi, psi, p, t) = semiclassical.dmcwf_h_dynamic!(dpsi, fquantum, fclassical!, nothing, psi, t) +fmcwf_h_dynamic!(dstate, state, p, t) = semiclassical.dmcwf_h_dynamic!(dstate, fquantum, fclassical, nothing, state, t) prob_mcwf_h_dynamic! = ODEProblem(fmcwf_h_dynamic!, ψ0, (t0, t1)) @test_nowarn sol_mcwf_h_dynamic! = solve(prob_mcwf_h_dynamic!, DP5(); save_everystep=false) From 2280d97e4fa78667f8c073b9e90244798181da5a Mon Sep 17 00:00:00 2001 From: apkille Date: Sun, 11 Aug 2024 05:13:07 -0400 Subject: [PATCH 4/4] add using OrdinaryDiffEq --- test/test_timeevolution_master.jl | 1 + test/test_timeevolution_mcwf.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/test/test_timeevolution_master.jl b/test/test_timeevolution_master.jl index cf0802f5..df96a264 100644 --- a/test/test_timeevolution_master.jl +++ b/test/test_timeevolution_master.jl @@ -1,6 +1,7 @@ using Test using QuantumOptics using LinearAlgebra +using OrdinaryDiffEq @testset "master" begin diff --git a/test/test_timeevolution_mcwf.jl b/test/test_timeevolution_mcwf.jl index cbba11bf..979ad208 100644 --- a/test/test_timeevolution_mcwf.jl +++ b/test/test_timeevolution_mcwf.jl @@ -1,6 +1,7 @@ using Test using QuantumOptics using Random, LinearAlgebra +using OrdinaryDiffEq @testset "mcwf" begin