Quantum Equation Addition : propmachinery.jl

qe_propmachinery.jl

##################################################
#####                                        #####
#####            propmachinery.jl            #####
#####                                        #####
##################################################

abstract QuPropagatorMethod

abstract QuEquation

immutable SchrodingerEquation{H<:QuBase.AbstractQuMatrix} <: QuEquation
    hamiltonian::H
    SchrodingerEquation(hamiltonian) = new(hamiltonian)
end

SchrodingerEquation{H<:QuBase.AbstractQuMatrix}(hamiltonian::H) = SchrodingerEquation{H}(hamiltonian)

immutable LiouvillevonNeumannEquation{H<:QuBase.AbstractQuMatrix} <: QuEquation
    liouvillian::H
    LiouvillevonNeumannEquation(superoperator) = new(superoperator)
end

LiouvillevonNeumannEquation{H<:QuBase.AbstractQuMatrix}(liouvillian::H) = LiouvillevonNeumannEquation{H}(liouvillian)

immutable QuPropagator{QPM<:QuPropagatorMethod, QVM<:Union(QuBase.AbstractQuVector,QuBase.AbstractQuMatrix)}
    eq
    init_state::QVM
    tlist
    method::QPM
    QuPropagator(eq, init_state, tlist, method) = new(eq, init_state, tlist, method)
end

QuPropagator{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(eq::SchrodingerEquation, init_state::QV, tlist, method::QPM) = QuPropagator{QPM,QV}(eq, init_state, tlist, method)

QuPropagator{QPM<:QuPropagatorMethod, QV<:QuBase.AbstractQuVector}(param::QuBase.AbstractQuMatrix, init_state::QV,  tlist, method::QPM) = QuPropagator{QPM,QV}(SchrodingerEquation(param),init_state, tlist, method)

QuPropagator{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(eq::LiouvillevonNeumannEquation, init_state::QM, tlist, method::QPM) = QuPropagator{QPM,QM}(eq, init_state, tlist, method)

QuPropagator{QPM<:QuPropagatorMethod, QM<:QuBase.AbstractQuMatrix}(param::QuBase.AbstractQuMatrix, init_state::QM,  tlist, method::QPM) = QuPropagator{QPM,QM}(LiouvillevonNeumannEquation(liouvillian_op(param)),init_state, tlist, method)

immutable QuPropagatorState
    psi
    t
    t_state
end

function Base.start(prob::QuPropagator)
    init_state = prob.init_state
    t_state = start(prob.tlist)
    t,t_state = next(prob.tlist,t_state)
    return QuPropagatorState(init_state,t,t_state)
end

function operator(qu_eq::LiouvillevonNeumannEquation)
    return qu_eq.liouvillian
end

function operator(qu_eq::SchrodingerEquation)
    return qu_eq.hamiltonian
end


function Base.next{QPM<:QuPropagatorMethod}(prob::QuPropagator{QPM}, qustate::QuPropagatorState)
    current_qustate = qustate.psi
    current_t = qustate.t
    t,t_state = next(prob.tlist, qustate.t_state)
    next_qustate = propagate(prob.method, prob.eq, t, current_t, current_qustate)
    return (t, next_qustate), QuPropagatorState(next_qustate, t, t_state)
end

Base.done(prob::QuPropagator, qustate::QuPropagatorState) = done(prob.tlist, qustate.t_state)

function liouvillian_op(h::QuBase.AbstractQuMatrix)
    # number of basis functions
    nb = size(coeffs(h), 1)
    # construct Lindblad superoperator
    SI = Array(Int,0)
    SJ = Array(Int,0)
    Lvals = Array(Complex128,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)
                    if j==n
                        lv = lv - im*h[m,i]
                    end
                    if i==m
                        lv = lv + im*h[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 QuArray(sparse(SI, SJ, Lvals, nb*nb, nb*nb))
end

function liouvillian_tensor(h::QuBase.AbstractQuMatrix)
    return -im*(tensor(eye(h), h) - tensor(h',eye(h)))
end