Use FillArrays.jl

Project.toml

@@ -6,6 +6,7 @@ authors = ["Yi-Te Huang"]
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 FastExpm = "7868e603-8603-432e-a1a1-694bd70b01f2"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -30,14 +31,15 @@ HierarchicalEOM_CUDAExt = "CUDA"
 CUDA = "5"
 DiffEqCallbacks = "2 - 4"
 FastExpm = "1.1"
+FillArrays = "1"
 IncompleteLU = "0.2"
 LinearAlgebra = "1"
 LinearSolve = "2.4.2 - 3"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "1"
 ProgressMeter = "1.11.0"
-QuantumToolbox = "0.38.1"
+QuantumToolbox = "0.39"
 Reexport = "1"
 SciMLBase = "2.105"
 SciMLOperators = "1.4"

docs/make.jl

@@ -6,7 +6,8 @@ using Documenter, HierarchicalEOM
 
 DocMeta.setdocmeta!(HierarchicalEOM, :DocTestSetup, :(using HierarchicalEOM); recursive = true)
 
-const DRAFT = false # set `true` to disable cell evaluation
+# some options for makedocs
+const DRAFT = get(ENV, "DRAFT", false) == "true"  # DRAFT   = true  disables cell evaluation
 
 const MathEngine = MathJax3(
     Dict(

docs/src/bath_boson/bosonic_bath_intro.md

@@ -62,7 +62,7 @@ c_list = correlation_function(bath, tlist)
 ```
 Here, `c_list` is a list which contains the value of ``C(t)`` corresponds to the given time series `tlist`.
 
-## Exponent
+## Exponential terms of the Bosonic Bath
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i`` and ``\gamma_i``:
 ```julia
 e = bath[2] # the 2nd-term

docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md

@@ -47,7 +47,7 @@ cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
-## Exponent
+## Exponential terms of the Bosonic Bath (RWA)
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i^\nu`` and ``\gamma_i^\nu``:
 ```julia
 e = bath[2] # the 2nd-term

docs/src/bath_fermion/fermionic_bath_intro.md

@@ -50,7 +50,7 @@ cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
-## Exponent
+## Exponential terms of the Fermionic Bath
 `HierarchicalEOM.jl` also supports users to access the specific exponential term with brackets `[]`. This returns an [`Exponent`](@ref) object, which contains the corresponding value of ``\eta_i^\nu`` and ``\gamma_i^\nu``:
 ```julia
 e = bath[2] # the 2nd-term

docs/src/libraryAPI.md

@@ -72,12 +72,14 @@ M_Fermion
 M_Fermion(Hsys::QuantumObject, tier::Int, Bath::Vector{FermionBath}, parity::AbstractParity=EVEN; threshold::Real=0.0, verbose::Bool=true)
 M_Boson_Fermion
 M_Boson_Fermion(Hsys::QuantumObject, tier_b::Int, tier_f::Int, Bath_b::Vector{BosonBath}, Bath_f::Vector{FermionBath}, parity::AbstractParity=EVEN; threshold::Real=0.0, verbose::Bool=true)
-size(M::HEOMSuperOp)
-size(M::HEOMSuperOp, dim::Int)
-size(M::AbstractHEOMLSMatrix)
-size(M::AbstractHEOMLSMatrix, dim::Int)
-eltype(M::HEOMSuperOp)
-eltype(M::AbstractHEOMLSMatrix)
+Base.size(M::HEOMSuperOp)
+Base.size(M::HEOMSuperOp, dim::Int)
+Base.size(M::AbstractHEOMLSMatrix)
+Base.size(M::AbstractHEOMLSMatrix, dim::Int)
+Base.eltype(M::HEOMSuperOp)
+Base.eltype(M::AbstractHEOMLSMatrix)
+iscached
+isconstant
 Propagator
 addBosonDissipator
 addFermionDissipator
@@ -88,8 +90,8 @@ addTerminator
 ```@docs
 ADOs
 ADOs(V::AbstractVector, N::Int)
-length(A::ADOs)
-eltype(A::ADOs)
+Base.length(A::ADOs)
+Base.eltype(A::ADOs)
 getRho
 getADO
 QuantumToolbox.expect
@@ -106,6 +108,7 @@ getIndexEnsemble
 ## [Time Evolution](@id lib-Time-Evolution)
 There are two function definitions of `HEOMsolve`, which depend on different methods to solve the time evolution:
 ```@docs
+HEOMsolveProblem
 HEOMsolve
 HEOMsolve_map
 TimeEvolutionHEOMSol

docs/src/time_evolution.md

@@ -66,7 +66,7 @@ The first method is implemented by solving the ordinary differential equation (O
 See the docstring of this method:  
 
 ```@docs
-HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, alg::OrdinaryDiffEqAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), progress_bar::Union{Val,Bool} = Val(true), inplace::Union{Val,Bool} = Val(true), kwargs...) where {T_state<:Union{QuantumObject,ADOs}}
+HEOMsolve(M::AbstractHEOMLSMatrix, ρ0::T_state, tlist::AbstractVector; e_ops::Union{Nothing,AbstractVector} = nothing, alg::AbstractODEAlgorithm = DP5(), H_t::Union{Nothing,Function} = nothing, params::NamedTuple = NamedTuple(), progress_bar::Union{Val,Bool} = Val(true), inplace::Union{Val,Bool} = Val(true), kwargs...) where {T_state<:Union{QuantumObject,ADOs}}
 ```
 
 ```julia

src/HeomBase.jl

@@ -1,5 +1,10 @@
 export AbstractHEOMLSMatrix
 
+@doc raw"""
+    abstract type AbstractHEOMLSMatrix{T}
+
+Abstract type for all HEOM liouvillian superoperators
+"""
 abstract type AbstractHEOMLSMatrix{T} end
 
 function Base.getproperty(M::AbstractHEOMLSMatrix, key::Symbol)

src/HierarchicalEOM.jl

@@ -4,6 +4,9 @@ module HierarchicalEOM
 using LinearAlgebra
 using SparseArrays
 
+import Base.Threads: @threads, nthreads, Channel
+import Pkg
+
 # Re-export QuantumToolbox
 import Reexport: @reexport
 @reexport using QuantumToolbox
@@ -17,11 +20,13 @@ import QuantumToolbox:
     _spost,
     _sprepost,
     _liouvillian,
+    _sum_lindblad_dissipators,
     _gen_dimensions,
     _get_dims_string,
     dimensions_to_dims,
     _save_func,
     _merge_saveat,
+    _merge_tstops,
     _merge_kwargs_with_callback,
     _get_expvals,
     _se_me_map_prob_func,
@@ -47,19 +52,24 @@ import SciMLBase:
     EnsembleThreads,
     FullSpecialize,
     CallbackSet,
-    NullParameters
+    NullParameters,
+    AbstractODEAlgorithm
 import SciMLOperators:
-    AbstractSciMLOperator, MatrixOperator, ScaledOperator, AddedOperator, update_coefficients!, concretize
-import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
+    SciMLOperators,
+    AbstractSciMLOperator,
+    MatrixOperator,
+    ScaledOperator,
+    AddedOperator,
+    update_coefficients!,
+    concretize
 import OrdinaryDiffEqLowOrderRK: DP5
 import DiffEqCallbacks: FunctionCallingCallback, TerminateSteadyState
 import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, KrylovJL_GMRES
 
 # other dependencies (in alphabetical order)
-import Base.Threads: @threads, nthreads, Channel
 import FastExpm: fastExpm
+import FillArrays: Eye
 import IncompleteLU: ilu
-import Pkg
 import ProgressMeter: Progress, next!
 
 # Basic functions

src/bath/BosonBath.jl

@@ -240,8 +240,7 @@ function bosonReal(op::QuantumObject, η_real::Vector{Ti}, γ_real::Vector{Tj})
     if N_exp_term != length(γ_real)
         error("The length of \'η_real\' and \'γ_real\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    return bosonReal(_spre(_op.data, Id_cache) - _spost(_op.data, Id_cache), _op.dimensions, η_real, γ_real, N_exp_term)
+    return bosonReal(_spre(_op.data) - _spost(_op.data), _op.dimensions, η_real, γ_real, N_exp_term)
 end
 
 @doc raw"""
@@ -285,9 +284,8 @@ function bosonImag(op::QuantumObject, η_imag::Vector{Ti}, γ_imag::Vector{Tj})
     if N_exp_term != length(γ_imag)
         error("The length of \'η_imag\' and \'γ_imag\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    spreQ = _spre(_op.data, Id_cache)
-    spostQ = _spost(_op.data, Id_cache)
+    spreQ = _spre(_op.data)
+    spostQ = _spost(_op.data)
     return bosonImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_imag, γ_imag, N_exp_term)
 end
 
@@ -340,9 +338,8 @@ function bosonRealImag(
     if (N_exp_term != length(η_imag)) || (N_exp_term != length(γ))
         error("The length of \'η_real\', \'η_imag\' and \'γ\' should be the same.")
     end
-    Id_cache = I(size(_op, 1))
-    spreQ = _spre(_op.data, Id_cache)
-    spostQ = _spost(_op.data, Id_cache)
+    spreQ = _spre(_op.data)
+    spostQ = _spost(_op.data)
     return bosonRealImag(spreQ - spostQ, spreQ + spostQ, _op.dimensions, η_real, η_imag, γ, N_exp_term)
 end
 
@@ -437,11 +434,10 @@ function bosonAbsorb(
     if (N_exp_term != length(γ_absorb)) || (N_exp_term != length(η_emit))
         error("The length of \'η_absorb\', \'γ_absorb\' and \'η_emit\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return bosonAbsorb(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data) - _spost(adjoint(_op).data),
         _op.dimensions,
         η_absorb,
         γ_absorb,
@@ -501,11 +497,10 @@ function bosonEmit(
     if (N_exp_term != length(γ_emit)) || (N_exp_term != length(η_absorb))
         error("The length of \'η_emit\', \'γ_emit\' and \'η_absorb\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return bosonEmit(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache) - _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data) - _spost(adjoint(_op).data),
         _op.dimensions,
         η_emit,
         γ_emit,

src/bath/FermionBath.jl

@@ -136,12 +136,11 @@ function fermionAbsorb(
     if (N_exp_term != length(γ_absorb)) || (N_exp_term != length(η_emit))
         error("The length of \'η_absorb\', \'γ_absorb\' and \'η_emit\' should all be the same.")
     end
-    Id_cache = I(size(_op, 1))
     return fermionAbsorb(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache),
-        _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data),
+        _spost(adjoint(_op).data),
         _op.dimensions,
         η_absorb,
         γ_absorb,
@@ -204,12 +203,11 @@ function fermionEmit(
         error("The length of \'η_emit\', \'γ_emit\' and \'η_absorb\' should all be the same.")
     end
 
-    Id_cache = I(size(_op, 1))
     return fermionEmit(
-        _spre(_op.data, Id_cache),
-        _spost(_op.data, Id_cache),
-        _spre(adjoint(_op).data, Id_cache),
-        _spost(adjoint(_op).data, Id_cache),
+        _spre(_op.data),
+        _spost(_op.data),
+        _spre(adjoint(_op).data),
+        _spost(adjoint(_op).data),
         _op.dimensions,
         η_emit,
         γ_emit,

src/density_of_states.jl

@@ -53,10 +53,8 @@ Calculate density of states for the fermionic system in frequency domain.
 
     # Handle d_op
     _tr = _Tr(M)
-    Id_sys = I(prod(d_op.dimensions))
-    Id_HEOM = I(M.N)
-    d_normal = HEOMSuperOp(spre(d_op, Id_sys), ODD, M; Id_cache = Id_HEOM)
-    d_dagger = HEOMSuperOp(spre(d_op', Id_sys), ODD, M; Id_cache = Id_HEOM)
+    d_normal = HEOMSuperOp(spre(d_op), ODD, M)
+    d_dagger = HEOMSuperOp(spre(d_op'), ODD, M)
     b_m = _HandleVectorType(M, (d_normal * ados).data)
     b_p = _HandleVectorType(M, (d_dagger * ados).data)
     _tr_d_normal = _HandleTraceVectorType(M, adjoint(d_normal.data) * _tr) # another adjoint will be applied in dot function later
@@ -80,7 +78,7 @@ Calculate density of states for the fermionic system in frequency domain.
         QuantumToolbox.settings.ProgressMeterKWARGS...,
     )
     i = convert(ElType, 1im)
-    I_total = I(size(M, 1))
+    I_total = Eye(size(M, 1))
     Iω1 = i * ωList[1] * I_total
     cache_m = init(LinearProblem(M.data.A - Iω1, b_m), alg, kwargs...)
     cache_p = init(LinearProblem(M.data.A + Iω1, b_p), alg, kwargs...)

src/deprecated.jl

@@ -56,14 +56,9 @@ _HEOMSuperOp_deprecated_method_error() = error(
     "Specifying `mul_basis = \"L\", \"R\", or \"LR\"` to `HEOMSuperOp` is deprecated, try to construct system SuperOperator manually by using `spre`, `spost`, or `sprepost`.",
 )
 
-HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int, mul_basis::AbstractString; Id_cache = I(N)) =
+HEOMSuperOp(op, opParity::AbstractParity, dims, N::Int, mul_basis::AbstractString) =
     _HEOMSuperOp_deprecated_method_error()
-HEOMSuperOp(
-    op,
-    opParity::AbstractParity,
-    refHEOMLS::AbstractHEOMLSMatrix,
-    mul_basis::AbstractString;
-    Id_cache = I(refHEOMLS.N),
-) = HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N, mul_basis; Id_cache = Id_cache)
-HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString; Id_cache = I(refADOs.N)) =
-    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N, mul_basis; Id_cache = Id_cache)
+HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix, mul_basis::AbstractString) =
+    HEOMSuperOp(op, opParity, refHEOMLS.dimensions, refHEOMLS.N, mul_basis)
+HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString) =
+    HEOMSuperOp(op, opParity, refADOs.dimensions, refADOs.N, mul_basis)