SE, LVN, LME types and integration

prop_eq_types_3.jl

using QuBase
using Expokit
import Base: start, next, done

abstract QuPropagatorMethod

abstract QuantumEquation

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

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

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

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

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(lindbladsuperop(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.superoperator
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)

immutable QuEuler <: QuPropagatorMethod
end

function propagate(prob::QuEuler, qu_eq::QuantumEquation, t, current_t, current_qustate)
    dt = t - current_t
    return (eye(operator(qu_eq))-im*operator(qu_eq)*dt)*current_qustate
end

immutable QuExpm_Expo <: QuPropagatorMethod
    options::Dict{Symbol, Any}
end

QuExpm_Expo() = QuExpm_Expo(Dict())

function propagate(prob::QuExpm_Expo, qu_eq::QuantumEquation, t, current_t, current_qustate)
    dt = t - current_t
    next_state = Expokit.expmv(dt, -im*coeffs(operator(qu_eq)), coeffs(current_qustate), m = get(prob.options, :m, 30), tol = get(prob.options, :tol, 1e-7))
    return QuArray(next_state)
end

function lindbladsuperop(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

qv = normalize!(QuArray([0.5+0.1im, 0.2+0.2im]))
t = 0.:0.1:2*pi
lvn = LiouvillevonNeumannEquation(lindbladsuperop(sigmax))
qexpo = QuPropagator(sigmax, qv*qv', t, QuEuler())
qexpo_m = QuPropagator(lvn, qv*qv', t, QuExpm_Expo())