Merge pull request #46

CompatHelper: bump compat for SteadyStateDiffEq to 2, (keep existing compat)

.github/workflows/Runtests.yml

@@ -8,6 +8,7 @@ on:
     branches:
       - 'main'
     types:
+      - 'opened'
       - 'synchronize'
       - 'ready_for_review'
 
@@ -21,7 +22,6 @@ jobs:
         # for Core functionalities
         julia-version:
           - '1'
-          - '1.8'
         os:
           - ubuntu-latest
           - macOS-latest

.github/workflows/documentation.yml

@@ -10,6 +10,7 @@ on:
     branches:
       - 'main'
     types:
+      - 'opened'
       - 'synchronize'
       - 'ready_for_review'
 

Project.toml

@@ -1,10 +1,11 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 authors = ["Yi-Te Huang <[email protected]>"]
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 Crayons = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -28,19 +29,20 @@ HierarchicalEOM_QOExt = "QuantumOptics"
 
 [compat]
 Crayons = "4.1"
+DiffEqBase = "6.140"
 FastExpm = "1.1.0"
 JLD2 = "0.4.31"
-LinearAlgebra = "1.8"
+LinearAlgebra = "1.9"
 LinearSolve = "2.4.2"
 OrdinaryDiffEq = "6.53.4"
-Pkg = "1.8"
+Pkg = "1.9"
 PrecompileTools = "1"
 ProgressMeter = "1.7"
 Reexport = "1"
 SciMLOperators = "0.3.5"
-SparseArrays = "1.8"
-SteadyStateDiffEq = "1.16.0"
-julia = "1.8"
+SparseArrays = "1.9"
+SteadyStateDiffEq = "2.0"
+julia = "1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -42,10 +42,10 @@ Pkg.add("HierarchicalEOM")
 ```
 Alternatively, this can also be done in Julia's [Pkg REPL](https://julialang.github.io/Pkg.jl/v1/getting-started/) by pressing the key `]` in the REPL to use the package mode, and then type the following command:
 ```julia-REPL
-(1.8) pkg> add HierarchicalEOM
+(1.9) pkg> add HierarchicalEOM
 ```
 More information about `Julia`'s package manager can be found at [`Pkg.jl`](https://julialang.github.io/Pkg.jl/v1/).  
-`HierarchicalEOM.jl` requires Julia 1.8 or higher. Installing it on an older version of Julia will result in many errors.
+`HierarchicalEOM.jl` now requires Julia 1.9 or higher. Installing it on an older version of Julia will result in many errors.
 
 To load the package and check the version information, use the command:
 ```julia

docs/src/install.md

@@ -8,11 +8,11 @@ Pkg.add("HierarchicalEOM")
 ```
 Alternatively, this can also be done in Julia's [Pkg REPL](https://julialang.github.io/Pkg.jl/v1/getting-started/) by pressing the key `]` in the REPL to use the package mode, and then type the following command:
 ```julia-REPL
-(1.8) pkg> add HierarchicalEOM
+(1.9) pkg> add HierarchicalEOM
 ```
 More information about `Julia`'s package manager can be found at [`Pkg.jl`](https://julialang.github.io/Pkg.jl/v1/).  
-!!! note "Julia 1.8"
-    `HierarchicalEOM.jl` requires Julia 1.8 or higher
+!!! note "Julia 1.9"
+    `HierarchicalEOM.jl` requires Julia 1.9 or higher (we dropped Julia 1.8 since ver.1.12.1)
 
 To load the package and check the version information, use the command:
 ```julia

docs/src/stationary_state.md

@@ -26,7 +26,7 @@ ados_steady = SteadyState(M)
     This method does not require an initial condition ``\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(0)``. Although this method works for most of the cases, it does not guarantee that one can obtain a physical (or unique) solution. If there is any problem within the solution, please try the second method which solves with an initial condition, as shown below.
 
 ## Solve with [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/)
-The second method is implemented by solving the ordinary differential equation (ODE)
+The second method is implemented by solving the ordinary differential equation (ODE) method : [`SteadyStateDiffEq.jl`](https://github.com/SciML/SteadyStateDiffEq.jl)
 ```math
 \partial_{t}\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)=\hat{\mathcal{M}}\rho^{(m,n,p)}_{\textbf{j} \vert \textbf{q}}(t)
 ```
@@ -39,7 +39,7 @@ until finding a stationary solution.
 See the docstring of this method:  
 
 ```@docs
-SteadyState(M::AbstractHEOMLSMatrix, ρ0; solver = DP5(), reltol::Real = 1.0e-6, abstol::Real = 1.0e-8, maxiters::Real = 1e5, save_everystep::Bool=false, verbose::Bool = true, SOLVEROptions...)
+SteadyState(M::AbstractHEOMLSMatrix, ρ0, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
 ```
 
 ```julia
@@ -55,7 +55,7 @@ ados_steady = SteadyState(M, ρ0)
 ### Given the initial state as Auxiliary Density Operators
 See the docstring of this method:  
 ```@docs
-SteadyState(M::AbstractHEOMLSMatrix, ados::ADOs; solver = DP5(), reltol::Real = 1.0e-6, abstol::Real = 1.0e-8, maxiters::Real = 1e5, save_everystep::Bool=false, verbose::Bool = true, SOLVEROptions...)
+SteadyState(M::AbstractHEOMLSMatrix, ados::ADOs, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
 ```
 
 ```julia

src/HierarchicalEOM.jl

@@ -62,6 +62,7 @@ module HierarchicalEOM
         # for solving steady state
         import LinearSolve: LinearProblem, init, solve!, UMFPACKFactorization
         import OrdinaryDiffEq: SteadyStateProblem, solve
+        import DiffEqBase: NormTerminationMode
         import SteadyStateDiffEq: DynamicSS
 
         export

src/SteadyState.jl

@@ -25,7 +25,7 @@ For more details about solvers and extra options, please refer to [`LinearSolve.
 
     # solving x where A * x = b
     if verbose
-        print("Solving steady state for auxiliary density operators...")
+        print("Solving steady state for ADOs by linear-solve method...")
         flush(stdout)
     end
     cache = init(LinearProblem(A, _HandleVectorType(typeof(M.data), b)), solver, SOLVEROptions...)
@@ -39,79 +39,67 @@ For more details about solvers and extra options, please refer to [`LinearSolve.
 end
 
 @doc raw"""
-    SteadyState(M, ρ0; solver, reltol, abstol, maxiters, save_everystep, verbose, SOLVEROptions...)
-Solve the steady state of the auxiliary density operators based on time evolution (ordinary differential equations)
-with initial state is given in the type of density-matrix (`ρ0`).
+    SteadyState(M, ρ0, tspan; solver, termination_condition, verbose, SOLVEROptions...)
+Solve the steady state of the auxiliary density operators based on time evolution (ordinary differential equations; `SteadyStateDiffEq.jl`) with initial state is given in the type of density-matrix (`ρ0`).
 
 # Parameters
 - `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model, where the parity should be `EVEN`.
 - `ρ0` : system initial state (density matrix)
+- `tspan::Number` : the time limit to find stationary state. Default to `Inf`
 - `solver` : The ODE solvers in package `DifferentialEquations.jl`. Default to `DP5()`.
-- `reltol::Real` : Relative tolerance in adaptive timestepping. Default to `1.0e-6`.
-- `abstol::Real` : Absolute tolerance in adaptive timestepping. Default to `1.0e-8`.
-- `maxiters::Real` : Maximum number of iterations before stopping. Default to `1e5`.
-- `save_everystep::Bool` : Saves the result at every step. Defaults to `false`.
+- `termination_condition` : The stationary state terminate condition in `DiffEqBase.jl`. Default to `NormTerminationMode()`.
 - `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`.
 - `SOLVEROptions` : extra options for solver
 
-For more details about solvers and extra options, please refer to [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/)
+For more details about solvers, termination condition, and extra options, please refer to [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/)
 
 # Returns
 - `::ADOs` : The steady state of auxiliary density operators.
 """
 function SteadyState(
         M::AbstractHEOMLSMatrix, 
-        ρ0;
+        ρ0,
+        tspan::Number = Inf;
         solver = DP5(),
-        reltol::Real = 1.0e-6,
-        abstol::Real = 1.0e-8,
-        maxiters::Real = 1e5,
-        save_everystep::Bool=false,
+        termination_condition = NormTerminationMode(),
         verbose::Bool = true,
         SOLVEROptions...
     )
     return SteadyState(
         M, 
-        ADOs(ρ0, M.N, M.parity);
+        ADOs(ρ0, M.N, M.parity),
+        tspan;
         solver = solver,
-        reltol = reltol,
-        abstol = abstol,
-        maxiters = maxiters,
-        save_everystep = save_everystep,
+        termination_condition = termination_condition,
         verbose = verbose,
         SOLVEROptions...
     )
 end
 
 @doc raw"""
-    SteadyState(M, ados; solver, reltol, abstol, maxiters, save_everystep, verbose, SOLVEROptions...)
-Solve the steady state of the auxiliary density operators based on time evolution (ordinary differential equations)
-with initial state is given in the type of `ADOs`.
+    SteadyState(M, ados, tspan; solver, termination_condition, verbose, SOLVEROptions...)
+Solve the steady state of the auxiliary density operators based on time evolution (ordinary differential equations; `SteadyStateDiffEq.jl`) with initial state is given in the type of `ADOs`.
 
 # Parameters
 - `M::AbstractHEOMLSMatrix` : the matrix given from HEOM model, where the parity should be `EVEN`.
 - `ados::ADOs` : initial auxiliary density operators
+- `tspan::Number` : the time limit to find stationary state. Default to `Inf`
 - `solver` : The ODE solvers in package `DifferentialEquations.jl`. Default to `DP5()`.
-- `reltol::Real` : Relative tolerance in adaptive timestepping. Default to `1.0e-3`.
-- `abstol::Real` : Absolute tolerance in adaptive timestepping. Default to `1.0e-6`.
-- `maxiters::Real` : Maximum number of iterations before stopping. Default to `1e5`.
-- `save_everystep::Bool` : Saves the result at every step. Defaults to `false`.
+- `termination_condition` : The stationary state terminate condition in `DiffEqBase.jl`. Default to `NormTerminationMode()`.
 - `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`.
 - `SOLVEROptions` : extra options for solver
 
-For more details about solvers and extra options, please refer to [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/)
+For more details about solvers, termination condition, and extra options, please refer to [`DifferentialEquations.jl`](https://diffeq.sciml.ai/stable/)
 
 # Returns
 - `::ADOs` : The steady state of auxiliary density operators.
 """
 @noinline function SteadyState(
         M::AbstractHEOMLSMatrix, 
-        ados::ADOs;
+        ados::ADOs,
+        tspan::Number = Inf;
         solver = DP5(),
-        reltol = 1.0e-6,
-        abstol = 1.0e-8,
-        maxiters = 1e5,
-        save_everystep::Bool = false,
+        termination_condition = NormTerminationMode(),
         verbose::Bool = true,
         SOLVEROptions...
     )
@@ -121,20 +109,17 @@ For more details about solvers and extra options, please refer to [`Differential
 
     # problem: dρ(t)/dt = L * ρ(t)
     L = MatrixOperator(M.data)
-    ElType = eltype(M)
-    tspan  = Tuple{ElType, ElType}((0, Inf))
-    prob   = ODEProblem(L, _HandleVectorType(typeof(M.data), ados.data), tspan)
+    prob = SteadyStateProblem(L, _HandleVectorType(typeof(M.data), ados.data))
 
     # solving steady state of the ODE problem
     if verbose
-        print("Solving steady state for auxiliary density operators...")
+        print("Solving steady state for ADOs by Ordinary Differential Equations method...")
         flush(stdout)
     end
     sol = solve(
-        SteadyStateProblem(prob), 
-        DynamicSS(solver; abstol = _HandleFloatType(eltype(M), abstol), reltol = _HandleFloatType(eltype(M), reltol));
-        maxiters = maxiters,
-        save_everystep = save_everystep,
+        prob, 
+        DynamicSS(solver, _HandleFloatType(eltype(M), tspan));
+        termination_condition = termination_condition,
         SOLVEROptions...
     )
 

src/evolution.jl

@@ -107,7 +107,7 @@ For more details, please refer to [`FastExpm.jl`](https://github.com/fmentink/Fa
     # start solving
     ρvec = copy(ados.data)
     if verbose
-        print("Solving time evolution for auxiliary density operators...\n")
+        print("Solving time evolution for ADOs by propagator method...\n")
         flush(stdout)
         prog = Progress(steps + 1; start=1, desc="Progress : ", PROGBAR_OPTIONS...)
     end
@@ -259,7 +259,7 @@ For more details about solvers and extra options, please refer to [`Differential
 
     # start solving ode
     if verbose
-        print("Solving time evolution for auxiliary density operators...\n")
+        print("Solving time evolution for ADOs by Ordinary Differential Equations method...\n")
         flush(stdout)
         prog = Progress(length(Tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...)
     end
@@ -427,7 +427,7 @@ For more details about solvers and extra options, please refer to [`Differential
 
     # start solving ode
     if verbose
-        print("Solving time evolution for auxiliary density operators with time-dependent Hamiltonian...\n")
+        print("Solving time evolution for ADOs with time-dependent Hamiltonian by Ordinary Differential Equations method...\n")
         flush(stdout)
         prog = Progress(length(Tlist); start=1, desc="Progress : ", PROGBAR_OPTIONS...)
     end

src/precompile.jl

@@ -3,9 +3,11 @@ import PrecompileTools
 PrecompileTools.@setup_workload begin
     # Putting some things in `setup` can reduce the size of the
     # precompile file and potentially make loading faster.
-    op = [0.1 0.2; 0.2 0.4]
-    bB = Boson_DrudeLorentz_Pade(op, 1, 1., 1., 3)
-    fB = Fermion_Lorentz_Pade(op, 1., 1., 1., 1., 2)
+    d  = [0 1; 0 0]
+    Hs = d' * d
+    ρ0 = [1 0; 0 0]
+    bB = Boson_DrudeLorentz_Pade(d' * d, 1, 1., 0.05, 3)
+    fB = Fermion_Lorentz_Pade(d, 1., 0, 1., 0.05, 3)
 
     PrecompileTools.@compile_workload begin
         # all calls in this block will be precompiled, regardless of whether
@@ -20,22 +22,23 @@ PrecompileTools.@setup_workload begin
 
         # precompile HEOM matrices
         @info "Precompiling HEOM Liouvillian superoperator matrices..."
-        Ms   = M_S(op; verbose=false)
-        Mb   = M_Boson(op, 2, bB; verbose=false, threshold=1e-1)
-        Mfo  = M_Fermion(op, 2, fB, ODD; verbose=false, threshold=1e-1)
-        Mbfe = M_Boson_Fermion(op, 2, 2, bB, fB; verbose=false, threshold=1e-1)
+        Ms   = M_S(Hs; verbose=false)
+        Mb   = M_Boson(Hs, 2, bB; verbose=false, threshold=1e-6)
+        Mfe  = M_Fermion(Hs, 5, fB, EVEN; verbose=false, threshold=1e-8)
+        Mfo  = M_Fermion(Hs, 5, fB, ODD;  verbose=false, threshold=1e-8)
+        Mbfe = M_Boson_Fermion(Hs, 2, 2, bB, fB; verbose=false, threshold=1e-6)
 
         # precompile Steadystate
         @info "Precompiling steady state solver..."
-        ados1 = SteadyState(Mb; verbose=false)
-        ados1 = SteadyState(Mb, [1. 0.; 0. 0.]; verbose=false, reltol=1e-1, abstol=1e-3)
-        E1 = Expect(op, ados1)
+        ados1 = SteadyState(Mfe; verbose=false)
+        ados1 = SteadyState(Mfe, ρ0; verbose=false)
+        E1 = Expect(Hs, ados1)
 
         # precompile evolution
         @info "Precompiling time evolution solver..."
-        ados2  = evolution(Mb, [1. 0.; 0. 0.], 0:1:1; verbose=false)
-        ados2  = evolution(Mb, [1. 0.; 0. 0.], 1, 1; verbose=false)
-        E2 = Expect(op, ados2)
+        ados2  = evolution(Mb, ρ0, 0:1:1; verbose=false)
+        ados2  = evolution(Mb, ρ0, 1, 1; verbose=false)
+        E2 = Expect(Hs, ados2)
 
         # precompile ADOs for the support of Base functions 
         @info "Precompiling Auxiliary Density Operators (ADOs)..."
@@ -49,8 +52,8 @@ PrecompileTools.@setup_workload begin
 
         # precompile Spectrum functions
         @info "Precompiling solvers for calculating spectrum..."
-        psd = PowerSpectrum(Mb,    [1. 0.; 0. 0.], op, [1]; verbose=false)
-        dos = DensityOfStates(Mfo, [1. 0.; 0. 0.], op, [1]; verbose=false)
+        psd = PowerSpectrum(  Mfo, ados1, d, [1]; verbose=false)
+        dos = DensityOfStates(Mfo, ados1, d, [1]; verbose=false)
     end
     @info "HierarchicalEOM precompilation complete"
 end

test/stationary_state.jl

@@ -1,20 +1,25 @@
 # System Hamiltonian and initial state
-Hsys = 0.25 * [1 0; 0 -1] + 0.5 * [0 1; 1 0]
+d    = [0 1; 0 0]
+Hsys = d' * d
 ρ0   = [1 0; 0 0]
 
 # Bath properties:
-λ  = 0.1
-W  = 0.5
-kT = 0.5
-N  = 2
-Q = [1 0; 0 -1]  # System-bath coupling operator
-bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N)
+Γ  =  1.
+W  =  1
+kT =  0.025851991
+μL =  1.
+μR = -1.
+N  =  3
+baths = [
+    Fermion_Lorentz_Pade(d, Γ, μL, W, kT, N), 
+    Fermion_Lorentz_Pade(d, Γ, μR, W, kT, N)
+]
 
 # HEOM Liouvillian superoperator matrix
 tier = 5
-L = M_Boson(Hsys, tier, bath; verbose=false)
+L = M_Fermion(Hsys, tier, baths; verbose=false)
 
-ados = SteadyState(L, ρ0; verbose=false, reltol=1e-2, abstol=1e-4)
+ados = SteadyState(L, ρ0; verbose=false)
 ρs   = getRho(ados)
 O = [1 0.5; 0.5 1]
 @test Expect(O, ados) ≈ real(tr(O * ρs))
@@ -31,4 +36,4 @@ bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2)
 @test_throws ErrorException SteadyState(L, mat2)
 @test_throws ErrorException SteadyState(L, ADOs(zeros(8), 2))
 @test_throws ErrorException SteadyState(L, ADOs(ados.data, ados.N, ODD))
-@test_throws ErrorException SteadyState(L, HEOMSuperOp(Q, ODD, L) * ados)
+@test_throws ErrorException SteadyState(L, HEOMSuperOp(d, ODD, L) * ados)