Stochastic Schrödinger equation

Project.toml

@@ -6,6 +6,7 @@ version = "0.12.0"
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
@@ -18,6 +19,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -29,6 +31,7 @@ QuantumToolboxCUDAExt = "CUDA"
 ArrayInterface = "6, 7"
 CUDA = "5"
 DiffEqCallbacks = "2 - 3.1, 3.8"
+DiffEqNoiseProcess = "5"
 FFTW = "1.5"
 Graphs = "1.7"
 IncompleteLU = "0.2"
@@ -41,6 +44,7 @@ Random = "<0.0.1, 1"
 Reexport = "1"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "2"
+StochasticDiffEq = "6"
 Test = "<0.0.1, 1"
 julia = "1.7"
 

docs/src/api.md

@@ -179,11 +179,14 @@ sesolveProblem
 mesolveProblem
 lr_mesolveProblem
 mcsolveProblem
+ssesolveProblem
 mcsolveEnsembleProblem
+ssesolveEnsembleProblem
 sesolve
 mesolve
 lr_mesolve
 mcsolve
+ssesolve
 dfd_mesolve
 dsf_mesolve
 dsf_mcsolve

src/QuantumToolbox.jl

@@ -7,6 +7,7 @@ import Reexport: @reexport
 @reexport using LinearAlgebra
 @reexport using SparseArrays
 @reexport using OrdinaryDiffEq
+@reexport using StochasticDiffEq
 @reexport using LinearSolve
 
 # other functions in LinearAlgebra
@@ -24,6 +25,7 @@ end
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
 import DiffEqCallbacks: DiscreteCallback, PeriodicCallback, PresetTimeCallback, TerminateSteadyState
+import DiffEqNoiseProcess: RealWienerProcess, RealWienerProcess!
 import FFTW: fft, fftshift
 import Graphs: connected_components, DiGraph
 import IncompleteLU: ilu
@@ -60,6 +62,7 @@ include("time_evolution/mesolve.jl")
 include("time_evolution/lr_mesolve.jl")
 include("time_evolution/sesolve.jl")
 include("time_evolution/mcsolve.jl")
+include("time_evolution/ssesolve.jl")
 include("time_evolution/time_evolution_dynamical.jl")
 
 # Others

src/qobj/operator_sum.jl

@@ -18,7 +18,7 @@ struct OperatorSum{CT<:Vector{<:Number},OT<:AbstractVector} <: AbstractQuantumOb
         mapreduce(x -> size(x) == size_1, &, operators) ||
             throw(DimensionMismatch("All the operators must have the same dimensions."))
         T = promote_type(
-            mapreduce(x -> eltype(x.data), promote_type, operators),
+            mapreduce(x -> eltype(x), promote_type, operators),
             mapreduce(eltype, promote_type, coefficients),
         )
         coefficients2 = T.(coefficients)

src/time_evolution/mcsolve.jl

@@ -86,7 +86,7 @@ function _mcsolve_output_func(sol, i)
     return (sol, false)
 end
 
-function _mcsolve_generate_statistics(sol, i, times, states, expvals_all, jump_times, jump_which)
+function _mcsolve_generate_statistics!(sol, i, times, states, expvals_all, jump_times, jump_which)
     sol_i = sol[:, i]
     !isempty(sol_i.prob.kwargs[:saveat]) ?
     states[i] = [QuantumObject(sol_i.u[i], dims = sol_i.prob.p.Hdims) for i in 1:length(sol_i.u)] : nothing
@@ -524,7 +524,7 @@ function mcsolve(
     jump_which = Vector{Vector{Int16}}(undef, length(sol))
 
     foreach(
-        i -> _mcsolve_generate_statistics(sol, i, times, states, expvals_all, jump_times, jump_which),
+        i -> _mcsolve_generate_statistics!(sol, i, times, states, expvals_all, jump_times, jump_which),
         eachindex(sol),
     )
     expvals = dropdims(sum(expvals_all, dims = 1), dims = 1) ./ length(sol)