Improve efficiency of generation of Bloch-Redfield tensor and fix documentation typos

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

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

src/time_evolution/brmesolve.jl

@@ -37,18 +37,16 @@ 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(real.(rst.values)), Operator(), H.dimensions)
+    c_ops_new = isnothing(c_ops) ? nothing : map(x -> getVal(fock_basis) ? x : U' * x * U, c_ops)
+    R = 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))
+    isempty(a_ops) || (R += sum(x -> _brterm(rst, x[1], x[2], sec_cutoff, fock_basis), a_ops))
 
-    SU = sprepost(U, U') # transformation matrix from eigen basis back to fock basis
     if getVal(fock_basis)
-        return R0 + SU * R1 * SU'
+        return R
     else
-        return SU' * R0 * SU + R1, U
+        return R, U
     end
 end
 
@@ -99,7 +97,7 @@ function _brterm(
     a_op::T,
     spectra::F,
     sec_cutoff::Real,
-    fock_basis::Union{Val{true},Val{false}},
+    fock_basis::Union{Bool,Val},
 ) where {T<:QuantumObject{Operator},F<:Function}
     _check_br_spectra(spectra)
 
@@ -114,6 +112,13 @@ function _brterm(
     ac_term = (A_mat .* spectrum) * A_mat
     bd_term = A_mat * (A_mat .* trans(spectrum))
 
+    # Remove small values before passing in the Liouville space
+    if settings.auto_tidyup
+        tidyup!(A_mat)
+        tidyup!(ac_term)
+        tidyup!(bd_term)
+    end
+
     if sec_cutoff != -1
         m_cut = similar(skew)
         map!(x -> abs(x) < sec_cutoff, m_cut, skew)

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