ParticleNumberOperator (#255)

* ParticleNumberOperator

* use recent julia version for docs generation

---------

Co-authored-by: Joachim Brand <[email protected]>

.github/workflows/docs.yml

@@ -15,7 +15,7 @@ jobs:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
-          version: '1.7'
+          version: '1'
       - uses: julia-actions/cache@v1
       - name: Install dependencies
         run: |

docs/src/hamiltonians.md

@@ -102,6 +102,7 @@ Observables are [`AbstractHamiltonian`](@ref)s that represent a physical
 observable. Their ground state expectation values can be sampled by passing
 them into [`AllOverlaps`](@ref).
 ```@docs
+ParticleNumberOperator
 G2MomCorrelator
 G2RealCorrelator
 SuperfluidCorrelator

src/Hamiltonians/Hamiltonians.jl

@@ -36,6 +36,7 @@ Other
 - [`Stoquastic`](@ref)
 
 ## [Observables](#Observables)
+- [`ParticleNumberOperator`](@ref)
 - [`G2MomCorrelator`](@ref)
 - [`G2RealCorrelator`](@ref)
 - [`DensityMatrixDiagonal`](@ref)
@@ -83,6 +84,7 @@ export Stoquastic
 export Transcorrelated1D
 export hubbard_dispersion, continuum_dispersion
 export FroehlichPolaron
+export ParticleNumberOperator
 
 export G2MomCorrelator, G2RealCorrelator, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
 export StringCorrelator
@@ -125,6 +127,7 @@ include("Transcorrelated1D.jl")
 include("correlation_functions.jl")
 include("DensityMatrixDiagonal.jl")
 include("Momentum.jl")
+include("particle_number.jl")
 
 include("HOCartesianContactInteractions.jl")
 include("HOCartesianEnergyConservedPerDim.jl")

src/Hamiltonians/Momentum.jl

@@ -27,7 +27,6 @@ Momentum(component=0; fold=true) = Momentum{component}(fold)
 Base.show(io::IO, mom::Momentum{C}) where {C} = print(io, "Momentum($C; fold=$(mom.fold))")
 
 LOStructure(::Type{<:Momentum}) = IsDiagonal()
-num_offdiagonals(ham::Momentum, add) = 0
 
 @inline function _momentum(add::SingleComponentFockAddress, fold)
     M = num_modes(add)

src/Hamiltonians/particle_number.jl

@@ -0,0 +1,36 @@
+"""
+    ParticleNumberOperator([address]) <: AbstractHamiltonian
+
+The number operator in Fock space. This operator is diagonal in the Fock basis and
+returns the number of particles in the Fock state. Passing an address is optional.
+
+```jldoctest
+julia> h = FroehlichPolaron(fs"|0 0⟩{}"; mode_cutoff=5, v=3); bsr = BasisSetRep(h);
+
+julia> gs = DVec(zip(bsr.basis, eigen(Matrix(bsr)).vectors[:,1])); # ground state
+
+julia> dot(gs, ParticleNumberOperator(), gs) # particle number expectation value
+2.8823297252925917
+```
+
+See also [`AbstractHamiltonian`](@ref).
+"""
+struct ParticleNumberOperator{A} <: AbstractHamiltonian{Float64}
+    address::A
+end
+ParticleNumberOperator() = ParticleNumberOperator(BoseFS(1,))
+
+function Base.show(io::IO, n::ParticleNumberOperator)
+    io = IOContext(io, :compact => true)
+    print(io, "ParticleNumberOperator(")
+    n.address === BoseFS(1,) || show(io, n.address) # suppress if default
+    print(io, ")")
+end
+
+LOStructure(::Type{<:ParticleNumberOperator}) = IsDiagonal()
+starting_address(n::ParticleNumberOperator) = n.address
+
+function diagonal_element(::ParticleNumberOperator, addr::AbstractFockAddress)
+    return float(num_particles(addr))
+end
+allowed_address_type(::ParticleNumberOperator) = AbstractFockAddress

src/Interfaces/hamiltonians.jl

@@ -252,6 +252,12 @@ LOStructure(op) = LOStructure(typeof(op))
 LOStructure(::Type) = AdjointUnknown()
 LOStructure(::AbstractMatrix) = AdjointKnown()
 
+# diagonal matrices have zero offdiagonal elements
+function num_offdiagonals(h::H, addr) where {H<:AbstractHamiltonian}
+    return num_offdiagonals(LOStructure(H), h, addr)
+end
+num_offdiagonals(::IsDiagonal, _, _) = 0
+
 """
     has_adjoint(op)
 

test/Hamiltonians.jl

@@ -196,7 +196,9 @@ end
         HOCartesianCentralImpurity(BoseFS((1,0,0,0,0))),
 
         FroehlichPolaron(OccupationNumberFS(1,1,1)),
-        FroehlichPolaron(OccupationNumberFS(1,1,1); momentum_cutoff = 10.0)
+        FroehlichPolaron(OccupationNumberFS(1,1,1); momentum_cutoff = 10.0),
+
+        ParticleNumberOperator(OccupationNumberFS(1, 1, 1))
     )
         test_hamiltonian_interface(H)
     end