SuperfluidCorrelator and StringCorrelator

docs/src/hamiltonians.md

@@ -103,6 +103,8 @@ them into [`AllOverlaps`](@ref).
 ```@docs
 G2MomCorrelator
 G2RealCorrelator
+SuperfluidCorrelator
+StringCorrelator
 DensityMatrixDiagonal
 Momentum
 AxialAngularMomentumHO

src/Hamiltonians/Hamiltonians.jl

@@ -74,7 +74,8 @@ export Stoquastic
 export Transcorrelated1D
 export hubbard_dispersion, continuum_dispersion
 
-export G2MomCorrelator, G2RealCorrelator, DensityMatrixDiagonal, Momentum
+export G2MomCorrelator, G2RealCorrelator, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
+export StringCorrelator
 
 export LatticeGeometry, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
 export num_neighbours, neighbour_site, num_dimensions

src/Hamiltonians/correlation_functions.jl

@@ -236,3 +236,140 @@ function get_offdiagonal(
     gd = exp(-im*g.d*Δp*2π/M)*gamma
     return new_add, ComplexF64(gd/M)
 end
+
+"""
+    SuperfluidCorrelator(d::Int) <: AbstractHamiltonian{Float64}
+
+Operator for extracting superfluid correlation between sites separated by a distance `d` with `0 ≤ d < M`:
+
+```math
+    \\hat{C}_{\\text{SF}}(d) = \\frac{1}{M} \\sum_{i}^{M} a_{i}^{\\dagger} a_{i + d}
+```
+Assumes a one-dimensional lattice with ``M`` sites and periodic boundary conditions. ``M`` is also the number of modes in the Fock state address.
+
+# Usage
+Superfluid correlations can be extracted from a Monte Carlo calculation by wrapping `SuperfluidCorrelator` with 
+[`AllOverlaps`](@ref) and passing into [`lomc!`](@ref) with the `replica` keyword argument. For an example with a
+similar use of [`G2RealCorrelator`](@ref) see 
+[G2 Correlator Example](https://joachimbrand.github.io/Rimu.jl/previews/PR227/generated/G2-example.html).
+
+
+See also [`HubbardReal1D`](@ref), [`G2RealCorrelator`](@ref), [`AbstractHamiltonian`](@ref),
+and [`AllOverlaps`](@ref).
+"""
+struct SuperfluidCorrelator{D} <: AbstractHamiltonian{Float64}
+end
+
+SuperfluidCorrelator(d::Int) = SuperfluidCorrelator{d}()
+
+function Base.show(io::IO, ::SuperfluidCorrelator{D}) where {D}
+    print(io, "SuperfluidCorrelator($D)")
+end
+
+function num_offdiagonals(::SuperfluidCorrelator, add::SingleComponentFockAddress)
+    return num_occupied_modes(add)
+end
+
+function get_offdiagonal(::SuperfluidCorrelator{D}, add::SingleComponentFockAddress, chosen) where {D}
+    src = find_occupied_mode(add, chosen)
+    dst = find_mode(add, mod1(src.mode + D, num_modes(add)))
+    address, value = excitation(add, (dst,), (src,))
+    return address, value / num_modes(add)
+end
+
+function diagonal_element(::SuperfluidCorrelator{0}, add::SingleComponentFockAddress)
+    return num_particles(add) / num_modes(add)
+end
+function diagonal_element(::SuperfluidCorrelator{D}, add::SingleComponentFockAddress) where {D}
+    return 0.0
+end
+
+
+"""
+    StringCorrelator(d::Int) <: AbstractHamiltonian{Float64}
+
+Operator for extracting string correlation between lattice sites on a one-dimensional Hubbard lattice 
+separated by a distance `d` with `0 ≤ d < M`
+
+```math
+    \\hat{C}_{\\text{string}}(d) = \\frac{1}{M} \\sum_{j}^{M} \\delta n_j (e^{i \\pi \\sum_{j \\leq k < j + d} \\delta n_k}) \\delta n_{j+d}
+```
+Here, ``\\delta \\hat{n}_j = \\hat{n}_j - \\bar{n}`` is the boson number deviation from the mean filling 
+number and ``\\bar{n} = N/M`` is the mean filling number of lattice sites with ``N`` particles and 
+``M`` lattice sites (or modes). 
+
+Assumes a one-dimensional lattice with periodic boundary conditions. For usage
+see [`SuperfluidCorrelator`](@ref) and [`AllOverlaps`](@ref).
+
+See also [`HubbardReal1D`](@ref), [`G2RealCorrelator`](@ref), [`SuperfluidCorrelator`](@ref),
+[`AbstractHamiltonian`](@ref), and [`AllOverlaps`](@ref).
+"""
+struct StringCorrelator{D} <: AbstractHamiltonian{Float64}
+end
+
+StringCorrelator(d::Int) = StringCorrelator{d}()
+
+function Base.show(io::IO, ::StringCorrelator{D}) where {D}
+    print(io, "StringCorrelator($D)")
+end
+
+LOStructure(::Type{<:StringCorrelator}) = IsDiagonal()
+
+function diagonal_element(::StringCorrelator{0}, add::SingleComponentFockAddress)
+    M = num_modes(add)
+    N = num_particles(add)
+    n̄ = N/M
+    v = onr(add)
+
+    result = 0.0
+    for i in eachindex(v)
+        result += (v[i] - n̄)^2
+    end
+
+    return result / M
+end
+
+num_offdiagonals(::StringCorrelator, ::SingleComponentFockAddress) = 0
+
+function diagonal_element(::StringCorrelator{D}, add::SingleComponentFockAddress) where {D}
+    M = num_modes(add)
+    N = num_particles(add)
+    d = mod(D, M)
+
+    if !ismissing(N) && iszero(N % M)
+        return _string_diagonal_real(d, add)
+    else
+        return _string_diagonal_complex(d, add)
+    end
+end
+
+function _string_diagonal_complex(d, add)
+    M = num_modes(add)
+    N = num_particles(add)
+    n̄ = N/M
+    v = onr(add)
+
+    result = ComplexF64(0)
+    for i in eachindex(v)
+        phase_sum = sum((v[mod1(k, M)] - n̄) for k in i:1:(i+d-1))
+
+        result += (v[i] - n̄) * exp(pi * im * phase_sum) * (v[mod1(i + d, M)] - n̄)
+    end
+
+    return result / M
+end
+function _string_diagonal_real(d, add)
+    M = num_modes(add)
+    N = num_particles(add)
+    n̄ = N ÷ M
+    v = onr(add)
+
+    result = 0.0
+    for i in eachindex(v)
+        phase_sum = sum((v[mod1(k, M)] - n̄) for k in i:1:(i+d-1))
+
+        result += (v[i] - n̄) * (-1)^phase_sum * (v[mod1(i + d, M)] - n̄)
+    end
+
+    return result / M
+end

test/Hamiltonians.jl

@@ -4,6 +4,7 @@ using Random
 using Rimu
 using Test
 using DataFrames
+using Suppressor
 
 function exact_energy(ham)
     dv = DVec(starting_address(ham) => 1.0)
@@ -837,6 +838,87 @@ end
     # offdiagonals
     @test num_offdiagonals(G2RealCorrelator(0), add1) == 0
     @test num_offdiagonals(G2RealCorrelator(0), comp) == 0
+
+    # Test show method
+    d = 5
+    output = @capture_out print(G2RealCorrelator(d))
+    @test output == "G2RealCorrelator($d)"
+end
+
+@testset "SuperfluidCorrelator" begin
+    m = 6
+    n1 = 4
+    n2 = m
+    add1 = BoseFS((n1,0,0,0,0,0))
+    add2 = near_uniform(BoseFS{n2,m})
+
+    # localised state
+    @test @inferred diagonal_element(SuperfluidCorrelator(0), add1) == n1/m
+    @test @inferred diagonal_element(SuperfluidCorrelator(1), add1) == 0.
+
+    # constant density state
+    @test diagonal_element(SuperfluidCorrelator(0), add2) == n2/m
+    @test diagonal_element(SuperfluidCorrelator(1), add2) == 0.
+
+    # offdiagonals
+    @test num_offdiagonals(SuperfluidCorrelator(0), add1) == 1
+    @test num_offdiagonals(SuperfluidCorrelator(0), add2) == 6
+
+    # get_offdiagonal
+    @test get_offdiagonal(SuperfluidCorrelator(0), add1, 1) == (add1, n1/m)
+    @test get_offdiagonal(SuperfluidCorrelator(1), add1, 1) == (BoseFS((3,1,0,0,0,0)), sqrt(n1)/m)
+    @test get_offdiagonal(SuperfluidCorrelator(0), add2, 1) == (add2, 1/m)
+    @test get_offdiagonal(SuperfluidCorrelator(1), add2, 1) == (BoseFS((0,2,1,1,1,1)), sqrt(2)/m)
+
+    # Test show method
+    d = 5
+    output = @capture_out print(SuperfluidCorrelator(d))
+    @test output == "SuperfluidCorrelator($d)"
+end
+
+@testset "StringCorrelator" begin
+    m = 6
+    n1 = 4
+    n2 = m
+
+    # unital refers to n̄=1
+    non_unital_localised_state = BoseFS((n1,0,0,0,0,0))
+    non_unital_uniform_state = near_uniform(non_unital_localised_state)
+
+    localised_state = BoseFS((n2,0,0,0,0,0))
+    uniform_state = near_uniform(BoseFS{n2,m})
+
+    S0 = StringCorrelator(0)
+    S1 = StringCorrelator(1)
+    S2 = StringCorrelator(2)
+
+    @test num_offdiagonals(S0, localised_state) == 0
+
+    # non unital localised state
+    @test @inferred diagonal_element(S0, non_unital_localised_state) ≈ 20/9
+    @test @inferred diagonal_element(S1, non_unital_localised_state) ≈ (-4/9)*exp(im * -2pi/3)
+
+    # non unital near uniform state
+    @test @inferred diagonal_element(S0, non_unital_uniform_state) ≈ 2/9
+
+    # constant density localised state
+    @test @inferred diagonal_element(S0, localised_state) == 5.
+    @test @inferred diagonal_element(S1, localised_state) ≈ 1
+    @test @inferred diagonal_element(S2, localised_state) ≈ -1
+
+    # constant density uniform state
+    @test @inferred diagonal_element(S0, uniform_state) == 0
+    @test @inferred diagonal_element(S2, uniform_state) == 0
+
+    # Test return type for integer, and non-integer filling
+    @test @inferred diagonal_element(S0, localised_state) isa Float64
+    @test @inferred diagonal_element(S1, non_unital_localised_state) isa ComplexF64
+
+    # Test show method
+    d = 5
+    output = @capture_out print(StringCorrelator(d))
+    @test output == "StringCorrelator($d)"
+
 end
 
 @testset "Momentum" begin