update

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

src/QuantumOptics.jl

@@ -12,19 +12,19 @@ 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,
                 position, momentum, potentialoperator, transform,
         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

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...)

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)

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)
 
 """

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)

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)

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

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

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
+