From 37125a69e2feae092a38a616aa0f674595960f5d Mon Sep 17 00:00:00 2001 From: acroy Date: Sun, 10 Aug 2014 22:27:57 +0200 Subject: [PATCH 1/6] Initial commit for MCWF method + some fixes for QuState. --- examples/mcwf_damped_ho.jl | 57 ++++++++++++++++++++ src/Qudos.jl | 7 ++- src/mcwf.jl | 105 +++++++++++++++++++++++++++++++++++++ src/qme.jl | 6 ++- 4 files changed, 172 insertions(+), 3 deletions(-) create mode 100644 examples/mcwf_damped_ho.jl create mode 100644 src/mcwf.jl 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 762efc2..a63e662 100644 --- a/src/Qudos.jl +++ b/src/Qudos.jl @@ -143,7 +143,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) @@ -159,7 +159,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 @@ -420,4 +420,7 @@ include("propagate.jl") # QME types and functionality include("qme.jl") +# MCWF method (aka quantum jump) +include("mcwf.jl") + end # module diff --git a/src/mcwf.jl b/src/mcwf.jl new file mode 100644 index 0000000..a4f3528 --- /dev/null +++ b/src/mcwf.jl @@ -0,0 +1,105 @@ +# MCWF method +# +# inputs: qme::AbstractLindbladQME +# psi0::QuStateVec +# 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(psi0::QuStateVec, + 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) + + rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) + + for traj=1:ntraj + println("starting trajectory $traj/$ntraj") + + eps = rand() # draw random number from [0,1) + psi = copy(psi0) + + 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 ) + 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 + psi = 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 From 934433474018503466646cc7c3c2b2fadcd28508 Mon Sep 17 00:00:00 2001 From: acroy Date: Mon, 11 Aug 2014 17:32:11 +0200 Subject: [PATCH 2/6] Support for QuState as initial state in mcwfpropagate. --- src/Qudos.jl | 2 +- src/mcwf.jl | 64 ++++++++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 58 insertions(+), 8 deletions(-) diff --git a/src/Qudos.jl b/src/Qudos.jl index a63e662..fde3a7c 100644 --- a/src/Qudos.jl +++ b/src/Qudos.jl @@ -6,7 +6,7 @@ module QuDOS # exports -export QuStateVec, QuState +export AbstractQuState, QuStateVec, QuState export QuFixedStepPropagator, QuKrylovPropagator diff --git a/src/mcwf.jl b/src/mcwf.jl index a4f3528..d2e1fbe 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -1,7 +1,51 @@ + +# 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 -# psi0::QuStateVec +# state::AbstractQuState - initial state # ntraj::Int - number of trajectories # dt0::Float64 - time-step # nsteps::Int - number of time-steps @@ -11,7 +55,7 @@ # exvals - nsteps*nops matrix with expectation values # for each time-step and each op in ops # -function mcwfpropagate(psi0::QuStateVec, +function mcwfpropagate(state::AbstractQuState, qme::AbstractLindbladQME, ntraj::Int, dt0::Float64, nsteps::Int, ops=[]) @@ -29,13 +73,19 @@ function mcwfpropagate(psi0::QuStateVec, heff = eff_hamiltonian(qme) c = linbladops(qme) - rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) + # 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") - eps = rand() # draw random number from [0,1) - psi = copy(psi0) + # initial state vec for trajectory + psi = draw(qse) + + eps = rand() # draw random number from [0,1) for step=1:nsteps @@ -70,7 +120,8 @@ function mcwfpropagate(psi0::QuStateVec, end # now we have a jump - psi = QuStateVec( cop*psi1.elem ) + # psi = QuStateVec( cop*psi1.elem ) + psi = applyop(psi1, cop) normalize!(psi) eps = rand() # draw new random number from [0,1) @@ -81,7 +132,6 @@ function mcwfpropagate(psi0::QuStateVec, dt = dt/2 # jump in the current interval else psi = copy(psi1) # no jump - psi = psi1 # no jump accdt = accdt + dt dt = dt0 - accdt end From 8a91c0197ed26fae392e59fecbfd81b4bd11535e Mon Sep 17 00:00:00 2001 From: Xiaodong Qi Date: Mon, 11 Aug 2014 12:06:52 -0600 Subject: [PATCH 3/6] Initialize mcsolve.jl and the corresponding example. --- examples/mcsolve_example.jl | 57 ++++++++++++++++++++ src/mcsolve.jl | 105 ++++++++++++++++++++++++++++++++++++ 2 files changed, 162 insertions(+) create mode 100644 examples/mcsolve_example.jl create mode 100644 src/mcsolve.jl diff --git a/examples/mcsolve_example.jl b/examples/mcsolve_example.jl new file mode 100644 index 0000000..d0fbddf --- /dev/null +++ b/examples/mcsolve_example.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 ) # This is an error here. +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/mcsolve.jl b/src/mcsolve.jl new file mode 100644 index 0000000..a4f3528 --- /dev/null +++ b/src/mcsolve.jl @@ -0,0 +1,105 @@ +# MCWF method +# +# inputs: qme::AbstractLindbladQME +# psi0::QuStateVec +# 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(psi0::QuStateVec, + 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) + + rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) + + for traj=1:ntraj + println("starting trajectory $traj/$ntraj") + + eps = rand() # draw random number from [0,1) + psi = copy(psi0) + + 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 ) + 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 + psi = 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 From 1559a9146d091ca4495e01ce3ac3917d428093d6 Mon Sep 17 00:00:00 2001 From: Xiaodong Qi Date: Tue, 12 Aug 2014 23:43:40 -0600 Subject: [PATCH 4/6] Basic version of the mcsolve function and example. Need to double check if it actually works, and use parallel. --- examples/mcsolve_example.jl | 19 +-- src/Qudos.jl | 5 +- src/mcsolve.jl | 237 +++++++++++++++++++++++------------- 3 files changed, 164 insertions(+), 97 deletions(-) diff --git a/examples/mcsolve_example.jl b/examples/mcsolve_example.jl index d0fbddf..3e86f73 100644 --- a/examples/mcsolve_example.jl +++ b/examples/mcsolve_example.jl @@ -1,4 +1,4 @@ -# quantum dynamics of a damped harmonic oscillator +# Quantum dynamics of a damped harmonic oscillator using the mcsolve.jl. using QuDOS import Winston: plot, oplot, xlabel, ylabel @@ -8,10 +8,10 @@ gamma = 0.1 nsteps = 50 dt = 0.25 -ntraj = 5 # number of trajectories +ntraj = 10 # number of trajectories # initialize operators etc -ad = QuDOS.creationop( n ) # This is an error here. +ad = QuDOS.creationop( n ) h = ad*ad' x = QuDOS.positionop( n ) @@ -42,14 +42,17 @@ pex[1] = real(QuDOS.expectationval(rho, p)) # 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) -rho, exvals = QuDOS.mcwfpropagate(cs, qme, ntraj, dt, nsteps, {x, p}) - -xex[2:end] = real(exvals[:,1]) -pex[2:end] = real(exvals[:,2]) +xex[2:end] = real(exvals[1,1,:]) +pex[2:end] = real(exvals[1,2,:]) # plots results using Winston -plot(dt*[0:nsteps], xex, "b-", dt*[0:nsteps], pex, "r-") +plot([0,tlist], xex, "b-", [0,tlist], 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--") diff --git a/src/Qudos.jl b/src/Qudos.jl index fde3a7c..2c4254f 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 @@ -423,4 +423,7 @@ 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/mcsolve.jl b/src/mcsolve.jl index a4f3528..f9e5c93 100644 --- a/src/mcsolve.jl +++ b/src/mcsolve.jl @@ -1,105 +1,166 @@ -# MCWF method -# -# inputs: qme::AbstractLindbladQME -# psi0::QuStateVec -# 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(psi0::QuStateVec, - 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 = [] +# mcsolve method for Monte Carlo Wavefunctions. It has a similar constructure as the QuTIP mcsolve function. + +function mcsolve(H::Union(SparseMatrixCSC{Complex{Float64},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::AbstractMatrix{Complex128}. - The Hamiltonian for the system. + # psi0::Vector{Complex128}. - The initial state of the system. + # tlist::Vector. - Time points as a vector. + # ops::Vector{AbstractMatrix{Complex128}}. - The Lindblad operators or jump operators as a vector with matrix elements. + # measms::Vector{AbstractMatrix{Complex128}}. - The measurement operators. + # ntraj::Vector{Int}=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. + # + + # BEtter to uniformality define new types: 1. QuDOSResult{K<:Float64,OP<:Complex128,N<:Int64}-->StateVec,Ops,Expect + # 2. QOperators + + #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 Hamiltonians and jump operators. + # c = Vector{AbstractMatrix{Complex128}} # Initialize the jump operators. + # c = ops #linbladops(qme) # Jump operators. + heff = H # Initialize the effective Hamiltonians. + for i=1:nops + heff = heff - 0.5im*hbar*ops[i]'*ops[i] #eff_hamiltonian(qme) end - # get information of QME - heff = eff_hamiltonian(qme) - c = linbladops(qme) - - rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) - - for traj=1:ntraj - println("starting trajectory $traj/$ntraj") - - eps = rand() # draw random number from [0,1) - psi = copy(psi0) - - 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])) + # rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) + + # Loop over all trajectories. + maxntraj=ntraj[end] + trajind=1 # Counting readout index in the ntraj vector. + for traj=1:maxntraj # Loop over all trajectories and including all given stopping points. + println("Starting trajectory $traj/$maxntraj:") + # The initial state and measurement records. + eps = rand() # Draw a random number from [0,1) to represent the probability that a quantum jump occurs. + # psi = Vector{AbstractVector{Complex128}}(nsteps) # Initialize the state for each trajectory. + psi = copy(psi0) # Initialize the quantum state vector after each evolution step. + psi1 = copy(psi0) # Initialize the quantum state vector before evolution. + psi_norm = 1 # Initialize the norm of the state vector. + t0=tlist[1] # Assme 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 # Step over all time steps. + dt = tlist[step]-t0 + if psi_norm > eps # No jump. + # Propagate one time-step until the jumping condition satisfies. + #eff_prop = QuFixedStepPropagator( -1im*heff/hbar, dt) + psi1 = expmv(psi,dt,-1im*heff/hbar) # propagate(eff_prop, psi) + println("t = $(tlist[step]), norm = $(norm(psi1)), evolved from $(norm(psi)).") + psi=copy(psi1) # State for the next step after the evolution time dt. + else # Otherwise, jump happens. + println("Jump at t=$(tlist[step]), norm=$(norm(psi1)), 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 probability of jumps. + # dp = 0. # Initialize the accumulated probability of jumps. + for cind = 1:nops + Pi[cind] = norm(ops[cind]*psi1).^2 # real(psi1'*ops[cind]'*ops[cind]*psi1) + PnS += Pi[cind] # Total probability of jumps without normalization. end - - Pn = 0. - cop = speye(size(heff,1)) - for cind = 1:length(c) - Pn += real(expectationval(psi1, c[cind]'*c[cind]))/PnS + # Judge which jump. + Pn = 0. # Initialize the normalized accumulated probability of jumps. + cop = eye(size(heff,1)) # Initialize the jump operator. + for cind = 1:nops + Pn += Pi[cind]/PnS println("Pn($cind)=$Pn") if Pn >= eps - println("using operator $cind") - cop=c[cind] + println("Using jump operator $cind") + cop=ops[cind] break end end - - # now we have a jump - psi = QuStateVec( cop*psi1.elem ) - normalize!(psi) - + # New state after the jump of cop. + psi1 = cop*psi #QuStateVec( cop*psi1.elem ) + # psi = QuStateVec(psi1) + psi = psi1./norm(psi1) # normalize!(psi1) + # Generate a new random number after the jump. 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 - psi = psi1 # no jump - accdt = accdt + dt - dt = dt0 - accdt end + psi_norm=norm(psi) # Norm after this step of evolution. + 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 - println("dt = $dt, norm = $(norm(psi1)), $(norm(psi))") - end + end # Time-steps - for nop = 1:nops - exvals[step, nop] = exvals[step, nop] + expectationval(normalize(psi), ops[nop])/ntraj - end + # Averaging expectation values for requested quantum trajectory numbers. + if traj==ntraj[trajind] + exvals[trajind,:,:] = exvals_temp[:,:]./traj + trajind += 1 + end + end # trajectories - tj = tj + dt0 - end # time-steps + return exvals +end - normalize!(psi) - rho = rho + QuState(psi)/ntraj +###################################################################################### +# 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 - end # trajectories +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. - return rho, exvals + psi_init = vec(full(psi0.elem)) + exvals = mcsolve(H, psi_init, tlist, ops, measms, ntraj) + return exvals end From 67720e766635f60ea4f84114291d2f4a8738c075 Mon Sep 17 00:00:00 2001 From: Xiaodong Qi Date: Wed, 13 Aug 2014 14:02:50 -0600 Subject: [PATCH 5/6] Make a better documentation and performance opt. --- examples/mcsolve_example.jl | 6 +-- src/mcsolve.jl | 77 ++++++++++++++++++------------------- 2 files changed, 40 insertions(+), 43 deletions(-) diff --git a/examples/mcsolve_example.jl b/examples/mcsolve_example.jl index 3e86f73..10e8184 100644 --- a/examples/mcsolve_example.jl +++ b/examples/mcsolve_example.jl @@ -8,7 +8,7 @@ gamma = 0.1 nsteps = 50 dt = 0.25 -ntraj = 10 # number of trajectories +ntraj = 100 # number of trajectories # initialize operators etc ad = QuDOS.creationop( n ) @@ -53,8 +53,8 @@ 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/2.)*dt*[0:0.5:nsteps]),"k--") -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)*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/src/mcsolve.jl b/src/mcsolve.jl index f9e5c93..cbe34ce 100644 --- a/src/mcsolve.jl +++ b/src/mcsolve.jl @@ -1,26 +1,28 @@ # mcsolve method for Monte Carlo Wavefunctions. It has a similar constructure as the QuTIP mcsolve function. -function mcsolve(H::Union(SparseMatrixCSC{Complex{Float64},Int64},AbstractMatrix{Complex128}), +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::AbstractMatrix{Complex128}. - The Hamiltonian for the system. - # psi0::Vector{Complex128}. - The initial state of the system. - # tlist::Vector. - Time points as a vector. + # 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::Vector{AbstractMatrix{Complex128}}. - The measurement operators. - # ntraj::Vector{Int}=500. - Number of trajectories for averaging. It can be a vector with increasing integers. + # 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. # - - # BEtter to uniformality define new types: 1. QuDOSResult{K<:Float64,OP<:Complex128,N<:Int64}-->StateVec,Ops,Expect - # 2. QOperators + # 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? @@ -33,53 +35,48 @@ function mcsolve(H::Union(SparseMatrixCSC{Complex{Float64},Int64},AbstractMatrix 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 Hamiltonians and jump operators. - # c = Vector{AbstractMatrix{Complex128}} # Initialize the jump operators. - # c = ops #linbladops(qme) # Jump operators. + # 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 - - # rho = QuState( zeros(Complex128, psi0.nb, psi0.nb), psi0.subnb ) + 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 and including all given stopping points. + for traj=1:maxntraj # Loop over all trajectories including all given recording points. println("Starting trajectory $traj/$maxntraj:") - # The initial state and measurement records. + # 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 = Vector{AbstractVector{Complex128}}(nsteps) # Initialize the state for each trajectory. - psi = copy(psi0) # Initialize the quantum state vector after each evolution step. - psi1 = copy(psi0) # Initialize the quantum state vector before evolution. - psi_norm = 1 # Initialize the norm of the state vector. - t0=tlist[1] # Assme tlist has a length larger than 1. + 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 # Step over all 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. - # Propagate one time-step until the jumping condition satisfies. - #eff_prop = QuFixedStepPropagator( -1im*heff/hbar, dt) - psi1 = expmv(psi,dt,-1im*heff/hbar) # propagate(eff_prop, psi) - println("t = $(tlist[step]), norm = $(norm(psi1)), evolved from $(norm(psi)).") - psi=copy(psi1) # State for the next step after the evolution time dt. + 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=$(norm(psi1)), eps=$eps .") + 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 probability of jumps. - # dp = 0. # Initialize the accumulated probability of jumps. + PnS = 0. # Initialize the unnormalized accumulated (sum) probability of jumps. for cind = 1:nops - Pi[cind] = norm(ops[cind]*psi1).^2 # real(psi1'*ops[cind]'*ops[cind]*psi1) + 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. - cop = eye(size(heff,1)) # Initialize the jump operator. for cind = 1:nops Pn += Pi[cind]/PnS println("Pn($cind)=$Pn") @@ -91,23 +88,23 @@ function mcsolve(H::Union(SparseMatrixCSC{Complex{Float64},Int64},AbstractMatrix end # New state after the jump of cop. psi1 = cop*psi #QuStateVec( cop*psi1.elem ) - # psi = QuStateVec(psi1) + # Normalize the state for the next time step. psi = psi1./norm(psi1) # normalize!(psi1) - # Generate a new random number after the jump. + 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 - psi_norm=norm(psi) # Norm after this step of evolution. + 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] - exvals[trajind,:,:] = exvals_temp[:,:]./traj + for j=1:nsteps, i=1:nmeasm + exvals[trajind,i,j] = copy(exvals_temp[i,j]./traj) + end trajind += 1 end end # trajectories From bb999e2523f301966a35bf4a157a7527faf16ea7 Mon Sep 17 00:00:00 2001 From: Xiaodong Qi Date: Mon, 25 Aug 2014 13:30:47 -0600 Subject: [PATCH 6/6] Make creationop work. --- src/operators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators.jl b/src/operators.jl index 407556e..9adfa2f 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -9,7 +9,7 @@ # function creationop( nb::Int ) - sparse( [2:nb], [1:nb], sqrt(linspace( 1, nb-1, nb-1)), nb, nb ) + sparse( [2:nb], [1:nb-1], sqrt(linspace( 1, nb-1, nb-1)), nb, nb ) end