Merge branch 'master' into generic-operator

Project.toml

@@ -20,7 +20,7 @@ WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
 QuantumOpticsBase = "0.2"
 Arpack = "0.4"
 DiffEqCallbacks = "2"
-FFTW = "1.2"
+FFTW = "1"
 OrdinaryDiffEq = "5"
 RecursiveArrayTools = "2"
 Reexport = "0.2"

src/bloch_redfield_master.jl

@@ -23,36 +23,41 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 
 
     #Define function for transforming to Hamiltonian eigenbasis
-    function to_Heb(op)
+    function to_Heb(op, U)
         #Copy oper
         oper = copy(op)
         #Transform underlying array
-        oper.data = inv(transf_mat) * oper.data * transf_mat
+        oper.data = inv(U) * oper.data * U
         return oper
     end
 
     N = length(H_evals)
     K = length(a_ops)
 
-    #If only Lindblad collapse terms
+    # Calculate Liouvillian for Lindblad terms (unitary part + dissipation from J (if given)):
+    Heb = to_Heb(H, transf_mat)
+    #Use annon function
+    f = (x->to_Heb(x, transf_mat))
+    L = liouvillian(Heb, f.(J))
+    L = sparse(L)
+
+    #If only Lindblad collapse terms (no a_ops given)
     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{Float64}}(undef, N, N, K)
     for k in 1:K
-        A[:, :, k] = to_Heb(a_ops[k][1]).data
+        A[:, :, k] = to_Heb(a_ops[k][1], transf_mat).data
     end
 
     # Trasition frequencies between eigenstates
     W = transpose(H_evals) .- H_evals
 
     #Array for spectral functions evaluated at transition frequencies
-    Jw = Array{Complex}(undef, N, N, K)
+    Jw = Array{Complex{Float64}}(undef, N, N, K)
+    # Jw = zeros(Complex{Float64}, N, N, K)
     for k in 1:K
        # do explicit loops here
        for n in 1:N
@@ -74,23 +79,17 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
         Iabs[I, 2:3] = [indices[I].I...]
     end
 
-    # 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[]
-    cols = Int[]
-    data = Complex[]
+    # ALTERNATIVE DENSE METHOD - Main Bloch-Redfield operators part
+    data = zeros(ComplexF64, N^2, N^2)
     Is = view(Iabs, :, 1)
     As = view(Iabs, :, 2)
     Bs = view(Iabs, :, 3)
 
     for (I, a, b) in zip(Is, As, Bs)
 
         if use_secular
-            Jcds = Int.(zeros(size(Iabs)))
+            Jcds = zeros(Int, size(Iabs))
             for (row, (I2, a2, b2)) in enumerate(zip(Is, As, Bs))
                 if abs.(W[a, b] - W[a2, b2]) < dw_min * secular_cutoff
                     Jcds[row, :] = [I2 a2 b2]
@@ -114,36 +113,24 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 
         for (J, c, d) in zip(Js, Cs, Ds)
 
-            elem = 0.5 * sum(view(A, c, a, :) .* view(A, b, d, :) .* (view(Jw, c, a, :) + view(Jw, d, b, :) ))
+            sum!(   view(data, I, J),   view(A, c, a, :) .* view(A, b, d, :) .* (view(Jw, c, a, :) + view(Jw, d, b, :) )  )
+            # data[I, J] = 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, :, :) ))
+                data[I, J] -= sum( view(A, a, :, :) .* view(A, c, :, :) .* view(Jw, c, :, :) )
             end
 
             if a == c
-            elem -= 0.5 * sum(transpose(view(A, d, :, :)) .* transpose(view(A, b, :, :)) .* transpose(view(Jw, d, :, :) ))
-            end
-
-            #Add element
-            if abs(elem) != 0.0
-                push!(rows, I)
-                push!(cols, J)
-                push!(data, elem)
+                data[I, J] -= sum( view(A, d, :, :) .* view(A, b, :, :) .* view(Jw, d, :, :) )
             end
-        end
-    end
 
+            # data[I, J] *= 0.5
 
-    """
-    Need to be careful here since sparse function is happy to create rectangular arrays but we dont want this behaviour so need to make square explicitly.
-    """
-    if !any(rows .== N*N) || !any(cols .== N*N) #Check if row or column list has (N*N)th element, if not then add one
-        push!(rows, N*N)
-        push!(cols, N*N)
-        push!(data, 0.0)
+        end
     end
 
-    R = sparse(rows, cols, data) #Careful with rows/cols ordering...
+    data *= 0.5
+    R = sparse(data) # Removes any zero values and converts to sparse array
 
     #Add Bloch-Redfield part to Linblad Liouvillian calculated earlier
     L.data = L.data + R
@@ -153,6 +140,7 @@ function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secu
 end #Function
 
 
+
 """
     timeevolution.master_bloch_redfield(tspan, rho0, R, H; <keyword arguments>)
 

src/mcwf.jl

@@ -246,6 +246,7 @@ function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
                         display_beforeevent=false, display_afterevent=false,
                         display_jumps=false,
                         save_everystep=false, callback=nothing,
+                        saveat=tspan,
                         alg=OrdinaryDiffEq.DP5(),
                         kwargs...) where T
 

src/stochastic_base.jl

@@ -15,7 +15,7 @@ Integrate using StochasticDiffEq
 """
 function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
             state::T, dstate::T, fout::Function, n::Int;
-            save_everystep = false, callback=nothing,
+            save_everystep = false, callback=nothing, saveat=tspan,
             alg::StochasticDiffEq.StochasticDiffEqAlgorithm=StochasticDiffEq.EM(),
             noise_rate_prototype = nothing,
             noise_prototype_classical = nothing,
@@ -60,7 +60,7 @@ function integrate_stoch(tspan::Vector, df::Function, dg::Function, x0::Vector,
 
     out = DiffEqCallbacks.SavedValues(eltype(tspan),out_type)
 
-    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan,
+    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=saveat,
                                          save_everystep=save_everystep,
                                          save_start = false)
 

src/timeevolution_base.jl

@@ -14,7 +14,7 @@ Integrate using OrdinaryDiffEq
 function integrate(tspan, df::Function, x0::X,
             state::T, dstate::T, fout::Function;
             alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
-            steady_state = false, tol = 1e-3, save_everystep = false,
+            steady_state = false, tol = 1e-3, save_everystep = false, saveat=tspan,
             callback = nothing, kwargs...) where {T,X}
 
     function df_(dx::T, x::T, p, t) where T
@@ -32,7 +32,7 @@ function integrate(tspan, df::Function, x0::X,
 
     out = DiffEqCallbacks.SavedValues(eltype(tspan),out_type)
 
-    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan,
+    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=saveat,
                                          save_everystep=save_everystep,
                                          save_start = false)