xtalforge_core/
transforms.rs

1//! Point-set transforms, finite cells and displacement ellipsoids (Task 1.9).
2
3#![allow(clippy::indexing_slicing)] // small fixed-size matrix and vector math
4
5use crate::{Error, Site, Structure, StructureBuilder};
6
7/// A Cartesian plane through a point.
8#[derive(Debug, Clone, Copy, PartialEq)]
9pub struct Plane {
10    /// Unit normal.
11    pub normal: [f64; 3],
12    /// A point on the plane.
13    pub point: [f64; 3],
14}
15
16/// Displacement ellipsoid axes and semi-axis lengths.
17#[derive(Debug, Clone, Copy, PartialEq)]
18pub struct Ellipsoid {
19    /// Three orthonormal axes.
20    pub axes: [[f64; 3]; 3],
21    /// Semi-axis lengths in Angstrom.
22    pub lengths: [f64; 3],
23    /// Probability level used for the lengths.
24    pub probability: f64,
25}
26
27/// Translates Cartesian points.
28#[must_use]
29pub fn translate(points: &[[f64; 3]], delta: [f64; 3]) -> Vec<[f64; 3]> {
30    points
31        .iter()
32        .map(|point| {
33            [
34                point[0] + delta[0],
35                point[1] + delta[1],
36                point[2] + delta[2],
37            ]
38        })
39        .collect()
40}
41
42/// Rotates Cartesian points about an axis through `center`.
43///
44/// # Errors
45///
46/// Returns [`Error::InvalidInput`] for a zero axis or non-finite values.
47pub fn rotate(
48    points: &[[f64; 3]],
49    axis: [f64; 3],
50    angle_degrees: f64,
51    center: [f64; 3],
52) -> Result<Vec<[f64; 3]>, Error> {
53    if !axis.iter().all(|value| value.is_finite()) || !angle_degrees.is_finite() {
54        return Err(Error::InvalidInput(
55            "rotation axis and angle must be finite".to_owned(),
56        ));
57    }
58    let axis = normalize(axis)?;
59    let angle = angle_degrees.to_radians();
60    let (sine, cosine) = angle.sin_cos();
61    Ok(points
62        .iter()
63        .map(|point| {
64            let relative = subtract(*point, center);
65            let rotated = add(
66                add(scale(relative, cosine), scale(cross(axis, relative), sine)),
67                scale(axis, dot(axis, relative) * (1.0 - cosine)),
68            );
69            add(rotated, center)
70        })
71        .collect())
72}
73
74/// Mirrors Cartesian points through a plane.
75#[must_use]
76pub fn mirror(points: &[[f64; 3]], plane: Plane) -> Vec<[f64; 3]> {
77    points
78        .iter()
79        .map(|point| {
80            let distance = dot(subtract(*point, plane.point), plane.normal);
81            subtract(*point, scale(plane.normal, 2.0 * distance))
82        })
83        .collect()
84}
85
86/// Inverts Cartesian points through a centre.
87#[must_use]
88pub fn invert(points: &[[f64; 3]], center: [f64; 3]) -> Vec<[f64; 3]> {
89    points
90        .iter()
91        .map(|point| {
92            [
93                2.0 * center[0] - point[0],
94                2.0 * center[1] - point[1],
95                2.0 * center[2] - point[2],
96            ]
97        })
98        .collect()
99}
100
101/// Rotates points so `from` aligns with `to`, about `center`.
102///
103/// # Errors
104///
105/// Returns [`Error::InvalidInput`] for zero or non-finite vectors.
106pub fn align_vector(
107    points: &[[f64; 3]],
108    from: [f64; 3],
109    to: [f64; 3],
110    center: [f64; 3],
111) -> Result<Vec<[f64; 3]>, Error> {
112    let from = normalize(from)?;
113    let to = normalize(to)?;
114    let axis = cross(from, to);
115    let sine = norm(axis);
116    if sine < 1e-12 {
117        if dot(from, to) > 0.0 {
118            return Ok(points.to_vec());
119        }
120        let perpendicular = if from[0].abs() < 0.9 {
121            normalize(cross(from, [1.0, 0.0, 0.0]))?
122        } else {
123            normalize(cross(from, [0.0, 1.0, 0.0]))?
124        };
125        return rotate(points, perpendicular, 180.0, center);
126    }
127    let angle = sine.atan2(dot(from, to)).to_degrees();
128    rotate(points, axis, angle, center)
129}
130
131/// Fits a plane through at least three non-collinear points.
132///
133/// # Errors
134///
135/// Returns [`Error::InvalidInput`] for fewer than three points or collinear
136/// points.
137pub fn fit_plane(points: &[[f64; 3]]) -> Result<Plane, Error> {
138    if points.len() < 3 {
139        return Err(Error::InvalidInput(
140            "a plane needs at least three points".to_owned(),
141        ));
142    }
143    let centroid = centroid(points);
144    let mut first = [0.0; 3];
145    for point in points {
146        let candidate = subtract(*point, centroid);
147        if norm(candidate) > 1e-9 {
148            first = candidate;
149            break;
150        }
151    }
152    let mut second = [0.0; 3];
153    for point in points {
154        let candidate = subtract(*point, centroid);
155        let normal = cross(first, candidate);
156        if norm(normal) > 1e-9 {
157            second = candidate;
158            break;
159        }
160    }
161    let normal = normalize(cross(first, second))?;
162    Ok(Plane {
163        normal,
164        point: centroid,
165    })
166}
167
168/// Rotates points so one fitted plane aligns with another.
169///
170/// # Errors
171///
172/// Returns the errors of [`fit_plane`] or [`align_vector`].
173pub fn align_planes(
174    points: &[[f64; 3]],
175    from_normal: [f64; 3],
176    to_normal: [f64; 3],
177    center: [f64; 3],
178) -> Result<Vec<[f64; 3]>, Error> {
179    align_vector(points, from_normal, to_normal, center)
180}
181
182/// Converts anisotropic displacement parameters to an ellipsoid.
183///
184/// `uij` is `[U11, U22, U33, U12, U13, U23]` in Angstrom squared. The
185/// probability is one of 0.5, 0.9, 0.95 or 0.99; other values are interpolated
186/// from the chi-square table for three degrees of freedom.
187///
188/// # Errors
189///
190/// Returns [`Error::InvalidInput`] for non-finite or negative-definite tensors.
191pub fn displacement_ellipsoid(uij: [f64; 6], probability: f64) -> Result<Ellipsoid, Error> {
192    if !uij.iter().all(|value| value.is_finite()) || !probability.is_finite() {
193        return Err(Error::InvalidInput(
194            "displacement parameters must be finite".to_owned(),
195        ));
196    }
197    if !(0.0 < probability && probability < 1.0) {
198        return Err(Error::InvalidInput(
199            "ellipsoid probability must be in (0, 1)".to_owned(),
200        ));
201    }
202    let matrix = [
203        [uij[0], uij[3], uij[4]],
204        [uij[3], uij[1], uij[5]],
205        [uij[4], uij[5], uij[2]],
206    ];
207    let (values, vectors) = jacobi_eigen(matrix)?;
208    let mut lengths = [0.0; 3];
209    for index in 0..3 {
210        if values[index] < -1e-12 {
211            return Err(Error::InvalidInput(
212                "displacement tensor is not positive semidefinite".to_owned(),
213            ));
214        }
215        lengths[index] = values[index].max(0.0).sqrt() * chi2_scale(probability);
216    }
217    Ok(Ellipsoid {
218        axes: vectors,
219        lengths,
220        probability,
221    })
222}
223
224impl Structure {
225    /// Atoms of one finite cell, including images on faces, edges and corners.
226    ///
227    /// # Errors
228    ///
229    /// Returns [`Error::InvalidInput`] for a negative tolerance.
230    pub fn finite_cell(&self, tolerance: f64) -> Result<Self, Error> {
231        if !tolerance.is_finite() || tolerance < 0.0 {
232            return Err(Error::InvalidInput(
233                "finite-cell tolerance must be non-negative and finite".to_owned(),
234            ));
235        }
236        let p1 = self.p1()?;
237        if p1.atoms.is_empty() {
238            return Ok(self.clone());
239        }
240        let lattice = p1.lattice;
241        let mut sites: Vec<Site> = Vec::new();
242        for atom in &p1.atoms {
243            let shifts = atom.frac.map(|value| {
244                let mut out = vec![0];
245                if value <= tolerance {
246                    out.push(1);
247                }
248                if value >= 1.0 - tolerance {
249                    out.push(-1);
250                }
251                out
252            });
253            for x in &shifts[0] {
254                for y in &shifts[1] {
255                    for z in &shifts[2] {
256                        let frac = [
257                            atom.frac[0] + f64::from(*x),
258                            atom.frac[1] + f64::from(*y),
259                            atom.frac[2] + f64::from(*z),
260                        ];
261                        if frac
262                            .iter()
263                            .all(|value| (-tolerance..=1.0 + tolerance).contains(value))
264                        {
265                            let duplicate = sites.iter().any(|existing| {
266                                existing.element() == atom.element
267                                    && existing.frac().iter().zip(frac.iter()).all(
268                                        |(first, second)| {
269                                            (first - second).abs() <= tolerance.max(1e-9)
270                                        },
271                                    )
272                            });
273                            if !duplicate {
274                                sites.push(Site::new(
275                                    format!("{}{}", atom.element, sites.len() + 1),
276                                    &atom.element,
277                                    frac,
278                                )?);
279                            }
280                        }
281                    }
282                }
283            }
284        }
285        let mut builder = StructureBuilder::new(self.name(), lattice);
286        for site in sites {
287            builder = builder.site(site);
288        }
289        builder.build()
290    }
291}
292
293fn centroid(points: &[[f64; 3]]) -> [f64; 3] {
294    let mut sum = [0.0; 3];
295    for point in points {
296        sum[0] += point[0];
297        sum[1] += point[1];
298        sum[2] += point[2];
299    }
300    let count = points.len() as f64;
301    [sum[0] / count, sum[1] / count, sum[2] / count]
302}
303
304fn chi2_scale(probability: f64) -> f64 {
305    let table = [(0.5, 1.538), (0.9, 2.500), (0.95, 2.796), (0.99, 3.368)];
306    if probability <= table[0].0 {
307        return table[0].1;
308    }
309    let mut previous = table[0];
310    for (level, value) in table {
311        if probability <= level {
312            let fraction = (probability - previous.0) / (level - previous.0);
313            return previous.1 + fraction * (value - previous.1);
314        }
315        previous = (level, value);
316    }
317    table[table.len() - 1].1
318}
319
320fn jacobi_eigen(matrix: [[f64; 3]; 3]) -> Result<([f64; 3], [[f64; 3]; 3]), Error> {
321    let mut a = matrix;
322    let mut v = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
323    for _ in 0..32 {
324        let (p, q) = largest_off_diagonal(a);
325        if a[p][q].abs() < 1e-12 {
326            break;
327        }
328        let angle = 0.5 * (2.0 * a[p][q]).atan2(a[q][q] - a[p][p]);
329        let (sine, cosine) = angle.sin_cos();
330        let app = a[p][p];
331        let aqq = a[q][q];
332        let apq = a[p][q];
333        a[p][p] = cosine * cosine * app - 2.0 * sine * cosine * apq + sine * sine * aqq;
334        a[q][q] = sine * sine * app + 2.0 * sine * cosine * apq + cosine * cosine * aqq;
335        a[p][q] = 0.0;
336        a[q][p] = 0.0;
337        for index in 0..3 {
338            if index != p && index != q {
339                let aip = a[index][p];
340                let aiq = a[index][q];
341                a[index][p] = cosine * aip - sine * aiq;
342                a[p][index] = a[index][p];
343                a[index][q] = sine * aip + cosine * aiq;
344                a[q][index] = a[index][q];
345            }
346            let vpi = v[index][p];
347            let vqi = v[index][q];
348            v[index][p] = cosine * vpi - sine * vqi;
349            v[index][q] = sine * vpi + cosine * vqi;
350        }
351    }
352    let mut values = [a[0][0], a[1][1], a[2][2]];
353    let mut axes = [
354        [v[0][0], v[1][0], v[2][0]],
355        [v[0][1], v[1][1], v[2][1]],
356        [v[0][2], v[1][2], v[2][2]],
357    ];
358    for index in 0..3 {
359        if values[index] < -1e-9 {
360            return Err(Error::InvalidInput(
361                "displacement tensor is not positive semidefinite".to_owned(),
362            ));
363        }
364        values[index] = values[index].max(0.0);
365        axes[index] = normalize(axes[index]).unwrap_or([0.0, 0.0, 0.0]);
366    }
367    let mut order = [0usize, 1, 2];
368    order.sort_by(|first, second| {
369        values[*first]
370            .partial_cmp(&values[*second])
371            .unwrap_or(std::cmp::Ordering::Equal)
372    });
373    Ok((
374        [values[order[0]], values[order[1]], values[order[2]]],
375        [axes[order[0]], axes[order[1]], axes[order[2]]],
376    ))
377}
378
379fn largest_off_diagonal(matrix: [[f64; 3]; 3]) -> (usize, usize) {
380    let candidates = [(0, 1), (0, 2), (1, 2)];
381    let mut best = candidates[0];
382    for candidate in candidates {
383        if matrix[candidate.0][candidate.1].abs() > matrix[best.0][best.1].abs() {
384            best = candidate;
385        }
386    }
387    best
388}
389
390fn add(first: [f64; 3], second: [f64; 3]) -> [f64; 3] {
391    [
392        first[0] + second[0],
393        first[1] + second[1],
394        first[2] + second[2],
395    ]
396}
397
398fn subtract(first: [f64; 3], second: [f64; 3]) -> [f64; 3] {
399    [
400        first[0] - second[0],
401        first[1] - second[1],
402        first[2] - second[2],
403    ]
404}
405
406fn scale(vector: [f64; 3], factor: f64) -> [f64; 3] {
407    [vector[0] * factor, vector[1] * factor, vector[2] * factor]
408}
409
410fn cross(first: [f64; 3], second: [f64; 3]) -> [f64; 3] {
411    [
412        first[1] * second[2] - first[2] * second[1],
413        first[2] * second[0] - first[0] * second[2],
414        first[0] * second[1] - first[1] * second[0],
415    ]
416}
417
418fn dot(first: [f64; 3], second: [f64; 3]) -> f64 {
419    first[0] * second[0] + first[1] * second[1] + first[2] * second[2]
420}
421
422fn norm(vector: [f64; 3]) -> f64 {
423    dot(vector, vector).sqrt()
424}
425
426fn normalize(vector: [f64; 3]) -> Result<[f64; 3], Error> {
427    let length = norm(vector);
428    if length < 1e-12 {
429        return Err(Error::InvalidInput("the vector has zero length".to_owned()));
430    }
431    Ok(scale(vector, 1.0 / length))
432}
xtalforge_core/transforms.rs

xtalforge_core/
transforms.rs
