Implement new geometries

src/Hamiltonians/Hamiltonians.jl

@@ -86,7 +86,7 @@ export G2RealCorrelator, G2RealSpace, SuperfluidCorrelator, DensityMatrixDiagona
 export SingleParticleExcitation, TwoParticleExcitation, ReducedDensityMatrix
 export StringCorrelator, G2MomCorrelator
 
-export CubicGrid, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries
+export CubicGrid, PeriodicBoundaries, HardwallBoundaries, LadderBoundaries, HoneycombLattice, HexagonalLattice
 
 export HOCartesianContactInteractions, HOCartesianEnergyConservedPerDim, HOCartesianCentralImpurity
 export AxialAngularMomentumHO

src/Hamiltonians/HubbardRealSpace.jl

@@ -140,11 +140,14 @@ is produced if `address`is incompatible with the interaction parameters `u`.
 
 ## Geometries
 
-Implemented [`CubicGrid`](@ref)s for keyword `geometry`
+Implemented [`Geometry`](@ref)s for keyword `geometry`
 
+* [`CubicGrid`](@ref)
 * [`PeriodicBoundaries`](@ref)
 * [`HardwallBoundaries`](@ref)
 * [`LadderBoundaries`](@ref)
+* [`HoneycombLattice`](@ref)
+* [`HexagonalLattice`](@ref)
 
 Default is `geometry=PeriodicBoundaries(M,)`, i.e. a one-dimensional lattice with the
 number of sites `M` inferred from the number of modes in `address`.
@@ -163,7 +166,7 @@ number of sites `M` inferred from the number of modes in `address`.
 struct HubbardRealSpace{
     C, # components
     A<:AbstractFockAddress,
-    G<:CubicGrid,
+    G<:Geometry,
     D, # dimension
     # The following need to be type params.
     T<:SVector{C,Float64},
@@ -181,7 +184,7 @@ end
 
 function HubbardRealSpace(
     address::AbstractFockAddress;
-    geometry::CubicGrid=PeriodicBoundaries((num_modes(address),)),
+    geometry::Geometry=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))
@@ -314,7 +317,7 @@ struct HubbardRealSpaceCompOffdiagonals{G,A} <: AbstractOffdiagonals{A,Float64}
 end
 
 function offdiagonals(h::HubbardRealSpace, comp, add)
-    neighbours = 2 * num_dimensions(h.geometry)
+    neighbours = num_neighbors(h.geometry)
     return HubbardRealSpaceCompOffdiagonals(
         h.geometry, add, h.t[comp], num_occupied_modes(add) * neighbours
     )
@@ -323,7 +326,7 @@ end
 Base.size(o::HubbardRealSpaceCompOffdiagonals) = (o.length,)
 
 @inline function Base.getindex(o::HubbardRealSpaceCompOffdiagonals, chosen)
-    neighbours = 2 * num_dimensions(o.geometry)
+    neighbours = num_neighbors(o.geometry)
     particle, neigh = fldmod1(chosen, neighbours)
     src_index = find_occupied_mode(o.address, particle)
     neigh = neighbor_site(o.geometry, src_index.mode, neigh)

src/Hamiltonians/geometry.jl

@@ -1,16 +1,20 @@
 """
-    CubicGrid(dims::NTuple{D,Int}, fold::NTuple{D,Bool})
+    abstract type Geometry{D}
 
-Represents a `D`-dimensional grid. Used to define a cubic lattice and boundary conditions
-for some [`AbstractHamiltonian`](@ref)s, e.g. with the keyword argument `geometry` when
-constructing a [`HubbardRealSpace`](@ref). 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).
+## Interface:
 
-* `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).
+* `Base.size(::Geometry)`: return the size of the lattice.
+* [`periodic_dimensions`](@ref)`(::Geometry)`: return a `D`-tuple of `Bool`s that signal which dimensions
+  of the geometry are periodic.
+* [`num_neighbors`](@ref)`(::Geometry)`: return the number of sites each site is connected
+  to.
+* [`neighbor_site`](@ref)`(::Geometry, site, chosen)`: pick a neighbor of site `site`.
+
+See also [`CubicGrid`](@ref), [`HoneycombLattice`](@ref), and [`HexagonalLattice`](@ref) for
+concrete implementations and [`HubbardRealSpace`](@ref) for a Hamiltonian that uses
+[`Geometry`](@ref).
+
+## Example
 
 ```julia
 julia> geo = CubicGrid((2,3), (true,false))
@@ -37,12 +41,104 @@ julia> geo[(3,3)]
 julia> geo[(3,4)] # returns 0 if out of bounds
 0
 ```
+"""
+abstract type Geometry{D} end
+
+Base.size(g::Geometry, i) = size(g)[i]
+Base.length(g::Geometry) = prod(size(g))
+
+"""
+    num_neighbors(::Geometry)
+
+Return the number of neighbors each site has in the geometry.
+"""
+num_neighbors
+
+"""
+    neighbor_site(::Geometry, site, i)
+
+Find the `i`-th neighbor of `site` in the geometry. If the move is illegal, return 0.
+"""
+neighbor_site
+
+"""
+    periodic_dimensions(::Geometry{D})
+
+Return a `D`-tuple of `Bool`s signaling which dimensions of the [`Geometry`](@ref) are
+periodic.
+"""
+periodic_dimensions
+
+"""
+    num_dimensions(geom::Geometry)
+
+Return the number of dimensions of the lattice in this geometry.
+"""
+num_dimensions(::Geometry{D}) where {D} = D
+
+function Base.getindex(g::Geometry{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::Geometry, i::Int) = SVector(Tuple(CartesianIndices(size(g))[i]))
+
+"""
+    fold_vec(g::Geometry{D}, vec::SVector{D,Int}) -> SVector{D,Int}
+
+Use the [`Geometry`](@ref) to fold the `vec` in each dimension. If folding is disabled in a
+dimension, the vector is allowed to go out of bounds.
+
+```julia
+julia> geo = CubicGrid((2,3), (true,false))
+CubicGrid{2}((2, 3), (true, false))
+
+julia> fold_vec(geo, (3,1))
+(1, 1)
+
+julia> fold_vec(geo, (3,4))
+(1, 4)
+```
+"""
+function fold_vec(g::Geometry{D}, vec::SVector{D,Int}) where {D}
+    (_fold_vec(Tuple(vec), periodic_dimensions(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
+
+"""
+    CubicGrid(dims::NTuple{D,Int}, fold::NTuple{D,Bool})
+
+Represents a `D`-dimensional grid. Used to define a cubic lattice and boundary conditions
+for some [`AbstractHamiltonian`](@ref)s, e.g. with the keyword argument `geometry` when
+constructing a [`HubbardRealSpace`](@ref). 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).
+
+* `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).
 
 See also [`PeriodicBoundaries`](@ref), [`HardwallBoundaries`](@ref) and
-[`LadderBoundaries`](@ref) for special-case constructors. See also
-[`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref).
+[`LadderBoundaries`](@ref) for special-case constructors and [`HoneycombLattice`](@ref), and
+[`HexagonalLattice`](@ref) for alternative lattice geometries.
+
+See also [`HubbardRealSpace`](@ref) and [`G2RealSpace`](@ref).
+
+A 3×2 `CubicGrid` is indexed as follows.
+```
+ │ │
+─1─4─
+ │ │
+─2─5─
+ │ │
+─3─6─
+ │ │
+```
 """
-struct CubicGrid{D,Dims,Fold}
+struct CubicGrid{D,Dims,Fold} <: Geometry{D}
     function CubicGrid(
         dims::NTuple{D,Int}, fold::NTuple{D,Bool}=ntuple(Returns(true), Val(D))
     ) where {D}
@@ -94,49 +190,7 @@ function Base.show(io::IO, g::CubicGrid{<:Any,Dims,Fold}) where {Dims,Fold}
 end
 
 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
-
-"""
-    num_dimensions(geom::LatticeCubicGrid)
-
-Return the number of dimensions of the lattice in this geometry.
-
-See also [`CubicGrid`](@ref).
-"""
-num_dimensions(::CubicGrid{D}) where {D} = D
-
-"""
-    fold_vec(g::CubicGrid{D}, vec::SVector{D,Int}) -> SVector{D,Int}
-
-Use the [`CubicGrid`](@ref) 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
-julia> geo = CubicGrid((2,3), (true,false))
-CubicGrid{2}((2, 3), (true, false))
-
-julia> fold_vec(geo, (3,1))
-(1, 1)
-
-julia> fold_vec(geo, (3,4))
-(1, 4)
-```
-"""
-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
-
-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]))
+periodic_dimensions(g::CubicGrid{<:Any,<:Any,Fold}) where {Fold} = Fold
 
 """
     Directions(D) <: AbstractVector{SVector{D,Int}}
@@ -161,7 +215,7 @@ See also [`CubicGrid`](@ref).
 struct Directions{D} <: AbstractVector{SVector{D,Int}} end
 
 Directions(D) = Directions{D}()
-Directions(::CubicGrid{D}) where {D} = Directions{D}()
+Directions(::Geometry{D}) where {D} = Directions{D}()
 
 Base.size(::Directions{D}) where {D} = (2D,)
 
@@ -222,21 +276,108 @@ Base.size(off::Displacements) = (length(off.geometry),)
     end
 end
 
-"""
-    neighbor_site(geom::CubicGrid, site, i)
-
-Find the `i`-th neighbor of `site` in the geometry. If the move is illegal, return 0.
-
-See also [`CubicGrid`](@ref).
-"""
 function neighbor_site(g::CubicGrid{D}, mode, chosen) where {D}
     # TODO: reintroduce LadderBoundaries small dimensions
     return g[g[mode] + Directions(D)[chosen]]
 end
 
+num_neighbors(::CubicGrid{D}) where {D} = 2D
+
 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
+
+"""
+    HoneycombLattice((height, width), fold=(true, true))
+
+A honeycomb lattice where each site has three neighbors. If periodic, each dimension of the
+lattice must be divisible by 2.
+
+A 4×4 `HoneycombLattice` is indexed as follows.
+```
+  ╲    ╱    ╲    ╱
+   1──5      9─13
+  ╱    ╲    ╱    ╲
+─2      6─10     14─
+  ╲    ╱    ╲    ╱
+   3──7     11─15
+  ╱    ╲    ╱    ╲
+─4      8─12     16─
+  ╲    ╱    ╲    ╱
+```
+"""
+struct HoneycombLattice{Dims,Fold} <: Geometry{2}
+    function HoneycombLattice(dims::Tuple{Int,Int}, fold=(true, true))
+        if fold[1] && isodd(dims[1])
+            throw(ArgumentError("if `fold[1] == true`, the lattice height must be even"))
+        end
+        if fold[2] && isodd(dims[2])
+            throw(ArgumentError("if `fold[2] == true`, the lattice width must be even"))
+        end
+        if any(<(1), dims)
+            throw(ArgumentError("lattice dimensions must be positive!"))
+        end
+        return new{dims,fold}()
+    end
+end
+HoneycombLattice(h, w, fold=(true, true)) = HoneycombLattice((h, w), fold)
+
+Base.size(::HoneycombLattice{Dims}) where {Dims} = Dims
+periodic_dimensions(::HoneycombLattice{<:Any,Fold}) where {Fold} = Fold
+
+function neighbor_site(geom::HoneycombLattice, mode, chosen)
+    i, j = geom[mode]
+    if chosen ≤ 2
+        target = SVector(i + ifelse(chosen == 1, 1, -1), j)
+    else
+        target = SVector(i, j + ifelse(isodd(i + j), -1, 1))
+    end
+    source = geom[mode]
+    return geom[target]
+end
+
+num_neighbors(::HoneycombLattice) = 3
+
+"""
+    HexagonalLattice((height, width), fold=(true, true))
+
+A hexagonal lattice where each site has 6 neighbors.
+
+A 3×2 `HexagonalLattice` is indexed as follows.
+```
+╲ │ ╲ │ ╲
+─ 1 ─ 4 ─
+╲ │ ╲ │ ╲
+─ 2 ─ 5 ─
+╲ │ ╲ │ ╲
+─ 3 ─ 6 ─
+╲ │ ╲ │ ╲
+```
+"""
+struct HexagonalLattice{Dims,Fold} <: Geometry{2}
+    function HexagonalLattice(dims::Tuple{Int,Int}, fold::Tuple{Bool,Bool}=(true,true))
+        if any(<(1), dims)
+            throw(ArgumentError("lattice dimensions must be positive!"))
+        end
+        return new{dims,fold}()
+    end
+end
+HexagonalLattice(h, w, fold=(true, true)) = HexagonalLattice((h, w), fold)
+
+Base.size(::HexagonalLattice{Dims}) where {Dims} = Dims
+periodic_dimensions(::HexagonalLattice{<:Any,Fold}) where {Fold} = Fold
+
+function neighbor_site(geom::HexagonalLattice, mode, chosen)
+    if chosen ≤ 4
+        # same as CubicGrid{2}
+        offset = Directions(2)[chosen]
+    else
+        offset = SVector(1, 1) * ifelse(chosen == 5, 1, -1)
+    end
+    return geom[geom[mode] + offset]
+end
+
+num_neighbors(::HexagonalLattice) = 6