Make HEOMLS matrices' type parametric

docs/src/libraryAPI.md

@@ -60,6 +60,8 @@ M_Fermion(Hsys, tier::Int, Bath::Vector{FermionBath}, parity::Symbol=:even; thre
 M_Boson_Fermion
 M_Boson_Fermion(Hsys, tier_b::Int, tier_f::Int, Bath_b::Vector{BosonBath}, Bath_f::Vector{FermionBath}, parity::Symbol=:even; threshold::Real=0.0, verbose::Bool=true)
 size(M::AbstractHEOMLSMatrix)
+size(M::AbstractHEOMLSMatrix, dim::Int)
+eltype(M::AbstractHEOMLSMatrix)
 Propagator
 addBosonDissipator
 addFermionDissipator
@@ -71,6 +73,7 @@ addTerminator
 ADOs
 ADOs(V::AbstractVector, N::Int)
 length(A::ADOs)
+eltype(A::ADOs)
 getRho
 getADO
 Expect

src/ADOs.jl

@@ -39,7 +39,7 @@
 """
 function ADOs(V::AbstractVector, N::Int, parity::Symbol=:even)
     # check the dimension of V
-    d,  = size(V)
+    d  = size(V, 1)
     dim = √(d / N)
     if isinteger(dim)
         if (parity == :even) || (parity == :odd)
@@ -64,6 +64,12 @@
 """
 length(A::ADOs) = A.N
 
+@doc raw"""
+    eltype(A::ADOs)
+Returns the elements' type of the Auxiliary Density Operators (ADOs)
+"""
+eltype(A::ADOs) = eltype(A.data)
+
 lastindex(A::ADOs) = length(A)
 
 function getindex(A::ADOs, i::Int)

src/Bath.jl

@@ -155,7 +155,7 @@
 
 function _check_bosonic_coupling_operator(op)
     _op = HandleMatrixType(op, 0, "op (coupling operator)")
-    N,  = size(_op)
+    N   = size(_op, 1)
     if !ishermitian(_op)
         @warn "The system-bosonic-bath coupling operator \"op\" should be Hermitian operator."
     end
@@ -164,7 +164,7 @@
 
 function _check_bosonic_RWA_coupling_operator(op)
     _op = HandleMatrixType(op, 0, "op (coupling operator)")
-    N,  = size(_op)
+    N  = size(_op, 1)
     return N, _op
 end
 
@@ -183,7 +183,7 @@
 
 function _check_fermionic_coupling_operator(op)
     _op = HandleMatrixType(op, 0, "op (coupling operator)")
-    N,  = size(_op)
+    N  = size(_op, 1)
     return N, _op
 end
 

src/HierarchicalEOM.jl

@@ -48,7 +48,7 @@ module HierarchicalEOM
     # sub-module HeomAPI for HierarchicalEOM
     module HeomAPI
         using ..Bath
-        import Base: ==, show, length, size, getindex, keys, setindex!, lastindex, iterate, checkbounds, hash, copy
+        import Base: ==, show, length, size, getindex, keys, setindex!, lastindex, iterate, checkbounds, hash, copy, eltype
         import Base.Threads: @threads, threadid, nthreads, lock, unlock, SpinLock
         import LinearAlgebra: I, kron, tr
         import SparseArrays: sparse, spzeros, sparsevec, reshape, SparseVector, SparseMatrixCSC, AbstractSparseMatrix

src/Spectrum.jl

@@ -44,7 +44,7 @@ function spectrum(
         SOLVEROptions...
     )
 
-    Size, = size(M)
+    Size = size(M, 1)
 
     # check ρ
     if typeof(ρ) == ADOs  # ρ::ADOs
@@ -106,7 +106,7 @@ end
         end
     end
 
-    Size, = size(M)
+    Size = size(M, 1)
     I_total = sparse(I, Size, Size)
     I_heom  = sparse(I, M.N, M.N)
 
@@ -175,7 +175,7 @@ end
         end
     end
 
-    Size, = size(M)
+    Size = size(M, 1)
     I_total = sparse(I, Size, Size)
     I_heom  = sparse(I, M.N, M.N)
 

src/SteadyState.jl

@@ -20,7 +20,7 @@ For more details about solvers and extra options, please refer to [`LinearSolve.
     end    
 
     A = copy(M.data)
-    S, = size(A)
+    S = size(A, 1)
     A[1,1:S] .= 0
 
     # sparse(row_idx, col_idx, values, row_dims, col_dims)

src/heom_matrices/M_Boson.jl

@@ -1,9 +1,9 @@
 @doc raw"""
-    struct M_Boson <: AbstractHEOMLSMatrix
+    struct M_Boson{T} <: AbstractHEOMLSMatrix
 HEOM Liouvillian superoperator matrix for bosonic bath
 
 # Fields
-- `data` : the sparse matrix of HEOM Liouvillian superoperator
+- `data::T` : the sparse matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the bosonic hierarchy
 - `dim` : the dimension of system
 - `N` : the number of total ADOs
@@ -12,8 +12,8 @@ HEOM Liouvillian superoperator matrix for bosonic bath
 - `bath::Vector{BosonBath}` : the vector which stores all `BosonBath` objects
 - `hierarchy::HierarchyDict`: the object which contains all dictionaries for boson-bath-ADOs hierarchy.
 """
-struct M_Boson <: AbstractHEOMLSMatrix
-    data::SparseMatrixCSC{ComplexF64, Int64}
+struct M_Boson{T} <: AbstractHEOMLSMatrix
+    data::T
     tier::Int
     dim::Int
     N::Int
@@ -59,7 +59,7 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, 0, "Hsys (system Hamiltonian)")
-    Nsys,   = size(_Hsys)
+    Nsys    = size(_Hsys, 1)
     sup_dim = Nsys ^ 2
     I_sup   = sparse(I, sup_dim, sup_dim)
 
@@ -147,5 +147,5 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
         println("[DONE]")
         flush(stdout)
     end
-    return M_Boson(L_he, tier, Nsys, Nado, sup_dim, parity, Bath, hierarchy)
+    return M_Boson{SparseMatrixCSC{ComplexF64, Int64}}(L_he, tier, Nsys, Nado, sup_dim, parity, Bath, hierarchy)
 end

src/heom_matrices/M_Boson_Fermion.jl

@@ -1,9 +1,9 @@
 @doc raw"""
-    struct M_Boson_Fermion <: AbstractHEOMLSMatrix
+    struct M_Boson_Fermion{T} <: AbstractHEOMLSMatrix
 HEOM Liouvillian superoperator matrix for mixtured (bosonic and fermionic) bath 
 
 # Fields
-- `data` : the sparse matrix of HEOM Liouvillian superoperator
+- `data::T` : the sparse matrix of HEOM Liouvillian superoperator
 - `Btier` : the tier (cutoff level) for bosonic hierarchy
 - `Ftier` : the tier (cutoff level) for fermionic hierarchy
 - `dim` : the dimension of system
@@ -14,8 +14,8 @@ HEOM Liouvillian superoperator matrix for mixtured (bosonic and fermionic) bath
 - `Fbath::Vector{FermionBath}` : the vector which stores all `FermionBath` objects
 - `hierarchy::MixHierarchyDict`: the object which contains all dictionaries for mixed-bath-ADOs hierarchy.
 """
-struct M_Boson_Fermion <: AbstractHEOMLSMatrix
-    data::SparseMatrixCSC{ComplexF64, Int64}
+struct M_Boson_Fermion{T} <: AbstractHEOMLSMatrix
+    data::T
     Btier::Int
     Ftier::Int
     dim::Int
@@ -76,7 +76,7 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, 0, "Hsys (system Hamiltonian)")
-    Nsys,   = size(_Hsys)
+    Nsys    = size(_Hsys, 1)
     sup_dim = Nsys ^ 2
     I_sup   = sparse(I, sup_dim, sup_dim)
 
@@ -196,5 +196,5 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
         println("[DONE]") 
         flush(stdout)
     end
-    return M_Boson_Fermion(L_he, Btier, Ftier, Nsys, Nado, sup_dim, parity, Bbath, Fbath, hierarchy)
+    return M_Boson_Fermion{SparseMatrixCSC{ComplexF64, Int64}}(L_he, Btier, Ftier, Nsys, Nado, sup_dim, parity, Bbath, Fbath, hierarchy)
 end

src/heom_matrices/M_Fermion.jl

@@ -1,9 +1,9 @@
 @doc raw"""
-    struct M_Fermion <: AbstractHEOMLSMatrix
+    struct M_Fermion{T} <: AbstractHEOMLSMatrix
 HEOM Liouvillian superoperator matrix for fermionic bath
 
 # Fields
-- `data` : the sparse matrix of HEOM Liouvillian superoperator
+- `data::T` : the sparse matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the fermionic hierarchy
 - `dim` : the dimension of system
 - `N` : the number of total ADOs
@@ -12,8 +12,8 @@ HEOM Liouvillian superoperator matrix for fermionic bath
 - `bath::Vector{FermionBath}` : the vector which stores all `FermionBath` objects
 - `hierarchy::HierarchyDict`: the object which contains all dictionaries for fermion-bath-ADOs hierarchy.
 """
-struct M_Fermion <: AbstractHEOMLSMatrix
-    data::SparseMatrixCSC{ComplexF64, Int64}
+struct M_Fermion{T} <: AbstractHEOMLSMatrix
+    data::T
     tier::Int
     dim::Int
     N::Int
@@ -57,7 +57,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, 0, "Hsys (system Hamiltonian)")
-    Nsys,   = size(_Hsys)
+    Nsys    = size(_Hsys, 1)
     sup_dim = Nsys ^ 2
     I_sup   = sparse(I, sup_dim, sup_dim)
 
@@ -144,5 +144,5 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
         println("[DONE]")
         flush(stdout)
     end
-    return M_Fermion(L_he, tier, Nsys, Nado, sup_dim, parity, Bath, hierarchy)
+    return M_Fermion{SparseMatrixCSC{ComplexF64, Int64}}(L_he, tier, Nsys, Nado, sup_dim, parity, Bath, hierarchy)
 end

src/heom_matrices/M_S.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    struct M_S <: AbstractHEOMLSMatrix
+    struct M_S{T} <: AbstractHEOMLSMatrix
 HEOM Liouvillian superoperator matrix with cutoff level of the hierarchy equals to `0`.  
 This corresponds to the standard Schrodinger (Liouville-von Neumann) equation, namely
 ```math
@@ -8,15 +8,15 @@ M[\cdot]=-i \left[H_{sys}, \cdot \right]_-,
 where ``[\cdot, \cdot]_-`` stands for commutator.
 
 # Fields
-- `data` : the sparse matrix of HEOM Liouvillian superoperator
+- `data::T` : the sparse matrix of HEOM Liouvillian superoperator
 - `tier` : the tier (cutoff level) for the hierarchy, which equals to `0` in this case
 - `dim` : the dimension of system
 - `N` : the number of total ADOs, which equals to `1` (only the reduced density operator) in this case
 - `sup_dim` : the dimension of system superoperator
 - `parity` : the parity label of the fermionic system (usually `:even`, only set as `:odd` for calculating spectrum of fermionic system)
 """
-struct M_S <: AbstractHEOMLSMatrix
-    data::SparseMatrixCSC{ComplexF64, Int64}
+struct M_S{T} <: AbstractHEOMLSMatrix
+    data::T
     tier::Int
     dim::Int
     N::Int
@@ -53,7 +53,7 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
 
     # check for system dimension
     _Hsys = HandleMatrixType(Hsys, 0, "Hsys (system Hamiltonian)")
-    Nsys,   = size(_Hsys)
+    Nsys    = size(_Hsys, 1)
     sup_dim = Nsys ^ 2
 
     # the Liouvillian operator for free Hamiltonian
@@ -66,5 +66,5 @@ Note that the parity only need to be set as `:odd` when the system contains ferm
         println("[DONE]")
         flush(stdout)
     end
-    return M_S(Lsys, 0, Nsys, 1, sup_dim, parity)
+    return M_S{SparseMatrixCSC{ComplexF64, Int64}}(Lsys, 0, Nsys, 1, sup_dim, parity)
 end

src/heom_matrices/heom_matrix_base.jl

@@ -10,15 +10,27 @@
 """
 size(M::AbstractHEOMLSMatrix) = size(M.data)
 
+@doc raw"""
+    size(M::AbstractHEOMLSMatrix, dim::Int)
+Returns the specified dimension of the HEOM Liouvillian superoperator matrix
+"""
+size(M::AbstractHEOMLSMatrix, dim::Int) = size(M.data, dim)
+
+@doc raw"""
+    eltype(M::AbstractHEOMLSMatrix)
+Returns the elements' type of the HEOM Liouvillian superoperator matrix
+"""
+eltype(M::AbstractHEOMLSMatrix) = eltype(M.data)
+
 getindex(M::AbstractHEOMLSMatrix, i::Ti, j::Tj) where {Ti, Tj <: Any} = M.data[i, j]
 
 function show(io::IO, M::AbstractHEOMLSMatrix)
     T = typeof(M)
-    if T == M_S
+    if T <: M_S
         type = "Schrodinger Eq."
-    elseif T == M_Boson
+    elseif T <: M_Boson
         type = "Boson"
-    elseif T == M_Fermion
+    elseif T <: M_Fermion
         type = "Fermion"
     else
         type = "Boson-Fermion"
@@ -85,11 +97,11 @@
             L += spre(_J) * spost(_J') - 0.5 * (spre(_J' * _J) + spost(_J' * _J))
         end
 
-        if T == M_S
+        if T <: M_S
             return M_S(M.data + L, M.tier, M.dim, M.N, M.sup_dim, M.parity)
-        elseif T == M_Boson
+        elseif T <: M_Boson
             return M_Boson(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
-        elseif T == M_Fermion
+        elseif T <: M_Fermion
             return M_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
         else
             return M_Boson_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.Btier, M.Ftier, M.dim, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy)
@@ -128,11 +140,11 @@
             _J = HandleMatrixType(J, M.dim, "in jumpOP")
             L += ((-1) ^ parity) * spre(_J) * spost(_J') - 0.5 * (spre(_J' * _J) + spost(_J' * _J))
         end
-        if T == M_S
+        if T <: M_S
             return M_S(M.data + L, M.tier, M.dim, M.N, M.sup_dim, M.parity)
-        elseif T == M_Boson
+        elseif T <: M_Boson
             return M_Boson(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
-        elseif T == M_Fermion
+        elseif T <: M_Fermion
             return M_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
         else
             return M_Boson_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.Btier, M.Ftier, M.dim, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy)
@@ -163,11 +175,11 @@
 """
 function addTerminator(M::Mtype, Bath::Union{BosonBath, FermionBath}) where Mtype <: AbstractHEOMLSMatrix
     Btype = typeof(Bath)
-    if (Btype == BosonBath) && (Mtype == M_Fermion)
+    if (Btype == BosonBath) && (Mtype <: M_Fermion)
         error("For $(Btype), the type of HEOMLS matrix should be either M_Boson or M_Boson_Fermion.")
-    elseif (Btype == FermionBath) && (Mtype == M_Boson)
+    elseif (Btype == FermionBath) && (Mtype <: M_Boson)
         error("For $(Btype), the type of HEOMLS matrix should be either M_Fermion or M_Boson_Fermion.")
-    elseif Mtype == M_S
+    elseif Mtype <: M_S
         error("The type of input HEOMLS matrix does not support this functionality.")
     end
 
@@ -185,9 +197,9 @@
             spre(J) * spost(J') - 0.5 * (spre(J' * J) + spost(J' * J))
         )
 
-        if Mtype == M_Boson
+        if Mtype <: M_Boson
             return M_Boson(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
-        elseif Mtype == M_Fermion
+        elseif Mtype <: M_Fermion
             return M_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy)
         else
             return M_Boson_Fermion(M.data + kron(sparse(I, M.N, M.N), L), M.Btier, M.Ftier, M.dim, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy)

test/heom_api.jl

@@ -11,6 +11,7 @@
     ados = ados_list[end]
     @test ados.dim == L.dim
     @test length(ados) == L.N
+    @test eltype(L) == eltype(ados)
     ρ0 = ados[1]
     @test getRho(ados) == ρ0
     ρ1 = [
@@ -76,6 +77,7 @@ J = [0 0.1450 - 0.7414im; 0.1450 + 0.7414im 0]
     ados = SteadyState(L; verbose=false)
     @test ados.dim == L.dim
     @test length(ados) == L.N
+    @test eltype(L) == eltype(ados)
     ρ0 = ados[1]
     @test getRho(ados) == ρ0
     ρ1 = [
@@ -147,6 +149,7 @@ end
     ados = SteadyState(L; verbose=false)
     @test ados.dim == L.dim
     @test length(ados) == L.N
+    @test eltype(L) == eltype(ados)
     ρ0 = ados[1]
     @test getRho(ados) == ρ0
     ρ1 = [
@@ -213,6 +216,7 @@ end
     ados = SteadyState(L; verbose=false)
     @test ados.dim == L.dim
     @test length(ados) == L.N
+    @test eltype(L) == eltype(ados)
     ρ0 = ados[1]
     @test getRho(ados) == ρ0
     ρ1 = [