From ccd349d86212e58d6b723e16d06117e0129e405e Mon Sep 17 00:00:00 2001
From: lch <lchopn@gmail.com>
Date: Sat, 28 Jun 2025 02:12:11 +0800
Subject: [PATCH 1/6] Unified ModeMap Type (#330)

1. Renamed the type OccupiedModeMap to ModeMap and replaced all occurrences.
2. Added the occupied_mode_map function as a constructor for ModeMap.
3. Added the unoccupied_mode_map function as a constructor for FermiFS
and corresponding iterator type FermiUnoccupiedModes.
---
 src/BitStringAddresses/BitStringAddresses.jl  |   4 +-
 src/BitStringAddresses/bitstring.jl           |  39 ++++++
 src/BitStringAddresses/fermifs.jl             |  44 +++++++
 src/BitStringAddresses/fockaddress.jl         | 120 +++++++++++++-----
 src/BitStringAddresses/occupationnumberfs.jl  |   2 +-
 src/Hamiltonians/ExtendedHubbardMom1D.jl      |   8 +-
 src/Hamiltonians/ExtendedHubbardReal1D.jl     |   4 +-
 src/Hamiltonians/G2MomCorrelator.jl           |   2 +-
 .../HOCartesianCentralImpurity.jl             |   8 +-
 .../HOCartesianContactInteractions.jl         |   4 +-
 .../HOCartesianEnergyConservedPerDim.jl       |  22 ++--
 src/Hamiltonians/HubbardMom1D.jl              |  26 ++--
 src/Hamiltonians/HubbardMom1DEP.jl            |  14 +-
 src/Hamiltonians/Momentum.jl                  |   2 +-
 src/Hamiltonians/Transcorrelated1D.jl         |   8 +-
 src/Hamiltonians/angular_momentum.jl          |   2 +-
 src/Hamiltonians/excitations.jl               |  26 ++--
 src/Hamiltonians/ho-cart-tools.jl             |   4 +-
 test/BitStringAddresses.jl                    |   2 +-
 test/Hamiltonians.jl                          |  20 +--
 20 files changed, 247 insertions(+), 114 deletions(-)

diff --git a/src/BitStringAddresses/BitStringAddresses.jl b/src/BitStringAddresses/BitStringAddresses.jl
index 0ece599f6..3ec8dbb8c 100644
--- a/src/BitStringAddresses/BitStringAddresses.jl
+++ b/src/BitStringAddresses/BitStringAddresses.jl
@@ -18,8 +18,8 @@ export OccupationNumberFS
 export BoseFSIndex, FermiFSIndex
 export BitString, SortedParticleList
 export num_particles, num_modes, num_components
-export find_occupied_mode, find_mode, occupied_modes, is_occupied, num_occupied_modes
-export excitation, near_uniform, OccupiedModeMap, OccupiedPairsMap
+export find_occupied_mode, find_mode, occupied_modes, unoccupied_modes, is_occupied, num_occupied_modes, num_unoccupied_modes
+export excitation, near_uniform, ModeMap, OccupiedPairsMap, occupied_mode_map, unoccupied_mode_map
 export onr, occupation_number_representation
 export hopnextneighbour, bose_hubbard_interaction
 export @fs_str
diff --git a/src/BitStringAddresses/bitstring.jl b/src/BitStringAddresses/bitstring.jl
index 7255eb87d..86d7634f0 100644
--- a/src/BitStringAddresses/bitstring.jl
+++ b/src/BitStringAddresses/bitstring.jl
@@ -835,3 +835,42 @@ function LinearAlgebra.dot(
 ) where {S<:BitString}
     return count_ones(occ_a.storage & occ_b.storage)
 end
+
+
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString})
+    c = 0
+    chunk = o.storage.chunks[end]
+    while iszero(chunk)
+        c += 1
+        chunk = o.storage.chunks[end-c]
+    end
+    zeros = trailing_zeros(chunk % Int)
+    return iterate(o, (chunk >> (zeros % UInt64), c * 64 + zeros, c))
+end
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString}, st)
+    chunk, index, c = st
+    while iszero(chunk)
+        c += 1
+        c == num_chunks(o.storage) && return nothing
+        chunk = o.storage.chunks[end-c]
+        index = c * 64
+    end
+    zeros = trailing_zeros(chunk % Int)
+    index += zeros
+    chunk >>= zeros
+    return FermiFSIndex(0, index + 1, index), (chunk >> 1, index + 1, c)
+end
+
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{<:Any,1,T}}) where {T}
+    chunk = o.storage.chunks[end]
+    zeros = trailing_zeros(chunk % Int)
+    return iterate(o, (chunk >> (zeros % T), zeros))
+end
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{<:Any,1,T}}, st) where {T}
+    chunk, index = st
+    iszero(chunk) && return nothing
+    chunk >>= 0x1
+    index += 1
+    zeros = trailing_zeros(chunk % Int)
+    return FermiFSIndex(0, index, index - 1), (chunk >> (zeros % T), index + zeros)
+end
\ No newline at end of file
diff --git a/src/BitStringAddresses/fermifs.jl b/src/BitStringAddresses/fermifs.jl
index 94a548292..184e950af 100644
--- a/src/BitStringAddresses/fermifs.jl
+++ b/src/BitStringAddresses/fermifs.jl
@@ -131,6 +131,50 @@ Base.:(==)(a::FermiFS, b::FermiFS) = a.bs == b.bs
 num_occupied_modes(::FermiFS{N}) where {N} = N
 occupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiOccupiedModes{N,S}(a.bs)
 
+num_unoccupied_modes(::FermiFS{N,M}) where {N,M} = M - N
+unoccupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiUnoccupiedModes{N,S}(~a.bs)
+
+"""
+    unoccupied_mode_map(addr::FermiFS) <: AbstractVector
+    
+Get a map of unoccupied modes in [`FermiFS`](@ref) address as an `AbstractVector`
+of indices compatible with [`excitation`](@ref).
+
+`unoccupied_mode_map(addr)[i]` contains the index for the `i`-th unoccupied mode.
+This is useful because unoccupied modes is required in some cases.
+`unoccupied_mode_map(addr)` is an eager version of the iterator returned by
+[`unoccupied_modes`](@ref). It is similar to [`onr`](@ref) but contains more information.
+
+# Example
+
+```jldoctest
+julia> f = FermiFS(1,1,0,0)
+FermiFS{2,4}(1, 1, 0, 0)
+
+julia> mf = unoccupied_mode_map(f)
+2-element ModeMap{2, FermiFSIndex}:
+ FermiFSIndex(occnum=0, mode=3, offset=2)
+ FermiFSIndex(occnum=0, mode=4, offset=3)
+ 
+julia> mf == collect(unoccupied_modes(f))
+true
+
+```
+See also [`SingleComponentFockAddress`](@ref).
+"""
+function unoccupied_mode_map(addr::FermiFS{N,M}) where {N,M}
+    modes = unoccupied_modes(addr)
+    T = eltype(modes)
+    L = num_unoccupied_modes(addr)
+    indices = MVector{L,T}(undef)
+    i = 0
+    for index in modes
+        i += 1
+        @inbounds indices[i] = index
+    end
+    return ModeMap(SVector(indices), i)
+end
+
 function near_uniform(::Type{FermiFS{N,M}}) where {N,M}
     return FermiFS([fill(1, N); fill(0, M - N)])
 end
diff --git a/src/BitStringAddresses/fockaddress.jl b/src/BitStringAddresses/fockaddress.jl
index c9fea6cb4..45819a9e3 100644
--- a/src/BitStringAddresses/fockaddress.jl
+++ b/src/BitStringAddresses/fockaddress.jl
@@ -50,7 +50,7 @@ Implemented subtypes: [`BoseFS`](@ref), [`FermiFS`](@ref).
 * [`find_occupied_mode`](@ref)
 * [`num_occupied_modes`](@ref)
 * [`occupied_modes`](@ref): Lazy iterator.
-* [`OccupiedModeMap`](@ref): `AbstractVector` with eager construction.
+* [`occupied_mode_map`](@ref): `AbstractVector` with eager construction.
 * [`excitation`](@ref): Create a new address.
 * [`BoseFSIndex`](@ref) and [`FermiFSIndex`](@ref) for indexing.
 
@@ -109,7 +109,7 @@ julia> find_occupied_mode(BoseFS(1, 0, 2), 1, 2)
 BoseFSIndex(occnum=2, mode=3, offset=3)
 ```
 
-See also [`occupied_modes`](@ref), [`OccupiedModeMap`](@ref),
+See also [`occupied_modes`](@ref), [`occupied_mode_map`](@ref),
 [`SingleComponentFockAddress`](@ref).
 """
 function find_occupied_mode(b::SingleComponentFockAddress, index::Integer, n=1)
@@ -149,7 +149,7 @@ num_occupied_modes
 
 Return a lazy iterator over all occupied modes in an address. Iterates over
 [`BoseFSIndex`](@ref)s for [`BoseFS`](@ref), and over [`FermiFSIndex`](@ref)s for
-[`FermiFS`](@ref). See [`OccupiedModeMap`](@ref) for an eager version.
+[`FermiFS`](@ref). See [`occupied_mode_map`](@ref) for an eager version.
 
 # Example
 
@@ -176,6 +176,27 @@ See also [`find_occupied_mode`](@ref),
 """
 occupied_modes
 
+"""
+    unoccupied_modes(::FermiFS)
+
+Return a lazy iterator over all unoccupied modes in an Fermi type address. Iterates over
+over [`FermiFSIndex`](@ref)s for [`FermiFS`](@ref). 
+See [`unoccupied_mode_map`](@ref) for an eager version.
+
+# Example
+
+```jldoctest
+julia> f = FermiFS((1,1,0,1,0,0,1));
+
+julia> foreach(println, unoccupied_modes(f))
+FermiFSIndex(occnum=0, mode=3, offset=2)
+FermiFSIndex(occnum=0, mode=5, offset=4)
+FermiFSIndex(occnum=0, mode=6, offset=5)
+```
+See also [`find_occupied_mode`](@ref), [`SingleComponentFockAddress`](@ref).
+"""
+unoccupied_modes
+
 """
     excitation(addr::SingleComponentFockAddress, creations::NTuple, destructions::NTuple)
 
@@ -209,16 +230,39 @@ See [`SingleComponentFockAddress`](@ref).
 """
 excitation
 
+
 """
-    OccupiedModeMap(addr) <: AbstractVector
+    ModeMap <: AbstractVector
+
+A unified storage structure for indices of `SingleComponentFockAddress`. 
+It stores the FSIndex of corresponding address as an `AbstractVector` compatible with
+[`excitation`](@ref) - [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref).
 
+This struct is not intended to construct directly. Use [`occupied_mode_map`](@ref) or 
+[`unoccupied_mode_map`](@ref) to obtain an instance.
+
+See also [`dot`](@ref Main.Hamiltonians.dot), [`SingleComponentFockAddress`](@ref).
+"""
+struct ModeMap{N,T} <: AbstractVector{T}
+    indices::SVector{N,T} # N = min(N, M)
+    length::Int
+end
+
+Base.eltype(::ModeMap{N,T}) where {N,T} = T
+
+Base.@deprecate OccupiedModeMap(addr) occupied_mode_map(addr)
+
+
+"""
+    occupied_mode_map(addr) <: AbstractVector
+    
 Get a map of occupied modes in address as an `AbstractVector` of indices compatible with
 [`excitation`](@ref) - [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref).
 
-`OccupiedModeMap(addr)[i]` contains the index for the `i`-th occupied mode.
+`occupied_mode_map(addr)[i]` contains the index for the `i`-th occupied mode.
 This is useful because repeatedly looking for occupied modes with
 [`find_occupied_mode`](@ref) can be time-consuming.
-`OccupiedModeMap(addr)` is an eager version of the iterator returned by
+`occupied_mode_map(addr)` is an eager version of the iterator returned by
 [`occupied_modes`](@ref). It is similar to [`onr`](@ref) but contains more information.
 
 # Example
@@ -227,8 +271,8 @@ This is useful because repeatedly looking for occupied modes with
 julia> b = BoseFS(10, 0, 0, 0, 2, 0, 1)
 BoseFS{13,7}(10, 0, 0, 0, 2, 0, 1)
 
-julia> mb = OccupiedModeMap(b)
-3-element OccupiedModeMap{7, BoseFSIndex}:
+julia> mb = occupied_mode_map(b)
+3-element ModeMap{7, BoseFSIndex}:
  BoseFSIndex(occnum=10, mode=1, offset=0)
  BoseFSIndex(occnum=2, mode=5, offset=14)
  BoseFSIndex(occnum=1, mode=7, offset=18)
@@ -236,8 +280,8 @@ julia> mb = OccupiedModeMap(b)
 julia> f = FermiFS(1,1,1,1,0,0,1,0,0)
 FermiFS{5,9}(1, 1, 1, 1, 0, 0, 1, 0, 0)
 
-julia> mf = OccupiedModeMap(f)
-5-element OccupiedModeMap{5, FermiFSIndex}:
+julia> mf = occupied_mode_map(f)
+5-element ModeMap{5, FermiFSIndex}:
  FermiFSIndex(occnum=1, mode=1, offset=0)
  FermiFSIndex(occnum=1, mode=2, offset=1)
  FermiFSIndex(occnum=1, mode=3, offset=2)
@@ -255,12 +299,7 @@ julia> dot(mf, 1:20)
 ```
 See also [`dot`](@ref Main.Hamiltonians.dot), [`SingleComponentFockAddress`](@ref).
 """
-struct OccupiedModeMap{N,T} <: AbstractVector{T}
-    indices::SVector{N,T} # N = min(N, M)
-    length::Int
-end
-
-function OccupiedModeMap(addr::SingleComponentFockAddress{N,M}) where {N,M}
+function occupied_mode_map(addr::SingleComponentFockAddress{N,M}) where {N,M}
     modes = occupied_modes(addr)
     T = eltype(modes)
     # There are at most N occupied modes. This could be also @generated for cases where N ≫ M
@@ -271,36 +310,36 @@ function OccupiedModeMap(addr::SingleComponentFockAddress{N,M}) where {N,M}
         i += 1
         @inbounds indices[i] = index
     end
-    return OccupiedModeMap(SVector(indices), i)
+    return ModeMap(SVector(indices), i)
 end
 
-Base.size(om::OccupiedModeMap) = (om.length,)
-function Base.getindex(om::OccupiedModeMap, i)
+Base.size(om::ModeMap) = (om.length,)
+function Base.getindex(om::ModeMap, i)
     @boundscheck 1 ≤ i ≤ om.length || throw(BoundsError(om, i))
     return om.indices[i]
 end
 
 """
-    abstract type OccupiedModeIterator
+    abstract type ModeIterator
 
-Iterator over occupied modes with `eltype` [`BoseFSIndex`](@ref) or
+Iterator over modes with `eltype` [`BoseFSIndex`](@ref) or
 [`FermiFSIndex`](@ref). A subtype of this should be returned when calling
 [`occupied_modes`](@ref) on a Fock state.
 """
-abstract type OccupiedModeIterator end
+abstract type ModeIterator end
 
 """
-    dot(map::OccupiedModeMap, vec::AbstractVector)
-    dot(map1::OccupiedModeMap, map2::OccupiedModeMap)
-Dot product extracting mode occupation numbers from an [`OccupiedModeMap`](@ref) similar
+    dot(map::ModeMap, vec::AbstractVector)
+    dot(map1::ModeMap, map2::ModeMap)
+Dot product extracting mode occupation numbers from an [`ModeMap`](@ref) similar
 to [`onr`](@ref).
 
 ```jldoctest
 julia> b = BoseFS(10, 0, 0, 0, 2, 0, 1)
 BoseFS{13,7}(10, 0, 0, 0, 2, 0, 1)
 
-julia> mb = OccupiedModeMap(b)
-3-element OccupiedModeMap{7, BoseFSIndex}:
+julia> mb = occupied_mode_map(b)
+3-element ModeMap{7, BoseFSIndex}:
  BoseFSIndex(occnum=10, mode=1, offset=0)
  BoseFSIndex(occnum=2, mode=5, offset=14)
  BoseFSIndex(occnum=1, mode=7, offset=18)
@@ -313,18 +352,18 @@ true
 ```
 See also [`SingleComponentFockAddress`](@ref).
 """
-function LinearAlgebra.dot(map::OccupiedModeMap, vec::AbstractVector)
+function LinearAlgebra.dot(map::ModeMap, vec::AbstractVector)
     value = zero(eltype(vec))
     for index in map
         value += vec[index.mode] * index.occnum
     end
     return value
 end
-LinearAlgebra.dot(vec::AbstractVector, map::OccupiedModeMap) = dot(map, vec)
+LinearAlgebra.dot(vec::AbstractVector, map::ModeMap) = dot(map, vec)
 
 # Defined for consistency. Could also be used to compute cross-component interactions in
 # real space.
-function LinearAlgebra.dot(map1::OccupiedModeMap, map2::OccupiedModeMap)
+function LinearAlgebra.dot(map1::ModeMap, map2::ModeMap)
     i = j = 1
     value = 0
     while i ≤ length(map1) && j ≤ length(map2)
@@ -535,7 +574,7 @@ See [`occupied_modes`](@ref).
 Defining `Base.length` and `Base.iterate` for this struct is a part of the interface for an
 underlying storage format used by [`BoseFS`](@ref).
 """
-struct BoseOccupiedModes{N,M,S}<:OccupiedModeIterator
+struct BoseOccupiedModes{N,M,S} <: ModeIterator
     storage::S
 end
 Base.eltype(::BoseOccupiedModes) = BoseFSIndex
@@ -644,13 +683,24 @@ Base.show(io::IO, ::MIME"text/plain", i::FermiFSIndex) = show(io, i)
 
 Iterator over occupied modes in address. `N` is the number of fermions. See [`occupied_modes`](@ref).
 """
-struct FermiOccupiedModes{N,S}<:OccupiedModeIterator
+struct FermiOccupiedModes{N,S} <: ModeIterator
     storage::S
 end
 
 Base.length(::FermiOccupiedModes{N}) where {N} = N
 Base.eltype(::FermiOccupiedModes) = FermiFSIndex
 
+"""
+    FermiUnoccupiedModes{N,S<:BitString}
+
+Iterator over unoccupied modes in address. `N` is the number of fermions. See [`unoccupied_modes`](@ref).
+"""
+struct FermiUnoccupiedModes{N,S} <: ModeIterator
+    storage::S
+end
+Base.length(::FermiUnoccupiedModes{N}) where {N} = N
+Base.eltype(::FermiUnoccupiedModes) = FermiFSIndex
+
 """
     from_fermi_onr(::Type{B}, onr) -> B
 
@@ -689,7 +739,7 @@ fermi_excitation
 ###
 ### General
 ###
-function LinearAlgebra.dot(occ_a::OccupiedModeIterator, occ_b::OccupiedModeIterator)
+function LinearAlgebra.dot(occ_a::ModeIterator, occ_b::ModeIterator)
     (n_a, i_a, _), st_a = iterate(occ_a)
     (n_b, i_b, _), st_b = iterate(occ_b)
 
@@ -757,7 +807,7 @@ julia> excitation(addr, pairs[2], pairs[4])
 (BoseFS{13,7}(9, 0, 0, 0, 4, 0, 0), 10.954451150103322)
 ```
 
-See also [`OccupiedModeMap`](@ref).
+See also [`occupied_mode_map`](@ref).
 """
 struct OccupiedPairsMap{N,T} <: AbstractVector{T}
     pairs::SVector{N,T}
@@ -765,7 +815,7 @@ struct OccupiedPairsMap{N,T} <: AbstractVector{T}
 end
 
 function OccupiedPairsMap(addr::SingleComponentFockAddress{N}) where {N}
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     T = eltype(omm)
     P = N * (N - 1) ÷ 2
     pairs = MVector{P,Tuple{T,T}}(undef)
diff --git a/src/BitStringAddresses/occupationnumberfs.jl b/src/BitStringAddresses/occupationnumberfs.jl
index 33a2349de..2ddd0cc72 100644
--- a/src/BitStringAddresses/occupationnumberfs.jl
+++ b/src/BitStringAddresses/occupationnumberfs.jl
@@ -251,4 +251,4 @@ function Base.iterate(bom::BoseOccupiedModes{<:Any,<:Any,<:OccupationNumberFS},
 end
 
 # find_occupied_modes provided by generic implementation
-# OccupiedModeMap provided by generic implementation
+# occupied_mode_map provided by generic implementation
diff --git a/src/Hamiltonians/ExtendedHubbardMom1D.jl b/src/Hamiltonians/ExtendedHubbardMom1D.jl
index 40ed633a0..2af72ad62 100644
--- a/src/Hamiltonians/ExtendedHubbardMom1D.jl
+++ b/src/Hamiltonians/ExtendedHubbardMom1D.jl
@@ -91,25 +91,25 @@ end
 end
 
 @inline function diagonal_element(h::ExtendedHubbardMom1D{<:Any,M,A}, address::A) where {M,A<:SingleComponentFockAddress}
-    map = OccupiedModeMap(address)
+    map = occupied_mode_map(address)
     return (dot(h.kes, map) + (h.u/ 2M) * momentum_transfer_diagonal(map) 
         + (h.v/ M) * extended_momentum_transfer_diagonal(map, 2π / M))
 end
 
 @inline function diagonal_element(h::ExtendedHubbardMom1D{<:Any,M,A}, address::A) where {M,A<:FermiFS}
-    map = OccupiedModeMap(address)
+    map = occupied_mode_map(address)
     return dot(h.kes, map) + (h.v/ M) * extended_momentum_transfer_diagonal(map, 2π / M)
 end
 
 @inline function get_offdiagonal(
-    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=OccupiedModeMap(address)
+    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=occupied_mode_map(address)
 ) where {M,A<:SingleComponentFockAddress}
     address, onproduct,_,_,q = momentum_transfer_excitation(address, chosen, map)
     return address, ham.u * onproduct / 2M + ham.v * cos(q * 2π / M) * onproduct / M
 end
 
 @inline function get_offdiagonal(
-    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=OccupiedModeMap(address)
+    ham::ExtendedHubbardMom1D{<:Any,M,A}, address::A, chosen, map=occupied_mode_map(address)
 ) where {M,A<:FermiFS}
     address, onproduct,_,_,q = momentum_transfer_excitation(address, chosen, map)
     return address, -ham.v * onproduct * cos(q * 2π / M) / M
diff --git a/src/Hamiltonians/ExtendedHubbardReal1D.jl b/src/Hamiltonians/ExtendedHubbardReal1D.jl
index 34f86c61e..357b4491d 100644
--- a/src/Hamiltonians/ExtendedHubbardReal1D.jl
+++ b/src/Hamiltonians/ExtendedHubbardReal1D.jl
@@ -117,7 +117,7 @@ where ``f_{j-i} = 1`` for nearest neighbors (power = nothing) and
 See [`ExtendedHubbardReal1D`](@ref) and [`hopnextneighbour`](@ref).
 """
 function extended_hubbard_interaction(h::ExtendedHubbardReal1D, b::SingleComponentFockAddress, ::Nothing)
-    omm = OccupiedModeMap(b)
+    omm = occupied_mode_map(b)
 
     prev = zero(eltype(omm))
     ext_result = 0
@@ -139,7 +139,7 @@ function extended_hubbard_interaction(h::ExtendedHubbardReal1D, b::SingleCompone
 end
 
 function extended_hubbard_interaction(h::ExtendedHubbardReal1D, b::SingleComponentFockAddress, power::Number)
-    omm = OccupiedModeMap(b)
+    omm = occupied_mode_map(b)
     M = num_modes(b)
     ext_result = 0
     reg_result = 0
diff --git a/src/Hamiltonians/G2MomCorrelator.jl b/src/Hamiltonians/G2MomCorrelator.jl
index 08534d389..500d26497 100644
--- a/src/Hamiltonians/G2MomCorrelator.jl
+++ b/src/Hamiltonians/G2MomCorrelator.jl
@@ -71,7 +71,7 @@ function get_offdiagonal(
     chosen,
 )::Tuple{A,ComplexF64} where {A<:SingleComponentFockAddress}
     M = num_modes(addr)
-    new_add, gamma, Δp = momentum_transfer_excitation(addr, chosen, OccupiedModeMap(addr))
+    new_add, gamma, Δp = momentum_transfer_excitation(addr, chosen, occupied_mode_map(addr))
     gd = exp(-im * g.d * Δp * 2π / M) * gamma
     return new_add, ComplexF64(gd / M)
 end
diff --git a/src/Hamiltonians/HOCartesianCentralImpurity.jl b/src/Hamiltonians/HOCartesianCentralImpurity.jl
index 0e48ff46d..d84d5c53a 100644
--- a/src/Hamiltonians/HOCartesianCentralImpurity.jl
+++ b/src/Hamiltonians/HOCartesianCentralImpurity.jl
@@ -148,10 +148,10 @@ LOStructure(::Type{<:HOCartesianCentralImpurity}) = IsHermitian()
         dot(indices, aspect) + sum(aspect)/2
     end
 end
-noninteracting_energy(H::HOCartesianCentralImpurity, addr) = noninteracting_energy(H.S, H.aspect, OccupiedModeMap(addr))
+noninteracting_energy(H::HOCartesianCentralImpurity, addr) = noninteracting_energy(H.S, H.aspect, occupied_mode_map(addr))
 
 @inline function diagonal_element(H::HOCartesianCentralImpurity, addr)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     u = H.g / sqrt(prod(H.aspect))
     result = u * sum(p -> ho_delta_potential(H.S, p.mode, p.mode; vals = H.vtable), omm)
     if !H.impurity_only
@@ -173,7 +173,7 @@ num_offdiagonals(H::HOCartesianCentralImpurity, addr) = (num_modes(addr) - 1) *
 Specialized [`AbstractOffdiagonals`](@ref) for [`HOCartesianCentralImpurity`](@ref).
 """
 struct HOCartImpurityOffdiagonals{
-    A<:BoseFS,O<:OccupiedModeMap,H<:HOCartesianCentralImpurity
+    A<:BoseFS,O<:ModeMap,H<:HOCartesianCentralImpurity
 } <: AbstractOffdiagonals{A,Float64}
     ham::H
     address::A
@@ -183,7 +183,7 @@ struct HOCartImpurityOffdiagonals{
 end
 
 function offdiagonals(H::HOCartesianCentralImpurity, addr::SingleComponentFockAddress)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     num_offs = num_offdiagonals(H, addr)
     u = H.g / sqrt(prod(H.aspect))
     return HOCartImpurityOffdiagonals(H, addr, omm, u, num_offs)
diff --git a/src/Hamiltonians/HOCartesianContactInteractions.jl b/src/Hamiltonians/HOCartesianContactInteractions.jl
index 6f1808b18..19c5fe699 100644
--- a/src/Hamiltonians/HOCartesianContactInteractions.jl
+++ b/src/Hamiltonians/HOCartesianContactInteractions.jl
@@ -272,14 +272,14 @@ end
 
 noninteracting_energy(h::HOCartesianContactInteractions, omm::BoseOccupiedModeMap) = dot(h.energies, omm)
 @inline function noninteracting_energy(h::HOCartesianContactInteractions, addr::BoseFS)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     return noninteracting_energy(h, omm)
 end
 # fast method for finding blocks
 noninteracting_energy(h::HOCartesianContactInteractions, t::Union{Vector{Int64},NTuple{N,Int64}}) where {N} = sum(h.energies[j] for j in t)
 
 @inline function diagonal_element(h::HOCartesianContactInteractions, addr::BoseFS)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     return noninteracting_energy(h, omm) + energy_transfer_diagonal(h, omm)
 end
 
diff --git a/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl b/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
index 75b561e06..163a650ef 100644
--- a/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
+++ b/src/Hamiltonians/HOCartesianEnergyConservedPerDim.jl
@@ -20,15 +20,15 @@ function largest_two_point_box(i::Int, j::Int, dims::NTuple{D,Int}) where {D}
 end
 
 """
-    find_chosen_pair_moves(omm::OccupiedModeMap, c, S) -> p_i, p_j, c, box_ranges
+    find_chosen_pair_moves(omm::ModeMap, c, S) -> p_i, p_j, c, box_ranges
 
 Find size of valid moves for chosen pair of indices in `omm`. Returns two valid indices `p_i` and `p_j`
 for the initial pair and a tuple of ranges `box_ranges` defining the subbox of valid moves.
 The index for the chosen move `c` is updated to be valid for this box.
-Arguments are the `OccupiedModeMap` `omm`, the chosen move index `c`
+Arguments are the `ModeMap` `omm`, the chosen move index `c`
 and the size of the basis grid `S`.
 """
-function find_chosen_pair_moves(omm::OccupiedModeMap, c, S::Tuple)
+function find_chosen_pair_moves(omm::ModeMap, c, S::Tuple)
     for i in eachindex(omm)
         p_i = omm[i]
         if p_i.occnum > 1
@@ -196,14 +196,14 @@ end
 
 noninteracting_energy(h::HOCartesianEnergyConservedPerDim, omm::BoseOccupiedModeMap) = dot(h.energies, omm)
 @inline function noninteracting_energy(h::HOCartesianEnergyConservedPerDim, addr::BoseFS)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     return noninteracting_energy(h, omm)
 end
 # fast method for finding blocks
 noninteracting_energy(h::HOCartesianEnergyConservedPerDim, t::Union{Vector{Int64},NTuple{N,Int64}}) where {N} = sum(h.energies[j] for j in t)
 
 @inline function diagonal_element(h::HOCartesianEnergyConservedPerDim, addr::BoseFS)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     return noninteracting_energy(h, omm) + energy_transfer_diagonal(h, omm)
 end
 
@@ -213,7 +213,7 @@ end
 # To-Do: optimise these out for FCIQMC
 function num_offdiagonals(h::HOCartesianEnergyConservedPerDim, addr::BoseFS)
     S = h.S
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     noffs = 0
 
     for i in eachindex(omm)
@@ -233,7 +233,7 @@ end
 
 
 """
-    energy_transfer_offdiagonal(S, addr, chosen, omm = OccupiedModeMap(addr))
+    energy_transfer_offdiagonal(S, addr, chosen, omm = occupied_mode_map(addr))
         -> new_add, val, mode_i, mode_j, mode_l
 
 Return the new address `new_add`, the prefactor `val`, the initial particle modes
@@ -244,7 +244,7 @@ function energy_transfer_offdiagonal(
         S::Tuple,
         addr::BoseFS,
         chosen::Int,
-        omm::BoseOccupiedModeMap = OccupiedModeMap(addr)
+    omm::BoseOccupiedModeMap=occupied_mode_map(addr)
     )
     # find size of valid moves for each pair
     particle_i, particle_j, valid_box_ranges, chosen = find_chosen_pair_moves(omm, chosen, S)
@@ -271,7 +271,7 @@ function get_offdiagonal(
         h::HOCartesianEnergyConservedPerDim{D,A},
         addr::A,
         chosen::Int,
-        omm::BoseOccupiedModeMap = OccupiedModeMap(addr)
+    omm::BoseOccupiedModeMap=occupied_mode_map(addr)
     ) where {D,A}
 
     S = h.S
@@ -306,7 +306,7 @@ end
 Specialized [`AbstractOffdiagonals`](@ref) for [`HOCartesianEnergyConservedPerDim`](@ref).
 """
 struct HOCartSeparableOffdiagonals{
-    A<:BoseFS,T,H<:AbstractHamiltonian{T},O<:OccupiedModeMap
+    A<:BoseFS,T,H<:AbstractHamiltonian{T},O<:ModeMap
 } <: AbstractOffdiagonals{A,T}
     hamiltonian::H
     address::A
@@ -315,7 +315,7 @@ struct HOCartSeparableOffdiagonals{
 end
 
 function offdiagonals(h::HOCartesianEnergyConservedPerDim, addr::BoseFS)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     num = num_offdiagonals(h, addr)
     return HOCartSeparableOffdiagonals(h, addr, num, omm)
 end
diff --git a/src/Hamiltonians/HubbardMom1D.jl b/src/Hamiltonians/HubbardMom1D.jl
index 56a66c1c9..2d9d3b954 100644
--- a/src/Hamiltonians/HubbardMom1D.jl
+++ b/src/Hamiltonians/HubbardMom1D.jl
@@ -144,7 +144,7 @@ end
 end
 
 """
-    momentum_transfer_diagonal(H, map::OccupiedModeMap)
+    momentum_transfer_diagonal(H, map::ModeMap)
 
 Compute diagonal interaction energy term.
 
@@ -156,7 +156,7 @@ BoseFS{6,5}(1, 2, 3, 0, 0)
 
 julia> H = HubbardMom1D(a);
 
-julia> Hamiltonians.momentum_transfer_diagonal(H, OccupiedModeMap(a))
+julia> Hamiltonians.momentum_transfer_diagonal(H, occupied_mode_map(a))
 5.2
 ```
 """
@@ -172,30 +172,30 @@ end
 end
 
 @inline function diagonal_element(h::HubbardMom1D, address::SingleComponentFockAddress)
-    map = OccupiedModeMap(address)
+    map = occupied_mode_map(address)
     return dot(h.kes, map) + momentum_transfer_diagonal(h, map)
 end
 @inline function diagonal_element(h::HubbardMom1D, address::FermiFS)
-    map = OccupiedModeMap(address)
+    map = occupied_mode_map(address)
     return dot(h.kes, map)
 end
 @inline function diagonal_element(h::HubbardMom1D, address::FermiFS2C)
-    map_a = OccupiedModeMap(address.components[1])
-    map_b = OccupiedModeMap(address.components[2])
+    map_a = occupied_mode_map(address.components[1])
+    map_b = occupied_mode_map(address.components[2])
     return dot(h.kes, map_a) + dot(h.kes, map_b) +
         momentum_transfer_diagonal(h, map_a, map_b)
 end
 
 @inline function get_offdiagonal(
-    ham::HubbardMom1D{<:Any,M,A}, address::A, chosen, map=OccupiedModeMap(address)
+    ham::HubbardMom1D{<:Any,M,A}, address::A, chosen, map=occupied_mode_map(address)
 ) where {M,A<:SingleComponentFockAddress}
     address, onproduct = momentum_transfer_excitation(address, chosen, map)
     return address, ham.u/(2*M)*onproduct
 end
 @inline function get_offdiagonal(
     ham::HubbardMom1D{<:Any,M,A}, address::A, chosen,
-    map_a=OccupiedModeMap(address.components[1]),
-    map_b=OccupiedModeMap(address.components[2])
+    map_a=occupied_mode_map(address.components[1]),
+    map_b=occupied_mode_map(address.components[2])
 ) where {M,A<:FermiFS2C}
     add_a, add_b = address.components
     new_add_a, new_add_b, onproduct = momentum_transfer_excitation(
@@ -214,7 +214,7 @@ Specialized [`AbstractOffdiagonals`](@ref) that keeps track of singly and doubly
 sites in current address.
 """
 struct OffdiagonalsBoseMom1D{
-    A<:SingleComponentFockAddress,T,H<:AbstractHamiltonian{T},O<:OccupiedModeMap
+    A<:SingleComponentFockAddress,T,H<:AbstractHamiltonian{T},O<:ModeMap
 } <: AbstractOffdiagonals{A,T}
     hamiltonian::H
     address::A
@@ -223,7 +223,7 @@ struct OffdiagonalsBoseMom1D{
 end
 
 function offdiagonals(h::HubbardMom1D, a::SingleComponentFockAddress)
-    map = OccupiedModeMap(a)
+    map = occupied_mode_map(a)
     singlies = length(map)
     doublies = count(i -> i.occnum ≥ 2, map)
     num = num_offdiagonals(h, a, singlies, doublies)
@@ -252,8 +252,8 @@ end
 
 function offdiagonals(h::HubbardMom1D, f::FermiFS2C)
     comp_a, comp_b = f.components
-    map_a = OccupiedModeMap(comp_a)
-    map_b = OccupiedModeMap(comp_b)
+    map_a = occupied_mode_map(comp_a)
+    map_b = occupied_mode_map(comp_b)
     num = num_offdiagonals(h, f)
     return OffdiagonalsFermiMom1D2C(h, f, num, map_a, map_b)
 end
diff --git a/src/Hamiltonians/HubbardMom1DEP.jl b/src/Hamiltonians/HubbardMom1DEP.jl
index 44cb28154..0d2ba8f3e 100644
--- a/src/Hamiltonians/HubbardMom1DEP.jl
+++ b/src/Hamiltonians/HubbardMom1DEP.jl
@@ -146,14 +146,14 @@ end
 end
 
 @inline function diagonal_element(h::HubbardMom1DEP, address::SingleComponentFockAddress)
-    map = OccupiedModeMap(address)
+    map = occupied_mode_map(address)
     return dot(h.kes, map) + momentum_transfer_diagonal(h, map) +
         momentum_external_potential_diagonal(h.ep, address, map)
 end
 @inline function diagonal_element(h::HubbardMom1DEP, address::FermiFS2C)
     c1, c2 = address.components
-    map_a = OccupiedModeMap(c1)
-    map_b = OccupiedModeMap(c2)
+    map_a = occupied_mode_map(c1)
+    map_b = occupied_mode_map(c2)
     return dot(h.kes, map_a) + dot(h.kes, map_b) +
         momentum_transfer_diagonal(h, map_a, map_b) +
         momentum_external_potential_diagonal(h.ep, c1, map_a) +
@@ -164,7 +164,7 @@ end
 ### offdiagonals
 ###
 struct OffdiagonalsBoseMom1DEP{
-    A<:SingleComponentFockAddress,T,H<:AbstractHamiltonian{T},O<:OccupiedModeMap
+    A<:SingleComponentFockAddress,T,H<:AbstractHamiltonian{T},O<:ModeMap
 } <: AbstractOffdiagonals{A,T}
     hamiltonian::H
     address::A
@@ -175,7 +175,7 @@ end
 
 function offdiagonals(h::HubbardMom1DEP, a::SingleComponentFockAddress)
     M = num_modes(a)
-    map = OccupiedModeMap(a)
+    map = occupied_mode_map(a)
     singlies = length(map)
     doublies = count(i -> i.occnum ≥ 2, map)
     num_mom = singlies * (singlies - 1) * (M - 2) + doublies * (M - 1)
@@ -219,8 +219,8 @@ function offdiagonals(h::HubbardMom1DEP, a::FermiFS2C)
     comp_a, comp_b = a.components
     N1 = num_particles(comp_a)
     N2 = num_particles(comp_b)
-    map_a = OccupiedModeMap(comp_a)
-    map_b = OccupiedModeMap(comp_b)
+    map_a = occupied_mode_map(comp_a)
+    map_b = occupied_mode_map(comp_b)
     num_mom = N1 * N2 * (M - 1)
     num_ep_a = N1 * (M - 1)
     num_ep_b = N2 * (M - 1)
diff --git a/src/Hamiltonians/Momentum.jl b/src/Hamiltonians/Momentum.jl
index 464bb8a50..d9fa2e0a7 100644
--- a/src/Hamiltonians/Momentum.jl
+++ b/src/Hamiltonians/Momentum.jl
@@ -33,7 +33,7 @@ end
 
 @inline function _momentum(address::SingleComponentFockAddress, fold)
     M = num_modes(address)
-    momentum = float(dot(-cld(M,2)+1:fld(M,2), OccupiedModeMap(address)))
+    momentum = float(dot(-cld(M, 2)+1:fld(M, 2), occupied_mode_map(address)))
     if fold
         return mod1(momentum + cld(M,2), M) - cld(M,2)
     else
diff --git a/src/Hamiltonians/Transcorrelated1D.jl b/src/Hamiltonians/Transcorrelated1D.jl
index 514b3eead..fbf046182 100644
--- a/src/Hamiltonians/Transcorrelated1D.jl
+++ b/src/Hamiltonians/Transcorrelated1D.jl
@@ -249,8 +249,8 @@ end
 
 function diagonal_element(h::Transcorrelated1D{<:Any,F}, add::F) where {F}
     c1, c2 = add.components
-    map1 = OccupiedModeMap(c1)
-    map2 = OccupiedModeMap(c2)
+    map1 = occupied_mode_map(c1)
+    map2 = occupied_mode_map(c2)
 
     value = dot(h.kes, map1) + dot(h.kes, map2) +
         momentum_transfer_diagonal(h, map1, map2)
@@ -277,8 +277,8 @@ end
 function offdiagonals(h::Transcorrelated1D{M,F}, add::F) where {M,F}
     offdiags = Tuple{F,Float64}[]
     c1, c2 = add.components
-    map1 = OccupiedModeMap(c1)
-    map2 = OccupiedModeMap(c2)
+    map1 = occupied_mode_map(c1)
+    map2 = occupied_mode_map(c2)
     N1 = length(map1)
     N2 = length(map2)
     n_mom = N1 * N2 * (M - 1)
diff --git a/src/Hamiltonians/angular_momentum.jl b/src/Hamiltonians/angular_momentum.jl
index 900aae7ee..e7488377d 100644
--- a/src/Hamiltonians/angular_momentum.jl
+++ b/src/Hamiltonians/angular_momentum.jl
@@ -88,7 +88,7 @@ num_offdiagonals(::AxialAngularMomentumHO, addr::SingleComponentFockAddress) = 2
 function get_offdiagonal(L::AxialAngularMomentumHO{<:Any,D}, addr::SingleComponentFockAddress, chosen::Int) where {D}
     S = L.S
     states = CartesianIndices(S)
-    omm = OccupiedModeMap(addr)
+    omm = occupied_mode_map(addr)
     x, y, z = L.xyz
 
     # mode selects current mode, b = 0,1 selects left or right branch
diff --git a/src/Hamiltonians/excitations.jl b/src/Hamiltonians/excitations.jl
index 14c27f2e5..849ef35e9 100644
--- a/src/Hamiltonians/excitations.jl
+++ b/src/Hamiltonians/excitations.jl
@@ -1,8 +1,8 @@
 using Rimu
 using Rimu.Hamiltonians: hopnextneighbour
 
-const BoseOccupiedModeMap{N} = OccupiedModeMap{N,BitStringAddresses.BoseFSIndex}
-const FermiOccupiedModeMap{N} = OccupiedModeMap{N,BitStringAddresses.FermiFSIndex}
+const BoseOccupiedModeMap{N} = ModeMap{N,BitStringAddresses.BoseFSIndex}
+const FermiOccupiedModeMap{N} = ModeMap{N,BitStringAddresses.FermiFSIndex}
 
 """
     momentum_transfer_excitation(add, chosen, map; fold=true) -> nadd, α, p, q, -k
@@ -17,14 +17,14 @@ a^†_{p + k} a^†_{q - k} a_q a_p |\\mathtt{add}⟩
 
 The `fold` argument controls whether the terms `p + k` and `q - k` are done modulo M. If
 not, zero is returned when either of those terms is less than 1 or larger than M.
-It is expected that `map == OccupiedModeMap(add)`.
+It is expected that `map == occupied_mode_map(add)`.
 
 Return the new address(es), the value, modes `p` and `q`, and the momentum change `-k`.
 
-See [`excitation`](@ref), [`OccupiedModeMap`](@ref).
+See [`excitation`](@ref), [`occupied_mode_map`](@ref).
 """
 @inline function momentum_transfer_excitation(
-    add::SingleComponentFockAddress, chosen, map::OccupiedModeMap; fold=true
+    add::SingleComponentFockAddress, chosen, map::ModeMap; fold=true
 )
     M = num_modes(add)
     singlies = length(map) # number of at least singly occupied modes
@@ -142,7 +142,7 @@ end
 The diagonal part of nearest neighbour term [`momentum_transfer_excitation`](@ref) in [`ExtendedHubbardMom1D`](@ref).
 Where `step` is the separation of single-particle momenta in the momentum grid.
 """
-function extended_momentum_transfer_diagonal(map::OccupiedModeMap, step)
+function extended_momentum_transfer_diagonal(map::ModeMap, step)
     onproduct = 0
     for i in 1:length(map)
         occ_i = map[i].occnum
@@ -189,7 +189,7 @@ a^†_{p+k,↑} a^†_{q+l,↑} a^†_{s-k-l,↓} a_{s,↓} a_{q,↑} a_{p,↑}
 ```
 
 The index `i` enumerates the possible non-zero terms and determines ``p, q, s, k, l``.
-It is expected that `map↑, map↓ == OccupiedModeMap(add↑), OccupiedModeMap(add↓)`.
+It is expected that `map↑, map↓ == occupied_mode_map(add↑), occupied_mode_map(add↓)`.
 
 Return new addresses, prefactor `value`, `k`, and `l`. Note: If either `k` or `l` are zero,
 or the excitation is diagonal, the function returns `value == 0`.
@@ -238,7 +238,7 @@ function transcorrelated_three_body_excitation(add_a, add_b, i, map_a, map_b)
 end
 
 """
-    momentum_external_potential_excitation(ep, add, i, map::OccupiedModeMap) -> nadd, α
+    momentum_external_potential_excitation(ep, add, i, map::ModeMap) -> nadd, α
 
 The momentum space version of an external potential. `ep` may be a discrete Fourier
 transform of a real-space potential.
@@ -249,12 +249,12 @@ transform of a real-space potential.
 
 Return the new address `nadd` and value `α` of the `i`th of the
 `(num_modes(add)-1) * num_occupied_modes(add)` terms in the sum, excluding the diagonal
-term `∝ |add⟩`. It is expected that `map == OccupiedModeMap(add)`.
+term `∝ |add⟩`. It is expected that `map == occupied_mode_map(add)`.
 
-See [`momentum_external_potential_diagonal`](@ref), [`OccupiedModeMap`](@ref),
+See [`momentum_external_potential_diagonal`](@ref), [`occupied_mode_map`](@ref),
 [`num_occupied_modes`](@ref), [`num_modes`](@ref).
 """
-function momentum_external_potential_excitation(ep, add, i, map::OccupiedModeMap)
+function momentum_external_potential_excitation(ep, add, i, map::ModeMap)
     M = num_modes(add)
     p, q = fldmod1(i, M - 1) # i == (p-1)*(M-1) + q
     p_index = map[p] # p-th occupied mode in add
@@ -267,11 +267,11 @@ function momentum_external_potential_excitation(ep, add, i, map::OccupiedModeMap
 end
 
 """
-    momentum_external_potential_diagonal(ep, add, map::OccupiedModeMap)
+    momentum_external_potential_diagonal(ep, add, map::ModeMap)
 
 The diagonal part of [`momentum_external_potential_excitation`](@ref).
 """
-function momentum_external_potential_diagonal(ep, add, map::OccupiedModeMap)
+function momentum_external_potential_diagonal(ep, add, map::ModeMap)
     onproduct = sum(map) do index
         index.occnum
     end
diff --git a/src/Hamiltonians/ho-cart-tools.jl b/src/Hamiltonians/ho-cart-tools.jl
index ecb58529f..a59a267cd 100644
--- a/src/Hamiltonians/ho-cart-tools.jl
+++ b/src/Hamiltonians/ho-cart-tools.jl
@@ -201,7 +201,7 @@ function get_all_blocks_parity(h::HOCartesianContactInteractions{D,<:Any,B}) whe
     # rough upper bound assuming `BasisSetRepresentation` will filter by energy cutoff
     df = DataFrame()
     for (block_id, addr) in enumerate(parity_block_seed_addresses(h))
-        t = vcat(map(p -> [p.mode for _ in 1:p.occnum], OccupiedModeMap(addr))...)
+        t = vcat(map(p -> [p.mode for _ in 1:p.occnum], occupied_mode_map(addr))...)
         block_E0 = noninteracting_energy(h, t)
         push!(df, (; block_id, block_E0, block_size = bs_estimate, addr, indices = Tuple(t), t_basis = 0.0))
     end
@@ -227,7 +227,7 @@ function fock_to_cart(
 
     cart = vcat(map(
         p -> [Tuple(states[p.mode]) .- Int(zero_index) for _ in 1:p.occnum],
-        OccupiedModeMap(addr))...)
+        occupied_mode_map(addr))...)
 
     return SVector{N,NTuple{D,Int}}(cart)
 end
diff --git a/test/BitStringAddresses.jl b/test/BitStringAddresses.jl
index 666977baa..072c386c6 100644
--- a/test/BitStringAddresses.jl
+++ b/test/BitStringAddresses.jl
@@ -549,7 +549,7 @@ end
         lfs = OccupationNumberFS{6}([1 0 0; 1 1 0])
         @test onr(lfs, LadderBoundaries(2, 3)) == [1 0 0; 1 1 0]
         @test num_occupied_modes(lfs) == length(occupied_modes(lfs)) == 3
-        @test OccupiedModeMap(lfs) == collect(occupied_modes(lfs))
+        @test occupied_mode_map(lfs) == collect(occupied_modes(lfs))
         b1, b2 = BoseFS(1, 6), BoseFS(3, 4)
         o1, o2 = OccupationNumberFS(b1), OccupationNumberFS(b2)
         @test (o1 < o2) == (b1 < b2)
diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl
index 9e9733a20..41b4a4ebe 100644
--- a/test/Hamiltonians.jl
+++ b/test/Hamiltonians.jl
@@ -105,32 +105,32 @@ using Rimu.Hamiltonians: momentum_transfer_excitation
         add1 = BoseFS((0,1,1,0))
         add2 = BoseFS((1,0,0,1))
         for i in 1:4
-            ex = momentum_transfer_excitation(add1, i, OccupiedModeMap(add1); fold=true)
+            ex = momentum_transfer_excitation(add1, i, occupied_mode_map(add1); fold=true)
             @test ex[1] == add2
             @test ex[2] == 1
 
-            ex = momentum_transfer_excitation(add1, i, OccupiedModeMap(add1); fold=false)
+            ex = momentum_transfer_excitation(add1, i, occupied_mode_map(add1); fold=false)
             @test ex[1] == add2
             @test ex[2] == 1
         end
 
         add3 = BoseFS((1,1,0,0))
         for i in 1:4
-            ex = momentum_transfer_excitation(add3, i, OccupiedModeMap(add3); fold=true)
+            ex = momentum_transfer_excitation(add3, i, occupied_mode_map(add3); fold=true)
             @test ex[2] == 1
 
-            ex = momentum_transfer_excitation(add3, i, OccupiedModeMap(add3); fold=false)
+            ex = momentum_transfer_excitation(add3, i, occupied_mode_map(add3); fold=false)
             @test ex[2] == 0
         end
 
         add4 = BoseFS((0,3,0))
         add5 = BoseFS((1,1,1))
         for i in 1:2
-            ex = momentum_transfer_excitation(add4, i, OccupiedModeMap(add4); fold=false)
+            ex = momentum_transfer_excitation(add4, i, occupied_mode_map(add4); fold=false)
             @test ex[1] == add5
             @test ex[2] ≈ √6
 
-            ex = momentum_transfer_excitation(add4, i, OccupiedModeMap(add4); fold=true)
+            ex = momentum_transfer_excitation(add4, i, occupied_mode_map(add4); fold=true)
             @test ex[1] == add5
             @test ex[2] ≈ √6
         end
@@ -138,8 +138,8 @@ using Rimu.Hamiltonians: momentum_transfer_excitation
     @testset "FermiFS" begin
         add1 = FermiFS((0,0,1,0))
         add2 = FermiFS((0,1,0,0))
-        occ1 = OccupiedModeMap(add1)
-        occ2 = OccupiedModeMap(add2)
+        occ1 = occupied_mode_map(add1)
+        occ2 = occupied_mode_map(add2)
         for i in 1:3
             ex = momentum_transfer_excitation(add1, add2, i, occ1, occ2; fold=true)
             @test ex[3] == 1
@@ -150,8 +150,8 @@ using Rimu.Hamiltonians: momentum_transfer_excitation
 
         add3 = FermiFS((1,0,0,0))
         add4 = FermiFS((0,1,0,0))
-        occ3 = OccupiedModeMap(add3)
-        occ4 = OccupiedModeMap(add4)
+        occ3 = occupied_mode_map(add3)
+        occ4 = occupied_mode_map(add4)
         for i in 1:3
             ex = momentum_transfer_excitation(add3, add4, i, occ3, occ4; fold=true)
             @test ex[3] == 1

From aed71ab9a45398d5ac297d0cd878bc936b9bf937 Mon Sep 17 00:00:00 2001
From: lch <9247311+lch@users.noreply.github.com>
Date: Fri, 4 Jul 2025 21:54:01 +0800
Subject: [PATCH 2/6] Unified ModeMap Type (#330)

1. Apply suggestions from code review.
2. Import ModeMap in Hamiltonians module explicitly.

Co-authored-by: Joachim Brand <joachim.brand@gmail.com>
---
 src/BitStringAddresses/BitStringAddresses.jl |  2 +-
 src/BitStringAddresses/fermifs.jl            |  7 ++++---
 src/BitStringAddresses/fockaddress.jl        | 13 +++++++------
 src/Hamiltonians/Hamiltonians.jl             |  1 +
 4 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/BitStringAddresses/BitStringAddresses.jl b/src/BitStringAddresses/BitStringAddresses.jl
index 3ec8dbb8c..ed562317f 100644
--- a/src/BitStringAddresses/BitStringAddresses.jl
+++ b/src/BitStringAddresses/BitStringAddresses.jl
@@ -19,7 +19,7 @@ export BoseFSIndex, FermiFSIndex
 export BitString, SortedParticleList
 export num_particles, num_modes, num_components
 export find_occupied_mode, find_mode, occupied_modes, unoccupied_modes, is_occupied, num_occupied_modes, num_unoccupied_modes
-export excitation, near_uniform, ModeMap, OccupiedPairsMap, occupied_mode_map, unoccupied_mode_map
+export excitation, near_uniform, OccupiedPairsMap, occupied_mode_map, unoccupied_mode_map
 export onr, occupation_number_representation
 export hopnextneighbour, bose_hubbard_interaction
 export @fs_str
diff --git a/src/BitStringAddresses/fermifs.jl b/src/BitStringAddresses/fermifs.jl
index 184e950af..1d811d511 100644
--- a/src/BitStringAddresses/fermifs.jl
+++ b/src/BitStringAddresses/fermifs.jl
@@ -143,7 +143,8 @@ of indices compatible with [`excitation`](@ref).
 `unoccupied_mode_map(addr)[i]` contains the index for the `i`-th unoccupied mode.
 This is useful because unoccupied modes is required in some cases.
 `unoccupied_mode_map(addr)` is an eager version of the iterator returned by
-[`unoccupied_modes`](@ref). It is similar to [`onr`](@ref) but contains more information.
+[`unoccupied_modes`](@ref). It is similar to 
+[`onr`](@ref) but contains more information.
 
 # Example
 
@@ -152,7 +153,7 @@ julia> f = FermiFS(1,1,0,0)
 FermiFS{2,4}(1, 1, 0, 0)
 
 julia> mf = unoccupied_mode_map(f)
-2-element ModeMap{2, FermiFSIndex}:
+2-element Rimu.BitStringAddresses.ModeMap{2, FermiFSIndex}:
  FermiFSIndex(occnum=0, mode=3, offset=2)
  FermiFSIndex(occnum=0, mode=4, offset=3)
  
@@ -160,7 +161,7 @@ julia> mf == collect(unoccupied_modes(f))
 true
 
 ```
-See also [`SingleComponentFockAddress`](@ref).
+See also [`occupied_mode_map`](@ref).
 """
 function unoccupied_mode_map(addr::FermiFS{N,M}) where {N,M}
     modes = unoccupied_modes(addr)
diff --git a/src/BitStringAddresses/fockaddress.jl b/src/BitStringAddresses/fockaddress.jl
index 45819a9e3..5f3e02ce0 100644
--- a/src/BitStringAddresses/fockaddress.jl
+++ b/src/BitStringAddresses/fockaddress.jl
@@ -193,7 +193,7 @@ FermiFSIndex(occnum=0, mode=3, offset=2)
 FermiFSIndex(occnum=0, mode=5, offset=4)
 FermiFSIndex(occnum=0, mode=6, offset=5)
 ```
-See also [`find_occupied_mode`](@ref), [`SingleComponentFockAddress`](@ref).
+See also [`find_occupied_mode`](@ref), [`occupied_modes`](@ref).
 """
 unoccupied_modes
 
@@ -241,7 +241,7 @@ It stores the FSIndex of corresponding address as an `AbstractVector` compatible
 This struct is not intended to construct directly. Use [`occupied_mode_map`](@ref) or 
 [`unoccupied_mode_map`](@ref) to obtain an instance.
 
-See also [`dot`](@ref Main.Hamiltonians.dot), [`SingleComponentFockAddress`](@ref).
+See also [`SingleComponentFockAddress`](@ref).
 """
 struct ModeMap{N,T} <: AbstractVector{T}
     indices::SVector{N,T} # N = min(N, M)
@@ -272,7 +272,7 @@ julia> b = BoseFS(10, 0, 0, 0, 2, 0, 1)
 BoseFS{13,7}(10, 0, 0, 0, 2, 0, 1)
 
 julia> mb = occupied_mode_map(b)
-3-element ModeMap{7, BoseFSIndex}:
+3-element Rimu.BitStringAddresses.ModeMap{7, BoseFSIndex}:
  BoseFSIndex(occnum=10, mode=1, offset=0)
  BoseFSIndex(occnum=2, mode=5, offset=14)
  BoseFSIndex(occnum=1, mode=7, offset=18)
@@ -281,7 +281,7 @@ julia> f = FermiFS(1,1,1,1,0,0,1,0,0)
 FermiFS{5,9}(1, 1, 1, 1, 0, 0, 1, 0, 0)
 
 julia> mf = occupied_mode_map(f)
-5-element ModeMap{5, FermiFSIndex}:
+5-element Rimu.BitStringAddresses.ModeMap{5, FermiFSIndex}:
  FermiFSIndex(occnum=1, mode=1, offset=0)
  FermiFSIndex(occnum=1, mode=2, offset=1)
  FermiFSIndex(occnum=1, mode=3, offset=2)
@@ -297,7 +297,8 @@ julia> dot(mf, mb)
 julia> dot(mf, 1:20)
 17
 ```
-See also [`dot`](@ref Main.Hamiltonians.dot), [`SingleComponentFockAddress`](@ref).
+See also [`dot`](@ref Main.Hamiltonians.dot), [`unoccupied_mode_map`](@ref),
+[`SingleComponentFockAddress`](@ref).
 """
 function occupied_mode_map(addr::SingleComponentFockAddress{N,M}) where {N,M}
     modes = occupied_modes(addr)
@@ -339,7 +340,7 @@ julia> b = BoseFS(10, 0, 0, 0, 2, 0, 1)
 BoseFS{13,7}(10, 0, 0, 0, 2, 0, 1)
 
 julia> mb = occupied_mode_map(b)
-3-element ModeMap{7, BoseFSIndex}:
+3-element Rimu.BitStringAddresses.ModeMap{7, BoseFSIndex}:
  BoseFSIndex(occnum=10, mode=1, offset=0)
  BoseFSIndex(occnum=2, mode=5, offset=14)
  BoseFSIndex(occnum=1, mode=7, offset=18)
diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl
index 0cc398843..599a93718 100644
--- a/src/Hamiltonians/Hamiltonians.jl
+++ b/src/Hamiltonians/Hamiltonians.jl
@@ -62,6 +62,7 @@ using StaticArrays: StaticArrays, SA, SMatrix, SVector, SArray, setindex
 using TupleTools: TupleTools
 
 using ..BitStringAddresses
+import ..BitStringAddresses: ModeMap
 using ..Interfaces
 using ..Interfaces: sum_mutating!
 import ..Interfaces: diagonal_element, num_offdiagonals, get_offdiagonal, starting_address,

From d91dbdd1a032ac3e98f866431216ec341a0c8c2b Mon Sep 17 00:00:00 2001
From: lch <9247311+lch@users.noreply.github.com>
Date: Mon, 7 Jul 2025 14:53:37 +0800
Subject: [PATCH 3/6] Unified ModeMap Type (#330)

1. Add tests to new methods according to suggestion

Co-authored-by: mtsch <matijacufar@gmail.com>
---
 src/BitStringAddresses/fockaddress.jl |  2 +-
 test/BitStringAddresses.jl            | 16 +++++++++++++++-
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/src/BitStringAddresses/fockaddress.jl b/src/BitStringAddresses/fockaddress.jl
index 5f3e02ce0..75c2df707 100644
--- a/src/BitStringAddresses/fockaddress.jl
+++ b/src/BitStringAddresses/fockaddress.jl
@@ -238,7 +238,7 @@ A unified storage structure for indices of `SingleComponentFockAddress`.
 It stores the FSIndex of corresponding address as an `AbstractVector` compatible with
 [`excitation`](@ref) - [`BoseFSIndex`](@ref) or [`FermiFSIndex`](@ref).
 
-This struct is not intended to construct directly. Use [`occupied_mode_map`](@ref) or 
+This struct is not intended to be constructed directly. Use [`occupied_mode_map`](@ref) or 
 [`unoccupied_mode_map`](@ref) to obtain an instance.
 
 See also [`SingleComponentFockAddress`](@ref).
diff --git a/test/BitStringAddresses.jl b/test/BitStringAddresses.jl
index 072c386c6..91d5cc5f9 100644
--- a/test/BitStringAddresses.jl
+++ b/test/BitStringAddresses.jl
@@ -2,7 +2,7 @@ using Rimu
 using Rimu.BitStringAddresses
 using Rimu.BitStringAddresses: num_chunks, chunks
 using Rimu.BitStringAddresses: remove_ghost_bits, has_ghost_bits
-using Rimu.BitStringAddresses: occupied_modes, update_component
+using Rimu.BitStringAddresses: occupied_modes, occupied_mode_map, unoccupied_modes, unoccupied_mode_map, update_component
 using Rimu.BitStringAddresses: parse_address
 using Rimu.BitStringAddresses: destroy, create
 using Random
@@ -229,6 +229,11 @@ end
 
         @test collect(occupied_modes(middle_empty)) == [[2, 1, 0], [8, 127, 128]]
     end
+    @testset "occupied_mode_map" begin
+        for b in (middle_full, middle_empty, two_full)
+            @test occupied_mode_map(b) == collect(occupied_modes(b))
+        end
+    end
     @testset "Randomized tests" begin
         # Note: the random number for these tests will be the same everytime. This is still
         # an ok way to look for errors.
@@ -354,9 +359,18 @@ end
         for o in (small, big, giant)
             f = FermiFS(o)
             sites = collect(occupied_modes(f))
+            @test occupied_mode_map(f) == sites
             @test getproperty.(sites, :mode) == findall(≠(0), onr(f))
         end
     end
+    @testset "unoccupied_modes" begin
+        for o in (small, big, giant)
+            f = FermiFS(o)
+            sites = collect(unoccupied_modes(f))
+            @test unoccupied_mode_map(f) == sites
+            @test getproperty.(sites, :mode) == findall(==(0), onr(f))
+        end
+    end
     @testset "Randomized Tests" begin
         function rand_onr_fermi(N, M)
             return SVector{M}(shuffle([ones(Int,N); zeros(Int,M - N)]))

From a3604f84cc2a6f6a4c698da7605484bc2490a956 Mon Sep 17 00:00:00 2001
From: lch <lchopn@gmail.com>
Date: Wed, 9 Jul 2025 20:56:33 +0800
Subject: [PATCH 4/6] Unified ModeMap Type (#330)

1. Fix bug in unoccupied_modes method
---
 src/BitStringAddresses/fermifs.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/BitStringAddresses/fermifs.jl b/src/BitStringAddresses/fermifs.jl
index 1d811d511..ffbae95bc 100644
--- a/src/BitStringAddresses/fermifs.jl
+++ b/src/BitStringAddresses/fermifs.jl
@@ -132,7 +132,7 @@ num_occupied_modes(::FermiFS{N}) where {N} = N
 occupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiOccupiedModes{N,S}(a.bs)
 
 num_unoccupied_modes(::FermiFS{N,M}) where {N,M} = M - N
-unoccupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiUnoccupiedModes{N,S}(~a.bs)
+unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S} = FermiUnoccupiedModes{M - N,S}(~a.bs)
 
 """
     unoccupied_mode_map(addr::FermiFS) <: AbstractVector

From d22beb26fdcdb242ccd6b55bcb8cfd556f428e28 Mon Sep 17 00:00:00 2001
From: lch <lchopn@gmail.com>
Date: Thu, 10 Jul 2025 03:07:20 +0800
Subject: [PATCH 5/6] Unified ModeMap Type (#330)

1. Resolve unoccupied_modes method test failed when FermiFS is
   constructed with SortedParticleList
---
 src/BitStringAddresses/fermifs.jl            |  3 ++-
 src/BitStringAddresses/fockaddress.jl        |  2 +-
 src/BitStringAddresses/sortedparticlelist.jl | 16 ++++++++++++++++
 3 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/src/BitStringAddresses/fermifs.jl b/src/BitStringAddresses/fermifs.jl
index ffbae95bc..e8b11d00b 100644
--- a/src/BitStringAddresses/fermifs.jl
+++ b/src/BitStringAddresses/fermifs.jl
@@ -132,7 +132,8 @@ num_occupied_modes(::FermiFS{N}) where {N} = N
 occupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiOccupiedModes{N,S}(a.bs)
 
 num_unoccupied_modes(::FermiFS{N,M}) where {N,M} = M - N
-unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S} = FermiUnoccupiedModes{M - N,S}(~a.bs)
+unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S<:BitString} = FermiUnoccupiedModes{M - N,typeof(a.bs)}(~a.bs)
+unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S<:SortedParticleList} = FermiUnoccupiedModes{M - N,typeof(a.bs)}(a.bs)
 
 """
     unoccupied_mode_map(addr::FermiFS) <: AbstractVector
diff --git a/src/BitStringAddresses/fockaddress.jl b/src/BitStringAddresses/fockaddress.jl
index 75c2df707..8d06722a3 100644
--- a/src/BitStringAddresses/fockaddress.jl
+++ b/src/BitStringAddresses/fockaddress.jl
@@ -694,7 +694,7 @@ Base.eltype(::FermiOccupiedModes) = FermiFSIndex
 """
     FermiUnoccupiedModes{N,S<:BitString}
 
-Iterator over unoccupied modes in address. `N` is the number of fermions. See [`unoccupied_modes`](@ref).
+Iterator over unoccupied modes in address. `N` is the number of unoccupied orbitals. See [`unoccupied_modes`](@ref).
 """
 struct FermiUnoccupiedModes{N,S} <: ModeIterator
     storage::S
diff --git a/src/BitStringAddresses/sortedparticlelist.jl b/src/BitStringAddresses/sortedparticlelist.jl
index d5ac18551..33886c6ff 100644
--- a/src/BitStringAddresses/sortedparticlelist.jl
+++ b/src/BitStringAddresses/sortedparticlelist.jl
@@ -277,3 +277,19 @@ function Base.iterate(fom::FermiOccupiedModes{<:Any,<:SortedParticleList}, i=1)
         return FermiFSIndex(res...), i
     end
 end
+
+Base.length(::FermiUnoccupiedModes{N,<:SortedParticleList}) where {N} = N
+function Base.iterate(fum::FermiUnoccupiedModes{<:Any,<:SortedParticleList{N,M}}, state=(1, 1)) where {N,M}
+    if isnothing(state)
+        return nothing
+    else
+        i, j = state
+        while i ∈ fum.storage.storage
+            i += 1
+        end
+        if i > M
+            return nothing
+        end
+        return FermiFSIndex(0, i, i - 1), (i + 1, j + 1)
+    end
+end
\ No newline at end of file

From 82045157affc510a15513cdfc18dd67b3cb0ebe9 Mon Sep 17 00:00:00 2001
From: lch <lchopn@gmail.com>
Date: Thu, 10 Jul 2025 23:04:41 +0800
Subject: [PATCH 6/6] Unified ModeMap Type (#330)

1. Rewrite the FermiUnoccupiedModes iterator based on the suggestions.
2. Improve the docstring.
---
 src/BitStringAddresses/bitstring.jl          | 64 +++++++++-----------
 src/BitStringAddresses/fermifs.jl            |  8 +--
 src/BitStringAddresses/fockaddress.jl        |  2 +-
 src/BitStringAddresses/sortedparticlelist.jl | 19 +++---
 4 files changed, 45 insertions(+), 48 deletions(-)

diff --git a/src/BitStringAddresses/bitstring.jl b/src/BitStringAddresses/bitstring.jl
index 86d7634f0..da772546b 100644
--- a/src/BitStringAddresses/bitstring.jl
+++ b/src/BitStringAddresses/bitstring.jl
@@ -836,41 +836,37 @@ function LinearAlgebra.dot(
     return count_ones(occ_a.storage & occ_b.storage)
 end
 
-
-function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString})
-    c = 0
-    chunk = o.storage.chunks[end]
-    while iszero(chunk)
-        c += 1
-        chunk = o.storage.chunks[end-c]
-    end
-    zeros = trailing_zeros(chunk % Int)
-    return iterate(o, (chunk >> (zeros % UInt64), c * 64 + zeros, c))
-end
-function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString}, st)
-    chunk, index, c = st
-    while iszero(chunk)
-        c += 1
-        c == num_chunks(o.storage) && return nothing
-        chunk = o.storage.chunks[end-c]
-        index = c * 64
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{N,NC,T}}, st=(1, 1)) where {N,NC,T}
+    l_index, c_index = st
+    chunk_size = sizeof(T) * 8
+    while c_index <= NC
+        chunk = o.storage.chunks[end+1-c_index]
+        while (chunk >> (l_index - 1) & 1 ≠ 0) && l_index <= chunk_size
+            l_index += 1
+        end
+        if l_index <= chunk_size
+            g_index = l_index + (c_index - 1) * chunk_size
+            if g_index <= N
+                return FermiFSIndex(0, g_index, g_index - 1), (l_index + 1, c_index)
+            else
+                return nothing
+            end
+        else
+            l_index = 1
+            c_index += 1
+        end
     end
-    zeros = trailing_zeros(chunk % Int)
-    index += zeros
-    chunk >>= zeros
-    return FermiFSIndex(0, index + 1, index), (chunk >> 1, index + 1, c)
+    return nothing
 end
 
-function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{<:Any,1,T}}) where {T}
-    chunk = o.storage.chunks[end]
-    zeros = trailing_zeros(chunk % Int)
-    return iterate(o, (chunk >> (zeros % T), zeros))
-end
-function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{<:Any,1,T}}, st) where {T}
-    chunk, index = st
-    iszero(chunk) && return nothing
-    chunk >>= 0x1
-    index += 1
-    zeros = trailing_zeros(chunk % Int)
-    return FermiFSIndex(0, index, index - 1), (chunk >> (zeros % T), index + zeros)
+function Base.iterate(o::FermiUnoccupiedModes{<:Any,<:BitString{N,1,T}}, st=1) where {N,T}
+    l_index = st
+    chunk_size = sizeof(T) * 8
+    while (o.storage.chunks[end] >> (l_index - 1) & 1 ≠ 0) && l_index <= chunk_size
+        l_index += 1
+    end
+    if l_index <= chunk_size && l_index <= N
+        return FermiFSIndex(0, l_index, l_index - 1), l_index + 1
+    end
+    return nothing
 end
\ No newline at end of file
diff --git a/src/BitStringAddresses/fermifs.jl b/src/BitStringAddresses/fermifs.jl
index e8b11d00b..6953c6d7e 100644
--- a/src/BitStringAddresses/fermifs.jl
+++ b/src/BitStringAddresses/fermifs.jl
@@ -132,8 +132,7 @@ num_occupied_modes(::FermiFS{N}) where {N} = N
 occupied_modes(a::FermiFS{N,<:Any,S}) where {N,S} = FermiOccupiedModes{N,S}(a.bs)
 
 num_unoccupied_modes(::FermiFS{N,M}) where {N,M} = M - N
-unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S<:BitString} = FermiUnoccupiedModes{M - N,typeof(a.bs)}(~a.bs)
-unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S<:SortedParticleList} = FermiUnoccupiedModes{M - N,typeof(a.bs)}(a.bs)
+unoccupied_modes(a::FermiFS{N,M,S}) where {N,M,S} = FermiUnoccupiedModes{M - N,S}(a.bs)
 
 """
     unoccupied_mode_map(addr::FermiFS) <: AbstractVector
@@ -144,8 +143,9 @@ of indices compatible with [`excitation`](@ref).
 `unoccupied_mode_map(addr)[i]` contains the index for the `i`-th unoccupied mode.
 This is useful because unoccupied modes is required in some cases.
 `unoccupied_mode_map(addr)` is an eager version of the iterator returned by
-[`unoccupied_modes`](@ref). It is similar to 
-[`onr`](@ref) but contains more information.
+[`unoccupied_modes`](@ref). It is similar to [`onr`](@ref) but contains more information.
+
+Note that this function is only implemented for addresses of type [`FermiFS`](@ref).
 
 # Example
 
diff --git a/src/BitStringAddresses/fockaddress.jl b/src/BitStringAddresses/fockaddress.jl
index 8d06722a3..b65ebe83a 100644
--- a/src/BitStringAddresses/fockaddress.jl
+++ b/src/BitStringAddresses/fockaddress.jl
@@ -692,7 +692,7 @@ Base.length(::FermiOccupiedModes{N}) where {N} = N
 Base.eltype(::FermiOccupiedModes) = FermiFSIndex
 
 """
-    FermiUnoccupiedModes{N,S<:BitString}
+    FermiUnoccupiedModes{N}
 
 Iterator over unoccupied modes in address. `N` is the number of unoccupied orbitals. See [`unoccupied_modes`](@ref).
 """
diff --git a/src/BitStringAddresses/sortedparticlelist.jl b/src/BitStringAddresses/sortedparticlelist.jl
index 33886c6ff..4d4d6b3d4 100644
--- a/src/BitStringAddresses/sortedparticlelist.jl
+++ b/src/BitStringAddresses/sortedparticlelist.jl
@@ -278,18 +278,19 @@ function Base.iterate(fom::FermiOccupiedModes{<:Any,<:SortedParticleList}, i=1)
     end
 end
 
-Base.length(::FermiUnoccupiedModes{N,<:SortedParticleList}) where {N} = N
-function Base.iterate(fum::FermiUnoccupiedModes{<:Any,<:SortedParticleList{N,M}}, state=(1, 1)) where {N,M}
+function Base.iterate(fum::FermiUnoccupiedModes{<:Any,<:SortedParticleList{N,M}}, state=(1, 1, 1)) where {N,M}
     if isnothing(state)
         return nothing
     else
-        i, j = state
-        while i ∈ fum.storage.storage
-            i += 1
-        end
-        if i > M
-            return nothing
+        mode, idx, count = state
+        while mode <= M
+            if idx <= N && mode == fum.storage.storage[idx]
+                mode += 1
+                idx += 1
+            else
+                return FermiFSIndex(0, mode, mode - 1), (mode + 1, idx, count + 1)
+            end
         end
-        return FermiFSIndex(0, i, i - 1), (i + 1, j + 1)
+        return nothing
     end
 end
\ No newline at end of file