Rework Geometry, add new G2 correlator (#254)

* New Geometry, new G2 correlator

* fix exports

* Move size and fold into type params

* cleanup

* Fix tests, G2RealCorrelator

* use SVectors in Geometry

* add some tests, doc updates

* Geometry tests

* summing over components

* fix doctest

* update docs

* improve test coverage

* Fix tests and docstrings

* is_zero -> iszero

* include *Boundaries in docs

* Small docstring tweaks

* clean up tests

* Apply suggestions from code review

Co-authored-by: Joachim Brand <j.brand@massey.ac.nz>

* Apply suggestions

* Rename Geometry and UnitVectors

* relax test to isapprox

* fix Direction tests

* refix

* fix interaction between MPI and array operators

* Use ref instead of vector in Allreduce

* Rename Offsets to Displacements

* bugfixes and documentation updates

---------

Co-authored-by: Joachim Brand <j.brand@massey.ac.nz>

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

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)

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

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)

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

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
 

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