xtalforge_core/
surface.rs

1#![allow(clippy::indexing_slicing, clippy::needless_range_loop)] // fixed-size numeric kernels
2
3//! Crystal surfaces: Wulff construction and adsorption sites (Task 9.4).
4//!
5//! The Wulff shape is the intersection of the half-spaces whose normals are
6//! the reciprocal-lattice directions and whose offsets are the surface
7//! energies. Adsorption sites are the top, bridge and hollow positions of the
8//! uppermost layer of a slab.
9
10use serde::{Deserialize, Serialize};
11
12use crate::polyhedron::{HalfSpace, Polyhedron};
13use crate::{Error, Lattice, Site, Structure, StructureBuilder};
14
15/// Surface energy of one Miller direction.
16#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
17pub struct SurfaceEnergy {
18    /// Miller index.
19    pub hkl: [i32; 3],
20    /// Surface energy in J/m2; only ratios matter for the shape.
21    pub energy: f64,
22}
23
24/// One facet of the Wulff shape.
25#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
26pub struct WulffFacet {
27    /// Miller index of the family that produced the facet.
28    pub hkl: [i32; 3],
29    /// Outward unit normal in Cartesian coordinates.
30    pub normal: [f64; 3],
31    /// Facet area in square Angstrom.
32    pub area: f64,
33    /// Fraction of the total surface area.
34    pub fraction: f64,
35    /// Surface energy of the facet.
36    pub energy: f64,
37}
38
39/// The Wulff polyhedron.
40#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
41pub struct WulffShape {
42    /// Vertices in Cartesian coordinates.
43    pub vertices: Vec<[f64; 3]>,
44    /// Facets with their areas and fractions.
45    pub facets: Vec<WulffFacet>,
46    /// Volume in cubic Angstrom.
47    pub volume: f64,
48    /// Total surface area in square Angstrom.
49    pub surface_area: f64,
50}
51
52/// Builds the Wulff shape from a lattice and a list of surface energies.
53///
54/// Directions are reciprocal-lattice vectors of `hkl`; include every symmetry
55/// equivalent direction (the CLI expands them with the point group).
56///
57/// # Errors
58///
59/// Returns [`Error::InvalidInput`] for an empty list, a non-positive energy or
60/// a degenerate intersection.
61pub fn wulff(lattice: &Lattice, energies: &[SurfaceEnergy]) -> Result<WulffShape, Error> {
62    if energies.is_empty() {
63        return Err(Error::InvalidInput(
64            "the Wulff construction needs at least one surface energy".to_owned(),
65        ));
66    }
67    let reciprocal = lattice.reciprocal(false)?;
68    let matrix = reciprocal.matrix();
69    let mut planes = Vec::with_capacity(energies.len());
70    let mut normals = Vec::with_capacity(energies.len());
71    for surface in energies {
72        if !surface.energy.is_finite() || surface.energy <= 0.0 {
73            return Err(Error::InvalidInput(format!(
74                "surface energy for {:?} must be positive and finite",
75                surface.hkl
76            )));
77        }
78        let [h, k, l] = surface.hkl.map(f64::from);
79        let normal = [
80            h * matrix[0][0] + k * matrix[1][0] + l * matrix[2][0],
81            h * matrix[0][1] + k * matrix[1][1] + l * matrix[2][1],
82            h * matrix[0][2] + k * matrix[1][2] + l * matrix[2][2],
83        ];
84        let length = (normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]).sqrt();
85        if length < 1e-12 {
86            return Err(Error::InvalidInput(format!(
87                "direction {:?} has no reciprocal vector",
88                surface.hkl
89            )));
90        }
91        planes.push(HalfSpace::new(normal, surface.energy)?);
92        normals.push((
93            [normal[0] / length, normal[1] / length, normal[2] / length],
94            surface.hkl,
95            surface.energy,
96        ));
97    }
98    let polyhedron = Polyhedron::from_half_spaces(&planes)?;
99    let total = polyhedron.surface_area();
100    if total <= 0.0 {
101        return Err(Error::InvalidInput(
102            "the Wulff shape has no surface".to_owned(),
103        ));
104    }
105    let mut facets = Vec::new();
106    for (face, _) in polyhedron.faces.iter().enumerate() {
107        let Some(normal) = polyhedron.face_normal(face) else {
108            continue;
109        };
110        let mut best: Option<(f64, [i32; 3], f64)> = None;
111        for (candidate, hkl, energy) in &normals {
112            let cosine =
113                normal[0] * candidate[0] + normal[1] * candidate[1] + normal[2] * candidate[2];
114            if cosine > best.map_or(0.99, |(value, _, _)| value) {
115                best = Some((cosine, *hkl, *energy));
116            }
117        }
118        let Some((_, hkl, energy)) = best else {
119            continue;
120        };
121        let area = polyhedron.face_area(face);
122        facets.push(WulffFacet {
123            hkl,
124            normal,
125            area,
126            fraction: area / total,
127            energy,
128        });
129    }
130    facets.sort_by(|first, second| {
131        second
132            .area
133            .partial_cmp(&first.area)
134            .unwrap_or(std::cmp::Ordering::Equal)
135            .then_with(|| first.hkl.cmp(&second.hkl))
136    });
137    let volume = polyhedron.volume();
138    Ok(WulffShape {
139        vertices: polyhedron.vertices,
140        facets,
141        volume,
142        surface_area: total,
143    })
144}
145
146/// Type of an adsorption site.
147#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
148#[serde(rename_all = "snake_case")]
149pub enum AdsorptionKind {
150    /// On top of one surface atom.
151    Top,
152    /// Above the midpoint of two bonded surface atoms.
153    Bridge,
154    /// Above the centroid of three or four bonded surface atoms.
155    Hollow,
156}
157
158/// One adsorption site on a slab.
159#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
160pub struct AdsorptionSite {
161    /// Site type.
162    pub kind: AdsorptionKind,
163    /// Cartesian position on the surface.
164    pub position: [f64; 3],
165    /// Fractional coordinates in the slab cell.
166    pub frac: [f64; 3],
167    /// P1 indices of the surface atoms that define the site.
168    pub atoms: Vec<usize>,
169}
170
171/// Enumerates top, bridge and hollow sites of the uppermost slab layer.
172///
173/// # Errors
174///
175/// Returns [`Error::InvalidInput`] for an empty slab or a non-positive
176/// tolerance, or the errors of the P1 expansion.
177pub fn adsorption_sites(slab: &Structure, tolerance: f64) -> Result<Vec<AdsorptionSite>, Error> {
178    if !tolerance.is_finite() || tolerance <= 0.0 {
179        return Err(Error::InvalidInput(
180            "adsorption tolerance must be positive and finite".to_owned(),
181        ));
182    }
183    let p1 = slab.p1()?;
184    if p1.atoms.is_empty() {
185        return Err(Error::InvalidInput("the slab has no atoms".to_owned()));
186    }
187    let top = p1
188        .atoms
189        .iter()
190        .map(|atom| atom.frac[2])
191        .fold(f64::NEG_INFINITY, f64::max);
192    let surface: Vec<usize> = p1
193        .atoms
194        .iter()
195        .enumerate()
196        .filter(|(_, atom)| (atom.frac[2] - top).abs() <= tolerance)
197        .map(|(index, _)| index)
198        .collect();
199    let mut sites: Vec<AdsorptionSite> = Vec::new();
200    for &index in &surface {
201        let Some(atom) = p1.atoms.get(index) else {
202            continue;
203        };
204        let frac = [atom.frac[0], atom.frac[1], top];
205        sites.push(AdsorptionSite {
206            kind: AdsorptionKind::Top,
207            position: p1.lattice.to_cart(frac),
208            frac,
209            atoms: vec![index],
210        });
211    }
212    let mut bonds = Vec::new();
213    for (position, &first) in surface.iter().enumerate() {
214        for &second in surface.iter().skip(position + 1) {
215            let (Some(a), Some(b)) = (p1.atoms.get(first), p1.atoms.get(second)) else {
216                continue;
217            };
218            let distance = p1.lattice.distance(a.frac, b.frac);
219            let radius = crate::elements::element_by_symbol(&a.element)?.covalent_radius
220                + crate::elements::element_by_symbol(&b.element)?.covalent_radius;
221            if distance <= 1.25 * radius {
222                bonds.push((first, second, distance));
223                let frac = [
224                    (a.frac[0] + b.frac[0]) / 2.0,
225                    (a.frac[1] + b.frac[1]) / 2.0,
226                    top,
227                ];
228                sites.push(AdsorptionSite {
229                    kind: AdsorptionKind::Bridge,
230                    position: p1.lattice.to_cart(frac),
231                    frac,
232                    atoms: vec![first, second],
233                });
234            }
235        }
236    }
237    for (position, &first) in surface.iter().enumerate() {
238        for (other, &second) in surface.iter().enumerate().skip(position + 1) {
239            for &third in surface.iter().skip(other + 1) {
240                let has = |a: usize, b: usize| {
241                    bonds.iter().any(|entry| {
242                        (entry.0 == a && entry.1 == b) || (entry.0 == b && entry.1 == a)
243                    })
244                };
245                if !(has(first, second) && has(second, third) && has(first, third)) {
246                    continue;
247                }
248                let (Some(a), Some(b), Some(c)) = (
249                    p1.atoms.get(first),
250                    p1.atoms.get(second),
251                    p1.atoms.get(third),
252                ) else {
253                    continue;
254                };
255                let frac = [
256                    (a.frac[0] + b.frac[0] + c.frac[0]) / 3.0,
257                    (a.frac[1] + b.frac[1] + c.frac[1]) / 3.0,
258                    top,
259                ];
260                sites.push(AdsorptionSite {
261                    kind: AdsorptionKind::Hollow,
262                    position: p1.lattice.to_cart(frac),
263                    frac,
264                    atoms: vec![first, second, third],
265                });
266            }
267        }
268    }
269    deduplicate(&mut sites, tolerance);
270    Ok(sites)
271}
272
273/// Places a molecule above an adsorption site.
274///
275/// The molecule centroid is moved to `site.position` plus `offset` Angstrom
276/// along the slab normal (the `c` axis of a slab cell), and its atoms are
277/// appended in fractional coordinates of the slab.
278///
279/// # Errors
280///
281/// Returns [`Error::InvalidInput`] for an empty molecule or a non-finite
282/// offset, or a builder error.
283pub fn place_molecule(
284    slab: &Structure,
285    molecule: &Structure,
286    site: &AdsorptionSite,
287    offset: f64,
288) -> Result<Structure, Error> {
289    if !offset.is_finite() {
290        return Err(Error::InvalidInput(
291            "the adsorption offset must be finite".to_owned(),
292        ));
293    }
294    let molecule_p1 = molecule.p1()?;
295    if molecule_p1.atoms.is_empty() {
296        return Err(Error::InvalidInput("the molecule has no atoms".to_owned()));
297    }
298    let normal = slab.lattice().plane_normal([0.0, 0.0, 1.0])?;
299    let target = [
300        site.position[0] + normal[0] * offset,
301        site.position[1] + normal[1] * offset,
302        site.position[2] + normal[2] * offset,
303    ];
304    let mut centroid = [0.0; 3];
305    for atom in &molecule_p1.atoms {
306        let cart = molecule_p1.lattice.to_cart(atom.frac);
307        centroid[0] += cart[0];
308        centroid[1] += cart[1];
309        centroid[2] += cart[2];
310    }
311    let count = molecule_p1.atoms.len() as f64;
312    centroid = [
313        centroid[0] / count,
314        centroid[1] / count,
315        centroid[2] / count,
316    ];
317    let mut builder = StructureBuilder::new(slab.name(), *slab.lattice()).hall(slab.hall());
318    for operation in slab.operations() {
319        builder = builder.operation(*operation);
320    }
321    for site in slab.sites() {
322        builder = builder.site(site.clone());
323    }
324    for atom in &molecule_p1.atoms {
325        let cart = molecule_p1.lattice.to_cart(atom.frac);
326        let moved = [
327            cart[0] - centroid[0] + target[0],
328            cart[1] - centroid[1] + target[1],
329            cart[2] - centroid[2] + target[2],
330        ];
331        let frac = slab.lattice().to_frac(moved);
332        builder = builder.site(
333            Site::new(atom.label.clone(), &atom.element, frac)?.with_occupancy(atom.occupancy)?,
334        );
335    }
336    for (key, value) in slab.metadata() {
337        builder = builder.metadata(key.clone(), value.clone());
338    }
339    builder.build()
340}
341
342fn deduplicate(sites: &mut Vec<AdsorptionSite>, tolerance: f64) {
343    let mut kept: Vec<AdsorptionSite> = Vec::new();
344    for site in sites.drain(..) {
345        let duplicate = kept.iter().any(|existing| {
346            existing.kind == site.kind
347                && existing.atoms.iter().all(|atom| site.atoms.contains(atom))
348                && distance(existing.position, site.position) <= tolerance
349        });
350        if !duplicate {
351            kept.push(site);
352        }
353    }
354    kept.sort_by(|first, second| {
355        kind_order(first.kind)
356            .cmp(&kind_order(second.kind))
357            .then_with(|| {
358                first
359                    .position
360                    .partial_cmp(&second.position)
361                    .unwrap_or(std::cmp::Ordering::Equal)
362            })
363    });
364    *sites = kept;
365}
366
367fn kind_order(kind: AdsorptionKind) -> u8 {
368    match kind {
369        AdsorptionKind::Top => 0,
370        AdsorptionKind::Bridge => 1,
371        AdsorptionKind::Hollow => 2,
372    }
373}
374
375fn distance(first: [f64; 3], second: [f64; 3]) -> f64 {
376    let dx = first[0] - second[0];
377    let dy = first[1] - second[1];
378    let dz = first[2] - second[2];
379    (dx * dx + dy * dy + dz * dz).sqrt()
380}
xtalforge_core/surface.rs

xtalforge_core/
surface.rs
