From d4d48e7390dcab4521c0de1db3467b6459e52267 Mon Sep 17 00:00:00 2001 From: Yi-Te Huang Date: Fri, 22 Nov 2024 11:16:02 +0900 Subject: [PATCH] fix runtests for `HEOMSuperOp` and spectra functions --- src/HeomBase.jl | 4 ++-- src/density_of_states.jl | 2 +- src/deprecated.jl | 13 +++++++++++++ src/power_spectrum.jl | 5 ++--- test/HEOMSuperOp.jl | 14 ++++++++------ test/density_of_states.jl | 2 +- test/power_spectrum.jl | 4 ++-- test/stationary_state.jl | 2 +- test/time_evolution.jl | 2 +- 9 files changed, 31 insertions(+), 17 deletions(-) diff --git a/src/HeomBase.jl b/src/HeomBase.jl index 99861842..770e54a8 100644 --- a/src/HeomBase.jl +++ b/src/HeomBase.jl @@ -2,8 +2,8 @@ export AbstractHEOMLSMatrix abstract type AbstractHEOMLSMatrix{T} end -QuantumToolbox._FType(::AbstractHEOMLSMatrix{<:AbstractArray{T}}) where {T<:Number} = _FType(T) -QuantumToolbox._CType(::AbstractHEOMLSMatrix{<:AbstractArray{T}}) where {T<:Number} = _CType(T) +QuantumToolbox._FType(M::AbstractHEOMLSMatrix) = _FType(eltype(M)) +QuantumToolbox._CType(M::AbstractHEOMLSMatrix) = _CType(eltype(M)) _get_SciML_matrix_wrapper(M::AbstractArray) = QuantumToolbox.get_typename_wrapper(M){eltype(M)} _get_SciML_matrix_wrapper(M::MatrixOperator) = _get_SciML_matrix_wrapper(M.A) diff --git a/src/density_of_states.jl b/src/density_of_states.jl index 7660d9c9..b3515e37 100644 --- a/src/density_of_states.jl +++ b/src/density_of_states.jl @@ -54,7 +54,7 @@ Calculate density of states for the fermionic system in frequency domain. # Handle d_op MType = _get_SciML_matrix_wrapper(M) _tr = transpose(_Tr(M)) - Id_sys = I(prod(M.dims)) + Id_sys = I(prod(d_op.dims)) Id_HEOM = I(M.N) d_normal = HEOMSuperOp(spre(d_op, Id_sys), ODD, M; Id_cache = Id_HEOM) d_dagger = HEOMSuperOp(spre(d_op', Id_sys), ODD, M; Id_cache = Id_HEOM) diff --git a/src/deprecated.jl b/src/deprecated.jl index f1e7f6ee..fd46f42f 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -19,3 +19,16 @@ export evolution, SteadyState SteadyState(args...; kwargs...) = error("`SteadyState` has been deprecated, please use `steadystate` instead.") evolution(args...; kwargs...) = error("`evolution` has been deprecated, please use `HEOMsolve` instead.") + +_HEOMSuperOp_deprecated_method_error() = error("Specifying `mul_basis = \"L\", \"R\", or \"LR\"` to `HEOMSuperOp` has been deprecated, try to construct system SuperOperator manually by using `spre`, `spost`, or `sprepost`.") + +HEOMSuperOp(op, opParity::AbstractParity, dims::SVector, N::Int, mul_basis::AbstractString; Id_cache = I(N)) = _HEOMSuperOp_deprecated_method_error() +HEOMSuperOp( + op, + opParity::AbstractParity, + refHEOMLS::AbstractHEOMLSMatrix, + mul_basis::AbstractString; + Id_cache = I(refHEOMLS.N), +) = HEOMSuperOp(op, opParity, refHEOMLS.dims, refHEOMLS.N, mul_basis; Id_cache = Id_cache) +HEOMSuperOp(op, opParity::AbstractParity, refADOs::ADOs, mul_basis::AbstractString; Id_cache = I(refADOs.N)) = + HEOMSuperOp(op, opParity, refADOs.dims, refADOs.N, mul_basis; Id_cache = Id_cache) diff --git a/src/power_spectrum.jl b/src/power_spectrum.jl index 29e6a277..186d4ed5 100644 --- a/src/power_spectrum.jl +++ b/src/power_spectrum.jl @@ -85,7 +85,6 @@ remember to set the parameters: end _check_sys_dim_and_ADOs_num(M, ados) - Id_sys = I(prod(M.dims)) Id_HEOM = I(M.N) # Handle P_op @@ -93,7 +92,7 @@ remember to set the parameters: _check_sys_dim_and_ADOs_num(M, P_op) _P = P_op else - _P = HEOMSuperOp(spre(P_op, Id_sys), EVEN, M; Id_cache = Id_HEOM) + _P = HEOMSuperOp(spre(P_op), EVEN, M; Id_cache = Id_HEOM) end MType = _get_SciML_matrix_wrapper(M) _tr_P = transpose(_Tr(M)) * MType(_P).data @@ -109,7 +108,7 @@ remember to set the parameters: else new_parity = !ados.parity end - _Q_ados = HEOMSuperOp(spre(Q_op, Id_sys), new_parity, M; Id_cache = Id_HEOM) * ados + _Q_ados = HEOMSuperOp(spre(Q_op), new_parity, M; Id_cache = Id_HEOM) * ados end b = _HandleVectorType(M, _Q_ados.data) diff --git a/test/HEOMSuperOp.jl b/test/HEOMSuperOp.jl index 34e13c4a..c2863fba 100644 --- a/test/HEOMSuperOp.jl +++ b/test/HEOMSuperOp.jl @@ -45,7 +45,7 @@ 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") + J_me = HEOMSuperOp(γ_eR * spre(d), ODD, M_me) * HEOMSuperOp(γ_eR * spost(d'), ODD, M_me) @test size(M_me) == size(J_me) @test size(M_me, 1) == size(J_me, 1) @test eltype(M_me) == eltype(J_me) @@ -53,7 +53,7 @@ @test J_me[4, 1] ≈ 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 + (HEOMSuperOp(√(2) * γ_eR * spre(d), ODD, M_me) * HEOMSuperOp(√(2) * γ_eR * spost(d'), ODD, M_me)).data @test show(devnull, MIME("text/plain"), J_me) == nothing # create HEOM Liouvuillian superoperator without quantum jump @@ -77,7 +77,7 @@ 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") + J_heom = HEOMSuperOp(γ_eR * spre(d), ODD, M_heom) * HEOMSuperOp(γ_eR * spost(d'), ODD, M_heom) # create HEOM Liouvuillian superoperator without quantum jump M1 = M_heom - J_heom @@ -93,9 +93,9 @@ end ## test for error - J_wrong1 = HEOMSuperOp(γ_eR * d, ODD, M_me, "L") - J_wrong2 = HEOMSuperOp(γ_eR * d, EVEN, M_heom, "L") - @test_throws ErrorException HEOMSuperOp(d, EVEN, M_heom, "LR") + J_wrong1 = HEOMSuperOp(γ_eR * spre(d), ODD, M_me) + J_wrong2 = HEOMSuperOp(γ_eR * spre(d), EVEN, M_heom) + @test_throws ErrorException HEOMSuperOp(sprepost(d), EVEN, M_heom) @test_throws ErrorException J_me * J_wrong2 @test_throws ErrorException J_me + J_wrong1 @test_throws ErrorException J_me + J_wrong2 @@ -105,4 +105,6 @@ @test_throws ErrorException J_me - J_wrong2 @test_throws ErrorException M_me - J_wrong1 @test_throws ErrorException M_me - J_wrong2 + @test_throws ErrorException HEOMSuperOp(d, ODD, M_heom, "L") + @test_throws ErrorException HEOMSuperOp(d, ODD, ados_s1, "L") end diff --git a/test/density_of_states.jl b/test/density_of_states.jl index 62c2d050..9f39cc65 100644 --- a/test/density_of_states.jl +++ b/test/density_of_states.jl @@ -30,7 +30,7 @@ if isfile("DOS.txt") rm("DOS.txt") end - d_up_normal = HEOMSuperOp(d_up, ODD, Le) + d_up_normal = HEOMSuperOp(spre(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) .+ diff --git a/test/power_spectrum.jl b/test/power_spectrum.jl index d103371b..260628c1 100644 --- a/test/power_spectrum.jl +++ b/test/power_spectrum.jl @@ -51,8 +51,8 @@ mat = Qobj(spzeros(ComplexF64, 2, 2)) mat2 = Qobj(spzeros(ComplexF64, 3, 3)) - a_even = HEOMSuperOp(a, EVEN, L) - a_odd = HEOMSuperOp(a, ODD, L) + a_even = HEOMSuperOp(spre(a), EVEN, L) + a_odd = HEOMSuperOp(spre(a), ODD, L) bathf = Fermion_Lorentz_Pade(mat, 1, 1, 1, 1, 2) @test_throws ErrorException PowerSpectrum(L, ados_s, a, ωlist; verbose = false, filename = "PSD") @test_throws ErrorException PowerSpectrum(L, ados_s, a_even, a_odd, ωlist; verbose = false) diff --git a/test/stationary_state.jl b/test/stationary_state.jl index e3e9b857..1d5a7de4 100644 --- a/test/stationary_state.jl +++ b/test/stationary_state.jl @@ -37,5 +37,5 @@ @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(d, ODD, L) * ados) + @test_throws ErrorException steadystate(L, HEOMSuperOp(spre(d), ODD, L) * ados) end diff --git a/test/time_evolution.jl b/test/time_evolution.jl index fda52665..217363d4 100644 --- a/test/time_evolution.jl +++ b/test/time_evolution.jl @@ -32,7 +32,7 @@ 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] + ados_wrong4 = HEOMSuperOp(spre(Q), ODD, ados_list[end]) * ados_list[end] ρ_list_p = getRho.(ados_list) @test show(devnull, MIME("text/plain"), sol_p) === nothing @test length(sol_p.ados) == 1