add QuPhysExt and corresponding runtests

Project.toml

@@ -22,10 +22,12 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c"
+QuPhys = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
 
 [extensions]
 HierarchicalEOM_CUDAExt = "CUDA"
 HierarchicalEOM_QOExt = "QuantumOptics"
+HierarchicalEOM_QuPhysExt = "QuPhys"
 
 [compat]
 Crayons = "4.1"

ext/HierarchicalEOM_QuPhysExt.jl

@@ -0,0 +1,34 @@
+module HierarchicalEOM_QuPhysExt
+
+import QuPhys: QuantumObject
+using HierarchicalEOM
+import HierarchicalEOM.HeomBase: HandleMatrixType
+
+@doc raw"""
+    QuantumObject(data::AbstractMatrix, refOP::QuantumObject)
+Return the operator under the type of `QuPhys.QuantumObject` from a given matrix and reference operator (type, dims).
+
+# Parameters
+- `data` : The matrix-type operator.
+- `refOP<:QuPhys.QuantumObject` : the reference operator from `QuPhys`.
+"""
+QuantumObject(data::AbstractMatrix, refOP::QuantumObject) = QuantumObject(data, refOP.type, refOP.dims)
+
+function HandleMatrixType(M::QuantumObject, dim::Int=0, MatrixName::String="")
+    if dim > 0
+        if size(M.data) == (dim, dim)
+            return copy(M.data)
+        else
+            error("The size of matrix $(MatrixName) should be: ($(dim), $(dim)).")
+        end
+    elseif dim == 0
+        N1, N2 = size(M.data)
+        if N1 == N2
+            return copy(M.data)
+        else
+            error("The size of matrix $(MatrixName) should be squared matrix.")
+        end
+    end
+end
+
+end

test/QuPhys/Project.toml

@@ -0,0 +1,3 @@
+[deps]
+HierarchicalEOM = "a62dbcb7-80f5-4d31-9a88-8b19fd92b128"
+QuPhys = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"

test/QuPhys/QuPhysExt.jl

@@ -0,0 +1,74 @@
+import QuPhys
+
+# System Hamiltonian and initial state
+σx = QuPhys.sigmax()
+σz = QuPhys.sigmaz()
+σm = QuPhys.sigmam()
+I2 = QuPhys.eye(2)
+
+Hsys = 0.25 * σz + 0.5 * σx
+ρ0 = QuPhys.ket2dm(QuPhys.basis(2, 0))
+
+λ  = 0.1
+W  = 0.5
+kT = 0.5
+N  = 2
+Q = σz  # System-bath coupling operator
+tier = 2
+bath = Boson_DrudeLorentz_Pade(Q, λ, W, kT, N)
+L = M_Boson(Hsys, tier, bath; verbose=false)
+
+# SteadyState
+ados = SteadyState(L, ρ0; verbose=false, reltol=1e-2, abstol=1e-4)
+
+# Evolution
+tier = 5
+O = I2 + 0.5 * σx
+Δt    = 10
+steps = 1
+tlist = 0:Δt:(Δt * steps)
+ados_list = evolution(L, ρ0, Δt, steps; verbose=false)
+ados_list = evolution(L, ρ0, tlist; verbose=false)
+@test Expect(O, ados_list[end]) ≈ QuPhys.expect(O, QuPhys.QuantumObject(ados_list[end][1], Hsys))
+
+# Power spectral density
+a = QuPhys.sigmam()
+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]
+@test PowerSpectrum(L, ados_s, a, ωlist; verbose=false)[1] ≈ 0.0008880367286438112
+
+# Density of states
+e = -5
+U = 10
+d_up = kron( σm, I2)
+d_dn = kron(-σz, σm)
+iden = kron( I2, I2)
+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 = [0]
+@test DensityOfStates(Lo, ados_s, d_up, ωlist; verbose=false)[1] ≈ 0.17217519700362002

test/QuantumOptics/QOExt.jl

@@ -1,14 +1,14 @@
-using QuantumOptics
+import QuantumOptics
 
 # System Hamiltonian and initial state
-basis = SpinBasis(1//2)
-σx = sigmax(basis)
-σz = sigmaz(basis)
-σm = sigmam(basis)
-I2 = identityoperator(basis)
+basis = QuantumOptics.SpinBasis(1//2)
+σx = QuantumOptics.sigmax(basis)
+σz = QuantumOptics.sigmaz(basis)
+σm = QuantumOptics.sigmam(basis)
+I2 = QuantumOptics.identityoperator(basis)
 
 Hsys = 0.25 * σz + 0.5 * σx
-ρ0 = dm(Ket(basis, [1, 0]))
+ρ0 = QuantumOptics.dm(QuantumOptics.Ket(basis, [1, 0]))
 
 λ  = 0.1
 W  = 0.5
@@ -30,10 +30,10 @@ steps = 1
 tlist = 0:Δt:(Δt * steps)
 ados_list = evolution(L, ρ0, Δt, steps; verbose=false)
 ados_list = evolution(L, ρ0, tlist; verbose=false)
-@test Expect(O, ados_list[end]) ≈ real(tr(O * Operator(ados_list[end][1], Hsys)))
+@test Expect(O, ados_list[end]) ≈ QuantumOptics.expect(O, QuantumOptics.Operator(ados_list[end][1], Hsys))
 
 # Power spectral density
-a = sigmam(basis)
+a = QuantumOptics.sigmam(basis)
 Hsys = a' * a
 λ  = 1e-4
 W  = 2e-1
@@ -51,9 +51,9 @@ ados_s = SteadyState(L; verbose=false)
 # Density of states
 e = -5
 U = 10
-d_up =  σm ⊗ I2
-d_dn = -σz ⊗ σm
-iden =  I2 ⊗ I2
+d_up = QuantumOptics.tensor( σm, I2)
+d_dn = QuantumOptics.tensor(-σz, σm)
+iden = QuantumOptics.tensor( I2, I2)
 H0 = e * (d_up' * d_up + d_dn' * d_dn)
 H1 = U * (d_up' * d_up * d_dn' * d_dn)
 Hsys = H0 + H1

test/runtests.jl

@@ -84,4 +84,14 @@ if (GROUP == "All") || (GROUP == "HierarchicalEOM_QOExt")
     @testset "QuantumOptics Extension" begin
         include("QuantumOptics/QOExt.jl")
     end
+end
+
+if (GROUP == "All") || (GROUP == "HierarchicalEOM_QuPhysExt")
+    Pkg.activate("QuPhys")
+    Pkg.add(url="https://github.com/albertomercurio/QuPhys.jl.git")
+    Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
+    Pkg.instantiate()
+    @testset "QuPhys Extension" begin
+        include("QuPhys/QuPhysExt.jl")
+    end
 end