diff --git a/docs/src/hamiltonians.md b/docs/src/hamiltonians.md
index e674ba248..e79c02316 100644
--- a/docs/src/hamiltonians.md
+++ b/docs/src/hamiltonians.md
@@ -103,8 +103,9 @@ observable. Their ground state expectation values can be sampled by passing
 them into [`AllOverlaps`](@ref).
 ```@docs
 ParticleNumberOperator
-G2MomCorrelator
 G2RealCorrelator
+G2RealSpace
+G2MomCorrelator
 SuperfluidCorrelator
 StringCorrelator
 DensityMatrixDiagonal
@@ -156,20 +157,21 @@ Hamiltonians.number_conserving_fermi_dimension
 ```
 
 ## Geometry
-Lattices in higher dimensions are defined here for [`HubbardRealSpace`](@ref).
+
+Lattices in higher dimensions are defined here for [`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref).
 
 ```@docs
-LatticeGeometry
+CubicGrid
+Hamiltonians.Directions
+Hamiltonians.Displacements
+Hamiltonians.neighbor_site
 PeriodicBoundaries
 HardwallBoundaries
 LadderBoundaries
-num_neighbours
-num_dimensions
-neighbour_site
 ```
 
 ## Harmonic Oscillator
-Useful utilities for harmonic oscillator in Cartesian basis, see [`HOCartesianContactInteractions`](@ref) 
+Useful utilities for harmonic oscillator in Cartesian basis, see [`HOCartesianContactInteractions`](@ref)
 and [`HOCartesianEnergyConservedPerDim`](@ref).
 ```@docs
 get_all_blocks
diff --git a/src/DictVectors/abstractdvec.jl b/src/DictVectors/abstractdvec.jl
index 6065af5dd..085e83c44 100644
--- a/src/DictVectors/abstractdvec.jl
+++ b/src/DictVectors/abstractdvec.jl
@@ -296,7 +296,8 @@ Internal function evaluates the 3-argument `dot()` function in order from right
 to left.
 """
 function dot_from_right(w, op, v::AbstractDVec)
-    result = zero(promote_type(valtype(w), eltype(op), valtype(v)))
+    T = typeof(zero(valtype(w)) * zero(eltype(op)) * zero(valtype(v)))
+    result = zero(T)
     for (key, val) in pairs(v)
         result += conj(w[key]) * diagonal_element(op, key) * val
         for (add, elem) in offdiagonals(op, key)
diff --git a/src/DictVectors/communicators.jl b/src/DictVectors/communicators.jl
index 62273dccc..a702fc75c 100644
--- a/src/DictVectors/communicators.jl
+++ b/src/DictVectors/communicators.jl
@@ -47,7 +47,7 @@ is_distributed(::Communicator) = true
 
 Merge the results of reductions over MPI. By default, it uses `MPI.Allreduce`.
 """
-merge_remote_reductions(c::Communicator, op, x) = MPI.Allreduce(x, op, mpi_comm(c))
+merge_remote_reductions(c::Communicator, op, x) = MPI.Allreduce!(Ref(x), op, mpi_comm(c))[]
 
 """
     total_num_segments(c::Communicator, n) -> Int
diff --git a/src/DictVectors/pdvec.jl b/src/DictVectors/pdvec.jl
index 82015a7dc..e4ad2620f 100644
--- a/src/DictVectors/pdvec.jl
+++ b/src/DictVectors/pdvec.jl
@@ -769,7 +769,7 @@ end
 function LinearAlgebra.dot(
     ::IsDiagonal, t::PDVec, op::AbstractHamiltonian, u::PDVec, w=nothing
 )
-    T = promote_type(valtype(t), eltype(op), valtype(u))
+    T = typeof(zero(valtype(t)) * zero(valtype(u)) * zero(eltype(op)))
     return sum(pairs(u); init=zero(T)) do (k, v)
         T(conj(t[k]) * diagonal_element(op, k) * v)
     end
@@ -809,7 +809,8 @@ function LinearAlgebra.dot(
 end
 
 function dot_from_right(target, op, source::PDVec)
-    T = promote_type(valtype(target), valtype(source), eltype(op))
+    T = typeof(zero(valtype(target)) * zero(valtype(source)) * zero(eltype(op)))
+
     result = sum(pairs(source); init=zero(T)) do (k, v)
         res = conj(target[k]) * diagonal_element(op, k) * v
         for (k_off, v_off) in offdiagonals(op, k)
diff --git a/src/Hamiltonians/GutzwillerSampling.jl b/src/Hamiltonians/GutzwillerSampling.jl
index 411679023..b5c88cf80 100644
--- a/src/Hamiltonians/GutzwillerSampling.jl
+++ b/src/Hamiltonians/GutzwillerSampling.jl
@@ -135,7 +135,7 @@ function TransformUndoer(k::GutzwillerSampling, op::Union{Nothing,AbstractHamilt
     if isnothing(op)
         T = eltype(k)
     else
-        T = promote_type(eltype(k), eltype(op))
+        T = typeof(zero(eltype(k)) * zero(eltype(op)))
     end
     return TransformUndoer{T,typeof(k),typeof(op)}(k, op)
 end
diff --git a/src/Hamiltonians/Hamiltonians.jl b/src/Hamiltonians/Hamiltonians.jl
index 022b6b533..b2ee064d6 100644
--- a/src/Hamiltonians/Hamiltonians.jl
+++ b/src/Hamiltonians/Hamiltonians.jl
@@ -37,8 +37,9 @@ Other
 
 ## [Observables](#Observables)
 - [`ParticleNumberOperator`](@ref)
-- [`G2MomCorrelator`](@ref)
 - [`G2RealCorrelator`](@ref)
+- [`G2RealSpace`](@ref)
+- [`G2MomCorrelator`](@ref)
 - [`DensityMatrixDiagonal`](@ref)
 - [`Momentum`](@ref)
 - [`AxialAngularMomentumHO`](@ref)
@@ -59,7 +60,7 @@ using Parameters: Parameters, @unpack
 using Setfield: Setfield
 using SparseArrays: SparseArrays, nnz, nzrange, sparse
 using SpecialFunctions: SpecialFunctions, gamma
-using StaticArrays: StaticArrays, SA, SMatrix, SVector
+using StaticArrays: StaticArrays, SA, SMatrix, SVector, SArray, setindex
 using TupleTools: TupleTools
 
 using ..BitStringAddresses
@@ -86,11 +87,10 @@ export hubbard_dispersion, continuum_dispersion
 export FroehlichPolaron
 export ParticleNumberOperator
 
-export G2MomCorrelator, G2RealCorrelator, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
+export G2MomCorrelator, G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagonal, Momentum
 export StringCorrelator
 
-export LatticeGeometry, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
-export num_neighbours, neighbour_site, num_dimensions
+export CubicGrid, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
 
 export sparse # from SparseArrays
 
diff --git a/src/Hamiltonians/HubbardRealSpace.jl b/src/Hamiltonians/HubbardRealSpace.jl
index f3c156062..9e25c64a6 100644
--- a/src/Hamiltonians/HubbardRealSpace.jl
+++ b/src/Hamiltonians/HubbardRealSpace.jl
@@ -18,10 +18,10 @@ See also [`BoseFS`](@ref), [`FermiFS`](@ref), [`CompositeFS`](@ref).
 local_interaction(b::SingleComponentFockAddress, u) = u * bose_hubbard_interaction(b) / 2
 local_interaction(f::FermiFS, _) = 0
 function local_interaction(a::SingleComponentFockAddress, b::SingleComponentFockAddress, u)
-    u * dot(occupied_modes(a), occupied_modes(b))
+    return u * dot(occupied_modes(a), occupied_modes(b))
 end
 function local_interaction(fs::CompositeFS, u)
-    _interactions(fs.components, u)
+    return _interactions(fs.components, u)
 end
 
 """
@@ -140,7 +140,7 @@ is produced if `address`is incompatible with the interaction parameters `u`.
 
 ## Geometries
 
-Implemented [`LatticeGeometry`](@ref)s for keyword `geometry`
+Implemented [`CubicGrid`](@ref)s for keyword `geometry`
 
 * [`PeriodicBoundaries`](@ref)
 * [`HardwallBoundaries`](@ref)
@@ -163,7 +163,7 @@ number of sites `M` inferred from the number of modes in `address`.
 struct HubbardRealSpace{
     C, # components
     A<:AbstractFockAddress,
-    G<:LatticeGeometry,
+    G<:CubicGrid,
     D, # dimension
     # The following need to be type params.
     T<:SVector{C,Float64},
@@ -181,7 +181,7 @@ end
 
 function HubbardRealSpace(
     address::AbstractFockAddress;
-    geometry::LatticeGeometry=PeriodicBoundaries((num_modes(address),)),
+    geometry::CubicGrid=PeriodicBoundaries((num_modes(address),)),
     t=ones(num_components(address)),
     u=ones(num_components(address), num_components(address)),
     v=zeros(num_components(address), num_dimensions(geometry))
@@ -313,7 +313,7 @@ struct HubbardRealSpaceCompOffdiagonals{G,A} <: AbstractOffdiagonals{A,Float64}
 end
 
 function offdiagonals(h::HubbardRealSpace, comp, add)
-    neighbours = num_neighbours(h.geometry)
+    neighbours = 2 * num_dimensions(h.geometry)
     return HubbardRealSpaceCompOffdiagonals(
         h.geometry, add, h.t[comp], num_occupied_modes(add) * neighbours
     )
@@ -322,10 +322,10 @@ end
 Base.size(o::HubbardRealSpaceCompOffdiagonals) = (o.length,)
 
 @inline function Base.getindex(o::HubbardRealSpaceCompOffdiagonals, chosen)
-    neighbours = num_neighbours(o.geometry)
+    neighbours = 2 * num_dimensions(o.geometry)
     particle, neigh = fldmod1(chosen, neighbours)
     src_index = find_occupied_mode(o.address, particle)
-    neigh = neighbour_site(o.geometry, src_index.mode, neigh)
+    neigh = neighbor_site(o.geometry, src_index.mode, neigh)
 
     if neigh == 0
         return o.address, 0.0
diff --git a/src/Hamiltonians/correlation_functions.jl b/src/Hamiltonians/correlation_functions.jl
index 0f4392065..2d7699164 100644
--- a/src/Hamiltonians/correlation_functions.jl
+++ b/src/Hamiltonians/correlation_functions.jl
@@ -1,4 +1,36 @@
-# TO-DO: add geometry for higher dimensions.
+"""
+    circshift_dot(arr1, arr2, inds)
+
+Fast, non-allocation version of
+
+```julia
+dot(arr1, circshift(arr2, inds))
+```
+"""
+function circshift_dot(arr1, arr2, inds)
+    _circshift_dot!(arr1, (), arr2, (), axes(arr2), Tuple(inds))
+end
+
+# The following is taken from Julia's implementation of circshift.
+@inline function _circshift_dot!(
+    dst, rdest, src, rsrc,
+    inds::Tuple{AbstractUnitRange,Vararg{Any}},
+    shiftamt::Tuple{Integer,Vararg{Any}}
+)
+    ind1, d = inds[1], shiftamt[1]
+    s = mod(d, length(ind1))
+    sf, sl = first(ind1)+s, last(ind1)-s
+    r1, r2 = first(ind1):sf-1, sf:last(ind1)
+    r3, r4 = first(ind1):sl, sl+1:last(ind1)
+    tinds, tshiftamt = Base.tail(inds), Base.tail(shiftamt)
+
+    return _circshift_dot!(dst, (rdest..., r1), src, (rsrc..., r4), tinds, tshiftamt) +
+        _circshift_dot!(dst, (rdest..., r2), src, (rsrc..., r3), tinds, tshiftamt)
+end
+@inline function _circshift_dot!(dst, rdest, src, rsrc, inds, shiftamt)
+    return dot(view(dst, rdest...), view(src, rsrc...))
+end
+
 """
     G2RealCorrelator(d::Int) <: AbstractHamiltonian{Float64}
 
@@ -28,6 +60,7 @@ equivalent to stacking all components into a single Fock state.
 # See also
 
 * [`HubbardReal1D`](@ref)
+* [`G2RealSpace`](@ref)
 * [`G2MomCorrelator`](@ref)
 * [`AbstractHamiltonian`](@ref)
 * [`AllOverlaps`](@ref)
@@ -52,11 +85,7 @@ function diagonal_element(::G2RealCorrelator{D}, add::SingleComponentFockAddress
     M = num_modes(add)
     d = mod(D, M)
     v = onr(add)
-    result = 0
-    for i in eachindex(v)
-        result += v[i] * v[mod1(i + d, M)]
-    end
-    return result / M
+    return circshift_dot(v, v, (d,)) / M
 end
 
 function diagonal_element(::G2RealCorrelator{0}, add::CompositeFS)
@@ -68,22 +97,133 @@ function diagonal_element(::G2RealCorrelator{D}, add::CompositeFS) where {D}
     M = num_modes(add)
     d = mod(D, M)
     v = sum(map(onr, add.components))
-    result = 0
-    for i in eachindex(v)
-        result += v[i] * v[mod1(i + d, M)]
-    end
-    return result / M
+    return circshift_dot(v, v, (d,)) / M
 end
 
 num_offdiagonals(::G2RealCorrelator, ::SingleComponentFockAddress) = 0
 num_offdiagonals(::G2RealCorrelator, ::CompositeFS) = 0
 
-# not needed:
-# get_offdiagonal(::G2RealCorrelator, add)
-# starting_address(::G2RealCorrelator)
+"""
+    G2RealSpace(geometry::CubicGrid, σ=1, τ=1; sum_components=false) <: AbstractHamiltonian{SArray}
+
+Two-body operator for density-density correlation for all [`Displacements`](@ref) ``d⃗`` in the specified `geometry`.
+
+```math
+    \\hat{G}^{(2)}_{σ,τ}(d⃗) = \\frac{1}{M} ∑_{i⃗} n̂_{σ,i⃗} (n̂_{τ,i⃗+d⃗} - δ_{0⃗,d⃗}δ_{σ,τ}).
+```
+
+For multicomponent addresses, `σ` and `τ` control the components involved. Alternatively,
+`sum_components` can be set to `true`, which treats all particles as belonging to the same
+component.
+
+# Examples
+
+```jldoctest
+julia> geom = CubicGrid(2, 2);
+
+julia> g2 = G2RealSpace(geom)
+G2RealSpace(CubicGrid((2, 2), (true, true)), 1,1)
+
+julia> diagonal_element(g2, BoseFS(2,0,1,1))
+2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
+ 0.5  1.0
+ 0.5  1.0
+
+julia> g2_cross = G2RealSpace(geom, 1, 2)
+G2RealSpace(CubicGrid((2, 2), (true, true)), 1,2)
+
+julia> g2_sum = G2RealSpace(geom, sum_components=true)
+G2RealSpace(CubicGrid((2, 2), (true, true)); sum_components=true)
+
+julia> diagonal_element(g2, fs"|⇅⋅↓↑⟩")
+2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
+ 0.0  0.0
+ 0.0  0.5
+
+julia> diagonal_element(g2_cross, fs"|⇅⋅↓↑⟩")
+2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
+ 0.25  0.25
+ 0.25  0.25
+
+julia> diagonal_element(g2_sum, fs"|⇅⋅↓↑⟩")
+2×2 StaticArraysCore.SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
+ 0.5  1.0
+ 0.5  1.0
+```
+
+# See also
+
+* [`CubicGrid`](@ref)
+* [`HubbardRealSpace`](@ref)
+* [`G2RealCorrelator`](@ref)
+* [`G2MomCorrelator`](@ref)
+* [`AbstractHamiltonian`](@ref)
+* [`AllOverlaps`](@ref)
+"""
+struct G2RealSpace{A,B,G<:CubicGrid,S} <: AbstractHamiltonian{S}
+    geometry::G
+    init::S
+end
+function G2RealSpace(geometry::CubicGrid, σ::Int=1, τ::Int=σ; sum_components=false)
+    if σ < 1 || τ < 1
+        throw(ArgumentError("`σ` and `τ` must be positive integers"))
+    end
+    if sum_components
+        if σ ≠ 1 || τ ≠ 1
+            throw(ArgumentError("`σ` or `τ` can't be set if `sum_components=true`"))
+        end
+        σ = τ = 0
+    end
+
+    init = zeros(SArray{Tuple{size(geometry)...}})
+    return G2RealSpace{σ,τ,typeof(geometry),typeof(init)}(geometry, init)
+end
+
+function Base.show(io::IO, g2::G2RealSpace{A,B}) where {A,B}
+    print(io, "G2RealSpace($(g2.geometry), $A,$B)")
+end
+function Base.show(io::IO, g2::G2RealSpace{0,0})
+    print(io, "G2RealSpace($(g2.geometry); sum_components=true)")
+end
+
+LOStructure(::Type{<:G2RealSpace}) = IsDiagonal()
+
+num_offdiagonals(g2::G2RealSpace, _) = 0
+
+@inline function _g2_diagonal_element(
+    g2::G2RealSpace{A,B}, onr1::SArray, onr2::SArray
+) where {A, B}
+    geo = g2.geometry
+    result = g2.init
+
+    @inbounds for i in eachindex(result)
+        res_i = 0.0
+        displacement = Displacements(geo)[i]
+
+        # Case of n_i(n_i - 1) on the same component
+        if A == B && iszero(displacement)
+            onr1_minus_1 = max.(onr1 .- 1, 0)
+            result = setindex(result, dot(onr2, onr1_minus_1), i)
+        else
+            result = setindex(result, circshift_dot(onr2, onr1, displacement), i)
+        end
+    end
+    return result ./ length(geo)
+end
+function diagonal_element(g2::G2RealSpace{A,A}, addr::SingleComponentFockAddress) where {A}
+    onr1 = onr(addr, g2.geometry)
+    return _g2_diagonal_element(g2, onr1, onr1)
+end
+function diagonal_element(g2::G2RealSpace{A,B}, addr::CompositeFS) where {A,B}
+    onr1 = onr(addr.components[A], g2.geometry)
+    onr2 = onr(addr.components[B], g2.geometry)
+    return _g2_diagonal_element(g2, onr1, onr2)
+end
+function diagonal_element(g2::G2RealSpace{0,0}, addr::CompositeFS)
+    onr1 = sum(x -> onr(x, g2.geometry), addr.components)
+    return _g2_diagonal_element(g2, onr1, onr1)
+end
 
-# Methods that need to be implemented:
-#   •  starting_address(::AbstractHamiltonian) - not needed
 """
     G2MomCorrelator(d::Int,c=:cross) <: AbstractHamiltonian{ComplexF64}
 
@@ -125,6 +265,7 @@ and let it be the default value.
 * [`BoseHubbardMom1D2C`](@ref)
 * [`BoseFS2C`](@ref)
 * [`G2RealCorrelator`](@ref)
+* [`G2RealSpace`](@ref)
 * [`AbstractHamiltonian`](@ref)
 * [`AllOverlaps`](@ref)
 """
@@ -165,7 +306,6 @@ function num_offdiagonals(g::G2MomCorrelator{3}, add::BoseFS2C)
     sa = num_occupied_modes(add.bsa)
     sb = num_occupied_modes(add.bsb)
     return sa*(m-1)*sb
-    # number of excitations that can be made
 end
 
 function num_offdiagonals(g::G2MomCorrelator, add::SingleComponentFockAddress)
@@ -251,9 +391,9 @@ Operator for extracting superfluid correlation between sites separated by a dist
 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 
+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 
+similar use of [`G2RealCorrelator`](@ref) see
 [G2 Correlator Example](https://joachimbrand.github.io/Rimu.jl/previews/PR227/generated/G2-example.html).
 
 
@@ -291,15 +431,15 @@ end
 """
     StringCorrelator(d::Int) <: AbstractHamiltonian{Float64}
 
-Operator for extracting string correlation between lattice sites on a one-dimensional Hubbard lattice 
+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). 
+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).
diff --git a/src/Hamiltonians/geometry.jl b/src/Hamiltonians/geometry.jl
index 0de29feb5..5bef66d14 100644
--- a/src/Hamiltonians/geometry.jl
+++ b/src/Hamiltonians/geometry.jl
@@ -1,256 +1,236 @@
 """
-    abstract type LatticeGeometry{D}
+    CubicGrid(dims::NTuple{D,Int}, fold::NTuple{D,Bool})
 
-A `LatticeGeometry` controls which sites in an [`AbstractFockAddress`](@ref) are considered
-to be neighbours.
+Represents a `D`-dimensional grid. Used to define a cubic lattice and boundary conditions
+for some [`AbstractHamiltonian`](@ref)s. The type instance can be used to convert between
+cartesian vector indices (tuples or `SVector`s) and linear indices (integers). When indexed
+with vectors, it folds them back into the grid if the out-of-bounds dimension is periodic and
+0 otherwise (see example below).
 
-Currently only supported by [`HubbardRealSpace`](@ref).
+* `dims` controls the size of the grid in each dimension.
+* `fold` controls whether the boundaries in each dimension are periodic (or folded in the
+  case of momentum space).
 
-# Available implementations
+```julia
+julia> geo = CubicGrid((2,3), (true,false))
+CubicGrid{2}((2, 3), (true, false))
 
-* [`PeriodicBoundaries`](@ref)
-* [`HardwallBoundaries`](@ref)
-* [`LadderBoundaries`](@ref)
+julia> geo[1]
+(1, 1)
 
-# Interface to implement
+julia> geo[2]
+(2, 1)
 
-* `Base.size`: return the lattice size.
-* [`neighbour_site(::LatticeGeometry, ::Int, ::Int)`](@ref)
-* [`num_dimensions(::LatticeGeometry)`](@ref)
-* [`num_neighbours(::LatticeGeometry)`](@ref)
-"""
-abstract type LatticeGeometry{D} end
+julia> geo[3]
+(1, 2)
+
+julia> geo[(1,2)]
+3
+
+julia> geo[(3,2)] # 3 is folded back into 1
+3
+
+julia> geo[(3,3)]
+5
+
+julia> geo[(3,4)] # returns 0 if out of bounds
+0
+```
 
-function Base.show(io::IO, geom::LatticeGeometry)
-    print(io, nameof(typeof(geom)), size(geom))
+See also [`PeriodicBoundaries`](@ref), [`HardwallBoundaries`](@ref) and
+[`LadderBoundaries`](@ref) for special-case constructors.
+"""
+struct CubicGrid{D,Dims,Fold}
+    function CubicGrid(
+        dims::NTuple{D,Int}, fold::NTuple{D,Bool}=ntuple(Returns(true), Val(D))
+    ) where {D}
+        if any(≤(1), dims)
+            throw(ArgumentError("All dimensions must be at least 2 in size"))
+        end
+        return new{D,dims,fold}()
+    end
 end
+CubicGrid(args::Vararg{Int}) = CubicGrid(args)
 
 """
-    onr(add::AbstractFockAddress, geom::LatticeGeometry)
-Returns the occupation number representation of a Fock state address as an `SArray`
-with the shape of the lattice geometry `geom`. For composite addresses, a tuple
-of `onr`s is returned.
+    PeriodicBoundaries(dims...) -> CubicGrid
+    PeriodicBoundaries(dims) -> CubicGrid
+
+Return `CubicGrid` with all dimensions periodic. Equivalent to `CubicGrid(dims)`.
 """
-function BitStringAddresses.onr(add, geom::LatticeGeometry)
-    return reshape(onr(add), size(geom))
-end
-function BitStringAddresses.onr(add::CompositeFS, geom::LatticeGeometry)
-    return map(fs -> onr(fs, geom), add.components)
+function PeriodicBoundaries(dims::NTuple{D,Int}) where {D}
+    return CubicGrid(dims, ntuple(Returns(true), Val(D)))
 end
+PeriodicBoundaries(dims::Vararg{Int}) = PeriodicBoundaries(dims)
 
 """
-    neighbour_site(geom::LatticeGeometry, site, i)
+    HardwallBoundaries(dims...) -> CubicGrid
+    HardwallBoundaries(dims) -> CubicGrid
 
-Find the `i`-th neighbour of `site` in the geometry. If the move is illegal, return 0.
+Return `CubicGrid` with all dimensions non-periodic. Equivalent to
+`CubicGrid(dims, (false, false, ...))`.
 """
-neighbour_site
+function HardwallBoundaries(dims::NTuple{D,Int}) where {D}
+    return CubicGrid(dims, ntuple(Returns(false), Val(D)))
+end
+HardwallBoundaries(dims::Vararg{Int}) = HardwallBoundaries(dims)
 
 """
-    num_dimensions(geom::LatticeGeometry)
+    LadderBoundaries(dims...) -> CubicGrid
+    LadderBoundaries(dims) -> CubicGrid
 
-Return the number of dimensions of the lattice in this geometry.
+Return `CubicGrid` where the first dimension is dimensions non-periodic and the rest are
+periodic. Equivalent to `CubicGrid(dims, (true, false, ...))`.
 """
-num_dimensions(::LatticeGeometry{D}) where D = D
+function LadderBoundaries(dims::NTuple{D,Int}) where {D}
+    return CubicGrid(dims, ntuple(>(1), Val(D)))
+end
+LadderBoundaries(dims::Vararg{Int}) = LadderBoundaries(dims)
 
-"""
-    num_neighbours(geom::LatticeGeometry)
+function Base.show(io::IO, g::CubicGrid{<:Any,Dims,Fold}) where {Dims,Fold}
+    print(io, "CubicGrid($Dims, $Fold)")
+end
 
-Return the number of neighbours each lattice site has in this geometry.
+Base.size(g::CubicGrid{<:Any,Dims}) where {Dims} = Dims
+Base.size(g::CubicGrid{<:Any,Dims}, i) where {Dims} = Dims[i]
+Base.length(g::CubicGrid) = prod(size(g))
+fold(g::CubicGrid{<:Any,<:Any,Fold}) where {Fold} = Fold
 
-Note that for efficiency reasons, all sites are expected to have the same number of
-neighbours. If some of the neighbours are invalid, this is handled by having
-[`neighbour_site`](@ref) return 0.
 """
-num_neighbours
-
+    num_dimensions(geom::LatticeCubicGrid)
 
+Return the number of dimensions of the lattice in this geometry.
 """
-    PeriodicBoundaries(size...) <: LatticeGeometry
-
-Rectangular lattice with periodic boundary conditions of size `size`.
-
-The dimension of the lattice is controlled by the number of arguments given to its
-constructor.
-
-This is the default geometry used by [`HubbardRealSpace`](@ref).
-
-## Example
+num_dimensions(::CubicGrid{D}) where {D} = D
 
-```
-julia> lattice = PeriodicBoundaries(5, 4) # 2D lattice of size 5 × 4
-PeriodicBoundaries(5, 4)
-
-julia> num_neighbours(lattice)
-4
+"""
+    fold_vec(g::CubicGrid{D}, vec::SVector{D,Int}) -> SVector{D,Int}
 
-julia> neighbour_site(lattice, 1, 1)
-2
+Use the CubicGrid to fold the `vec` in each dimension. If folding is disabled in a
+dimension, and the vector is allowed to go out of bounds.
 
-julia> neighbour_site(lattice, 1, 2)
-5
+```julia
+julia> geo = CubicGrid((2,3), (true,false))
+CubicGrid{2}((2, 3), (true, false))
 
-julia> neighbour_site(lattice, 1, 3)
-6
+julia> fold_vec(geo, (3,1))
+(1, 1)
 
-julia> neighbour_site(lattice, 1, 4)
-16
+julia> fold_vec(geo, (3,4))
+(1, 4)
 ```
-
-## See also
-
-* [`LatticeGeometry`](@ref)
-* [`HardwallBoundaries`](@ref)
-* [`LadderBoundaries`](@ref)
-* [`num_neighbours`](@ref)
-* [`neighbour_site`](@ref)
 """
-struct PeriodicBoundaries{D} <: LatticeGeometry{D}
-    size::NTuple{D,Int}
+function fold_vec(g::CubicGrid{D}, vec::SVector{D,Int}) where {D}
+    (_fold_vec(Tuple(vec), fold(g), size(g)))
+end
+@inline _fold_vec(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()
+@inline function _fold_vec((x, xs...), (f, fs...), (d, ds...))
+    x = f ? mod1(x, d) : x
+    return (x, _fold_vec(xs, fs, ds)...)
 end
 
-PeriodicBoundaries(args::Vararg{Int}) = PeriodicBoundaries(args)
-
-Base.size(geom::PeriodicBoundaries) = geom.size
-
-num_neighbours(::PeriodicBoundaries{D}) where {D} = 2D
-
-function neighbour_site(geom::PeriodicBoundaries{D}, site, i) where {D}
-    # Neighbour indexing
-    # i |  x  y  z …
-    # 0 | +1  0  0
-    # 1 | -1  0  0
-    # 2 |  0 +1  0
-    # 3 |  0 -1  0
-    # 4 |  0  0 +1
-    #         ⋮    ⋱
-    i -= 1
-    cart_indices = CartesianIndices(size(geom))
-    cart_index = Tuple(cart_indices[site])
-    offset = ntuple(Val(D)) do k
-        ifelse(2(k - 1) ≤ i < 2k, ifelse(iseven(i), 1, -1), 0)
-    end
-    new_index = CartesianIndex(mod1.(cart_index .+ offset, size(geom) .% UInt))
-    return LinearIndices(cart_indices)[new_index]
+function Base.getindex(g::CubicGrid{D}, vec::Union{NTuple{D,Int},SVector{D,Int}}) where {D}
+    return get(LinearIndices(size(g)), fold_vec(g, SVector(vec)), 0)
 end
+Base.getindex(g::CubicGrid, i::Int) = SVector(Tuple(CartesianIndices(size(g))[i]))
 
 """
-    HardwallBoundaries
+    Directions(D) <: AbstractVector{SVector{D,Int}}
+    Directions(geometry::CubicGrid) <: AbstractVector{SVector{D,Int}}
 
-Rectangular lattice with hard wall boundary conditions of size `size`.
-[`neighbour_site()`](@ref) will return 0 for some neighbours of boundary sites.
+Iterate over axis-aligned direction vectors in `D` dimensions.
 
-The dimension of the lattice is controlled by the number of arguments given to its
-constructor.
-
-## Example
-
-```
-julia> lattice = HardwallBoundaries(5) # 1D lattice of size 5
-HardwallBoundaries(5)
-
-julia> neighbour_site(lattice, 1, 1)
-2
-
-julia> neighbour_site(lattice, 1, 2)
-0
+```jldoctest; setup=:(using Rimu.Hamiltonians: Directions)
+julia> Directions(3)
+6-element Directions{3}:
+ [1, 0, 0]
+ [0, 1, 0]
+ [0, 0, 1]
+ [-1, 0, 0]
+ [0, -1, 0]
+ [0, 0, -1]
 
 ```
 
-## See also
-
-* [`LatticeGeometry`](@ref)
-* [`PeriodicBoundaries`](@ref)
-* [`LadderBoundaries`](@ref)
+See also [`CubicGrid`](@ref).
 """
-struct HardwallBoundaries{D} <: LatticeGeometry{D}
-    size::NTuple{D,Int}
-end
+struct Directions{D} <: AbstractVector{SVector{D,Int}} end
 
-HardwallBoundaries(args::Vararg{Int}) = HardwallBoundaries(args)
+Directions(D) = Directions{D}()
+Directions(::CubicGrid{D}) where {D} = Directions{D}()
 
-Base.size(geom::HardwallBoundaries) = geom.size
+Base.size(::Directions{D}) where {D} = (2D,)
 
-num_neighbours(::HardwallBoundaries{D}) where {D} = 2D
-
-function neighbour_site(geom::HardwallBoundaries{D}, site, i) where {D}
-    i -= 1
-    cart_indices = CartesianIndices(size(geom))
-    cart_index = Tuple(cart_indices[site])
-    offset = ntuple(Val(D)) do k
-        ifelse(2(k - 1) ≤ i < 2k, ifelse(iseven(i), 1, -1), 0)
-    end
-    new_index = CartesianIndex(cart_index .+ offset)
-    if new_index in cart_indices
-        return LinearIndices(cart_indices)[new_index]
+function Base.getindex(uv::Directions{D}, i) where {D}
+    @boundscheck 0 < i ≤ length(uv) || throw(BoundsError(uv, i))
+    if i ≤ D
+        return SVector(_unit_vec(Val(D), i, 1))
     else
-        return 0
+        return SVector(_unit_vec(Val(D), i - D, -1))
     end
 end
 
-"""
-    LadderBoundaries(size...; subgeometry=PeriodicBoundaries) <: LatticeGeometry
-
-Lattice geometry where the first dimension is of size 2 and has hardwall boundary conditions.
-Using this geometry is more efficient than using [`HardwallBoundaries`](@ref) with a size of
-2, as it does not generate rejected neighbours.
-
-In other dimensions, it behaves like its subgeometry, which can be any
-[`LatticeGeometry`](@ref).
-
-## Example
-
-```
-julia> lattice = LadderBoundaries(2, 3, 4) # 3D lattice of size 2 × 3 × 4
-LadderBoundaries(2, 3, 4)
+@inline _unit_vec(::Val{0}, _, _) = ()
+@inline function _unit_vec(::Val{I}, i, x) where {I}
+    val = ifelse(i == I, x, 0)
+    return (_unit_vec(Val(I-1), i, x)..., val)
+end
 
-julia> num_neighbours(lattice)
-5
+"""
+    Displacements(geometry::CubicGrid) <: AbstractVector{SVector{D,Int}}
 
-julia> neighbour_site(lattice, 1, 1)
-2
+Return all valid offset vectors in a [`CubicGrid`](@ref). If `center=true` the (0,0) displacement is
+placed at the centre of the array.
 
-julia> neighbour_site(lattice, 1, 2)
-3
+```jldoctest; setup=:(using Rimu.Hamiltonians: Displacements)
+julia> geometry = CubicGrid((3,4));
 
-julia> neighbour_site(lattice, 1, 3)
-5
+julia> reshape(Displacements(geometry), (3,4))
+3×4 reshape(::Displacements{2}, 3, 4) with eltype StaticArraysCore.SVector{2, Int64}:
+ [0, 0]  [0, 1]  [0, 2]  [0, 3]
+ [1, 0]  [1, 1]  [1, 2]  [1, 3]
+ [2, 0]  [2, 1]  [2, 2]  [2, 3]
 
-julia> neighbour_site(lattice, 1, 4)
-7
+julia> reshape(Displacements(geometry; center=true), (3,4))
+3×4 reshape(::Displacements{2}, 3, 4) with eltype StaticArraysCore.SVector{2, Int64}:
+ [-1, -1]  [-1, 0]  [-1, 1]  [-1, 2]
+ [0, -1]   [0, 0]   [0, 1]   [0, 2]
+ [1, -1]   [1, 0]   [1, 1]   [1, 2]
 
-julia> neighbour_site(lattice, 1, 5)
-19
 ```
-
-## See also
-
-* [`LatticeGeometry`](@ref)
-* [`PeriodicBoundaries`](@ref)
-* [`HardwallBoundaries`](@ref)
 """
-struct LadderBoundaries{D,G<:LatticeGeometry} <: LatticeGeometry{D}
-    subgeometry::G
+struct Displacements{D} <: AbstractVector{SVector{D,Int}}
+    geometry::CubicGrid{D}
+    center::Bool
 end
+Displacements(geometry; center=false) = Displacements(geometry, center)
+
+Base.size(off::Displacements) = (length(off.geometry),)
 
-function LadderBoundaries(arg::Int, args::Vararg{Int}; subgeometry=PeriodicBoundaries)
-    if arg ≠ 2
-        throw(ArgumentError("First dimension must be of size 2"))
+@inline function Base.getindex(off::Displacements{D}, i) where {D}
+    @boundscheck 0 < i ≤ length(off) || throw(BoundsError(off, i))
+    geo = off.geometry
+    vec = geo[i]
+    if !off.center
+        return vec - ones(SVector{D,Int})
+    else
+        return vec - SVector(ntuple(i -> cld(size(geo, i), 2), Val(D)))
     end
-    D = length(args) + 1
-    return LadderBoundaries{D,subgeometry{D-1}}(subgeometry(args...))
 end
 
-Base.size(geom::LadderBoundaries) = (2, size(geom.subgeometry)...)
+"""
+    neighbor_site(geom::CubicGrid, site, i)
 
-num_neighbours(geom::LadderBoundaries) = num_neighbours(geom.subgeometry) + 1
+Find the `i`-th neighbor of `site` in the geometry. If the move is illegal, return 0.
+"""
+function neighbor_site(g::CubicGrid{D}, mode, chosen) where {D}
+    # TODO: reintroduce LadderBoundaries small dimensions
+    return g[g[mode] + Directions(D)[chosen]]
+end
 
-function neighbour_site(geom::LadderBoundaries, site, i)
-    if i == 1
-        # Make odd site even, or make even site odd.
-        return site + site % 2 - (site + 1) % 2
-    else
-        i -= 1
-        subneighbour = neighbour_site(geom.subgeometry, cld(site, 2), i)
-        return ifelse(iszero(subneighbour), 0, 2 * subneighbour - site % 2)
-    end
+function BitStringAddresses.onr(address, geom::CubicGrid{<:Any,S}) where {S}
+    return SArray{Tuple{S...}}(onr(address))
+end
+function BitStringAddresses.onr(address::CompositeFS, geom::CubicGrid)
+    return map(fs -> onr(fs, geom), address.components)
 end
diff --git a/test/Hamiltonians.jl b/test/Hamiltonians.jl
index 7bd356035..2c0c8fa77 100644
--- a/test/Hamiltonians.jl
+++ b/test/Hamiltonians.jl
@@ -5,6 +5,7 @@ using Rimu
 using Test
 using DataFrames
 using Suppressor
+using StaticArrays
 
 function exact_energy(ham)
     dv = DVec(starting_address(ham) => 1.0)
@@ -149,6 +150,60 @@ using Rimu.Hamiltonians: momentum_transfer_excitation
     end
 end
 
+using Rimu.Hamiltonians: Directions, Displacements
+
+@testset "CubicGrid" begin
+    @testset "construtors and basic properties" begin
+        @test PeriodicBoundaries(3, 3) == CubicGrid(3, 3)
+        @test HardwallBoundaries(3, 4, 5) == CubicGrid((3, 4, 5), (false, false, false))
+        @test LadderBoundaries(2, 3, 4) == CubicGrid((2, 3, 4), (false, true, true))
+
+        for (dims, fold) in (
+            ((4,), (false,)), ((2, 5), (true, false)), ((5, 6, 7), (true, true, false))
+        )
+            geom = CubicGrid(dims, fold)
+            @test size(geom) == dims
+            @test length(geom) == prod(dims)
+            @test Rimu.Hamiltonians.fold(geom) == fold
+            @test eval(Meta.parse(repr(geom))) == geom
+        end
+    end
+
+    @testset "getindex" begin
+        g = CubicGrid((2,3,4), (false,true,false))
+        for i in 1:length(g)
+            v = SVector(Tuple(CartesianIndices((2,3,4))[i])...)
+            @test g[i] == v
+            @test g[v] == i
+        end
+
+        @test g[SVector(3, 1, 1)] == 0
+        @test g[SVector(1,4,1)] == 1
+        @test g[SVector(2,0,4)] == 24
+        @test g[SVector(2,3,0)] == 0
+    end
+
+    @testset "Directions" begin
+        @test Directions(1) == [[1], [-1]]
+        @test Directions(CubicGrid(2,3)) == [[1,0], [0,1], [-1,0], [0,-1]]
+        @test Directions(3) == [[1,0,0], [0,1,0], [0,0,1], [-1,0,0], [0,-1,0], [0,0,-1]]
+        @test_throws BoundsError Directions(3)[0]
+        @test_throws BoundsError Directions(2)[5]
+        @test_throws BoundsError Directions(1)[15]
+    end
+
+    @testset "Displacements" begin
+        @test collect(Displacements(CubicGrid(3), center=false)) == [[0], [1], [2]]
+        @test collect(Displacements(CubicGrid(3), center=true)) == [[-1], [0], [1]]
+
+        @test collect(Displacements(CubicGrid(2,2))) == [[0,0], [1,0], [0,1], [1,1]]
+        @test collect(Displacements(CubicGrid(2,3))) == [[0,0], [1,0], [0,1], [1,1], [0,2], [1,2]]
+        @test collect(Displacements(CubicGrid(2,3); center=true)) == [
+            [0,-1], [1,-1], [0,0], [1,0], [0,1], [1,1]
+        ]
+    end
+end
+
 @testset "Interface tests" begin
     for H in (
         HubbardReal1D(BoseFS((1, 2, 3, 4)); u=1.0, t=2.0),
@@ -340,19 +395,6 @@ end
         od_nonzeros = filter(!iszero, od_values)
         @test length(od_values) == 12
         @test length(od_nonzeros) == 4
-
-        H = HubbardRealSpace(f, geometry=LadderBoundaries(2, 6))
-        od_values = last.(offdiagonals(H, f))
-        od_nonzeros = filter(!iszero, od_values)
-        @test length(od_values) == 9
-        @test length(od_nonzeros) == 5
-
-        hard_ladder = LadderBoundaries(2, 6, subgeometry=HardwallBoundaries)
-        H = HubbardRealSpace(f, geometry=hard_ladder)
-        od_values = last.(offdiagonals(H, f))
-        od_nonzeros = filter(!iszero, od_values)
-        @test length(od_values) == 9
-        @test length(od_nonzeros) == 3
     end
     @testset "1D Bosons (single)" begin
         H1 = HubbardReal1D(BoseFS((1, 1, 1, 1, 1, 0)); u=2, t=3)
@@ -512,24 +554,19 @@ end
             )
             @test exact_energy(H3) < -16
         end
-        @testset "hardwall and ladder" begin
-            geom1 = LadderBoundaries(2, 3, subgeometry=HardwallBoundaries)
-            geom2 = HardwallBoundaries(2, 3)
-            geom3 = HardwallBoundaries(3, 2)
+        @testset "Hardwall" begin
+            geom1 = HardwallBoundaries(2, 3)
+            geom2 = HardwallBoundaries(3, 2)
             bose = BoseFS((1, 1, 1, 0, 0, 0))
             fermi = FermiFS((1, 0, 0, 0, 1, 0))
 
             H1 = HubbardRealSpace(bose, geometry=geom1)
             H2 = HubbardRealSpace(bose, geometry=geom2)
-            H3 = HubbardRealSpace(bose, geometry=geom3)
-            @test exact_energy(H1) == exact_energy(H2)
-            @test exact_energy(H1) ≈ exact_energy(H3)
+            @test exact_energy(H1) ≈ exact_energy(H2)
 
             H1 = HubbardRealSpace(fermi, geometry=geom1)
             H2 = HubbardRealSpace(fermi, geometry=geom2)
-            H3 = HubbardRealSpace(fermi, geometry=geom3)
-            @test exact_energy(H1) == exact_energy(H2)
-            @test exact_energy(H1) ≈ exact_energy(H3)
+            @test exact_energy(H1) ≈ exact_energy(H2)
         end
     end
 end
@@ -775,205 +812,304 @@ end
     @test f.shift ≈ a.shift
 end
 
-@testset "G2MomCorrelator" begin
-    # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1)
-    bfs1=BoseFS([0,2,0])
-    bfs2=BoseFS([0,1,0])
-    aIni = BoseFS2C(bfs1,bfs2)
-    v0 = DVec(
-        BoseFS2C((0, 2, 0), (0, 1, 0)) => 0.9999389545691221,
-        BoseFS2C((1, 1, 0), (0, 0, 1)) => -0.007812695959057453,
-        BoseFS2C((0, 1, 1), (1, 0, 0)) => -0.007812695959057453,
-        BoseFS2C((2, 0, 0), (1, 0, 0)) => 4.046694762039993e-5,
-        BoseFS2C((0, 0, 2), (0, 0, 1)) => 4.046694762039993e-5,
-        BoseFS2C((1, 0, 1), (0, 1, 0)) => 8.616127793651117e-5,
-    )
-    g0 = G2MomCorrelator(0)
-    g1 = G2MomCorrelator(1)
-    g2 = G2MomCorrelator(2)
-    g3 = G2MomCorrelator(3)
-    @test imag(dot(v0,g0,v0)) == 0 # should be strictly real
-    @test abs(imag(dot(v0,g3,v0))) < 1e-10
-    @test dot(v0,g0,v0) ≈ 0.65 rtol=0.01
-    @test dot(v0,g1,v0) ≈ 0.67 rtol=0.01
-    @test dot(v0,g2,v0) ≈ 0.67 rtol=0.01
-    @test dot(v0,g3,v0) ≈ 0.65 rtol=0.01
-    @test num_offdiagonals(g0,aIni) == 2
-
-    # on first component
-    g0f = G2MomCorrelator(0,:first)
-    g1f = G2MomCorrelator(1,:first)
-    @test imag(dot(v0,g0f,v0)) == 0 # should be strictly real
-    @test dot(v0,g0f,v0) ≈ 1.33 rtol=0.01
-    @test dot(v0,g1f,v0) ≈ 1.33 + 7.08e-5im rtol=0.01
-    # on second component
-    g0s = G2MomCorrelator(0,:second)
-    g1s = G2MomCorrelator(1,:second)
-    #@test_throws ErrorException("invalid ONR") get_offdiagonal(g0s,aIni,1) # should fail due to invalid ONR
-    @test dot(v0,g0s,v0) ≈ 1/3
-    @test dot(v0,g1s,v0) ≈ 1/3
-    # test against BoseFS
-    ham1 = HubbardMom1D(bfs1)
-    ham2 = HubbardMom1D(bfs2)
-    @test num_offdiagonals(g0f,aIni) == num_offdiagonals(ham1,bfs1)
-    @test num_offdiagonals(g0s,aIni) == num_offdiagonals(ham2,bfs2)
-    aIni = BoseFS2C(bfs2,bfs1) # flip bfs1 and bfs2
-    @test get_offdiagonal(g0s,aIni,1) == (BoseFS2C(BoseFS{1,3}((0, 1, 0)),BoseFS{2,3}((1, 0, 1))), 0.47140452079103173)
-    # test on BoseFS
-    @test diagonal_element(g0s,bfs1) == 4/3
-    @test diagonal_element(g0s,bfs2) == 1/3
-end
+using Rimu.Hamiltonians: circshift_dot
 
-@testset "G2RealCorrelator" begin
-    m = 6
-    n1 = 4
-    n2 = m
-    add1 = BoseFS((n1,0,0,0,0,0))
-    add2 = near_uniform(BoseFS{n2,m})
-
-    # localised state
-    @test diagonal_element(G2RealCorrelator(0), add1) == n1 * (n1 - 1) / m
-    @test diagonal_element(G2RealCorrelator(1), add1) == 0.
-
-    # constant density state
-    @test diagonal_element(G2RealCorrelator(0), add2) == (n2/m) * ((n2/m) - 1)
-    @test diagonal_element(G2RealCorrelator(1), add2) == (n2/m)^2
-
-    # local-local
-    comp = CompositeFS(add1,add1)
-    @test diagonal_element(G2RealCorrelator(0), comp) == 2n1 * (2n1 - 1) / m
-    @test diagonal_element(G2RealCorrelator(1), comp) == 0.
-
-    # local-uniform (assuming unit filling)
-    comp = CompositeFS(add1,add2)
-    @test diagonal_element(G2RealCorrelator(0), comp) == (n1 + 1) * n1 / m
-    @test diagonal_element(G2RealCorrelator(1), comp) == (2 * (n1 + 1) + (m - 2)) / m
-
-    # uniform-uniform
-    comp = CompositeFS(add2,add2)
-    @test diagonal_element(G2RealCorrelator(0), comp) == (2n2 / m) * (2 * (n2 / m) - 1)
-    @test diagonal_element(G2RealCorrelator(1), comp) == (2n2 / m)^2
-
-    # 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 "Correlation functions" begin
+    @testset "circhshift_dot" begin
+        for i in 1:10
+            A = rand(3, 4, 5)
+            B = rand(3, 4, 5)
+            inds = (rand(0:3), rand(0:4), rand(0:5))
+            @test circshift_dot(A, B, inds) ≈ dot(A, circshift(B, inds))
+        end
+    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 "G2RealCorrelator" begin
+        m = 6
+        n1 = 4
+        n2 = m
+        add1 = BoseFS((n1,0,0,0,0,0))
+        add2 = near_uniform(BoseFS{n2,m})
+
+        # localised state
+        @test diagonal_element(G2RealCorrelator(0), add1) == n1 * (n1 - 1) / m
+        @test diagonal_element(G2RealCorrelator(1), add1) == 0.0
+
+        # constant density state
+        @test diagonal_element(G2RealCorrelator(0), add2) == (n2/m) * ((n2/m) - 1)
+        @test diagonal_element(G2RealCorrelator(1), add2) == (n2/m)^2
+
+        # local-local
+        comp = CompositeFS(add1,add1)
+        @test diagonal_element(G2RealCorrelator(0), comp) == 2n1 * (2n1 - 1) / m
+        @test diagonal_element(G2RealCorrelator(1), comp) == 0.0
+
+        # local-uniform (assuming unit filling)
+        comp = CompositeFS(add1,add2)
+        @test diagonal_element(G2RealCorrelator(0), comp) == (n1 + 1) * n1 / m
+        @test diagonal_element(G2RealCorrelator(1), comp) == (2 * (n1 + 1) + (m - 2)) / m
+
+        # uniform-uniform
+        comp = CompositeFS(add2,add2)
+        @test diagonal_element(G2RealCorrelator(0), comp) == (2n2 / m) * (2 * (n2 / m) - 1)
+        @test diagonal_element(G2RealCorrelator(1), comp) == (2n2 / m)^2
+
+        # 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 "StringCorrelator" begin
-    m = 6
-    n1 = 4
-    n2 = m
+    @testset "G2RealSpace" begin
+        @testset "1D G2RealCorrelator comparison" begin
+            @testset "constructors" begin
+                g2_1 = G2RealSpace(CubicGrid(2, 2, 3), 1, 3)
+                g2_2 = G2RealSpace(CubicGrid(2, 2), 2)
+                g2_3 = G2RealSpace(CubicGrid(2, 2); sum_components=true)
+                @test g2_1 isa G2RealSpace{1,3}
+                @test g2_2 isa G2RealSpace{2,2}
+                @test g2_3 isa G2RealSpace{0,0}
+
+                @test eval(Meta.parse(repr(g2_1))) == g2_1
+                @test eval(Meta.parse(repr(g2_2))) == g2_2
+                @test eval(Meta.parse(repr(g2_3))) == g2_3
+
+                @test_throws ArgumentError G2RealSpace(CubicGrid(3), 1, 0)
+                @test_throws ArgumentError G2RealSpace(CubicGrid(2, 2), 0, 0)
+                @test_throws ArgumentError G2RealSpace(CubicGrid(1, 2, 3), -1, 2)
+                @test_throws ArgumentError G2RealSpace(CubicGrid(12), 3; sum_components=true)
+            end
+            @testset "single components" begin
+                addr = near_uniform(BoseFS{6,6})
+                H = HubbardReal1D(addr)
+                v = normalize!(H * (H * (H * DVec(addr => 1.0))))
+
+                g2 = dot(v, G2RealSpace(CubicGrid(6)), v)
+                for d in 0:5
+                    @test g2[d + 1] ≈ dot(v, G2RealCorrelator(d), v)
+                end
+            end
 
-    # 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)
+            @testset "sum of components" begin
+                addr = CompositeFS(BoseFS(4, 1=>1), BoseFS(4, 1=>1), BoseFS(4, 1=>1))
+                H = HubbardRealSpace(addr)
+                v = normalize!(H * (H * (H * DVec(addr => 1.0))))
 
-    localised_state = BoseFS((n2,0,0,0,0,0))
-    uniform_state = near_uniform(BoseFS{n2,m})
+                g2 = dot(v, G2RealSpace(CubicGrid(4); sum_components=true), v)
+                for d in 0:3
+                    @test g2[d + 1] ≈ dot(v, G2RealCorrelator(d), v)
+                end
 
-    S0 = StringCorrelator(0)
-    S1 = StringCorrelator(1)
-    S2 = StringCorrelator(2)
+                g2_nosum = dot(v, G2RealSpace(CubicGrid(4)), v)
+                @test iszero(g2_nosum)
+            end
+        end
 
-    @test num_offdiagonals(S0, localised_state) == 0
+        @testset "2-component, 2D" begin
+            c1 = BoseFS(1, 0, 0, 0, 0, 0)
+            c2 = BoseFS(1, 2, 3, 4, 5, 6)
+            addr = CompositeFS(c1, c2)
+            geom = CubicGrid((3, 2))
+
+            g2_11 = diagonal_element(G2RealSpace(geom, 1, 1), addr)
+            g2_12 = diagonal_element(G2RealSpace(geom, 1, 2), addr)
+            g2_21 = diagonal_element(G2RealSpace(geom, 2, 1), addr)
+            g2_22 = diagonal_element(G2RealSpace(geom, 2, 2), addr)
+
+            # Sum is N1 * (N2 - δ_12) / M
+            @test sum(g2_11) == 0
+            @test sum(g2_12) == 3.5
+            @test sum(g2_21) == 3.5
+            @test sum(g2_22) == 70
+
+            # Swapping components flips axes
+            @test g2_11 == [0 0; 0 0; 0 0]
+            @test g2_12 == [1 4; 2 5; 3 6] ./ 6
+            @test g2_21 == [1 4; 3 6; 2 5] ./ 6
+            @test g2_22[1, 1] == sum(n -> n * (n-1), onr(c2)) / 6
+            @test g2_22[3, 2] == dot(onr(c2, geom), circshift(onr(c2, geom), (2, 1))) / 6
+        end
 
-    # 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)
+        @testset "G2 is symmetric for translationally invariant ground states" begin
+            addr = near_uniform(BoseFS{3,18})
+            geom = CubicGrid((2,3,3), (false, true, true))
+            H = HubbardRealSpace(addr; geometry=geom)
+            bsr = BasisSetRep(H)
+            v0 = PDVec(zip(bsr.basis, eigen(Matrix(bsr)).vectors[:,1]))
 
-    # non unital near uniform state
-    @test @inferred diagonal_element(S0, non_unital_uniform_state) ≈ 2/9
+            g2 = dot(v0, G2RealSpace(geom), v0)
 
-    # 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
+            @test sum(g2) ≈ 3 * 2 / 18
+            @test g2[:,2,:] ≈ g2[:,3,:]
+            @test g2[:,:,2] ≈ g2[:,:,3]
+            @test minimum(g2) == first(g2)
+        end
+    end
 
-    # constant density uniform state
-    @test @inferred diagonal_element(S0, uniform_state) == 0
-    @test @inferred diagonal_element(S2, uniform_state) == 0
+    @testset "G2MomCorrelator" begin
+        # v0 is the exact ground state from BoseHubbardMom1D2C(aIni;ua=0,ub=0,v=0.1)
+        bfs1 = BoseFS([0, 2, 0])
+        bfs2 = BoseFS([0, 1, 0])
+        aIni = BoseFS2C(bfs1,bfs2)
+        v0 = DVec(
+            BoseFS2C((0, 2, 0), (0, 1, 0)) => 0.9999389545691221,
+            BoseFS2C((1, 1, 0), (0, 0, 1)) => -0.007812695959057453,
+            BoseFS2C((0, 1, 1), (1, 0, 0)) => -0.007812695959057453,
+            BoseFS2C((2, 0, 0), (1, 0, 0)) => 4.046694762039993e-5,
+            BoseFS2C((0, 0, 2), (0, 0, 1)) => 4.046694762039993e-5,
+            BoseFS2C((1, 0, 1), (0, 1, 0)) => 8.616127793651117e-5,
+        )
+        g0 = G2MomCorrelator(0)
+        g1 = G2MomCorrelator(1)
+        g2 = G2MomCorrelator(2)
+        g3 = G2MomCorrelator(3)
+        @test imag(dot(v0,g0,v0)) == 0 # should be strictly real
+        @test abs(imag(dot(v0,g3,v0))) < 1e-10
+        @test dot(v0,g0,v0) ≈ 0.65 rtol=0.01
+        @test dot(v0,g1,v0) ≈ 0.67 rtol=0.01
+        @test dot(v0,g2,v0) ≈ 0.67 rtol=0.01
+        @test dot(v0,g3,v0) ≈ 0.65 rtol=0.01
+        @test num_offdiagonals(g0,aIni) == 2
+
+        # on first component
+        g0f = G2MomCorrelator(0,:first)
+        g1f = G2MomCorrelator(1,:first)
+        @test imag(dot(v0,g0f,v0)) == 0 # should be strictly real
+        @test dot(v0,g0f,v0) ≈ 1.33 rtol=0.01
+        @test dot(v0,g1f,v0) ≈ 1.33 + 7.08e-5im rtol=0.01
+        # on second component
+        g0s = G2MomCorrelator(0,:second)
+        g1s = G2MomCorrelator(1,:second)
+        #@test_throws ErrorException("invalid ONR") get_offdiagonal(g0s,aIni,1) # should fail due to invalid ONR
+        @test dot(v0,g0s,v0) ≈ 1/3
+        @test dot(v0,g1s,v0) ≈ 1/3
+        # test against BoseFS
+        ham1 = HubbardMom1D(bfs1)
+        ham2 = HubbardMom1D(bfs2)
+        @test num_offdiagonals(g0f,aIni) == num_offdiagonals(ham1,bfs1)
+        @test num_offdiagonals(g0s,aIni) == num_offdiagonals(ham2,bfs2)
+        aIni = BoseFS2C(bfs2,bfs1) # flip bfs1 and bfs2
+        @test get_offdiagonal(g0s,aIni,1) == (BoseFS2C(BoseFS{1,3}((0, 1, 0)),BoseFS{2,3}((1, 0, 1))), 0.47140452079103173)
+        # test on BoseFS
+        @test diagonal_element(g0s,bfs1) == 4/3
+        @test diagonal_element(g0s,bfs2) == 1/3
+    end
 
-    # 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
+    @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
 
-    # Test show method
-    d = 5
-    output = @capture_out print(StringCorrelator(d))
-    @test output == "StringCorrelator($d)"
+    @testset "StringCorrelator" begin
+        m = 6
+        n1 = 4
+        n2 = m
 
-end
+        # 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)
 
-@testset "Momentum" begin
-    @test diagonal_element(Momentum(), BoseFS((0,0,2,1,3))) ≡ 2.0
-    @test diagonal_element(Momentum(fold=false), BoseFS((0,0,2,1,3))) ≡ 7.0
-    @test diagonal_element(Momentum(1), BoseFS((1,0,0,0))) ≡ -1.0
-    @test_throws MethodError diagonal_element(Momentum(2), BoseFS((0,1,0)))
+        localised_state = BoseFS((n2,0,0,0,0,0))
+        uniform_state = near_uniform(BoseFS{n2,m})
 
-    for add in (BoseFS2C((0,1,2,3,0), (1,2,3,4,5)), FermiFS2C((1,0,0,1), (0,0,1,0)))
-        @test diagonal_element(Momentum(1), add) + diagonal_element(Momentum(2), add) ≡
-            diagonal_element(Momentum(0), add)
-    end
+        S0 = StringCorrelator(0)
+        S1 = StringCorrelator(1)
+        S2 = StringCorrelator(2)
 
-    @test num_offdiagonals(Momentum(), BoseFS((0,1,0))) == 0
-    @test LOStructure(Momentum(2; fold=true)) == IsDiagonal()
-    @test Momentum(1)' === Momentum(1)
-end
+        @test num_offdiagonals(S0, localised_state) == 0
 
-@testset "DensityMatrixDiagonal" begin
-    @test diagonal_element(DensityMatrixDiagonal(5), FermiFS((0,1,0,1,0,1,0))) == 0
-    @test diagonal_element(DensityMatrixDiagonal(2; component=1), BoseFS((1,5,1,0))) == 5
+        # 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)
 
-    for add in (
-        CompositeFS(BoseFS((1,2,3,4,5)), BoseFS((5,4,3,2,1))),
-        BoseFS2C((1,2,3,4,5), (5,4,3,2,1))
-    )
-        for i in 1:5
-            @test diagonal_element(DensityMatrixDiagonal(i, component=1), add) == i
-            @test diagonal_element(DensityMatrixDiagonal(i, component=2), add) == 6 - i
-            @test diagonal_element(DensityMatrixDiagonal(i), add) == 6
+        # 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
+        @test diagonal_element(Momentum(), BoseFS((0,0,2,1,3))) ≡ 2.0
+        @test diagonal_element(Momentum(fold=false), BoseFS((0,0,2,1,3))) ≡ 7.0
+        @test diagonal_element(Momentum(1), BoseFS((1,0,0,0))) ≡ -1.0
+        @test_throws MethodError diagonal_element(Momentum(2), BoseFS((0,1,0)))
+
+        for add in (BoseFS2C((0,1,2,3,0), (1,2,3,4,5)), FermiFS2C((1,0,0,1), (0,0,1,0)))
+            @test diagonal_element(Momentum(1), add) + diagonal_element(Momentum(2), add) ≡
+                diagonal_element(Momentum(0), add)
         end
+
+        @test num_offdiagonals(Momentum(), BoseFS((0,1,0))) == 0
+        @test LOStructure(Momentum(2; fold=true)) == IsDiagonal()
+        @test Momentum(1)' === Momentum(1)
     end
 
-    @test num_offdiagonals(DensityMatrixDiagonal(1), BoseFS((0,1,0))) == 0
-    @test LOStructure(DensityMatrixDiagonal(2)) == IsDiagonal()
-    @test DensityMatrixDiagonal(15)' === DensityMatrixDiagonal(15)
+    @testset "DensityMatrixDiagonal" begin
+        @test diagonal_element(DensityMatrixDiagonal(5), FermiFS((0,1,0,1,0,1,0))) == 0
+        @test diagonal_element(DensityMatrixDiagonal(2; component=1), BoseFS((1,5,1,0))) == 5
+
+        for add in (
+            CompositeFS(BoseFS((1,2,3,4,5)), BoseFS((5,4,3,2,1))),
+            BoseFS2C((1,2,3,4,5), (5,4,3,2,1))
+            )
+            for i in 1:5
+                @test diagonal_element(DensityMatrixDiagonal(i, component=1), add) == i
+                @test diagonal_element(DensityMatrixDiagonal(i, component=2), add) == 6 - i
+                @test diagonal_element(DensityMatrixDiagonal(i), add) == 6
+            end
+        end
+
+        @test num_offdiagonals(DensityMatrixDiagonal(1), BoseFS((0,1,0))) == 0
+        @test LOStructure(DensityMatrixDiagonal(2)) == IsDiagonal()
+        @test DensityMatrixDiagonal(15)' === DensityMatrixDiagonal(15)
+    end
 end
 
 @testset "HubbardReal1DEP" begin
diff --git a/test/mpi_runtests.jl b/test/mpi_runtests.jl
index d445c33f7..526472739 100644
--- a/test/mpi_runtests.jl
+++ b/test/mpi_runtests.jl
@@ -266,6 +266,7 @@ end
             add = BoseFS((0,0,10,0,0))
             H = HubbardMom1D(add)
             D = DensityMatrixDiagonal(1)
+            G2 = G2RealSpace(PeriodicBoundaries(5))
 
             # Need to seed here to get the same random vectors on all ranks.
             Random.seed!(1)
@@ -300,6 +301,12 @@ end
                 @test norm(u, 2) ≈ norm(pu, 2)
                 @test norm(u, Inf) ≈ norm(pu, Inf)
             end
+            # dot only for G2
+            @test dot(v, G2, w) ≈ dot(pv, G2, pw)
+            @test dot(w, G2, v) ≈ dot(pw, G2, pv)
+
+            @test dot(v, G2, w) ≈ dot(pv, G2, pw, wm)
+            @test dot(w, G2, v) ≈ dot(pw, G2, pv, wm)
 
             @test dot(pv, (H, D), pw, wm) == (dot(pv, H, pw), dot(pv, D, pw))
             @test dot(pv, (H, D), pw) == (dot(pv, H, pw), dot(pv, D, pw))