Merge pull request #100 from ytdHuang/fix/error

Fix operator types

Project.toml

@@ -1,7 +1,7 @@
 name = "HierarchicalEOM"
 uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
 authors = ["Yi-Te Huang <[email protected]>"]
-version = "2.1.2"
+version = "2.1.3"
 
 [deps]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"

docs/make.jl

@@ -73,10 +73,7 @@ const PAGES = Any[
         "Stationary State"=>"stationary_state.md",
         "Spectrum"=>"spectrum.md",
         "Examples"=>EX_output_files,
-        "Solvers Lists"=>Any[
-            "ODE_solvers.md",
-            "LS_solvers.md",
-        ],
+        "Solvers Lists"=>Any["ODE_solvers.md", "LS_solvers.md"],
         "Extensions"=>Any["CUDA.jl"=>"extensions/CUDA.md"],
     ],
     "Library API"=>"libraryAPI.md",

docs/src/bath_boson/bosonic_bath_intro.md

@@ -52,7 +52,7 @@ One can check the information of the [`BosonBath`](@ref) by the `print` function
 print(bath)
 ```
 ```
-BosonBath object with (system) dims = [2] and 4 exponential-expansion terms
+BosonBath object with 4 exponential-expansion terms
 ```
 
 ## Calculate the correlation function

docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md

@@ -36,7 +36,7 @@ One can check the information of the [`BosonBath`](@ref) by the `print` function
 print(bath)
 ```
 ```
-BosonBath object with (system) dim = 2 and 4 exponential-expansion terms
+BosonBath object with 4 exponential-expansion terms
 ```
 Note that [`BosonBath`](@ref) under RWA always have even number of exponential terms (half for ``C^{\nu=+}`` and half for ``C^{\nu=-}``)
 

docs/src/bath_fermion/fermionic_bath_intro.md

@@ -39,7 +39,7 @@ One can check the information of the [`FermionBath`](@ref) by the `print` functi
 print(bath)
 ```
 ```
-FermionBath object with (system) dim = 2 and 4 exponential-expansion terms
+FermionBath object with 4 exponential-expansion terms
 ```
 Note that [`FermionBath`](@ref) always have even number of exponential terms (half for ``C^{\nu=+}`` and half for ``C^{\nu=-}``)
 

docs/src/spectrum.md

@@ -73,7 +73,7 @@ Finially, one can obtain the value of the power spectrum for specific ``\omega``
 
 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...)
+PowerSpectrum(M::AbstractHEOMLSMatrix, ρ::Union{QuantumObject,ADOs}, P_op::Union{QuantumObject,HEOMSuperOp}, Q_op::Union{QuantumObject,HEOMSuperOp}, ωlist::AbstractVector, reverse::Bool = false; solver = UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
 ```
 
 ```julia

src/Bath.jl

@@ -9,7 +9,7 @@ An object which describes a single exponential-expansion term (naively, an excit
 The expansion of a bath correlation function can be expressed as : ``C(t) = \sum_i \eta_i \exp(-\gamma_i t)``.
 
 # Fields
-- `op` : The system coupling operator according to system-bath interaction.
+- `op::QuantumObject` : The system coupling operator according to system-bath interaction.
 - `η::Number` : the coefficient ``\eta_i`` in bath correlation function.
 - `γ::Number` : the coefficient ``\gamma_i`` in bath correlation function.
 - `types::String` : The type-tag of the exponent.
@@ -24,7 +24,7 @@ The different types of the Exponent:
 - `\"fE\"` : from emission fermionic correlation function ``C^{\nu=-}(t)``
 """
 struct Exponent
-    op::Any
+    op::QuantumObject
     η::Number
     γ::Number
     types::String
@@ -33,16 +33,13 @@ end
 show(io::IO, E::Exponent) = print(io, "Bath Exponent with types = \"$(E.types)\", η = $(E.η), γ = $(E.γ).\n")
 show(io::IO, m::MIME"text/plain", E::Exponent) = show(io, E)
 
-show(io::IO, B::AbstractBath) =
-    print(io, "$(typeof(B)) object with (system) dims = $(B.dims) and $(B.Nterm) exponential-expansion terms\n")
+show(io::IO, B::AbstractBath) = print(io, "$(typeof(B)) object with $(B.Nterm) exponential-expansion terms\n")
 show(io::IO, m::MIME"text/plain", B::AbstractBath) = show(io, B)
 
-show(io::IO, B::AbstractBosonBath) =
-    print(io, "$(typeof(B))-type bath with (system) dims = $(B.dims) and $(B.Nterm) exponential-expansion terms\n")
+show(io::IO, B::AbstractBosonBath) = print(io, "$(typeof(B))-type bath with $(B.Nterm) exponential-expansion terms\n")
 show(io::IO, m::MIME"text/plain", B::AbstractBosonBath) = show(io, B)
 
-show(io::IO, B::AbstractFermionBath) =
-    print(io, "$(typeof(B))-type bath with (system) dims = $(B.dims) and $(B.Nterm) exponential-expansion terms\n")
+show(io::IO, B::AbstractFermionBath) = print(io, "$(typeof(B))-type bath with $(B.Nterm) exponential-expansion terms\n")
 show(io::IO, m::MIME"text/plain", B::AbstractFermionBath) = show(io, B)
 
 checkbounds(B::AbstractBath, i::Int) =
@@ -199,7 +196,6 @@ An object which describes the interaction between system and bosonic bath
 # Fields
 - `bath` : the different boson-bath-type objects which describes the interaction between system and bosonic bath
 - `op` : The system coupling operator, must be Hermitian and, for fermionic systems, even-parity to be compatible with charge conservation.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
 - `Nterm` : the number of exponential-expansion term of correlation functions
 - `δ` : The approximation discrepancy which is used for adding the terminator to HEOM matrix (see function: addTerminator)
 
@@ -218,7 +214,6 @@ end
 struct BosonBath <: AbstractBath
     bath::Vector{AbstractBosonBath}
     op::QuantumObject
-    dims::SVector
     Nterm::Int
     δ::Number
 end
@@ -257,7 +252,7 @@ function BosonBath(
     else
         bRI = bosonRealImag(_op, real.(η), imag.(η), γ)
     end
-    return BosonBath(AbstractBosonBath[bRI], _op, _op.dims, bRI.Nterm, δ)
+    return BosonBath(AbstractBosonBath[bRI], _op, bRI.Nterm, δ)
 end
 
 @doc raw"""
@@ -281,7 +276,7 @@ C(\tau)=\sum_{u=\textrm{R},\textrm{I}}(\delta_{u, \textrm{R}} + i\delta_{u, \tex
 where ``\delta`` is the Kronecker delta function and ``C^{u}(\tau)=\sum_i \eta_i^u \exp(-\gamma_i^u \tau)``
 
 # Parameters
-- `op` : The system coupling operator, must be Hermitian and, for fermionic systems, even-parity to be compatible with charge conservation.
+- `op::QuantumObject` : The system coupling operator, must be Hermitian and, for fermionic systems, even-parity to be compatible with charge conservation.
 - `η_real::Vector{Ti<:Number}` : the coefficients ``\eta_i`` in real part of bath correlation function ``C^{u=\textrm{R}}``.
 - `γ_real::Vector{Tj<:Number}` : the coefficients ``\gamma_i`` in real part of bath correlation function ``C^{u=\textrm{R}}``.
 - `η_imag::Vector{Tk<:Number}` : the coefficients ``\eta_i`` in imaginary part of bath correlation function ``C^{u=\textrm{I}}``.
@@ -290,7 +285,7 @@ where ``\delta`` is the Kronecker delta function and ``C^{u}(\tau)=\sum_i \eta_i
 - `combine::Bool` : Whether to combine the exponential-expansion terms with the same frequency. Defaults to `true`.
 """
 function BosonBath(
-    op,
+    op::QuantumObject,
     η_real::Vector{Ti},
     γ_real::Vector{Tj},
     η_imag::Vector{Tk},
@@ -342,12 +337,12 @@ function BosonBath(
             error("Conflicts occur in combining real and imaginary parts of bath correlation function.")
         end
 
-        return BosonBath(AbstractBosonBath[bR, bI, bRI], _op, _op.dims, Nterm_new, δ)
+        return BosonBath(AbstractBosonBath[bR, bI, bRI], _op, Nterm_new, δ)
 
     else
         bR = bosonReal(_op, η_real, γ_real)
         bI = bosonImag(_op, η_imag, γ_imag)
-        return BosonBath(AbstractBosonBath[bR, bI], _op, _op.dims, bR.Nterm + bI.Nterm, δ)
+        return BosonBath(AbstractBosonBath[bR, bI], _op, bR.Nterm + bI.Nterm, δ)
     end
 end
 
@@ -524,7 +519,7 @@ function BosonBathRWA(
 
     bA = bosonAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
     bE = bosonEmit(_op, η_emit, γ_emit, η_absorb)
-    return BosonBath(AbstractBosonBath[bA, bE], _op, _op.dims, bA.Nterm + bE.Nterm, δ)
+    return BosonBath(AbstractBosonBath[bA, bE], _op, bA.Nterm + bE.Nterm, δ)
 end
 
 @doc raw"""
@@ -656,7 +651,6 @@ An object which describes the interaction between system and fermionic bath
 # Fields
 - `bath` : the different fermion-bath-type objects which describes the interaction
 - `op` : The system \"emission\" operator according to the system-fermionic-bath interaction.
-- `dims` : the dimension list of the coupling operator (should be equal to the system dims).
 - `Nterm` : the number of exponential-expansion term of correlation functions
 - `δ` : The approximation discrepancy which is used for adding the terminator to HEOM matrix (see function: addTerminator)
 
@@ -675,7 +669,6 @@ end
 struct FermionBath <: AbstractBath
     bath::Vector{AbstractFermionBath}
     op::QuantumObject
-    dims::SVector
     Nterm::Int
     δ::Number
 end
@@ -717,7 +710,7 @@ function FermionBath(
     _op = HandleMatrixType(op, "op (coupling operator)")
     fA = fermionAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
     fE = fermionEmit(_op, η_emit, γ_emit, η_absorb)
-    return FermionBath(AbstractFermionBath[fA, fE], _op, _op.dims, fA.Nterm + fE.Nterm, δ)
+    return FermionBath(AbstractFermionBath[fA, fE], _op, fA.Nterm + fE.Nterm, δ)
 end
 
 @doc raw"""

src/HeomBase.jl

@@ -12,7 +12,7 @@ function HandleMatrixType(M::QuantumObject, MatrixName::String = ""; type::Quant
     if M.type == type
         return M
     else
-        error("The matrix $(MatrixName) should be an $(Type).")
+        error("The matrix $(MatrixName) should be an $(type).")
     end
 end
 function HandleMatrixType(M::QuantumObject, dims::SVector, MatrixName::String = ""; type::QuantumObjectType = Operator)

src/density_of_states.jl

@@ -8,8 +8,8 @@ Calculate density of states for the fermionic system in frequency domain.
 
 # Parameters
 - `M::AbstractHEOMLSMatrix` : the HEOMLS matrix which acts on `ODD`-parity operators.
-- `ρ` :  the system density matrix or the auxiliary density operators.
-- `d_op` : The annihilation operator (``d`` as shown above) acting on the fermionic system.
+- `ρ::Union{QuantumObject,ADOs}` :  the system density matrix or the auxiliary density operators.
+- `d_op::QuantumObject` : The annihilation operator (``d`` as shown above) acting on the fermionic system.
 - `ωlist::AbstractVector` : the specific frequency points to solve.
 - `solver::SciMLLinearSolveAlgorithm` : solver in package `LinearSolve.jl`. Default to `UMFPACKFactorization()`.
 - `verbose::Bool` : To display verbose output and progress bar during the process or not. Defaults to `true`.
@@ -24,7 +24,7 @@ Calculate density of states for the fermionic system in frequency domain.
 """
 @noinline function DensityOfStates(
     M::AbstractHEOMLSMatrix,
-    ρ,
+    ρ::Union{QuantumObject,ADOs},
     d_op::QuantumObject,
     ωlist::AbstractVector;
     solver::SciMLLinearSolveAlgorithm = UMFPACKFactorization(),

src/heom_matrices/HierarchyDict.jl

@@ -98,7 +98,7 @@ end
     if T == BosonBath
         baths = AbstractBosonBath[]
         for (α, b) in enumerate(B)
-            if b.dims != dims
+            if b.op.dims != dims
                 error("The matrix size of the bosonic bath coupling operators are not consistent.")
             end
             push!(baths, b.bath...)
@@ -112,7 +112,7 @@ end
     elseif T == FermionBath
         baths = AbstractFermionBath[]
         for (α, b) in enumerate(B)
-            if b.dims != dims
+            if b.op.dims != dims
                 error("The matrix size of the fermionic bath coupling operators are not consistent.")
             end
             push!(baths, b.bath...)
@@ -182,7 +182,7 @@ end
     bosonPtr = Tuple[]
     baths_b = AbstractBosonBath[]
     for (α, b) in enumerate(bB)
-        if b.dims != dims
+        if b.op.dims != dims
             error("The matrix size of the bosonic bath coupling operators are not consistent.")
         end
         push!(baths_b, b.bath...)
@@ -199,7 +199,7 @@ end
     fermionPtr = Tuple[]
     baths_f = AbstractFermionBath[]
     for (α, b) in enumerate(fB)
-        if b.dims != dims
+        if b.op.dims != dims
             error("The matrix size of the fermionic bath coupling operators are not consistent.")
         end
         push!(baths_f, b.bath...)

src/heom_matrices/heom_matrix_base.jl

@@ -358,8 +358,8 @@ function addTerminator(M::Mtype, Bath::Union{BosonBath,FermionBath}) where {Mtyp
         error("The type of input HEOMLS matrix does not support this functionality.")
     end
 
-    if M.dims != Bath.dims
-        error("The system dimension between the HEOMLS matrix and Bath are not consistent.")
+    if M.dims != Bath.op.dims
+        error("The system dims between the HEOMLS matrix and Bath coupling operator are not consistent.")
     end
 
     if Bath.δ == 0

src/power_spectrum.jl

@@ -5,31 +5,28 @@ Calculate power spectrum for the system in frequency domain where `P_op` will be
 This function is equivalent to:
 `PowerSpectrum(M, ρ, Q_op', Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...)`
 """
-function PowerSpectrum(
+PowerSpectrum(
     M::AbstractHEOMLSMatrix,
-    ρ,
-    Q_op,
+    ρ::Union{QuantumObject,ADOs},
+    Q_op::QuantumObject,
     ωlist::AbstractVector,
     reverse::Bool = false;
     solver::SciMLLinearSolveAlgorithm = UMFPACKFactorization(),
     verbose::Bool = true,
     filename::String = "",
     SOLVEROptions...,
+) = PowerSpectrum(
+    M,
+    ρ,
+    Q_op',
+    Q_op,
+    ωlist,
+    reverse;
+    solver = solver,
+    verbose = verbose,
+    filename = filename,
+    SOLVEROptions...,
 )
-    _Q_op = HandleMatrixType(Q_op, M.dims)
-    return PowerSpectrum(
-        M,
-        ρ,
-        _Q_op',
-        _Q_op,
-        ωlist,
-        reverse;
-        solver = solver,
-        verbose = verbose,
-        filename = filename,
-        SOLVEROptions...,
-    )
-end
 
 @doc raw"""
     PowerSpectrum(M, ρ, P_op, Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...)
@@ -49,9 +46,9 @@ remember to set the parameters:
 
 # Parameters
 - `M::AbstractHEOMLSMatrix` : the HEOMLS matrix.
-- `ρ` :  the system density matrix or the auxiliary density operators.
-- `P_op`: the system operator (or `HEOMSuperOp`) ``P`` acting on the system.
-- `Q_op`: the system operator (or `HEOMSuperOp`) ``Q`` acting on the system.
+- `ρ::Union{QuantumObject,ADOs}` :  the system density matrix or the auxiliary density operators.
+- `P_op::Union{QuantumObject,HEOMSuperOp}`: the system operator (or `HEOMSuperOp`) ``P`` acting on the system.
+- `Q_op::Union{QuantumObject,HEOMSuperOp}`: the system operator (or `HEOMSuperOp`) ``Q`` acting on the system.
 - `ωlist::AbstractVector` : the specific frequency points to solve.
 - `reverse::Bool` : If `true`, calculate ``\langle P(-t)Q(0) \rangle = \langle P(0)Q(t) \rangle = \langle P(t)Q(0) \rangle^*`` instead of ``\langle P(t) Q(0) \rangle``. Default to `false`.
 - `solver::SciMLLinearSolveAlgorithm` : solver in package `LinearSolve.jl`. Default to `UMFPACKFactorization()`.
@@ -67,9 +64,9 @@ remember to set the parameters:
 """
 @noinline function PowerSpectrum(
     M::AbstractHEOMLSMatrix,
-    ρ,
-    P_op,
-    Q_op,
+    ρ::Union{QuantumObject,ADOs},
+    P_op::Union{QuantumObject,HEOMSuperOp},
+    Q_op::Union{QuantumObject,HEOMSuperOp},
     ωlist::AbstractVector,
     reverse::Bool = false;
     solver::SciMLLinearSolveAlgorithm = UMFPACKFactorization(),

test/HEOMSuperOp.jl

@@ -95,6 +95,7 @@
     ## test for ERRORs
     J_wrong1 = HEOMSuperOp(γ_eR * d, ODD, M_me, "L")
     J_wrong2 = HEOMSuperOp(γ_eR * d, EVEN, M_heom, "L")
+    @test_throws ErrorException HEOMSuperOp(d, EVEN, M_heom, "LR")
     @test_throws ErrorException J_me * J_wrong2
     @test_throws ErrorException J_me + J_wrong1
     @test_throws ErrorException J_me + J_wrong2

test/power_spectrum.jl

@@ -54,7 +54,6 @@
     a_odd = HEOMSuperOp(a, ODD, L)
     bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2)
     @test_throws ErrorException PowerSpectrum(L, ados_s, a, ωlist; verbose = false, filename = "PSD")
-    @test_throws ErrorException PowerSpectrum(L, ados_s, a_even, ωlist; verbose = false)
     @test_throws ErrorException PowerSpectrum(L, ados_s, a_even, a_odd, ωlist; verbose = false)
     @test_throws ErrorException PowerSpectrum(L, ados_s, mat2, ωlist; verbose = false)
     @test_throws ErrorException PowerSpectrum(L, ADOs(zeros(8), 2), a, ωlist; verbose = false)