Liouvillian from Lindblad constructs

liouvillian_lindblad.jl

function lindblad_op(hamiltonian::QuBase.AbstractQuMatrix, collapse_ops::Vector)
    nb = size(coeffs(hamiltonian), 1)
    nlop = length(collapse_ops)
    heff = zeros(typeof(im*one(eltype(hamiltonian))), size(coeffs(hamiltonian)))
    for l=1:nlop
        heff = heff + 0.5*coeffs(collapse_ops[l])'*coeffs(collapse_ops[l])
    end
    SI = Array(Int,0)
    SJ = Array(Int,0)
    Lvals = Array(typeof(im*one(eltype(hamiltonian))),0)
    for m=1:nb
        for n=1:nb
            for i=1:nb
                for j=1:nb
                    sm = (n-1)*nb + m
                    sj = (j-1)*nb + i
                    lv = zero(Complex128)
                    for l=1:nlop
                        lv = lv + coeffs(collapse_ops[l])[m,i]*coeffs(collapse_ops[l]')[n,j]
                    end
                    if j==n
                        lv = lv - im*coeffs(hamiltonian)[m,i] - heff[m,i]
                    end
                    if i==m
                        lv = lv + im*coeffs(hamiltonian)[j,n] - heff[j,n]
                    end
                    if real(lv)!=0 || imag(lv)!=0
                        push!(SI, sm)
                        push!(SJ, sj)
                        push!(Lvals, lv)
                    end
                end
            end
        end
    end
    return QuBase.similar_type(hamiltonian)(sparse(SI, SJ, Lvals, nb*nb, nb*nb))
end

liouvillian_op(h::QuBase.AbstractQuMatrix) = lindblad_op(h, [])