Merge pull request #44 Develop HEOMSuperOp

Develop HEOM superoperator for users

.github/workflows/Runtests.yml

@@ -7,6 +7,9 @@ on:
   pull_request:
     branches:
       - 'main'
+    types:
+      - 'synchronize'
+      - 'ready_for_review'
 
 jobs:
   test:

.github/workflows/documentation.yml

@@ -9,6 +9,9 @@ on:
   pull_request:
     branches:
       - 'main'
+    types:
+      - 'synchronize'
+      - 'ready_for_review'
 
 jobs:
   build:

docs/make.jl

@@ -70,6 +70,7 @@ const PAGES = Any[
             "HEOMLS for Bosonic and Fermionic Bath" => "heom_matrix/M_Boson_Fermion.md",
             "HEOMLS for Master Equation" => "heom_matrix/master_eq.md",
         ],
+        "Parity Support" => "Parity.md",
         "Hierarchy Dictionary" => "hierarchy_dictionary.md",
         "Time Evolution" => "time_evolution.md",
         "Stationary State" => "stationary_state.md",

docs/src/Parity.md

@@ -0,0 +1,41 @@
+# [Parity Support](@id doc-Parity)
+## Introduction
+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 which ``\hat{\mathcal{M}}`` is acting on. This dependence intuitively originates from the properties of partial traces over composite fermionic spaces, where operators do not necessarily commute. 
+
+As an example, for an environment made out of a single fermion, the reduced matrix elements ``\langle{i}|\rho_\textrm{s}^p|{j}\rangle`` (in a basis  labeled by ``\langle i|`` and ``|{j}\rangle``) involve the perturbative sum of expressions of the form ``\langle i| (c \tilde{\rho}_\textrm{e} \tilde{\rho}_\textrm{s}^p c^\dagger+\tilde{\rho}_\textrm{e} \tilde{\rho}_\textrm{s}^p)|{j}\rangle`` (in terms of environmental operators ``\tilde{\rho}_{\textrm{e}}``, system operators ``\tilde{\rho}_\textrm{s}^p`` with parity ``p``, and the environment-annihilation operator ``c``). These quantities depend on the commutator between ``\tilde{\rho}_\textrm{s}^p`` and ``c``, which is trivial only for `EVEN`-parity (``p=+``). In the `ODD`-parity (``p=-``) case, the partial trace over the environment requires further anti-commutations, ultimately resulting in extra minus signs in the expression for the effective propagator describing the reduced dynamics. 
+
+It is important to explicitly note that, here, by `parity` we do not refer to the presence of an odd or even number of fermions in the system but, rather, to the number of fermionic (annihilation or creation) operators needed to represent ``\rho_\textrm{s}^p``. The reduced density matrix of the system should be an `EVEN`-parity operator and can be expressed as ``\rho_{\textrm{s}}^{p=+}(t)``. However, there are some situations (for example, [calculating density of states for fermionic systems](@ref doc-DOS)) where ``\hat{\mathcal{M}}`` is acting on `ODD`-parity ADOs, e.g., ``\rho_{\textrm{s}}^{p=-}(t)=d_{\textrm{s}}\rho_{\textrm{s}}^{+}(t)`` or ``\rho_{\textrm{s}}^{p=-}(t)=d_{\textrm{s}}^\dagger\rho_{\textrm{s}}^{+}(t)``, where ``d_{\textrm{s}}`` is an annihilation operator acting on fermionic systems.
+
+## Parity support for HEOMLS
+One can specify the parameter `parity::AbstractParity` in the function of constructing ``\hat{\mathcal{M}}`` which describes the dynamics of [`EVEN`](@ref)- or [`ODD`](@ref)-parity auxiliary density operators (ADOs). The default value of the parameter is `parity=EVEN`.
+```julia
+Hs::AbstractMatrix  # system Hamiltonian
+Bbath::BosonBath    # bosonic   bath object
+Fbath::FermionBath  # fermionic bath object
+Btier::Int          # bosonic   truncation level 
+Ftier::Int          # fermionic truncation level 
+
+# create HEOMLS matrix in EVEN or ODD parity
+M_even = M_S(Hs, EVEN)
+M_odd  = M_S(Hs, ODD)
+
+M_even = M_Boson(Hs, Btier, Bbath, EVEN) 
+M_odd  = M_Boson(Hs, Btier, Bbath, ODD) 
+
+M_even = M_Fermion(Hs, Ftier, Fbath, EVEN) 
+M_odd  = M_Fermion(Hs, Ftier, Fbath, ODD)
+
+M_even = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, EVEN) 
+M_odd  = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, ODD) 
+```
+
+## Base functions support
+### Multiplication between Parity labels
+```julia
+EVEN * EVEN # gives EVEN
+EVEN * ODD  # gives ODD
+ODD  * EVEN # gives ODD
+ODD  * ODD  # gives EVEN
+!EVEN       # gives ODD
+!ODD        # gives EVEN
+```

docs/src/heom_matrix/HEOMLS_intro.md

@@ -49,31 +49,6 @@ M = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath; threshold=1e-7)
 !!! note "Default value of importance threshold"
     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 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
-Hs::AbstractMatrix  # system Hamiltonian
-Bbath::BosonBath    # bosonic   bath object
-Fbath::FermionBath  # fermionic bath object
-Btier::Int          # bosonic   truncation level 
-Ftier::Int          # fermionic truncation level 
-
-# create HEOMLS matrix in EVEN or ODD parity
-M_even = M_S(Hs, EVEN)
-M_odd  = M_S(Hs, ODD)
-
-M_even = M_Boson(Hs, Btier, Bbath, EVEN) 
-M_odd  = M_Boson(Hs, Btier, Bbath, ODD) 
-
-M_even = M_Fermion(Hs, Ftier, Fbath, EVEN) 
-M_odd  = M_Fermion(Hs, Ftier, Fbath, ODD)
-
-M_even = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, EVEN) 
-M_odd  = M_Boson_Fermion(Hs, Btier, Ftier, Bbath, Fbath, ODD) 
-```
-
 ## Methods
 All of the HEOMLS matrices supports the following two `Base` Functions :
  - `size(M::AbstractHEOMLSMatrix)` : Returns the size of the HEOMLS matrix.

docs/src/libraryAPI.md

@@ -59,6 +59,10 @@ ODD
 
 ## HEOM Liouvillian superoperator matrices
 ```@docs
+HEOMSuperOp
+HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix, mul_basis::AbstractString="L")
+HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString="L")
+HEOMSuperOp(op, opParity::AbstractParity, dim::Int, N::Int, mul_basis::AbstractString)
 M_S
 M_S(Hsys, parity::AbstractParity=EVEN; verbose::Bool=true)
 M_Boson
@@ -67,8 +71,11 @@ M_Fermion
 M_Fermion(Hsys, tier::Int, Bath::Vector{FermionBath}, parity::AbstractParity=EVEN; threshold::Real=0.0, verbose::Bool=true)
 M_Boson_Fermion
 M_Boson_Fermion(Hsys, 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)
 Propagator
 addBosonDissipator

docs/src/spectrum.md

@@ -123,7 +123,7 @@ Finially, one can obtain the density of states for specific ``\omega``, namely
 
 See also the docstring : 
 ```@docs
-DensityOfStates(M::AbstractHEOMLSMatrix, ρ, op, ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
+DensityOfStates(M::AbstractHEOMLSMatrix, ρ, d_op, ωlist::AbstractVector; solver=UMFPACKFactorization(), verbose::Bool = true, filename::String = "", SOLVEROptions...)
 ```
 
 ```julia

ext/HierarchicalEOM_CUDAExt.jl

@@ -1,8 +1,7 @@
 module HierarchicalEOM_CUDAExt
 
 using HierarchicalEOM
-import HierarchicalEOM.HeomAPI:  _HandleVectorType, _HandleSteadyStateMatrix
-import HierarchicalEOM.Spectrum: _HandleIdentityType
+import HierarchicalEOM.HeomBase: _HandleVectorType, _HandleSteadyStateMatrix, _HandleIdentityType
 import CUDA
 import CUDA: cu, CuArray
 import CUDA.CUSPARSE: CuSparseMatrixCSC

src/ADOs.jl

@@ -48,6 +48,20 @@ function ADOs(V::AbstractVector, N::Int, parity::AbstractParity=EVEN)
     end    
 end
 
+@doc raw"""
+    ADOs(ρ, N, parity)
+Gernerate the object of auxiliary density operators for HEOM model.
+
+# Parameters
+- `ρ` : the reduced density operator
+- `N::Int` : the number of auxiliary density operators.
+- `parity::AbstractParity` : the parity label (`EVEN` or `ODD`). Default to `EVEN`.
+"""
+function ADOs(ρ, N::Int=1, parity::AbstractParity=EVEN)
+    _ρ = sparsevec(HandleMatrixType(ρ, 0, "ρ"))
+    return ADOs(sparsevec(_ρ.nzind, _ρ.nzval, N * length(_ρ)), N, parity)
+end
+
 function checkbounds(A::ADOs, i::Int)
     if (i > A.N) || (i < 1)
         error("Attempt to access $(A.N)-element ADOs at index [$(i)]")
@@ -158,9 +172,13 @@ where ``O`` is the operator and ``\rho`` is the reduced density operator in the
 """
 function Expect(op, ados::ADOs; take_real=true)
 
-    _op = HandleMatrixType(op, ados.dim, "op (observable)")
-
-    exp_val = tr(_op * getRho(ados))
+    if typeof(op) == HEOMSuperOp
+        _check_sys_dim_and_ADOs_num(op, ados)
+        exp_val = Tr(ados.dim, ados.N) * (op * ados).data
+    else
+        _op = HandleMatrixType(op, ados.dim, "op (observable)")
+        exp_val = tr(_op * getRho(ados))
+    end
 
     if take_real
         return real(exp_val)
@@ -188,37 +206,26 @@ where ``O`` is the operator and ``\rho`` is the reduced density operator in one
 function Expect(op, ados_list::Vector{ADOs}; take_real=true)
 
     dim = ados_list[1].dim
+    N   = ados_list[1].N
     for i in 2:length(ados_list)
-        if ados_list[i].dim != dim
-            error("The dimension of the elements in `ados_list` should be consistent.")
-        end
+        _check_sys_dim_and_ADOs_num(ados_list[1], ados_list[i])
     end
 
-    _op = HandleMatrixType(op, dim, "op (observable)")
+    if typeof(op) == HEOMSuperOp
+        _check_sys_dim_and_ADOs_num(op, ados_list[1])
+        _op = op
+    else
+        _op = HEOMSuperOp(op, EVEN, dim, N, "L")
+    end
+    tr_op = Tr(dim, N) * _op.data
 
-    ρlist = getRho.(ados_list)
     exp_val = [
-        tr(_op * ρlist[i]) for i in 1:length(ρlist)
+        (tr_op * ados.data) for ados in ados_list
     ]
 
     if take_real
         return real.(exp_val)
     else
         return exp_val
     end
-end
-
-# for changing a `Vector` back to `ADOs`
-function _HandleVectorType(V::T, cp::Bool=true) where T <: Vector
-    if cp
-        return Vector{ComplexF64}(V)
-    else
-        return V
-    end
-end
-
-# for changing the type of `ADOs` to match the type of HEOMLS matrix 
-function _HandleVectorType(MatrixType::Type{TM}, V::SparseVector) where TM <: SparseMatrixCSC
-    TE = eltype(MatrixType)
-    return Vector{TE}(V)
 end

src/Bath.jl

@@ -30,9 +30,6 @@ struct Exponent
     types::String
 end
 
-spre(q::AbstractMatrix)  = sparse(kron(Matrix(I, size(q)[1], size(q)[1]), q))
-spost(q::AbstractMatrix) = sparse(kron(transpose(q), Matrix(I, size(q)[1], size(q)[1])))
-
 function show(io::IO, E::Exponent)
     print(io, 
         "Bath Exponent with types = \"$(E.types)\", operator size = $(size(E.op)), η = $(E.η), γ = $(E.γ).\n"

src/HeomBase.jl

@@ -1,5 +1,13 @@
+abstract type AbstractHEOMLSMatrix end
+
 const PROGBAR_OPTIONS = Dict(:barlen=>20, :color=>:green, :showspeed=>true)
 
+spre(q::AbstractMatrix)  = sparse(kron(Matrix(I, size(q)[1], size(q)[1]), q))
+spost(q::AbstractMatrix) = sparse(kron(transpose(q), Matrix(I, size(q)[1], size(q)[1])))
+
+# equal to : transpose(sparse(vec(system_identity_matrix)))
+Tr(dim::Int, N::Int)  = transpose(SparseVector(N * dim ^ 2, [1 + n * (dim + 1) for n in 0:(dim - 1)], ones(ComplexF64, dim)))
+
 function HandleMatrixType(M, dim::Int=0, MatrixName::String="")
     error("HierarchicalEOM doesn't support matrix $(MatrixName) with type : $(typeof(M))")
 end
@@ -16,7 +24,7 @@ function HandleMatrixType(M::AbstractMatrix, dim::Int=0, MatrixName::String="")
         if N1 == N2
             return copy(M)
         else
-            error("The size of matrix $(MatrixName) should be squared matrix.")
+            error("The matrix $(MatrixName) should be squared matrix.")
         end
     end
 end
@@ -38,6 +46,52 @@ function _HandleFloatType(ElType::Type{T}, V::Any) where T <: Number
     end
 end
 
+# for changing a `Vector` back to `ADOs`
+function _HandleVectorType(V::T, cp::Bool=true) where T <: Vector
+    if cp
+        return Vector{ComplexF64}(V)
+    else
+        return V
+    end
+end
+
+# for changing the type of `ADOs` to match the type of HEOMLS matrix 
+function _HandleVectorType(MatrixType::Type{TM}, V::SparseVector) where TM <: SparseMatrixCSC
+    TE = eltype(MatrixType)
+    return Vector{TE}(V)
+end
+
+function _HandleSteadyStateMatrix(MatrixType::Type{TM}, M::AbstractHEOMLSMatrix, S::Int) where TM <: SparseMatrixCSC
+    ElType = eltype(M)
+    A = copy(M.data)
+    A[1,1:S] .= 0
+
+    # sparse(row_idx, col_idx, values, row_dims, col_dims)
+    A += sparse(ones(ElType, M.dim), [(n - 1) * (M.dim + 1) + 1 for n in 1:(M.dim)], ones(ElType, M.dim), S, S)
+    return A
+end
+
+function _HandleIdentityType(MatrixType::Type{TM}, S::Int) where TM <: SparseMatrixCSC
+    ElType = eltype(MatrixType)
+    return sparse(one(ElType) * I, S, S)
+end
+
+function _check_sys_dim_and_ADOs_num(A, B)
+    if (A.dim != B.dim)
+        error("Inconsistent system dimension (\"dim\").")
+    end
+
+    if (A.N != B.N)
+        error("Inconsistent number of ADOs (\"N\").")
+    end
+end
+
+function _check_parity(A, B)
+    if typeof(A.parity) != typeof(B.parity)
+        error("Inconsistent parity.")
+    end
+end
+
 function _get_pkg_version(pkg_name::String)
     D = Pkg.dependencies()
     for uuid in keys(D)