diff --git a/CHANGELOG.md b/CHANGELOG.md index ac58b3409..71fe85710 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main) +- Improve efficiency of `bloch_redfield_tensor` by avoiding unnecessary conversions. ([#509]) + ## [v0.33.0] Release date: 2025-07-22 @@ -267,3 +269,4 @@ Release date: 2024-11-13 [#504]: https://github.com/qutip/QuantumToolbox.jl/issues/504 [#506]: https://github.com/qutip/QuantumToolbox.jl/issues/506 [#507]: https://github.com/qutip/QuantumToolbox.jl/issues/507 +[#509]: https://github.com/qutip/QuantumToolbox.jl/issues/509 diff --git a/docs/src/users_guide/settings.md b/docs/src/users_guide/settings.md index be8d914e9..c2c530223 100644 --- a/docs/src/users_guide/settings.md +++ b/docs/src/users_guide/settings.md @@ -13,7 +13,8 @@ Here, we list out each setting along with the specific functions that will use i - `tidyup_tol::Float64 = 1e-14` : tolerance for [`tidyup`](@ref) and [`tidyup!`](@ref). - `auto_tidyup::Bool = true` : Automatically tidyup during the following situations: - * Solving for eigenstates, including [`eigenstates`](@ref), [`eigsolve`](@ref), and [`eigsolve_al`](@ref). + * Solving for eigenstates, including [`eigenstates`](@ref), [`eigsolve`](@ref), [`eigsolve_al`](@ref) + * Creating [`bloch_redfield_tensor`](@ref) or [`brterm`](@ref), and solving [`brmesolve`](@ref). - (to be announced) ## Change default settings diff --git a/docs/src/users_guide/time_evolution/brmesolve.md b/docs/src/users_guide/time_evolution/brmesolve.md index c7a46e612..681463d60 100644 --- a/docs/src/users_guide/time_evolution/brmesolve.md +++ b/docs/src/users_guide/time_evolution/brmesolve.md @@ -162,7 +162,7 @@ H = -Δ/2.0 * sigmax() - ε0/2 * sigmaz() ohmic_spectrum(ω) = (ω == 0.0) ? γ1 : γ1 / 2 * (ω / (2 * π)) * (ω > 0.0) -R, U = bloch_redfield_tensor(H, [(sigmax(), ohmic_spectrum)]) +R, U = bloch_redfield_tensor(H, ((sigmax(), ohmic_spectrum), )) R ``` @@ -183,6 +183,8 @@ H = - Δ/2.0 * sigmax() - ϵ0/2.0 * sigmaz() ohmic_spectrum(ω) = (ω == 0.0) ? γ1 : γ1 / 2 * (ω / (2 * π)) * (ω > 0.0) a_ops = ((sigmax(), ohmic_spectrum),) +R = bloch_redfield_tensor(H, a_ops; fock_basis = Val(true)) + e_ops = [sigmax(), sigmay(), sigmaz()] # same initial random ket state in QuTiP doc @@ -192,7 +194,7 @@ e_ops = [sigmax(), sigmay(), sigmaz()] ]) tlist = LinRange(0, 15.0, 1000) -sol = brmesolve(H, ψ0, tlist, a_ops, e_ops=e_ops) +sol = mesolve(R, ψ0, tlist, e_ops=e_ops) expt_list = real(sol.expect) # plot the evolution of state on Bloch sphere diff --git a/src/time_evolution/brmesolve.jl b/src/time_evolution/brmesolve.jl index 85ffb2594..1ef3720aa 100644 --- a/src/time_evolution/brmesolve.jl +++ b/src/time_evolution/brmesolve.jl @@ -37,18 +37,28 @@ function bloch_redfield_tensor( U = QuantumObject(rst.vectors, Operator(), H.dimensions) sec_cutoff = float(sec_cutoff) - # in fock basis - R0 = liouvillian(H, c_ops) + H_new = getVal(fock_basis) ? H : QuantumObject(Diagonal(rst.values), Operator(), H.dimensions) + c_ops_new = isnothing(c_ops) ? nothing : map(x -> getVal(fock_basis) ? x : U' * x * U, c_ops) + L0 = liouvillian(H_new, c_ops_new) - # set fock_basis=Val(false) and change basis together at the end - R1 = 0 - isempty(a_ops) || (R1 += mapreduce(x -> _brterm(rst, x[1], x[2], sec_cutoff, Val(false)), +, a_ops)) + # Check whether we can rotate the terms to the eigenbasis directly in the Hamiltonian space + fock_basis_hamiltonian = getVal(fock_basis) && sec_cutoff == -1 - SU = sprepost(U, U') # transformation matrix from eigen basis back to fock basis + R = isempty(a_ops) ? 0 : sum(x -> _brterm(rst, x[1], x[2], sec_cutoff, fock_basis_hamiltonian), a_ops) + + # If in fock basis, we need to transform the terms back to the fock basis + # Note: we can transform the terms in the Hamiltonian space only if sec_cutoff is -1 + # otherwise, we need to use the SU superoperator below to transform the entire Liouvillian + # at the end, due to the action of M_cut if getVal(fock_basis) - return R0 + SU * R1 * SU' + if fock_basis_hamiltonian + return L0 + R # Already rotated in the Hamiltonian space + else + SU = sprepost(U, U') + return L0 + SU * R * SU' + end else - return SU' * R0 * SU + R1, U + return L0 + R, U end end @@ -86,11 +96,21 @@ function brterm( fock_basis::Union{Bool,Val} = Val(false), ) rst = eigenstates(H) - term = _brterm(rst, a_op, spectra, sec_cutoff, makeVal(fock_basis)) + U = QuantumObject(rst.vectors, Operator(), H.dimensions) + + # Check whether we can rotate the terms to the eigenbasis directly in the Hamiltonian space + fock_basis_hamiltonian = getVal(fock_basis) && sec_cutoff == -1 + + term = _brterm(rst, a_op, spectra, sec_cutoff, fock_basis_hamiltonian) if getVal(fock_basis) - return term + if fock_basis_hamiltonian + return term # Already rotated in the Hamiltonian space + else + SU = sprepost(U, U') + return SU * term * SU' + end else - return term, Qobj(rst.vectors, Operator(), rst.dimensions) + return term, U end end @@ -99,7 +119,7 @@ function _brterm( a_op::T, spectra::F, sec_cutoff::Real, - fock_basis::Union{Val{true},Val{false}}, + fock_basis_hamiltonian::Union{Bool,Val}, ) where {T<:QuantumObject{Operator},F<:Function} _check_br_spectra(spectra) @@ -110,9 +130,11 @@ function _brterm( spectrum = spectra.(skew) A_mat = U' * a_op.data * U + A_mat_spec = A_mat .* spectrum + A_mat_spec_t = A_mat .* transpose(spectrum) - ac_term = (A_mat .* spectrum) * A_mat - bd_term = A_mat * (A_mat .* trans(spectrum)) + ac_term = A_mat_spec * A_mat + bd_term = A_mat * A_mat_spec_t if sec_cutoff != -1 m_cut = similar(skew) @@ -124,20 +146,29 @@ function _brterm( M_cut = @. abs(vec_skew - vec_skew') < sec_cutoff end - out = - 0.5 * ( - + _sprepost(A_mat .* trans(spectrum), A_mat) + _sprepost(A_mat, A_mat .* spectrum) - _spost(ac_term, Id) - - _spre(bd_term, Id) - ) + # Rotate the terms to the eigenbasis if possible + if getVal(fock_basis_hamiltonian) + A_mat = U * A_mat * U' + A_mat_spec = U * A_mat_spec * U' + A_mat_spec_t = U * A_mat_spec_t * U' + ac_term = U * ac_term * U' + bd_term = U * bd_term * U' + end + + # Remove small values before passing in the Liouville space + if settings.auto_tidyup + tidyup!(A_mat) + tidyup!(A_mat_spec) + tidyup!(A_mat_spec_t) + tidyup!(ac_term) + tidyup!(bd_term) + end + + out = (_sprepost(A_mat_spec_t, A_mat) + _sprepost(A_mat, A_mat_spec) - _spost(ac_term, Id) - _spre(bd_term, Id)) / 2 (sec_cutoff != -1) && (out .*= M_cut) - if getVal(fock_basis) - SU = _sprepost(U, U') - return QuantumObject(SU * out * SU', SuperOperator(), rst.dimensions) - else - return QuantumObject(out, SuperOperator(), rst.dimensions) - end + return QuantumObject(out, SuperOperator(), rst.dimensions) end @doc raw""" diff --git a/test/core-test/brmesolve.jl b/test/core-test/brmesolve.jl index 63b6ba02c..2ac87a0b1 100644 --- a/test/core-test/brmesolve.jl +++ b/test/core-test/brmesolve.jl @@ -5,8 +5,9 @@ A_op = a+a' spectra(x) = (x>0) * 0.5 for sec_cutoff in [0, 0.1, 1, 3, -1] - R = bloch_redfield_tensor(H, [(A_op, spectra)], [a^2], sec_cutoff = sec_cutoff, fock_basis = true) - R_eig, evecs = bloch_redfield_tensor(H, [(A_op, spectra)], [a^2], sec_cutoff = sec_cutoff, fock_basis = false) + R = bloch_redfield_tensor(H, ((A_op, spectra),), [a^2], sec_cutoff = sec_cutoff, fock_basis = Val(true)) + R_eig, evecs = + bloch_redfield_tensor(H, ((A_op, spectra),), [a^2], sec_cutoff = sec_cutoff, fock_basis = Val(false)) @test isa(R, QuantumObject) @test isa(R_eig, QuantumObject) @test isa(evecs, QuantumObject) @@ -27,7 +28,7 @@ end # this test applies for limited cutoff lindblad = lindblad_dissipator(a) - computation = brterm(H, A_op, spectra, sec_cutoff = 1.5, fock_basis = true) + computation = brterm(H, A_op, spectra, sec_cutoff = 1.5, fock_basis = Val(true)) @test isapprox(lindblad, computation, atol = 1e-15) end @@ -38,8 +39,8 @@ end A_op = a+a' spectra(x) = x>0 for sec_cutoff in [0, 0.1, 1, 3, -1] - R = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = true) - R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = false) + R = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = Val(true)) + R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = sec_cutoff, fock_basis = Val(false)) @test isa(R, QuantumObject) @test isa(R_eig, QuantumObject) @test isa(evecs, QuantumObject) @@ -76,8 +77,8 @@ end a = destroy(N) + destroy(N)^2/2 A_op = a+a' for spectra in spectra_list - R = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = true) - R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = false) + R = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = Val(true)) + R_eig, evecs = brterm(H, A_op, spectra, sec_cutoff = 0.1, fock_basis = Val(false)) @test isa(R, QuantumObject) @test isa(R_eig, QuantumObject) @test isa(evecs, QuantumObject) @@ -112,4 +113,4 @@ end @test all(me.expect .== brme.expect) end -end; +end