readcon_core/
correspondence.rs

1//! Pair atoms across frames by column-5 `atom_id`.
2//!
3//! CON stores atoms in type-grouped order. The id is the index from before
4//! that grouping, so a reactant and a product that a writer grouped
5//! differently still name the same atom. These helpers are the map, the
6//! reorder, the band check, and the minimum-image displacement of the
7//! matched pairs.
8//!
9//! The cell is [`FrameHeader::cell_matrix`] and [`FrameHeader::pbc_effective`].
10//! The matrix is inverted here; lengths and angles are not rebuilt.
11//!
12//! ```
13//! use readcon_core::types::ConFrameBuilder;
14//! use readcon_core::{
15//!     IdentityCheck, atom_id_permutation, check_atom_id_consistency,
16//!     minimum_image_displacements, reorder_like,
17//! };
18//!
19//! fn frame(rows: &[(u64, &str, f64)]) -> readcon_core::types::ConFrame {
20//!     let mut builder = ConFrameBuilder::new([10.0, 10.0, 10.0], [90.0, 90.0, 90.0]);
21//!     for (id, symbol, x) in rows {
22//!         builder.add_atom(symbol, *x, 0.0, 0.0, [false; 3], *id, 12.0);
23//!     }
24//!     builder.build().unwrap()
25//! }
26//!
27//! let reactant = frame(&[(1, "C", 0.0), (2, "O", 1.0)]);
28//! let product = frame(&[(2, "O", 1.2), (1, "C", 0.1)]);
29//! let perm = atom_id_permutation(&reactant, &product, IdentityCheck::Strict).unwrap();
30//! assert_eq!(perm, vec![1, 0]);
31//! let ordered = reorder_like(&reactant, &product, IdentityCheck::Strict).unwrap();
32//! assert_eq!(ordered.atom_data[0].atom_id, 1);
33//! assert_eq!(ordered.atom_data[1].atom_id, 2);
34//! check_atom_id_consistency(&[reactant.clone(), product.clone()], IdentityCheck::Strict)
35//!     .unwrap();
36//! let disp = minimum_image_displacements(&reactant, &product, IdentityCheck::Strict).unwrap();
37//! assert_eq!(disp.len(), 2);
38//! assert!((disp[0][0] - 0.1).abs() < 1e-12);
39//! ```
40
41use crate::storage_dtype::{FloatArray1, FloatArray2};
42use crate::types::{Bond, ConFrame, ConFrameBuilder};
43use rustc_hash::FxHashMap;
44
45/// How strictly two atoms that share an `atom_id` must agree.
46///
47/// Symbol mismatches, duplicate ids, and missing or extra ids fail in both
48/// modes. [`IdentityCheck::Strict`] also requires the mass and the fixed
49/// mask to agree. Masses within [`ConFrameBuilder::TYPE_MASS_ABS_TOL`] agree.
50#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
51pub enum IdentityCheck {
52    /// Symbol, mass, and fixed mask agree.
53    #[default]
54    Strict,
55    /// Symbol agrees. Mass and fixed mask may differ.
56    Lenient,
57}
58
59/// Why two frames (or a later frame in a sequence) do not correspond.
60#[derive(Debug, Clone, PartialEq)]
61pub enum CorrespondenceError {
62    /// `check_atom_id_consistency` was given no frames.
63    EmptySequence,
64    /// The same id occurs twice in one frame.
65    Duplicate {
66        frame_index: usize,
67        atom_id: u64,
68        /// Index of the second occurrence.
69        atom_index: usize,
70        first_index: usize,
71    },
72    /// `atom_id` is in the reference frame and absent from this frame.
73    Missing {
74        frame_index: usize,
75        atom_id: u64,
76        /// Index of the atom in the reference frame.
77        reference_index: usize,
78    },
79    /// `atom_id` is in this frame and absent from the reference frame.
80    Extra {
81        frame_index: usize,
82        atom_id: u64,
83        atom_index: usize,
84    },
85    Symbol {
86        frame_index: usize,
87        atom_id: u64,
88        atom_index: usize,
89        reference_symbol: String,
90        found_symbol: String,
91    },
92    Mass {
93        frame_index: usize,
94        atom_id: u64,
95        atom_index: usize,
96        reference_mass: f64,
97        found_mass: f64,
98    },
99    Fixed {
100        frame_index: usize,
101        atom_id: u64,
102        atom_index: usize,
103        reference_fixed: [bool; 3],
104        found_fixed: [bool; 3],
105    },
106    /// Reference cell matrix cannot be inverted.
107    SingularCell {
108        determinant: f64,
109    },
110    NonFiniteDelta,
111    NonFiniteCoordinate {
112        frame_index: usize,
113        atom_id: u64,
114        atom_index: usize,
115    },
116    SectionLength {
117        section: &'static str,
118        rows: usize,
119        atoms: usize,
120    },
121    BondEndpoint {
122        index: u32,
123        atoms: usize,
124    },
125}
126
127impl CorrespondenceError {
128    /// Frame that failed. `None` when the failure is not a frame.
129    pub fn frame_index(&self) -> Option<usize> {
130        match self {
131            Self::EmptySequence | Self::SingularCell { .. } | Self::NonFiniteDelta => None,
132            Self::Duplicate { frame_index, .. }
133            | Self::Missing { frame_index, .. }
134            | Self::Extra { frame_index, .. }
135            | Self::Symbol { frame_index, .. }
136            | Self::Mass { frame_index, .. }
137            | Self::Fixed { frame_index, .. }
138            | Self::NonFiniteCoordinate { frame_index, .. } => Some(*frame_index),
139            Self::SectionLength { .. } | Self::BondEndpoint { .. } => None,
140        }
141    }
142
143    /// Atom id involved, when the failure names one atom.
144    pub fn atom_id(&self) -> Option<u64> {
145        match self {
146            Self::Duplicate { atom_id, .. }
147            | Self::Missing { atom_id, .. }
148            | Self::Extra { atom_id, .. }
149            | Self::Symbol { atom_id, .. }
150            | Self::Mass { atom_id, .. }
151            | Self::Fixed { atom_id, .. }
152            | Self::NonFiniteCoordinate { atom_id, .. } => Some(*atom_id),
153            _ => None,
154        }
155    }
156
157    /// Storage index of that atom.
158    ///
159    /// For [`CorrespondenceError::Missing`] this is the index in the
160    /// reference frame. Otherwise it is the index in [`Self::frame_index`].
161    pub fn atom_index(&self) -> Option<usize> {
162        match self {
163            Self::Duplicate { atom_index, .. }
164            | Self::Extra { atom_index, .. }
165            | Self::Symbol { atom_index, .. }
166            | Self::Mass { atom_index, .. }
167            | Self::Fixed { atom_index, .. }
168            | Self::NonFiniteCoordinate { atom_index, .. } => Some(*atom_index),
169            Self::Missing {
170                reference_index, ..
171            } => Some(*reference_index),
172            _ => None,
173        }
174    }
175}
176
177impl std::fmt::Display for CorrespondenceError {
178    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
179        match self {
180            Self::EmptySequence => {
181                write!(f, "atom-id consistency check needs at least one frame")
182            }
183            Self::Duplicate {
184                frame_index,
185                atom_id,
186                atom_index,
187                first_index,
188            } => write!(
189                f,
190                "duplicate atom_id {atom_id} in frame {frame_index} at index {atom_index} (first index {first_index})"
191            ),
192            Self::Missing {
193                frame_index,
194                atom_id,
195                reference_index,
196            } => write!(
197                f,
198                "frame {frame_index} is missing atom_id {atom_id} (reference index {reference_index})"
199            ),
200            Self::Extra {
201                frame_index,
202                atom_id,
203                atom_index,
204            } => write!(
205                f,
206                "frame {frame_index} atom_id {atom_id} at index {atom_index} is not in the reference frame"
207            ),
208            Self::Symbol {
209                frame_index,
210                atom_id,
211                atom_index,
212                reference_symbol,
213                found_symbol,
214            } => write!(
215                f,
216                "frame {frame_index} atom_id {atom_id} at index {atom_index} has symbol {found_symbol}; reference symbol is {reference_symbol}"
217            ),
218            Self::Mass {
219                frame_index,
220                atom_id,
221                atom_index,
222                reference_mass,
223                found_mass,
224            } => write!(
225                f,
226                "frame {frame_index} atom_id {atom_id} at index {atom_index} has mass {found_mass}; reference mass is {reference_mass} (tolerance {})",
227                ConFrameBuilder::TYPE_MASS_ABS_TOL
228            ),
229            Self::Fixed {
230                frame_index,
231                atom_id,
232                atom_index,
233                reference_fixed,
234                found_fixed,
235            } => write!(
236                f,
237                "frame {frame_index} atom_id {atom_id} at index {atom_index} has fixed mask {found_fixed:?}; reference fixed mask is {reference_fixed:?}"
238            ),
239            Self::SingularCell { determinant } => write!(
240                f,
241                "minimum-image displacement: reference cell matrix is singular (determinant {determinant})"
242            ),
243            Self::NonFiniteDelta => {
244                write!(
245                    f,
246                    "minimum-image displacement: cartesian delta is non-finite"
247                )
248            }
249            Self::NonFiniteCoordinate {
250                frame_index,
251                atom_id,
252                atom_index,
253            } => write!(
254                f,
255                "frame {frame_index} atom_id {atom_id} at index {atom_index} has a non-finite coordinate"
256            ),
257            Self::SectionLength {
258                section,
259                rows,
260                atoms,
261            } => write!(f, "{section} has {rows} rows for {atoms} atoms"),
262            Self::BondEndpoint { index, atoms } => {
263                write!(
264                    f,
265                    "bond endpoint {index} is outside the atom list ({atoms} atoms)"
266                )
267            }
268        }
269    }
270}
271
272impl std::error::Error for CorrespondenceError {}
273
274/// `perm[i]` is the index in `other` of the atom whose `atom_id` equals
275/// `reference.atom_data[i].atom_id`.
276///
277/// Frame index 0 in an error is `reference`. Frame index 1 is `other`.
278pub fn atom_id_permutation(
279    reference: &ConFrame,
280    other: &ConFrame,
281    check: IdentityCheck,
282) -> Result<Vec<usize>, CorrespondenceError> {
283    atom_id_permutation_at(reference, other, 1, check)
284}
285
286fn atom_id_permutation_at(
287    reference: &ConFrame,
288    other: &ConFrame,
289    other_index: usize,
290    check: IdentityCheck,
291) -> Result<Vec<usize>, CorrespondenceError> {
292    index_ids(&reference.atom_data, 0)?;
293    let other_index_map = index_ids(&other.atom_data, other_index)?;
294    let ref_mass = masses_of(reference)?;
295    let other_mass = masses_of(other)?;
296    let mut perm = Vec::with_capacity(reference.atom_data.len());
297    for (ri, atom) in reference.atom_data.iter().enumerate() {
298        let Some(&bi) = other_index_map.get(&atom.atom_id) else {
299            return Err(CorrespondenceError::Missing {
300                frame_index: other_index,
301                atom_id: atom.atom_id,
302                reference_index: ri,
303            });
304        };
305        let found = &other.atom_data[bi];
306        if found.symbol != atom.symbol {
307            return Err(CorrespondenceError::Symbol {
308                frame_index: other_index,
309                atom_id: atom.atom_id,
310                atom_index: bi,
311                reference_symbol: atom.symbol.to_string(),
312                found_symbol: found.symbol.to_string(),
313            });
314        }
315        if check == IdentityCheck::Strict {
316            let reference_mass = ref_mass[ri];
317            let found_mass = other_mass[bi];
318            if !reference_mass.is_finite()
319                || !found_mass.is_finite()
320                || (reference_mass - found_mass).abs() > ConFrameBuilder::TYPE_MASS_ABS_TOL
321            {
322                return Err(CorrespondenceError::Mass {
323                    frame_index: other_index,
324                    atom_id: atom.atom_id,
325                    atom_index: bi,
326                    reference_mass,
327                    found_mass,
328                });
329            }
330            if found.fixed() != atom.fixed() {
331                return Err(CorrespondenceError::Fixed {
332                    frame_index: other_index,
333                    atom_id: atom.atom_id,
334                    atom_index: bi,
335                    reference_fixed: atom.fixed(),
336                    found_fixed: found.fixed(),
337                });
338            }
339        }
340        perm.push(bi);
341    }
342    if other.atom_data.len() != reference.atom_data.len() {
343        let mut seen = vec![false; other.atom_data.len()];
344        for &bi in &perm {
345            seen[bi] = true;
346        }
347        for (i, atom) in other.atom_data.iter().enumerate() {
348            if !seen[i] {
349                return Err(CorrespondenceError::Extra {
350                    frame_index: other_index,
351                    atom_id: atom.atom_id,
352                    atom_index: i,
353                });
354            }
355        }
356    }
357    Ok(perm)
358}
359
360fn index_ids(
361    atoms: &[crate::types::AtomDatum],
362    frame_index: usize,
363) -> Result<FxHashMap<u64, usize>, CorrespondenceError> {
364    let mut map = FxHashMap::with_capacity_and_hasher(atoms.len(), Default::default());
365    for (i, atom) in atoms.iter().enumerate() {
366        if let Some(first) = map.insert(atom.atom_id, i) {
367            return Err(CorrespondenceError::Duplicate {
368                frame_index,
369                atom_id: atom.atom_id,
370                atom_index: i,
371                first_index: first,
372            });
373        }
374    }
375    Ok(map)
376}
377
378fn masses_of(frame: &ConFrame) -> Result<Vec<f64>, CorrespondenceError> {
379    let n = frame.atom_data.len();
380    if frame.masses.len() != n {
381        return Err(CorrespondenceError::SectionLength {
382            section: "masses",
383            rows: frame.masses.len(),
384            atoms: n,
385        });
386    }
387    Ok((0..n).map(|i| frame.masses.get_f64(i)).collect())
388}
389
390/// `other` with atoms stored in `reference` order.
391///
392/// Positions, velocities, forces, and the other per-atom sections move
393/// with the atom. Bond endpoints, which are `atom_data` indices, are
394/// remapped. The cell, the metadata, and the per-frame scalars stay
395/// those of `other`. Type counts are rebuilt from the new order so the
396/// frame still writes as CON when that order is type-grouped.
397pub fn reorder_like(
398    reference: &ConFrame,
399    other: &ConFrame,
400    check: IdentityCheck,
401) -> Result<ConFrame, CorrespondenceError> {
402    let perm = atom_id_permutation(reference, other, check)?;
403    apply_permutation(other, &perm)
404}
405
406fn apply_permutation(frame: &ConFrame, perm: &[usize]) -> Result<ConFrame, CorrespondenceError> {
407    let n = perm.len();
408    let mut out = frame.clone();
409    let mut atoms = Vec::with_capacity(n);
410    for &j in perm {
411        let mut atom = frame.atom_data[j].clone();
412        let p = position_of(frame, j)?;
413        atom.set_x(p[0]);
414        atom.set_y(p[1]);
415        atom.set_z(p[2]);
416        overlay_sections(&mut atom, frame, j)?;
417        atoms.push(atom);
418    }
419    out.atom_data = atoms;
420    out.positions = permute_rows(&frame.positions, perm, "positions")?;
421    out.velocities = permute_rows(&frame.velocities, perm, "velocities")?.into();
422    out.forces = permute_rows(&frame.forces, perm, "forces")?.into();
423    out.magmoms = permute_rows(&frame.magmoms, perm, "magmoms")?.into();
424    out.displacements = permute_rows(&frame.displacements, perm, "displacements")?.into();
425    out.spreads = permute_rows(&frame.spreads, perm, "spreads")?.into();
426    out.atom_energies = permute_col(&frame.atom_energies, perm, "energies")?.into();
427    out.charges = permute_col(&frame.charges, perm, "charges")?.into();
428    out.spins = permute_col(&frame.spins, perm, "spins")?.into();
429    out.masses = permute_col(&frame.masses, perm, "masses")?;
430    let ids: Vec<u64> = perm.iter().map(|&j| frame.atom_data[j].atom_id).collect();
431    out.atom_ids = ndarray::Array1::from_vec(ids).into_shared();
432    // Positions were copied above. When the source stored coordinates only
433    // on the atoms, the permuted array is empty and the atoms carry xyz.
434    if frame.positions.nrows() == n {
435        for (i, atom) in out.atom_data.iter().enumerate() {
436            out.positions.set_f64_row(i, [atom.x(), atom.y(), atom.z()]);
437        }
438    }
439    refresh_type_runs(&mut out);
440    remap_bonds(&mut out, perm)?;
441    Ok(out)
442}
443
444fn position_of(frame: &ConFrame, i: usize) -> Result<[f64; 3], CorrespondenceError> {
445    let n = frame.atom_data.len();
446    if frame.positions.nrows() == n {
447        return Ok(frame.positions.as_f64_row(i));
448    }
449    if frame.positions.nrows() == 0 {
450        let atom = &frame.atom_data[i];
451        return Ok([atom.x(), atom.y(), atom.z()]);
452    }
453    Err(CorrespondenceError::SectionLength {
454        section: "positions",
455        rows: frame.positions.nrows(),
456        atoms: n,
457    })
458}
459
460fn overlay_sections(
461    atom: &mut crate::types::AtomDatum,
462    frame: &ConFrame,
463    j: usize,
464) -> Result<(), CorrespondenceError> {
465    let n = frame.atom_data.len();
466    if frame.velocities.nrows() == n {
467        atom.set_velocity(Some(frame.velocities.as_f64_row(j)));
468    } else if frame.velocities.nrows() != 0 {
469        return Err(CorrespondenceError::SectionLength {
470            section: "velocities",
471            rows: frame.velocities.nrows(),
472            atoms: n,
473        });
474    }
475    if frame.forces.nrows() == n {
476        atom.set_force(Some(frame.forces.as_f64_row(j)));
477    } else if frame.forces.nrows() != 0 {
478        return Err(CorrespondenceError::SectionLength {
479            section: "forces",
480            rows: frame.forces.nrows(),
481            atoms: n,
482        });
483    }
484    if frame.atom_energies.len() == n {
485        atom.set_energy(Some(frame.atom_energies.get_f64(j)));
486    } else if !frame.atom_energies.is_empty() {
487        return Err(CorrespondenceError::SectionLength {
488            section: "energies",
489            rows: frame.atom_energies.len(),
490            atoms: n,
491        });
492    }
493    if frame.charges.len() == n {
494        atom.set_charge(Some(frame.charges.get_f64(j)));
495    } else if !frame.charges.is_empty() {
496        return Err(CorrespondenceError::SectionLength {
497            section: "charges",
498            rows: frame.charges.len(),
499            atoms: n,
500        });
501    }
502    if frame.spins.len() == n {
503        atom.set_spin(Some(frame.spins.get_f64(j)));
504    } else if !frame.spins.is_empty() {
505        return Err(CorrespondenceError::SectionLength {
506            section: "spins",
507            rows: frame.spins.len(),
508            atoms: n,
509        });
510    }
511    if frame.magmoms.nrows() == n {
512        atom.set_magmom(Some(frame.magmoms.as_f64_row(j)));
513    } else if frame.magmoms.nrows() != 0 {
514        return Err(CorrespondenceError::SectionLength {
515            section: "magmoms",
516            rows: frame.magmoms.nrows(),
517            atoms: n,
518        });
519    }
520    if frame.displacements.nrows() == n {
521        atom.set_displacement(Some(frame.displacements.as_f64_row(j)));
522    } else if frame.displacements.nrows() != 0 {
523        return Err(CorrespondenceError::SectionLength {
524            section: "displacements",
525            rows: frame.displacements.nrows(),
526            atoms: n,
527        });
528    }
529    if frame.spreads.nrows() == n {
530        atom.set_spread(Some(frame.spreads.as_f64_row(j)));
531    } else if frame.spreads.nrows() != 0 {
532        return Err(CorrespondenceError::SectionLength {
533            section: "spreads",
534            rows: frame.spreads.nrows(),
535            atoms: n,
536        });
537    }
538    Ok(())
539}
540
541fn permute_rows(
542    src: &FloatArray2,
543    perm: &[usize],
544    section: &'static str,
545) -> Result<FloatArray2, CorrespondenceError> {
546    let n = perm.len();
547    if src.nrows() == 0 {
548        return Ok(FloatArray2::zeros(src.kind(), 0, 3));
549    }
550    if src.nrows() != n {
551        return Err(CorrespondenceError::SectionLength {
552            section,
553            rows: src.nrows(),
554            atoms: n,
555        });
556    }
557    let mut out = FloatArray2::zeros(src.kind(), n, src.ncols().max(3));
558    for (i, &j) in perm.iter().enumerate() {
559        out.set_f64_row(i, src.as_f64_row(j));
560    }
561    Ok(out)
562}
563
564fn permute_col(
565    src: &FloatArray1,
566    perm: &[usize],
567    section: &'static str,
568) -> Result<FloatArray1, CorrespondenceError> {
569    let n = perm.len();
570    if src.is_empty() {
571        return Ok(FloatArray1::zeros(src.kind(), 0));
572    }
573    if src.len() != n {
574        return Err(CorrespondenceError::SectionLength {
575            section,
576            rows: src.len(),
577            atoms: n,
578        });
579    }
580    let mut out = FloatArray1::zeros(src.kind(), n);
581    for (i, &j) in perm.iter().enumerate() {
582        out.set_f64(i, src.get_f64(j));
583    }
584    Ok(out)
585}
586
587fn refresh_type_runs(frame: &mut ConFrame) {
588    let n = frame.atom_data.len();
589    if n == 0 {
590        frame.header.natm_types = 0;
591        frame.header.natms_per_type.clear();
592        frame.header.masses_per_type.clear();
593        return;
594    }
595    let mut counts = Vec::new();
596    let mut masses = Vec::new();
597    let mut run_symbol = frame.atom_data[0].symbol.clone();
598    let mut run_count = 0usize;
599    let mut run_mass = frame.masses.get_f64(0);
600    for i in 0..n {
601        if frame.atom_data[i].symbol != run_symbol {
602            counts.push(run_count);
603            masses.push(run_mass);
604            run_symbol = frame.atom_data[i].symbol.clone();
605            run_count = 0;
606            run_mass = frame.masses.get_f64(i);
607        }
608        run_count += 1;
609    }
610    counts.push(run_count);
611    masses.push(run_mass);
612    frame.header.natm_types = counts.len();
613    frame.header.natms_per_type = counts;
614    frame.header.masses_per_type = masses;
615}
616
617fn remap_bonds(frame: &mut ConFrame, perm: &[usize]) -> Result<(), CorrespondenceError> {
618    let bonds = frame.bonds();
619    if bonds.is_empty() {
620        return Ok(());
621    }
622    let n = perm.len();
623    let mut inverse = vec![0usize; n];
624    for (i, &j) in perm.iter().enumerate() {
625        inverse[j] = i;
626    }
627    let mut mapped = Vec::with_capacity(bonds.len());
628    for bond in bonds {
629        if bond.i as usize >= n {
630            return Err(CorrespondenceError::BondEndpoint {
631                index: bond.i,
632                atoms: n,
633            });
634        }
635        if bond.j as usize >= n {
636            return Err(CorrespondenceError::BondEndpoint {
637                index: bond.j,
638                atoms: n,
639            });
640        }
641        mapped.push(Bond {
642            i: inverse[bond.i as usize] as u32,
643            j: inverse[bond.j as usize] as u32,
644            order: bond.order,
645        });
646    }
647    frame.header.set_bonds(&mapped);
648    Ok(())
649}
650
651/// Walk `frames` in order. Frame 0 is the reference.
652///
653/// The error is the first frame that disagrees with frame 0, and the
654/// atom involved. A one-frame sequence checks that its ids are unique.
655pub fn check_atom_id_consistency(
656    frames: &[ConFrame],
657    check: IdentityCheck,
658) -> Result<(), CorrespondenceError> {
659    let refs: Vec<&ConFrame> = frames.iter().collect();
660    check_atom_id_consistency_refs(&refs, check)
661}
662
663pub(crate) fn check_atom_id_consistency_refs(
664    frames: &[&ConFrame],
665    check: IdentityCheck,
666) -> Result<(), CorrespondenceError> {
667    let Some(reference) = frames.first() else {
668        return Err(CorrespondenceError::EmptySequence);
669    };
670    if frames.len() == 1 {
671        index_ids(&reference.atom_data, 0)?;
672        return Ok(());
673    }
674    for (i, frame) in frames.iter().enumerate().skip(1) {
675        atom_id_permutation_at(reference, frame, i, check)?;
676    }
677    Ok(())
678}
679
680/// Minimum-image displacement of each reference atom toward the atom in
681/// `other` with the same `atom_id`.
682///
683/// Row `i` is the wrapped `r_other - r_reference` for reference atom `i`.
684/// The cell and the periodic flags are the reference frame's
685/// [`FrameHeader::cell_matrix`] and [`FrameHeader::pbc_effective`].
686pub fn minimum_image_displacements(
687    reference: &ConFrame,
688    other: &ConFrame,
689    check: IdentityCheck,
690) -> Result<Vec<[f64; 3]>, CorrespondenceError> {
691    let perm = atom_id_permutation(reference, other, check)?;
692    let cell = reference.header.cell_matrix();
693    let pbc = reference.header.pbc_effective();
694    let mut out = Vec::with_capacity(perm.len());
695    for (i, &j) in perm.iter().enumerate() {
696        let ra = position_of(reference, i)?;
697        let rb = position_of(other, j)?;
698        if !finite3(ra) {
699            return Err(CorrespondenceError::NonFiniteCoordinate {
700                frame_index: 0,
701                atom_id: reference.atom_data[i].atom_id,
702                atom_index: i,
703            });
704        }
705        if !finite3(rb) {
706            return Err(CorrespondenceError::NonFiniteCoordinate {
707                frame_index: 1,
708                atom_id: other.atom_data[j].atom_id,
709                atom_index: j,
710            });
711        }
712        let delta = [rb[0] - ra[0], rb[1] - ra[1], rb[2] - ra[2]];
713        out.push(minimum_image_delta(cell, pbc, delta)?);
714    }
715    Ok(out)
716}
717
718fn finite3(v: [f64; 3]) -> bool {
719    v.iter().all(|x| x.is_finite())
720}
721
722/// Wrap one cartesian delta into the reference cell.
723///
724/// `cell` is a [`FrameHeader::cell_matrix`]: each row is a lattice vector.
725/// Fractional coordinates are `delta * cell^{-1}` (row vector on the left).
726/// A periodic axis subtracts the nearest integer. `f64::round` breaks a
727/// half-integer toward ±infinity, away from zero.
728pub fn minimum_image_delta(
729    cell: [[f64; 3]; 3],
730    pbc: [bool; 3],
731    delta: [f64; 3],
732) -> Result<[f64; 3], CorrespondenceError> {
733    if !finite3(delta) {
734        return Err(CorrespondenceError::NonFiniteDelta);
735    }
736    let (inv, _) = invert_cell(cell).ok_or(CorrespondenceError::SingularCell {
737        determinant: det3(cell),
738    })?;
739    let mut frac = mul_row(delta, inv);
740    for axis in 0..3 {
741        if pbc[axis] {
742            frac[axis] -= frac[axis].round();
743        }
744    }
745    Ok(mul_row(frac, cell))
746}
747
748fn det3(m: [[f64; 3]; 3]) -> f64 {
749    m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
750        - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
751        + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])
752}
753
754fn invert_cell(m: [[f64; 3]; 3]) -> Option<([[f64; 3]; 3], f64)> {
755    let det = det3(m);
756    if !det.is_finite() || det.abs() < 1e-18 {
757        return None;
758    }
759    let inv_det = 1.0 / det;
760    let c00 = m[1][1] * m[2][2] - m[1][2] * m[2][1];
761    let c01 = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
762    let c02 = m[1][0] * m[2][1] - m[1][1] * m[2][0];
763    let c10 = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
764    let c11 = m[0][0] * m[2][2] - m[0][2] * m[2][0];
765    let c12 = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
766    let c20 = m[0][1] * m[1][2] - m[0][2] * m[1][1];
767    let c21 = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
768    let c22 = m[0][0] * m[1][1] - m[0][1] * m[1][0];
769    let inv = [
770        [c00 * inv_det, c10 * inv_det, c20 * inv_det],
771        [c01 * inv_det, c11 * inv_det, c21 * inv_det],
772        [c02 * inv_det, c12 * inv_det, c22 * inv_det],
773    ];
774    Some((inv, det))
775}
776
777/// Row vector times a row-major matrix.
778fn mul_row(v: [f64; 3], m: [[f64; 3]; 3]) -> [f64; 3] {
779    [
780        v[0] * m[0][0] + v[1] * m[1][0] + v[2] * m[2][0],
781        v[0] * m[0][1] + v[1] * m[1][1] + v[2] * m[2][1],
782        v[0] * m[0][2] + v[1] * m[1][2] + v[2] * m[2][2],
783    ]
784}
785
786#[cfg(test)]
787mod tests {
788    use super::*;
789    use crate::types::cell_vectors;
790
791    fn close(a: [f64; 3], b: [f64; 3]) {
792        for k in 0..3 {
793            assert!(
794                (a[k] - b[k]).abs() < 1e-9,
795                "component {k}: {} vs {}",
796                a[k],
797                b[k]
798            );
799        }
800    }
801
802    #[test]
803    fn orthogonal_wrap_uses_the_cell_matrix() {
804        let cell = cell_vectors([10.0, 12.0, 14.0], [90.0, 90.0, 90.0]);
805        let wrapped = minimum_image_delta(cell, [true, true, true], [9.0, 0.0, 0.0]).unwrap();
806        close(wrapped, [-1.0, 0.0, 0.0]);
807        let open = minimum_image_delta(cell, [false, true, true], [9.0, 9.0, 0.0]).unwrap();
808        close(open, [9.0, -3.0, 0.0]);
809    }
810
811    #[test]
812    fn a_lattice_vector_wraps_to_zero() {
813        let cell = cell_vectors([10.0, 12.0, 14.0], [80.0, 90.0, 70.0]);
814        let along_a = minimum_image_delta(cell, [true, true, true], cell[0]).unwrap();
815        close(along_a, [0.0, 0.0, 0.0]);
816        let along_b = minimum_image_delta(cell, [true, true, true], cell[1]).unwrap();
817        close(along_b, [0.0, 0.0, 0.0]);
818    }
819
820    #[test]
821    fn a_flat_cell_is_rejected() {
822        let err = minimum_image_delta([[0.0; 3]; 3], [true; 3], [1.0, 0.0, 0.0]).unwrap_err();
823        assert!(matches!(err, CorrespondenceError::SingularCell { .. }));
824        assert!(err.to_string().contains("singular"));
825    }
826}
readcon_core/correspondence.rs

readcon_core/
correspondence.rs
