Type-stable bloch_redfield_tensor (#262)

* Made bloch_redfield_tensor function type-stable * Quick patch failing tests
qojulia · Nov 26, 2019 · 194bff1 · 194bff1
1 parent a3b24fa
commit 194bff1
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 4 deletions.
diff --git a/src/bloch_redfield_master.jl b/src/bloch_redfield_master.jl
@@ -19,7 +19,8 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 
     # use the eigenbasis
     H_evals, transf_mat = eigen(DenseOperator(H).data)
-    H_ekets = [Ket(H.basis_l, transf_mat[:, i]) for i in 1:length(H_evals)]
+    H_ekets = [Ket(H.basis_l, transf_mat[:, i]) for i in 1:length(H_evals)]::Array{Ket{typeof(H.basis_l), Array{Complex{Float64}, 1}}, 1}
+
 
     #Define function for transforming to Hamiltonian eigenbasis
     function to_Heb(op)
@@ -37,11 +38,12 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     if K==0
         Heb = to_Heb(H)
         L = liouvillian(Heb, to_Heb.(J))
+        L = sparse(L)
         return L, H_ekets
     end
 
     #Transform interaction operators to Hamiltonian eigenbasis
-    A = Array{Complex}(undef, N, N, K)
+    A = Array{Complex{Float64}}(undef, N, N, K)
     for k in 1:K
         A[:, :, k] = to_Heb(a_ops[k][1]).data
     end
@@ -75,6 +77,7 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
     # Calculate Liouvillian for Lindblad temrs (unitary part + dissipation from J (if given)):
     Heb = to_Heb(H)
     L = liouvillian(Heb, to_Heb.(J))
+    L = sparse(L)
 
     # Main Bloch-Redfield operators part
     rows = Int[]
@@ -114,7 +117,7 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
             elem = 0.5 * sum(view(A, c, a, :) .* view(A, b, d, :) .* (view(Jw, c, a, :) + view(Jw, d, b, :) ))
 
             if b == d
-            elem -= 0.5 * sum(transpose(view(A, a, :, :)) .* transpose(view(A, c, :, :)) .* transpose(view(Jw, c, :, :) ))
+                elem -= 0.5 * sum(transpose(view(A, a, :, :)) .* transpose(view(A, c, :, :)) .* transpose(view(Jw, c, :, :) ))
             end
 
             if a == c

diff --git a/test/test_timecorrelations.jl b/test/test_timecorrelations.jl
@@ -52,7 +52,7 @@ omega_sample = mod(n, 2) == 0 ? [-n/2:n/2-1;] : [-(n-1)/2:(n-1)/2;]
 omega_sample .*= 2pi/tspan[end]
 omega, S = timecorrelations.spectrum(omega_sample, H, J, op; rho_ss=ρ₀)
 
-omega2, S2 = timecorrelations.spectrum(H, J, op)
+omega2, S2 = timecorrelations.spectrum(H, J, op; tol=1e-3)
 @test length(omega2) == length(S2)
 
 omegaFFT, SFFT = timecorrelations.correlation2spectrum(tspan, exp_values)