Merge pull request #102 from ytdHuang/compat/QT

Make `HierarchicalEOM` more compatible with `QuantumToolbox`

Project.toml

@@ -1,7 +1,7 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 authors = ["Yi-Te Huang <[email protected]>"]
-version = "2.1.3"
+version = "2.2.0"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
@@ -12,7 +12,6 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -36,7 +35,6 @@ LinearSolve = "2.4.2 - 2"
 OrdinaryDiffEqCore = "1"
 OrdinaryDiffEqLowOrderRK = "1"
 Pkg = "<0.0.1, 1"
-ProgressMeter = "1.7"
 QuantumToolbox = "0.13"
 Reexport = "1"
 SciMLBase = "2"

docs/src/LS_solvers.md

@@ -1,6 +1,6 @@
 # [LinearSolve solvers](@id LS-solvers)
 
-In this page, we list several recommended solvers provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving SteadyState and spectrum in hierarchical equations of motion approach.  
+In this page, we list several recommended solvers provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving [`steadystate`](@ref) and spectrum in hierarchical equations of motion approach.  
 
 Remember to import `LinearSolve.jl`
 

docs/src/extensions/CUDA.md

@@ -65,11 +65,11 @@ M_odd_cpu  = M_Fermion(Hsys, tier, bath_list, ODD)
 M_odd_gpu  = cu(M_odd_cpu)
 
 # solve steady state with CPU
-ados_ss = SteadyState(M_even_cpu);
+ados_ss = steadystate(M_even_cpu);
 ```
 
 !!! note "Note"
-    This extension does not support for solving [`SteadyState`](@ref doc-Stationary-State) on GPU since it is not efficient and might get wrong solutions. If you really want to obtain the stationary state with GPU, you can repeatedly solve the [`evolution`](@ref doc-Time-Evolution) until you find it.
+    This extension does not support for solving [stationary state](@ref doc-Stationary-State) on GPU since it is not efficient and might get wrong solutions. If you really want to obtain the stationary state with GPU, you can repeatedly solve the [`evolution`](@ref doc-Time-Evolution) until you find it.
 
 ### Solving time evolution with CPU
 

docs/src/libraryAPI.md

@@ -108,10 +108,10 @@ There are six function definitions of `evolution`, which depend on different inp
 evolution
 ```
 
-## Steady State
-There are three function definitions of `SteadyState`, which depend on different input types and methods to solve the stationary state:
+## Stationary State
+There are three function definitions of `steadystate`, which depend on different input types and methods to solve the stationary state:
 ```@docs
-SteadyState
+steadystate
 ```
 
 ## Spectrum

docs/src/spectrum.md

@@ -81,7 +81,7 @@ M::AbstractHEOMLSMatrix
 
 # the input state can be in either type (but usually ADOs):
 ρ::QuantumObject # the reduced density operator
-ρ::ADOs # the ADOs solved from "evolution" or "SteadyState"
+ρ::ADOs # the ADOs solved from "evolution" or "steadystate"
 
 P::QuantumObject 
 Q::QuantumObject
@@ -136,7 +136,7 @@ M_even = M_Fermion(Hs, tier, bath)
 M_odd  = M_Fermion(Hs, tier, bath, ODD) 
 
 # the input state can be in either type of density operator matrix or ADOs (but usually ADOs):
-ados = SteadyState(M_even)
+ados = steadystate(M_even)
 
 # the (usually annihilation) operator "d" as shown above
 d::QuantumObject 

docs/src/stationary_state.md

@@ -1,7 +1,7 @@
 # [Stationary State](@id doc-Stationary-State)
 `HierarchicalEOM.jl` implements two different ways to calculate stationary states of all [Auxiliary Density Operators (ADOs)](@ref doc-ADOs).
 
-To solve the stationary state of the reduced state and also all the [ADOs](@ref doc-ADOs), you only need to call [`SteadyState`](@ref). Different methods are implemented with different input parameters of the function which makes it easy to switch between different methods. The output of the function [`SteadyState`](@ref) for each methods will always be in the type of the auxiliary density operators [`ADOs`](@ref).
+To solve the stationary state of the reduced state and also all the [ADOs](@ref doc-ADOs), you only need to call [`steadystate`](@ref). Different methods are implemented with different input parameters of the function which makes it easy to switch between different methods. The output of the function [`steadystate`](@ref) for each methods will always be in the type of the auxiliary density operators [`ADOs`](@ref).
 
 ## Solve with [LinearSolve.jl](http://linearsolve.sciml.ai/stable/)
 The first method is implemented by solving the linear problem
@@ -14,19 +14,19 @@ The first method is implemented by solving the linear problem
 See the docstring of this method:  
 
 ```@docs
-SteadyState(M::AbstractHEOMLSMatrix; solver=UMFPACKFactorization(), verbose::Bool=true, SOLVEROptions...)
+steadystate(M::AbstractHEOMLSMatrix; solver=UMFPACKFactorization(), verbose::Bool=true, SOLVEROptions...)
 ```
 
 ```julia
 # the HEOMLS matrix
 M::AbstractHEOMLSMatrix  
-ados_steady = SteadyState(M)
+ados_steady = steadystate(M)
 ```
 !!! warning "Unphysical solution"
     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) method : [`SteadyStateDiffEq.jl`](https://github.com/SciML/SteadyStateDiffEq.jl)
+The second method is implemented by solving the ordinary differential equation (ODE) method :
 ```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::QuantumObject, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
+steadystate(M::AbstractHEOMLSMatrix, ρ0::QuantumObject, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
 ```
 
 ```julia
@@ -49,13 +49,13 @@ M::AbstractHEOMLSMatrix
 # the initial state of the system density operator
 ρ0::QuantumObject
 
-ados_steady = SteadyState(M, ρ0)
+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, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
+steadystate(M::AbstractHEOMLSMatrix, ados::ADOs, tspan::Number = Inf; solver = DP5(), termination_condition = NormTerminationMode(), verbose::Bool = true, SOLVEROptions...)
 ```
 
 ```julia
@@ -65,5 +65,5 @@ M::AbstractHEOMLSMatrix
 # the initial state of the ADOs
 ados::ADOs
 
-ados_steady = SteadyState(M, ados)
+ados_steady = steadystate(M, ados)
 ```

examples/SIAM.jl

@@ -54,7 +54,7 @@ M_odd = M_Fermion(Hsys, tier, bath_list, ODD)
 
 # ## Solve stationary state of ADOs
 # (see also [Stationary State](@ref doc-Stationary-State))
-ados_s = SteadyState(M_even)
+ados_s = steadystate(M_even)
 
 # ## Calculate density of states (DOS)
 # (see also [Spectrum](@ref doc-Spectrum))

examples/benchmark_LS_solvers.jl

@@ -1,6 +1,6 @@
 # # [LinearSolve solvers](@id benchmark-LS-solvers)
 # 
-# In this page, we will benchmark several solvers provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving SteadyState and spectrum in hierarchical equations of motion approach.  
+# In this page, we will benchmark several solvers provided by [LinearSolve.jl](https://docs.sciml.ai/LinearSolve/stable/) for solving steadystate and spectrum in hierarchical equations of motion approach.  
 
 using LinearSolve
 using BenchmarkTools
@@ -30,7 +30,7 @@ bath_dn = Fermion_Lorentz_Pade(d_dn, Γ, μ, W, kT, N)
 bath_list = [bath_up, bath_dn]
 M_even = M_Fermion(Hsys, tier, bath_list)
 M_odd = M_Fermion(Hsys, tier, bath_list, ODD)
-ados_s = SteadyState(M_even);
+ados_s = steadystate(M_even);
 
 # ## LinearSolve Solver List
 # (click [here](https://docs.sciml.ai/LinearSolve/stable/solvers/solvers/) to see the full solver list provided by `LinearSolve.jl`)
@@ -57,19 +57,19 @@ MKLPardisoIterate();
 # ## Solving Stationary State
 # Since we are using [`BenchmarkTools`](https://juliaci.github.io/BenchmarkTools.jl/stable/) (`@benchmark`) in the following, we set `verbose = false` to disable the output message.
 # ### UMFPACKFactorization (Default solver)
-@benchmark SteadyState(M_even; verbose = false)
+@benchmark steadystate(M_even; verbose = false)
 
 # ### KLUFactorization
-@benchmark SteadyState(M_even; solver = KLUFactorization(), verbose = false)
+@benchmark steadystate(M_even; solver = KLUFactorization(), verbose = false)
 
 # ### KrylovJL_BICGSTAB
-@benchmark SteadyState(M_even; solver = KrylovJL_BICGSTAB(rtol = 1e-10, atol = 1e-12), verbose = false)
+@benchmark steadystate(M_even; solver = KrylovJL_BICGSTAB(rtol = 1e-10, atol = 1e-12), verbose = false)
 
 # ### MKLPardisoFactorize
-@benchmark SteadyState(M_even; solver = MKLPardisoFactorize(), verbose = false)
+@benchmark steadystate(M_even; solver = MKLPardisoFactorize(), verbose = false)
 
 # ### MKLPardisoIterate
-@benchmark SteadyState(M_even; solver = MKLPardisoIterate(), verbose = false)
+@benchmark steadystate(M_even; solver = MKLPardisoIterate(), verbose = false)
 
 # ## Calculate Spectrum
 # ### UMFPACKFactorization (Default solver)

examples/cavityQED.jl

@@ -105,7 +105,7 @@ evo_H = evolution(M_Heom, ψ0, t_list);
 
 # ## Solve stationary state of ADOs
 # (see also [Stationary State](@ref doc-Stationary-State))
-steady_H = SteadyState(M_Heom);
+steady_H = steadystate(M_Heom);
 
 # ## Expectation values
 # observable of atom: $\sigma_z$
@@ -164,7 +164,7 @@ M_master = addBosonDissipator(M_master, jump_op)
 evo_M = evolution(M_master, ψ0, t_list);
 
 ## steady
-steady_M = SteadyState(M_master);
+steady_M = steadystate(M_master);
 
 ## expectation value of σz
 σz_evo_M = expect(σz, evo_M)

examples/electronic_current.jl

@@ -58,7 +58,7 @@ ados_evolution = evolution(M, ψ0, tlist);
 
 # ## Solve stationary state of ADOs
 # (see also [Stationary State](@ref doc-Stationary-State))
-ados_steady = SteadyState(M);
+ados_steady = steadystate(M);
 
 # ## Calculate current
 # Within the influence functional approach, the expectation value of the electronic current from the $\alpha$-fermionic bath into the system can be written in terms of the first-level-fermionic ($n=1$) auxiliary density operators, namely

examples/quick_start.jl

@@ -13,10 +13,10 @@
 # Here are the functions in `HierarchicalEOM.jl` that we will use in this tutorial (Quick Start):
 
 import HierarchicalEOM
-import HierarchicalEOM: Boson_DrudeLorentz_Pade, M_Boson, evolution, SteadyState, getRho, BosonBath
+import HierarchicalEOM: Boson_DrudeLorentz_Pade, M_Boson, evolution, getRho, BosonBath
 
 # Note that you can also type `using HierarchicalEOM` to import everything you need in `HierarchicalEOM.jl`.
-# To check the versions of dependencies of `HierarchicalEOM.jl` , run the following function
+# To check the versions of dependencies of `HierarchicalEOM.jl`, run the following function
 
 HierarchicalEOM.versioninfo()
 
@@ -33,7 +33,7 @@ HierarchicalEOM.versioninfo()
 # #### System Hamiltonian and initial state
 # You must construct system hamiltonian, initial state, and coupling operators by [`QuantumToolbox`](https://github.com/qutip/QuantumToolbox.jl) framework. It provides many useful functions to create arbitrary quantum states and operators which can be combined in all the expected ways.
 
-import QuantumToolbox: Qobj, sigmaz, sigmax, basis, ket2dm, expect
+import QuantumToolbox: Qobj, sigmaz, sigmax, basis, ket2dm, expect, steadystate
 
 ## The system Hamiltonian
 ϵ = 0.5 # energy of 2-level system
@@ -76,10 +76,10 @@ ados_list = evolution(L, ρ0, tlist);
 # To learn more about `evolution`, please refer to [Time Evolution](@ref doc-Time-Evolution).
 
 # ### Stationary State
-# We can also solve the stationary state of the auxiliary density operators (ADOs) by calling [`SteadyState`](@ref).
-ados_steady = SteadyState(L)
+# We can also solve the stationary state of the auxiliary density operators (ADOs) by calling [`steadystate`](@ref).
+ados_steady = steadystate(L)
 
-# To learn more about `SteadyState`, please refer to [Stationary State](@ref doc-Stationary-State).
+# To learn more about `steadystate`, please refer to [Stationary State](@ref doc-Stationary-State).
 
 # ### Reduced Density Operator
 # To obtain the reduced density operator, one can either access the first element of auxiliary density operator (`ADOs`) or call [`getRho`](@ref):
@@ -96,7 +96,7 @@ ados_steady = SteadyState(L)
 # state but also easily take high-order terms into account without struggling in finding the indices (see [Auxiliary Density Operators](@ref doc-ADOs) and [Hierarchy Dictionary](@ref doc-Hierarchy-Dictionary) for more details).
 
 # ### Expectation Value
-# We can now compare the results obtained from `evolution` and `SteadyState`:
+# We can now compare the results obtained from `evolution` and `steadystate`:
 
 ## Define the operators that measure the populations of the two
 ## system states:

src/HeomBase.jl

@@ -1,7 +1,5 @@
 abstract type AbstractHEOMLSMatrix end
 
-const PROGBAR_OPTIONS = Dict(:barlen => 20, :color => :green, :showspeed => true)
-
 # equal to : transpose(sparse(vec(system_identity_matrix)))
 function _Tr(dims::SVector, N::Int)
     D = prod(dims)

src/HierarchicalEOM.jl

@@ -12,7 +12,6 @@ module HeomBase
     import QuantumToolbox: QuantumObject, QuantumObjectType, Operator, SuperOperator
 
     export _Tr,
-        PROGBAR_OPTIONS,
         AbstractHEOMLSMatrix,
         _check_sys_dim_and_ADOs_num,
         _check_parity,
@@ -100,12 +99,12 @@ module HeomAPI
         sprepost,
         expect,
         ket2dm,
-        lindblad_dissipator
-    import ProgressMeter: Progress, next!
+        lindblad_dissipator,
+        steadystate,
+        ProgressBar,
+        next!
     import FastExpm: fastExpm
-
-    # solving time evolution and steady state
-    import SciMLBase: solve, solve!, init, step!, ODEProblem
+    import SciMLBase: init, solve, solve!, step!, ODEProblem
     import SciMLOperators: MatrixOperator
     import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm
     import OrdinaryDiffEqLowOrderRK: DP5
@@ -137,7 +136,9 @@ module HeomAPI
         addFermionDissipator,
         addTerminator,
         evolution,
-        SteadyState
+        SteadyState, # has been deprecated, throws error only
+        PowerSpectrum,
+        DensityOfStates
 
     include("Parity.jl")
     include("ADOs.jl")
@@ -151,24 +152,10 @@ module HeomAPI
     include("heom_matrices/M_Boson_Fermion.jl")
 
     include("evolution.jl")
-    include("SteadyState.jl")
-end
-@reexport using .HeomAPI
-
-# sub-module Spectrum for HierarchicalEOM
-module Spectrum
-    using ..HeomBase
-    import ..HeomAPI: HEOMSuperOp, ADOs, EVEN, ODD
-    import SciMLBase: init, solve!
-    import LinearSolve: LinearProblem, SciMLLinearSolveAlgorithm, UMFPACKFactorization
-    import ProgressMeter: Progress, next!
-    import QuantumToolbox: QuantumObject
-
-    export PowerSpectrum, DensityOfStates
-
+    include("steadystate.jl")
     include("power_spectrum.jl")
     include("density_of_states.jl")
 end
-@reexport using .Spectrum
+@reexport using .HeomAPI
 
 end

src/density_of_states.jl

@@ -75,7 +75,7 @@ Calculate density of states for the fermionic system in frequency domain.
     if verbose
         print("Calculating density of states in frequency domain...\n")
         flush(stdout)
-        prog = Progress(Length; desc = "Progress : ", PROGBAR_OPTIONS...)
+        prog = ProgressBar(Length)
     end
     i = convert(ElType, 1im)
     I_total = _HandleIdentityType(typeof(M.data), Size)

src/evolution.jl

@@ -108,7 +108,7 @@ For more details, please refer to [`FastExpm.jl`](https://github.com/fmentink/Fa
     if verbose
         print("Solving time evolution for ADOs by propagator method...\n")
         flush(stdout)
-        prog = Progress(steps + 1; start = 1, desc = "Progress : ", PROGBAR_OPTIONS...)
+        prog = ProgressBar(steps + 1)
     end
     for n in 1:steps
         ρvec = exp_Mt * ρvec
@@ -264,7 +264,7 @@ with initial state is given in the type of `ADOs`.
     if verbose
         print("Solving time evolution for ADOs by Ordinary Differential Equations method...\n")
         flush(stdout)
-        prog = Progress(length(Tlist); start = 1, desc = "Progress : ", PROGBAR_OPTIONS...)
+        prog = ProgressBar(length(Tlist))
     end
     idx = 1
     dt_list = diff(Tlist)
@@ -426,7 +426,7 @@ with initial state is given in the type of `ADOs`.
             "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...)
+        prog = ProgressBar(length(Tlist))
     end
 
     parameters = (H = H, param = param, dims = M.dims, N = M.N, parity = M.parity, Id_cache = Id)
@@ -452,7 +452,7 @@ with initial state is given in the type of `ADOs`.
             "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...)
+        prog = ProgressBar(length(Tlist))
     end
     idx = 1
     dt_list = diff(Tlist)

src/heom_matrices/M_Boson.jl

@@ -89,7 +89,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
     if verbose
         println("Preparing block matrices for HEOM Liouvillian superoperator (using $(Nthread) threads)...")
         flush(stdout)
-        prog = Progress(Nado; desc = "Processing: ", PROGBAR_OPTIONS...)
+        prog = ProgressBar(Nado)
     end
     @threads for idx in 1:Nado
         tID = threadid()