Add smesolve solver

.typos.toml

@@ -1,2 +1,4 @@
 [default.extend-words]
-ket = "ket"
+ket = "ket"
+sme = "sme"
+SME = "SME"

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Fix CUDA `sparse_to_dense`. ([#386])
 - Improve pseudo inverse spectrum solver. ([#388])
+- Add `smesolve` function for stochastic master equation. ([#389])
 
 ## [v0.25.2]
 Release date: 2025-02-02
@@ -109,3 +110,4 @@ Release date: 2024-11-13
 [#383]: https://github.com/qutip/QuantumToolbox.jl/issues/383
 [#386]: https://github.com/qutip/QuantumToolbox.jl/issues/386
 [#388]: https://github.com/qutip/QuantumToolbox.jl/issues/388
+[#389]: https://github.com/qutip/QuantumToolbox.jl/issues/389

docs/src/resources/api.md

@@ -196,10 +196,13 @@ mcsolveProblem
 mcsolveEnsembleProblem
 ssesolveProblem
 ssesolveEnsembleProblem
+smesolveProblem
+smesolveEnsembleProblem
 sesolve
 mesolve
 mcsolve
 ssesolve
+smesolve
 dfd_mesolve
 liouvillian
 liouvillian_generalized

src/QuantumToolbox.jl

@@ -106,6 +106,7 @@ include("time_evolution/lr_mesolve.jl")
 include("time_evolution/sesolve.jl")
 include("time_evolution/mcsolve.jl")
 include("time_evolution/ssesolve.jl")
+include("time_evolution/smesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
 
 # Others

src/qobj/superoperators.jl

@@ -164,22 +164,34 @@
 ) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, Id_cache)
     if !(c_ops isa Nothing)
-        # sum all the (time-independent) c_ops first
-        c_ops_ti = filter(op -> isa(op, QuantumObject), c_ops)
-        if !isempty(c_ops_ti)
-            L += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_ti)
-        end
-
-        # sum rest of the QobjEvo together
-        c_ops_td = filter(op -> isa(op, QuantumObjectEvolution), c_ops)
-        if !isempty(c_ops_td)
-            L += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_td)
-        end
+        L += _sum_lindblad_dissipators(c_ops, Id_cache)
     end
     return L
 end
 
+liouvillian(H::Nothing, c_ops::Union{AbstractVector,Tuple}, Id_cache::Diagonal = I(prod(c_ops[1].dims))) =
+    _sum_lindblad_dissipators(c_ops, Id_cache)
+
+liouvillian(H::Nothing, c_ops::Nothing) = 0
+
 liouvillian(H::AbstractQuantumObject{OperatorQuantumObject}, Id_cache::Diagonal = I(prod(H.dimensions))) =
     get_typename_wrapper(H)(_liouvillian(H.data, Id_cache), SuperOperator, H.dimensions)
 
 liouvillian(H::AbstractQuantumObject{SuperOperatorQuantumObject}, Id_cache::Diagonal) = H
+
+function _sum_lindblad_dissipators(c_ops, Id_cache::Diagonal)
+    D = 0
+    # sum all the (time-independent) c_ops first
+    c_ops_ti = filter(op -> isa(op, QuantumObject), c_ops)
+    if !isempty(c_ops_ti)
+        D += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_ti)
+    end
+
+    # sum rest of the QobjEvo together
+    c_ops_td = filter(op -> isa(op, QuantumObjectEvolution), c_ops)
+    if !isempty(c_ops_td)
+        D += mapreduce(op -> lindblad_dissipator(op, Id_cache), +, c_ops_td)
+    end
+
+    return D
+end

src/time_evolution/mcsolve.jl

@@ -12,11 +12,6 @@ function _mcsolve_prob_func(prob, i, repeat, global_rng, seeds, tlist)
     return remake(prob, f = f, callback = cb)
 end
 
-function _mcsolve_dispatch_prob_func(rng, ntraj, tlist)
-    seeds = map(i -> rand(rng, UInt64), 1:ntraj)
-    return (prob, i, repeat) -> _mcsolve_prob_func(prob, i, repeat, rng, seeds, tlist)
-end
-
 # Standard output function
 function _mcsolve_output_func(sol, i)
     idx = _mc_get_jump_callback(sol).affect!.jump_times_which_idx[]
@@ -25,43 +20,6 @@ function _mcsolve_output_func(sol, i)
     return (sol, false)
 end
 
-# Output function with progress bar update
-function _mcsolve_output_func_progress(sol, i, progr)
-    next!(progr)
-    return _mcsolve_output_func(sol, i)
-end
-
-# Output function with distributed channel update for progress bar
-function _mcsolve_output_func_distributed(sol, i, channel)
-    put!(channel, true)
-    return _mcsolve_output_func(sol, i)
-end
-
-function _mcsolve_dispatch_output_func(::ET, progress_bar, ntraj) where {ET<:Union{EnsembleSerial,EnsembleThreads}}
-    if getVal(progress_bar)
-        progr = ProgressBar(ntraj, enable = getVal(progress_bar))
-        f = (sol, i) -> _mcsolve_output_func_progress(sol, i, progr)
-        return (f, progr, nothing)
-    else
-        return (_mcsolve_output_func, nothing, nothing)
-    end
-end
-function _mcsolve_dispatch_output_func(
-    ::ET,
-    progress_bar,
-    ntraj,
-) where {ET<:Union{EnsembleSplitThreads,EnsembleDistributed}}
-    if getVal(progress_bar)
-        progr = ProgressBar(ntraj, enable = getVal(progress_bar))
-        progr_channel::RemoteChannel{Channel{Bool}} = RemoteChannel(() -> Channel{Bool}(1))
-
-        f = (sol, i) -> _mcsolve_output_func_distributed(sol, i, progr_channel)
-        return (f, progr, progr_channel)
-    else
-        return (_mcsolve_output_func, nothing, nothing)
-    end
-end
-
 function _normalize_state!(u, dims, normalize_states)
     getVal(normalize_states) && normalize!(u)
     return QuantumObject(u, type = Ket, dims = dims)
@@ -280,9 +238,10 @@ function mcsolveEnsembleProblem(
     output_func::Union{Tuple,Nothing} = nothing,
     kwargs...,
 ) where {TJC<:LindbladJumpCallbackType}
-    _prob_func = prob_func isa Nothing ? _mcsolve_dispatch_prob_func(rng, ntraj, tlist) : prob_func
+    _prob_func = isnothing(prob_func) ? _ensemble_dispatch_prob_func(rng, ntraj, tlist, _mcsolve_prob_func) : prob_func
     _output_func =
-        output_func isa Nothing ? _mcsolve_dispatch_output_func(ensemble_method, progress_bar, ntraj) : output_func
+        output_func isa Nothing ?
+        _ensemble_dispatch_output_func(ensemble_method, progress_bar, ntraj, _mcsolve_output_func) : output_func
 
     prob_mc = mcsolveProblem(
         H,
@@ -431,46 +390,14 @@ function mcsolve(
     return mcsolve(ens_prob_mc, alg, ntraj, ensemble_method, normalize_states)
 end
 
-function _mcsolve_solve_ens(
-    ens_prob_mc::TimeEvolutionProblem,
-    alg::OrdinaryDiffEqAlgorithm,
-    ensemble_method::ET,
-    ntraj::Int,
-) where {ET<:Union{EnsembleSplitThreads,EnsembleDistributed}}
-    sol = nothing
-
-    @sync begin
-        @async while take!(ens_prob_mc.kwargs.channel)
-            next!(ens_prob_mc.kwargs.progr)
-        end
-
-        @async begin
-            sol = solve(ens_prob_mc.prob, alg, ensemble_method, trajectories = ntraj)
-            put!(ens_prob_mc.kwargs.channel, false)
-        end
-    end
-
-    return sol
-end
-
-function _mcsolve_solve_ens(
-    ens_prob_mc::TimeEvolutionProblem,
-    alg::OrdinaryDiffEqAlgorithm,
-    ensemble_method,
-    ntraj::Int,
-)
-    sol = solve(ens_prob_mc.prob, alg, ensemble_method, trajectories = ntraj)
-    return sol
-end
-
 function mcsolve(
     ens_prob_mc::TimeEvolutionProblem,
     alg::OrdinaryDiffEqAlgorithm = Tsit5(),
     ntraj::Int = 1,
     ensemble_method = EnsembleThreads(),
     normalize_states = Val(true),
 )
-    sol = _mcsolve_solve_ens(ens_prob_mc, alg, ensemble_method, ntraj)
+    sol = _ensemble_dispatch_solve(ens_prob_mc, alg, ensemble_method, ntraj)
 
     dims = ens_prob_mc.dimensions
     _sol_1 = sol[:, 1]

src/time_evolution/mesolve.jl

@@ -1,7 +1,9 @@
 export mesolveProblem, mesolve
 
-_mesolve_make_L_QobjEvo(H::QuantumObject, c_ops) = QobjEvo(liouvillian(H, c_ops); type = SuperOperator)
+_mesolve_make_L_QobjEvo(H::Union{QuantumObject,Nothing}, c_ops) = QobjEvo(liouvillian(H, c_ops); type = SuperOperator)
 _mesolve_make_L_QobjEvo(H::Union{QuantumObjectEvolution,Tuple}, c_ops) = liouvillian(QobjEvo(H), c_ops)
+_mesolve_make_L_QobjEvo(H::Nothing, c_ops::Nothing) = throw(ArgumentError("Both H and
+c_ops are Nothing. You are probably running the wrong function."))
 
 @doc raw"""
     mesolveProblem(
@@ -31,7 +33,7 @@
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`KetQuantumObject`](@ref) or a [`OperatorQuantumObject`](@ref) for mixed states.
 - `tlist`: List of times at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `e_ops`: List of operators for which to calculate expectation values. It can be either a `Vector` or a `Tuple`.
@@ -119,7 +121,7 @@
 # Arguments
 
 - `H`: Hamiltonian of the system ``\hat{H}``. It can be either a [`QuantumObject`](@ref), a [`QuantumObjectEvolution`](@ref), or a `Tuple` of operator-function pairs.
-- `ψ0`: Initial state of the system ``|\psi(0)\rangle``.
+- `ψ0`: Initial state of the system ``|\psi(0)\rangle``. It can be either a [`KetQuantumObject`](@ref) or a [`OperatorQuantumObject`](@ref) for mixed states.
 - `tlist`: List of times at which to save either the state or the expectation values of the system.
 - `c_ops`: List of collapse operators ``\{\hat{C}_n\}_n``. It can be either a `Vector` or a `Tuple`.
 - `alg`: The algorithm for the ODE solver. The default value is `Tsit5()`.