diff --git a/REQUIRE b/REQUIRE index aea0520e..2c007aa4 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,6 +1,6 @@ julia 0.6 -Compat 0.52.0 -OrdinaryDiffEq 3.8.0 -DiffEqCallbacks 1.0 -StochasticDiffEq 4.2.4 +Compat 0.64.0 +OrdinaryDiffEq 3.19.1 +DiffEqCallbacks 1.1 +StochasticDiffEq 4.4.5 RecursiveArrayTools diff --git a/src/QuantumOptics.jl b/src/QuantumOptics.jl index cacd9c58..c1d594d8 100644 --- a/src/QuantumOptics.jl +++ b/src/QuantumOptics.jl @@ -12,11 +12,11 @@ export bases, Basis, GenericBasis, CompositeBasis, basis, operators_lazysum, LazySum, operators_lazyproduct, LazyProduct, operators_lazytensor, LazyTensor, - randstate, randoperator, superoperators, SuperOperator, DenseSuperOperator, SparseSuperOperator, spre, spost, liouvillian, fock, FockBasis, number, destroy, create, fockstate, coherentstate, displace, + randstate, randoperator, thermalstate, coherentthermalstate, phase_average, passive_state, spin, SpinBasis, sigmax, sigmay, sigmaz, sigmap, sigmam, spinup, spindown, subspace, SubspaceBasis, projector, particle, PositionBasis, MomentumBasis, samplepoints, gaussianstate, @@ -24,7 +24,7 @@ export bases, Basis, GenericBasis, CompositeBasis, basis, nlevel, NLevelBasis, transition, nlevelstate, manybody, ManyBodyBasis, fermionstates, bosonstates, manybodyoperator, onebodyexpect, occupation, - phasespace, qfunc, wigner, + phasespace, qfunc, wigner, coherentspinstate, metrics, tracenorm, tracenorm_h, tracenorm_nh, tracedistance, tracedistance_h, tracedistance_nh, entropy_vn, fidelity, ptranspose, PPT, @@ -48,10 +48,10 @@ include("operators_sparse.jl") include("operators_lazysum.jl") include("operators_lazyproduct.jl") include("operators_lazytensor.jl") -include("random.jl") include("superoperators.jl") include("spin.jl") include("fock.jl") +include("state_definitions.jl") include("subspace.jl") include("particle.jl") include("nlevel.jl") @@ -74,10 +74,12 @@ include("timecorrelations.jl") include("spectralanalysis.jl") include("semiclassical.jl") module stochastic + include("stochastic_definitions.jl") include("stochastic_schroedinger.jl") include("stochastic_master.jl") include("stochastic_semiclassical.jl") using .stochastic_schroedinger, .stochastic_master, .stochastic_semiclassical + using .stochastic_definitions end include("printing.jl") @@ -89,10 +91,10 @@ using .operators_sparse using .operators_lazysum using .operators_lazyproduct using .operators_lazytensor -using .random using .superoperators using .spin using .fock +using .state_definitions using .subspace using .particle using .nlevel diff --git a/src/master.jl b/src/master.jl index 950a82bc..be7bffe4 100644 --- a/src/master.jl +++ b/src/master.jl @@ -181,15 +181,15 @@ master_nh_dynamic(tspan, psi0::Ket, f::Function; kwargs...) = master_nh_dynamic( # Recasting needed for the ODE solver is just providing the underlying data -function recast!(x::Vector{Complex128}, rho::DenseOperator) - rho.data = reshape(x, size(rho.data)) +function recast!(x::Array{Complex128, 2}, rho::DenseOperator) + rho.data = x end -recast!(rho::DenseOperator, x::Vector{Complex128}) = nothing +recast!(rho::DenseOperator, x::Array{Complex128, 2}) = nothing function integrate_master(tspan, df::Function, rho0::DenseOperator, fout::Union{Void, Function}; kwargs...) tspan_ = convert(Vector{Float64}, tspan) - x0 = reshape(rho0.data, length(rho0)) + x0 = rho0.data state = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data) dstate = DenseOperator(rho0.basis_l, rho0.basis_r, rho0.data) integrate(tspan_, df, x0, state, dstate, fout; kwargs...) diff --git a/src/metrics.jl b/src/metrics.jl index 009671c5..a977bde2 100644 --- a/src/metrics.jl +++ b/src/metrics.jl @@ -5,7 +5,7 @@ export tracenorm, tracenorm_h, tracenorm_nh, entropy_vn, fidelity, ptranspose, PPT, negativity, logarithmic_negativity -using ..bases, ..operators, ..operators_dense +using ..bases, ..operators, ..operators_dense, ..states """ tracenorm(rho) @@ -147,8 +147,17 @@ S(ρ) = -Tr(ρ \\log(ρ)) = -\\sum_n λ_n\\log(λ_n) where ``λ_n`` are the eigenvalues of the density matrix ``ρ``, ``\\log`` is the natural logarithm and ``\\log(0) ≡ 0``. + +# Arguments +* `rho`: Density operator of which to calculate Von Neumann entropy. +* `tol=1e-15`: Tolerance for rounding errors in the computed eigenvalues. """ -entropy_vn(rho::DenseOperator) = sum([d == 0 ? 0 : -d*log(d) for d=eigvals(rho.data)]) +function entropy_vn(rho::DenseOperator; tol::Float64=1e-15) + evals = eigvals(rho.data) + evals[abs.(evals) .< tol] = 0.0 + sum([d == 0 ? 0 : -d*log(d) for d=evals]) +end +entropy_vn(psi::StateVector; kwargs...) = entropy_vn(dm(psi); kwargs...) """ fidelity(rho, sigma) diff --git a/src/operators.jl b/src/operators.jl index fb6b3ed0..4cdfadb3 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -4,7 +4,7 @@ export Operator, length, basis, dagger, ishermitian, tensor, embed, trace, ptrace, normalize, normalize!, expect, variance, expm, permutesystems, identityoperator -import Base: ==, +, -, *, /, length, trace, one, ishermitian, expm +import Base: ==, +, -, *, /, length, trace, one, ishermitian, expm, conj, conj! import ..bases: basis, tensor, ptrace, permutesystems, samebases, check_samebases, multiplicable import ..states: dagger, normalize, normalize! @@ -50,6 +50,9 @@ basis(a::Operator) = (check_samebases(a); a.basis_l) dagger(a::Operator) = arithmetic_unary_error("Hermitian conjugate", a) +conj(a::Operator) = arithmetic_unary_error("Complex conjugate", a) +conj!(a::Operator) = conj(a::Operator) + """ ishermitian(op::Operator) @@ -155,7 +158,8 @@ Expectation value of the given operator `op` for the specified `state`. `state` can either be a (density) operator or a ket. """ -expect(op::Operator, state::Ket) = dagger(state)*(op*state) + +expect(op::Operator, state::Ket) = dagger(state) * op * state expect(op::Operator, state::Operator) = trace(op*state) """ diff --git a/src/operators_dense.jl b/src/operators_dense.jl index f79d2909..b964bbbb 100644 --- a/src/operators_dense.jl +++ b/src/operators_dense.jl @@ -85,6 +85,10 @@ operators.dagger(x::DenseOperator) = DenseOperator(x.basis_r, x.basis_l, x.data' operators.ishermitian(A::DenseOperator) = ishermitian(A.data) operators.tensor(a::DenseOperator, b::DenseOperator) = DenseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(b.data, a.data)) + +operators.conj(a::DenseOperator) = DenseOperator(a.basis_l, a.basis_r, conj(a.data)) +operators.conj!(a::DenseOperator) = conj!(a.data) + """ tensor(x::Ket, y::Bra) @@ -121,6 +125,12 @@ end operators.normalize!(op::DenseOperator) = scale!(op.data, 1./trace(op)) +function operators.expect(op::DenseOperator, state::Ket)# where T <: Union{Ket, Bra} + check_samebases(op.basis_r, state.basis) + check_samebases(op.basis_l, state.basis) + state.data' * op.data * state.data +end + function operators.expect(op::DenseOperator, state::Operator) check_samebases(op.basis_r, state.basis_l) check_samebases(op.basis_l, state.basis_r) diff --git a/src/operators_sparse.jl b/src/operators_sparse.jl index 47fc0cd6..a165d03a 100644 --- a/src/operators_sparse.jl +++ b/src/operators_sparse.jl @@ -75,6 +75,9 @@ operators.tensor(a::SparseOperator, b::DenseOperator) = SparseOperator(tensor(a. operators.trace(op::SparseOperator) = (check_samebases(op); trace(op.data)) +operators.conj(op::SparseOperator) = SparseOperator(op.basis_l, op.basis_r, conj(op.data)) +operators.conj!(op::SparseOperator) = conj!(op.data) + function operators.ptrace(op::SparseOperator, indices::Vector{Int}) operators.check_ptrace_arguments(op, indices) shape = [op.basis_l.shape; op.basis_r.shape] @@ -84,6 +87,13 @@ function operators.ptrace(op::SparseOperator, indices::Vector{Int}) SparseOperator(b_l, b_r, data) end +function operators.expect(op::SparseOperator, state::Ket)# where T <: Union{Ket, Bra} + check_samebases(op.basis_r, state.basis) + check_samebases(op.basis_l, state.basis) + state.data' * op.data * state.data +end + + function operators.expect(op::SparseOperator, state::DenseOperator) check_samebases(op.basis_r, state.basis_l) check_samebases(op.basis_l, state.basis_r) diff --git a/src/phasespace.jl b/src/phasespace.jl index 4519aa9e..3f79910f 100644 --- a/src/phasespace.jl +++ b/src/phasespace.jl @@ -1,8 +1,8 @@ module phasespace -export qfunc, wigner +export qfunc, wigner, coherentspinstate -using ..bases, ..states, ..operators, ..operators_dense, ..fock +using ..bases, ..states, ..operators, ..operators_dense, ..fock, ..spin """ @@ -204,4 +204,21 @@ function _clenshaw(L::Int, abs2_2α::Float64, ρ::Matrix{Complex128}) end end +############################ coherent spin state ############################## +function coherentspinstate(b::SpinBasis, theta::Number, phi::Number, + result = Ket(b, Vector{Complex128}(length(b)))) + #theta = real(theta) + #phi = real(phi) + data = result.data + N = length(b)-1 + @inbounds for n=0:N + data[n+1] = + sqrt(factorial(BigInt(N))/(factorial(BigInt(n)) * + factorial(BigInt(N-n)))) * + (sin(theta/2.)*exp(1im*phi/2.))^n * + (cos(theta/2.)*exp(-1im*phi/2.))^(N-n) + end + return result +end + end #module diff --git a/src/random.jl b/src/random.jl deleted file mode 100644 index 1c4ea874..00000000 --- a/src/random.jl +++ /dev/null @@ -1,27 +0,0 @@ -module random - -export randstate, randoperator - -using ..bases, ..states, ..operators_dense - - -""" - randstate(basis) - -Calculate a random normalized ket state. -""" -function randstate(b::Basis) - psi = Ket(b, rand(Complex128, length(b))) - normalize!(psi) - psi -end - -""" - randoperator(b1[, b2]) - -Calculate a random unnormalized dense operator. -""" -randoperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, rand(Complex128, length(b1), length(b2))) -randoperator(b::Basis) = randoperator(b, b) - -end #module diff --git a/src/spectralanalysis.jl b/src/spectralanalysis.jl index bfeb8043..be90450c 100644 --- a/src/spectralanalysis.jl +++ b/src/spectralanalysis.jl @@ -18,6 +18,10 @@ about the way the calculation is done is needed, use the functions directly. More details can be found at [http://docs.julialang.org/en/stable/stdlib/linalg/]. +NOTE: Especially for small systems full diagonalization with Julia's `eig` +function is often more desirable. You can convert a sparse operator `A` to a +dense one using `full(A)`. + If the given operator is non-hermitian a warning is given. This behavior can be turned off using the keyword `warning=false`. """ @@ -41,15 +45,15 @@ end """ For sparse operators by default it only returns the 6 lowest eigenvalues. """ -function eigenstates(op::SparseOperator, n::Int=length(basis(op)); warning=true) +function eigenstates(op::SparseOperator, n::Int=6; warning::Bool=true, + info::Bool=true, kwargs...) b = basis(op) - if ishermitian(op) - data = Hermitian(op.data) - else - warning && warn(nonhermitian_warning) - data = op.data - end - D, V = eigs(data; nev=n, which=:SR) + # TODO: Change to sparese-Hermitian specific algorithm if more efficient + ishermitian(op) || (warning && warn(nonhermitian_warning)) + info && println("INFO: Defaulting to sparse diagonalization. + If storing the full operator is possible, it might be faster to do + eigenstates(full(op)). Set info=false to turn off this message.") + D, V = eigs(op.data; which=:SR, nev=n, kwargs...) states = [Ket(b, V[:, k]) for k=1:length(D)] D, states end @@ -84,7 +88,7 @@ end """ For sparse operators by default it only returns the 6 lowest eigenvalues. """ -eigenenergies(op::SparseOperator, n::Int=6; warning=true) = eigenstates(op, n; warning=warning)[1] +eigenenergies(op::SparseOperator, n::Int=6; kwargs...) = eigenstates(op, n; kwargs...)[1] arithmetic_unary_error = operators.arithmetic_unary_error diff --git a/src/state_definitions.jl b/src/state_definitions.jl new file mode 100644 index 00000000..b0c98dd7 --- /dev/null +++ b/src/state_definitions.jl @@ -0,0 +1,64 @@ +module state_definitions + +export randstate, randoperator, thermalstate, coherentthermalstate, phase_average, passive_state + +using ..bases, ..states, ..operators, ..operators_dense, ..fock + + +""" + randstate(basis) + +Calculate a random normalized ket state. +""" +function randstate(b::Basis) + psi = Ket(b, rand(Complex128, length(b))) + normalize!(psi) + psi +end + +""" + randoperator(b1[, b2]) + +Calculate a random unnormalized dense operator. +""" +randoperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, rand(Complex128, length(b1), length(b2))) +randoperator(b::Basis) = randoperator(b, b) + +""" + thermalstate(H,T) + +Thermal state ``exp(-H/T)/Tr[exp(-H/T)]``. +""" +function thermalstate(H::Operator,T::Real) + return normalize(expm(-full(H)/T)) +end + +""" + coherentthermalstate(basis::FockBasis,H,T,alpha) + +Coherent thermal state ``D(α)exp(-H/T)/Tr[exp(-H/T)]D^†(α)``. +""" +function coherentthermalstate(basis::FockBasis,H::Operator,T::Real,alpha::Number) + return displace(basis,alpha)*thermalstate(H,T)*dagger(displace(basis,alpha)) +end + +""" + phase_average(rho) + +Returns the phase-average of ``ρ`` containing only the diagonal elements. +""" +function phase_average(rho::DenseOperator) + return DenseOperator(basis(rho),diagm(diag(rho.data))) +end + +""" + passive_state(rho,IncreasingEigenenergies::Bool=true) + +Passive state ``π`` of ``ρ``. IncreasingEigenenergies=true means that higher indices correspond to higher energies. +""" +function passive_state(rho::DenseOperator,IncreasingEigenenergies::Bool=true) + return DenseOperator(basis(rho),diagm(sort(abs.(eigvals(rho.data)),rev=IncreasingEigenenergies))) +end + +end #module + diff --git a/src/stochastic_definitions.jl b/src/stochastic_definitions.jl new file mode 100644 index 00000000..ed4ad63d --- /dev/null +++ b/src/stochastic_definitions.jl @@ -0,0 +1,73 @@ +module stochastic_definitions + +export homodyne_carmichael + +using ...operators, ...states + +""" + stochastic.homodyne_carmichael(H0, C, theta) + +Helper function that defines the functions needed to compute homodyne detection +trajectories according to Carmichael with `stochastic.schroedinger_dynamic`. + +# Arguments +* `H0`: The deterministic, time-independent system Hamiltonian. +* `C`: Collapse operator (or vector of operators) of the detected output channel(s). +* `theta`: The phase difference between the local oscillator and the signal field. + Defines the operator of the measured quadrature as + ``X_\\theta = C e^{-i\\theta} + C^\\dagger e^{i\\theta}``. Needs to be a + vector of the same length as `C` if `C` is a vector. +* `normalize_expect=true`: Specifiy whether or not to normalize the state vector + when the expectation value in the nonlinear term is calculated. NOTE: + should only be set to `false` if the state is guaranteed to be normalized, + e.g. by setting `normalize_state=true` in `stochastic.schroedinger_dynamic`. + +Returns `(fdeterm, fstoch)`, where `fdeterm(t, psi) -> H` and +`fstoch(t, psi) -> Hs` are functions returning the deterministic and stochastic +part of the Hamiltonian required for calling `stochastic.schroedinger_dynamic`. + +The deterministic and stochastic parts of the Hamiltonian are constructed as + +```math +H_{det} = H_0 + H_{nl}, +``` + +where + +```math +H_{nl} = iCe^{-i\\theta} \\langle X_\\theta \\rangle - \\frac{i}{2} C^\\dagger C, +``` + +and + +```math +H_s = iCe^{-i\\theta}. +``` +""" +function homodyne_carmichael(H0::Operator, C::Vector{T}, theta::Vector{R}; + normalize_expect::Bool=true) where {T <: Operator, R <: Real} + n = length(C) + @assert n == length(theta) + Hs = 1.0im*C .* exp.(-1.0im .* theta) + X = C .* exp.(-1.0im .* theta) + dagger.(C) .* exp.(1.0im .* theta) + CdagC = -0.5im .* dagger.(C) .* C + + fstoch(t::Float64, psi::StateVector) = Hs + if normalize_expect + function H_nl_n(psi::StateVector) + psi_n = normalize(psi) + sum(expect(X[i], psi_n)*Hs[i] + CdagC[i] for i=1:n) + end + fdeterm_n(t::Float64, psi::StateVector) = H0 + H_nl_n(psi) + return fdeterm_n, fstoch + else + H_nl_un(psi::StateVector) = + sum(expect(X[i], psi)*Hs[i] + CdagC[i] for i=1:n) + fdeterm_un(t::Float64, psi::StateVector) = H0 + H_nl_un(psi) + return fdeterm_un, fstoch + end +end +homodyne_carmichael(H0::Operator, C::Operator, theta::Real; kwargs...) = + homodyne_carmichael(H0, [C], [theta]; kwargs...) + +end # module diff --git a/src/stochastic_master.jl b/src/stochastic_master.jl index 267b0ee4..76b7a646 100644 --- a/src/stochastic_master.jl +++ b/src/stochastic_master.jl @@ -12,7 +12,7 @@ const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void} const DiffArray = Union{Vector{Complex128}, Array{Complex128, 2}} """ - stochastic.master(tspan, rho0, H, J, Js; ) + stochastic.master(tspan, rho0, H, J, C; ) Time-evolution according to a stochastic master equation. @@ -27,65 +27,38 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster. * `H`: Deterministic part of the Hamiltonian. * `J`: Vector containing all deterministic jump operators which can be of any arbitrary operator type. -* `Js`: Vector containing the stochastic jump operators for a superoperator - describing a measurement which has the form of the standard linear - stochastic master equation, `Js[i]*rho + rho*Jsdagger[i]`. -* `Hs=nothing`: Vector containing additional stochastic terms of the Hamiltonian. +* `C`: Vector containing the stochastic operators for a superoperator + of the form `C[i]*rho + rho*Cdagger[i]`. * `rates=nothing`: Vector or matrix specifying the coefficients (decay rates) for the jump operators. If nothing is specified all rates are assumed to be 1. -* `rates_s=nothing`: Vector specifying the coefficients (decay rates) - for the stochastic jump operators. If nothing is specified all rates - are assumed to be 1. * `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump operators. If they are not given they are calculated automatically. -* `Jsdagger=dagger.(Js)`: Vector containing the hermitian conjugates of the - stochastic jump operators. * `fout=nothing`: If given, this function `fout(t, rho)` is called every time an output should be displayed. ATTENTION: The given state rho is not permanent! It is still in use by the ode solver and therefore must not be changed. -* `nonlinear=true`: Specify whether or not to include the nonlinear term - `expect(Js[i] + Jsdagger[i],rho)*rho` in the equation. This ensures - the trace of `rho` is conserved. * `kwargs...`: Further arguments are passed on to the ode solver. """ function master(tspan, rho0::DenseOperator, H::Operator, - J::Vector, Js::Vector; Hs::Union{Void, Vector}=nothing, - rates::DecayRates=nothing, rates_s::DecayRates=nothing, - Jdagger::Vector=dagger.(J), Jsdagger::Vector=dagger.(Js), + J::Vector, C::Vector; + rates::DecayRates=nothing, + Jdagger::Vector=dagger.(J), Cdagger::Vector=dagger.(C), fout::Union{Function,Void}=nothing, - nonlinear::Bool=true, kwargs...) tmp = copy(rho0) - if isa(rates_s, Matrix{Float64}) - throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates. - You may want to use diagonaljumps.")) - end - - n = length(Js) + (isa(Hs, Void) ? 0 : length(Hs)) + n = length(C) - if nonlinear - dmaster_stoch_nl(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stochastic_nl(dx, rho, Hs, rates_s, Js, Jsdagger, drho, n) - else - dmaster_stoch_lin(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stochastic(dx, rho, Hs, rates_s, Js, Jsdagger, drho, n) - end + dmaster_stoch(dx::DiffArray, t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = + dmaster_stochastic(dx, rho, C, Cdagger, drho, n) isreducible = check_master(rho0, H, J, Jdagger, rates) - check_master(rho0, H, Js, Jsdagger, rates_s) if !isreducible - dmaster_h_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) - if nonlinear - integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch_nl, rho0, fout, n; kwargs...) - else - integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch_lin, rho0, fout, n; kwargs...) - end + dmaster_h_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = + dmaster_h(rho, H, rates, J, Jdagger, drho, tmp) + integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...) else Hnh = copy(H) if typeof(rates) == Matrix{Float64} @@ -103,18 +76,15 @@ function master(tspan, rho0::DenseOperator, H::Operator, end Hnhdagger = dagger(Hnh) - dmaster_nh_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) - if nonlinear - integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch_nl, rho0, fout, n; kwargs...) - else - integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch_lin, rho0, fout, n; kwargs...) - end + dmaster_nh_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = + dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp) + integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...) end end master(tspan, psi0::Ket, args...; kwargs...) = master(tspan, dm(psi0), args...; kwargs...) """ - stochastic.master_dynamic(tspan, rho0, f, fs; ) + stochastic.master_dynamic(tspan, rho0, fdeterm, fstoch; ) Time-evolution according to a stochastic master equation with a dynamic Hamiltonian and J. @@ -126,20 +96,11 @@ dynamic Hamiltonian and J. * `fdeterm`: Function `f(t, rho) -> (H, J, Jdagger)` or `f(t, rho) -> (H, J, Jdagger, rates)` giving the deterministic part of the master equation. -* `fstoch`: Function `f(t, rho) -> (Js, Jsdagger)` or - `f(t, rho) -> (Js, Jsdagger, rates)` giving the stochastic superoperator - of the form `Js[i]*rho + rho*Jsdagger[i]`. +* `fstoch`: Function `f(t, rho) -> (C, Cdagger)` giving the stochastic + superoperator of the form `C[i]*rho + rho*Cdagger[i]`. * `rates=nothing`: Vector or matrix specifying the coefficients (decay rates) for the jump operators. If nothing is specified all rates are assumed to be 1. -* `rates_s=nothing`: Vector or matrix specifying the coefficients (decay rates) - for the stochastic jump operators. If nothing is specified all rates are assumed - to be 1. -* `fstoch_H=nothing`: Function `f(t, rho) -> Hs` providing a vector of operators - that correspond to stochastic terms of the Hamiltonian. -* `fstoch_J=nothing`: Function `f(t, rho) -> (J, Jdagger)` or - `f(t, rho) -> (J, Jdagger, rates)` giving a stochastic - Lindblad term. * `fout=nothing`: If given, this function `fout(t, rho)` is called every time an output should be displayed. ATTENTION: The given state rho is not permanent! It is still in use by the ode solver and therefore must not @@ -149,342 +110,56 @@ dynamic Hamiltonian and J. returned by `fstoch` and all optional functions. If unset, the number is calculated automatically from the function outputs. NOTE: Set this number if you want to avoid an initial calculation of function outputs! -* `nonlinear=true`: Specify whether or not to include the nonlinear term - `expect(Js[i] + Jsdagger[i],rho)*rho` in the equation. This ensures - the trace of `rho` is conserved. * `kwargs...`: Further arguments are passed on to the ode solver. """ function master_dynamic(tspan::Vector{Float64}, rho0::DenseOperator, fdeterm::Function, fstoch::Function; - fstoch_H::Union{Function, Void}=nothing, fstoch_J::Union{Function, Void}=nothing, - rates::DecayRates=nothing, rates_s::DecayRates=nothing, + rates::DecayRates=nothing, fout::Union{Function,Void}=nothing, - noise_processes::Int=0, nonlinear::Bool=true, + noise_processes::Int=0, kwargs...) tmp = copy(rho0) - if isa(rates_s, Matrix{Float64}) - throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates. - You may want to use diagonaljumps.")) - end - if noise_processes == 0 fs_out = fstoch(0, rho0) n = length(fs_out[1]) - if isa(fstoch_H, Function) - n += length(fstoch_H(0, rho0)) - end - if isa(fstoch_J, Function) - n += length(fstoch_J(0, rho0)[1]) - end else n = noise_processes end dmaster_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp) - if isa(fstoch_H, Void) && isa(fstoch_J, Void) - if nonlinear - dmaster_stoch_std_nl(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stoch_dynamic_nl(dx, t, rho, fstoch, rates_s, drho, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std_nl, rho0, fout, n; kwargs...) - else - dmaster_stoch_std(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stoch_dynamic(dx, t, rho, fstoch, rates_s, drho, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std, rho0, fout, n; kwargs...) - end - else - if nonlinear - dmaster_stoch_gen_nl(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stoch_dynamic_general_nl(dx, t, rho, fstoch, fstoch_H, fstoch_J, - rates, rates_s, drho, tmp, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen_nl, rho0, fout, n; kwargs...) - else - dmaster_stoch_gen(dx::DiffArray, - t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = - dmaster_stoch_dynamic_general(dx, t, rho, fstoch, fstoch_H, fstoch_J, - rates, rates_s, drho, tmp, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen, rho0, fout, n; kwargs...) - end - end + dmaster_stoch(dx::DiffArray, t::Float64, rho::DenseOperator, drho::DenseOperator, n::Int) = + dmaster_stoch_dynamic(dx, t, rho, fstoch, drho, n) + integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; kwargs...) end master_dynamic(tspan::Vector{Float64}, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...) -# Terms in SME -function dneumann(rho::DenseOperator, H::Operator, drho::DenseOperator) - operators.gemm!(-1.0im, H, rho, 0.0, drho) - operators.gemm!(1.0im, rho, H, 1.0, drho) -end - -function dlindblad(rho::DenseOperator, rates::Void, J::Vector, Jdagger::Vector, - drho::DenseOperator, tmp::DenseOperator, i::Int) - operators.gemm!(1, J[i], rho, 0, tmp) - operators.gemm!(1, tmp, Jdagger[i], 1, drho) -end -function dlindblad(rho::DenseOperator, rates::Vector{Float64}, J::Vector, - Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator, i::Int) - operators.gemm!(rates[i], J[i], rho, 0, tmp) - operators.gemm!(1, tmp, Jdagger[i], 1, drho) -end - -function dwiseman(rho::DenseOperator, rates::Void, J::Vector, - Jdagger::Vector, drho::DenseOperator, i::Int) - operators.gemm!(1, J[i], rho, 0, drho) - operators.gemm!(1, rho, Jdagger[i], 1, drho) -end -function dwiseman(rho::DenseOperator, rates::Vector{Float64}, J::Vector, - Jdagger::Vector, drho::DenseOperator, i::Int) - operators.gemm!(rates[i], J[i], rho, 0, drho) - operators.gemm!(rates[i], rho, Jdagger[i], 1, drho) -end - -function dwiseman_nl(rho::DenseOperator, rates::Void, J::Vector, - Jdagger::Vector, drho::DenseOperator, i::Int) - operators.gemm!(1, J[i], rho, 0, drho) - operators.gemm!(1, rho, Jdagger[i], 1, drho) - drho.data .-= (expect(J[i], rho) + expect(Jdagger[i], rho))*rho.data -end -function dwiseman_nl(rho::DenseOperator, rates::Vector{Float64}, J::Vector, - Jdagger::Vector, drho::DenseOperator, i::Int) - operators.gemm!(rates[i], J[i], rho, 0, drho) - operators.gemm!(rates[i], rho, Jdagger[i], 1, drho) - drho.data .-= rates[i]*(expect(J[i], rho) + expect(Jdagger[i], rho))*rho.data -end - # Derivative functions -function dmaster_stochastic(dx::Vector{Complex128}, rho::DenseOperator, H::Void, rates::DecayRates, - J::Vector, Jdagger::Vector, drho::DenseOperator, ::Int) +function dmaster_stochastic(dx::Vector{Complex128}, rho::DenseOperator, + C::Vector, Cdagger::Vector, drho::DenseOperator, ::Int) recast!(dx, drho) - dwiseman(rho, rates, J, Jdagger, drho, 1) - recast!(drho, dx) + operators.gemm!(1, C[1], rho, 0, drho) + operators.gemm!(1, rho, Cdagger[1], 1, drho) + drho.data .-= trace(drho)*rho.data end -function dmaster_stochastic(dx::Array{Complex128, 2}, rho::DenseOperator, H::Void, rates::DecayRates, - J::Vector, Jdagger::Vector, drho::DenseOperator, n::Int) +function dmaster_stochastic(dx::Array{Complex128, 2}, rho::DenseOperator, + C::Vector, Cdagger::Vector, drho::DenseOperator, n::Int) for i=1:n dx_i = @view dx[:, i] recast!(dx_i, drho) - dwiseman(rho, rates, J, Jdagger, drho, i) - recast!(drho, dx_i) - end -end -function dmaster_stochastic(dx::Array{Complex128, 2}, rho::DenseOperator, H::Vector, - rates::DecayRates, J::Vector, Jdagger::Vector, drho::DenseOperator, n::Int) - m = length(H) - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dneumann(rho, H[i-n+m], drho) + operators.gemm!(1, C[i], rho, 0, drho) + operators.gemm!(1, rho, Cdagger[i], 1, drho) + drho.data .-= trace(drho)*rho.data recast!(drho, dx_i) end - dmaster_stochastic(dx, rho, nothing, rates, J, Jdagger, drho, n-m) -end - -function dmaster_stochastic_nl(dx::Vector{Complex128}, rho::DenseOperator, H::Void, rates::DecayRates, - J::Vector, Jdagger::Vector, drho::DenseOperator, ::Int) - recast!(dx, drho) - dwiseman_nl(rho, rates, J, Jdagger, drho, 1) - recast!(drho, dx) -end -function dmaster_stochastic_nl(dx::Array{Complex128, 2}, rho::DenseOperator, H::Void, rates::DecayRates, - J::Vector, Jdagger::Vector, drho::DenseOperator, n::Int) - for i=1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dwiseman_nl(rho, rates, J, Jdagger, drho, i) - recast!(drho, dx_i) - end -end -function dmaster_stochastic_nl(dx::Array{Complex128, 2}, rho::DenseOperator, H::Vector, - rates::DecayRates, J::Vector, Jdagger::Vector, drho::DenseOperator, n::Int) - m = length(H) - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dneumann(rho, H[i-n+m], drho) - recast!(drho, dx_i) - end - dmaster_stochastic_nl(dx, rho, nothing, rates, J, Jdagger, drho, n-m) end function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, rho::DenseOperator, - f::Function, rates::DecayRates, drho::DenseOperator, n::Int) + f::Function, drho::DenseOperator, n::Int) result = f(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_ = rates - else - J, Jdagger, rates_ = result - end - dmaster_stochastic(dx, rho, nothing, rates_, J, Jdagger, drho, n) -end -function dmaster_stoch_dynamic_nl(dx::DiffArray, t::Float64, rho::DenseOperator, - f::Function, rates::DecayRates, drho::DenseOperator, n::Int) - result = f(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_ = rates - else - J, Jdagger, rates_ = result - end - dmaster_stochastic_nl(dx, rho, nothing, rates_, J, Jdagger, drho, n) -end - -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - H = fstoch_H(t, rho) - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_ = rates_s - else - J, Jdagger, rates_ = result - end - dmaster_stochastic(dx, rho, H, rates_, J, Jdagger, drho, n) -end -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - result_J = fstoch_J(t, rho) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J_stoch, J_stoch_dagger = result_J - rates_ = rates - else - J_stoch, J_stoch_dagger, rates_ = result_J - end - l = length(J_stoch) - - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_s_ = rates_s - else - J, Jdagger, rates_s_ = result - end - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, i-n+l) - recast!(drho, dx_i) - end - dmaster_stochastic(dx, rho, nothing, rates_s_, J, Jdagger, drho, n-l) -end -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - H = fstoch_H(t, rho) - - result_J = fstoch_J(t, rho) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J_stoch, J_stoch_dagger = result_J - rates_ = rates - else - J_stoch, J_stoch_dagger, rates_ = result_J - end - l = length(J_stoch) - - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_s_ = rates_s - else - J, Jdagger, rates_s_ = result - end - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, i-n+l) - recast!(drho, dx_i) - end - dmaster_stochastic(dx, rho, H, rates_s_, J, Jdagger, drho, n-l) -end - -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - H = fstoch_H(t, rho) - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_ = rates_s - else - J, Jdagger, rates_ = result - end - dmaster_stochastic_nl(dx, rho, H, rates_, J, Jdagger, drho, n) -end -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - result_J = fstoch_J(t, rho) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J_stoch, J_stoch_dagger = result_J - rates_ = rates - else - J_stoch, J_stoch_dagger, rates_ = result_J - end - l = length(J_stoch) - - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_s_ = rates_s - else - J, Jdagger, rates_s_ = result - end - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, i-n+l) - recast!(drho, dx_i) - end - dmaster_stochastic_nl(dx, rho, nothing, rates_s_, J, Jdagger, drho, n-l) -end -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, rho::DenseOperator, fstoch::Function, - fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - drho::DenseOperator, tmp::DenseOperator, n::Int) - H = fstoch_H(t, rho) - - result_J = fstoch_J(t, rho) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J_stoch, J_stoch_dagger = result_J - rates_ = rates - else - J_stoch, J_stoch_dagger, rates_ = result_J - end - l = length(J_stoch) - - result = fstoch(t, rho) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - J, Jdagger = result - rates_s_ = rates_s - else - J, Jdagger, rates_s_ = result - end - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, drho) - dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, i-n+l) - recast!(drho, dx_i) - end - dmaster_stochastic_nl(dx, rho, H, rates_s_, J, Jdagger, drho, n-l) + @assert 2 == length(result) + C, Cdagger = result + dmaster_stochastic(dx, rho, C, Cdagger, drho, n) end function integrate_master_stoch(tspan, df::Function, dg::Function, @@ -498,9 +173,11 @@ function integrate_master_stoch(tspan, df::Function, dg::Function, integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...) end -function recast!(x::SubArray{Complex128, 1}, rho::DenseOperator) +# TODO: Speed up by recasting to n-d arrays, remove vector methods +function recast!(x::Union{Vector{Complex128}, SubArray{Complex128, 1}}, rho::DenseOperator) rho.data = reshape(x, size(rho.data)) end recast!(state::DenseOperator, x::SubArray{Complex128, 1}) = (x[:] = state.data) +recast!(state::DenseOperator, x::Vector{Complex128}) = nothing end # module diff --git a/src/stochastic_schroedinger.jl b/src/stochastic_schroedinger.jl index c2bf265a..ca85dddb 100644 --- a/src/stochastic_schroedinger.jl +++ b/src/stochastic_schroedinger.jl @@ -8,6 +8,8 @@ using ...timeevolution import ...timeevolution: integrate_stoch, recast! import ...timeevolution.timeevolution_schroedinger: dschroedinger, dschroedinger_dynamic, check_schroedinger +import DiffEqCallbacks + const DiffArray = Union{Vector{Complex128}, Array{Complex128, 2}} """ @@ -25,10 +27,14 @@ Integrate stochastic Schrödinger equation. * `fout=nothing`: If given, this function `fout(t, state)` is called every time an output should be displayed. ATTENTION: The given state is neither normalized nor permanent! +* `normalize_state=false`: Specify whether or not to normalize the state after + each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector; fout::Union{Function,Void}=nothing, + normalize_state::Bool=false, + calback=nothing, kwargs...) tspan_ = convert(Vector{Float64}, tspan) @@ -42,7 +48,18 @@ function schroedinger(tspan, psi0::Ket, H::Operator, Hs::Vector; dschroedinger_stoch(dx::DiffArray, t::Float64, psi::Ket, dpsi::Ket, n::Int) = dschroedinger_stochastic(dx, psi, Hs, dpsi, n) - integrate_stoch(tspan_, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...) + if normalize_state + norm_func(u::Vector{Complex128}, t::Float64, integrator) = normalize!(u) + ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func; + func_everystep=true, + func_start=false) + else + ncb = nothing + end + + integrate_stoch(tspan_, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n; + ncb=ncb, + kwargs...) end schroedinger(tspan, psi0::Ket, H::Operator, Hs::Operator; kwargs...) = schroedinger(tspan, psi0, H, [Hs]; kwargs...) @@ -68,10 +85,13 @@ Integrate stochastic Schrödinger equation with dynamic Hamiltonian. from the function output. NOTE: Set this number if you want to avoid an initial calculation of the function output! +* `normalize_state=false`: Specify whether or not to normalize the state after + each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Function; fout::Union{Function,Void}=nothing, noise_processes::Int=0, + normalize_state::Bool=false, kwargs...) tspan_ = convert(Vector{Float64}, tspan) @@ -91,15 +111,27 @@ function schroedinger_dynamic(tspan, psi0::Ket, fdeterm::Function, fstoch::Funct t::Float64, psi::Ket, dpsi::Ket, n::Int) = dschroedinger_stochastic(dx, t, psi, fstoch, dpsi, n) - integrate_stoch(tspan, dschroedinger_determ, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...) + if normalize_state + norm_func(u::Vector{Complex128}, t::Float64, integrator) = normalize!(u) + ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func; + func_everystep=true, + func_start=false) + else + ncb = nothing + end + + integrate_stoch(tspan, dschroedinger_determ, dschroedinger_stoch, x0, state, + dstate, fout, n; + ncb=ncb, + kwargs...) end + function dschroedinger_stochastic(dx::Vector{Complex128}, psi::Ket, Hs::Vector{T}, dpsi::Ket, index::Int) where T <: Operator check_schroedinger(psi, Hs[index]) recast!(dx, dpsi) dschroedinger(psi, Hs[index], dpsi) - recast!(dpsi, dx) end function dschroedinger_stochastic(dx::Array{Complex128, 2}, psi::Ket, Hs::Vector{T}, dpsi::Ket, n::Int) where T <: Operator diff --git a/src/stochastic_semiclassical.jl b/src/stochastic_semiclassical.jl index 1bb4c7bd..cabf97de 100644 --- a/src/stochastic_semiclassical.jl +++ b/src/stochastic_semiclassical.jl @@ -10,7 +10,8 @@ using ...timeevolution import ...timeevolution: integrate_stoch import ...timeevolution.timeevolution_schroedinger: dschroedinger, dschroedinger_dynamic using ...stochastic -import ...stochastic.stochastic_master: dneumann, dwiseman, dwiseman_nl, dlindblad + +import DiffEqCallbacks const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void} const DiffArray = Union{Vector{Complex128}, Array{Complex128, 2}} @@ -38,12 +39,19 @@ Integrate time-dependent Schrödinger equation coupled to a classical system. * `fout=nothing`: If given, this function `fout(t, state)` is called every time an output should be displayed. ATTENTION: The given state is neither normalized nor permanent! -* `noise_processes=0`: Number of distinct white-noise processes in the equation. +* `noise_processes=0`: Number of distinct quantum noise processes in the equation. This number has to be equal to the total number of noise operators - returned by `fstoch_quantum`. Add 1 if `fstoch_classical` is specificed. - If unset, the number is automatically calculated from function outputs. + returned by `fstoch`. If unset, the number is calculated automatically + from the function output. NOTE: Set this number if you want to avoid an initial calculation of - function outputs! + the function output! +* `noise_prototype_classical=nothing`: The equivalent of the optional argument + `noise_rate_prototype` in `StochasticDiffEq` for the classical + stochastic function `fstoch_classical` only. Must be set for + non-diagonal classical noise or combinations of quantum and classical + noise. See the documentation for details. +* `normalize_state=false`: Specify whether or not to normalize the state after + each time step taken by the solver. * `kwargs...`: Further arguments are passed on to the ode solver. """ function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Function, @@ -51,6 +59,8 @@ function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Functio fstoch_classical::Union{Void, Function}=nothing, fout::Union{Function,Void}=nothing, noise_processes::Int=0, + noise_prototype_classical=nothing, + normalize_state::Bool=false, kwargs...) tspan_ = convert(Vector{Float64}, tspan) dschroedinger_det(t::Float64, state::State{Ket}, dstate::State{Ket}) = semiclassical.dschroedinger_dynamic(t, state, fquantum, fclassical, dstate) @@ -70,17 +80,35 @@ function schroedinger_semiclassical(tspan, state0::State{Ket}, fquantum::Functio fs_out = fstoch_quantum(0.0, state0.quantum, state0.classical) n += length(fs_out) end - if isa(fstoch_classical, Function) - n += 1 - end else n = noise_processes end + if n > 0 && isa(fstoch_classical, Function) + if isa(noise_prototype_classical, Void) + throw(ArgumentError("noise_prototype_classical must be set for combinations of quantum and classical noise!")) + end + end + + if normalize_state + len_q = length(state0.quantum) + function norm_func(u::Vector{Complex128}, t::Float64, integrator) + u .= [normalize!(u[1:len_q]), u[len_q+1:end];] + end + ncb = DiffEqCallbacks.FunctionCallingCallback(norm_func; + func_everystep=true, + func_start=false) + else + ncb = nothing + end + dschroedinger_stoch(dx::DiffArray, - t::Float64, state::State{Ket}, dstate::State{Ket}, index::Int) = - dschroedinger_stochastic(dx, t, state, fstoch_quantum, fstoch_classical, dstate, index) - integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n; kwargs...) + t::Float64, state::State{Ket}, dstate::State{Ket}, n::Int) = + dschroedinger_stochastic(dx, t, state, fstoch_quantum, fstoch_classical, dstate, n) + integrate_stoch(tspan_, dschroedinger_det, dschroedinger_stoch, x0, state, dstate, fout, n; + noise_prototype_classical = noise_prototype_classical, + ncb=ncb, + kwargs...) end """ @@ -101,54 +129,44 @@ non-hermitian Hamiltonian and then calls master_nh which is slightly faster. part of the master equation. * `fclassical`: Function `f(t, rho, u, du)` that calculates the classical derivatives `du`. -* `fstoch_quantum=nothing`: Function `f(t, rho, u) -> Js, Jsdagger` or - `f(t, psi, u) -> Js, Jsdagger, rates_s` that returns the stochastic - operator for the superoperator of the form `Js[i]*rho + rho*Jsdagger[i]`. +* `fstoch_quantum=nothing`: Function `f(t, rho, u) -> C, Cdagger` + that returns the stochastic operator for the superoperator of the form + `C[i]*rho + rho*Cdagger[i]`. * `fstoch_classical=nothing`: Function `f(t, rho, u, du)` that calculates the stochastic terms of the derivative `du`. -* `fstoch_H=nothing`: Function `f(t, rho, u) -> Hs` providing a vector of operators - that correspond to stochastic terms of the Hamiltonian. -* `fstoch_J=nothing`: Function `f(t, rho, u) -> (J, Jdagger)` or - `f(t, rho, u) -> (J, Jdagger, rates)` giving a stochastic - Lindblad term. * `rates=nothing`: Vector or matrix specifying the coefficients (decay rates) for the jump operators. If nothing is specified all rates are assumed to be 1. -* `rates_s=nothing`: Vector or matrix specifying the coefficients (decay rates) - for the stochastic jump operators. If nothing is specified all rates are assumed - to be 1. * `fout=nothing`: If given, this function `fout(t, rho)` is called every time an output should be displayed. ATTENTION: The given state rho is not permanent! It is still in use by the ode solver and therefore must not be changed. -* `noise_processes=0`: Number of distinct white-noise processes in the equation. +* `noise_processes=0`: Number of distinct quantum noise processes in the equation. This number has to be equal to the total number of noise operators - returned by all stochastic functions. Add 1 for classical noise. - If unset, the number is calculated automatically from the function outputs. + returned by `fstoch`. If unset, the number is calculated automatically + from the function output. NOTE: Set this number if you want to avoid an initial calculation of - function outputs! -* `nonlinear=true`: Specify whether or not to include the nonlinear term - `expect(Js[i] + Jsdagger[i],rho)*rho` in the equation. This ensures - the trace of `rho` is conserved. + the function output! +* `noise_prototype_classical=nothing`: The equivalent of the optional argument + `noise_rate_prototype` in `StochasticDiffEq` for the classical + stochastic function `fstoch_classical` only. Must be set for + non-diagonal classical noise or combinations of quantum and classical + noise. See the documentation for details. * `kwargs...`: Further arguments are passed on to the ode solver. """ function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator}, fquantum::Function, fclassical::Function; fstoch_quantum::Union{Function, Void}=nothing, fstoch_classical::Union{Function, Void}=nothing, - fstoch_H::Union{Function, Void}=nothing, fstoch_J::Union{Function, Void}=nothing, - rates::DecayRates=nothing, rates_s::DecayRates=nothing, + rates::DecayRates=nothing, fout::Union{Function,Void}=nothing, - noise_processes::Int=0, nonlinear::Bool=true, + noise_processes::Int=0, + noise_prototype_classical=nothing, + nonlinear::Bool=true, kwargs...) tmp = copy(rho0.quantum) - - if isa(rates_s, Matrix{Float64}) - throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates. - You may want to set them as ones or use diagonaljumps.")) - end - if isa(fstoch_quantum, Void) && isa(fstoch_classical, Void) && isa(fstoch_H, Void) && isa(fstoch_J, Void) + if isa(fstoch_quantum, Void) && isa(fstoch_classical, Void) throw(ArgumentError("No stochastic functions provided!")) end @@ -158,56 +176,31 @@ function master_semiclassical(tspan::Vector{Float64}, rho0::State{DenseOperator} fq_out = fstoch_quantum(0, rho0.quantum, rho0.classical) n += length(fq_out[1]) end - if isa(fstoch_classical, Function) - n += 1 - end - if isa(fstoch_H, Function) - n += length(fstoch_H(0, rho0.quantum, rho0.classical)) - end - if isa(fstoch_J, Function) - n += length(fstoch_J(0, rho0.quantum, rho0.classical)[1]) - end else n = noise_processes end - dmaster_determ(t::Float64, rho::State{DenseOperator}, drho::State{DenseOperator}) = - dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp) - if isa(fstoch_H, Void) && isa(fstoch_J, Void) - if nonlinear - dmaster_stoch_std_nl(dx::DiffArray, t::Float64, rho::State{DenseOperator}, - drho::State{DenseOperator}, n::Int) = - dmaster_stoch_dynamic_nl(dx, t, rho, fstoch_quantum, fstoch_classical, - rates_s, drho, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std_nl, rho0, fout, n; kwargs...) - else - dmaster_stoch_std(dx::DiffArray, t::Float64, rho::State{DenseOperator}, - drho::State{DenseOperator}, n::Int) = - dmaster_stoch_dynamic(dx, t, rho, fstoch_quantum, fstoch_classical, - rates_s, drho, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std, rho0, fout, n; kwargs...) - end - else - if nonlinear - dmaster_stoch_gen_nl(dx::DiffArray, t::Float64, rho::State{DenseOperator}, - drho::State{DenseOperator}, n::Int) = - dmaster_stoch_dynamic_general_nl(dx, t, rho, fstoch_quantum, - fstoch_classical, fstoch_H, fstoch_J, rates, rates_s, - drho, tmp, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen_nl, rho0, fout, n; kwargs...) - else - dmaster_stoch_gen(dx::DiffArray, t::Float64, rho::State{DenseOperator}, - drho::State{DenseOperator}, n::Int) = - dmaster_stoch_dynamic_general(dx, t, rho, fstoch_quantum, - fstoch_classical, fstoch_H, fstoch_J, rates, rates_s, - drho, tmp, n) - integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen, rho0, fout, n; kwargs...) + if n > 0 && isa(fstoch_classical, Function) + if isa(noise_prototype_classical, Void) + throw(ArgumentError("noise_prototype_classical must be set for combinations of quantum and classical noise!")) end end + + dmaster_determ(t::Float64, rho::State{DenseOperator}, drho::State{DenseOperator}) = + dmaster_h_dynamic(t, rho, fquantum, fclassical, rates, drho, tmp) + + dmaster_stoch(dx::DiffArray, t::Float64, rho::State{DenseOperator}, + drho::State{DenseOperator}, n::Int) = + dmaster_stoch_dynamic(dx, t, rho, fstoch_quantum, fstoch_classical, drho, n) + + integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch, rho0, fout, n; + noise_prototype_classical=noise_prototype_classical, + kwargs...) end master_semiclassical(tspan::Vector{Float64}, psi0::State{Ket}, args...; kwargs...) = master_semiclassical(tspan, dm(psi0), args...; kwargs...) +# Derivative functions function dschroedinger_stochastic(dx::Vector{Complex128}, t::Float64, state::State{Ket}, fstoch_quantum::Function, fstoch_classical::Void, dstate::State{Ket}, ::Int) @@ -227,270 +220,60 @@ function dschroedinger_stochastic(dx::Array{Complex128, 2}, recast!(dstate, dx_i) end end -function dschroedinger_stochastic(dx::Vector{Complex128}, t::Float64, - state::State{Ket}, fstoch_quantum::Void, fstoch_classical::Function, - dstate::State{Ket}, ::Int) - recast!(dx, dstate) - fstoch_classical(t, state.quantum, state.classical, dstate.classical) - recast!(dstate, dx) +function dschroedinger_stochastic(dx::DiffArray, t::Float64, + state::State{Ket}, fstoch_quantum::Void, fstoch_classical::Function, + dstate::State{Ket}, ::Int) + dclassical = @view dx[length(state.quantum)+1:end, :] + fstoch_classical(t, state.quantum, state.classical, dclassical) end function dschroedinger_stochastic(dx::Array{Complex128, 2}, t::Float64, state::State{Ket}, fstoch_quantum::Function, fstoch_classical::Function, dstate::State{Ket}, n::Int) - dschroedinger_stochastic(dx, t, state, fstoch_quantum, nothing, dstate, n-1) + dschroedinger_stochastic(dx, t, state, fstoch_quantum, nothing, dstate, n) - dx_i = @view dx[:, end] - recast!(dx_i, dstate) - fstoch_classical(t, state.quantum, state.classical, dstate.classical) - recast!(dstate, dx_i) + dx_i = @view dx[length(state.quantum)+1:end, n+1:end] + fstoch_classical(t, state.quantum, state.classical, dx_i) end function dmaster_stoch_dynamic(dx::Vector{Complex128}, t::Float64, state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Void, - rates_s::DecayRates, dstate::State{DenseOperator}, ::Int) + fstoch_classical::Void, dstate::State{DenseOperator}, ::Int) result = fstoch_quantum(t, state.quantum, state.classical) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - Js, Jsdagger = result - rates_s_ = rates_s - else - Js, Jsdagger, rates_s_ = result - end + @assert length(result) == 2 + C, Cdagger = result recast!(dx, dstate) - dwiseman(state.quantum, rates_s_, Js, Jsdagger, dstate.quantum, 1) + operators.gemm!(1, C[1], state.quantum, 0, dstate.quantum) + operators.gemm!(1, state.quantum, Cdagger[1], 1, dstate.quantum) + dstate.quantum.data .-= trace(dstate.quantum)*state.quantum.data recast!(dstate, dx) end function dmaster_stoch_dynamic(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Void, - rates_s::DecayRates, dstate::State{DenseOperator}, n::Int) + fstoch_classical::Void, dstate::State{DenseOperator}, n::Int) result = fstoch_quantum(t, state.quantum, state.classical) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - Js, Jsdagger = result - rates_s_ = rates_s - else - Js, Jsdagger, rates_s_ = result - end + @assert length(result) == 2 + C, Cdagger = result for i=1:n dx_i = @view dx[:, i] recast!(dx_i, dstate) - dwiseman(state.quantum, rates_s_, Js, Jsdagger, dstate.quantum, i) + operators.gemm!(1, C[i], state.quantum, 0, dstate.quantum) + operators.gemm!(1, state.quantum, Cdagger[i], 1, dstate.quantum) + dstate.quantum.data .-= trace(dstate.quantum)*state.quantum.data recast!(dstate, dx_i) end end -function dmaster_stoch_dynamic(dx::Vector{Complex128}, t::Float64, +function dmaster_stoch_dynamic(dx::DiffArray, t::Float64, state::State{DenseOperator}, fstoch_quantum::Void, - fstoch_classical::Function, - rates_s::DecayRates, dstate::State{DenseOperator}, ::Int) - recast!(dx, dstate) - fstoch_classical(t, state.quantum, state.classical, dstate.classical) - recast!(dstate, dx) + fstoch_classical::Function, dstate::State{DenseOperator}, ::Int) + dclassical = @view dx[length(state.quantum)+1:end, :] + fstoch_classical(t, state.quantum, state.classical, dclassical) end function dmaster_stoch_dynamic(dx::Array{Complex128, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Void, - fstoch_classical::Function, - rates_s::DecayRates, dstate::State{DenseOperator}, i::Int) - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - fstoch_classical(t, state.quantum, state.classical, dstate.classical) - recast!(dstate, dx_i) -end -function dmaster_stoch_dynamic(dx::Array{Complex128, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Function, - rates_s::DecayRates, dstate::State{DenseOperator}, n::Int) - dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, nothing, rates_s, dstate, n-1) - dmaster_stoch_dynamic(dx, t, state, nothing, fstoch_classical, rates_s, dstate, n) -end -dmaster_stoch_dynamic(dx::DiffArray, t::Float64, state::State{DenseOperator}, ::Void, ::Void, args...) = nothing - -function dmaster_stoch_dynamic_nl(dx::Vector{Complex128}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Void, - rates_s::DecayRates, dstate::State{DenseOperator}, ::Int) - result = fstoch_quantum(t, state.quantum, state.classical) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - Js, Jsdagger = result - rates_s_ = rates_s - else - Js, Jsdagger, rates_s_ = result - end - recast!(dx, dstate) - dwiseman_nl(state.quantum, rates_s_, Js, Jsdagger, dstate.quantum, 1) - recast!(dstate, dx) -end -function dmaster_stoch_dynamic_nl(dx::Array{Complex128, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Void, - rates_s::DecayRates, dstate::State{DenseOperator}, n::Int) - result = fstoch_quantum(t, state.quantum, state.classical) - @assert 2 <= length(result) <= 3 - if length(result) == 2 - Js, Jsdagger = result - rates_s_ = rates_s - else - Js, Jsdagger, rates_s_ = result - end - for i=1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dwiseman_nl(state.quantum, rates_s_, Js, Jsdagger, dstate.quantum, i) - recast!(dstate, dx_i) - end -end -function dmaster_stoch_dynamic_nl(dx::DiffArray, - t::Float64, state::State{DenseOperator}, fstoch_quantum::Void, - fstoch_classical::Function, rates_s::DecayRates, dstate::State{DenseOperator}, n::Int) - dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, fstoch_classical, rates_s, dstate, n) -end -function dmaster_stoch_dynamic_nl(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, fstoch_quantum::Function, - fstoch_classical::Function, - rates_s::DecayRates, dstate::State{DenseOperator}, n::Int) - dmaster_stoch_dynamic_nl(dx, t, state, fstoch_quantum, nothing, rates_s, dstate, n-1) - dmaster_stoch_dynamic_nl(dx, t, state, nothing, fstoch_classical, rates_s, dstate, 1) -end -dmaster_stoch_dynamic_nl(dx::DiffArray, t::Float64, state::State{DenseOperator}, ::Void, ::Void, args...) = nothing - - -function dmaster_stoch_dynamic_general(dx::Vector{Complex128}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Void, fstoch_classical::Void, - fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, ::Int) - H = fstoch_H(t, state.quantum, state.classical) - recast!(dx, dstate) - dneumann(state.quantum, H[1], dstate.quantum) - recast!(dstate, dx) -end -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void}, - fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - H = fstoch_H(t, state.quantum, state.classical) - m = length(H) - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dneumann(state.quantum, H[i-n+m], dstate.quantum) - end - dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, fstoch_classical, - rates_s, dstate, n-m) -end -function dmaster_stoch_dynamic_general(dx::Vector{Complex128}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Void, fstoch_classical::Void, - fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, ::Int) - result_J = fstoch_J(t, state.quantum, state.classical) - - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J, Jdagger = result_J - rates_ = rates - else - J, Jdagger, rates_ = result_J - end - recast!(dx, dstate) - dlindblad(state.quantum, rates_, J, Jdagger, dstate.quantum, tmp, 1) - recast!(dstate, dx) -end -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void}, - fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - result_J = fstoch_J(t, state.quantum, state.classical) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J, Jdagger = result_J - rates_ = rates - else - J, Jdagger, rates_ = result_J - end - l = length(J) - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dlindblad(state.quantum, rates_, J, Jdagger, dstate.quantum, tmp, i-n+l) - recast!(dstate, dx_i) - end - dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, fstoch_classical, - rates_s, dstate, n-l) -end -function dmaster_stoch_dynamic_general(dx::Array{Complex128, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Union{Function, Void}, - fstoch_classical::Union{Function, Void}, - fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - H = fstoch_H(t, state.quantum, state.classical) - m = length(H) - - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dneumann(state.quantum, H[i-n+m], dstate.quantum) - recast!(dstate, dx_i) - end - dmaster_stoch_dynamic_general(dx, t, state, fstoch_quantum, fstoch_classical, - nothing, fstoch_J, rates, rates_s, dstate, tmp, n-m) -end + fstoch_classical::Function, dstate::State{DenseOperator}, n::Int) + dmaster_stoch_dynamic(dx, t, state, fstoch_quantum, nothing, dstate, n) -dmaster_stoch_dynamic_general_nl(dx::Vector{Complex128}, args...; kwargs...) = - dmaster_stoch_dynamic_general(dx, args...; kwargs...) -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void}, - fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - H = fstoch_H(t, state.quantum, state.classical) - m = length(H) - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dneumann(state.quantum, H[i-n+m], dstate.quantum) - end - dmaster_stoch_dynamic_nl(dx, t, state, fstoch_quantum, fstoch_classical, - rates_s, dstate, n-m) -end -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, state::State{DenseOperator}, - fstoch_quantum::Union{Function, Void}, fstoch_classical::Union{Function, Void}, - fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - result_J = fstoch_J(t, state.quantum, state.classical) - @assert 2 <= length(result_J) <= 3 - if length(result_J) == 2 - J, Jdagger = result_J - rates_ = rates - else - J, Jdagger, rates_ = result_J - end - l = length(J) - - for i=n-l+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dlindblad(state.quantum, rates_, J, Jdagger, dstate.quantum, tmp, i-n+l) - recast!(dstate, dx_i) - end - dmaster_stoch_dynamic_nl(dx, t, state, fstoch_quantum, fstoch_classical, - rates_s, dstate, n-l) -end -function dmaster_stoch_dynamic_general_nl(dx::Array{Complex128, 2}, t::Float64, - state::State{DenseOperator}, fstoch_quantum::Union{Function, Void}, - fstoch_classical::Union{Function, Void}, - fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates, - dstate::State{DenseOperator}, tmp::DenseOperator, n::Int) - H = fstoch_H(t, state.quantum, state.classical) - m = length(H) - - for i=n-m+1:n - dx_i = @view dx[:, i] - recast!(dx_i, dstate) - dneumann(state.quantum, H[i-n+m], dstate.quantum) - recast!(dstate, dx_i) - end - dmaster_stoch_dynamic_general_nl(dx, t, state, fstoch_quantum, fstoch_classical, - nothing, fstoch_J, rates, rates_s, dstate, tmp, n-m) + dx_i = @view dx[length(state.quantum)+1:end, n+1:end] + fstoch_classical(t, state.quantum, state.classical, dx_i) end function integrate_master_stoch(tspan, df::Function, dg::Function, diff --git a/src/superoperators.jl b/src/superoperators.jl index f46c2862..1f386c18 100644 --- a/src/superoperators.jl +++ b/src/superoperators.jl @@ -104,6 +104,11 @@ function *(a::SuperOperator, b::DenseOperator) data = a.data*reshape(b.data, length(b.data)) return DenseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2]))) end +function *(a::SuperOperator, b::SparseOperator) + check_multiplicable(a, b) + data = a.data*reshape(b.data, length(b.data)) + return SparseOperator(a.basis_l[1], a.basis_l[2], reshape(data, length(a.basis_l[1]), length(a.basis_l[2]))) +end function *(a::SuperOperator, b::SuperOperator) check_multiplicable(a, b) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index d44b43ae..ebb1b981 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -2,6 +2,8 @@ using ..metrics import OrdinaryDiffEq, DiffEqCallbacks, StochasticDiffEq +const DiffArray = Union{Vector{Complex128}, Array{Complex128, 2}} + function recast! end """ @@ -10,19 +12,19 @@ function recast! end Integrate using OrdinaryDiffEq """ -function integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, +function integrate(tspan::Vector{Float64}, df::Function, x0::DiffArray, state::T, dstate::T, fout::Function; alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(), steady_state = false, tol = 1e-3, save_everystep = false, callback = nothing, kwargs...) where T - function df_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t) + function df_(dx::DiffArray, x::DiffArray, p, t) recast!(x, state) recast!(dx, dstate) df(t, state, dstate) recast!(dstate, dx) end - function fout_(x::Vector{Complex128}, t::Float64, integrator) + function fout_(x::DiffArray, t::Float64, integrator) recast!(x, state) fout(t, state) end @@ -64,7 +66,7 @@ function integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, out.t,out.saveval end -function integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128}, +function integrate(tspan::Vector{Float64}, df::Function, x0::DiffArray, state::T, dstate::T, ::Void; kwargs...) where T function fout(t::Float64, state::T) copy(state) @@ -96,8 +98,11 @@ Integrate using StochasticDiffEq function integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Function, x0::Vector{Complex128}, state::T, dstate::T, fout::Function, n::Int; save_everystep = false, callback=nothing, - alg = nothing, noise_rate_prototype = Array{Complex128}(length(state), n), - noise=StochasticDiffEq.RealWienerProcess(0.0, randn(n)), + alg::StochasticDiffEq.StochasticDiffEqAlgorithm=StochasticDiffEq.EM(), + noise_rate_prototype = nothing, + noise_prototype_classical = nothing, + noise=nothing, + ncb=nothing, kwargs...) where T function df_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t) @@ -118,6 +123,18 @@ function integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Function, x0: fout(t, state) end + nc = isa(noise_prototype_classical, Void) ? 0 : size(noise_prototype_classical)[2] + if isa(noise, Void) && n > 0 + noise_ = StochasticDiffEq.RealWienerProcess!(0.0, randn(n + nc)) + else + noise_ = noise + end + if isa(noise_rate_prototype, Void) + if n > 1 || nc > 1 || (n > 0 && nc > 0) + noise_rate_prototype = zeros(Complex128, length(x0), n + nc) + end + end + out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)}) out = DiffEqCallbacks.SavedValues(Float64,out_type) @@ -126,32 +143,17 @@ function integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Function, x0: save_everystep=save_everystep, save_start = false) - full_cb = OrdinaryDiffEq.CallbackSet(callback, scb) + full_cb = OrdinaryDiffEq.CallbackSet(callback, ncb, scb) - if n == 1 - prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]), - noise=noise) - else - prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]), - noise=noise, noise_rate_prototype=noise_rate_prototype) - end - - if isa(alg, Void) - if n == 1 - alg_ = StochasticDiffEq.RKMil(interpretation=:Stratonovich) - else - alg_ = StochasticDiffEq.LambaEulerHeun() - end - else - @assert isa(alg, StochasticDiffEq.StochasticDiffEqAlgorithm) - alg_ = alg - end + prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]), + noise=noise_, + noise_rate_prototype=noise_rate_prototype) sol = StochasticDiffEq.solve( prob, - alg_; - reltol = 1.0e-4, - abstol = 1.0e-4, + alg; + reltol = 1.0e-3, + abstol = 1.0e-3, save_everystep = false, save_start = false, save_end = false, callback=full_cb, kwargs...) diff --git a/test/runtests.jl b/test/runtests.jl index 99bd0fa0..10735496 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,7 @@ names = [ "test_manybody.jl", "test_nlevel.jl", "test_subspace.jl", + "test_state_definitions.jl", "test_phasespace.jl", "test_transformations.jl", @@ -40,6 +41,7 @@ names = [ "test_semiclassical.jl", + "test_stochastic_definitions.jl", "test_stochastic_schroedinger.jl", "test_stochastic_master.jl", "test_stochastic_semiclassical.jl", diff --git a/test/test_metrics.jl b/test/test_metrics.jl index ced476ab..9d6dfe84 100644 --- a/test/test_metrics.jl +++ b/test/test_metrics.jl @@ -56,6 +56,8 @@ rho = spinup(b1) ⊗ dagger(coherentstate(b2, 0.1)) # entropy rho_mix = full(identityoperator(b1))/2. @test entropy_vn(rho_mix)/log(2) ≈ 1 +psi = coherentstate(FockBasis(20), 2.0) +@test isapprox(entropy_vn(psi), 0.0, atol=1e-8) # fidelity rho = tensor(psi1, dagger(psi1)) diff --git a/test/test_operators.jl b/test/test_operators.jl index 976c8ec1..0c98438a 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -68,4 +68,6 @@ op_test3 = test_operators(b1 ⊗ b2, b2 ⊗ b1, randoperator(b, b).data) @test one(b1).data == diagm(ones(b1.shape[1])) @test one(op1).data == diagm(ones(b1.shape[1])) +@test_throws ArgumentError conj!(op_test) + end # testset diff --git a/test/test_operators_dense.jl b/test/test_operators_dense.jl index 1f26132e..73c4448e 100644 --- a/test/test_operators_dense.jl +++ b/test/test_operators_dense.jl @@ -72,6 +72,10 @@ xbra2 = Bra(b_l, rand(Complex128, length(b_l))) @test 1e-12 > D((op1 + op2)*dagger(0.3*op3), 0.3*op1*dagger(op3) + 0.3*op2*dagger(op3)) @test 1e-12 > D(0.3*(op1*dagger(op2)), op1*(0.3*dagger(op2))) +tmp = copy(op1) +conj!(tmp) +@test tmp == conj(op1) && conj(tmp.data) == op1.data + # Internal layout b1 = GenericBasis(2) b2 = GenericBasis(3) diff --git a/test/test_operators_sparse.jl b/test/test_operators_sparse.jl index 350669d5..cc90aa96 100644 --- a/test/test_operators_sparse.jl +++ b/test/test_operators_sparse.jl @@ -87,6 +87,11 @@ xbra2 = dagger(randstate(b_l)) # Test division @test 1e-14 > D(op1/7, op1_/7) +# Conjugation +tmp = copy(op1) +conj!(tmp) +@test tmp == conj(op1) && conj(tmp.data) == op1.data + # Test identityoperator Idense = identityoperator(DenseOperator, b_r) I = identityoperator(SparseOperator, b_r) diff --git a/test/test_phasespace.jl b/test/test_phasespace.jl index cca71aba..e5db5a93 100644 --- a/test/test_phasespace.jl +++ b/test/test_phasespace.jl @@ -85,4 +85,15 @@ for (i,x)=enumerate(X), (j,y)=enumerate(Y) @test 1e-14 > abs(qfunc(rho, c) - q_rho) end +# Test SU(2) phasespace +b = SpinBasis(50) +theta = π*rand() +phi =2π*rand() +css = coherentspinstate(b,theta,phi) +sx = expect(sigmax(b)/2,css); # eigenstate of jx operator +sy = expect(sigmay(b)/2,css); # eigenstate of jy +sz = expect(sigmaz(b)/2,css); # eigenstate of jz operator +ssq = sx^2 + sy^2 + sz^2 + +@test ssq ≈ float(b.spinnumber)^2 end # testset diff --git a/test/test_spectralanalysis.jl b/test/test_spectralanalysis.jl index 4411aeb2..ebff6831 100644 --- a/test/test_spectralanalysis.jl +++ b/test/test_spectralanalysis.jl @@ -1,8 +1,7 @@ using Base.Test using QuantumOptics -mutable struct SpectralanalysisTestOperator <: Operator -end +mutable struct SpectralanalysisTestOperator <: Operator end @testset "spectralanalysis" begin @@ -26,15 +25,15 @@ d = [-3, -2.6, -0.1, 0.0, 0.6] D = DenseOperator(b, diagm(d)) Dsp = sparse(D) @test eigenenergies(D) ≈ d -@test eigenenergies(Dsp, 3) ≈ d[1:3] +@test eigenenergies(Dsp, 3, info=false) ≈ d[1:3] R = U*D*dagger(U) Rsp = sparse(R) @test eigenenergies((R+dagger(R))/2) ≈ d -@test eigenenergies((Rsp+dagger(Rsp))/2, 3) ≈ d[1:3] +@test eigenenergies((Rsp+dagger(Rsp))/2, 3; info=false) ≈ d[1:3] E, states = eigenstates((R+dagger(R))/2, 3) -Esp, states_sp = eigenstates((Rsp+dagger(Rsp))/2, 3) +Esp, states_sp = eigenstates((Rsp+dagger(Rsp))/2, 3, info=false) for i=1:3 @test E[i] ≈ d[i] @test Esp[i] ≈ d[i] @@ -54,15 +53,15 @@ d = [-3+0.2im, -2.6-0.1im, -0.1+0.5im, 0.0, 0.6+0.3im] D = DenseOperator(b, diagm(d)) Dsp = sparse(D) @test eigenenergies(D; warning=false) ≈ d -@test eigenenergies(Dsp, 3; warning=false) ≈ d[1:3] +@test eigenenergies(Dsp, 3; warning=false, info=false) ≈ d[1:3] R = U*D*dagger(U) Rsp = sparse(R) @test eigenenergies(R; warning=false) ≈ d -@test eigenenergies(Rsp, 3; warning=false) ≈ d[1:3] +@test eigenenergies(Rsp, 3; warning=false, info=false) ≈ d[1:3] E, states = eigenstates(R, 3; warning=false) -Esp, states_sp = eigenstates(Rsp, 3; warning=false) +Esp, states_sp = eigenstates(Rsp, 3; warning=false, info=false) for i=1:3 @test E[i] ≈ d[i] @test Esp[i] ≈ d[i] diff --git a/test/test_state_definitions.jl b/test/test_state_definitions.jl new file mode 100644 index 00000000..58a173d9 --- /dev/null +++ b/test/test_state_definitions.jl @@ -0,0 +1,40 @@ +using Base.Test +using QuantumOptics + +@testset "state_definitions" begin + +n=30 +b=FockBasis(n) +omega=40.3 +T=2.3756 +r=thermalstate(omega*number(b),T) +for k=1:n-1 + @test isapprox(r.data[k+1,k+1]/r.data[k,k],exp(-omega/T)) +end +S=entropy_vn(r) +z=sum(exp.(-[0:n;]*omega)) +s=expect(omega*number(b),r)/T+log(z) +isapprox(S,s) + +alpha=rand()+im*rand() +r=coherentthermalstate(b,omega*number(b),T,alpha) +@test isapprox(expect(number(b),r),abs(alpha)^2+1/(exp(omega/T)-1)) +@test isapprox(expect(destroy(b),r),alpha) +@test isapprox(entropy_vn(r),S) + +rp=phase_average(r) +@test isapprox(expect(number(b),rp),abs(alpha)^2+1/(exp(omega/T)-1)) +@test isapprox(expect(destroy(b),rp),0) +for k=1:n + @test isapprox(rp.data[k,k],r.data[k,k]) +end + +rpas=passive_state(r) +for k=1:n-1 + @test real(rpas.data[k+1,k+1])