1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
use super::offsets::hypercubic;
/// What a forward bond that crosses the edge of the cell along one dimension becomes.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Boundary {
Periodic,
/// The bond is removed.
Open,
/// The bond changes sign each time it crosses.
Antiperiodic,
}
impl TryFrom<&str> for Boundary {
type Error = String;
fn try_from(s: &str) -> Result<Self, Self::Error> {
match s {
"periodic" => Ok(Self::Periodic),
"open" => Ok(Self::Open),
"antiperiodic" => Ok(Self::Antiperiodic),
_ => Err(format!(
"unknown boundary '{s}', expected 'periodic', 'open' or 'antiperiodic'"
)),
}
}
}
/// Periodic Bravais lattice with precomputed neighbor table.
///
/// Sites are indexed in row-major (C) order. Couplings are stored in a flat
/// array of length `n_spins * n_neighbors`, where element `i * n_neighbors + d`
/// is the coupling on the bond from site `i` to its forward neighbor in
/// direction `d`.
pub struct Lattice {
/// Extent along each dimension (e.g. `[8, 8, 8]`).
pub shape: Vec<usize>,
/// Row-major strides: `strides[d] = product of shape[d+1..]`.
pub strides: Vec<usize>,
/// Total number of sites (`shape.iter().product()`).
pub n_spins: usize,
/// Number of spatial dimensions (`shape.len()`).
pub n_dims: usize,
/// Number of forward neighbor directions per site.
pub n_neighbors: usize,
/// Forward neighbors, length `n_spins * n_neighbors`.
/// Layout: `fwd_neighbors[i * n_neighbors + d]`.
fwd_neighbors: Vec<u32>,
/// Backward neighbors, length `n_spins * n_neighbors`.
/// Layout: `bwd_neighbors[i * n_neighbors + d]`.
bwd_neighbors: Vec<u32>,
offsets: Vec<Vec<isize>>,
canonical_square_shape: Option<(usize, usize)>,
canonical_cubic_shape: Option<[usize; 3]>,
}
impl Lattice {
/// Create a hypercubic lattice with the given shape (e.g. `vec![16, 16]`).
pub fn new(shape: Vec<usize>) -> Self {
let n_dims = shape.len();
let canonical_square_shape = (n_dims == 2).then(|| (shape[0], shape[1]));
let canonical_cubic_shape = (n_dims == 3).then(|| [shape[0], shape[1], shape[2]]);
let mut lattice = Self::with_offsets(shape, hypercubic(n_dims));
lattice.canonical_square_shape = canonical_square_shape;
lattice.canonical_cubic_shape = canonical_cubic_shape;
lattice
}
/// Create a lattice with arbitrary neighbor offsets.
///
/// Each offset is a vector of length `n_dims` specifying a displacement in
/// lattice coordinates. The backward neighbor is the negation of the offset.
/// Periodic boundary conditions are applied via `rem_euclid`.
pub fn with_offsets(shape: Vec<usize>, offsets: Vec<Vec<isize>>) -> Self {
let n_dims = shape.len();
let n_neighbors = offsets.len();
let n_spins: usize = shape.iter().product();
for (idx, off) in offsets.iter().enumerate() {
assert_eq!(
off.len(),
n_dims,
"offset {idx} has length {}, expected {n_dims}",
off.len(),
);
}
let mut strides = vec![1usize; n_dims];
for d in (0..n_dims.saturating_sub(1)).rev() {
strides[d] = strides[d + 1] * shape[d + 1];
}
let mut fwd_neighbors = vec![0u32; n_spins * n_neighbors];
let mut bwd_neighbors = vec![0u32; n_spins * n_neighbors];
for i in 0..n_spins {
let coords: Vec<usize> = (0..n_dims).map(|d| (i / strides[d]) % shape[d]).collect();
for (d, off) in offsets.iter().enumerate() {
for (sign, table) in [(1isize, &mut fwd_neighbors), (-1isize, &mut bwd_neighbors)] {
let mut flat = 0usize;
for dim in 0..n_dims {
let c = (coords[dim] as isize + sign * off[dim])
.rem_euclid(shape[dim] as isize)
as usize;
flat += c * strides[dim];
}
table[i * n_neighbors + d] = flat as u32;
}
}
}
Self {
shape,
strides,
n_spins,
n_dims,
n_neighbors,
fwd_neighbors,
bwd_neighbors,
offsets,
canonical_square_shape: None,
canonical_cubic_shape: None,
}
}
#[inline]
pub(crate) fn square_shape(&self) -> Option<(usize, usize)> {
self.canonical_square_shape
}
/// Shape of a lattice built by [`Lattice::new`] with three dimensions, whose
/// directions have strides `L1 * L2`, `L2` and 1.
#[inline]
pub(crate) fn cubic_shape(&self) -> Option<[usize; 3]> {
self.canonical_cubic_shape
}
/// How many times the forward bond from `site` along offset `d` crosses the edge of
/// the cell along `dim`, with the sign of the crossing.
pub fn forward_wraps(&self, site: usize, d: usize, dim: usize) -> isize {
let extent = self.shape[dim] as isize;
let coord = ((site / self.strides[dim]) % self.shape[dim]) as isize;
(coord + self.offsets[d][dim]).div_euclid(extent)
}
/// Per forward bond, the factor that turns periodic couplings into ones with the
/// given `boundaries` (one per dimension): 1 for a bond kept as is, -1 for one
/// that changes sign, 0 for one removed. Laid out like the couplings.
///
/// The factors are in {-1, 0, 1}, so unit couplings stay unit couplings, and a
/// pair of runs with periodic and antiperiodic boundaries gives the domain-wall
/// free energy. Wrapping observables are zero with open boundaries.
pub fn boundary_factors(&self, boundaries: &[Boundary]) -> Result<Vec<i8>, String> {
if boundaries.len() != self.n_dims {
return Err(format!(
"expected one boundary per dimension ({}), got {}",
self.n_dims,
boundaries.len()
));
}
let mut factors = vec![1i8; self.n_spins * self.n_neighbors];
for site in 0..self.n_spins {
for d in 0..self.n_neighbors {
for (dim, &boundary) in boundaries.iter().enumerate() {
let wraps = self.forward_wraps(site, d, dim);
let factor = &mut factors[site * self.n_neighbors + d];
match boundary {
Boundary::Periodic => {}
Boundary::Open if wraps != 0 => *factor = 0,
Boundary::Open => {}
Boundary::Antiperiodic if wraps % 2 != 0 => *factor = -*factor,
Boundary::Antiperiodic => {}
}
}
}
}
Ok(factors)
}
/// Whether a neighbor offset maps every site to itself, as an extent of 1 along
/// it does. That bond is a constant of the energy, but the single-spin kernels
/// would read it as a field.
pub fn has_self_neighbors(&self) -> bool {
// Offsets act by translation, so site 0 decides for all sites.
self.n_spins > 0 && (0..self.n_neighbors).any(|d| self.neighbor_fwd(0, d) == 0)
}
#[inline]
pub fn neighbor_fwd(&self, flat_idx: usize, dim: usize) -> usize {
self.fwd_neighbors[flat_idx * self.n_neighbors + dim] as usize
}
#[inline]
pub fn neighbor_bwd(&self, flat_idx: usize, dim: usize) -> usize {
self.bwd_neighbors[flat_idx * self.n_neighbors + dim] as usize
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_2d_neighbors() {
// 3x4 lattice
let lat = Lattice::new(vec![3, 4]);
assert_eq!(lat.n_spins, 12);
assert_eq!(lat.strides, vec![4, 1]);
// Spin 0 = (0,0): forward in dim 0 -> (1,0)=4, forward in dim 1 -> (0,1)=1
assert_eq!(lat.neighbor_fwd(0, 0), 4);
assert_eq!(lat.neighbor_fwd(0, 1), 1);
// Spin 0 = (0,0): backward in dim 0 -> (2,0)=8 (wrap), backward in dim 1 -> (0,3)=3 (wrap)
assert_eq!(lat.neighbor_bwd(0, 0), 8);
assert_eq!(lat.neighbor_bwd(0, 1), 3);
// Spin 11 = (2,3): forward in dim 0 -> (0,3)=3 (wrap), forward in dim 1 -> (2,0)=8 (wrap)
assert_eq!(lat.neighbor_fwd(11, 0), 3);
assert_eq!(lat.neighbor_fwd(11, 1), 8);
}
#[test]
fn test_3d_neighbors() {
let lat = Lattice::new(vec![2, 3, 4]);
assert_eq!(lat.n_spins, 24);
assert_eq!(lat.strides, vec![12, 4, 1]);
// Spin 0 = (0,0,0)
assert_eq!(lat.neighbor_fwd(0, 0), 12); // (1,0,0)
assert_eq!(lat.neighbor_fwd(0, 1), 4); // (0,1,0)
assert_eq!(lat.neighbor_fwd(0, 2), 1); // (0,0,1)
}
#[test]
fn test_triangular_neighbors() {
use super::super::offsets::triangular;
// 4x4 triangular lattice: offsets [1,0], [0,1], [1,-1]
let lat = Lattice::with_offsets(vec![4, 4], triangular());
assert_eq!(lat.n_neighbors, 3);
assert_eq!(lat.n_spins, 16);
// Site 0 = (0,0)
// offset [1,0] -> (1,0) = 4
assert_eq!(lat.neighbor_fwd(0, 0), 4);
// offset [0,1] -> (0,1) = 1
assert_eq!(lat.neighbor_fwd(0, 1), 1);
// offset [1,-1] -> (1, -1 mod 4) = (1,3) = 7
assert_eq!(lat.neighbor_fwd(0, 2), 7);
// backward of [1,0] from (0,0) -> (-1 mod 4, 0) = (3,0) = 12
assert_eq!(lat.neighbor_bwd(0, 0), 12);
// backward of [0,1] from (0,0) -> (0, -1 mod 4) = (0,3) = 3
assert_eq!(lat.neighbor_bwd(0, 1), 3);
// backward of [1,-1] from (0,0) -> (-1 mod 4, 1 mod 4) = (3,1) = 13
assert_eq!(lat.neighbor_bwd(0, 2), 13);
// Site 5 = (1,1)
// offset [1,-1] -> (2, 0) = 8
assert_eq!(lat.neighbor_fwd(5, 2), 8);
// backward of [1,-1] from (1,1) -> (0, 2) = 2
assert_eq!(lat.neighbor_bwd(5, 2), 2);
// Site 15 = (3,3): all forward neighbors wrap
// offset [1,0] -> (0,3) = 3
assert_eq!(lat.neighbor_fwd(15, 0), 3);
// offset [0,1] -> (3,0) = 12
assert_eq!(lat.neighbor_fwd(15, 1), 12);
// offset [1,-1] -> (0, 2) = 2
assert_eq!(lat.neighbor_fwd(15, 2), 2);
}
#[test]
fn boundary_factors_mark_the_bonds_that_cross_the_cell_edge() {
// 3x4: dimension 0 has stride 4, dimension 1 stride 1.
let lattice = Lattice::new(vec![3, 4]);
let factors = lattice
.boundary_factors(&[Boundary::Antiperiodic, Boundary::Open])
.unwrap();
for site in 0..12 {
let (row, col) = (site / 4, site % 4);
assert_eq!(factors[site * 2], if row == 2 { -1 } else { 1 }, "{site}");
assert_eq!(
factors[site * 2 + 1],
if col == 3 { 0 } else { 1 },
"{site}"
);
}
// A diagonal offset crosses both edges from the corner, so both rules apply.
let tri = Lattice::with_offsets(vec![3, 3], vec![vec![1, 1]]);
let factors = tri
.boundary_factors(&[Boundary::Antiperiodic, Boundary::Antiperiodic])
.unwrap();
assert_eq!(factors[8], 1, "two sign changes cancel");
assert_eq!(factors[2], -1);
assert!(lattice.boundary_factors(&[Boundary::Open]).is_err());
}
#[test]
fn detects_offsets_that_map_sites_to_themselves() {
assert!(Lattice::new(vec![4, 1]).has_self_neighbors());
assert!(Lattice::with_offsets(vec![2, 3], vec![vec![2, 0]]).has_self_neighbors());
// Extent 2 gives a double bond to a distinct site, which the kernels handle.
assert!(!Lattice::new(vec![2, 2]).has_self_neighbors());
assert!(
!Lattice::with_offsets(vec![4, 4], vec![vec![1, 0], vec![1, -1]]).has_self_neighbors()
);
}
#[test]
fn square_shape_uses_construction_time_topology() {
let mut square = Lattice::new(vec![5, 7]);
square.shape = vec![1, 1];
square.n_dims = 3;
assert_eq!(square.square_shape(), Some((5, 7)));
let custom = Lattice::with_offsets(vec![5, 7], vec![vec![0, 1], vec![1, 0]]);
assert_eq!(custom.square_shape(), None);
}
}