diff --git a/.github/workflows/Runtests.yml b/.github/workflows/Runtests.yml index d7041196..a1758844 100644 --- a/.github/workflows/Runtests.yml +++ b/.github/workflows/Runtests.yml @@ -7,6 +7,9 @@ on: pull_request: branches: - 'main' + types: + - 'synchronize' + - 'ready_for_review' jobs: test: diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index bf181bd8..88caf936 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -9,6 +9,9 @@ on: pull_request: branches: - 'main' + types: + - 'synchronize' + - 'ready_for_review' jobs: build: diff --git a/docs/make.jl b/docs/make.jl index cf2d1570..e0216717 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/Parity.md b/docs/src/Parity.md new file mode 100644 index 00000000..c36a4f79 --- /dev/null +++ b/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 +``` \ No newline at end of file diff --git a/docs/src/heom_matrix/HEOMLS_intro.md b/docs/src/heom_matrix/HEOMLS_intro.md index e9420f7a..9b879300 100644 --- a/docs/src/heom_matrix/HEOMLS_intro.md +++ b/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. diff --git a/docs/src/libraryAPI.md b/docs/src/libraryAPI.md index 8f324be5..510e5cce 100644 --- a/docs/src/libraryAPI.md +++ b/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 diff --git a/docs/src/spectrum.md b/docs/src/spectrum.md index b92c879a..8a095a93 100644 --- a/docs/src/spectrum.md +++ b/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 diff --git a/ext/HierarchicalEOM_CUDAExt.jl b/ext/HierarchicalEOM_CUDAExt.jl index f428a247..fd2f5c73 100644 --- a/ext/HierarchicalEOM_CUDAExt.jl +++ b/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 diff --git a/src/ADOs.jl b/src/ADOs.jl index 64a8cbaa..7d835e32 100644 --- a/src/ADOs.jl +++ b/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,17 +206,21 @@ 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 @@ -206,19 +228,4 @@ function Expect(op, ados_list::Vector{ADOs}; take_real=true) 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 \ No newline at end of file diff --git a/src/Bath.jl b/src/Bath.jl index 99f5419b..db5c8a1a 100644 --- a/src/Bath.jl +++ b/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" diff --git a/src/HeomBase.jl b/src/HeomBase.jl index ac36a3e6..c5eae7e8 100644 --- a/src/HeomBase.jl +++ b/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) diff --git a/src/HierarchicalEOM.jl b/src/HierarchicalEOM.jl index d12bc46f..62996b5e 100644 --- a/src/HierarchicalEOM.jl +++ b/src/HierarchicalEOM.jl @@ -7,26 +7,34 @@ module HierarchicalEOM # sub-module HeomBase for HierarchicalEOM module HeomBase import Pkg - import LinearAlgebra: BLAS + import LinearAlgebra: BLAS, kron + import SparseArrays: I, sparse, SparseVector, SparseMatrixCSC import Crayons: Crayon + export + spre, spost, Tr, + PROGBAR_OPTIONS, + AbstractHEOMLSMatrix, + _check_sys_dim_and_ADOs_num, _check_parity, + HandleMatrixType, _HandleFloatType, _HandleVectorType, _HandleSteadyStateMatrix, _HandleIdentityType + include("HeomBase.jl") end import .HeomBase.versioninfo as versioninfo import .HeomBase.print_logo as print_logo + @reexport using .HeomBase: AbstractHEOMLSMatrix, spre, spost # sub-module Bath for HierarchicalEOM module Bath + using ..HeomBase import Base: show, length, getindex, lastindex, iterate, checkbounds - import LinearAlgebra: I, kron, ishermitian, eigvals - import SparseArrays: sparse, SparseMatrixCSC - import ..HeomBase: HandleMatrixType - + import LinearAlgebra: ishermitian, eigvals + import SparseArrays: SparseMatrixCSC + export AbstractBath, BosonBath, BosonBathRWA, FermionBath, Exponent, C, AbstractBosonBath, bosonReal, bosonImag, bosonRealImag, bosonAbsorb, bosonEmit, AbstractFermionBath, fermionAbsorb, fermionEmit, - spre, spost, Boson_DrudeLorentz_Matsubara, Boson_DrudeLorentz_Pade, Fermion_Lorentz_Matsubara, Fermion_Lorentz_Pade @@ -37,14 +45,14 @@ module HierarchicalEOM # sub-module HeomAPI for HierarchicalEOM module HeomAPI - using ..Bath - import Base: ==, show, length, size, getindex, keys, setindex!, lastindex, iterate, checkbounds, hash, copy, eltype + using ..HeomBase + using ..Bath + 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 + import SparseArrays: sparse, sparsevec, spzeros, SparseVector, SparseMatrixCSC import ProgressMeter: Progress, next! import FastExpm: fastExpm - import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType, _HandleFloatType # for solving time evolution import SciMLOperators: MatrixOperator @@ -57,29 +65,40 @@ module HierarchicalEOM import SteadyStateDiffEq: DynamicSS export - AbstractHEOMLSMatrix, M_S, M_Boson, M_Fermion, M_Boson_Fermion, AbstractParity, OddParity, EvenParity, value, ODD, EVEN, ADOs, getRho, getADO, Expect, Nvec, AbstractHierarchyDict, HierarchyDict, MixHierarchyDict, getIndexEnsemble, + HEOMSuperOp, M_S, M_Boson, M_Fermion, M_Boson_Fermion, Propagator, addBosonDissipator, addFermionDissipator, addTerminator, evolution, SteadyState - include("heom_api.jl") + include("Parity.jl") + include("ADOs.jl") + + include("heom_matrices/heom_matrix_base.jl") + include("heom_matrices/Nvec.jl") + include("heom_matrices/HierarchyDict.jl") + include("heom_matrices/M_S.jl") + include("heom_matrices/M_Boson.jl") + include("heom_matrices/M_Fermion.jl") + include("heom_matrices/M_Boson_Fermion.jl") + + include("evolution.jl") + include("SteadyState.jl") end @reexport using .HeomAPI # sub-module Spectrum for HierarchicalEOM module Spectrum - import ..HeomAPI: AbstractHEOMLSMatrix, OddParity, ADOs, spre, _HandleVectorType - import LinearAlgebra: I, kron - import SparseArrays: sparse, sparsevec, SparseMatrixCSC + using ..HeomBase + import ..HeomAPI: HEOMSuperOp, ADOs, EVEN, ODD import LinearSolve: LinearProblem, init, solve!, UMFPACKFactorization import ProgressMeter: Progress, next! - import ..HeomBase: PROGBAR_OPTIONS, HandleMatrixType, _HandleFloatType export spectrum, PowerSpectrum, DensityOfStates - include("Spectrum.jl") + include("power_spectrum.jl") + include("density_of_states.jl") end @reexport using .Spectrum diff --git a/src/Parity.jl b/src/Parity.jl new file mode 100644 index 00000000..ac58663e --- /dev/null +++ b/src/Parity.jl @@ -0,0 +1,39 @@ +abstract type AbstractParity end + +@doc raw""" + struct OddParity <: AbstractParity +""" +struct OddParity <: AbstractParity end + +@doc raw""" + struct EvenParity <: AbstractParity +""" +struct EvenParity <: AbstractParity end + +value(p::OddParity) = 1 +value(p::EvenParity) = 0 + +!(p::OddParity) = EVEN +!(p::EvenParity) = ODD +*(p1::TP1, p2::TP2) where {TP1, TP2 <: AbstractParity} = TP1 == TP2 ? EVEN : ODD + +function show(io::IO, p::OddParity) + print(io, "odd-parity") +end + +function show(io::IO, p::EvenParity) + print(io, "even-parity") +end + +# Parity label +@doc raw""" + const ODD = OddParity() +Label of odd-parity +""" +const ODD = OddParity() + +@doc raw""" + const EVEN = EvenParity() +Label of even-parity +""" +const EVEN = EvenParity() \ No newline at end of file diff --git a/src/Spectrum.jl b/src/Spectrum.jl deleted file mode 100644 index f9f8119d..00000000 --- a/src/Spectrum.jl +++ /dev/null @@ -1,320 +0,0 @@ -@doc raw""" - spectrum(M, ρ, op, ωlist; solver, verbose, filename, SOLVEROptions...) -!!! warning "Warning" - This function has been deprecated start from `HierarchicalEOM v1.1`, use `PowerSpectrum` or `DensityOfStates` instead. -""" -function spectrum( - M::AbstractHEOMLSMatrix, - ρ, - op, - ωlist::AbstractVector; - solver=UMFPACKFactorization(), - verbose::Bool = true, - filename::String = "", - SOLVEROptions... - ) - error("This function has been deprecated start from \`HierarchicalEOM v1.1\`, use \`PowerSpectrum\` or \`DensityOfStates\` instead.") -end - -@doc raw""" - PowerSpectrum(M, ρ, Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...) -Calculate power spectrum for the system where `P_op` will be automatically set as the adjoint of `Q_op`. - -This function is equivalent to: -`PowerSpectrum(M, ρ, Q_op', Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...)` -""" -function PowerSpectrum( - M::AbstractHEOMLSMatrix, - ρ, - Q_op, - ωlist::AbstractVector, - reverse::Bool = false; - solver=UMFPACKFactorization(), - verbose::Bool = true, - filename::String = "", - SOLVEROptions... - ) - 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...) -Calculate power spectrum for the system. - -```math -\pi S(\omega)=\textrm{Re}\left\{\int_0^\infty dt \langle P(t) Q(0)\rangle e^{-i\omega t}\right\}, -``` - -# To calculate spectrum when input operator `Q_op` has `EVEN`-parity: -remember to set the parameters: -- `M::AbstractHEOMLSMatrix`: should be `EVEN` parity - -# To calculate spectrum when input operator `Q_op` has `ODD`-parity: -remember to set the parameters: -- `M::AbstractHEOMLSMatrix`: should be `ODD` parity - -# Parameters -- `M::AbstractHEOMLSMatrix` : the HEOMLS matrix. -- `ρ` : the system density matrix or the auxiliary density operators. -- `P_op`: the operator ``P`` acting on the system. -- `Q_op`: the operator ``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` : 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`. -- `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process. -- `SOLVEROptions` : extra options for solver - -For more details about solvers and extra options, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/) - -# Returns -- `spec::AbstractVector` : the spectrum list corresponds to the specified `ωlist` -""" -function PowerSpectrum( - M::AbstractHEOMLSMatrix, - ρ, - P_op, - Q_op, - ωlist::AbstractVector, - reverse::Bool = false; - solver = UMFPACKFactorization(), - verbose::Bool = true, - filename::String = "", - SOLVEROptions... - ) - - Size = size(M, 1) - - # check ρ - if typeof(ρ) == ADOs # ρ::ADOs - if (M.dim != ρ.dim) - error("The system dimension between M and ρ are not consistent.") - end - if (M.N != ρ.N) - error("The ADOs number \"N\" between M and ados are not consistent.") - end - if (typeof(ρ.parity) == OddParity) - error("The parity of ρ or the ADOs must be `EVEN`.") - end - - ados_vec = ρ.data - - else - _ρ = HandleMatrixType(ρ, M.dim, "ρ (state)") - v = sparsevec(_ρ) - ados_vec = sparsevec(v.nzind, v.nzval, Size) - end - - # check dimension of P_op and Q_op - _P = HandleMatrixType(P_op, M.dim, "P_op (operator)") - _Q = HandleMatrixType(Q_op, M.dim, "Q_op (operator)") - - SAVE::Bool = (filename != "") - if SAVE - FILENAME = filename * ".txt" - if isfile(FILENAME) - error("FILE: $(FILENAME) already exist.") - end - end - - I_heom = sparse(one(ComplexF64) * I, M.N, M.N) - I_total = _HandleIdentityType(typeof(M.data), Size) - - # equal to : transpose(sparse(vec(system_identity_matrix))) - I_dual_vec = transpose(sparsevec([1 + n * (M.dim + 1) for n in 0:(M.dim - 1)], ones(ComplexF64, M.dim), M.sup_dim)) - - # operator for calculating two-time correlation functions in frequency domain - P_sup = kron(I_heom, spre(_P)) - Q_sup = kron(I_heom, spre(_Q)) - X = _HandleVectorType(typeof(M.data), Q_sup * ados_vec) - - ElType = eltype(M) - ωList = _HandleFloatType(ElType, ωlist) - Length = length(ωList) - Sω = Vector{Float64}(undef, Length) - - if verbose - print("Calculating power spectrum...\n") - flush(stdout) - prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) - end - if reverse - i = convert(ElType, 1im) - else - i = convert(ElType, -1im) - end - Iω = i * ωList[1] * I_total - cache = init(LinearProblem(M.data + Iω, X), solver, SOLVEROptions...) - sol = solve!(cache) - @inbounds for (j, ω) in enumerate(ωList) - if j > 1 - Iω = i * ω * I_total - cache.A = M.data + Iω - sol = solve!(cache) - end - - # trace over the Hilbert space of system (expectation value) - Sω[j] = -1 * real(I_dual_vec * (P_sup * _HandleVectorType(sol.u, false))[1:(M.sup_dim)]) - - if SAVE - open(FILENAME, "a") do file - write(file, "$(Sω[j]),\n") - end - end - - if verbose - next!(prog) - end - end - if verbose - println("[DONE]") - end - - return Sω -end - -@doc raw""" - DensityOfStates(M, ρ, P_op, Q_op, ωlist; solver, verbose, filename, SOLVEROptions...) -Calculate density of states for the fermionic system. - -```math - \pi A(\omega)=\textrm{Re}\left\{\int_0^\infty dt \left[\langle d(t) d^\dagger(0)\rangle^* + \langle d^\dagger(t) d(0)\rangle \right] e^{-i\omega t}\right\}, -``` - -# Parameters -- `M::AbstractHEOMLSMatrix` : the HEOMLS matrix which acts on `ODD`-parity operators. -- `ρ` : the system density matrix or the auxiliary density operators. -- `op` : The annihilation operator (``d`` as shown above) acting on the fermionic system. -- `ωlist::AbstractVector` : the specific frequency points to solve. -- `solver` : 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`. -- `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process. -- `SOLVEROptions` : extra options for solver - -For more details about solvers and extra options, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/) - -# Returns -- `dos::AbstractVector` : the list of density of states corresponds to the specified `ωlist` -""" -@noinline function DensityOfStates( - M::AbstractHEOMLSMatrix, - ρ, - op, - ωlist::AbstractVector; - solver=UMFPACKFactorization(), - verbose::Bool = true, - filename::String = "", - SOLVEROptions... - ) - - Size = size(M, 1) - - # check M - if (typeof(M.parity) != OddParity) - error("The HEOMLS matrix M must be acting on `ODD`-parity operators.") - end - - # check ρ - if typeof(ρ) == ADOs # ρ::ADOs - if (M.dim != ρ.dim) - error("The system dimension between M and ρ are not consistent.") - end - if (M.N != ρ.N) - error("The ADOs number \"N\" between M and ados are not consistent.") - end - if (typeof(ρ.parity) == OddParity) - error("The parity of ρ or the ADOs must be `EVEN`.") - end - - ados_vec = ρ.data - - else - _ρ = HandleMatrixType(ρ, M.dim, "ρ (state)") - v = sparsevec(_ρ) - ados_vec = sparsevec(v.nzind, v.nzval, Size) - end - - # check dimension of op - _op = HandleMatrixType(op, M.dim, "op (operator)") - - SAVE::Bool = (filename != "") - if SAVE - FILENAME = filename * ".txt" - if isfile(FILENAME) - error("FILE: $(FILENAME) already exist.") - end - end - - I_heom = sparse(one(ComplexF64) * I, M.N, M.N) - I_total = _HandleIdentityType(typeof(M.data), Size) - - # equal to : transpose(sparse(vec(system_identity_matrix))) - I_dual_vec = transpose(sparsevec([1 + n * (M.dim + 1) for n in 0:(M.dim - 1)], ones(ComplexF64, M.dim), M.sup_dim)) - - # operators for calculating two-time correlation functions in frequency domain - d_normal = kron(I_heom, spre(_op)) - d_dagger = kron(I_heom, spre(_op')) - X_m = _HandleVectorType(typeof(M.data), d_normal * ados_vec) - X_p = _HandleVectorType(typeof(M.data), d_dagger * ados_vec) - - ElType = eltype(M) - ωList = _HandleFloatType(ElType, ωlist) - Length = length(ωList) - Aω = Vector{Float64}(undef, Length) - - if verbose - print("Calculating density of states...\n") - flush(stdout) - prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) - end - i = convert(ElType, 1im) - Iω = i * ωList[1] * I_total - cache_m = init(LinearProblem(M.data - Iω, X_m), solver, SOLVEROptions...) - cache_p = init(LinearProblem(M.data + Iω, X_p), solver, SOLVEROptions...) - sol_m = solve!(cache_m) - sol_p = solve!(cache_p) - @inbounds for (j, ω) in enumerate(ωList) - if j > 1 - Iω = i * ω * I_total - - cache_m.A = M.data - Iω - sol_m = solve!(cache_m) - - cache_p.A = M.data + Iω - sol_p = solve!(cache_p) - end - Cω_m = d_dagger * _HandleVectorType(sol_m.u, false) - Cω_p = d_normal * _HandleVectorType(sol_p.u, false) - - # trace over the Hilbert space of system (expectation value) - Aω[j] = -1 * ( - real(I_dual_vec * Cω_p[1:(M.sup_dim)]) + - real(I_dual_vec * Cω_m[1:(M.sup_dim)]) - ) - - if SAVE - open(FILENAME, "a") do file - write(file, "$(Aω[j]),\n") - end - end - - if verbose - next!(prog) - end - end - if verbose - println("[DONE]") - end - - return Aω -end - -function _HandleIdentityType(MatrixType::Type{TM}, S::Int) where TM <: SparseMatrixCSC - ElType = eltype(MatrixType) - return sparse(one(ElType) * I, S, S) -end \ No newline at end of file diff --git a/src/SteadyState.jl b/src/SteadyState.jl index 78450219..ce0a4150 100644 --- a/src/SteadyState.jl +++ b/src/SteadyState.jl @@ -70,14 +70,9 @@ function SteadyState( verbose::Bool = true, SOLVEROptions... ) - - _ρ0 = HandleMatrixType(ρ0, M.dim, "ρ0 (initial state)") - - # vectorize initial state - ρ1 = sparse(sparsevec(_ρ0)) - ados = ADOs(sparsevec(ρ1.nzind, ρ1.nzval, M.N * M.sup_dim), M.N, M.parity) - - return SteadyState(M, ados; + return SteadyState( + M, + ADOs(ρ0, M.N, M.parity); solver = solver, reltol = reltol, abstol = abstol, @@ -121,22 +116,8 @@ For more details about solvers and extra options, please refer to [`Differential SOLVEROptions... ) - # check parity - if typeof(M.parity) == OddParity - error("The parity of M should be \"EVEN\".") - end - - if (M.dim != ados.dim) - error("The system dimension between M and ados are not consistent.") - end - - if (M.N != ados.N) - error("The ADOs number \"N\" between M and ados are not consistent.") - end - - if (typeof(M.parity) != typeof(ados.parity)) - error("The parity between M and ados are not consistent.") - end + _check_sys_dim_and_ADOs_num(M, ados) + _check_parity(M, ados) # problem: dρ(t)/dt = L * ρ(t) L = MatrixOperator(M.data) @@ -163,14 +144,4 @@ For more details about solvers and extra options, please refer to [`Differential end return ADOs(_HandleVectorType(sol.u, false), M.dim, M.N, M.parity) -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 \ No newline at end of file diff --git a/src/density_of_states.jl b/src/density_of_states.jl new file mode 100644 index 00000000..4c9bff95 --- /dev/null +++ b/src/density_of_states.jl @@ -0,0 +1,119 @@ +@doc raw""" + DensityOfStates(M, ρ, d_op, ωlist; solver, verbose, filename, SOLVEROptions...) +Calculate density of states for the fermionic system in frequency domain. + +```math + \pi A(\omega)=\textrm{Re}\left\{\int_0^\infty dt \left[\langle d(t) d^\dagger(0)\rangle^* + \langle d^\dagger(t) d(0)\rangle \right] e^{-i\omega t}\right\}, +``` + +# 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. +- `ωlist::AbstractVector` : the specific frequency points to solve. +- `solver` : 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`. +- `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process. +- `SOLVEROptions` : extra options for solver + +For more details about solvers and extra options, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/) + +# Returns +- `dos::AbstractVector` : the list of density of states corresponds to the specified `ωlist` +""" +@noinline function DensityOfStates( + M::AbstractHEOMLSMatrix, + ρ, + d_op, + ωlist::AbstractVector; + solver=UMFPACKFactorization(), + verbose::Bool = true, + filename::String = "", + SOLVEROptions... + ) + + Size = size(M, 1) + + # check M + if M.parity == EVEN + error("The HEOMLS matrix M must be acting on `ODD`-parity operators.") + end + + # Handle ρ + if typeof(ρ) == ADOs # ρ::ADOs + ados = ρ + if ados.parity != EVEN + error("The parity of ρ must be `EVEN`.") + end + else + ados = ADOs(ρ, M.N) + end + _check_sys_dim_and_ADOs_num(M, ados) + + # Handle d_op + _tr = Tr(M.dim, M.N) + d_normal = HEOMSuperOp(d_op, ODD, M) + d_dagger = HEOMSuperOp(d_op', ODD, M) + b_m = _HandleVectorType(typeof(M.data), (d_normal * ados).data) + b_p = _HandleVectorType(typeof(M.data), (d_dagger * ados).data) + _tr_d_normal = _tr * d_normal.data + _tr_d_dagger = _tr * d_dagger.data + + SAVE::Bool = (filename != "") + if SAVE + FILENAME = filename * ".txt" + if isfile(FILENAME) + error("FILE: $(FILENAME) already exist.") + end + end + + ElType = eltype(M) + ωList = _HandleFloatType(ElType, ωlist) + Length = length(ωList) + Aω = Vector{Float64}(undef, Length) + + if verbose + print("Calculating density of states in frequency domain...\n") + flush(stdout) + prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) + end + i = convert(ElType, 1im) + I_total = _HandleIdentityType(typeof(M.data), Size) + Iω = i * ωList[1] * I_total + cache_m = init(LinearProblem(M.data - Iω, b_m), solver, SOLVEROptions...) + cache_p = init(LinearProblem(M.data + Iω, b_p), solver, SOLVEROptions...) + sol_m = solve!(cache_m) + sol_p = solve!(cache_p) + @inbounds for (j, ω) in enumerate(ωList) + if j > 1 + Iω = i * ω * I_total + + cache_m.A = M.data - Iω + sol_m = solve!(cache_m) + + cache_p.A = M.data + Iω + sol_p = solve!(cache_p) + end + + # trace over the Hilbert space of system (expectation value) + Aω[j] = -1 * ( + real(_tr_d_normal * _HandleVectorType(sol_p.u, false)) + + real(_tr_d_dagger * _HandleVectorType(sol_m.u, false)) + ) + + if SAVE + open(FILENAME, "a") do file + write(file, "$(Aω[j]),\n") + end + end + + if verbose + next!(prog) + end + end + if verbose + println("[DONE]") + end + + return Aω +end \ No newline at end of file diff --git a/src/evolution.jl b/src/evolution.jl index 95ac99c2..ac014eb7 100644 --- a/src/evolution.jl +++ b/src/evolution.jl @@ -30,18 +30,15 @@ function evolution( verbose::Bool = true, filename::String = "" ) - - _ρ0 = HandleMatrixType(ρ0, M.dim, "ρ0 (initial state)") - - # vectorize initial state - ρ1 = sparse(sparsevec(_ρ0)) - ados = ADOs(sparsevec(ρ1.nzind, ρ1.nzval, M.N * M.sup_dim), M.N, M.parity) - - return evolution(M, ados, Δt, steps; + return evolution( + M, + ADOs(ρ0, M.N, M.parity), + Δt, + steps; threshold = threshold, nonzero_tol = nonzero_tol, - verbose = verbose, - filename = filename + verbose = verbose, + filename = filename ) end @@ -78,17 +75,8 @@ For more details, please refer to [`FastExpm.jl`](https://github.com/fmentink/Fa filename::String = "" ) - if (M.dim != ados.dim) - error("The system dimension between M and ados are not consistent.") - end - - if (M.N != ados.N) - error("The number N between M and ados are not consistent.") - end - - if (typeof(M.parity) != typeof(ados.parity)) - error("The parity between M and ados are not consistent.") - end + _check_sys_dim_and_ADOs_num(M, ados) + _check_parity(M, ados) SAVE::Bool = (filename != "") if SAVE @@ -183,20 +171,16 @@ function evolution( filename::String = "", SOLVEROptions... ) - - _ρ0 = HandleMatrixType(ρ0, M.dim, "ρ0 (initial state)") - - # vectorize initial state - ρ1 = sparse(sparsevec(_ρ0)) - ados = ADOs(sparsevec(ρ1.nzind, ρ1.nzval, M.N * M.sup_dim), M.N, M.parity) - - return evolution(M, ados, tlist; + return evolution( + M, + ADOs(ρ0, M.N, M.parity), + tlist; solver = solver, reltol = reltol, abstol = abstol, maxiters = maxiters, save_everystep = save_everystep, - verbose = verbose, + verbose = verbose, filename = filename, SOLVEROptions... ) @@ -239,17 +223,8 @@ For more details about solvers and extra options, please refer to [`Differential SOLVEROptions... ) - if (M.dim != ados.dim) - error("The system dimension between M and ados are not consistent.") - end - - if (M.N != ados.N) - error("The ADOs number \"N\" between M and ados are not consistent.") - end - - if (typeof(M.parity) != typeof(ados.parity)) - error("The parity between M and ados are not consistent.") - end + _check_sys_dim_and_ADOs_num(M, ados) + _check_parity(M, ados) SAVE::Bool = (filename != "") if SAVE @@ -354,20 +329,18 @@ function evolution( filename::String = "", SOLVEROptions... ) - - _ρ0 = HandleMatrixType(ρ0, M.dim, "ρ0 (initial state)") - - # vectorize initial state - ρ1 = sparse(sparsevec(_ρ0)) - ados = ADOs(sparsevec(ρ1.nzind, ρ1.nzval, M.N * M.sup_dim), M.N, M.parity) - - return evolution(M, ados, tlist, H, param; + return evolution( + M, + ADOs(ρ0, M.N, M.parity), + tlist, + H, + param; solver = solver, reltol = reltol, abstol = abstol, maxiters = maxiters, save_everystep = save_everystep, - verbose = verbose, + verbose = verbose, filename = filename, SOLVEROptions... ) @@ -413,17 +386,8 @@ For more details about solvers and extra options, please refer to [`Differential SOLVEROptions... ) - if (M.dim != ados.dim) - error("The system dimension between M and ados are not consistent.") - end - - if (M.N != ados.N) - error("The ADOs number \"N\" between M and ados are not consistent.") - end - - if (typeof(M.parity) != typeof(ados.parity)) - error("The parity between M and ados are not consistent.") - end + _check_sys_dim_and_ADOs_num(M, ados) + _check_parity(M, ados) SAVE::Bool = (filename != "") if SAVE @@ -443,8 +407,8 @@ For more details about solvers and extra options, please refer to [`Differential Ht = H(param, Tlist[1]) _Ht = HandleMatrixType(Ht, M.dim, "H (Hamiltonian) at t=$(Tlist[1])") - Lt = kron(sparse(I, M.N, M.N), - 1im * (spre(_Ht) - spost(_Ht))) - L = MatrixOperator(M.data + Lt, update_func! = _update_L!) + Lt = HEOMSuperOp(minus_i_L_op(_Ht), M.parity, M, "LR") + L = MatrixOperator((M + Lt).data, update_func! = _update_L!) # problem: dρ/dt = L(t) * ρ(0) ## M.dim will check whether the returned time-dependent Hamiltonian has the correct dimension @@ -503,6 +467,6 @@ function _update_L!(L, u, p, t) _Ht = HandleMatrixType(Ht, M.dim, "H (Hamiltonian) at t=$(t)") # update the block diagonal terms of L - L .= M.data - kron(sparse(I, M.N, M.N), 1im * (spre(_Ht) - spost(_Ht))) + L .= (M + HEOMSuperOp(minus_i_L_op(_Ht), M.parity, M, "LR")).data nothing end \ No newline at end of file diff --git a/src/heom_api.jl b/src/heom_api.jl deleted file mode 100644 index dbb397b2..00000000 --- a/src/heom_api.jl +++ /dev/null @@ -1,11 +0,0 @@ -include("heom_matrices/Nvec.jl") -include("heom_matrices/HierarchyDict.jl") -include("heom_matrices/heom_matrix_base.jl") -include("heom_matrices/M_S.jl") -include("heom_matrices/M_Boson.jl") -include("heom_matrices/M_Fermion.jl") -include("heom_matrices/M_Boson_Fermion.jl") - -include("ADOs.jl") -include("evolution.jl") -include("SteadyState.jl") \ No newline at end of file diff --git a/src/heom_matrices/M_Boson.jl b/src/heom_matrices/M_Boson.jl index bbef6141..3433ee1c 100644 --- a/src/heom_matrices/M_Boson.jl +++ b/src/heom_matrices/M_Boson.jl @@ -59,7 +59,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim) # the Liouvillian operator for free Hamiltonian term - Lsys = -1im * (spre(_Hsys) - spost(_Hsys)) + Lsys = minus_i_L_op(_Hsys) # bosonic bath if verbose && (threshold > 0.0) @@ -111,7 +111,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_minus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, nvec_neigh) idx_neigh = nvec2idx[nvec_neigh] - op = _D_op(bB, k, n_k) + op = minus_i_D_op(bB, k, n_k) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_plus!(nvec_neigh, mode) @@ -122,7 +122,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_plus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, nvec_neigh) idx_neigh = nvec2idx[nvec_neigh] - op = _B_op(bB) + op = minus_i_B_op(bB) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_minus!(nvec_neigh, mode) diff --git a/src/heom_matrices/M_Boson_Fermion.jl b/src/heom_matrices/M_Boson_Fermion.jl index 5232d1d9..5ca02a49 100644 --- a/src/heom_matrices/M_Boson_Fermion.jl +++ b/src/heom_matrices/M_Boson_Fermion.jl @@ -76,7 +76,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim) # the Liouvillian operator for free Hamiltonian term - Lsys = -1im * (spre(_Hsys) - spost(_Hsys)) + Lsys = minus_i_L_op(_Hsys) # check for bosonic and fermionic bath if verbose && (threshold > 0.0) @@ -129,7 +129,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_minus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, (nvec_neigh, nvec_f)) idx_neigh = nvec2idx[(nvec_neigh, nvec_f)] - op = _D_op(bB, k, n_k) + op = minus_i_D_op(bB, k, n_k) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_plus!(nvec_neigh, mode) @@ -140,7 +140,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_plus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, (nvec_neigh, nvec_f)) idx_neigh = nvec2idx[(nvec_neigh, nvec_f)] - op = _B_op(bB) + op = minus_i_B_op(bB) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_minus!(nvec_neigh, mode) @@ -161,7 +161,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_minus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, (nvec_b, nvec_neigh)) idx_neigh = nvec2idx[(nvec_b, nvec_neigh)] - op = _C_op(fB, k, nvec_f.level, sum(nvec_neigh[1:(mode - 1)]), parity) + op = minus_i_C_op(fB, k, nvec_f.level, sum(nvec_neigh[1:(mode - 1)]), parity) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_plus!(nvec_neigh, mode) @@ -171,7 +171,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi Nvec_plus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, (nvec_b, nvec_neigh)) idx_neigh = nvec2idx[(nvec_b, nvec_neigh)] - op = _A_op(fB, nvec_f.level, sum(nvec_neigh[1:(mode - 1)]), parity) + op = minus_i_A_op(fB, nvec_f.level, sum(nvec_neigh[1:(mode - 1)]), parity) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_minus!(nvec_neigh, mode) diff --git a/src/heom_matrices/M_Fermion.jl b/src/heom_matrices/M_Fermion.jl index a70246dc..e76460ae 100644 --- a/src/heom_matrices/M_Fermion.jl +++ b/src/heom_matrices/M_Fermion.jl @@ -57,7 +57,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix I_sup = sparse(one(ComplexF64) * I, sup_dim, sup_dim) # the Liouvillian operator for free Hamiltonian term - Lsys = -1im * (spre(_Hsys) - spost(_Hsys)) + Lsys = minus_i_L_op(_Hsys) # fermionic bath if verbose && (threshold > 0.0) @@ -109,7 +109,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix Nvec_minus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, nvec_neigh) idx_neigh = nvec2idx[nvec_neigh] - op = _C_op(fB, k, nvec.level, sum(nvec_neigh[1:(mode - 1)]), parity) + op = minus_i_C_op(fB, k, nvec.level, sum(nvec_neigh[1:(mode - 1)]), parity) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_plus!(nvec_neigh, mode) @@ -119,7 +119,7 @@ Generate the fermion-type HEOM Liouvillian superoperator matrix Nvec_plus!(nvec_neigh, mode) if (threshold == 0.0) || haskey(nvec2idx, nvec_neigh) idx_neigh = nvec2idx[nvec_neigh] - op = _A_op(fB, nvec.level, sum(nvec_neigh[1:(mode - 1)]), parity) + op = minus_i_A_op(fB, nvec.level, sum(nvec_neigh[1:(mode - 1)]), parity) add_operator!(op, L_row[tID], L_col[tID], L_val[tID], Nado, idx, idx_neigh) end Nvec_minus!(nvec_neigh, mode) diff --git a/src/heom_matrices/M_S.jl b/src/heom_matrices/M_S.jl index 2be5d178..a0b2f96a 100644 --- a/src/heom_matrices/M_S.jl +++ b/src/heom_matrices/M_S.jl @@ -56,7 +56,7 @@ Note that the parity only need to be set as `ODD` when the system contains fermi println("Constructing Liouville-von Neumann superoperator...") flush(stdout) end - Lsys = -1im * (spre(_Hsys) - spost(_Hsys)) + Lsys = minus_i_L_op(_Hsys) if verbose println("[DONE]") flush(stdout) diff --git a/src/heom_matrices/heom_matrix_base.jl b/src/heom_matrices/heom_matrix_base.jl index bec30bdd..c971afa7 100644 --- a/src/heom_matrices/heom_matrix_base.jl +++ b/src/heom_matrices/heom_matrix_base.jl @@ -1,39 +1,102 @@ -abstract type AbstractHEOMLSMatrix end -abstract type AbstractParity end +@doc raw""" + struct HEOMSuperOp +General HEOM superoperator matrix. + +# Fields +- `data` : the HEOM superoperator matrix +- `dim` : the dimension of the system +- `N` : the number of auxiliary density operators +- `parity`: the parity label (`EVEN` or `ODD`). +""" +struct HEOMSuperOp + data::SparseMatrixCSC{ComplexF64, Int64} + dim::Int + N::Int + parity::AbstractParity +end @doc raw""" - struct OddParity <: AbstractParity + HEOMSuperOp(op, opParity, refHEOMLS, mul_basis="L") +Construct the HEOM superoperator matrix corresponding to the given system operator which acts on all `ADOs`. + +During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator. + +# Parameters +- `op` : The system operator which will act on all `ADOs`. +- `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`. +- `refHEOMLS::AbstractHEOMLSMatrix` : copy the system `dim` and number of `ADOs` (`N`) from this reference HEOMLS matrix +- `mul_basis::AbstractString` : this specifies the basis for `op` to multiply on all `ADOs`. Defaults to `"L"`. + +if `mul_basis` is specified as +- `"L"` : the matrix `op` has same dimension with the system and acts on left-hand side. +- `"R"` : the matrix `op` has same dimension with the system and acts on right-hand side. +- `"LR"` : the matrix `op` is a superoperator of the system. """ -struct OddParity <: AbstractParity end +HEOMSuperOp(op, opParity::AbstractParity, refHEOMLS::AbstractHEOMLSMatrix, mul_basis::AbstractString="L") = HEOMSuperOp(op, opParity, refHEOMLS.dim, refHEOMLS.N, mul_basis) @doc raw""" - struct EvenParity <: AbstractParity + HEOMSuperOp(op, opParity, refADOs, mul_basis="L") +Construct the HEOMLS matrix corresponding to the given system operator which multiplies on the "L"eft-hand ("R"ight-hand) side basis of all `ADOs`. + +During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator. + +# Parameters +- `op` : The system operator which will act on all `ADOs`. +- `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`. +- `refADOs::ADOs` : copy the system `dim` and number of `ADOs` (`N`) from this reference `ADOs` +- `mul_basis::AbstractString` : this specifies the basis for `op` to multiply on all `ADOs`. Defaults to `"L"`. + +if `mul_basis` is specified as +- `"L"` : the matrix `op` has same dimension with the system and acts on left-hand side. +- `"R"` : the matrix `op` has same dimension with the system and acts on right-hand side. +- `"LR"` : the matrix `op` is a superoperator of the system. """ -struct EvenParity <: AbstractParity end +HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString="L") = HEOMSuperOp(op, opParity, refADOs.dim, refADOs.N, mul_basis) -value(p::OddParity) = 1 -value(p::EvenParity) = 0 +@doc raw""" + HEOMSuperOp(op, opParity, dim, N, mul_basis) +Construct the HEOM superoperator matrix corresponding to the given system operator which acts on all `ADOs`. -function show(io::IO, p::OddParity) - print(io, "odd-parity") -end +During the multiplication on all the `ADOs`, the parity of the output `ADOs` might change depend on the parity of this HEOM superoperator. -function show(io::IO, p::EvenParity) - print(io, "even-parity") +# Parameters +- `op` : The system operator which will act on all `ADOs`. +- `opParity::AbstractParity` : the parity label of the given operator (`op`), should be `EVEN` or `ODD`. +- `dim::Int` : the system dimension. +- `N::Int` : the number of `ADOs`. +- `mul_basis::AbstractString` : this specifies the basis for `op` to multiply on all `ADOs`. + +if `mul_basis` is specified as +- `"L"` : the matrix `op` has same dimension with the system and acts on left-hand side. +- `"R"` : the matrix `op` has same dimension with the system and acts on right-hand side. +- `"LR"` : the matrix `op` is a superoperator of the system. +""" +function HEOMSuperOp(op, opParity::AbstractParity, dim::Int, N::Int, mul_basis::AbstractString) + if mul_basis == "L" + sup_op = spre(HandleMatrixType(op, dim, "op (operator)")) + elseif mul_basis == "R" + sup_op = spost(HandleMatrixType(op, dim, "op (operator)")) + elseif mul_basis == "LR" + sup_op = HandleMatrixType(op, dim ^ 2, "op (operator)") + else + error("The multiplication basis (mul_basis) can only be given as a string with either \"L\", \"R\", or \"LR\".") + end + + HEOMLS = kron(sparse(one(ComplexF64) * I, N, N), sup_op) + return HEOMSuperOp(HEOMLS, dim, N, opParity) end -# Parity label @doc raw""" - const ODD = OddParity() -Label of odd-parity + size(M::HEOMSuperOp) +Returns the size of the HEOM superoperator matrix """ -const ODD = OddParity() +size(M::HEOMSuperOp) = size(M.data) @doc raw""" - const EVEN = EvenParity() -Label of even-parity + size(M::HEOMSuperOp, dim::Int) +Returns the specified dimension of the HEOM superoperator matrix """ -const EVEN = EvenParity() +size(M::HEOMSuperOp, dim::Int) = size(M.data, dim) @doc raw""" size(M::AbstractHEOMLSMatrix) @@ -47,14 +110,31 @@ Returns the specified dimension of the HEOM Liouvillian superoperator matrix """ size(M::AbstractHEOMLSMatrix, dim::Int) = size(M.data, dim) +@doc raw""" + eltype(M::HEOMSuperOp) +Returns the elements' type of the HEOM superoperator matrix +""" +eltype(M::HEOMSuperOp) = eltype(M.data) + @doc raw""" eltype(M::AbstractHEOMLSMatrix) Returns the elements' type of the HEOM Liouvillian superoperator matrix """ eltype(M::AbstractHEOMLSMatrix) = eltype(M.data) +getindex(M::HEOMSuperOp, i::Ti, j::Tj) where {Ti, Tj <: Any} = M.data[i, j] getindex(M::AbstractHEOMLSMatrix, i::Ti, j::Tj) where {Ti, Tj <: Any} = M.data[i, j] +function show(io::IO, M::HEOMSuperOp) + print(io, + "$(M.parity) HEOM superoperator matrix acting on arbitrary-parity-ADOs\n", + "system dim = $(M.dim)\n", + "number of ADOs N = $(M.N)\n", + "data =\n" + ) + show(io, MIME("text/plain"), M.data) +end + function show(io::IO, M::AbstractHEOMLSMatrix) T = typeof(M) if T <: M_S @@ -76,8 +156,52 @@ function show(io::IO, M::AbstractHEOMLSMatrix) show(io, MIME("text/plain"), M.data) end +function show(io::IO, m::MIME"text/plain", M::HEOMSuperOp) show(io, M) end function show(io::IO, m::MIME"text/plain", M::AbstractHEOMLSMatrix) show(io, M) end +function *(Sup::HEOMSuperOp, ados::ADOs) + _check_sys_dim_and_ADOs_num(Sup, ados) + + return ADOs(Sup.data * ados.data, ados.dim, ados.N, Sup.parity * ados.parity) +end + +function *(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp) + _check_sys_dim_and_ADOs_num(Sup1, Sup2) + + return HEOMSuperOp(Sup1.data * Sup2.data, Sup1.dim, Sup1.N, Sup1.parity * Sup2.parity) +end + +*(n::Number, Sup::HEOMSuperOp) = HEOMSuperOp(n * Sup.data, Sup.dim, Sup.N, Sup.parity) +*(Sup::HEOMSuperOp, n::Number) = n * Sup + +function +(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp) + _check_sys_dim_and_ADOs_num(Sup1, Sup2) + _check_parity(Sup1, Sup2) + + return HEOMSuperOp(Sup1.data + Sup2.data, Sup1.dim, Sup1.N, Sup1.parity) +end + +function +(M::AbstractHEOMLSMatrix, Sup::HEOMSuperOp) + _check_sys_dim_and_ADOs_num(M, Sup) + _check_parity(M, Sup) + + return _reset_HEOMLS_data(M, M.data + Sup.data) +end + +function -(Sup1::HEOMSuperOp, Sup2::HEOMSuperOp) + _check_sys_dim_and_ADOs_num(Sup1, Sup2) + _check_parity(Sup1, Sup2) + + return HEOMSuperOp(Sup1.data - Sup2.data, Sup1.dim, Sup1.N, Sup1.parity) +end + +function -(M::AbstractHEOMLSMatrix, Sup::HEOMSuperOp) + _check_sys_dim_and_ADOs_num(M, Sup) + _check_parity(M, Sup) + + return _reset_HEOMLS_data(M, M.data - Sup.data) +end + @doc raw""" Propagator(M, Δt; threshold, nonzero_tol) Use `FastExpm.jl` to calculate the propagator matrix from a given HEOM Liouvillian superoperator matrix ``M`` with a specific time step ``\Delta t``. @@ -103,6 +227,18 @@ For more details, please refer to [`FastExpm.jl`](https://github.com/fmentink/Fa return fastExpm(M.data * Δt; threshold=threshold, nonzero_tol=nonzero_tol) end +function _reset_HEOMLS_data(M::T, new_data::SparseMatrixCSC{ComplexF64, Int64}) where T <: AbstractHEOMLSMatrix + if T <: M_S + return M_S(new_data, M.tier, M.dim, M.N, M.sup_dim, M.parity) + elseif T <: M_Boson + return M_Boson(new_data, M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy) + elseif T <: M_Fermion + return M_Fermion(new_data, M.tier, M.dim, M.N, M.sup_dim, M.parity, M.bath, M.hierarchy) + else + return M_Boson_Fermion(new_data, M.Btier, M.Ftier, M.dim, M.N, M.sup_dim, M.parity, M.Bbath, M.Fbath, M.hierarchy) + end +end + @doc raw""" addBosonDissipator(M, jumpOP) Adding bosonic dissipator to a given HEOMLS matrix which describes how the system dissipatively interacts with an extra bosonic environment. @@ -121,7 +257,7 @@ Note that if ``V`` is acting on fermionic systems, it should be even-parity to b # Return - `M_new::AbstractHEOMLSMatrix` : the new HEOM Liouvillian superoperator matrix """ -function addBosonDissipator(M::T, jumpOP::Vector=[]) where T <: AbstractHEOMLSMatrix +function addBosonDissipator(M::AbstractHEOMLSMatrix, jumpOP::Vector=[]) if length(jumpOP) > 0 L = spzeros(ComplexF64, M.sup_dim, M.sup_dim) for J in jumpOP @@ -129,15 +265,7 @@ function addBosonDissipator(M::T, jumpOP::Vector=[]) where T <: AbstractHEOMLSMa L += spre(_J) * spost(_J') - 0.5 * (spre(_J' * _J) + spost(_J' * _J)) end - 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 - return M_Boson(M.data + kron(sparse(one(ComplexF64) * 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 - return M_Fermion(M.data + kron(sparse(one(ComplexF64) * 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(one(ComplexF64) * 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) - end + return M + HEOMSuperOp(L, M.parity, M, "LR") end end @@ -173,15 +301,8 @@ function addFermionDissipator(M::T, jumpOP::Vector=[]) where T <: AbstractHEOMLS _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 - return M_S(M.data + L, M.tier, M.dim, M.N, M.sup_dim, M.parity) - elseif T <: M_Boson - return M_Boson(M.data + kron(sparse(one(ComplexF64) * 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 - return M_Fermion(M.data + kron(sparse(one(ComplexF64) * 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(one(ComplexF64) * 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) - end + + return M + HEOMSuperOp(L, M.parity, M, "LR") end end @@ -230,13 +351,7 @@ function addTerminator(M::Mtype, Bath::Union{BosonBath, FermionBath}) where Mtyp spre(J) * spost(J') - 0.5 * (spre(J' * J) + spost(J' * J)) ) - if Mtype <: M_Boson - return M_Boson(M.data + kron(sparse(one(ComplexF64) * 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 - return M_Fermion(M.data + kron(sparse(one(ComplexF64) * 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(one(ComplexF64) * 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) - end + return M + HEOMSuperOp(L, M.parity, M, "LR") end end @@ -303,67 +418,72 @@ function bath_sum_γ(nvec, baths::Vector{T}) where T <: Union{AbstractBosonBath, return sum_γ end +# commutator of system Hamiltonian +function minus_i_L_op(Hsys) + return -1.0im * (spre(Hsys) - spost(Hsys)) +end + # connect to bosonic (n-1)th-level for "Real & Imag combined operator" -function _D_op(bath::bosonRealImag, k, n_k) +function minus_i_D_op(bath::bosonRealImag, k, n_k) return n_k * ( - -1im * bath.η_real[k] * bath.Comm + - bath.η_imag[k] * bath.anComm + -1.0im * bath.η_real[k] * bath.Comm + + bath.η_imag[k] * bath.anComm ) end # connect to bosonic (n-1)th-level for (Real & Imag combined) operator "Real operator" -function _D_op(bath::bosonReal, k, n_k) - return -1im * n_k * bath.η[k] * bath.Comm +function minus_i_D_op(bath::bosonReal, k, n_k) + return -1.0im * n_k * bath.η[k] * bath.Comm end # connect to bosonic (n-1)th-level for "Imag operator" -function _D_op(bath::bosonImag, k, n_k) +function minus_i_D_op(bath::bosonImag, k, n_k) return n_k * bath.η[k] * bath.anComm end # connect to bosonic (n-1)th-level for "Absorption operator" -function _D_op(bath::bosonAbsorb, k, n_k) - return -1im * n_k * ( +function minus_i_D_op(bath::bosonAbsorb, k, n_k) + return -1.0im * n_k * ( bath.η[k] * bath.spre - conj(bath.η_emit[k]) * bath.spost ) end # connect to bosonic (n-1)th-level for "Emission operator" -function _D_op(bath::bosonEmit, k, n_k) - return -1im * n_k * ( +function minus_i_D_op(bath::bosonEmit, k, n_k) + return -1.0im * n_k * ( bath.η[k] * bath.spre - conj(bath.η_absorb[k]) * bath.spost ) end # connect to fermionic (n-1)th-level for "absorption operator" -function _C_op(bath::fermionAbsorb, k, n_exc, n_exc_before, parity) - return -1im * ((-1) ^ n_exc_before) * ( +function minus_i_C_op(bath::fermionAbsorb, k, n_exc, n_exc_before, parity) + return -1.0im * ((-1) ^ n_exc_before) * ( ((-1) ^ value(parity)) * bath.η[k] * bath.spre - (-1) ^ (n_exc - 1) * conj(bath.η_emit[k]) * bath.spost ) end # connect to fermionic (n-1)th-level for "emission operator" -function _C_op(bath::fermionEmit, k, n_exc, n_exc_before, parity) - return -1im * ((-1) ^ n_exc_before) * ( +function minus_i_C_op(bath::fermionEmit, k, n_exc, n_exc_before, parity) + return -1.0im * ((-1) ^ n_exc_before) * ( (-1) ^ (value(parity)) * bath.η[k] * bath.spre - (-1) ^ (n_exc - 1) * conj(bath.η_absorb[k]) * bath.spost ) end # connect to bosonic (n+1)th-level for real-and-imaginary-type bosonic bath -function _B_op(bath::T) where T <: Union{bosonReal, bosonImag, bosonRealImag} - return -1im * bath.Comm +function minus_i_B_op(bath::T) where T <: Union{bosonReal, bosonImag, bosonRealImag} + return -1.0im * bath.Comm end # connect to bosonic (n+1)th-level for absorption-and-emission-type bosonic bath -function _B_op(bath::T) where T <: Union{bosonAbsorb, bosonEmit} - return -1im * bath.CommD +function minus_i_B_op(bath::T) where T <: Union{bosonAbsorb, bosonEmit} + return -1.0im * bath.CommD end # connect to fermionic (n+1)th-level -function _A_op(bath::T, n_exc, n_exc_before, parity) where T <: AbstractFermionBath - return -1im * ((-1) ^ n_exc_before) * ( +function minus_i_A_op(bath::T, n_exc, n_exc_before, parity) where T <: AbstractFermionBath + return -1.0im * ((-1) ^ n_exc_before) * ( (-1) ^ (value(parity)) * bath.spreD + (-1) ^ (n_exc + 1) * bath.spostD ) diff --git a/src/power_spectrum.jl b/src/power_spectrum.jl new file mode 100644 index 00000000..a7933a74 --- /dev/null +++ b/src/power_spectrum.jl @@ -0,0 +1,179 @@ +@doc raw""" + spectrum(M, ρ, op, ωlist; solver, verbose, filename, SOLVEROptions...) +!!! warning "Warning" + This function has been deprecated start from `HierarchicalEOM v1.1`, use `PowerSpectrum` or `DensityOfStates` instead. +""" +function spectrum( + M::AbstractHEOMLSMatrix, + ρ, + op, + ωlist::AbstractVector; + solver=UMFPACKFactorization(), + verbose::Bool = true, + filename::String = "", + SOLVEROptions... + ) + error("This function has been deprecated start from \`HierarchicalEOM v1.1\`, use \`PowerSpectrum\` or \`DensityOfStates\` instead.") +end + +@doc raw""" + PowerSpectrum(M, ρ, Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...) +Calculate power spectrum for the system in frequency domain where `P_op` will be automatically set as the adjoint of `Q_op`. + +This function is equivalent to: +`PowerSpectrum(M, ρ, Q_op', Q_op, ωlist, reverse; solver, verbose, filename, SOLVEROptions...)` +""" +function PowerSpectrum( + M::AbstractHEOMLSMatrix, + ρ, + Q_op, + ωlist::AbstractVector, + reverse::Bool = false; + solver=UMFPACKFactorization(), + verbose::Bool = true, + filename::String = "", + SOLVEROptions... + ) + HandleMatrixType(Q_op, M.dim) + 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...) +Calculate power spectrum for the system in frequency domain. + +```math +\pi S(\omega)=\textrm{Re}\left\{\int_0^\infty dt \langle P(t) Q(0)\rangle e^{-i\omega t}\right\}, +``` + +# To calculate spectrum when input operator `Q_op` has `EVEN`-parity: +remember to set the parameters: +- `M::AbstractHEOMLSMatrix`: should be `EVEN` parity + +# To calculate spectrum when input operator `Q_op` has `ODD`-parity: +remember to set the parameters: +- `M::AbstractHEOMLSMatrix`: should be `ODD` parity + +# 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. +- `ω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` : 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`. +- `filename::String` : If filename was specified, the value of spectrum for each ω will be saved into the file "filename.txt" during the solving process. +- `SOLVEROptions` : extra options for solver + +For more details about solvers and extra options, please refer to [`LinearSolve.jl`](http://linearsolve.sciml.ai/stable/) + +# Returns +- `spec::AbstractVector` : the spectrum list corresponds to the specified `ωlist` +""" +@noinline function PowerSpectrum( + M::AbstractHEOMLSMatrix, + ρ, + P_op, + Q_op, + ωlist::AbstractVector, + reverse::Bool = false; + solver = UMFPACKFactorization(), + verbose::Bool = true, + filename::String = "", + SOLVEROptions... + ) + + Size = size(M, 1) + + # Handle ρ + if typeof(ρ) == ADOs # ρ::ADOs + ados = ρ + else + ados = ADOs(ρ, M.N) + end + _check_sys_dim_and_ADOs_num(M, ados) + + # Handle P_op + if typeof(P_op) == HEOMSuperOp + _check_sys_dim_and_ADOs_num(M, P_op) + _P = P_op + else + _P = HEOMSuperOp(P_op, EVEN, M) + end + _tr_P = Tr(M.dim, M.N) * _P.data + + # Handle Q_op + if typeof(Q_op) == HEOMSuperOp + _check_sys_dim_and_ADOs_num(M, Q_op) + _Q_ados = Q_op * ados + _check_parity(M, _Q_ados) + else + if M.parity == EVEN + _Q = HEOMSuperOp(Q_op, ados.parity, M) + else + _Q = HEOMSuperOp(Q_op, !ados.parity, M) + end + _Q_ados = _Q * ados + end + b = _HandleVectorType(typeof(M.data), _Q_ados.data) + + SAVE::Bool = (filename != "") + if SAVE + FILENAME = filename * ".txt" + if isfile(FILENAME) + error("FILE: $(FILENAME) already exist.") + end + end + + ElType = eltype(M) + ωList = _HandleFloatType(ElType, ωlist) + Length = length(ωList) + Sω = Vector{Float64}(undef, Length) + + if verbose + print("Calculating power spectrum in frequency domain...\n") + flush(stdout) + prog = Progress(Length; desc="Progress : ", PROGBAR_OPTIONS...) + end + + if reverse + i = convert(ElType, 1im) + else + i = convert(ElType, -1im) + end + I_total = _HandleIdentityType(typeof(M.data), Size) + Iω = i * ωList[1] * I_total + cache = init(LinearProblem(M.data + Iω, b), solver, SOLVEROptions...) + sol = solve!(cache) + @inbounds for (j, ω) in enumerate(ωList) + if j > 1 + Iω = i * ω * I_total + cache.A = M.data + Iω + sol = solve!(cache) + end + + # trace over the Hilbert space of system (expectation value) + Sω[j] = -1 * real(_tr_P * _HandleVectorType(sol.u, false)) + + if SAVE + open(FILENAME, "a") do file + write(file, "$(Sω[j]),\n") + end + end + + if verbose + next!(prog) + end + end + if verbose + println("[DONE]") + end + + return Sω +end \ No newline at end of file diff --git a/test/ADOs.jl b/test/ADOs.jl new file mode 100644 index 00000000..9f3a54b7 --- /dev/null +++ b/test/ADOs.jl @@ -0,0 +1,19 @@ +ados_b = ADOs(spzeros(Int64, 20), 5) +ados_f = ADOs(spzeros(Int64, 8), 2) +ados_bf = ADOs(spzeros(Int64, 40), 10) +@test show(devnull, MIME("text/plain"), ados_b) == nothing +@test show(devnull, MIME("text/plain"), ados_f) == nothing +@test show(devnull, MIME("text/plain"), ados_bf) == nothing +@test_throws ErrorException ADOs(zeros(8), 4) + +ρ_b = ados_b[:] +# check iteration +for (i, ado) in enumerate(ados_b) + @test ρ_b[i] == ado +end + +# expections for expect +ados_wrong = ADOs(spzeros(Int64, 18), 2) +@test_throws ErrorException Expect([0 0 0; 0 0 0; 0 0 0], ados_f) +@test_throws ErrorException Expect([0 0; 0 0], [ados_b, ados_wrong]) +@test_throws ErrorException Expect([0 0 0; 0 0 0; 0 0 0], [ados_b, ados_f]) \ No newline at end of file diff --git a/test/HEOMSuperOp.jl b/test/HEOMSuperOp.jl new file mode 100644 index 00000000..055cacab --- /dev/null +++ b/test/HEOMSuperOp.jl @@ -0,0 +1,103 @@ +# Take waiting time distribution as an example +WTD_ans = [ + 0.0, + 0.0015096807591478143, + 0.0019554827640519456, + 0.0019720219482432183, + 0.0018278123354472694, + 0.0016348729788834187, + 0.0014385384773227094, + 0.0012559122034857784, + 0.0010923035083891884, + 0.0009482348815344435, + 0.0008224107708326052, + 0.0007129580704171314, + 0.0006179332764127152, + 0.0005355140701828154, + 0.00046406230317384855, + 0.00040213314306247375, + 0.0003484637380485898, + 0.0003019551325507645, + 0.00026165305303984653, + 0.0002267297384879625, + 0.0001964675440128574 +] + +ϵ = -0.1 +Γ = 0.01 +μ = 0 +W = 1 +kT = 0.5 +N = 5 +tier = 3 +tlist = 0:50:1000 +γ_eL = γ_eR = 0.06676143449267714 # emission rate towards Left and Right lead +γ_aL = γ_aR = 0.0737827958503192 # absorption rate towards Left and Right lead + +d = [0 1; 0 0] +Hs = ϵ * d' * d + +## Only local master equation +M_me = M_S(Hs; verbose=false) +jumpOPs = [γ_eL * d, γ_eR * d, γ_aL * d', γ_aR * d'] +M_me = addFermionDissipator(M_me, jumpOPs) + +# create quantum jump superoperator (only the emission part) +J_me = HEOMSuperOp(γ_eR * d, ODD, M_me, "L") * HEOMSuperOp(γ_eR * d', ODD, M_me, "R") +@test size(M_me) == size(J_me) +@test size(M_me, 1) == size(J_me, 1) +@test eltype(M_me) == eltype(J_me) +@test nnz(J_me.data) == 1 +@test J_me[1, 4] ≈ 0.0044570891355200214 +@test J_me.data * 2 ≈ (J_me - (-1) * J_me).data +@test J_me.data * 2 ≈ (HEOMSuperOp(√(2) * γ_eR * d, ODD, M_me, "L") * HEOMSuperOp(√(2) * γ_eR * d', ODD, M_me, "R")).data +@test show(devnull, MIME("text/plain"), J_me) == nothing + +# create HEOM Liouvuillian superoperator without quantum jump +M0 = M_me - J_me +ados_s = SteadyState(M_me; verbose=false) +ados_list = evolution(M0, J_me * ados_s, tlist; verbose=false) + +# calculating waiting time distribution +WTD_me = [] +trJρ = Expect(J_me, ados_s) +WTD_me = Expect(J_me, ados_list) ./ trJρ +for i in 1:length(tlist) + @test WTD_me[i] ≈ WTD_ans[i] +end +@test Expect(J_me, ados_list[end]) / trJρ ≈ WTD_ans[end] + +## Left lead with HEOM method and Right lead with local master equation +bath = Fermion_Lorentz_Pade(d, Γ, μ, W, kT, N) +jumpOPs = [γ_eR * d, γ_aR * d'] +M_heom = M_Fermion(Hs, tier, bath; verbose=false) +M_heom = addFermionDissipator(M_heom, jumpOPs) + +# create quantum jump superoperator (only the emission part) +J_heom = HEOMSuperOp(γ_eR * d, ODD, M_heom, "L") * HEOMSuperOp(γ_eR * d', ODD, M_heom, "R") + +# create HEOM Liouvuillian superoperator without quantum jump +M1 = M_heom - J_heom +ados_s1 = SteadyState(M_heom; verbose=false) +ados_list1 = evolution(M1, J_heom * ados_s1, tlist; verbose=false) + +# calculating waiting time distribution +WTD_heom = [] +trJρ1 = Expect(J_heom, ados_s1) +WTD_heom = Expect(J_heom, ados_list1) ./ trJρ1 +for i in 1:length(tlist) + @test WTD_heom[i] ≈ WTD_ans[i] atol=1e-5 +end + +## test for ERRORs +J_wrong1 = HEOMSuperOp(γ_eR * d, ODD, M_me, "L") +J_wrong2 = HEOMSuperOp(γ_eR * d, EVEN, M_heom, "L") +@test_throws ErrorException J_me * J_wrong2 +@test_throws ErrorException J_me + J_wrong1 +@test_throws ErrorException J_me + J_wrong2 +@test_throws ErrorException M_me + J_wrong1 +@test_throws ErrorException M_me + J_wrong2 +@test_throws ErrorException J_me - J_wrong1 +@test_throws ErrorException J_me - J_wrong2 +@test_throws ErrorException M_me - J_wrong1 +@test_throws ErrorException M_me - J_wrong2 \ No newline at end of file diff --git a/test/M_Boson.jl b/test/M_Boson.jl new file mode 100644 index 00000000..17839be5 --- /dev/null +++ b/test/M_Boson.jl @@ -0,0 +1,65 @@ +# Test Boson-type HEOM Liouvillian superoperator matrix +λ = 0.1450 +W = 0.6464 +kT = 0.7414 +μ = 0.8787 +N = 5 +tier = 3 + +# System Hamiltonian +Hsys = [ + 0.6969 0.4364; + 0.4364 0.3215 +] + +# system-bath coupling operator +Q = [ + 0.1234 0.1357 + 0.2468im; + 0.1357 - 0.2468im 0.5678 +] +Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) + +# jump operator +J = [0 0.1450 - 0.7414im; 0.1450 + 0.7414im 0] + +L = M_Boson(Hsys, tier, Bbath; verbose=false) +@test show(devnull, MIME("text/plain"), L) == nothing +@test size(L) == (336, 336) +@test L.N == 84 +@test nnz(L.data) == 4422 +L = addBosonDissipator(L, J) +@test nnz(L.data) == 4760 +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 = [ + 0.4969521584882579 - 2.27831302340618e-13im -0.0030829715611090133 + 0.002534368458048467im; + -0.0030829715591718203 - 0.0025343684616701547im 0.5030478415140676 + 2.3661885315257474e-13im +] +@test _is_Matrix_approx(ρ0, ρ1) + +L = M_Boson(Hsys, tier, [Bbath, Bbath]; verbose=false) +@test size(L) == (1820, 1820) +@test L.N == 455 +@test nnz(L.data) == 27662 +L = addBosonDissipator(L, J) +@test nnz(L.data) == 29484 +ados = SteadyState(L; verbose=false) +@test ados.dim == L.dim +@test length(ados) == L.N +ρ0 = ados[1] +@test getRho(ados) == ρ0 +ρ1 = [ + 0.49406682844513267 + 9.89558173111355e-13im -0.005261234545120281 + 0.0059968903987593im; + -0.005261234550122085 - 0.005996890386139547im 0.5059331715578721 - 9.413847493320824e-13im +] +@test _is_Matrix_approx(ρ0, ρ1) + +## check exceptions +@test_throws BoundsError L[1, 1821] +@test_throws BoundsError L[1:1821, 336] +@test_throws ErrorException ados[L.N + 1] +@test_throws ErrorException M_Boson([0, 0], tier, Bbath; verbose=false) \ No newline at end of file diff --git a/test/M_Boson_Fermion.jl b/test/M_Boson_Fermion.jl new file mode 100644 index 00000000..b7b4818e --- /dev/null +++ b/test/M_Boson_Fermion.jl @@ -0,0 +1,85 @@ +# Test Boson-Fermion-type HEOM Liouvillian superoperator matrix +λ = 0.1450 +W = 0.6464 +kT = 0.7414 +μ = 0.8787 +N = 3 +tierb = 2 +tierf = 2 + +Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) +Fbath = Fermion_Lorentz_Pade(Q, λ, μ, W, kT, N) + +# System Hamiltonian +Hsys = [ + 0.6969 0.4364; + 0.4364 0.3215 +] + +# system-bath coupling operator +Q = [ + 0.1234 0.1357 + 0.2468im; + 0.1357 - 0.2468im 0.5678 +] +Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) +Fbath = Fermion_Lorentz_Pade(Q, λ, μ, W, kT, N) + +# jump operator +J = [0 0.1450 - 0.7414im; 0.1450 + 0.7414im 0] + +L = M_Boson_Fermion(Hsys, tierb, tierf, Bbath, Fbath; verbose=false) +@test show(devnull, MIME("text/plain"), L) == nothing +@test size(L) == (2220, 2220) +@test L.N == 555 +@test nnz(L.data) == 43368 +L = addBosonDissipator(L, J) +@test nnz(L.data) == 45590 +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 = [ + 0.496709+4.88415e-12im -0.00324048+0.00286376im; +-0.00324048-0.00286376im 0.503291-4.91136e-12im +] +@test _is_Matrix_approx(ρ0, ρ1) + +L = M_Boson_Fermion(Hsys, tierb, tierf, [Bbath, Bbath], Fbath; verbose=false) +@test size(L) == (6660, 6660) +@test L.N == 1665 +@test nnz(L.data) == 139210 +L = addFermionDissipator(L, J) +@test nnz(L.data) == 145872 +ados = SteadyState(L; verbose=false) +@test ados.dim == L.dim +@test length(ados) == L.N +ρ0 = ados[1] +@test getRho(ados) == ρ0 +ρ1 = [ + 0.493774+6.27624e-13im -0.00536526+0.00651746im; + -0.00536526-0.00651746im 0.506226-6.15855e-13im +] +@test _is_Matrix_approx(ρ0, ρ1) + +L = M_Boson_Fermion(Hsys, tierb, tierf, Bbath, [Fbath, Fbath]; verbose=false) +@test size(L) == (8220, 8220) +@test L.N == 2055 +@test nnz(L.data) == 167108 +L = addBosonDissipator(L, J) +@test nnz(L.data) == 175330 +ados = SteadyState(L; verbose=false) +@test ados.dim == L.dim +@test length(ados) == L.N +ρ0 = ados[1] +@test getRho(ados) == ρ0 +ρ1 = [ + 0.496468-4.32253e-12im -0.00341484+0.00316445im; +-0.00341484-0.00316445im 0.503532+4.32574e-12im +] +@test _is_Matrix_approx(ρ0, ρ1) + +## check exceptions +@test_throws ErrorException ados[L.N + 1] +@test_throws ErrorException M_Boson_Fermion([0, 0], tierb, tierf, Bbath, Fbath; verbose=false) \ No newline at end of file diff --git a/test/M_Boson_RWA.jl b/test/M_Boson_RWA.jl new file mode 100644 index 00000000..591975a2 --- /dev/null +++ b/test/M_Boson_RWA.jl @@ -0,0 +1,24 @@ +# Test Boson-type HEOM Liouvillian superoperator matrix under rotating wave approximation +ωq = 1.1 +Λ = 0.01 +Γ = 0.02 +Hsys_rwa = 0.5 * ωq * [1 0; 0 -1] +op_rwa = [0 0; 1 0] +ρ0 = 0.5 * [1 1; 1 1] + +tlist = 0:1:20 +d = 1im * √(Λ * (2 * Γ - Λ)) # non-Markov regime + +B_rwa = BosonBathRWA(op_rwa, [0], [Λ - 1im * ωq], [0.5 * Γ * Λ], [Λ + 1im * ωq]) +L = M_Boson(Hsys_rwa, tier, B_rwa; verbose=false) +ados_list = evolution(L, ρ0, tlist; reltol=1e-10, abstol=1e-12, verbose=false) + +for (i, t) in enumerate(tlist) + ρ_rwa = getRho(ados_list[i]) + + # analytical result + Gt = exp(-1 * (Λ / 2 + 1im * ωq) * t) * (cosh(t * d / 2) + Λ * sinh(t * d / 2) / d) + + @test ρ_rwa[1, 1] ≈ abs(Gt) ^ 2 * ρ0[1, 1] + @test ρ_rwa[1, 2] ≈ Gt * ρ0[1, 2] +end \ No newline at end of file diff --git a/test/M_Fermion.jl b/test/M_Fermion.jl new file mode 100644 index 00000000..10e31499 --- /dev/null +++ b/test/M_Fermion.jl @@ -0,0 +1,75 @@ +# Test Fermion-type HEOM Liouvillian superoperator matrix +λ = 0.1450 +W = 0.6464 +kT = 0.7414 +μ = 0.8787 +N = 5 +tier = 3 + +# System Hamiltonian +Hsys = [ + 0.6969 0.4364; + 0.4364 0.3215 +] + +# system-bath coupling operator +Q = [ + 0.1234 0.1357 + 0.2468im; + 0.1357 - 0.2468im 0.5678 +] +Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) +Fbath = Fermion_Lorentz_Pade(Q, λ, μ, W, kT, N) + +# jump operator +J = [0 0.1450 - 0.7414im; 0.1450 + 0.7414im 0] + +L = M_Fermion(Hsys, tier, Fbath; verbose=false) +@test show(devnull, MIME("text/plain"), L) == nothing +@test size(L) == (1196, 1196) +@test L.N == 299 +@test nnz(L.data) == 21318 +L = addFermionDissipator(L, J) +@test nnz(L.data) == 22516 +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 = [ + 0.49971864340781574 + 1.5063528845574697e-11im -0.00025004129095353573 + 0.00028356932981729176im; + -0.0002500413218393161 - 0.0002835693203755187im 0.5002813565929579 - 1.506436545359778e-11im +] +@test _is_Matrix_approx(ρ0, ρ1) + +L = M_Fermion(Hsys, tier, Fbath; threshold = 1e-8, verbose=false) +L = addFermionDissipator(L, J) +@test size(L) == (148, 148) +@test L.N == 37 +@test nnz(L.data) == 2054 +ados = SteadyState(L; verbose=false) +ρ2 = ados[1] +@test _is_Matrix_approx(ρ0, ρ2) + +L = M_Fermion(Hsys, tier, [Fbath, Fbath]; verbose=false) +@test size(L) == (9300, 9300) +@test L.N == 2325 +@test nnz(L.data) == 174338 +L = addFermionDissipator(L, J) +@test nnz(L.data) == 183640 +ados = SteadyState(L; verbose=false) +@test ados.dim == L.dim +@test length(ados) == L.N +ρ0 = ados[1] +@test getRho(ados) == ρ0 +ρ1 = [ + 0.4994229368103249 + 2.6656157051929377e-12im -0.0005219753638749304 + 0.0005685093274121244im; + -0.0005219753958601764 - 0.0005685092413099392im 0.5005770631903512 - 2.6376966158390854e-12im +] +@test _is_Matrix_approx(ρ0, ρ1) + +## check exceptions +@test_throws BoundsError L[1, 9301] +@test_throws BoundsError L[1:9301, 9300] +@test_throws ErrorException ados[L.N + 1] +@test_throws ErrorException M_Fermion([0, 0], tier, Fbath; verbose=false) \ No newline at end of file diff --git a/test/M_S.jl b/test/M_S.jl new file mode 100644 index 00000000..7e33a067 --- /dev/null +++ b/test/M_S.jl @@ -0,0 +1,39 @@ +# Test Schrodinger type HEOM Liouvillian superoperator matrix +t = 10 +Hsys = [0 1; 1 0] +L = M_S(Hsys; verbose=false) +@test show(devnull, MIME("text/plain"), L) == nothing +@test size(L) == (4, 4) +@test L.N == 1 +@test nnz(L.data) == 8 +ados_list = evolution(L, [1 0; 0 0], 0:1:t; reltol=1e-8, abstol=1e-10, verbose=false) +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 = [ + cos(t)^2 -0.5im*sin(2*t); + 0.5im*sin(2*t) sin(t)^2 +] +@test _is_Matrix_approx(ρ0, ρ1) + +L = addBosonDissipator(L, √(0.01) * [1 0; 0 -1]) +L = addFermionDissipator(L, √(0.01) * [1 0; 0 -1]) +@test nnz(L.data) == 10 +ados_list = evolution(L, [1 0; 0 0], 0:0.5:t; reltol=1e-8, abstol=1e-10, verbose=false) +ados = ados_list[end] +ρ0 = ados[1] +@test getRho(ados) == ρ0 +ρ1 = [ + 0.6711641119639493+0.0im 0.3735796014268062im; + -0.3735796014268062im 0.32883588803605024 +] +@test _is_Matrix_approx(ρ0, ρ1) + +## check exceptions +@test_throws BoundsError L[1, 5] +@test_throws BoundsError L[1:5, 2] +@test_throws ErrorException ados[L.N + 1] +@test_throws ErrorException M_S([0, 0]; verbose=false) \ No newline at end of file diff --git a/test/corr_func.jl b/test/bath_corr_func.jl similarity index 100% rename from test/corr_func.jl rename to test/bath_corr_func.jl diff --git a/test/density_of_states.jl b/test/density_of_states.jl new file mode 100644 index 00000000..f5fb1342 --- /dev/null +++ b/test/density_of_states.jl @@ -0,0 +1,75 @@ +e = -5 +U = 10 +d_up = kron( [0 1; 0 0], [1 0; 0 1]) +d_dn = kron(-1 * [1 0; 0 -1], [0 1; 0 0]) +iden = kron( [1 0; 0 1], [1 0; 0 1]) + +H0 = e * (d_up' * d_up + d_dn' * d_dn) +H1 = U * (d_up' * d_up * d_dn' * d_dn) +Hsys = H0 + H1 + +λ = 1 +μ_l = 1 +μ_r = -1 +W = 10 +kT = 0.5 +N = 5 +fuL = Fermion_Lorentz_Pade(d_up, λ, μ_l, W, kT, N) +fdL = Fermion_Lorentz_Pade(d_dn, λ, μ_l, W, kT, N) +fuR = Fermion_Lorentz_Pade(d_up, λ, μ_r, W, kT, N) +fdR = Fermion_Lorentz_Pade(d_dn, λ, μ_r, W, kT, N) + +tier = 2 +Le = M_Fermion(Hsys, tier, [fuL, fdL, fuR, fdR]; verbose=false) +Lo = M_Fermion(Hsys, tier, [fuL, fdL, fuR, fdR], ODD; verbose=false) + +ados_s = SteadyState(Le; verbose=false) +ωlist = -20:2:20 + +if isfile("DOS.txt") + rm("DOS.txt") +end +d_up_normal = HEOMSuperOp(d_up, ODD, Le) +dos1 = DensityOfStates(Lo, ados_s, d_up, ωlist; verbose=false, filename="DOS") +dos2 = PowerSpectrum(Lo, ados_s, d_up_normal, d_up', ωlist, true; verbose=false) .+ PowerSpectrum(Lo, ados_s, d_up', d_up_normal, ωlist, false; verbose=false) +dos3 = [ + 0.0007920428534358747, + 0.0012795202828027256, + 0.0022148985361417936, + 0.004203651086852703, + 0.009091831276694192, + 0.024145570990254425, + 0.09111929563034224, + 0.2984323182857298, + 0.1495897397559431, + 0.12433521300531171, + 0.17217519700362036, + 0.1243352130053117, + 0.14958973975594306, + 0.29843231828572964, + 0.09111929563034214, + 0.024145570990254432, + 0.009091831276694183, + 0.004203651086852702, + 0.0022148985361417923, + 0.0012795202828027236, + 0.0007920428534358735 +] +for i in 1:length(ωlist) + @test dos1[i] ≈ dos3[i] atol=1.0e-10 + @test dos2[i] ≈ dos3[i] atol=1.0e-10 +end + +mat = spzeros(ComplexF64, 2, 2) +mat2 = spzeros(ComplexF64, 3, 3) +bathb = Boson_DrudeLorentz_Pade(mat, 1, 1, 1, 2) +@test_throws ErrorException spectrum(Lo, ados_s, d_up, ωlist; verbose=false) +@test_throws ErrorException DensityOfStates(Lo, ados_s, mat2, ωlist; verbose=false) +@test_throws ErrorException DensityOfStates(Lo, ADOs(zeros(8), 2), d_up, ωlist; verbose=false) +@test_throws ErrorException DensityOfStates(Lo, ADOs(zeros(32), 2), d_up, ωlist; verbose=false) +@test_throws ErrorException DensityOfStates(Lo, ADOs(ados_s.data, ados_s.N, ODD), d_up, ωlist; verbose=false) +@test_throws ErrorException DensityOfStates(M_Boson(mat, 2, bathb; verbose=false), mat, mat, [0]) +@test_throws ErrorException DensityOfStates(M_Fermion(mat, 2, fuL; verbose=false), mat, mat, [0]) + +# remove all the temporary files +rm("DOS.txt") \ No newline at end of file diff --git a/test/heom_api.jl b/test/heom_api.jl deleted file mode 100644 index 2585b5d9..00000000 --- a/test/heom_api.jl +++ /dev/null @@ -1,347 +0,0 @@ -# Test Schrodinger type HEOM Liouvillian superoperator matrix -@testset "M_S" begin - t = 10 - Hsys = [0 1; 1 0] - L = M_S(Hsys; verbose=false) - @test show(devnull, MIME("text/plain"), L) == nothing - @test size(L) == (4, 4) - @test L.N == 1 - @test nnz(L.data) == 8 - ados_list = evolution(L, [1 0; 0 0], 0:1:t; reltol=1e-8, abstol=1e-10, verbose=false) - 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 = [ - cos(t)^2 -0.5im*sin(2*t); - 0.5im*sin(2*t) sin(t)^2 - ] - @test _is_Matrix_approx(ρ0, ρ1) - - L = addBosonDissipator(L, √(0.01) * [1 0; 0 -1]) - L = addFermionDissipator(L, √(0.01) * [1 0; 0 -1]) - @test nnz(L.data) == 10 - ados_list = evolution(L, [1 0; 0 0], 0:0.5:t; reltol=1e-8, abstol=1e-10, verbose=false) - ados = ados_list[end] - ρ0 = ados[1] - @test getRho(ados) == ρ0 - ρ1 = [ - 0.6711641119639493+0.0im 0.3735796014268062im; - -0.3735796014268062im 0.32883588803605024 - ] - @test _is_Matrix_approx(ρ0, ρ1) - - ## check exceptions - @test_throws BoundsError L[1, 5] - @test_throws BoundsError L[1:5, 2] - @test_throws ErrorException ados[L.N + 1] - @test_throws ErrorException @test_warn "HEOM doesn't support matrix type : Vector{Int64}" M_S([0, 0]; verbose=false) -end - -λ = 0.1450 -W = 0.6464 -kT = 0.7414 -μ = 0.8787 -N = 5 -tier = 3 - -# System Hamiltonian -Hsys = [ - 0.6969 0.4364; - 0.4364 0.3215 -] - -# system-bath coupling operator -Q = [ - 0.1234 0.1357 + 0.2468im; - 0.1357 - 0.2468im 0.5678 -] -Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) -Fbath = Fermion_Lorentz_Pade(Q, λ, μ, W, kT, N) - -# jump operator -J = [0 0.1450 - 0.7414im; 0.1450 + 0.7414im 0] - -# Test Boson-type HEOM Liouvillian superoperator matrix -@testset "M_Boson" begin - L = M_Boson(Hsys, tier, Bbath; verbose=false) - @test show(devnull, MIME("text/plain"), L) == nothing - @test size(L) == (336, 336) - @test L.N == 84 - @test nnz(L.data) == 4422 - L = addBosonDissipator(L, J) - @test nnz(L.data) == 4760 - 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 = [ - 0.4969521584882579 - 2.27831302340618e-13im -0.0030829715611090133 + 0.002534368458048467im; - -0.0030829715591718203 - 0.0025343684616701547im 0.5030478415140676 + 2.3661885315257474e-13im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - L = M_Boson(Hsys, tier, [Bbath, Bbath]; verbose=false) - @test size(L) == (1820, 1820) - @test L.N == 455 - @test nnz(L.data) == 27662 - L = addBosonDissipator(L, J) - @test nnz(L.data) == 29484 - ados = SteadyState(L; verbose=false) - @test ados.dim == L.dim - @test length(ados) == L.N - ρ0 = ados[1] - @test getRho(ados) == ρ0 - ρ1 = [ - 0.49406682844513267 + 9.89558173111355e-13im -0.005261234545120281 + 0.0059968903987593im; - -0.005261234550122085 - 0.005996890386139547im 0.5059331715578721 - 9.413847493320824e-13im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - ## check exceptions - @test_throws BoundsError L[1, 1821] - @test_throws BoundsError L[1:1821, 336] - @test_throws ErrorException ados[L.N + 1] - @test_throws ErrorException M_Boson([0, 0], tier, Bbath; verbose=false) -end - -@testset "M_Boson (RWA)" begin - ωq = 1.1 - Λ = 0.01 - Γ = 0.02 - Hsys_rwa = 0.5 * ωq * [1 0; 0 -1] - op_rwa = [0 0; 1 0] - ρ0 = 0.5 * [1 1; 1 1] - - tlist = 0:1:20 - d = 1im * √(Λ * (2 * Γ - Λ)) # non-Markov regime - - B_rwa = BosonBathRWA(op_rwa, [0], [Λ - 1im * ωq], [0.5 * Γ * Λ], [Λ + 1im * ωq]) - L = M_Boson(Hsys_rwa, tier, B_rwa; verbose=false) - ados_list = evolution(L, ρ0, tlist; reltol=1e-10, abstol=1e-12, verbose=false) - - for (i, t) in enumerate(tlist) - ρ_rwa = getRho(ados_list[i]) - - # analytical result - Gt = exp(-1 * (Λ / 2 + 1im * ωq) * t) * (cosh(t * d / 2) + Λ * sinh(t * d / 2) / d) - - @test ρ_rwa[1, 1] ≈ abs(Gt) ^ 2 * ρ0[1, 1] - @test ρ_rwa[1, 2] ≈ Gt * ρ0[1, 2] - end -end - -# Test Fermion-type HEOM Liouvillian superoperator matrix -@testset "M_Fermion" begin - L = M_Fermion(Hsys, tier, Fbath; verbose=false) - @test show(devnull, MIME("text/plain"), L) == nothing - @test size(L) == (1196, 1196) - @test L.N == 299 - @test nnz(L.data) == 21318 - L = addFermionDissipator(L, J) - @test nnz(L.data) == 22516 - 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 = [ - 0.49971864340781574 + 1.5063528845574697e-11im -0.00025004129095353573 + 0.00028356932981729176im; - -0.0002500413218393161 - 0.0002835693203755187im 0.5002813565929579 - 1.506436545359778e-11im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - L = M_Fermion(Hsys, tier, Fbath; threshold = 1e-8, verbose=false) - L = addFermionDissipator(L, J) - @test size(L) == (148, 148) - @test L.N == 37 - @test nnz(L.data) == 2054 - ados = SteadyState(L; verbose=false) - ρ2 = ados[1] - @test _is_Matrix_approx(ρ0, ρ2) - - L = M_Fermion(Hsys, tier, [Fbath, Fbath]; verbose=false) - @test size(L) == (9300, 9300) - @test L.N == 2325 - @test nnz(L.data) == 174338 - L = addFermionDissipator(L, J) - @test nnz(L.data) == 183640 - ados = SteadyState(L; verbose=false) - @test ados.dim == L.dim - @test length(ados) == L.N - ρ0 = ados[1] - @test getRho(ados) == ρ0 - ρ1 = [ - 0.4994229368103249 + 2.6656157051929377e-12im -0.0005219753638749304 + 0.0005685093274121244im; - -0.0005219753958601764 - 0.0005685092413099392im 0.5005770631903512 - 2.6376966158390854e-12im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - ## check exceptions - @test_throws BoundsError L[1, 9301] - @test_throws BoundsError L[1:9301, 9300] - @test_throws ErrorException ados[L.N + 1] - @test_throws ErrorException M_Fermion([0, 0], tier, Fbath; verbose=false) -end - -# Test Boson-Fermion-type HEOM Liouvillian superoperator matrix -@testset "M_Boson_Fermion" begin - # re-define the bath (make the matrix smaller) - λ = 0.1450 - W = 0.6464 - kT = 0.7414 - μ = 0.8787 - N = 3 - tierb = 2 - tierf = 2 - - Bbath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) - Fbath = Fermion_Lorentz_Pade(Q, λ, μ, W, kT, N) - - L = M_Boson_Fermion(Hsys, tierb, tierf, Bbath, Fbath; verbose=false) - @test show(devnull, MIME("text/plain"), L) == nothing - @test size(L) == (2220, 2220) - @test L.N == 555 - @test nnz(L.data) == 43368 - L = addBosonDissipator(L, J) - @test nnz(L.data) == 45590 - 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 = [ - 0.496709+4.88415e-12im -0.00324048+0.00286376im; - -0.00324048-0.00286376im 0.503291-4.91136e-12im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - L = M_Boson_Fermion(Hsys, tierb, tierf, [Bbath, Bbath], Fbath; verbose=false) - @test size(L) == (6660, 6660) - @test L.N == 1665 - @test nnz(L.data) == 139210 - L = addFermionDissipator(L, J) - @test nnz(L.data) == 145872 - ados = SteadyState(L; verbose=false) - @test ados.dim == L.dim - @test length(ados) == L.N - ρ0 = ados[1] - @test getRho(ados) == ρ0 - ρ1 = [ - 0.493774+6.27624e-13im -0.00536526+0.00651746im; - -0.00536526-0.00651746im 0.506226-6.15855e-13im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - L = M_Boson_Fermion(Hsys, tierb, tierf, Bbath, [Fbath, Fbath]; verbose=false) - @test size(L) == (8220, 8220) - @test L.N == 2055 - @test nnz(L.data) == 167108 - L = addBosonDissipator(L, J) - @test nnz(L.data) == 175330 - ados = SteadyState(L; verbose=false) - @test ados.dim == L.dim - @test length(ados) == L.N - ρ0 = ados[1] - @test getRho(ados) == ρ0 - ρ1 = [ - 0.496468-4.32253e-12im -0.00341484+0.00316445im; - -0.00341484-0.00316445im 0.503532+4.32574e-12im - ] - @test _is_Matrix_approx(ρ0, ρ1) - - ## check exceptions - @test_throws ErrorException ados[L.N + 1] - @test_throws ErrorException M_Boson_Fermion([0, 0], tierb, tierf, Bbath, Fbath; verbose=false) -end - -@testset "Hierarchy Dictionary" begin - Btier = 2 - Ftier = 2 - Nb = 3 - Nf = 3 - threshold = 1e-5 - - Γ = 0.0025 - Dα = 30 - Λ = 0.0025 - ωcα = 0.2 - μL = 0.5 - μR = -0.5 - kT = 0.025 - - Hsys = [ - 0 0 0 0; - 0 0.2 0 0; - 0 0 0.208 0.04; - 0 0 0.04 0.408 - ] - - ρ0 = [ - 1 0 0 0; - 0 0 0 0; - 0 0 0 0; - 0 0 0 0 - ] - - cop = [ - 0 1 0 0; - 1 0 0 0; - 0 0 0 1; - 0 0 1 0 - ] - - dop = [ - 0 0 1 0; - 0 0 0 1; - 0 0 0 0; - 0 0 0 0 - ] - - bbath = Boson_DrudeLorentz_Matsubara(cop, Λ, ωcα, kT, Nb) - fbath = [ - Fermion_Lorentz_Pade(dop, Γ, μL, Dα, kT, Nf), - Fermion_Lorentz_Pade(dop, Γ, μR, Dα, kT, Nf) - ] - - L = M_Boson_Fermion(Hsys, Btier, Ftier, bbath, fbath; threshold = threshold, verbose=false) - L = addTerminator(L, bbath) - - @test size(L) == (1696, 1696) - @test L.N == 106 - @test nnz(L.data) == 13536 - - ados = SteadyState(L; verbose=false) - @test Ic(ados, L, 1) ≈ 0.2883390125832726 - nvec_b, nvec_f = L.hierarchy.idx2nvec[1] - @test_throws ErrorException getIndexEnsemble(nvec_f, L.hierarchy.bosonPtr) - @test_throws ErrorException getIndexEnsemble(nvec_b, L.hierarchy.fermionPtr) -end - -@testset "Auxiliary density operators" begin - ados_b = ADOs(spzeros(Int64, 20), 5) - ados_f = ADOs(spzeros(Int64, 8), 2) - ados_bf = ADOs(spzeros(Int64, 40), 10) - @test show(devnull, MIME("text/plain"), ados_b) == nothing - @test show(devnull, MIME("text/plain"), ados_f) == nothing - @test show(devnull, MIME("text/plain"), ados_bf) == nothing - @test_throws ErrorException ADOs(zeros(8), 4) - - ρ_b = ados_b[:] - # check iteration - for (i, ado) in enumerate(ados_b) - @test ρ_b[i] == ado - end - - # expections for expect - ados_wrong = ADOs(spzeros(Int64, 18), 2) - @test_throws ErrorException Expect([0 0 0; 0 0 0; 0 0 0], ados_f) - @test_throws ErrorException Expect([0 0; 0 0], [ados_b, ados_wrong]) - @test_throws ErrorException Expect([0 0 0; 0 0 0; 0 0 0], [ados_b, ados_f]) -end \ No newline at end of file diff --git a/test/hierarchy_dictionary.jl b/test/hierarchy_dictionary.jl new file mode 100644 index 00000000..34615ae3 --- /dev/null +++ b/test/hierarchy_dictionary.jl @@ -0,0 +1,60 @@ +Btier = 2 +Ftier = 2 +Nb = 3 +Nf = 3 +threshold = 1e-5 + +Γ = 0.0025 +Dα = 30 +Λ = 0.0025 +ωcα = 0.2 +μL = 0.5 +μR = -0.5 +kT = 0.025 + +Hsys = [ + 0 0 0 0; + 0 0.2 0 0; + 0 0 0.208 0.04; + 0 0 0.04 0.408 +] + +ρ0 = [ + 1 0 0 0; + 0 0 0 0; + 0 0 0 0; + 0 0 0 0 +] + +cop = [ + 0 1 0 0; + 1 0 0 0; + 0 0 0 1; + 0 0 1 0 +] + +dop = [ + 0 0 1 0; + 0 0 0 1; + 0 0 0 0; + 0 0 0 0 +] + +bbath = Boson_DrudeLorentz_Matsubara(cop, Λ, ωcα, kT, Nb) +fbath = [ + Fermion_Lorentz_Pade(dop, Γ, μL, Dα, kT, Nf), + Fermion_Lorentz_Pade(dop, Γ, μR, Dα, kT, Nf) +] + +L = M_Boson_Fermion(Hsys, Btier, Ftier, bbath, fbath; threshold = threshold, verbose=false) +L = addTerminator(L, bbath) + +@test size(L) == (1696, 1696) +@test L.N == 106 +@test nnz(L.data) == 13536 + +ados = SteadyState(L; verbose=false) +@test Ic(ados, L, 1) ≈ 0.2883390125832726 +nvec_b, nvec_f = L.hierarchy.idx2nvec[1] +@test_throws ErrorException getIndexEnsemble(nvec_f, L.hierarchy.bosonPtr) +@test_throws ErrorException getIndexEnsemble(nvec_b, L.hierarchy.fermionPtr) \ No newline at end of file diff --git a/test/phys_analysis.jl b/test/phys_analysis.jl deleted file mode 100644 index abc054b9..00000000 --- a/test/phys_analysis.jl +++ /dev/null @@ -1,365 +0,0 @@ -# System Hamiltonian and initial state -Hsys = 0.25 * [1 0; 0 -1] + 0.5 * [0 1; 1 0] -ρ0 = [1 0; 0 0] - -# Bath properties: -λ = 0.1 -W = 0.5 -kT = 0.5 -N = 2 -Q = [1 0; 0 -1] # System-bath coupling operator -bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) - -# HEOM Liouvillian superoperator matrix -tier = 5 -L = M_Boson(Hsys, tier, bath; verbose=false) - -ados = SteadyState(L, ρ0; verbose=false, reltol=1e-2, abstol=1e-4) -ρs = getRho(ados) -@testset "Steady state" begin - O = [1 0.5; 0.5 1] - @test Expect(O, ados) ≈ real(tr(O * ρs)) - @test Expect(O, ados, take_real=false) ≈ tr(O * ρs) - - ρ1 = getRho(SteadyState(L; verbose=false)) - @test _is_Matrix_approx(ρ1, ρs) - - mat = spzeros(ComplexF64, 2, 2) - mat2 = spzeros(ComplexF64, 3, 3) - bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2) - @test_throws ErrorException SteadyState(M_Fermion(mat, 2, bathf, ODD; verbose=false)) - @test_throws ErrorException SteadyState(M_Fermion(mat, 2, bathf, ODD; verbose=false), mat) - @test_throws ErrorException SteadyState(L, mat2) - @test_throws ErrorException SteadyState(L, ADOs(zeros(8), 2)) - @test_throws ErrorException SteadyState(L, ADOs(ados.data, ados.N, ODD)) -end - -@testset "Time evolution" begin - bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) - L = M_Boson(Hsys, tier, bath; verbose=false) - ρ0 = [1 0; 0 0] - ρ_wrong = zeros(3, 3) - - Δt = 10 - steps = 10 - tlist = 0:Δt:(Δt * steps) - if isfile("evolution_p.jld2") - rm("evolution_p.jld2") - end - # using the method based on propagator - ados_list = evolution(L, ρ0, Δt, steps; verbose=false, filename="evolution_p") - ados_wrong1 = ADOs(zeros(8), 2) - ados_wrong2 = ADOs(zeros(32), 2) - ados_wrong3 = ADOs((ados_list[1]).data, (ados_list[1]).N, ODD) - ρ_list_p = getRho.(ados_list) - @test_throws ErrorException evolution(L, ρ0, Δt, steps; verbose=false, filename="evolution_p") - @test_throws ErrorException evolution(L, ρ_wrong, Δt, steps; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong1, Δt, steps; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong2, Δt, steps; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong3, Δt, steps; verbose=false) - - if isfile("evolution_o.jld2") - rm("evolution_o.jld2") - end - # using the method based on ODE solver - ρ_list_e = getRho.(evolution(L, ρ0, tlist; verbose=false, filename="evolution_o")) - @test_throws ErrorException evolution(L, ρ0, tlist; verbose=false, filename="evolution_o") - @test_throws ErrorException evolution(L, ρ_wrong, tlist; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong1, tlist; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong2, tlist; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong3, tlist; verbose=false) - - for i in 1:(steps + 1) - @test _is_Matrix_approx(ρ_list_p[i], ρ_list_e[i]) - end - @test _is_Matrix_approx(ρs, ρ_list_p[end]) - @test _is_Matrix_approx(ρs, ρ_list_e[end]) - - # time-dependent Hamiltonian - σz = [1 0; 0 -1] - P01 = [0 1; 0 0] - - H_sys = 0 * σz - ρ0 = [0.5 0.5; 0.5 0.5] - - bath = Boson_DrudeLorentz_Pade(σz, 0.0005, 0.005, 0.05, 3) - L = M_Boson(H_sys, 6, bath; verbose=false) - - function Ht(param, t) - amplitude, delay, integral = param - duration = integral / amplitude - period = duration + delay - - t = t % period - if t < duration - return amplitude * [0 1; 1 0] - else - return 0 * zeros(2, 2) - end - end - - tlist = 0:10:400 - if isfile("evolution_t.jld2") - rm("evolution_t.jld2") - end - fastDD_ados = evolution(L, ρ0, tlist, Ht, (0.50, 20, π/2); reltol=1e-12, abstol=1e-12, verbose=false, filename="evolution_t"); - @test_throws ErrorException evolution(L, ρ0, tlist, Ht, (0.50, 20, π/2); verbose=false, filename="evolution_t") - @test_throws ErrorException evolution(L, ρ_wrong, tlist, Ht; verbose=false) - fastBoFiN = [ - 0.4999999999999999, - 0.4972948770876402, - 0.48575366430956535, - 0.48623357850926063, - 0.4953126026707093, - 0.4956964571136953, - 0.4920357760578866, - 0.4791950252661603, - 0.4892461036798634, - 0.4928465660900118, - 0.4924652444726348, - 0.48681317030861687, - 0.4810852414658711, - 0.4895850463738865, - 0.4883063278114976, - 0.4891347504183582, - 0.4812529995344225, - 0.48397690318479186, - 0.48766721173648925, - 0.4863478506723415, - 0.4853757629933787, - 0.47619727575597826, - 0.4846116525810951, - 0.4838281128868905, - 0.4847126978506697, - 0.4809571657303493, - 0.47862209473563255, - 0.4832775448598134, - 0.4795481060217657, - 0.48232522269103417, - 0.47572829624493995, - 0.47970784844818903, - 0.48020739013048824, - 0.47927809262397914, - 0.47904939060891494, - 0.4728486106129371, - 0.47901940281020244, - 0.4755943953129941, - 0.47816314189739584, - 0.47479965067847246, - 0.47451220871416044 - ] - fastDD = Expect(P01, fastDD_ados) - @test typeof(fastDD) == Vector{Float64} - for i in 1:length(tlist) - @test fastDD[i] ≈ fastBoFiN[i] atol=1.0e-6 - end - - slowDD_ados = evolution(L, ρ0, tlist, Ht, (0.01, 20, π/2); reltol=1e-12, abstol=1e-12, verbose=false); - slowBoFiN = [ - 0.4999999999999999, - 0.4949826158957288, - 0.4808740017602084, - 0.4593383302633091, - 0.43252194834417496, - 0.402802685601788, - 0.3725210685299983, - 0.34375397372926303, - 0.3181509506463072, - 0.29684111346956915, - 0.2804084225914882, - 0.26892558062398203, - 0.26203207044640237, - 0.2590399640129812, - 0.25905156533518636, - 0.26107507001114955, - 0.2639313319708524, - 0.26368823254606255, - 0.25936367828335455, - 0.25399924039230276, - 0.2487144584444748, - 0.24360738475043517, - 0.23875449658047362, - 0.23419406109247376, - 0.2299206674661969, - 0.22588932666421335, - 0.22202631685357516, - 0.2182434870635807, - 0.21445292381955944, - 0.210579531547926, - 0.20656995199483505, - 0.20239715154725021, - 0.19806079107370503, - 0.19358407389379975, - 0.18856981068448642, - 0.18126256908840252, - 0.17215645949526584, - 0.16328689667616822, - 0.1552186906255167, - 0.14821389956195355, - 0.14240802098404504 - ] - slowDD = Expect(P01, slowDD_ados; take_real=false) - @test typeof(slowDD) == Vector{ComplexF64} - for i in 1:length(tlist) - @test slowDD[i] ≈ slowBoFiN[i] atol=1.0e-6 - end - - H_wrong1(param, t) = zeros(3, 3) - function H_wrong2(param, t) - if t == 0 - zeros(2, 2) - else - zeros(3, 3) - end - end - ados_wrong1 = ADOs(zeros(8), 2) - ados_wrong2 = ADOs(zeros(32), 2) - ados_wrong3 = ADOs((slowDD_ados[1]).data, (slowDD_ados[1]).N, ODD) - @test_throws ErrorException evolution(L, ρ0, tlist, H_wrong1; verbose=false) - @test_throws ErrorException evolution(L, ρ0, tlist, H_wrong2; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong1, tlist, Ht; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong2, tlist, Ht; verbose=false) - @test_throws ErrorException evolution(L, ados_wrong3, tlist, Ht; verbose=false) -end - -@testset "Power spectrum" begin - a = [0 1; 0 0] - - Hsys = a' * a - - λ = 1e-4 - W = 2e-1 - kT = 0.5 - N = 5 - bath = Boson_DrudeLorentz_Matsubara((a' + a), λ, W, kT, N) - - tier = 3 - L = M_Boson(Hsys, tier, bath; verbose=false) - L = addBosonDissipator(L, 1e-3 * a') - L = addTerminator(L, bath) - - ados_s = SteadyState(L; verbose=false) - ωlist = 0.9:0.01:1.1 - - if isfile("PSD.txt") - rm("PSD.txt") - end - psd1 = PowerSpectrum(L, ados_s, a, ωlist; verbose=false, filename="PSD") - psd2 = [ - 8.88036729e-04, - 1.06145358e-03, - 1.30081318e-03, - 1.64528197e-03, - 2.16857171e-03, - 3.02352743e-03, - 4.57224555e-03, - 7.85856353e-03, - 1.70441286e-02, - 6.49071303e-02, - 1.64934976e+02, - 6.60426108e-02, - 1.57205620e-02, - 6.74130995e-03, - 3.67130512e-03, - 2.27825666e-03, - 1.53535649e-03, - 1.09529424e-03, - 8.14605980e-04, - 6.25453434e-04, - 4.92451868e-04 - ] - for i in 1:length(ωlist) - @test psd1[i] ≈ psd2[i] atol=1.0e-6 - end - - mat = spzeros(ComplexF64, 2, 2) - mat2 = spzeros(ComplexF64, 3, 3) - bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2) - @test_throws ErrorException spectrum(L, ados_s, a, ωlist; verbose=false) - @test_throws ErrorException PowerSpectrum(L, ados_s, a, ωlist; verbose=false, filename="PSD") - @test_throws ErrorException PowerSpectrum(L, ados_s, mat2, ωlist; verbose=false) - @test_throws ErrorException PowerSpectrum(L, ADOs(zeros(8), 2), a, ωlist; verbose=false) - @test_throws ErrorException PowerSpectrum(L, ADOs(zeros(32), 2), a, ωlist; verbose=false) - @test_throws ErrorException PowerSpectrum(L, ADOs(ados_s.data, ados_s.N, ODD), a, ωlist; verbose=false) -end - -@testset "Density of states" begin - e = -5 - U = 10 - d_up = kron( [0 1; 0 0], [1 0; 0 1]) - d_dn = kron(-1 * [1 0; 0 -1], [0 1; 0 0]) - iden = kron( [1 0; 0 1], [1 0; 0 1]) - - H0 = e * (d_up' * d_up + d_dn' * d_dn) - H1 = U * (d_up' * d_up * d_dn' * d_dn) - Hsys = H0 + H1 - - λ = 1 - μ_l = 1 - μ_r = -1 - W = 10 - kT = 0.5 - N = 5 - fuL = Fermion_Lorentz_Pade(d_up, λ, μ_l, W, kT, N) - fdL = Fermion_Lorentz_Pade(d_dn, λ, μ_l, W, kT, N) - fuR = Fermion_Lorentz_Pade(d_up, λ, μ_r, W, kT, N) - fdR = Fermion_Lorentz_Pade(d_dn, λ, μ_r, W, kT, N) - - tier = 2 - Le = M_Fermion(Hsys, tier, [fuL, fdL, fuR, fdR]; verbose=false) - Lo = M_Fermion(Hsys, tier, [fuL, fdL, fuR, fdR], ODD; verbose=false) - - ados_s = SteadyState(Le; verbose=false) - ωlist = -20:2:20 - - if isfile("DOS.txt") - rm("DOS.txt") - end - dos1 = DensityOfStates(Lo, ados_s, d_up, ωlist; verbose=false, filename="DOS") - dos2 = PowerSpectrum(Lo, ados_s, d_up', ωlist, true; verbose=false) .+ PowerSpectrum(Lo, ados_s, d_up, ωlist, false; verbose=false) - dos3 = [ - 0.0007920428534358747, - 0.0012795202828027256, - 0.0022148985361417936, - 0.004203651086852703, - 0.009091831276694192, - 0.024145570990254425, - 0.09111929563034224, - 0.2984323182857298, - 0.1495897397559431, - 0.12433521300531171, - 0.17217519700362036, - 0.1243352130053117, - 0.14958973975594306, - 0.29843231828572964, - 0.09111929563034214, - 0.024145570990254432, - 0.009091831276694183, - 0.004203651086852702, - 0.0022148985361417923, - 0.0012795202828027236, - 0.0007920428534358735 - ] - for i in 1:length(ωlist) - @test dos1[i] ≈ dos3[i] atol=1.0e-10 - @test dos2[i] ≈ dos3[i] atol=1.0e-10 - end - - mat = spzeros(ComplexF64, 2, 2) - mat2 = spzeros(ComplexF64, 3, 3) - bathb = Boson_DrudeLorentz_Pade(mat, 1, 1, 1, 2) - @test_throws ErrorException spectrum(Lo, ados_s, d_up, ωlist; verbose=false) - @test_throws ErrorException DensityOfStates(Lo, ados_s, mat2, ωlist; verbose=false) - @test_throws ErrorException DensityOfStates(Lo, ADOs(zeros(8), 2), d_up, ωlist; verbose=false) - @test_throws ErrorException DensityOfStates(Lo, ADOs(zeros(32), 2), d_up, ωlist; verbose=false) - @test_throws ErrorException DensityOfStates(Lo, ADOs(ados_s.data, ados_s.N, ODD), d_up, ωlist; verbose=false) - @test_throws ErrorException DensityOfStates(M_Boson(mat, 2, bathb; verbose=false), mat, mat, [0]) - @test_throws ErrorException DensityOfStates(M_Fermion(mat, 2, fuL; verbose=false), mat, mat, [0]) -end - -# remove all the temporary files -rm("evolution_p.jld2") -rm("evolution_o.jld2") -rm("evolution_t.jld2") -rm("PSD.txt") -rm("DOS.txt") \ No newline at end of file diff --git a/test/power_spectrum.jl b/test/power_spectrum.jl new file mode 100644 index 00000000..4c5d44a9 --- /dev/null +++ b/test/power_spectrum.jl @@ -0,0 +1,64 @@ +a = [0 1; 0 0] + +Hsys = a' * a + +λ = 1e-4 +W = 2e-1 +kT = 0.5 +N = 5 +bath = Boson_DrudeLorentz_Matsubara((a' + a), λ, W, kT, N) + +tier = 3 +L = M_Boson(Hsys, tier, bath; verbose=false) +L = addBosonDissipator(L, 1e-3 * a') +L = addTerminator(L, bath) + +ados_s = SteadyState(L; verbose=false) +ωlist = 0.9:0.01:1.1 + +if isfile("PSD.txt") + rm("PSD.txt") +end +psd1 = PowerSpectrum(L, ados_s, a, ωlist; verbose=false, filename="PSD") +psd2 = [ + 8.88036729e-04, + 1.06145358e-03, + 1.30081318e-03, + 1.64528197e-03, + 2.16857171e-03, + 3.02352743e-03, + 4.57224555e-03, + 7.85856353e-03, + 1.70441286e-02, + 6.49071303e-02, + 1.64934976e+02, + 6.60426108e-02, + 1.57205620e-02, + 6.74130995e-03, + 3.67130512e-03, + 2.27825666e-03, + 1.53535649e-03, + 1.09529424e-03, + 8.14605980e-04, + 6.25453434e-04, + 4.92451868e-04 +] +for i in 1:length(ωlist) + @test psd1[i] ≈ psd2[i] atol=1.0e-6 +end + +mat = spzeros(ComplexF64, 2, 2) +mat2 = spzeros(ComplexF64, 3, 3) +a_even = HEOMSuperOp(a, EVEN, L) +a_odd = HEOMSuperOp(a, ODD, L) +bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2) +@test_throws ErrorException spectrum(L, ados_s, a, ωlist; verbose=false) +@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) +@test_throws ErrorException PowerSpectrum(L, ADOs(zeros(32), 2), a, ωlist; verbose=false) + +# remove all the temporary files +rm("PSD.txt") \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2a3e00fc..67989373 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,5 @@ -using Test, Pkg using HierarchicalEOM -using SparseArrays -using LinearAlgebra +using Test, Pkg, SparseArrays, LinearAlgebra const GROUP = get(ENV, "GROUP", "All") const HAS_EXTENSIONS = isdefined(Base, :get_extension) @@ -18,13 +16,57 @@ if GROUP == "All" || GROUP == "Core" include("bath.jl") end - @testset "Correlation functions" begin - include("corr_func.jl") + @testset "Bath correlation functions" begin + include("bath_corr_func.jl") end - include("heom_api.jl") + @testset "Auxiliary density operators" begin + include("ADOs.jl") + end + + @testset "HEOM superoperator" begin + include("HEOMSuperOp.jl") + end + + @testset "M_S" begin + include("M_S.jl") + end + + @testset "M_Boson" begin + include("M_Boson.jl") + end + + @testset "M_Boson (RWA)" begin + include("M_Boson_RWA.jl") + end + + @testset "M_Fermion" begin + include("M_Fermion.jl") + end - include("phys_analysis.jl") + @testset "M_Boson_Fermion" begin + include("M_Boson_Fermion.jl") + end + + @testset "Hierarchy Dictionary" begin + include("hierarchy_dictionary.jl") + end + + @testset "Stationary state" begin + include("stationary_state.jl") + end + + @testset "Time evolution" begin + include("time_evolution.jl") + end + + @testset "Power spectrum" begin + include("power_spectrum.jl") + end + + @testset "Density of states" begin + include("density_of_states.jl") + end end if GROUP == "HierarchicalEOM_CUDAExt" diff --git a/test/stationary_state.jl b/test/stationary_state.jl new file mode 100644 index 00000000..fe93e080 --- /dev/null +++ b/test/stationary_state.jl @@ -0,0 +1,34 @@ +# System Hamiltonian and initial state +Hsys = 0.25 * [1 0; 0 -1] + 0.5 * [0 1; 1 0] +ρ0 = [1 0; 0 0] + +# Bath properties: +λ = 0.1 +W = 0.5 +kT = 0.5 +N = 2 +Q = [1 0; 0 -1] # System-bath coupling operator +bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) + +# HEOM Liouvillian superoperator matrix +tier = 5 +L = M_Boson(Hsys, tier, bath; verbose=false) + +ados = SteadyState(L, ρ0; verbose=false, reltol=1e-2, abstol=1e-4) +ρs = getRho(ados) +O = [1 0.5; 0.5 1] +@test Expect(O, ados) ≈ real(tr(O * ρs)) +@test Expect(O, ados, take_real=false) ≈ tr(O * ρs) + +ρ1 = getRho(SteadyState(L; verbose=false)) +@test _is_Matrix_approx(ρ1, ρs) + +mat = spzeros(ComplexF64, 2, 2) +mat2 = spzeros(ComplexF64, 3, 3) +bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2) +@test_throws ErrorException SteadyState(M_Fermion(mat, 2, bathf, ODD; verbose=false)) +@test_throws ErrorException SteadyState(M_Fermion(mat, 2, bathf, ODD; verbose=false), mat) +@test_throws ErrorException SteadyState(L, mat2) +@test_throws ErrorException SteadyState(L, ADOs(zeros(8), 2)) +@test_throws ErrorException SteadyState(L, ADOs(ados.data, ados.N, ODD)) +@test_throws ErrorException SteadyState(L, HEOMSuperOp(Q, ODD, L) * ados) \ No newline at end of file diff --git a/test/time_evolution.jl b/test/time_evolution.jl new file mode 100644 index 00000000..3469ca4d --- /dev/null +++ b/test/time_evolution.jl @@ -0,0 +1,206 @@ +# System Hamiltonian and initial state +Hsys = 0.25 * [1 0; 0 -1] + 0.5 * [0 1; 1 0] +ρ0 = [1 0; 0 0] + +# Bath properties: +λ = 0.1 +W = 0.5 +kT = 0.5 +N = 2 +tier = 5 +Q = [1 0; 0 -1] # System-bath coupling operator + +bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N) + +L = M_Boson(Hsys, tier, bath; verbose=false) +ρ0 = [1 0; 0 0] +ρ_wrong = zeros(3, 3) + +Δt = 10 +steps = 10 +tlist = 0:Δt:(Δt * steps) +if isfile("evolution_p.jld2") + rm("evolution_p.jld2") +end +# using the method based on propagator +ados_list = evolution(L, ρ0, Δt, steps; verbose=false, filename="evolution_p") +ados_wrong1 = ADOs(zeros(8), 2) +ados_wrong2 = ADOs(zeros(32), 2) +ados_wrong3 = ADOs((ados_list[1]).data, (ados_list[1]).N, ODD) +ados_wrong4 = HEOMSuperOp(Q, ODD, ados_list[end]) * ados_list[end] +ρ_list_p = getRho.(ados_list) +@test_throws ErrorException evolution(L, ρ0, Δt, steps; verbose=false, filename="evolution_p") +@test_throws ErrorException evolution(L, ρ_wrong, Δt, steps; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong1, Δt, steps; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong2, Δt, steps; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong3, Δt, steps; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong4, Δt, steps; verbose=false) + +if isfile("evolution_o.jld2") + rm("evolution_o.jld2") +end +# using the method based on ODE solver +ρ_list_e = getRho.(evolution(L, ρ0, tlist; verbose=false, filename="evolution_o")) +@test_throws ErrorException evolution(L, ρ0, tlist; verbose=false, filename="evolution_o") +@test_throws ErrorException evolution(L, ρ_wrong, tlist; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong1, tlist; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong2, tlist; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong3, tlist; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong4, tlist; verbose=false) + +for i in 1:(steps + 1) + @test _is_Matrix_approx(ρ_list_p[i], ρ_list_e[i]) +end +@test _is_Matrix_approx(ρs, ρ_list_p[end]) +@test _is_Matrix_approx(ρs, ρ_list_e[end]) + +# time-dependent Hamiltonian +σz = [1 0; 0 -1] +P01 = [0 1; 0 0] + +H_sys = 0 * σz +ρ0 = [0.5 0.5; 0.5 0.5] + +bath = Boson_DrudeLorentz_Pade(σz, 0.0005, 0.005, 0.05, 3) +L = M_Boson(H_sys, 6, bath; verbose=false) + +function Ht(param, t) + amplitude, delay, integral = param + duration = integral / amplitude + period = duration + delay + + t = t % period + if t < duration + return amplitude * [0 1; 1 0] + else + return 0 * zeros(2, 2) + end +end + +tlist = 0:10:400 +if isfile("evolution_t.jld2") + rm("evolution_t.jld2") +end +fastDD_ados = evolution(L, ρ0, tlist, Ht, (0.50, 20, π/2); reltol=1e-12, abstol=1e-12, verbose=false, filename="evolution_t"); +@test_throws ErrorException evolution(L, ρ0, tlist, Ht, (0.50, 20, π/2); verbose=false, filename="evolution_t") +@test_throws ErrorException evolution(L, ρ_wrong, tlist, Ht; verbose=false) +fastBoFiN = [ + 0.4999999999999999, + 0.4972948770876402, + 0.48575366430956535, + 0.48623357850926063, + 0.4953126026707093, + 0.4956964571136953, + 0.4920357760578866, + 0.4791950252661603, + 0.4892461036798634, + 0.4928465660900118, + 0.4924652444726348, + 0.48681317030861687, + 0.4810852414658711, + 0.4895850463738865, + 0.4883063278114976, + 0.4891347504183582, + 0.4812529995344225, + 0.48397690318479186, + 0.48766721173648925, + 0.4863478506723415, + 0.4853757629933787, + 0.47619727575597826, + 0.4846116525810951, + 0.4838281128868905, + 0.4847126978506697, + 0.4809571657303493, + 0.47862209473563255, + 0.4832775448598134, + 0.4795481060217657, + 0.48232522269103417, + 0.47572829624493995, + 0.47970784844818903, + 0.48020739013048824, + 0.47927809262397914, + 0.47904939060891494, + 0.4728486106129371, + 0.47901940281020244, + 0.4755943953129941, + 0.47816314189739584, + 0.47479965067847246, + 0.47451220871416044 +] +fastDD = Expect(P01, fastDD_ados) +@test typeof(fastDD) == Vector{Float64} +for i in 1:length(tlist) + @test fastDD[i] ≈ fastBoFiN[i] atol=1.0e-6 +end + +slowDD_ados = evolution(L, ρ0, tlist, Ht, (0.01, 20, π/2); reltol=1e-12, abstol=1e-12, verbose=false); +slowBoFiN = [ + 0.4999999999999999, + 0.4949826158957288, + 0.4808740017602084, + 0.4593383302633091, + 0.43252194834417496, + 0.402802685601788, + 0.3725210685299983, + 0.34375397372926303, + 0.3181509506463072, + 0.29684111346956915, + 0.2804084225914882, + 0.26892558062398203, + 0.26203207044640237, + 0.2590399640129812, + 0.25905156533518636, + 0.26107507001114955, + 0.2639313319708524, + 0.26368823254606255, + 0.25936367828335455, + 0.25399924039230276, + 0.2487144584444748, + 0.24360738475043517, + 0.23875449658047362, + 0.23419406109247376, + 0.2299206674661969, + 0.22588932666421335, + 0.22202631685357516, + 0.2182434870635807, + 0.21445292381955944, + 0.210579531547926, + 0.20656995199483505, + 0.20239715154725021, + 0.19806079107370503, + 0.19358407389379975, + 0.18856981068448642, + 0.18126256908840252, + 0.17215645949526584, + 0.16328689667616822, + 0.1552186906255167, + 0.14821389956195355, + 0.14240802098404504 +] +slowDD = Expect(P01, slowDD_ados; take_real=false) +@test typeof(slowDD) == Vector{ComplexF64} +for i in 1:length(tlist) + @test slowDD[i] ≈ slowBoFiN[i] atol=1.0e-6 +end + +H_wrong1(param, t) = zeros(3, 3) +function H_wrong2(param, t) + if t == 0 + zeros(2, 2) + else + zeros(3, 3) + end +end +ados_wrong1 = ADOs(zeros(8), 2) +ados_wrong2 = ADOs(zeros(32), 2) +ados_wrong3 = ADOs((slowDD_ados[1]).data, (slowDD_ados[1]).N, ODD) +@test_throws ErrorException evolution(L, ρ0, tlist, H_wrong1; verbose=false) +@test_throws ErrorException evolution(L, ρ0, tlist, H_wrong2; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong1, tlist, Ht; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong2, tlist, Ht; verbose=false) +@test_throws ErrorException evolution(L, ados_wrong3, tlist, Ht; verbose=false) + +# remove all the temporary files +rm("evolution_p.jld2") +rm("evolution_o.jld2") +rm("evolution_t.jld2") \ No newline at end of file