Struct Cell 

pub struct Cell { /* private fields */ }

Periodic parallelepiped: columns of H, origin, and an ortho flag.

let cell = Cell::ortho(10.0, 11.0, 12.0)?;
assert!(cell.is_ortho());
let sheared = Cell::from_vectors(
    [10.0, 0.0, 0.0],
    [5.0, 8.66, 0.0],
    [0.0, 0.0, 10.0],
    [0.0, 0.0, 0.0],
)?;
assert!(!sheared.is_ortho());

Implementations

impl Cell

pub fn ortho(lx: f64, ly: f64, lz: f64) -> Result<Cell, Error>

Diagonal box with origin at zero.

Sets Self::is_ortho so Self::dist2 and Self::lattice_shift skip the two 3x3 matvecs.

Examples found in repository

examples/knn_scale.rs (line 28)

fn cubic(nside: usize) -> (Vec<[f64; 3]>, Cell) {
    let boxl = nside as f64 * 3.125;
    let cell = Cell::ortho(boxl, boxl, boxl).expect("box");
    let mut xyz = Vec::with_capacity(nside * nside * nside);
    let a = 3.125;
    for iz in 0..nside {
        for iy in 0..nside {
            for ix in 0..nside {
                xyz.push([ix as f64 * a, iy as f64 * a, iz as f64 * a]);
            }
        }
    }
    (xyz, cell)
}

examples/knn_hot.rs (line 59)

fn slab(n: usize) -> (Vec<[f64; 3]>, Cell) {
    let nside = nside_of(n);
    let film = nside as f64 * 3.125;
    let cell = Cell::ortho(film, film, film * 2.0).expect("slab");
    let mut xyz = Vec::with_capacity(n);
    for i in 0..n {
        let ix = i % nside;
        let iy = (i / nside) % nside;
        let iz = i / (nside * nside);
        xyz.push([
            (ix as f64 + 0.5) / nside as f64 * film,
            (iy as f64 + 0.5) / nside as f64 * film,
            (iz as f64 + 0.5) / nside as f64 * film,
        ]);
    }
    (xyz, cell)
}

/// Tall orthorhombic column at the cubic spacing. Points occupy a band
/// and the rest of z is vacuum, as in a seams nanotube trajectory.
fn column(n: usize) -> (Vec<[f64; 3]>, Cell) {
    let nside = nside_of(n);
    let side = nside as f64 * 3.125;
    let cell = Cell::ortho(side, side, side * 4.0).expect("column");
    let mut xyz = Vec::with_capacity(n);
    for i in 0..n {
        let ix = i % nside;
        let iy = (i / nside) % nside;
        let iz = i / (nside * nside);
        xyz.push([
            (ix as f64 + 0.5) / nside as f64 * side,
            (iy as f64 + 0.5) / nside as f64 * side,
            (iz as f64 + 0.5) / nside as f64 * side + side * 1.5,
        ]);
    }
    (xyz, cell)
}

fn system(shape: &str, n: usize) -> (Vec<[f64; 3]>, Cell) {
    match shape {
        "cubic" => {
            let boxl = 50.0;
            (
                lattice(n, boxl),
                Cell::ortho(boxl, boxl, boxl).expect("box"),
            )
        }
        "slab" => slab(n),
        "column" => column(n),
        "hex" => {
            // rgsaddle hexagonal cell, scaled to the cubic spacing.
            let cell = Cell::from_vectors(
                [50.0, 0.0, 0.0],
                [25.0, 43.30127018922193, 0.0],
                [0.0, 0.0, 55.55555555555556],
                [0.0, 0.0, 0.0],
            )
            .expect("hex");
            (frac_grid(&cell, n), cell)
        }
        "tilt" => {
            // genice sH dump (BOX BOUNDS xy xz yz), scaled to the cubic spacing.
            let s = 2.85;
            let cell = Cell::from_lammps(
                -6.2106057 * s,
                24.8424228 * s,
                0.0,
                10.7570846 * s,
                0.0,
                20.131291 * s,
                -6.2106057 * s,
                0.0,
                0.0,
            )
            .expect("tilt");
            (frac_grid(&cell, n), cell)
        }
        other => panic!("unknown KNN_HOT_SHAPE {other}"),
    }
}

examples/pairs_scale.rs (line 71)

fn main() {
    let n = env_usize("PAIRS_N", 4096);
    let reps = env_usize("PAIRS_REPS", 5).max(1);
    let half = std::env::var("PAIRS_HALF").ok().as_deref() == Some("1");
    let hint = std::env::var("PAIRS_HINT")
        .ok()
        .and_then(|s| s.parse().ok());
    let threads = std::env::var("RAYON_NUM_THREADS").unwrap_or_else(|_| "default".to_string());
    let mut xyz = fill(n);
    if let Some(seed) = std::env::var("PAIRS_SHUFFLE")
        .ok()
        .and_then(|s| s.parse().ok())
    {
        shuffle(&mut xyz, seed);
    }
    let cell = Cell::ortho(18.0, 18.0, 18.0).expect("box");
    linkcell::pop_engage();
    linkcell::pop_prepare();
    let mut pairs = Vec::new();
    let mut early = Vec::new();
    for _ in 0..3 {
        let t = Instant::now();
        pairs = pairs_within(&xyz, &cell, 4.0, None, hint, half).expect("pairs");
        early.push(t.elapsed().as_secs_f64() * 1e3);
    }
    linkcell::pop_reset();
    let t0 = Instant::now();
    let mut acc = 0usize;
    for _ in 0..reps {
        let got = pairs_within(&xyz, &cell, 4.0, None, hint, half).expect("pairs");
        acc += got.len();
    }
    let wall_ns = t0.elapsed().as_nanos() as u64;
    let ms = wall_ns as f64 / 1e6;
    println!(
        "n={} half={} hint={:?} pairs={} cold={:.3} call2={:.3} call3={:.3} reps={} ms_per={:.4} acc={}",
        xyz.len(),
        half as u8,
        hint,
        pairs.len(),
        early[0],
        early[1],
        early[2],
        reps,
        ms / reps as f64,
        acc
    );
    let useful = linkcell::pop_snapshot();
    let (lb, ce, pe) = linkcell::pop_efficiencies(&useful, wall_ns);
    let list = useful
        .iter()
        .map(|ns| ns.to_string())
        .collect::<Vec<_>>()
        .join(",");
    println!(
        "pop shape=pairs-{n}-h{} threads={threads} n={n} reps={reps} wall_ns={wall_ns} useful_ns={list} lb={lb:.4} ce={ce:.4} pe={pe:.4}",
        half as u8
    );
}

pub fn ortho_origin( lx: f64, ly: f64, lz: f64, origin: [f64; 3], ) -> Result<Cell, Error>

Diagonal box with an explicit dump-cell origin.

pub fn from_vectors( a: [f64; 3], b: [f64; 3], c: [f64; 3], origin: [f64; 3], ) -> Result<Cell, Error>

Parallelepiped from lattice vectors a, b, c and an origin.

Vectors are the columns of H. A non-diagonal H clears Self::is_ortho.

Inlined so the caller writes the cell once. An outlined return stores it and then copies it.

Examples found in repository

examples/knn_scale.rs (lines 45-50)

fn tilt(nside: usize) -> (Vec<[f64; 3]>, Cell) {
    let lx = nside as f64 * 3.125;
    let cell = Cell::from_vectors(
        [lx, 0.0, 0.0],
        [lx * 0.2, lx * 0.9, 0.0],
        [lx * 0.05, lx * -0.08, lx * 0.95],
        [0.0, 0.0, 0.0],
    )
    .expect("tilt");
    let mut xyz = Vec::with_capacity(nside * nside * nside);
    let nf = nside as f64;
    for iz in 0..nside {
        for iy in 0..nside {
            for ix in 0..nside {
                let s = [
                    (ix as f64 + 0.5) / nf,
                    (iy as f64 + 0.5) / nf,
                    (iz as f64 + 0.5) / nf,
                ];
                xyz.push(cell.cartesian(s));
            }
        }
    }
    (xyz, cell)
}

examples/knn_hot.rs (lines 107-112)

fn system(shape: &str, n: usize) -> (Vec<[f64; 3]>, Cell) {
    match shape {
        "cubic" => {
            let boxl = 50.0;
            (
                lattice(n, boxl),
                Cell::ortho(boxl, boxl, boxl).expect("box"),
            )
        }
        "slab" => slab(n),
        "column" => column(n),
        "hex" => {
            // rgsaddle hexagonal cell, scaled to the cubic spacing.
            let cell = Cell::from_vectors(
                [50.0, 0.0, 0.0],
                [25.0, 43.30127018922193, 0.0],
                [0.0, 0.0, 55.55555555555556],
                [0.0, 0.0, 0.0],
            )
            .expect("hex");
            (frac_grid(&cell, n), cell)
        }
        "tilt" => {
            // genice sH dump (BOX BOUNDS xy xz yz), scaled to the cubic spacing.
            let s = 2.85;
            let cell = Cell::from_lammps(
                -6.2106057 * s,
                24.8424228 * s,
                0.0,
                10.7570846 * s,
                0.0,
                20.131291 * s,
                -6.2106057 * s,
                0.0,
                0.0,
            )
            .expect("tilt");
            (frac_grid(&cell, n), cell)
        }
        other => panic!("unknown KNN_HOT_SHAPE {other}"),
    }
}

pub fn from_ase(rows: [[f64; 3]; 3]) -> Result<Cell, Error>

ASE-style 3x3 cell: rows are lattice vectors a, b, c.

Origin is zero. Matches atoms.cell in ASE.

pub fn from_con(rows: [[f64; 3]; 3]) -> Result<Cell, Error>

CON lattice: rows are lattice vectors a, b, c, origin zero.

Same storage as readcon lattice_vectors and an eOn CON header that already holds the 3x3 basis.

pub fn from_con_box(boxl: [f64; 3], angles_deg: [f64; 3]) -> Result<Cell, Error>

CON / crystallographic lengths and angles (degrees).

alpha is the angle at the origin between b and c, beta between a and c, gamma between a and b. The basis is the standard reduced triclinic frame used by ASE and eOn.

pub fn from_vesin(rows: [[f64; 3]; 3]) -> Result<Cell, Error>

vesin box: rows are lattice vectors a, b, c, origin zero.

pub fn from_lammps( xlo: f64, xhi: f64, ylo: f64, yhi: f64, zlo: f64, zhi: f64, xy: f64, xz: f64, yz: f64, ) -> Result<Cell, Error>

LAMMPS restricted-triclinic box from xlo xhi ylo yhi zlo zhi and tilt factors xy, xz, yz.

Recovers the same H and origin as a dump ITEM BOX BOUNDS line whose bound lo/hi already include the tilt padding.

Examples found in repository

examples/knn_hot.rs (lines 119-129)

fn system(shape: &str, n: usize) -> (Vec<[f64; 3]>, Cell) {
    match shape {
        "cubic" => {
            let boxl = 50.0;
            (
                lattice(n, boxl),
                Cell::ortho(boxl, boxl, boxl).expect("box"),
            )
        }
        "slab" => slab(n),
        "column" => column(n),
        "hex" => {
            // rgsaddle hexagonal cell, scaled to the cubic spacing.
            let cell = Cell::from_vectors(
                [50.0, 0.0, 0.0],
                [25.0, 43.30127018922193, 0.0],
                [0.0, 0.0, 55.55555555555556],
                [0.0, 0.0, 0.0],
            )
            .expect("hex");
            (frac_grid(&cell, n), cell)
        }
        "tilt" => {
            // genice sH dump (BOX BOUNDS xy xz yz), scaled to the cubic spacing.
            let s = 2.85;
            let cell = Cell::from_lammps(
                -6.2106057 * s,
                24.8424228 * s,
                0.0,
                10.7570846 * s,
                0.0,
                20.131291 * s,
                -6.2106057 * s,
                0.0,
                0.0,
            )
            .expect("tilt");
            (frac_grid(&cell, n), cell)
        }
        other => panic!("unknown KNN_HOT_SHAPE {other}"),
    }
}

pub fn from_lammps_bounds( xspan: f64, yspan: f64, zspan: f64, xy: f64, xz: f64, yz: f64, xlo_b: f64, ylo_b: f64, zlo_b: f64, ) -> Result<Cell, Error>

LAMMPS dump bound spans plus tilts and bound lo.

xspan, yspan, zspan are xhi_bound - xlo_bound etc. xy, xz, yz are the tilt factors. xlo_b, ylo_b, zlo_b are the bound lo. Inverse of xlo_bound = xlo + min(0, xy, xz, xy+xz).

pub fn is_ortho(&self) -> bool

True when H is diagonal. Distances and lattice shifts then skip the two 3x3 matvecs and use three independent wraps.

pub fn a(&self) -> [f64; 3]

Lattice vector a (first column of H).

pub fn b(&self) -> [f64; 3]

Lattice vector b.

pub fn c(&self) -> [f64; 3]

Lattice vector c.

pub fn origin(&self) -> [f64; 3]

Dump-cell origin.

pub fn widths(&self) -> [f64; 3]

Perpendicular widths of the three faces.

pub fn h(&self) -> [[f64; 3]; 3]

Columns of H.

pub fn hinv(&self) -> [[f64; 3]; 3]

Inverse of H, stored by columns.

pub fn fractional(&self, r: [f64; 3]) -> [f64; 3]

Fractional coordinates in [0, 1).

Orthorhombic boxes divide by the three widths. The general path is s = Hinv (r - origin), then wrap.

pub fn cartesian(&self, s: [f64; 3]) -> [f64; 3]

Cartesian from fractional: r = H s + origin.

Examples found in repository

examples/knn_hot.rs (line 45)

fn frac_grid(cell: &Cell, n: usize) -> Vec<[f64; 3]> {
    let nside = ((n as f64).cbrt().ceil() as usize).max(1);
    let mut xyz = Vec::with_capacity(n);
    for i in 0..n {
        let ix = i % nside;
        let iy = (i / nside) % nside;
        let iz = i / (nside * nside);
        let s = [
            (ix as f64 + 0.5) / nside as f64,
            (iy as f64 + 0.5) / nside as f64,
            (iz as f64 + 0.5) / nside as f64,
        ];
        xyz.push(cell.cartesian(s));
    }
    xyz
}

examples/knn_scale.rs (line 62)

fn tilt(nside: usize) -> (Vec<[f64; 3]>, Cell) {
    let lx = nside as f64 * 3.125;
    let cell = Cell::from_vectors(
        [lx, 0.0, 0.0],
        [lx * 0.2, lx * 0.9, 0.0],
        [lx * 0.05, lx * -0.08, lx * 0.95],
        [0.0, 0.0, 0.0],
    )
    .expect("tilt");
    let mut xyz = Vec::with_capacity(nside * nside * nside);
    let nf = nside as f64;
    for iz in 0..nside {
        for iy in 0..nside {
            for ix in 0..nside {
                let s = [
                    (ix as f64 + 0.5) / nf,
                    (iy as f64 + 0.5) / nf,
                    (iz as f64 + 0.5) / nf,
                ];
                xyz.push(cell.cartesian(s));
            }
        }
    }
    (xyz, cell)
}

pub fn lattice_shift(&self, na: i32, nb: i32, nc: i32) -> [f64; 3]

Cartesian translation by integer lattice counts (na, nb, nc).

pub fn dist2_shifted(&self, p: [f64; 3], q: [f64; 3], shift: [f64; 3]) -> f64

Squared Cartesian distance after applying a lattice shift to q.

This is not a minimum-image wrap. The caller supplies the image via Self::lattice_shift.

pub fn displacement(&self, p: [f64; 3], q: [f64; 3]) -> [f64; 3]

Minimum-image displacement from p to q: wrap(q - p).

Orthorhombic boxes wrap each axis independently. The general path is ds = wrap(Hinv (q - p)), then dr = H ds.

pub fn dist2(&self, p: [f64; 3], q: [f64; 3]) -> f64

Squared minimum-image distance.

pub fn displacement_cartesian(&self, p: [f64; 3], q: [f64; 3]) -> [f64; 3]

Euclidean MIC by checking the 27 nearest lattice images.

Self::displacement wraps in fractional coordinates. That is the LAMMPS / eOn / ASE convention and is the Cartesian nearest image on a reduced (restricted-triclinic) cell. A highly skewed H can make a neighbouring image shorter in Cartesian space. This walk is the check those codes omit.

pub fn is_restricted(&self) -> bool

True when a is along x and b lies in the xy plane.

This is the GROMACS / LAMMPS / HOOMD restricted-triclinic frame. Those engines refuse or rotate any other orientation before wrap.

pub fn tilts_reduced(&self) -> bool

True when the cell is restricted and the tilts sit inside the GROMACS correct_box / LAMMPS half-box limits |xy| <= lx/2, |xz| <= lx/2, |yz| <= ly/2.

On that domain the fractional wrap is the Euclidean MIC, which is why LAMMPS Domain::minimum_image and HOOMD BoxDim::minImage never search extra images.

pub fn reduce_tilts(&self) -> Result<Cell, Error>

Unimodular tilt reduction: GROMACS correct_box.

Subtracts integer lattice vectors so |xy|, |xz|, |yz| sit inside half the corresponding edge. Same Cartesian lattice, new basis. No-op when Self::tilts_reduced is already true. A general orientation is rotated into the restricted frame first (LAMMPS define_general_triclinic).

pub fn to_restricted(&self) -> Result<Cell, Error>

Rotate a general parallelepiped onto the restricted-triclinic frame (LAMMPS general-to-restricted). Identity when already restricted.

pub fn displacement_euclidean(&self, p: [f64; 3], q: [f64; 3]) -> [f64; 3]

Euclidean nearest image, including displacements longer than the cell.

An orthorhombic box is the per-axis wrap. Otherwise Smith, CCP5 Newsletter 1989: the engine wrap is nearest when its length is strictly below half the smallest face altitude. That is the cutoff regime of GROMACS, LAMMPS, and HOOMD. Linkcell k-NN has no cutoff. A hex-prism body diagonal is an engine wrap that is longer than another image.

Past that test the Selling superbasis of H is handed to McKilliam, Grant, and Clarkson, SIAM J. Discrete Math. 28, 1405 (2014). The superbasis is cached on the calling thread, keyed by H, and is not stored in the cell. A tie keeps the engine vector. Self::displacement stays that engine vector on the caller’s H.

pub fn dist2_shifted_many( &self, p: [f64; 3], qs: &[[f64; 3]], shift: [f64; 3], out: &mut [f64], ) -> Result<(), Error>

Squared distances from p to each q, with one lattice shift already applied to every q.

Rapaport’s linked-cell pair: |q + shift - p|^2. No wrap and no Hinv. out has one entry per q.

pub fn dist2_shifted_indexed( &self, p: [f64; 3], positions: &[[f64; 3]], indices: &[usize], shift: [f64; 3], out: &mut [f64], ) -> Result<(), Error>

Squared distances from p to positions named by indices, with one lattice shift already applied.

A linked-cell bin is an index list into the neighbour array, not a contiguous slice. This gathers a chunk of those positions and runs the same kernel as Self::dist2_shifted_many. The gather copies every point before the kernel runs, and on a long index list that copy was slower than the inlined subtract, so the production walk keeps the subtract.

out has one entry per index. An index past positions panics.

pub fn dist2_euclidean(&self, p: [f64; 3], q: [f64; 3]) -> f64

Squared Euclidean MIC distance.

pub fn fractional_matches_cartesian(&self) -> bool

True when fractional wrap matches the 27-image check on a grid of short Cartesian probes (the cutoff-scale displacements MD engines actually wrap).

A hex prism at the tilt limit can still make a body-diagonal fractional wrap longer than a neighbouring image. Production codes ignore that because the cutoff is well below half the long diagonal.

Trait Implementations

impl Clone for Cell

fn clone(&self) -> Cell

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Copy for Cell

impl Debug for Cell

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter.

impl PartialEq for Cell

fn eq(&self, other: &Cell) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl StructuralPartialEq for Cell