GitHub

Momentum space

LatticeCore's momentum-space layer is designed to be usable without thinking about whether your lattice is periodic or quasiperiodic. A periodic lattice produces a PeriodicMomentumLattice; a quasicrystal produces a BraggPeakSet. Both subtype AbstractMomentumLattice, and both are themselves AbstractLattice{D, T} — so every observable that walks num_sites / position on a real-space lattice transparently walks num_k_points / k_point on a momentum lattice.

For the physics background, see Reciprocal lattice and Brillouin zone.

The trait-dispatched entry point

Any AbstractLattice declares its k-space capability through the reciprocal_support trait:

The unified helper is momentum_lattice(lat), which dispatches on the trait and throws ArgumentError for NoReciprocal:

lat = SimpleSquareLattice(4, 4, PeriodicAxis())
reciprocal_support(lat)      # HasReciprocal()
ml = momentum_lattice(lat)   # PeriodicMomentumLattice{2, Float64}

obc = SimpleSquareLattice(4, 4, OpenAxis())
reciprocal_support(obc)      # NoReciprocal()
momentum_lattice(obc)        # throws ArgumentError

For a periodic lattice, momentum_lattice returns the same object as reciprocal_lattice(lat).

Building a mesh from a basis

If you know the reciprocal basis and you want to sample the BZ directly, use monkhorst_pack or gamma_centered.

using StaticArrays
B = SMatrix{2, 2, Float64}(2π, 0.0,
                           0.0, 2π)

mp = monkhorst_pack(B, (4, 4))          # 16 Monkhorst–Pack k-points
mp isa PeriodicMomentumLattice{2, Float64}

γ = gamma_centered(B, (4, 4))           # 16 Γ-inclusive k-points

Both return a PeriodicMomentumLattice{D, T} that behaves like any other lattice: num_sites(mp) is 16, position(mp, 1) is a k-vector, and boundary(mp) is a uniform periodic wrap (k-space repeats modulo the reciprocal basis).

The difference between the two is where the mesh sits relative to Γ:

  • gamma_centered((3, 3)) produces fractions {0, 1/3, 2/3} on each axis — Γ is one of the sampled points.
  • monkhorst_pack((3, 3)) produces fractions {(n + 0.5)/3 - 0.5} = {-1/3, 0, 1/3} — the mesh is centred on Γ but offset by half a step, which is the convention solid-state codes use when integrating smooth functions over the BZ.

Reference-lattice reciprocals

LineLattice and SimpleSquareLattice both implement basis_vectors (the unit-spacing identity) and reciprocal_lattice, which constructs $B = 2\pi A^{-\top}$ and feeds it through monkhorst_pack. On a square lattice that gives you $B = 2\pi \, \mathrm{I}_2$ and a mesh whose density matches the real-space sample:

lat = SimpleSquareLattice(4, 4, PeriodicAxis())
A = basis_vectors(lat)                   # SMatrix{2, 2, Float64}(1, 0, 0, 1)
ml = reciprocal_lattice(lat)
B = reciprocal_basis(ml)                  # 2π * I

# Orthogonality: a_i · b_j = 2π δ_ij
using LinearAlgebra
[dot(A[:, i], B[:, j]) for i in 1:2, j in 1:2]
# → [2π 0; 0 2π]

Non-periodic boundaries (OBC or cylinders) throw ArgumentError from reciprocal_lattice, matching the reciprocal_support trait.

Structure factor

structure_factor evaluates

\[S(\mathbf{k}) = \frac{1}{N}\, \left| \sum_i s_i \, e^{-i \mathbf{k} \cdot \mathbf{r}_i} \right|^2\]

for a scalar configuration state. The default per-k-point implementation is naive O(N); when a whole momentum lattice is passed the call is dispatched on reciprocal_support(lat):

  • HasReciprocal() (Bravais periodic): the FFT extension LatticeCoreFFTWExt (loaded automatically by using FFTW) takes the O(N log N) path on regular meshes whose dims match the lattice grid, falling back to naive otherwise.
  • HasFourierModule() (cut-and-project quasicrystals): the LatticeCoreNFFTExt extension reserves a NUFFT entry point (issue #28); for now it returns the naive result.
  • otherwise: naive.

Opting a custom lattice into the FFT fast path

The grid-layout hooks the FFT extension dispatches on live on LatticeCore itself, so downstream packages can register their own Bravais lattices without depending on FFTW. Add two methods next to the lattice definition:

using LatticeCore

struct MyLattice <: AbstractLattice{2,Float64}
    Lx::Int
    Ly::Int
    # ...
end

# Tell LatticeCore that `state[i]` reshapes onto the natural
# `(Lx, Ly)` grid (e.g. via `site_index(x, y) = (y - 1) * Lx + x`):
LatticeCore._has_known_grid(::MyLattice) = true
LatticeCore._reshape_state(::MyLattice, state, dims) = reshape(state, dims)

With those two lines in place, structure_factor(lat, state, ml) will use the FFT path whenever ml::PeriodicMomentumLattice matches the lattice dims and using FFTW has been issued; if the user hasn't loaded FFTW the call still works — it just stays on the naive helper.

Three canonical checks:

Ferromagnet: S(0) = N

lat = SimpleSquareLattice(4, 4, PeriodicAxis())
state = ones(Int8, num_sites(lat))
structure_factor(lat, state, SVector(0.0, 0.0))   # 16.0

Ferromagnet vanishes at BZ-interior k

# k = (π/2, 0) is inside the BZ, so Σ_x exp(-i(π/2)x) = 0 for x=1..4
structure_factor(lat, state, SVector(π / 2, 0.0))  # ≈ 0

Note that S does not vanish at (2π, 0) — that is a reciprocal-lattice vector, hence Γ-equivalent, and the ferromagnetic signal concentrates there too.

Néel antiferromagnet: S(π, π) = N

lat = SimpleSquareLattice(4, 4, PeriodicAxis())
N = num_sites(lat)
neel = Int8[
    (-1)^(Int(position(lat, i)[1]) + Int(position(lat, i)[2]))
    for i in 1:N
]

structure_factor(lat, neel, SVector(π, π))    # 16.0
structure_factor(lat, neel, SVector(0.0, 0.0)) # ≈ 0

Bipartiteness is exactly the condition that this momentum lies on the reciprocal lattice of the sample; the 4 × 4 square is bipartite (both axes even), so the signal is clean.

Sweeping a whole momentum lattice

You can pass a momentum lattice in place of a single k-vector and LatticeCore will evaluate the structure factor at every k-point in the mesh:

ml = reciprocal_lattice(lat)
Sks = structure_factor(lat, neel, ml)        # Vector{Float64}, length num_sites(ml)
argmax(Sks)                                   # the k-index where the peak lives

The result is a plain Vector{Float64} of the same length as num_k_points(ml), indexed the same way as k_point(ml, i).

Quasicrystal Fourier modules (preview)

Everything above generalises: a BraggPeakSet subtypes AbstractMomentumLattice, exposes num_k_points(bps) = length(bps.peaks), and answers position(bps, i) = bps.peaks[i]. That means the same structure_factor(lat, state, ml) code path works for a quasicrystal — there is no "quasi" branch in the observer. You iterate over the Bragg peaks and each one returns a plain Float64.

The generation algorithm — choosing an acceptance window, walking the higher-dimensional reciprocal lattice, computing the window Fourier transform — lives in QuasiCrystal.jl, not here. LatticeCore ships the type skeletons HyperReciprocalLattice and AcceptanceWindow so downstream packages have a stable interface to target.

See also