Improve efficiency of bloch_redfield_tensor

Author: albertomercurio

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
 

test/core-test/brmesolve.jl

@@ -5,8 +5,8 @@
     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 +27,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 +38,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 +76,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 +112,4 @@ end
 
         @test all(me.expect .== brme.expect)
     end
-end;
+end