hoomd_microstate/boundary/periodic/
cuboid.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement periodic boundary conditions for cuboids in cartesian space.
5
6use arrayvec::ArrayVec;
7
8use crate::{
9    boundary::{
10        Error, GenerateGhosts, MAX_GHOSTS, MaximumAllowableInteractionRange, Periodic, Wrap,
11    },
12    property::Position,
13};
14use hoomd_geometry::{IsPointInside, shape::Hypercuboid};
15use hoomd_utility::valid::PositiveReal;
16use hoomd_vector::Cartesian;
17
18impl<const N: usize> MaximumAllowableInteractionRange for Hypercuboid<N> {
19    /// The largest value that the maximum interaction range can take.
20    ///
21    /// For a cuboid, the maximum is
22    /// ```math
23    /// \frac{L_\mathrm{min}}{2}
24    /// ```
25    /// where $`L_\mathrm{min}`$ is the smallest edge length.
26    ///
27    /// # Example
28    ///
29    /// ```
30    /// use hoomd_geometry::shape::Hypercuboid;
31    /// use hoomd_microstate::boundary::MaximumAllowableInteractionRange;
32    ///
33    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
34    /// let rectangular_prism = Hypercuboid {
35    ///     edge_lengths: [2.0.try_into()?, 3.0.try_into()?, 9.0.try_into()?],
36    /// };
37    ///
38    /// assert_eq!(rectangular_prism.maximum_allowable_interaction_range(), 1.0);
39    /// # Ok(())
40    /// # }
41    /// ```
42    #[inline]
43    fn maximum_allowable_interaction_range(&self) -> f64 {
44        let minimum_l = self
45            .edge_lengths
46            .iter()
47            .map(PositiveReal::get)
48            .reduce(f64::min)
49            .expect("cuboid should have dimension 1 or greater");
50        minimum_l / 2.0
51    }
52}
53
54impl<const N: usize, P> Wrap<P> for Periodic<Hypercuboid<N>>
55where
56    P: Position<Position = Cartesian<N>>,
57{
58    /// Wrap any cartesian vector to the inside of the given cuboid.
59    ///
60    /// # Example
61    ///
62    /// ```
63    /// use hoomd_geometry::shape::Rectangle;
64    /// use hoomd_microstate::{
65    ///     boundary::{Periodic, Wrap},
66    ///     property::Point,
67    /// };
68    /// use hoomd_vector::Cartesian;
69    ///
70    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
71    /// let periodic =
72    ///     Periodic::new(2.5, Rectangle::with_equal_edges(10.0.try_into()?))?;
73    /// let point = Point::new(Cartesian::from([6.0, -15.0]));
74    ///
75    /// let wrapped_point = periodic.wrap(point)?;
76    /// assert_eq!(wrapped_point.position, [-4.0, -5.0].into());
77    /// # Ok(())
78    /// # }
79    /// ```
80    #[inline]
81    fn wrap(&self, properties: P) -> Result<P, Error> {
82        let mut properties = properties;
83        let r = properties.position_mut();
84
85        for (coordinate, edge_length) in r.coordinates.iter_mut().zip(self.shape.edge_lengths) {
86            let edge_length = edge_length.get();
87            let lambda = *coordinate / edge_length;
88            let lambda = lambda - lambda.round();
89            let lambda = if lambda == 0.5 { -0.5 } else { lambda };
90            *coordinate = lambda * edge_length;
91        }
92        Ok(properties)
93    }
94}
95
96impl<S> GenerateGhosts<S> for Periodic<Hypercuboid<2>>
97where
98    S: Position<Position = Cartesian<2>> + Copy + Default,
99{
100    #[inline]
101    fn maximum_interaction_range(&self) -> f64 {
102        self.maximum_interaction_range
103    }
104
105    /// Place periodic images of sites near the edge of the periodic boundary.
106    ///
107    /// For 2D cuboids, `generate_ghosts` places ghosts near the 4 edges and 4
108    /// vertices.
109    #[inline]
110    fn generate_ghosts(&self, site_properties: &S) -> ArrayVec<S, MAX_GHOSTS> {
111        let mut result = ArrayVec::new();
112
113        let r = site_properties.position();
114        let max = self.shape.maximal_extents();
115        let min = self.shape.minimal_extents();
116
117        if !self.shape.is_point_inside(r) {
118            return result;
119        }
120
121        let new_site = |x, y| {
122            let mut new_site = *site_properties;
123            new_site.position_mut()[0] += x * self.shape.edge_lengths[0].get();
124            new_site.position_mut()[1] += y * self.shape.edge_lengths[1].get();
125            new_site
126        };
127
128        let near_left = r[0] < min[0] + self.maximum_interaction_range;
129        let near_right = r[0] > max[0] - self.maximum_interaction_range;
130        let near_top = r[1] > max[1] - self.maximum_interaction_range;
131        let near_bottom = r[1] < min[1] + self.maximum_interaction_range;
132
133        if near_right {
134            result.push(new_site(-1.0, 0.0));
135        }
136        if near_left {
137            result.push(new_site(1.0, 0.0));
138        }
139        if near_top {
140            result.push(new_site(0.0, -1.0));
141        }
142        if near_bottom {
143            result.push(new_site(0.0, 1.0));
144        }
145        if near_right && near_top {
146            result.push(new_site(-1.0, -1.0));
147        }
148        if near_right && near_bottom {
149            result.push(new_site(-1.0, 1.0));
150        }
151        if near_left && near_top {
152            result.push(new_site(1.0, -1.0));
153        }
154        if near_left && near_bottom {
155            result.push(new_site(1.0, 1.0));
156        }
157
158        result
159    }
160}
161
162impl<S> GenerateGhosts<S> for Periodic<Hypercuboid<3>>
163where
164    S: Position<Position = Cartesian<3>> + Copy + Default,
165{
166    #[inline]
167    fn maximum_interaction_range(&self) -> f64 {
168        self.maximum_interaction_range
169    }
170
171    /// Place periodic images of sites near the edge of the periodic boundary.
172    ///
173    /// For 3D cuboids, `generate_ghosts` places ghosts near the 6 faces, 12 edges,
174    /// and 8 vertices.
175    #[inline]
176    fn generate_ghosts(&self, site_properties: &S) -> ArrayVec<S, MAX_GHOSTS> {
177        let mut result = ArrayVec::new();
178
179        let r = site_properties.position();
180        let max = self.shape.maximal_extents();
181        let min = self.shape.minimal_extents();
182
183        if !self.shape.is_point_inside(r) {
184            return result;
185        }
186
187        let new_site = |x, y, z| {
188            let mut new_site = *site_properties;
189            new_site.position_mut()[0] += x * self.shape.edge_lengths[0].get();
190            new_site.position_mut()[1] += y * self.shape.edge_lengths[1].get();
191            new_site.position_mut()[2] += z * self.shape.edge_lengths[2].get();
192            new_site
193        };
194
195        let near_left = r[0] < min[0] + self.maximum_interaction_range;
196        let near_right = r[0] > max[0] - self.maximum_interaction_range;
197        let near_top = r[1] > max[1] - self.maximum_interaction_range;
198        let near_bottom = r[1] < min[1] + self.maximum_interaction_range;
199        let near_front = r[2] > max[2] - self.maximum_interaction_range;
200        let near_back = r[2] < min[2] + self.maximum_interaction_range;
201
202        if near_right {
203            result.push(new_site(-1.0, 0.0, 0.0));
204        }
205        if near_left {
206            result.push(new_site(1.0, 0.0, 0.0));
207        }
208        if near_top {
209            result.push(new_site(0.0, -1.0, 0.0));
210        }
211        if near_bottom {
212            result.push(new_site(0.0, 1.0, 0.0));
213        }
214        if near_front {
215            result.push(new_site(0.0, 0.0, -1.0));
216        }
217        if near_back {
218            result.push(new_site(0.0, 0.0, 1.0));
219        }
220
221        if near_right && near_top {
222            result.push(new_site(-1.0, -1.0, 0.0));
223        }
224        if near_right && near_bottom {
225            result.push(new_site(-1.0, 1.0, 0.0));
226        }
227        if near_right && near_front {
228            result.push(new_site(-1.0, 0.0, -1.0));
229        }
230        if near_right && near_back {
231            result.push(new_site(-1.0, 0.0, 1.0));
232        }
233        if near_left && near_top {
234            result.push(new_site(1.0, -1.0, 0.0));
235        }
236        if near_left && near_bottom {
237            result.push(new_site(1.0, 1.0, 0.0));
238        }
239        if near_left && near_front {
240            result.push(new_site(1.0, 0.0, -1.0));
241        }
242        if near_left && near_back {
243            result.push(new_site(1.0, 0.0, 1.0));
244        }
245
246        if near_top && near_front {
247            result.push(new_site(0.0, -1.0, -1.0));
248        }
249        if near_bottom && near_front {
250            result.push(new_site(0.0, 1.0, -1.0));
251        }
252        if near_top && near_back {
253            result.push(new_site(0.0, -1.0, 1.0));
254        }
255        if near_bottom && near_back {
256            result.push(new_site(0.0, 1.0, 1.0));
257        }
258
259        if near_right && near_top && near_front {
260            result.push(new_site(-1.0, -1.0, -1.0));
261        }
262        if near_right && near_top && near_back {
263            result.push(new_site(-1.0, -1.0, 1.0));
264        }
265        if near_right && near_bottom && near_front {
266            result.push(new_site(-1.0, 1.0, -1.0));
267        }
268        if near_right && near_bottom && near_back {
269            result.push(new_site(-1.0, 1.0, 1.0));
270        }
271        if near_left && near_top && near_front {
272            result.push(new_site(1.0, -1.0, -1.0));
273        }
274        if near_left && near_top && near_back {
275            result.push(new_site(1.0, -1.0, 1.0));
276        }
277        if near_left && near_bottom && near_front {
278            result.push(new_site(1.0, 1.0, -1.0));
279        }
280        if near_left && near_bottom && near_back {
281            result.push(new_site(1.0, 1.0, 1.0));
282        }
283
284        result
285    }
286}
287
288#[cfg(test)]
289mod tests {
290    use super::*;
291    use crate::property::Point;
292
293    use approxim::assert_relative_eq;
294    use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
295
296    const N_SAMPLES: usize = 1024;
297
298    mod cuboid_2 {
299        use super::*;
300
301        #[test]
302        fn maximum_allowable() {
303            let cuboid = Hypercuboid {
304                edge_lengths: [
305                    10.0.try_into()
306                        .expect("hard-coded constant should be positive"),
307                    6.0.try_into()
308                        .expect("hard-coded constant should be positive"),
309                ],
310            };
311
312            assert_eq!(cuboid.maximum_allowable_interaction_range(), 3.0);
313
314            let cuboid = Hypercuboid {
315                edge_lengths: [
316                    4.0.try_into()
317                        .expect("hard-coded constant should be positive"),
318                    6.0.try_into()
319                        .expect("hard-coded constant should be positive"),
320                ],
321            };
322
323            assert_eq!(cuboid.maximum_allowable_interaction_range(), 2.0);
324
325            let cuboid = Hypercuboid {
326                edge_lengths: [
327                    100.0
328                        .try_into()
329                        .expect("hard-coded constant should be positive"),
330                    18.0.try_into()
331                        .expect("hard-coded constant should be positive"),
332                ],
333            };
334
335            assert_eq!(cuboid.maximum_allowable_interaction_range(), 9.0);
336        }
337
338        #[test]
339        fn wrap() {
340            let cuboid = Hypercuboid {
341                edge_lengths: [
342                    20.0.try_into()
343                        .expect("hard-coded constant should be positive"),
344                    20.0.try_into()
345                        .expect("hard-coded constant should be positive"),
346                ],
347            };
348
349            let periodic = Periodic::new(0.0, cuboid).expect("hard-coded range should be valid");
350
351            let point = Point::new([5.0, 3.0].into());
352            assert_eq!(periodic.wrap(point), Ok(point));
353
354            let point = Point::new([-10.0, -10.0].into());
355            assert_eq!(periodic.wrap(point), Ok(point));
356
357            let point = Point::new([10.0_f64.next_down(), 10.0_f64.next_down()].into());
358            assert_eq!(periodic.wrap(point), Ok(point));
359
360            let point = Point::new([10.0, 10.0].into());
361            assert_eq!(periodic.wrap(point), Ok(Point::new([-10.0, -10.0].into())));
362
363            let point = Point::new([20.0, 20.0].into());
364            assert_eq!(periodic.wrap(point), Ok(Point::new([0.0, 0.0].into())));
365
366            let point = Point::new([30.0, 30.0].into());
367            assert_eq!(periodic.wrap(point), Ok(Point::new([-10.0, -10.0].into())));
368
369            let point = Point::new([25.0, -35.0].into());
370            assert_eq!(periodic.wrap(point), Ok(Point::new([5.0, 5.0].into())));
371        }
372
373        #[test]
374        fn no_ghosts() {
375            let cuboid = Hypercuboid {
376                edge_lengths: [
377                    20.0.try_into()
378                        .expect("hard-coded constant should be positive"),
379                    10.0.try_into()
380                        .expect("hard-coded constant should be positive"),
381                ],
382            };
383
384            let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid");
385
386            let inner = Hypercuboid {
387                edge_lengths: [
388                    18.0.try_into()
389                        .expect("hard-coded constant should be positive"),
390                    8.0.try_into()
391                        .expect("hard-coded constant should be positive"),
392                ],
393            };
394
395            let mut rng = StdRng::seed_from_u64(1);
396            for _ in 0..N_SAMPLES {
397                let point = inner.sample(&mut rng);
398                let ghosts = periodic.generate_ghosts(&Point::new(point));
399                assert!(ghosts.is_empty());
400            }
401        }
402
403        #[test]
404        fn ghosts() {
405            let cuboid = Hypercuboid {
406                edge_lengths: [
407                    20.0.try_into()
408                        .expect("hard-coded constant should be positive"),
409                    10.0.try_into()
410                        .expect("hard-coded constant should be positive"),
411                ],
412            };
413
414            let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid");
415
416            // no ghosts for points outside the boundary
417            let ghosts = periodic.generate_ghosts(&Point::new([10.5, 0.0].into()));
418            assert!(ghosts.is_empty());
419            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 5.5].into()));
420            assert!(ghosts.is_empty());
421
422            // edges
423            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0].into()));
424            assert_eq!(ghosts.len(), 1);
425            assert_relative_eq!(ghosts[0].position, [-10.5, 0.0].into());
426
427            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0].into()));
428            assert_eq!(ghosts.len(), 1);
429            assert_relative_eq!(ghosts[0].position, [10.5, 0.0].into());
430
431            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5].into()));
432            assert_eq!(ghosts.len(), 1);
433            assert_relative_eq!(ghosts[0].position, [0.0, -5.5].into());
434
435            let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5].into()));
436            assert_eq!(ghosts.len(), 1);
437            assert_relative_eq!(ghosts[0].position, [0.0, 5.5].into());
438
439            // vertices
440            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5].into()));
441            assert_eq!(ghosts.len(), 3);
442            assert_relative_eq!(ghosts[0].position, [-10.5, 4.5].into());
443            assert_relative_eq!(ghosts[1].position, [9.5, -5.5].into());
444            assert_relative_eq!(ghosts[2].position, [-10.5, -5.5].into());
445
446            let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5].into()));
447            assert_eq!(ghosts.len(), 3);
448            assert_relative_eq!(ghosts[0].position, [-10.5, -4.5].into());
449            assert_relative_eq!(ghosts[1].position, [9.5, 5.5].into());
450            assert_relative_eq!(ghosts[2].position, [-10.5, 5.5].into());
451
452            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5].into()));
453            assert_eq!(ghosts.len(), 3);
454            assert_relative_eq!(ghosts[0].position, [10.5, 4.5].into());
455            assert_relative_eq!(ghosts[1].position, [-9.5, -5.5].into());
456            assert_relative_eq!(ghosts[2].position, [10.5, -5.5].into());
457
458            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5].into()));
459            assert_eq!(ghosts.len(), 3);
460            assert_relative_eq!(ghosts[0].position, [10.5, -4.5].into());
461            assert_relative_eq!(ghosts[1].position, [-9.5, 5.5].into());
462            assert_relative_eq!(ghosts[2].position, [10.5, 5.5].into());
463        }
464    }
465
466    mod cuboid_3 {
467        use super::*;
468
469        #[test]
470        fn maximum_allowable() {
471            let cuboid = Hypercuboid {
472                edge_lengths: [
473                    10.0.try_into()
474                        .expect("hard-coded constant should be positive"),
475                    6.0.try_into()
476                        .expect("hard-coded constant should be positive"),
477                    4.0.try_into()
478                        .expect("hard-coded constant should be positive"),
479                ],
480            };
481
482            assert_eq!(cuboid.maximum_allowable_interaction_range(), 2.0);
483
484            let cuboid = Hypercuboid {
485                edge_lengths: [
486                    6.0.try_into()
487                        .expect("hard-coded constant should be positive"),
488                    8.0.try_into()
489                        .expect("hard-coded constant should be positive"),
490                    10.0.try_into()
491                        .expect("hard-coded constant should be positive"),
492                ],
493            };
494
495            assert_eq!(cuboid.maximum_allowable_interaction_range(), 3.0);
496
497            let cuboid = Hypercuboid {
498                edge_lengths: [
499                    18.0.try_into()
500                        .expect("hard-coded constant should be positive"),
501                    8.0.try_into()
502                        .expect("hard-coded constant should be positive"),
503                    10.0.try_into()
504                        .expect("hard-coded constant should be positive"),
505                ],
506            };
507
508            assert_eq!(cuboid.maximum_allowable_interaction_range(), 4.0);
509        }
510
511        #[test]
512        fn wrap() {
513            let cuboid = Hypercuboid {
514                edge_lengths: [
515                    20.0.try_into()
516                        .expect("hard-coded constant should be positive"),
517                    20.0.try_into()
518                        .expect("hard-coded constant should be positive"),
519                    20.0.try_into()
520                        .expect("hard-coded constant should be positive"),
521                ],
522            };
523
524            let periodic = Periodic::new(0.0, cuboid).expect("hard-coded range should be valid");
525
526            let point = Point::new([5.0, 3.0, 8.0].into());
527            assert_eq!(periodic.wrap(point), Ok(point));
528
529            let point = Point::new([-10.0, -10.0, -10.0].into());
530            assert_eq!(periodic.wrap(point), Ok(point));
531
532            let point = Point::new(
533                [
534                    10.0_f64.next_down(),
535                    10.0_f64.next_down(),
536                    10.0_f64.next_down(),
537                ]
538                .into(),
539            );
540            assert_eq!(periodic.wrap(point), Ok(point));
541
542            let point = Point::new([10.0, 10.0, 10.0].into());
543            assert_eq!(
544                periodic.wrap(point),
545                Ok(Point::new([-10.0, -10.0, -10.0].into()))
546            );
547
548            let point = Point::new([20.0, 20.0, 20.0].into());
549            assert_eq!(periodic.wrap(point), Ok(Point::new([0.0, 0.0, 0.0].into())));
550
551            let point = Point::new([30.0, 30.0, 30.0].into());
552            assert_eq!(
553                periodic.wrap(point),
554                Ok(Point::new([-10.0, -10.0, -10.0].into()))
555            );
556
557            let point = Point::new([25.0, -35.0, 55.0].into());
558            assert_eq!(
559                periodic.wrap(point),
560                Ok(Point::new([5.0, 5.0, -5.0].into()))
561            );
562        }
563
564        #[test]
565        fn no_ghosts() {
566            let cuboid = Hypercuboid {
567                edge_lengths: [
568                    40.0.try_into()
569                        .expect("hard-coded constant should be positive"),
570                    20.0.try_into()
571                        .expect("hard-coded constant should be positive"),
572                    10.0.try_into()
573                        .expect("hard-coded constant should be positive"),
574                ],
575            };
576
577            let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid");
578
579            let inner = Hypercuboid {
580                edge_lengths: [
581                    38.0.try_into()
582                        .expect("hard-coded constant should be positive"),
583                    18.0.try_into()
584                        .expect("hard-coded constant should be positive"),
585                    8.0.try_into()
586                        .expect("hard-coded constant should be positive"),
587                ],
588            };
589
590            let mut rng = StdRng::seed_from_u64(1);
591            for _ in 0..N_SAMPLES {
592                let point = inner.sample(&mut rng);
593                let ghosts = periodic.generate_ghosts(&Point::new(point));
594                assert!(ghosts.is_empty());
595            }
596        }
597
598        #[expect(clippy::too_many_lines, reason = "There are many cases to test.")]
599        #[test]
600        fn ghosts() {
601            let cuboid = Hypercuboid {
602                edge_lengths: [
603                    20.0.try_into()
604                        .expect("hard-coded constant should be positive"),
605                    10.0.try_into()
606                        .expect("hard-coded constant should be positive"),
607                    40.0.try_into()
608                        .expect("hard-coded constant should be positive"),
609                ],
610            };
611
612            let periodic = Periodic::new(1.0, cuboid).expect("hard-coded range should be valid");
613
614            // no ghosts for points outside the boundary
615            let ghosts = periodic.generate_ghosts(&Point::new([10.5, 0.0, 0.0].into()));
616            assert!(ghosts.is_empty());
617            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 5.5, 0.0].into()));
618            assert!(ghosts.is_empty());
619
620            // faces
621            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, 0.0].into()));
622            assert_eq!(ghosts.len(), 1);
623            assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, 0.0].into());
624
625            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, 0.0].into()));
626            assert_eq!(ghosts.len(), 1);
627            assert_relative_eq!(ghosts[0].position, [10.5, 0.0, 0.0].into());
628
629            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, 0.0].into()));
630            assert_eq!(ghosts.len(), 1);
631            assert_relative_eq!(ghosts[0].position, [0.0, -5.5, 0.0].into());
632
633            let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, 0.0].into()));
634            assert_eq!(ghosts.len(), 1);
635            assert_relative_eq!(ghosts[0].position, [0.0, 5.5, 0.0].into());
636
637            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 0.0, 19.5].into()));
638            assert_eq!(ghosts.len(), 1);
639            assert_relative_eq!(ghosts[0].position, [0.0, 0.0, -20.5].into());
640
641            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 0.0, -19.5].into()));
642            assert_eq!(ghosts.len(), 1);
643            assert_relative_eq!(ghosts[0].position, [0.0, 0.0, 20.5].into());
644
645            // edges
646            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, 0.0].into()));
647            assert_eq!(ghosts.len(), 3);
648            assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, 0.0].into());
649            assert_relative_eq!(ghosts[1].position, [9.5, -5.5, 0.0].into());
650            assert_relative_eq!(ghosts[2].position, [-10.5, -5.5, 0.0].into());
651
652            let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, 0.0].into()));
653            assert_eq!(ghosts.len(), 3);
654            assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, 0.0].into());
655            assert_relative_eq!(ghosts[1].position, [9.5, 5.5, 0.0].into());
656            assert_relative_eq!(ghosts[2].position, [-10.5, 5.5, 0.0].into());
657
658            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, 0.0].into()));
659            assert_eq!(ghosts.len(), 3);
660            assert_relative_eq!(ghosts[0].position, [10.5, 4.5, 0.0].into());
661            assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, 0.0].into());
662            assert_relative_eq!(ghosts[2].position, [10.5, -5.5, 0.0].into());
663
664            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, 0.0].into()));
665            assert_eq!(ghosts.len(), 3);
666            assert_relative_eq!(ghosts[0].position, [10.5, -4.5, 0.0].into());
667            assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, 0.0].into());
668            assert_relative_eq!(ghosts[2].position, [10.5, 5.5, 0.0].into());
669
670            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, 19.5].into()));
671            assert_eq!(ghosts.len(), 3);
672            assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, 19.5].into());
673            assert_relative_eq!(ghosts[1].position, [9.5, 0.0, -20.5].into());
674            assert_relative_eq!(ghosts[2].position, [-10.5, 0.0, -20.5].into());
675
676            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, 19.5].into()));
677            assert_eq!(ghosts.len(), 3);
678            assert_relative_eq!(ghosts[0].position, [10.5, 0.0, 19.5].into());
679            assert_relative_eq!(ghosts[1].position, [-9.5, 0.0, -20.5].into());
680            assert_relative_eq!(ghosts[2].position, [10.5, 0.0, -20.5].into());
681
682            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, 19.5].into()));
683            assert_eq!(ghosts.len(), 3);
684            assert_relative_eq!(ghosts[0].position, [0.0, -5.5, 19.5].into());
685            assert_relative_eq!(ghosts[1].position, [0.0, 4.5, -20.5].into());
686            assert_relative_eq!(ghosts[2].position, [0.0, -5.5, -20.5].into());
687
688            let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, 19.5].into()));
689            assert_eq!(ghosts.len(), 3);
690            assert_relative_eq!(ghosts[0].position, [0.0, 5.5, 19.5].into());
691            assert_relative_eq!(ghosts[1].position, [0.0, -4.5, -20.5].into());
692            assert_relative_eq!(ghosts[2].position, [0.0, 5.5, -20.5].into());
693
694            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 0.0, -19.5].into()));
695            assert_eq!(ghosts.len(), 3);
696            assert_relative_eq!(ghosts[0].position, [-10.5, 0.0, -19.5].into());
697            assert_relative_eq!(ghosts[1].position, [9.5, 0.0, 20.5].into());
698            assert_relative_eq!(ghosts[2].position, [-10.5, 0.0, 20.5].into());
699
700            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 0.0, -19.5].into()));
701            assert_eq!(ghosts.len(), 3);
702            assert_relative_eq!(ghosts[0].position, [10.5, 0.0, -19.5].into());
703            assert_relative_eq!(ghosts[1].position, [-9.5, 0.0, 20.5].into());
704            assert_relative_eq!(ghosts[2].position, [10.5, 0.0, 20.5].into());
705
706            let ghosts = periodic.generate_ghosts(&Point::new([0.0, 4.5, -19.5].into()));
707            assert_eq!(ghosts.len(), 3);
708            assert_relative_eq!(ghosts[0].position, [0.0, -5.5, -19.5].into());
709            assert_relative_eq!(ghosts[1].position, [0.0, 4.5, 20.5].into());
710            assert_relative_eq!(ghosts[2].position, [0.0, -5.5, 20.5].into());
711
712            let ghosts = periodic.generate_ghosts(&Point::new([0.0, -4.5, -19.5].into()));
713            assert_eq!(ghosts.len(), 3);
714            assert_relative_eq!(ghosts[0].position, [0.0, 5.5, -19.5].into());
715            assert_relative_eq!(ghosts[1].position, [0.0, -4.5, 20.5].into());
716            assert_relative_eq!(ghosts[2].position, [0.0, 5.5, 20.5].into());
717
718            // vertices
719            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, 19.5].into()));
720            assert_eq!(ghosts.len(), 7);
721            assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, 19.5].into());
722            assert_relative_eq!(ghosts[1].position, [9.5, -5.5, 19.5].into());
723            assert_relative_eq!(ghosts[2].position, [9.5, 4.5, -20.5].into());
724            assert_relative_eq!(ghosts[3].position, [-10.5, -5.5, 19.5].into());
725            assert_relative_eq!(ghosts[4].position, [-10.5, 4.5, -20.5].into());
726            assert_relative_eq!(ghosts[5].position, [9.5, -5.5, -20.5].into());
727            assert_relative_eq!(ghosts[6].position, [-10.5, -5.5, -20.5].into());
728
729            let ghosts = periodic.generate_ghosts(&Point::new([9.5, 4.5, -19.5].into()));
730            assert_eq!(ghosts.len(), 7);
731            assert_relative_eq!(ghosts[0].position, [-10.5, 4.5, -19.5].into());
732            assert_relative_eq!(ghosts[1].position, [9.5, -5.5, -19.5].into());
733            assert_relative_eq!(ghosts[2].position, [9.5, 4.5, 20.5].into());
734            assert_relative_eq!(ghosts[3].position, [-10.5, -5.5, -19.5].into());
735            assert_relative_eq!(ghosts[4].position, [-10.5, 4.5, 20.5].into());
736            assert_relative_eq!(ghosts[5].position, [9.5, -5.5, 20.5].into());
737            assert_relative_eq!(ghosts[6].position, [-10.5, -5.5, 20.5].into());
738
739            let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, 19.5].into()));
740            assert_eq!(ghosts.len(), 7);
741            assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, 19.5].into());
742            assert_relative_eq!(ghosts[1].position, [9.5, 5.5, 19.5].into());
743            assert_relative_eq!(ghosts[2].position, [9.5, -4.5, -20.5].into());
744            assert_relative_eq!(ghosts[3].position, [-10.5, 5.5, 19.5].into());
745            assert_relative_eq!(ghosts[4].position, [-10.5, -4.5, -20.5].into());
746            assert_relative_eq!(ghosts[5].position, [9.5, 5.5, -20.5].into());
747            assert_relative_eq!(ghosts[6].position, [-10.5, 5.5, -20.5].into());
748
749            let ghosts = periodic.generate_ghosts(&Point::new([9.5, -4.5, -19.5].into()));
750            assert_eq!(ghosts.len(), 7);
751            assert_relative_eq!(ghosts[0].position, [-10.5, -4.5, -19.5].into());
752            assert_relative_eq!(ghosts[1].position, [9.5, 5.5, -19.5].into());
753            assert_relative_eq!(ghosts[2].position, [9.5, -4.5, 20.5].into());
754            assert_relative_eq!(ghosts[3].position, [-10.5, 5.5, -19.5].into());
755            assert_relative_eq!(ghosts[4].position, [-10.5, -4.5, 20.5].into());
756            assert_relative_eq!(ghosts[5].position, [9.5, 5.5, 20.5].into());
757            assert_relative_eq!(ghosts[6].position, [-10.5, 5.5, 20.5].into());
758
759            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, 19.5].into()));
760            assert_eq!(ghosts.len(), 7);
761            assert_relative_eq!(ghosts[0].position, [10.5, 4.5, 19.5].into());
762            assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, 19.5].into());
763            assert_relative_eq!(ghosts[2].position, [-9.5, 4.5, -20.5].into());
764            assert_relative_eq!(ghosts[3].position, [10.5, -5.5, 19.5].into());
765            assert_relative_eq!(ghosts[4].position, [10.5, 4.5, -20.5].into());
766            assert_relative_eq!(ghosts[5].position, [-9.5, -5.5, -20.5].into());
767            assert_relative_eq!(ghosts[6].position, [10.5, -5.5, -20.5].into());
768
769            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, 4.5, -19.5].into()));
770            assert_eq!(ghosts.len(), 7);
771            assert_relative_eq!(ghosts[0].position, [10.5, 4.5, -19.5].into());
772            assert_relative_eq!(ghosts[1].position, [-9.5, -5.5, -19.5].into());
773            assert_relative_eq!(ghosts[2].position, [-9.5, 4.5, 20.5].into());
774            assert_relative_eq!(ghosts[3].position, [10.5, -5.5, -19.5].into());
775            assert_relative_eq!(ghosts[4].position, [10.5, 4.5, 20.5].into());
776            assert_relative_eq!(ghosts[5].position, [-9.5, -5.5, 20.5].into());
777            assert_relative_eq!(ghosts[6].position, [10.5, -5.5, 20.5].into());
778
779            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, 19.5].into()));
780            assert_eq!(ghosts.len(), 7);
781            assert_relative_eq!(ghosts[0].position, [10.5, -4.5, 19.5].into());
782            assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, 19.5].into());
783            assert_relative_eq!(ghosts[2].position, [-9.5, -4.5, -20.5].into());
784            assert_relative_eq!(ghosts[3].position, [10.5, 5.5, 19.5].into());
785            assert_relative_eq!(ghosts[4].position, [10.5, -4.5, -20.5].into());
786            assert_relative_eq!(ghosts[5].position, [-9.5, 5.5, -20.5].into());
787            assert_relative_eq!(ghosts[6].position, [10.5, 5.5, -20.5].into());
788
789            let ghosts = periodic.generate_ghosts(&Point::new([-9.5, -4.5, -19.5].into()));
790            assert_eq!(ghosts.len(), 7);
791            assert_relative_eq!(ghosts[0].position, [10.5, -4.5, -19.5].into());
792            assert_relative_eq!(ghosts[1].position, [-9.5, 5.5, -19.5].into());
793            assert_relative_eq!(ghosts[2].position, [-9.5, -4.5, 20.5].into());
794            assert_relative_eq!(ghosts[3].position, [10.5, 5.5, -19.5].into());
795            assert_relative_eq!(ghosts[4].position, [10.5, -4.5, 20.5].into());
796            assert_relative_eq!(ghosts[5].position, [-9.5, 5.5, 20.5].into());
797            assert_relative_eq!(ghosts[6].position, [10.5, 5.5, 20.5].into());
798        }
799    }
800}