Excitation number restricted (ENR) state space implementation (#500)

Author: TendonFFF

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Implement `EnrSpace` and corresponding functionality. ([#500])
 - Check for orthogonality breakdown in `Lanczos` solver for `spectrum`. ([#501])
 
 ## [v0.32.1]
@@ -255,4 +256,5 @@ Release date: 2024-11-13
 [#487]: https://github.com/qutip/QuantumToolbox.jl/issues/487
 [#489]: https://github.com/qutip/QuantumToolbox.jl/issues/489
 [#494]: https://github.com/qutip/QuantumToolbox.jl/issues/494
+[#500]: https://github.com/qutip/QuantumToolbox.jl/issues/500
 [#501]: https://github.com/qutip/QuantumToolbox.jl/issues/501

docs/src/resources/api.md

@@ -18,6 +18,7 @@ end
 
 ```@docs
 Space
+EnrSpace
 Dimensions
 GeneralDimensions
 AbstractQuantumObject
@@ -122,6 +123,8 @@ coherent_dm
 thermal_dm
 maximally_mixed_dm
 rand_dm
+enr_fock
+enr_thermal_dm
 spin_state
 spin_coherent
 bell_state
@@ -152,6 +155,8 @@ QuantumToolbox.momentum
 phase
 fdestroy
 fcreate
+enr_destroy
+enr_identity
 tunneling
 qft
 eye
@@ -312,6 +317,7 @@ PhysicalConstants
 convert_unit
 row_major_reshape
 meshgrid
+enr_state_dictionaries
 ```
 
 ## [Visualization](@id doc-API:Visualization)

docs/src/users_guide/QuantumObject/QuantumObject.md

@@ -66,14 +66,16 @@ Manually specifying the data for each quantum object is inefficient. Even more s
 
 ### States
 - [`zero_ket`](@ref): zero ket vector
-- [`fock`](@ref) or [`basis`](@ref): fock state ket vector
-- [`fock_dm`](@ref): density matrix of a fock state
+- [`fock`](@ref) or [`basis`](@ref): Fock state ket vector
+- [`fock_dm`](@ref): density matrix of a Fock state
 - [`coherent`](@ref): coherent state ket vector 
 - [`rand_ket`](@ref): random ket vector
 - [`coherent_dm`](@ref): density matrix of a coherent state
 - [`thermal_dm`](@ref): density matrix of a thermal state
 - [`maximally_mixed_dm`](@ref): density matrix of a maximally mixed state
 - [`rand_dm`](@ref): random density matrix
+- [`enr_fock`](@ref): Fock state in the excitation number restricted (ENR) space
+- [`enr_thermal_dm`](@ref): thermal state in the excitation number restricted (ENR) space
 - [`spin_state`](@ref): spin state
 - [`spin_coherent`](@ref): coherent spin state
 - [`bell_state`](@ref): Bell state
@@ -108,6 +110,8 @@ Manually specifying the data for each quantum object is inefficient. Even more s
 - [`spin_J_set`](@ref): a set of Spin-`j` operators ``(S_x, S_y, S_z)``
 - [`fdestroy`](@ref): fermion destruction operator
 - [`fcreate`](@ref): fermion creation operator
+- [`enr_destroy`](@ref): destruction operator in the excitation number restricted (ENR) space
+- [`enr_identity`](@ref): identity operator in the excitation number restricted (ENR) space
 - [`commutator`](@ref): commutator or anti-commutator
 - [`tunneling`](@ref): tunneling operator
 - [`qft`](@ref): discrete quantum Fourier transform matrix

src/QuantumToolbox.jl

@@ -84,6 +84,7 @@ include("linear_maps.jl")
 
 # Quantum Object
 include("qobj/space.jl")
+include("qobj/energy_restricted.jl")
 include("qobj/dimensions.jl")
 include("qobj/quantum_object_base.jl")
 include("qobj/quantum_object.jl")

src/negativity.jl

@@ -71,18 +71,20 @@ Return the partial transpose of a density matrix ``\rho``, where `mask` is an ar
 - `ρ_pt::QuantumObject`: The density matrix with the selected subsystems transposed.
 """
 function partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
-    if length(mask) != length(ρ.dimensions)
+    any(s -> s isa EnrSpace, ρ.dimensions.to) && throw(ArgumentError("partial_transpose does not support EnrSpace"))
+
+    (length(mask) != length(ρ.dimensions)) &&
         throw(ArgumentError("The length of \`mask\` should be equal to the length of \`ρ.dims\`."))
-    end
-    return _partial_transpose(ρ, mask)
-end
 
-# for dense matrices
-function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     isa(ρ.dimensions, GeneralDimensions) &&
         (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
         throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
 
+    return _partial_transpose(ρ, mask)
+end
+
+# for dense matrices
+function _partial_transpose(ρ::QuantumObject{Operator}, mask::Vector{Bool})
     mask2 = [1 + Int(i) for i in mask]
     # mask2 has elements with values equal to 1 or 2
     #   1 - the subsystem don't need to be transposed
@@ -107,10 +109,6 @@ function _partial_transpose(
     ρ::QuantumObject{Operator,DimsType,<:AbstractSparseArray},
     mask::Vector{Bool},
 ) where {DimsType<:AbstractDimensions}
-    isa(ρ.dimensions, GeneralDimensions) &&
-        (get_dimensions_to(ρ) != get_dimensions_from(ρ)) &&
-        throw(ArgumentError("Invalid partial transpose for dims = $(_get_dims_string(ρ.dimensions))"))
-
     M, N = size(ρ)
     dimsTuple = Tuple(dimensions_to_dims(get_dimensions_to(ρ)))
     colptr = ρ.data.colptr

src/qobj/arithmetic_and_attributes.jl

@@ -505,6 +505,8 @@ Quantum Object:   type=Operator()   dims=[2]   size=(2, 2)   ishermitian=true
 ```
 """
 function ptrace(QO::QuantumObject{Ket}, sel::Union{AbstractVector{Int},Tuple})
+    any(s -> s isa EnrSpace, QO.dimensions.to) && throw(ArgumentError("ptrace does not support EnrSpace"))
+
     _non_static_array_warning("sel", sel)
 
     if length(sel) == 0 # return full trace for empty sel
@@ -527,6 +529,8 @@ end
 ptrace(QO::QuantumObject{Bra}, sel::Union{AbstractVector{Int},Tuple}) = ptrace(QO', sel)
 
 function ptrace(QO::QuantumObject{Operator}, sel::Union{AbstractVector{Int},Tuple})
+    any(s -> s isa EnrSpace, QO.dimensions.to) && throw(ArgumentError("ptrace does not support EnrSpace"))
+
     # TODO: support for special cases when some of the subsystems have same `to` and `from` space
     isa(QO.dimensions, GeneralDimensions) &&
         (get_dimensions_to(QO) != get_dimensions_from(QO)) &&
@@ -714,6 +718,8 @@ function SparseArrays.permute(
     A::QuantumObject{ObjType},
     order::Union{AbstractVector{Int},Tuple},
 ) where {ObjType<:Union{Ket,Bra,Operator}}
+    any(s -> s isa EnrSpace, A.dimensions.to) && throw(ArgumentError("permute does not support EnrSpace"))
+
     (length(order) != length(A.dimensions)) &&
         throw(ArgumentError("The order list must have the same length as the number of subsystems (A.dims)"))
 

src/qobj/dimensions.jl

@@ -17,7 +17,7 @@ struct Dimensions{N,T<:Tuple} <: AbstractDimensions{N,N}
     to::T
 
     # make sure the elements in the tuple are all AbstractSpace
-    Dimensions(to::NTuple{N,T}) where {N,T<:AbstractSpace} = new{N,typeof(to)}(to)
+    Dimensions(to::NTuple{N,AbstractSpace}) where {N} = new{N,typeof(to)}(to)
 end
 function Dimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Integer,N}
     _non_static_array_warning("dims", dims)
@@ -43,12 +43,11 @@ Dimensions(dims::Any) = throw(
 A structure that describes the left-hand side (`to`) and right-hand side (`from`) Hilbert [`Space`](@ref) of an [`Operator`](@ref).
 """
 struct GeneralDimensions{M,N,T1<:Tuple,T2<:Tuple} <: AbstractDimensions{M,N}
-    # note that the number `N` should be the same for both `to` and `from`
     to::T1   # space acting on the left
     from::T2 # space acting on the right
 
     # make sure the elements in the tuple are all AbstractSpace
-    GeneralDimensions(to::NTuple{M,T1}, from::NTuple{N,T2}) where {M,N,T1<:AbstractSpace,T2<:AbstractSpace} =
+    GeneralDimensions(to::NTuple{M,AbstractSpace}, from::NTuple{N,AbstractSpace}) where {M,N} =
         new{M,N,typeof(to),typeof(from)}(to, from)
 end
 function GeneralDimensions(dims::Union{AbstractVector{T},NTuple{N,T}}) where {T<:Union{AbstractVector,NTuple},N}
@@ -98,3 +97,8 @@ function _get_dims_string(dimensions::GeneralDimensions)
     return "[$(string(dims[1])), $(string(dims[2]))]"
 end
 _get_dims_string(::Nothing) = "nothing" # for EigsolveResult.dimensions = nothing
+
+Base.:(==)(dim1::Dimensions, dim2::Dimensions) = dim1.to == dim2.to
+Base.:(==)(dim1::GeneralDimensions, dim2::GeneralDimensions) = (dim1.to == dim2.to) && (dim1.from == dim2.from)
+Base.:(==)(dim1::Dimensions, dim2::GeneralDimensions) = false
+Base.:(==)(dim1::GeneralDimensions, dim2::Dimensions) = false