xtalforge_core/
edits.rs

1#![allow(clippy::indexing_slicing, clippy::needless_range_loop)] // fixed-size numeric kernels
2
3//! Structure edits: strain, rattle, site operations and optimal supercells.
4//!
5//! Every operation is additive and deterministic. The random operations take
6//! an explicit seed (AGENTS.md rule 6) and use a small local SplitMix64
7//! generator so the results do not depend on the platform.
8
9use std::collections::BTreeMap;
10
11use crate::{Error, Lattice, Site, Structure, StructureBuilder};
12
13/// Converts engineering Voigt strain `[exx, eyy, ezz, gyz, gxz, gxy]` to a
14/// symmetric tensor.
15///
16/// # Errors
17///
18/// Returns [`Error::InvalidInput`] for non-finite values.
19pub fn voigt_to_tensor(voigt: [f64; 6]) -> Result<[[f64; 3]; 3], Error> {
20    if !voigt.iter().all(|value| value.is_finite()) {
21        return Err(Error::InvalidInput(
22            "strain components must be finite".to_owned(),
23        ));
24    }
25    let [xx, yy, zz, yz, xz, xy] = voigt;
26    Ok([
27        [xx, 0.5 * xy, 0.5 * xz],
28        [0.5 * xy, yy, 0.5 * yz],
29        [0.5 * xz, 0.5 * yz, zz],
30    ])
31}
32
33/// Converts a symmetric strain tensor to engineering Voigt components.
34#[must_use]
35pub fn tensor_to_voigt(tensor: [[f64; 3]; 3]) -> [f64; 6] {
36    [
37        tensor[0][0],
38        tensor[1][1],
39        tensor[2][2],
40        2.0 * tensor[1][2],
41        2.0 * tensor[0][2],
42        2.0 * tensor[0][1],
43    ]
44}
45
46/// Applies a Cartesian strain tensor to the cell; fractional coordinates stay.
47///
48/// # Errors
49///
50/// Returns the errors of [`Lattice::strained`] or a builder error.
51pub fn strain(structure: &Structure, tensor: [[f64; 3]; 3]) -> Result<Structure, Error> {
52    let lattice = structure.lattice().strained(tensor)?;
53    copy_with_lattice(structure, lattice)
54}
55
56/// Applies a deformation gradient `F` (`x' = F x`).
57///
58/// # Errors
59///
60/// Returns [`Error::InvalidInput`] for a non-finite or singular gradient, or a
61/// builder error.
62pub fn deform(structure: &Structure, gradient: [[f64; 3]; 3]) -> Result<Structure, Error> {
63    if !gradient.iter().flatten().all(|value| value.is_finite()) {
64        return Err(Error::InvalidInput(
65            "deformation gradient must be finite".to_owned(),
66        ));
67    }
68    let determinant = determinant(gradient);
69    if determinant.abs() < 1e-9 || determinant <= 0.0 {
70        return Err(Error::InvalidInput(
71            "deformation gradient must have a positive determinant".to_owned(),
72        ));
73    }
74    // `Lattice::strained` computes M (I + eps)^T; with eps = F^T - I the
75    // result is M F, which is the row form of x' = F x.
76    let transpose = [
77        [gradient[0][0], gradient[1][0], gradient[2][0]],
78        [gradient[0][1], gradient[1][1], gradient[2][1]],
79        [gradient[0][2], gradient[1][2], gradient[2][2]],
80    ];
81    let eps = [
82        [transpose[0][0] - 1.0, transpose[0][1], transpose[0][2]],
83        [transpose[1][0], transpose[1][1] - 1.0, transpose[1][2]],
84        [transpose[2][0], transpose[2][1], transpose[2][2] - 1.0],
85    ];
86    strain(structure, eps)
87}
88
89/// Isotropically scales the cell to a target volume in cubic Angstrom.
90///
91/// # Errors
92///
93/// Returns the errors of [`Lattice::scaled_to_volume`] or a builder error.
94pub fn scale_to_volume(structure: &Structure, volume: f64) -> Result<Structure, Error> {
95    let lattice = structure.lattice().scaled_to_volume(volume)?;
96    copy_with_lattice(structure, lattice)
97}
98
99/// Isotropically scales the cell to a target density in g/cm3.
100///
101/// # Errors
102///
103/// Returns [`Error::InvalidInput`] for a non-positive density, or the errors
104/// of [`Structure::cell_mass`] and [`scale_to_volume`].
105pub fn scale_to_density(structure: &Structure, density: f64) -> Result<Structure, Error> {
106    if !density.is_finite() || density <= 0.0 {
107        return Err(Error::InvalidInput(format!(
108            "target density must be positive and finite, got {density}"
109        )));
110    }
111    let mass = structure.cell_mass()?;
112    let volume = mass / (density * 0.602_214_076_9);
113    scale_to_volume(structure, volume)
114}
115
116/// Options for [`rattle`].
117#[derive(Debug, Clone, PartialEq)]
118pub struct RattleOptions {
119    /// Gaussian standard deviation in Angstrom (or the maximum norm when
120    /// `bounded` is set).
121    pub amplitude: f64,
122    /// Explicit seed for the deterministic generator.
123    pub seed: u64,
124    /// Draw a uniform direction with a fixed norm instead of a Gaussian.
125    pub bounded: bool,
126    /// Per-element amplitude overrides in Angstrom.
127    pub per_element: BTreeMap<String, f64>,
128}
129
130impl Default for RattleOptions {
131    fn default() -> Self {
132        Self {
133            amplitude: 0.05,
134            seed: 0,
135            bounded: false,
136            per_element: BTreeMap::new(),
137        }
138    }
139}
140
141/// The rattled structure and its displacement statistics.
142#[derive(Debug, Clone, PartialEq)]
143pub struct RattleReport {
144    /// The displaced P1 structure.
145    pub structure: Structure,
146    /// Largest Cartesian displacement in Angstrom.
147    pub max_displacement: f64,
148    /// Root mean square Cartesian displacement in Angstrom.
149    pub rms_displacement: f64,
150    /// Number of displaced atoms.
151    pub n_atoms: usize,
152}
153
154/// Randomly displaces every atom of the P1 cell.
155///
156/// # Errors
157///
158/// Returns [`Error::InvalidInput`] for a non-finite amplitude, or a builder
159/// error.
160pub fn rattle(structure: &Structure, options: &RattleOptions) -> Result<RattleReport, Error> {
161    if !options.amplitude.is_finite() || options.amplitude < 0.0 {
162        return Err(Error::InvalidInput(
163            "rattle amplitude must be non-negative and finite".to_owned(),
164        ));
165    }
166    for (element, amplitude) in &options.per_element {
167        if !amplitude.is_finite() || *amplitude < 0.0 {
168            return Err(Error::InvalidInput(format!(
169                "rattle amplitude for {element} must be non-negative and finite"
170            )));
171        }
172    }
173    let p1 = structure.p1()?;
174    let mut rng = SplitMix64::new(options.seed);
175    let mut builder = StructureBuilder::new(structure.name(), p1.lattice);
176    let mut squared = 0.0;
177    let mut largest = 0.0_f64;
178    for atom in &p1.atoms {
179        let amplitude = options
180            .per_element
181            .get(&atom.element)
182            .copied()
183            .unwrap_or(options.amplitude);
184        let displacement = if options.bounded {
185            let direction = rng.unit_vector();
186            [
187                direction[0] * amplitude,
188                direction[1] * amplitude,
189                direction[2] * amplitude,
190            ]
191        } else {
192            let scale = amplitude * (rng.gaussian() / std::f64::consts::SQRT_2);
193            let direction = rng.unit_vector();
194            [
195                direction[0] * scale,
196                direction[1] * scale,
197                direction[2] * scale,
198            ]
199        };
200        let cart = p1.lattice.to_cart(atom.frac);
201        let moved = [
202            cart[0] + displacement[0],
203            cart[1] + displacement[1],
204            cart[2] + displacement[2],
205        ];
206        let frac = Lattice::wrap(p1.lattice.to_frac(moved));
207        let site =
208            Site::new(atom.label.clone(), &atom.element, frac)?.with_occupancy(atom.occupancy)?;
209        builder = builder.site(site);
210        let norm = (displacement[0] * displacement[0]
211            + displacement[1] * displacement[1]
212            + displacement[2] * displacement[2])
213            .sqrt();
214        squared += norm * norm;
215        largest = largest.max(norm);
216    }
217    Ok(RattleReport {
218        structure: builder.build()?,
219        max_displacement: largest,
220        rms_displacement: if p1.atoms.is_empty() {
221            0.0
222        } else {
223            (squared / p1.atoms.len() as f64).sqrt()
224        },
225        n_atoms: p1.atoms.len(),
226    })
227}
228
229mod sites;
230mod supercell;
231
232pub use sites::{
233    SiteSort, SpeciesReplacement, add_sites, move_site, remove_sites, replace_site,
234    replace_species, sort_sites, wrap_sites,
235};
236pub use supercell::{SupercellChoice, optimal_supercell};
237
238fn copy_with_lattice(structure: &Structure, lattice: Lattice) -> Result<Structure, Error> {
239    let mut builder = StructureBuilder::new(structure.name(), lattice).hall(structure.hall());
240    for operation in structure.operations() {
241        builder = builder.operation(*operation);
242    }
243    for site in structure.sites() {
244        builder = builder.site(site.clone());
245    }
246    for (key, value) in structure.metadata() {
247        builder = builder.metadata(key.clone(), value.clone());
248    }
249    builder.build()
250}
251
252/// Deterministic SplitMix64 generator.
253#[derive(Debug, Clone)]
254pub(crate) struct SplitMix64 {
255    state: u64,
256}
257
258impl SplitMix64 {
259    pub(crate) fn new(seed: u64) -> Self {
260        Self { state: seed }
261    }
262
263    pub(crate) fn next_u64(&mut self) -> u64 {
264        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
265        let mut z = self.state;
266        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
267        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
268        z ^ (z >> 31)
269    }
270
271    fn uniform(&mut self) -> f64 {
272        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
273    }
274
275    pub(crate) fn gaussian(&mut self) -> f64 {
276        let u1 = self.uniform().max(1e-300);
277        let u2 = self.uniform();
278        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
279    }
280
281    pub(crate) fn unit_vector(&mut self) -> [f64; 3] {
282        let z = 2.0 * self.uniform() - 1.0;
283        let phi = 2.0 * std::f64::consts::PI * self.uniform();
284        let radius = (1.0 - z * z).max(0.0).sqrt();
285        [radius * phi.cos(), radius * phi.sin(), z]
286    }
287
288    pub(crate) fn shuffle<T>(&mut self, values: &mut [T]) {
289        for index in (1..values.len()).rev() {
290            let other = (self.next_u64() % (index as u64 + 1)) as usize;
291            values.swap(index, other);
292        }
293    }
294}
295
296fn determinant(matrix: [[f64; 3]; 3]) -> f64 {
297    let [[a, b, c], [d, e, f], [g, h, i]] = matrix;
298    a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)
299}
xtalforge_core/edits.rs

xtalforge_core/
edits.rs
