Skip to content

Commit

Permalink
Merge pull request #130 from ytdHuang/sys-liouvillian
Browse files Browse the repository at this point in the history 
Support `SuperOperator` input for constructing HEOMLS
  • Loading branch information
@ytdHuang
ytdHuang authored Dec 16, 2024
2 parents ece438e + aad160a commit a62be5c
Show file tree
Hide file tree
Showing 10 changed files with 46 additions and 28 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "HierarchicalEOM"
uuid = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
authors = ["Yi-Te Huang <[email protected]>"]
version = "2.3.1"
version = "2.3.2"

[deps]
DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
Expand Down
32 changes: 24 additions & 8 deletions src/HeomBase.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,28 +39,44 @@ end
_Tr(M::AbstractHEOMLSMatrix) = _Tr(_get_SciML_matrix_wrapper(M), M.dims, M.N)
_Tr(M::Type{<:SparseMatrixCSC}, dims::SVector, N::Int) = _Tr(eltype(M), dims, N)

function HandleMatrixType(M::AbstractQuantumObject, MatrixName::String = ""; type::QuantumObjectType = Operator)
if M.type == type
return M
else
function HandleMatrixType(
M::AbstractQuantumObject,
MatrixName::String = "";
type::T = nothing,
) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}}
if (type isa Nothing)
(M.type isa OperatorQuantumObject) ||
(M.type isa SuperOperatorQuantumObject) ||
error("The matrix $(MatrixName) should be either an Operator or a SuperOperator.")
elseif M.type != type
error("The matrix $(MatrixName) should be an $(type).")
end
return M
end
function HandleMatrixType(
M::AbstractQuantumObject,
dims::SVector,
MatrixName::String = "";
type::QuantumObjectType = Operator,
)
type::T = nothing,
) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}}
if M.dims == dims
return HandleMatrixType(M, MatrixName; type = type)
else
error("The dims of $(MatrixName) should be: $(dims)")
end
end
HandleMatrixType(M, dims::SVector, MatrixName::String = ""; type::QuantumObjectType = Operator) =
HandleMatrixType(
M,
dims::SVector,
MatrixName::String = "";
type::T = nothing,
) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}} =
HandleMatrixType(M, MatrixName; type = type)
HandleMatrixType(M, MatrixName::String = ""; type::QuantumObjectType = Operator) =
HandleMatrixType(
M,
MatrixName::String = "";
type::T = nothing,
) where {T<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject}} =
error("HierarchicalEOM doesn't support matrix $(MatrixName) with type : $(typeof(M))")

# change the type of `ADOs` to match the type of HEOMLS matrix
Expand Down
10 changes: 5 additions & 5 deletions src/bath/BosonBath.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,14 @@ Base.show(io::IO, B::AbstractBosonBath) =
Base.show(io::IO, m::MIME"text/plain", B::AbstractBosonBath) = show(io, B)

function _check_bosonic_coupling_operator(op)
_op = HandleMatrixType(op, "op (coupling operator)")
_op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
if !ishermitian(_op)
@warn "The system-bosonic-bath coupling operator \"op\" should be Hermitian Operator."
end
return _op
end

_check_bosonic_RWA_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)")
_check_bosonic_RWA_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)", type = Operator)

function _combine_same_gamma(η::Vector{Ti}, γ::Vector{Tj}) where {Ti,Tj<:Number}
if length(η) != length(γ)
Expand Down Expand Up @@ -96,7 +96,7 @@ function BosonBath(
δ::Number = 0.0;
combine::Bool = true,
) where {Ti,Tj<:Number}
_op = HandleMatrixType(op, "op (coupling operator)")
_op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
if combine
ηnew, γnew = _combine_same_gamma(η, γ)
bRI = bosonRealImag(_op, real.(ηnew), imag.(ηnew), γnew)
Expand Down Expand Up @@ -144,7 +144,7 @@ function BosonBath(
δ::Tm = 0.0;
combine::Bool = true,
) where {Ti,Tj,Tk,Tl,Tm<:Number}
_op = HandleMatrixType(op, "op (coupling operator)")
_op = HandleMatrixType(op, "op (coupling operator)", type = Operator)

if combine
ηR, γR = _combine_same_gamma(η_real, γ_real)
Expand Down Expand Up @@ -366,7 +366,7 @@ function BosonBathRWA(
δ::Tm = 0.0,
) where {Ti,Tj,Tk,Tl,Tm<:Number}
_check_gamma_absorb_and_emit(γ_absorb, γ_emit)
_op = HandleMatrixType(op, "op (coupling operator)")
_op = HandleMatrixType(op, "op (coupling operator)", type = Operator)

bA = bosonAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
bE = bosonEmit(_op, η_emit, γ_emit, η_absorb)
Expand Down
4 changes: 2 additions & 2 deletions src/bath/FermionBath.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ Base.show(io::IO, B::AbstractFermionBath) =
print(io, "$(typeof(B))-type bath with $(B.Nterm) exponential-expansion terms\n")
Base.show(io::IO, m::MIME"text/plain", B::AbstractFermionBath) = show(io, B)

_check_fermionic_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)")
_check_fermionic_coupling_operator(op) = HandleMatrixType(op, "op (coupling operator)", type = Operator)

@doc raw"""
struct FermionBath <: AbstractBath
Expand Down Expand Up @@ -72,7 +72,7 @@ function FermionBath(
) where {Ti,Tj,Tk,Tl,Tm<:Number}
_check_gamma_absorb_and_emit(γ_absorb, γ_emit)

_op = HandleMatrixType(op, "op (coupling operator)")
_op = HandleMatrixType(op, "op (coupling operator)", type = Operator)
fA = fermionAbsorb(adjoint(_op), η_absorb, γ_absorb, η_emit)
fE = fermionEmit(_op, η_emit, γ_emit, η_absorb)
return FermionBath(AbstractFermionBath[fA, fE], _op, fA.Nterm + fE.Nterm, δ)
Expand Down
4 changes: 2 additions & 2 deletions src/heom_matrices/M_Boson.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ end
Generate the boson-type HEOM Liouvillian superoperator matrix
# Parameters
- `Hsys` : The time-independent system Hamiltonian
- `Hsys` : The time-independent system Hamiltonian or Liouvillian
- `tier::Int` : the tier (cutoff level) for the bosonic bath
- `Bath::Vector{BosonBath}` : objects for different bosonic baths
- `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
Expand All @@ -62,7 +62,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
)

# check for system dimension
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
sup_dim = prod(_Hsys.dims)^2
I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)

Expand Down
4 changes: 2 additions & 2 deletions src/heom_matrices/M_Boson_Fermion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ end
Generate the boson-fermion-type HEOM Liouvillian superoperator matrix
# Parameters
- `Hsys` : The time-independent system Hamiltonian
- `Hsys` : The time-independent system Hamiltonian or Liouvillian
- `Btier::Int` : the tier (cutoff level) for the bosonic bath
- `Ftier::Int` : the tier (cutoff level) for the fermionic bath
- `Bbath::Vector{BosonBath}` : objects for different bosonic baths
Expand All @@ -99,7 +99,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
)

# check for system dimension
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
sup_dim = prod(_Hsys.dims)^2
I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)

Expand Down
4 changes: 2 additions & 2 deletions src/heom_matrices/M_Fermion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ end
Generate the fermion-type HEOM Liouvillian superoperator matrix
# Parameters
- `Hsys` : The time-independent system Hamiltonian
- `Hsys` : The time-independent system Hamiltonian or Liouvillian
- `tier::Int` : the tier (cutoff level) for the fermionic bath
- `Bath::Vector{FermionBath}` : objects for different fermionic baths
- `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
Expand All @@ -60,7 +60,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix
)

# check for system dimension
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
sup_dim = prod(_Hsys.dims)^2
I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim)

Expand Down
4 changes: 2 additions & 2 deletions src/heom_matrices/M_S.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ M[\cdot]=-i \left[H_{sys}, \cdot \right]_-,
where ``[\cdot, \cdot]_-`` stands for commutator.
# Parameters
- `Hsys` : The time-independent system Hamiltonian
- `Hsys` : The time-independent system Hamiltonian or Liouvillian
- `parity::AbstractParity` : the parity label of the operator which HEOMLS is acting on (usually `EVEN`, only set as `ODD` for calculating spectrum of fermionic system).
- `verbose::Bool` : To display verbose output during the process or not. Defaults to `true`.
Expand All @@ -45,7 +45,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi
@noinline function M_S(Hsys::QuantumObject, parity::AbstractParity = EVEN; verbose::Bool = true)

# check for system dimension
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian)")
_Hsys = HandleMatrixType(Hsys, "Hsys (system Hamiltonian or Liouvillian)")
sup_dim = prod(_Hsys.dims)^2

# the Liouvillian operator for free Hamiltonian
Expand Down
2 changes: 1 addition & 1 deletion src/heom_matrices/heom_matrix_base.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -442,7 +442,7 @@ function bath_sum_γ(nvec, baths::Vector{T}) where {T<:Union{AbstractBosonBath,A
end

# commutator of system Hamiltonian
minus_i_L_op(Hsys::QuantumObject, Id = I(size(Hsys, 1))) = _liouvillian(Hsys.data, Id)
minus_i_L_op(Hsys::QuantumObject, Id = I(size(Hsys, 1))) = liouvillian(Hsys, Id).data

# connect to bosonic (n-1)th-level for "Real & Imag combined operator"
minus_i_D_op(bath::bosonRealImag, k, n_k) = n_k * (-1.0im * bath.η_real[k] * bath.Comm + bath.η_imag[k] * bath.anComm)
Expand Down
8 changes: 5 additions & 3 deletions test/M_S.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,13 @@
t = 10
Hsys = sigmax()
L = M_S(Hsys; verbose = false)
L_super = M_S(liouvillian(Hsys); verbose = false) # test if input is already Liouvillian
ψ0 = basis(2, 0)
@test show(devnull, MIME("text/plain"), L) === nothing
@test size(L) == (4, 4)
@test L.N == 1
@test nnz(L.data.A) == nnz(L(0)) == 8
@test size(L) == size(L_super) == (4, 4)
@test L.N == L_super.N == 1
@test nnz(L.data.A) == nnz(L_super.data.A) == nnz(L(0)) == nnz(L_super(0)) == 8
@test L.data == L_super.data
ados_list = HEOMsolve(L, ψ0, 0:1:t; reltol = 1e-8, abstol = 1e-10, verbose = false).ados
ados = ados_list[end]
@test ados.dims == L.dims
Expand Down

4 comments on commit a62be5c

@ytdHuang
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator registry

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Error while trying to register: Action not recognized: registry

@ytdHuang
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/121458

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v2.3.2 -m "<description of version>" a62be5c53bd95f1923dcbb98a9e35a86e8871023
git push origin v2.3.2

Please sign in to comment.