Merge pull request #43 Generalize spectrum calculation

Generalize spectrum calculation and optimize runtests

.buildkite/pipeline.yml

@@ -15,6 +15,6 @@ steps:
     env:
       GROUP: "HierarchicalEOM_CUDAExt"
       SECRET_CODECOV_TOKEN: "xPPS+qHx6R93hyaYefGl5Hj0ojomvcMvB0azWyK0LkOrhuHNT0OegnOvcDfckKkPAw3c3zDPZyO7ogBUSUw/R9JRXF8oQd3yOxySp1vcdjsNsMvazW8Nc9QS+Cbb7/yQmtRHtkPNc4pzI0HsGPqr6pKyX8GOZRdspQ6R6KCmP95AF2KPgvhW+t/I98Z0qz8ksKgE6OQvo0XMSLKlAQRthWIBmeYf1IgjSxpibhLZhb5OKYjOUElGvRWL79SFBmbVqIHAdAVOMaL8VVoRMfgqE89SvshYEiDLNBCenxNrLAYcOHEA4ElrhmUhlVDwsqYfq7ooxHRIEkaIZ5/sL348LA==;U2FsdGVkX1+pTkirGXYb1wl7BDHCxbkVqFTOLY0ciN8+3TKlwnrtGpP+PrBs+fqJTIch2eJndJfpvbAL2DAzyg=="
-    timeout_in_minutes: 30
+    timeout_in_minutes: 60
     # Don't run Buildkite if the commit message includes the text [skip tests]
     if: build.message !~ /\[skip tests\]/

Project.toml

@@ -40,15 +40,7 @@ SteadyStateDiffEq = "1.16.0"
 julia = "1.8"
 
 [extras]
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "Documenter", "LaTeXStrings", "Literate", "Pardiso", "Pkg", "Plots", "QuantumOptics", "Test"]
+test = ["Test"]

docs/src/extensions/CUDA.md

@@ -1,11 +1,11 @@
 # [Extension for CUDA.jl](@id doc-ext-CUDA)
 
-This is an extension to support GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) acceleration for solving the [time evolution](@ref doc-Time-Evolution) and [spectrum](@ref doc-Spectrum). This improves the execution time and memory usage especially when the HEOMLS matrix is super large.
+This is an extension to support GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) acceleration for solving the [time evolution](@ref doc-Time-Evolution) and [spectra](@ref doc-Spectrum). This improves the execution time and memory usage especially when the HEOMLS matrix is super large.
 
 !!! compat "Compat"
     The described feature requires `Julia 1.9+`.
 
-The functions [`evolution`](@ref doc-Time-Evolution) (only supports ODE method with time-independent system Hamiltonian) and [`spectrum`](@ref doc-Spectrum) will automatically choose to solve on CPU or GPU depend on the type of the sparse matrix in `M::AbstractHEOMLSMatrix` objects (i.e., the type of the field `M.data`). 
+The functions of calculating [time evolution](@ref doc-Time-Evolution) (only supports ODE method with time-independent system Hamiltonian) and [spectra](@ref doc-Spectrum) will automatically choose to solve on CPU or GPU depend on the type of the sparse matrix in `M::AbstractHEOMLSMatrix` objects (i.e., the type of the field `M.data`). 
 
 ```julia
 typeof(M.data) <:   SparseMatrixCSC # solve on CPU
@@ -85,11 +85,11 @@ ados_list_gpu = evolution(M_even_gpu, ρ0, tlist)
 ### Solving Spectrum with CPU
 
 ```julia
-dos_cpu = spectrum(M_odd_cpu, ados_ss, d_up, ωlist)
+dos_cpu = DensityOfStates(M_odd_cpu, ados_ss, d_up, ωlist)
 ```
 
 ### Solving Spectrum with GPU
 
 ```julia
-dos_gpu = spectrum(M_odd_gpu, ados_ss, d_up, ωlist; solver=KrylovJL_BICGSTAB(rtol=1f-10, atol=1f-12))
+dos_gpu = DensityOfStates(M_odd_gpu, ados_ss, d_up, ωlist; solver=KrylovJL_BICGSTAB(rtol=1f-10, atol=1f-12))
 ```

docs/src/heom_matrix/HEOMLS_intro.md

@@ -50,7 +50,7 @@ M = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath; threshold=1e-7)
     The full hierarchical equations can be recovered in the limiting case ``\mathcal{I}_\textrm{th}\rightarrow 0``, which is the default value of the parameter : `threshold=0.0`. This means that all of the ADOs will be taken into account by default.
 
 ## [Parity Support for HEOMLS Matrices](@id doc-Parity)
-When the system Hamiltonian contains fermionic systems, the HEOMLS matrix ``\hat{\mathcal{M}}`` might be constructed into a different one depend on the parity of the input operator (HEOMLS) it is acting on. Usually, it is acting on the reduced density operator and [auxiliary density operators (ADOs)](@ref doc-ADOs), which are all in `EVEN`-parity. However, there are some situations (for example, [calculating spectrum for fermionic systems](@ref doc-DOS)) where ``\hat{\mathcal{M}}`` is acting on ADOs with `ODD`-parity.
+When the system Hamiltonian contains fermionic systems, the HEOMLS matrix ``\hat{\mathcal{M}}`` might be constructed into a different one depend on the parity of the input operator (HEOMLS) it is acting on. Usually, it is acting on the reduced density operator and [auxiliary density operators (ADOs)](@ref doc-ADOs), which are all in `EVEN`-parity. However, there are some situations (for example, [calculating density of states for fermionic systems](@ref doc-DOS)) where ``\hat{\mathcal{M}}`` is acting on ADOs with `ODD`-parity.
 
 One can specify the parameter `parity::AbstractParity` in the function of constructing ``\hat{\mathcal{M}}`` to be [`EVEN`](@ref) or [`ODD`](@ref). The default value of the parameter is `parity=EVEN`.
 ```julia

docs/src/install.md

@@ -32,10 +32,10 @@ In order to get a better experience and take full advantage of `HierarchicalEOM`
     `HierarchicalEOM.jl` provides an extension to support `QuantumOptics`-type object, but this feature requires `Julia 1.9+` and `HierarchicalEOM 0.3+`. See [here](@ref doc-ext-QuantumOptics) for more details.
 
 ### [DifferentialEquations.jl](https://diffeq.sciml.ai/stable/)
-`DifferentialEquations` is needed to provide the low-level ODE solvers especially for solving time [`evolution`](@ref). For [low dependency usage](https://diffeq.sciml.ai/stable/features/low_dep/), users can use [`OrdinaryDiffEq.jl`](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) instead.
+`DifferentialEquations` is needed to provide the low-level ODE solvers especially for solving [time evolution](@ref doc-Time-Evolution). For [low dependency usage](https://diffeq.sciml.ai/stable/features/low_dep/), users can use [`OrdinaryDiffEq.jl`](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) instead.
 
 ### [LinearSolve.jl](http://linearsolve.sciml.ai/stable/)
-`LinearSolve` is a unified interface for the linear solving packages of Julia. It interfaces with other packages of the Julia ecosystem to make it easier to test alternative solver packages and pass small types to control algorithm swapping. It is needed to provide the solvers especially for solving [`SteadyState`](@ref) and [`spectrum`](@ref) for both bosonic and fermionic systems.
+`LinearSolve` is a unified interface for the linear solving packages of Julia. It interfaces with other packages of the Julia ecosystem to make it easier to test alternative solver packages and pass small types to control algorithm swapping. It is needed to provide the solvers especially for solving [stationary state](@ref doc-Stationary-State) and [spectra](@ref doc-Spectrum) for both bosonic and fermionic systems.
 
 ### [JLD2.jl](https://juliaio.github.io/JLD2.jl/stable/)
 `JLD2` saves and loads Julia data structures in a format comprising a subset of HDF5. Because the size of matrix in `HierarchicalEOM` is usually super large and leads to long time calculation, we support the functionality for saving and loading the `HierarchicalEOM`-type objects into files by `JLD2 >= 0.4.23`.

docs/src/libraryAPI.md

@@ -40,7 +40,7 @@ fermionEmit
 fermionEmit(op::AbstractMatrix, η_emit::Vector{Ti}, γ_emit::Vector{Tj}, η_absorb::Vector{Tk}) where {Ti, Tj, Tk <: Number}
 ```
 
-## Correlation Functions
+## Bath Correlation Functions
 
 ```@docs
 Boson_DrudeLorentz_Matsubara
@@ -110,6 +110,8 @@ SteadyState
 ## Spectrum
 ```@docs
 spectrum
+PowerSpectrum
+DensityOfStates
 ```
 
 ## Misc.

docs/src/spectrum.md

@@ -3,64 +3,98 @@
 ## Introduction
 We briefly summarize how to numerically compute the spectrum associated with the system degree of freedom. [Phys. Rev. Lett. 109, 266403 (2012)](https://link.aps.org/doi/10.1103/PhysRevLett.109.266403) showed that the spectrum can be evaluated either in time or frequency domain.
 
-`HierarchicalEOM.jl` provides a function [`spectrum`](@ref) which performs the calculation in frequency domain. There are two different methods (as shown below) which depends on the [parity](@ref doc-Parity) of the HEOMLS matrices ``\hat{\mathcal{M}}`` corresponds to different system degree of freedom. 
+`HierarchicalEOM.jl` provides the following listed functions which performs the calculation of spectrum in frequency domain.
 
-If you want to calculate the spectrum associated with
- - [bosonic systems (Power Spectral Density)](@ref doc-PSD) : you have to provide ``\hat{\mathcal{M}}`` constructed in `EVEN` parity.
- - [fermionic systems (Density of States)](@ref doc-DOS) : you have to provide ``\hat{\mathcal{M}}`` constructed in `ODD` parity.
-
-The function [`spectrum`](@ref) will automatically detect the [parity](@ref doc-Parity) of ``\hat{\mathcal{M}}`` by itself. Furthermore, the output of the function [`spectrum`](@ref) for both cases will always be in the type of `Vector{Float64}`, which contains the list of the spectrum values corresponding to the given `ω_list`.
+ - [Power Spectrum](@ref doc-PS)
+ - [Density of States](@ref doc-DOS)
 
 `HierarchicalEOM.jl` wraps some of the functions in [LinearSolve.jl](http://linearsolve.sciml.ai/stable/), which is a very rich numerical library for solving the linear problems and provides many solvers. It offers quite a few options for the user to tailor the solver to their specific needs. The default solver (and its corresponding settings) are chosen to suit commonly encountered problems and should work fine for most of the cases. If you require more specialized methods, such as the choice of algorithm, please refer to [benchmark for LinearSolve solvers](@ref benchmark-LS-solvers) and also the documentation of [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
 
 !!! compat "Extension for CUDA.jl"
     `HierarchicalEOM.jl` provides an extension to support GPU ([`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl)) acceleration for solving the spectrum, but this feature requires `Julia 1.9+` and `HierarchicalEOM 1.1+`. See [here](@ref doc-ext-CUDA) for more details.
 
-## [Power Spectral Density](@id doc-PSD)
-Start from the spectrum for bosonic systems (power spectral density) in the time-domain. We write the system two-time correlation function in terms of the propagator ``\hat{\mathcal{G}}(t)=\exp(\hat{\mathcal{M}} t)`` for ``t>0``. The power spectral density ``S(\omega)`` can be obtained as
+The output of the above listed functions will always be in the type of `Vector{Float64}`, which contains the list of the spectrum values corresponding to the given `ωlist`.
+
+## Common and optional parameters
+Furthermore, there are two common optional parameters for all the functions provided below:
+- `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`.
+- `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process.
+
+If the filename is specified, the function will automatically save (update) the value (together with a comma behind it) to a new line in the file (with ".txt" behind the filename) once it obtains the solution of each specified ``\omega``. For example, if you specify `filename="test"` and `ωlist=0:1:5`, you will obtain a file `test.txt` where each line in this file (as shown below) is the result of spectrum corresponding to the given `ωlist`:
+```
+# (the content inside test.txt) #
+0.4242990296334028,
+0.28617768129333854,
+0.21332961856387556,
+0.1751179183484055,
+0.15739257286986685,
+0.1518018484057393,
+
+```
+
+For your convenience, we add those commas (",") in the end of each line for the users to easily do "copy-and-paste" and load these results back into julia's kernel (construct a vector of the given results) again, namely
+
+```julia
+results = [
+0.4242990296334028,
+0.28617768129333854,
+0.21332961856387556,
+0.1751179183484055,
+0.15739257286986685,
+0.1518018484057393
+]
+```
+
+## [Power Spectrum](@id doc-PS)
+Start from the power spectrum in the time-domain. We write the system two-time correlation function in terms of the propagator ``\hat{\mathcal{G}}(t)=\exp(\hat{\mathcal{M}} t)`` for ``t>0``. The power spectrum ``\pi S(\omega)`` can be obtained as
 ```math
 \begin{aligned}
 \pi S(\omega) 
-&= \textrm{Re}\left\{\int_0^\infty dt \langle a^\dagger(t)a(0)\rangle e^{-i\omega t}\right\}\\
-&= \textrm{Re}\left\{\int_0^\infty dt \langle a^\dagger e^{\hat{\mathcal{M}} t}a\rangle e^{-i\omega t}\right\}\\
-&= -\textrm{Re}\left\{\langle a^\dagger (\hat{\mathcal{M}} -i\omega)^{-1} a\rangle\right\}\\
-&= -\textrm{Re}\left\{\textrm{Tr}\left[ a^\dagger (\hat{\mathcal{M}} -i\omega)^{-1} a\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}\right]\right\},
+&= \textrm{Re}\left\{\int_0^\infty dt \langle P(t)Q(0)\rangle e^{-i\omega t}\right\}\\
+&= \textrm{Re}\left\{\int_0^\infty dt \langle P e^{\hat{\mathcal{M}} t}Q\rangle e^{-i\omega t}\right\}\\
+&= -\textrm{Re}\left\{\langle P (\hat{\mathcal{M}} -i\omega)^{-1} Q\rangle\right\}\\
+&= -\textrm{Re}\left\{\textrm{Tr}\left[ P (\hat{\mathcal{M}} -i\omega)^{-1} Q\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}\right]\right\},
 \end{aligned}
 ```
 where a half-Fourier transform has been introduced in the third line. We note that only the reduced density operator (``m=n=0``) is considered when taking the final trace operation.
 
-The function [`spectrum`](@ref) solves the linear problem ``\textbf{A x}=\textbf{b}`` at a fixed frequency ``\omega`` where 
+This function solves the linear problem ``\textbf{A x}=\textbf{b}`` at a fixed frequency ``\omega`` where 
  - ``\textbf{A}=\hat{\mathcal{M}}-i\omega``
- - ``\textbf{b}=a\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` 
+ - ``\textbf{b}=Q\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` 
 using the package [LinearSolve.jl](http://linearsolve.sciml.ai/stable/).
 
-Finially, one can obtain the value of the power spectral density for specific ``\omega``, namely
+Finially, one can obtain the value of the power spectrum for specific ``\omega``, namely
 ```math
-\pi S(\omega) = -\textrm{Re}\left\{\textrm{Tr}\left[ a^\dagger \textbf{x}\right]\right\}.
+\pi S(\omega) = -\textrm{Re}\left\{\textrm{Tr}\left[ P \textbf{x}\right]\right\}.
 ```
 
-See also the docstring : [`spectrum`](@ref)
+!!! note "Odd-Parity for Power Spectrum"
+    When ``Q`` is an operator acting on fermionic systems and has `ODD`-parity, the HEOMLS matrix ``\hat{\mathcal{M}}`` is acting on the `ODD`-parity space because ``\textbf{b}=Q\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}``. Therefore, remember to construct ``\hat{\mathcal{M}}`` with `ODD` [parity](@ref doc-Parity) in this kind of cases.
+
+See also the docstring : 
+```@docs
+PowerSpectrum(M::AbstractHEOMLSMatrix, ρ, P_op, Q_op, ωlist::AbstractVector, reverse::Bool = false; solver = UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
+```
 
 ```julia
-M::AbstractHEOMLSMatrix # need to be in "EVEN" parity
+M::AbstractHEOMLSMatrix
 
 # the input state can be in either type (but usually ADOs):
 ρ::AbstractMatrix # the reduced density operator
 ρ::ADOs # the ADOs solved from "evolution" or "SteadyState"
 
-# the (usually annihilation) operator "a" as shown above
-a::AbstractMatrix 
+P::AbstractMatrix 
+Q::AbstractMatrix
 
 # the spectrum value for the specific frequency ω which need to be solved
-ω_list = 0:0.5:2 # [0.0, 0.5, 1.0, 1.5, 2.0]
+ωlist = 0:0.5:2 # [0.0, 0.5, 1.0, 1.5, 2.0]
 
-πSω = spectrum(M, ρ, a, ω_list)
+πSω = spectrum(M, ρ,    Q, ωlist) # P will automatically be considered as "the adjoint of Q operator"
+πSω = spectrum(M, ρ, P, Q, ωlist) # user specify both P and Q operator
 ```
-!!! note "Note"
-    To calculate power spectral density, remember to construct ``\hat{\mathcal{M}}`` with `EVEN` [parity](@ref doc-Parity).
 
 ## [Density of States](@id doc-DOS)
-Start from the spectrum for fermionic systems (density of states) in the time-domain. We write the system two-time correlation function in terms of the propagator ``\hat{\mathcal{G}}(t)=\exp(\hat{\mathcal{M}} t)`` for ``t>0``. The density of states ``A(\omega)`` can be obtained as
+Start from the density of states for fermionic systems in the time-domain. We write the system two-time correlation function in terms of the propagator ``\hat{\mathcal{G}}(t)=\exp(\hat{\mathcal{M}} t)`` for ``t>0``. The density of states ``\pi A(\omega)`` can be obtained as
 ```math
 \begin{aligned}
 \pi A(\omega) 
@@ -72,7 +106,7 @@ Start from the spectrum for fermionic systems (density of states) in the time-do
 ```
 where a half-Fourier transform has been introduced in the third line. We note that only the reduced density operator (``m=n=0``) is considered when taking the final trace operation.
 
-The function [`spectrum`](@ref) solves two linear problems ``\textbf{A}_+ \textbf{x}_+=\textbf{b}_+`` and ``\textbf{A}_- \textbf{x}_-=\textbf{b}_-`` at a fixed frequency ``\omega`` where 
+This functionsolves two linear problems ``\textbf{A}_+ \textbf{x}_+=\textbf{b}_+`` and ``\textbf{A}_- \textbf{x}_-=\textbf{b}_-`` at a fixed frequency ``\omega`` where 
  - ``\textbf{A}_+=\hat{\mathcal{M}}+i\omega``
  - ``\textbf{b}_+=d^\dagger\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` 
  - ``\textbf{A}_-=\hat{\mathcal{M}}-i\omega``
@@ -85,10 +119,12 @@ Finially, one can obtain the density of states for specific ``\omega``, namely
 ```
 
 !!! note "Odd-Parity for Density of States"
-    As shown above, the HEOMLS matrix ``\hat{\mathcal{M}}`` acts on the `ODD`-parity space, compatibly with the parity of both the operators ``d\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` and ``d^\dagger\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}``.  
-    Therefore, remember to construct ``\hat{\mathcal{M}}`` with `ODD` [parity](@ref doc-Parity) for solving spectrum of fermionic systems.
+    As shown above, the HEOMLS matrix ``\hat{\mathcal{M}}`` acts on the `ODD`-parity space, compatibly with the parity of both the operators ``\textbf{b}_-=d\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}`` and ``\textbf{b}_+=d^\dagger\rho^{(m,n,+)}_{\textbf{j} \vert \textbf{q}}``. Therefore, remember to construct ``\hat{\mathcal{M}}`` with `ODD` [parity](@ref doc-Parity) for solving spectrum of fermionic systems.
 
-See also the docstring : [`spectrum`](@ref)
+See also the docstring : 
+```@docs
+DensityOfStates(M::AbstractHEOMLSMatrix, ρ, op, ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
+```
 
 ```julia
 Hs::AbstractMatrix # system Hamiltonian
@@ -108,5 +144,5 @@ d::AbstractMatrix
 # the spectrum value for the specific frequency ω which need to be solved
 ω_list = 0:0.5:2 # [0.0, 0.5, 1.0, 1.5, 2.0]
 
-πAω = spectrum(M_odd, ados, d, ω_list)
+πAω = DensityOfStates(M_odd, ados, d, ω_list)
 ```

docs/src/stationary_state.md

@@ -13,7 +13,9 @@ The first method is implemented by solving the linear problem
 
 See the docstring of this method:  
 
-[`SteadyState(M::AbstractHEOMLSMatrix)`](@ref)
+```@docs
+SteadyState(M::AbstractHEOMLSMatrix; solver=UMFPACKFactorization(), verbose::Bool=true, SOLVEROptions...)
+```
 
 ```julia
 # the HEOMLS matrix
@@ -36,7 +38,9 @@ until finding a stationary solution.
 ### Given the initial state as Density Operator (`AbstractMatrix` type)
 See the docstring of this method:  
 
-[`SteadyState(M::AbstractHEOMLSMatrix, ρ0)`](@ref)
+```@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...)
+```
 
 ```julia
 # the HEOMLS matrix
@@ -50,7 +54,9 @@ ados_steady = SteadyState(M, ρ0)
 
 ### Given the initial state as Auxiliary Density Operators
 See the docstring of this method:  
-[`SteadyState(M::AbstractHEOMLSMatrix, ados::ADOs)`](@ref)
+```@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...)
+```
 
 ```julia
 # the HEOMLS matrix