diff --git a/examples/mcsolve_example.jl b/examples/mcsolve_example.jl new file mode 100644 index 0000000..10e8184 --- /dev/null +++ b/examples/mcsolve_example.jl @@ -0,0 +1,60 @@ +# Quantum dynamics of a damped harmonic oscillator using the mcsolve.jl. +using QuDOS +import Winston: plot, oplot, xlabel, ylabel + +# parameters +n = 20 +gamma = 0.1 +nsteps = 50 +dt = 0.25 + +ntraj = 100 # number of trajectories + +# initialize operators etc +ad = QuDOS.creationop( n ) +h = ad*ad' + +x = QuDOS.positionop( n ) +p = QuDOS.momentumop( n ) + +cs = QuDOS.coherentstatevec(n, 1.) +rho = QuDOS.QuState(cs) + +# construct Lindblad QME from Hamiltonian h and +# Lindblad operator(s) +qme = QuDOS.LindbladQME( h, {sqrt(gamma)*ad',}) + +# propagator for QME +# lb_prop = QuDOS.QuFixedStepPropagator( qme, dt) +# or +# lb_prop = QuDOS.QuKrylovPropagator( qme ) + +# propagation loop +xex = zeros(nsteps+1) +pex = zeros(nsteps+1) + +xex[1] = real(QuDOS.expectationval(rho, x)) +pex[1] = real(QuDOS.expectationval(rho, p)) + +# for step=1:nsteps +# rho = QuDOS.propagate(lb_prop, rho, tspan=[0., dt]) +# +# xex[step+1] = real(QuDOS.expectationval(rho, x)) +# pex[step+1] = real(QuDOS.expectationval(rho, p)) +# end +tlist=[1:nsteps]*dt +measms = Array(AbstractMatrix{Complex128},2) +measms[1]=complex(x) +measms[2]=complex(p) +exvals = QuDOS.mcsolve(qme, cs, tlist, measms, ntraj) + +xex[2:end] = real(exvals[1,1,:]) +pex[2:end] = real(exvals[1,2,:]) + +# plots results using Winston +plot([0,tlist], xex, "b-", [0,tlist], pex, "r-") +oplot(dt*[0:0.5:nsteps], sqrt(2.)*exp(-(gamma)*dt*[0:0.5:nsteps]),"k--") # Should the decay rate be Gamma or Gamma/2.0? +oplot(dt*[0:0.5:nsteps],-sqrt(2.)*exp(-(gamma)*dt*[0:0.5:nsteps]),"k--") # Should the decay rate be Gamma or Gamma/2.0? + +xlabel("time") +ylabel("\\langle x\\rangle (blue) and \\langle p\\rangle (red)") diff --git a/examples/mcwf_damped_ho.jl b/examples/mcwf_damped_ho.jl new file mode 100644 index 0000000..2fa6d75 --- /dev/null +++ b/examples/mcwf_damped_ho.jl @@ -0,0 +1,57 @@ +# quantum dynamics of a damped harmonic oscillator +using QuDOS +import Winston: plot, oplot, xlabel, ylabel + +# parameters +n = 20 +gamma = 0.1 +nsteps = 50 +dt = 0.25 + +ntraj = 5 # number of trajectories + +# initialize operators etc +ad = QuDOS.creationop( n ) +h = ad*ad' + +x = QuDOS.positionop( n ) +p = QuDOS.momentumop( n ) + +cs = QuDOS.coherentstatevec(n, 1.) +rho = QuDOS.QuState(cs) + +# construct Lindblad QME from Hamiltonian h and +# Lindblad operator(s) +qme = QuDOS.LindbladQME( h, {sqrt(gamma)*ad',}) + +# propagator for QME +# lb_prop = QuDOS.QuFixedStepPropagator( qme, dt) +# or +# lb_prop = QuDOS.QuKrylovPropagator( qme ) + +# propagation loop +xex = zeros(nsteps+1) +pex = zeros(nsteps+1) + +xex[1] = real(QuDOS.expectationval(rho, x)) +pex[1] = real(QuDOS.expectationval(rho, p)) + +# for step=1:nsteps +# rho = QuDOS.propagate(lb_prop, rho, tspan=[0., dt]) +# +# xex[step+1] = real(QuDOS.expectationval(rho, x)) +# pex[step+1] = real(QuDOS.expectationval(rho, p)) +# end + +rho, exvals = QuDOS.mcwfpropagate(cs, qme, ntraj, dt, nsteps, {x, p}) + +xex[2:end] = real(exvals[:,1]) +pex[2:end] = real(exvals[:,2]) + +# plots results using Winston +plot(dt*[0:nsteps], xex, "b-", dt*[0:nsteps], pex, "r-") +oplot(dt*[0:0.5:nsteps], sqrt(2.)*exp(-(gamma/2.)*dt*[0:0.5:nsteps]),"k--") +oplot(dt*[0:0.5:nsteps],-sqrt(2.)*exp(-(gamma/2.)*dt*[0:0.5:nsteps]),"k--") + +xlabel("time") +ylabel("\\langle x\\rangle (blue) and \\langle p\\rangle (red)") diff --git a/src/Qudos.jl b/src/Qudos.jl index 74e3968..c3663c6 100644 --- a/src/Qudos.jl +++ b/src/Qudos.jl @@ -13,7 +13,7 @@ export QuFixedStepPropagator, export propagate -export LindbladQME +export LindbladQME, mcsolve export stationary, superop, eff_hamiltonian @@ -146,7 +146,7 @@ function +(rho1::QuState, rho2::QuState) error("Quantum states are not compatible and cannot be added.") end - QuState( (rho1.elem+rho2.elem),copy(rho1.subnb) ) + QuState( reshape(rho1.elem+rho2.elem,rho1.nb,rho1.nb),copy(rho1.subnb) ) end function -(vec1::QuStateVec, vec2::QuStateVec) @@ -162,7 +162,7 @@ function -(rho1::QuState, rho2::QuState) error("Quantum states are not compatible and cannot be subtracted.") end - QuState( (rho1.elem-rho2.elem), copy(rho1.subnb) ) + QuState( reshape(rho1.elem-rho2.elem,rho1.nb,rho1.nb), copy(rho1.subnb) ) end @@ -423,4 +423,10 @@ include("propagate.jl") # QME types and functionality include("qme.jl") +# MCWF method (aka quantum jump) +include("mcwf.jl") + +# mcsolve method for MCWF (similar to the QuTIP's mcsolve function) +include("mcsolve.jl") + end # module diff --git a/src/expmv.jl b/src/expmv.jl index aacdd57..b137b69 100644 --- a/src/expmv.jl +++ b/src/expmv.jl @@ -23,9 +23,9 @@ # The work resulting from EXPOKIT has been published in ACM-Transactions # on Mathematical Software, 24(1):130-156, 1998. # -expmv{T}( vec::Vector{T}, t::Real, amat::AbstractMatrix, tol::Real=1e-7, m::Int=min(30,size(amat,1))) = expmv!(copy(vec), t, amat, tol, m) +expmv{T}( vec::Vector{T}, t::Number, amat::AbstractMatrix, tol::Real=1e-7, m::Int=min(30,size(amat,1))) = expmv!(copy(vec), t, amat, tol, m) -function expmv!{T}( vec::Vector{T}, t::Real, amat::AbstractMatrix, tol::Real=1e-7, m::Int=min(30,size(amat,1))) +function expmv!{T}( vec::Vector{T}, t::Number, amat::AbstractMatrix, tol::Real=1e-7, m::Int=min(30,size(amat,1))) if size(vec,1) != size(amat,2) error("dimension mismatch") @@ -53,7 +53,7 @@ function expmv!{T}( vec::Vector{T}, t::Real, amat::AbstractMatrix, tol::Real=1e- tf = abs(t) tsgn = sign(t) - tk = zero(t) + tk = zero(tf) v = vec mx = m diff --git a/src/mcsolve.jl b/src/mcsolve.jl new file mode 100644 index 0000000..cbe34ce --- /dev/null +++ b/src/mcsolve.jl @@ -0,0 +1,163 @@ +# mcsolve method for Monte Carlo Wavefunctions. It has a similar constructure as the QuTIP mcsolve function. + +function mcsolve(H::Union(SparseMatrixCSC{Complex128,Int64},AbstractMatrix{Complex128}), + psi0::Array{Complex128,1}, + tlist::Array{Float64,1}, ops::Vector{AbstractMatrix{Complex128}}, measms::Array{AbstractMatrix{Complex128},1}, + ntraj::Union(Int64,Array{Int64,1})=500 ) + # General mcsolve method with normal data type inputs. + # + # inputs: H::Union(SparseMatrixCSC{Complex128,Int64},AbstractMatrix{Complex128}). - The system Hamiltonian. + # psi0::Array{Complex128,1}. - The initial state of the system. + # tlist::Array{Float64,1}. - Time points as a vector. + # ops::Vector{AbstractMatrix{Complex128}}. - The Lindblad operators or jump operators as a vector with matrix elements. + # measms::Array{AbstractMatrix{Complex128},1}. - The measurement operators. + # ntraj::Union(Int64,Array{Int64,1})=500. - Number of trajectories for averaging. It can be a vector with increasing integers. + # + # outputs: + # exvals=Array(Float64, lengthtraj, nmeasm, nsteps) - Expectation values of the measurement operators. + # It is given in a matrix form, where lengthtraj is the number of the recorded trajectories, nmeasm is the + # measurement operators, and nsteps is the number of the tlist time points. + # + # To be improved: 1. Unstable to perform normalization operations on states. Have tried various ways, for example, on line 67 and so + # on. Suspect it is caused by the BLAS library. + # 2. hbar=1? + # 3. Full test and documentations. + # 4. Define more input patterns to at least remove the Unions, and use parametrized types of inputs. + + #Constants. + hbar=1 # Should be defined as the real value, right? + + nsteps=length(tlist) # Number of time steps. + nops = length(ops) # Number of jump operators. + nd = size(ops[1],1) # By default, all operators should be represented as square matrices. + nmeasm=length(measms) # Number of measurement operators. + lengthtraj=length(ntraj) # Number of averaging points for given trajectory numbers. + exvals = zeros(Float64, lengthtraj, nmeasm,nsteps) # Expection values at each time step for the given measurement operators for given number of trajectories. + exvals_temp = zeros(Float64, nmeasm,nsteps) # Expection values at each time step for the given measurement operators in accumulated trajectories. + + # Get information of the effective Hamiltonian operator. + heff = H # Initialize the effective Hamiltonians. + for i=1:nops + heff = heff - 0.5im*hbar*ops[i]'*ops[i] #eff_hamiltonian(qme) + end + cop = eye(size(heff,1)) # Initialize the jump operator chosen for each jump. + + # Loop over all trajectories. + maxntraj=ntraj[end] + trajind=1 # Counting readout index in the ntraj vector. + for traj=1:maxntraj # Loop over all trajectories including all given recording points. + println("Starting trajectory $traj/$maxntraj:") + # The initial state and temple measurement records. + eps = rand() # Draw a random number from [0,1) to represent the probability that a quantum jump occurs. + psi = copy(psi0) # Initialize the quantum state vector before each evolution step. + psi1 = copy(psi0) # Initialize the quantum state vector after each step of evolution. + psi_norm = 1 # Initialize the norm of the state vector. + t0=tlist[1] # Assumed tlist has a length larger than 1. + for nm=1:nmeasm + exvals_temp[nm,1] += real(psi'*measms[nm]*psi)[1] + end + # All the rest time steps. + for step=2:nsteps # Evolve over all time steps other than the initial one. + dt = tlist[step]-t0 + if psi_norm > eps # No jump case. + # Propagate one time-step until the jumping condition satisfies. + psi1 = expmv(psi,-dt/hbar,1im*heff) # Free evolution of the state. propagate(eff_prop, psi) + # No need to normalize the state here: psi1=psi1/norm(psi1) + psi = copy(psi1) # State for the next step after the evolution time dt. + psi_norm = real(psi'*psi)[1] #norm(psi) # Norm after this step of evolution. + println("t = $(tlist[step]), norm = $psi_norm.") + else # Otherwise, jump happens. + println("Jump at t=$(tlist[step]), norm=$psi_norm, eps=$eps .") + # Calculate the probability distribution for each jump operator. + Pi = zeros(Float64,nops) # Initialize the probability of each jump. + PnS = 0. # Initialize the unnormalized accumulated (sum) probability of jumps. + for cind = 1:nops + Pi[cind] = norm(ops[cind]*psi).^2 # Probability of jump i: real(psi1'*ops[cind]'*ops[cind]*psi1) + PnS += Pi[cind] # Total probability of jumps without normalization. + end + # Judge which jump. + Pn = 0. # Initialize the normalized accumulated probability of jumps. + for cind = 1:nops + Pn += Pi[cind]/PnS + println("Pn($cind)=$Pn") + if Pn >= eps + println("Using jump operator $cind") + cop=ops[cind] + break + end + end + # New state after the jump of cop. + psi1 = cop*psi #QuStateVec( cop*psi1.elem ) + # Normalize the state for the next time step. + psi = psi1./norm(psi1) # normalize!(psi1) + psi_norm = 1.0 # Norm of the state after this step of evolution. + # Generate a new random number for the next jump judgement. + eps = rand() # draw new random number from [0,1) + end + t0=tlist[step] # Record time to get dt for the next iteration. + # Expectation values of measurement operators. + for nm = 1:nmeasm + exvals_temp[nm,step] += real(psi'*measms[nm]*psi)[1] + end + end # Time-steps + # Averaging expectation values for requested quantum trajectory numbers. + if traj==ntraj[trajind] + for j=1:nsteps, i=1:nmeasm + exvals[trajind,i,j] = copy(exvals_temp[i,j]./traj) + end + trajind += 1 + end + end # trajectories + + return exvals +end + +###################################################################################### +# Other input type patterns. +###################################################################################3## +function mcsolve(qme::LindbladQME, psi0::AbstractQuState, + tlist::Array{Float64,1}, measms::Array{AbstractArray{Complex128,2},1}, + ntraj::Array{Int,1}=[500] ) + # mcsolve method with AbstractLindbladQME and AbstractQuState types of inputs. + + H = qme.hamop + psi_init = psi0.elem + ops = QuDOS.linbladops(qme) + exvals = mcsolve(H, psi0, tlist, ops, measms, ntraj) + return exvals +end + +function mcsolve(qme::LindbladQME, psi0::AbstractQuState, + tlist::Array{Float64,1}, measms::Array{AbstractArray{Complex128,2},1}, + ntraj::Int64=500 ) + # mcsolve method with AbstractLindbladQME and AbstractQuState types of inputs. + + H = qme.hamop + psi_init = vec(full(psi0.elem)) + ops = QuDOS.linbladops(qme) + exvals = mcsolve(H, psi_init, tlist, ops, measms, ntraj) + return exvals +end + +function mcsolve(qme::LindbladQME, psi0::AbstractQuState, + tlist::Array{Float64,1}, measms::Array{AbstractArray{Complex128,2},1}, + ntraj::Array{Int64,1}=[500] ) + # mcsolve method with AbstractLindbladQME and AbstractQuState types of inputs. + + H = qme.hamop + psi_init = vec(full(psi0.elem)) + ops = QuDOS.linbladops(qme) + exvals = mcsolve(H, psi_init, tlist, ops, measms, ntraj) + return exvals +end + +function mcsolve(H::SparseMatrixCSC{Complex{Float64},Int64}, psi0::QuStateVec, + tlist::Array{Float64,1}, ops::Array{AbstractArray{Complex{Float64},2},1}, + measms::Array{AbstractArray{Complex128,2},1}, + ntraj::Int64=500 ) + # mcsolve method with AbstractQuState type of inputs. + + psi_init = vec(full(psi0.elem)) + exvals = mcsolve(H, psi_init, tlist, ops, measms, ntraj) + return exvals +end diff --git a/src/mcwf.jl b/src/mcwf.jl new file mode 100644 index 0000000..d2e1fbe --- /dev/null +++ b/src/mcwf.jl @@ -0,0 +1,155 @@ + +# quantum state as an ensemble of pure states +type QuStateEnsemble{S<:AbstractQuState} + state::S + decomp + QuStateEnsemble{S}(rho::S, d) = new(rho, d) +end + +function QuStateEnsemble{S<:AbstractQuState}(rho::S, d) + return QuStateEnsemble{S}(rho, d) +end + +function QuStateEnsemble(rho::QuState) + println("QuStateEnsemble{QuState}") + QuStateEnsemble(rho, eigfact(reshape(full(rho.elem), rho.nb, rho.nb))) +end + +function QuStateEnsemble(psi::QuStateVec) + println("QuStateEnsemble{QuStateVec}") + QuStateEnsemble(psi, nothing) +end + +# for a mixed state we use the spectral decomposition +function draw(e::QuStateEnsemble{QuState}) + println("draw from QuStateEnsemble{QuState}") + + r = rand() # draw random number from [0,1) + pacc = 0. + for i=1:length(e.decomp[:values]) + pacc = pacc + e.decomp[:values][i] + if pacc >= r + return QuStateVec(reshape(complex(e.decomp[:vectors][:,i]), e.state.nb,1), e.state.subnb) + end + end +end + +# for a pure state drawing is easy +function draw(e::QuStateEnsemble{QuStateVec}) + println("draw from QuStateEnsemble{QuStateVec}") + + return copy(e.state) +end + + +# MCWF method +# +# inputs: qme::AbstractLindbladQME +# state::AbstractQuState - initial state +# ntraj::Int - number of trajectories +# dt0::Float64 - time-step +# nsteps::Int - number of time-steps +# ops - Array of operators for which +# expectation values are calcluated at each time-step +# outputs: rho::QuState - final density matrix after nsteps +# exvals - nsteps*nops matrix with expectation values +# for each time-step and each op in ops +# +function mcwfpropagate(state::AbstractQuState, + qme::AbstractLindbladQME, + ntraj::Int, dt0::Float64, nsteps::Int, + ops=[]) + + jtol = 1.e-6 # jump tolerance + + nops = length(ops) + if nops > 0 + exvals = Array(Complex128, nsteps, nops) + else + exvals = [] + end + + # get information of QME + heff = eff_hamiltonian(qme) + c = linbladops(qme) + + # final state + rho = QuState( zeros(Complex128, state.nb, state.nb), state.subnb ) + + # initial state -> ensemble + qse = QuStateEnsemble(state) + + for traj=1:ntraj + println("starting trajectory $traj/$ntraj") + + # initial state vec for trajectory + psi = draw(qse) + + eps = rand() # draw random number from [0,1) + + for step=1:nsteps + + dt = dt0 + accdt = 0. + tj = 0. + + while accdt < dt0 + + # propagate one time-step + eff_prop = QuFixedStepPropagator( -im*heff, dt) + psi1 = propagate(eff_prop, psi) + + if abs(norm(psi1) - eps) < jtol + println("jump at t=$(tj + accdt + dt), norm=$(norm(psi1)), eps=$eps") + + PnS = 0. + for cind = 1:length(c) + PnS += real(expectationval(psi1, c[cind]'*c[cind])) + end + + Pn = 0. + cop = speye(size(heff,1)) + for cind = 1:length(c) + Pn += real(expectationval(psi1, c[cind]'*c[cind]))/PnS + println("Pn($cind)=$Pn") + if Pn >= eps + println("using operator $cind") + cop=c[cind] + break + end + end + + # now we have a jump + # psi = QuStateVec( cop*psi1.elem ) + psi = applyop(psi1, cop) + normalize!(psi) + + eps = rand() # draw new random number from [0,1) + accdt = accdt + dt + dt = dt0 - accdt + + elseif norm(psi1) < eps + dt = dt/2 # jump in the current interval + else + psi = copy(psi1) # no jump + accdt = accdt + dt + dt = dt0 - accdt + end + + println("dt = $dt, norm = $(norm(psi1)), $(norm(psi))") + end + + for nop = 1:nops + exvals[step, nop] = exvals[step, nop] + expectationval(normalize(psi), ops[nop])/ntraj + end + + tj = tj + dt0 + end # time-steps + + normalize!(psi) + rho = rho + QuState(psi)/ntraj + + end # trajectories + + return rho, exvals +end diff --git a/src/qme.jl b/src/qme.jl index e5c4189..080c1ba 100644 --- a/src/qme.jl +++ b/src/qme.jl @@ -1,7 +1,7 @@ ## QuDOS ## interface for quantum master equations -# +# include("nulls.jl") # abstract super-type for all QMEs @@ -50,6 +50,10 @@ function eff_hamiltonian(qme::AbstractLindbladQME) return heff end +# return Lindblad operators +linbladops(qme::AbstractLindbladQME) = qme.lbops + + # find stationary state, starting from rho # TODO: for some Lindblad ops this doesn't converge # and it seems not to be reliable