Implement Stochastic Schrodinger equation (#238)

* Implement Stochastic Schrodinger equation

* Minor changes

* Minor changes on the import

* Fix wrong link

* Fix instabilities

* Avoid safety_copy

Author: albertomercurio

Project.toml

@@ -7,6 +7,7 @@ version = "0.14.1"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 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"
@@ -22,6 +23,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -34,6 +36,7 @@ ArrayInterface = "6, 7"
 CUDA = "5"
 DiffEqBase = "6"
 DiffEqCallbacks = "2 - 3.1, 3.8"
+DiffEqNoiseProcess = "5"
 FFTW = "1.5"
 Graphs = "1.7"
 IncompleteLU = "0.2"
@@ -49,6 +52,7 @@ SciMLOperators = "0.3"
 SparseArrays = "<0.0.1, 1"
 SpecialFunctions = "2"
 StaticArraysCore = "1"
+StochasticDiffEq = "6"
 Test = "<0.0.1, 1"
 julia = "1.10"
 

docs/src/api.md

@@ -175,18 +175,22 @@ qeye
 ```@docs
 TimeEvolutionSol
 TimeEvolutionMCSol
+TimeEvolutionSSESol
 sesolveProblem
 mesolveProblem
-lr_mesolveProblem
 mcsolveProblem
 mcsolveEnsembleProblem
+ssesolveProblem
+ssesolveEnsembleProblem
+lr_mesolveProblem
 sesolve
 mesolve
-lr_mesolve
 mcsolve
+ssesolve
 dfd_mesolve
 dsf_mesolve
 dsf_mcsolve
+lr_mesolve
 liouvillian
 liouvillian_generalized
 steadystate

src/QuantumToolbox.jl

@@ -22,18 +22,21 @@ import SciMLBase:
     remake,
     u_modified!,
     ODEProblem,
+    SDEProblem,
     EnsembleProblem,
     EnsembleThreads,
     FullSpecialize,
     CallbackSet,
     ContinuousCallback,
     DiscreteCallback
+import StochasticDiffEq: StochasticDiffEqAlgorithm, SRA1
 import SciMLOperators: MatrixOperator
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_MINRES, KrylovJL_GMRES
 import DiffEqBase: get_tstops
 import DiffEqCallbacks: PeriodicCallback, PresetTimeCallback, TerminateSteadyState
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
 import OrdinaryDiffEqTsit5: Tsit5
+import DiffEqNoiseProcess: RealWienerProcess
 
 # other dependencies (in alphabetical order)
 import ArrayInterface: allowed_getindex, allowed_setindex!
@@ -73,6 +76,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)