Skip to main content

lindhard/
material.rs

1//! Materials: elements combined by atom or mass fraction, with a mass density,
2//! a derived atom number density, and per-element displacement (`E_d`), lattice
3//! binding (`E_b`) and surface binding (`E_s`) energies.
4//!
5//! # Conventions (read this before trusting a default)
6//!
7//! * Fractions passed in are **relative weights**: they are normalised by
8//!   dividing by their sum, so stoichiometric counts (`Si 1, O 2`) are valid.
9//!   The sum must be finite and strictly positive and every entry finite and
10//!   non-negative; otherwise construction fails. Normalisation is exact
11//!   division, with no tolerance window.
12//! * `E_d`, `E_b`, `E_s` are **model parameters, not measured properties of
13//!   your sample**. Defaults come from [`crate::elements`] and are conventions
14//!   with citations there: `E_s` defaults to the elemental cohesive
15//!   (sublimation) energy where one is tabulated; `E_d` defaults to the ASTM
16//!   E521 convention value where one is recorded; `E_b` defaults to 0 eV, the
17//!   common choice for amorphous-target binary-collision runs. Where there is no
18//!   default (`None`) the getter returns [`MaterialError::EnergyNotSet`] so the
19//!   user has to decide. Energies are stored in eV (suffix `_ev`); see
20//!   [`crate::units`].
21//! * Elements keep the order they were given in, so every output derived from a
22//!   material is deterministic.
23//! * Atom number density follows `n_i = x_i * rho * N_A / sum_j(x_j * M_j)`
24//!   with `x` the atom fractions, `rho` the mass density and `M` the atomic
25//!   weights (CODATA 2022 `N_A`).
26
27use serde::{Deserialize, Serialize};
28
29use crate::constants::AVOGADRO;
30use crate::elements::{element, element_by_symbol};
31use crate::units::{g_cm3_to_kg_m3, kg_m3_to_g_cm3};
32
33/// Which per-element energy a message is about.
34#[derive(Debug, Clone, Copy, PartialEq, Eq)]
35pub enum EnergyKind {
36    /// Displacement energy `E_d`.
37    Displacement,
38    /// Lattice binding energy `E_b`.
39    LatticeBinding,
40    /// Surface binding energy `E_s`.
41    SurfaceBinding,
42}
43
44impl std::fmt::Display for EnergyKind {
45    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
46        f.write_str(match self {
47            Self::Displacement => "displacement energy E_d",
48            Self::LatticeBinding => "lattice binding energy E_b",
49            Self::SurfaceBinding => "surface binding energy E_s",
50        })
51    }
52}
53
54/// Errors from building or querying a [`Material`].
55#[derive(Debug, Clone, PartialEq, thiserror::Error)]
56pub enum MaterialError {
57    /// No elements were given.
58    #[error("material has no elements")]
59    Empty,
60    /// Atomic number outside 1..=92.
61    #[error("unknown atomic number {0}")]
62    UnknownAtomicNumber(u8),
63    /// Symbol not in the element table.
64    #[error("unknown element symbol {0:?}")]
65    UnknownSymbol(String),
66    /// The same element appeared twice.
67    #[error("element Z={0} listed more than once")]
68    DuplicateElement(u8),
69    /// A fraction was negative, NaN or infinite.
70    #[error("invalid fraction {value} for element Z={z}: must be finite and non-negative")]
71    InvalidFraction {
72        /// Atomic number.
73        z: u8,
74        /// Offending value.
75        value: f64,
76    },
77    /// Fractions sum to zero or to a non-finite value, so cannot be normalised.
78    #[error("fractions sum to {0}: cannot be normalised")]
79    NotNormalizable(f64),
80    /// A density was zero, negative or not finite.
81    #[error("invalid density {0} kg/m^3: must be finite and positive")]
82    InvalidDensity(f64),
83    /// No density given and none can be inferred (compound, or element with no
84    /// tabulated solid density).
85    #[error("no density given and none can be inferred for this composition")]
86    MissingDensity,
87    /// An energy was negative, NaN or infinite.
88    #[error("invalid {kind} {value} eV for element Z={z}: must be finite and non-negative")]
89    InvalidEnergy {
90        /// Atomic number.
91        z: u8,
92        /// Which energy.
93        kind: EnergyKind,
94        /// Offending value.
95        value: f64,
96    },
97    /// The element is not part of the material.
98    #[error("element Z={0} is not in this material")]
99    ElementNotInMaterial(u8),
100    /// No default exists for this energy and the user has not set one.
101    #[error("{kind} for element Z={z} has no default; set it explicitly")]
102    EnergyNotSet {
103        /// Atomic number.
104        z: u8,
105        /// Which energy.
106        kind: EnergyKind,
107    },
108    /// A serialized element entry was malformed (spec level problems).
109    #[error("invalid material specification: {0}")]
110    InvalidSpec(String),
111}
112
113/// One element within a material.
114#[derive(Debug, Clone, PartialEq)]
115pub struct Component {
116    z: u8,
117    atom_fraction: f64,
118    e_d_ev: Option<f64>,
119    e_b_ev: Option<f64>,
120    e_s_ev: Option<f64>,
121}
122
123impl Component {
124    /// Atomic number.
125    pub fn z(&self) -> u8 {
126        self.z
127    }
128    /// Normalised atom fraction.
129    pub fn atom_fraction(&self) -> f64 {
130        self.atom_fraction
131    }
132}
133
134/// A material: a list of elements, a mass density and per-element energies.
135///
136/// Serializes (serde) to the canonical form of [`MaterialSpec`] with atom
137/// fractions and explicit densities and energies; deserializes from any valid
138/// spec, so a TOML file may give atom or mass fractions.
139#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
140#[serde(try_from = "MaterialSpec", into = "MaterialSpec")]
141pub struct Material {
142    name: Option<String>,
143    components: Vec<Component>,
144    density_kg_m3: f64,
145}
146
147fn check_fraction(z: u8, v: f64) -> Result<(), MaterialError> {
148    if v.is_finite() && v >= 0.0 {
149        Ok(())
150    } else {
151        Err(MaterialError::InvalidFraction { z, value: v })
152    }
153}
154
155fn check_energy(z: u8, kind: EnergyKind, v: f64) -> Result<f64, MaterialError> {
156    if v.is_finite() && v >= 0.0 {
157        Ok(v)
158    } else {
159        Err(MaterialError::InvalidEnergy { z, kind, value: v })
160    }
161}
162
163/// Validate, normalise and (if mass fractions) convert to atom fractions.
164fn normalise(parts: &[(u8, f64)], by_mass: bool) -> Result<Vec<(u8, f64)>, MaterialError> {
165    if parts.is_empty() {
166        return Err(MaterialError::Empty);
167    }
168    let mut out: Vec<(u8, f64)> = Vec::with_capacity(parts.len());
169    for &(z, f) in parts {
170        let el = element(z).ok_or(MaterialError::UnknownAtomicNumber(z))?;
171        check_fraction(z, f)?;
172        if out.iter().any(|&(zz, _)| zz == z) {
173            return Err(MaterialError::DuplicateElement(z));
174        }
175        out.push((z, if by_mass { f / el.atomic_weight } else { f }));
176    }
177    let sum: f64 = out.iter().map(|&(_, f)| f).sum();
178    if !(sum.is_finite() && sum > 0.0) {
179        return Err(MaterialError::NotNormalizable(sum));
180    }
181    for p in &mut out {
182        p.1 /= sum;
183    }
184    Ok(out)
185}
186
187impl Material {
188    fn build(
189        parts: &[(u8, f64)],
190        by_mass: bool,
191        density_kg_m3: Option<f64>,
192    ) -> Result<Self, MaterialError> {
193        let norm = normalise(parts, by_mass)?;
194        let density = match density_kg_m3 {
195            Some(d) if d.is_finite() && d > 0.0 => d,
196            Some(d) => return Err(MaterialError::InvalidDensity(d)),
197            None if norm.len() == 1 => element(norm[0].0)
198                .and_then(|e| e.density_kg_m3())
199                .ok_or(MaterialError::MissingDensity)?,
200            None => return Err(MaterialError::MissingDensity),
201        };
202        let components = norm
203            .into_iter()
204            .map(|(z, x)| {
205                // `normalise` already checked the element exists.
206                let el = element(z).expect("validated");
207                Component {
208                    z,
209                    atom_fraction: x,
210                    e_d_ev: el.default_displacement_ev,
211                    e_b_ev: Some(0.0),
212                    e_s_ev: el.default_surface_binding_ev,
213                }
214            })
215            .collect();
216        Ok(Self {
217            name: None,
218            components,
219            density_kg_m3: density,
220        })
221    }
222
223    /// Build from `(Z, atom fraction)` pairs. Fractions are relative weights
224    /// (see the module docs). `density_kg_m3` may be `None` only for a single
225    /// element with a tabulated density.
226    pub fn from_atom_fractions(
227        parts: &[(u8, f64)],
228        density_kg_m3: Option<f64>,
229    ) -> Result<Self, MaterialError> {
230        Self::build(parts, false, density_kg_m3)
231    }
232
233    /// Build from `(Z, mass fraction)` pairs. Same rules as
234    /// [`Material::from_atom_fractions`].
235    pub fn from_mass_fractions(
236        parts: &[(u8, f64)],
237        density_kg_m3: Option<f64>,
238    ) -> Result<Self, MaterialError> {
239        Self::build(parts, true, density_kg_m3)
240    }
241
242    /// Set a display name (does not affect physics).
243    pub fn with_name(mut self, name: impl Into<String>) -> Self {
244        self.name = Some(name.into());
245        self
246    }
247
248    /// Display name, if set.
249    pub fn name(&self) -> Option<&str> {
250        self.name.as_deref()
251    }
252
253    /// Components in the order they were given.
254    pub fn components(&self) -> &[Component] {
255        &self.components
256    }
257
258    /// Mass density, kg/m³.
259    pub fn mass_density(&self) -> f64 {
260        self.density_kg_m3
261    }
262
263    /// Override the mass density, kg/m³.
264    pub fn set_mass_density(&mut self, density_kg_m3: f64) -> Result<(), MaterialError> {
265        if !(density_kg_m3.is_finite() && density_kg_m3 > 0.0) {
266            return Err(MaterialError::InvalidDensity(density_kg_m3));
267        }
268        self.density_kg_m3 = density_kg_m3;
269        Ok(())
270    }
271
272    /// Mean atomic weight `sum(x_j M_j)`, g/mol.
273    pub fn mean_atomic_weight(&self) -> f64 {
274        self.components
275            .iter()
276            .map(|c| c.atom_fraction * element(c.z).expect("validated").atomic_weight)
277            .sum()
278    }
279
280    /// Total atom number density, atoms per m³.
281    pub fn atom_number_density(&self) -> f64 {
282        // rho [kg/m^3] -> g/m^3 is 1e3; divide by g/mol for mol/m^3.
283        self.density_kg_m3 * 1.0e3 * AVOGADRO / self.mean_atomic_weight()
284    }
285
286    /// Number density of one element, atoms per m³.
287    pub fn number_density_of(&self, z: u8) -> Result<f64, MaterialError> {
288        let c = self.component(z)?;
289        Ok(c.atom_fraction * self.atom_number_density())
290    }
291
292    /// Atom fractions in component order.
293    pub fn atom_fractions(&self) -> Vec<(u8, f64)> {
294        self.components
295            .iter()
296            .map(|c| (c.z, c.atom_fraction))
297            .collect()
298    }
299
300    /// Mass fractions in component order.
301    pub fn mass_fractions(&self) -> Vec<(u8, f64)> {
302        let mean = self.mean_atomic_weight();
303        self.components
304            .iter()
305            .map(|c| {
306                let m = element(c.z).expect("validated").atomic_weight;
307                (c.z, c.atom_fraction * m / mean)
308            })
309            .collect()
310    }
311
312    fn component(&self, z: u8) -> Result<&Component, MaterialError> {
313        self.components
314            .iter()
315            .find(|c| c.z == z)
316            .ok_or(MaterialError::ElementNotInMaterial(z))
317    }
318
319    fn component_mut(&mut self, z: u8) -> Result<&mut Component, MaterialError> {
320        self.components
321            .iter_mut()
322            .find(|c| c.z == z)
323            .ok_or(MaterialError::ElementNotInMaterial(z))
324    }
325
326    /// Displacement energy `E_d` of element `z` in this material, eV.
327    pub fn displacement_energy_ev(&self, z: u8) -> Result<f64, MaterialError> {
328        self.component(z)?
329            .e_d_ev
330            .ok_or(MaterialError::EnergyNotSet {
331                z,
332                kind: EnergyKind::Displacement,
333            })
334    }
335
336    /// Lattice binding energy `E_b` of element `z`, eV.
337    pub fn lattice_binding_energy_ev(&self, z: u8) -> Result<f64, MaterialError> {
338        self.component(z)?
339            .e_b_ev
340            .ok_or(MaterialError::EnergyNotSet {
341                z,
342                kind: EnergyKind::LatticeBinding,
343            })
344    }
345
346    /// Surface binding energy `E_s` of element `z`, eV.
347    pub fn surface_binding_energy_ev(&self, z: u8) -> Result<f64, MaterialError> {
348        self.component(z)?
349            .e_s_ev
350            .ok_or(MaterialError::EnergyNotSet {
351                z,
352                kind: EnergyKind::SurfaceBinding,
353            })
354    }
355
356    /// Set `E_d` for element `z`, eV.
357    pub fn set_displacement_energy_ev(&mut self, z: u8, ev: f64) -> Result<(), MaterialError> {
358        let v = check_energy(z, EnergyKind::Displacement, ev)?;
359        self.component_mut(z)?.e_d_ev = Some(v);
360        Ok(())
361    }
362
363    /// Set `E_b` for element `z`, eV.
364    pub fn set_lattice_binding_energy_ev(&mut self, z: u8, ev: f64) -> Result<(), MaterialError> {
365        let v = check_energy(z, EnergyKind::LatticeBinding, ev)?;
366        self.component_mut(z)?.e_b_ev = Some(v);
367        Ok(())
368    }
369
370    /// Set `E_s` for element `z`, eV.
371    pub fn set_surface_binding_energy_ev(&mut self, z: u8, ev: f64) -> Result<(), MaterialError> {
372        let v = check_energy(z, EnergyKind::SurfaceBinding, ev)?;
373        self.component_mut(z)?.e_s_ev = Some(v);
374        Ok(())
375    }
376
377    /// Energies that still have no value, in component order: a run should
378    /// refuse to start while this is non-empty.
379    pub fn unset_energies(&self) -> Vec<(u8, EnergyKind)> {
380        let mut v = Vec::new();
381        for c in &self.components {
382            if c.e_d_ev.is_none() {
383                v.push((c.z, EnergyKind::Displacement));
384            }
385            if c.e_b_ev.is_none() {
386                v.push((c.z, EnergyKind::LatticeBinding));
387            }
388            if c.e_s_ev.is_none() {
389                v.push((c.z, EnergyKind::SurfaceBinding));
390            }
391        }
392        v
393    }
394}
395
396/// Serialized form of one element entry. Give exactly one of `symbol` / `z`
397/// and exactly one of `atom_fraction` / `mass_fraction`; every element of a
398/// material must use the same fraction kind.
399#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
400#[serde(deny_unknown_fields)]
401pub struct ElementSpec {
402    /// Element symbol, e.g. `"Si"`.
403    #[serde(default, skip_serializing_if = "Option::is_none")]
404    pub symbol: Option<String>,
405    /// Atomic number (alternative to `symbol`).
406    #[serde(default, skip_serializing_if = "Option::is_none")]
407    pub z: Option<u8>,
408    /// Relative atom fraction (or stoichiometric count).
409    #[serde(default, skip_serializing_if = "Option::is_none")]
410    pub atom_fraction: Option<f64>,
411    /// Relative mass fraction.
412    #[serde(default, skip_serializing_if = "Option::is_none")]
413    pub mass_fraction: Option<f64>,
414    /// Displacement energy override, eV.
415    #[serde(default, skip_serializing_if = "Option::is_none")]
416    pub e_d_ev: Option<f64>,
417    /// Lattice binding energy override, eV.
418    #[serde(default, skip_serializing_if = "Option::is_none")]
419    pub e_b_ev: Option<f64>,
420    /// Surface binding energy override, eV.
421    #[serde(default, skip_serializing_if = "Option::is_none")]
422    pub e_s_ev: Option<f64>,
423}
424
425/// Serialized form of a [`Material`] (what the CLI's TOML describes).
426#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
427#[serde(deny_unknown_fields)]
428pub struct MaterialSpec {
429    /// Display name.
430    #[serde(default, skip_serializing_if = "Option::is_none")]
431    pub name: Option<String>,
432    /// Mass density override, g/cm³. Required for compounds.
433    #[serde(default, skip_serializing_if = "Option::is_none")]
434    pub density_g_cm3: Option<f64>,
435    /// Elements, in order.
436    pub elements: Vec<ElementSpec>,
437}
438
439impl TryFrom<MaterialSpec> for Material {
440    type Error = MaterialError;
441
442    fn try_from(spec: MaterialSpec) -> Result<Self, MaterialError> {
443        let bad = |s: &str| MaterialError::InvalidSpec(s.to_string());
444        let mut parts = Vec::with_capacity(spec.elements.len());
445        let mut kind: Option<bool> = None; // Some(true) = by mass
446        for e in &spec.elements {
447            let z = match (&e.symbol, e.z) {
448                (Some(s), None) => {
449                    element_by_symbol(s)
450                        .ok_or_else(|| MaterialError::UnknownSymbol(s.clone()))?
451                        .z
452                }
453                (None, Some(z)) => z,
454                (Some(s), Some(z)) => {
455                    let el = element_by_symbol(s)
456                        .ok_or_else(|| MaterialError::UnknownSymbol(s.clone()))?;
457                    if el.z != z {
458                        return Err(bad("symbol and z disagree"));
459                    }
460                    z
461                }
462                (None, None) => return Err(bad("element needs `symbol` or `z`")),
463            };
464            let (by_mass, f) = match (e.atom_fraction, e.mass_fraction) {
465                (Some(f), None) => (false, f),
466                (None, Some(f)) => (true, f),
467                _ => return Err(bad("give exactly one of atom_fraction / mass_fraction")),
468            };
469            match kind {
470                None => kind = Some(by_mass),
471                Some(k) if k != by_mass => {
472                    return Err(bad("cannot mix atom_fraction and mass_fraction"))
473                }
474                _ => {}
475            }
476            parts.push((z, f));
477        }
478        let density = spec.density_g_cm3.map(g_cm3_to_kg_m3);
479        let mut m = Material::build(&parts, kind.unwrap_or(false), density)?;
480        m.name = spec.name;
481        for (e, &(z, _)) in spec.elements.iter().zip(&parts) {
482            if let Some(v) = e.e_d_ev {
483                m.set_displacement_energy_ev(z, v)?;
484            }
485            if let Some(v) = e.e_b_ev {
486                m.set_lattice_binding_energy_ev(z, v)?;
487            }
488            if let Some(v) = e.e_s_ev {
489                m.set_surface_binding_energy_ev(z, v)?;
490            }
491        }
492        Ok(m)
493    }
494}
495
496impl From<Material> for MaterialSpec {
497    fn from(m: Material) -> Self {
498        MaterialSpec {
499            name: m.name,
500            density_g_cm3: Some(kg_m3_to_g_cm3(m.density_kg_m3)),
501            elements: m
502                .components
503                .iter()
504                .map(|c| ElementSpec {
505                    symbol: Some(element(c.z).expect("validated").symbol.to_string()),
506                    z: None,
507                    atom_fraction: Some(c.atom_fraction),
508                    mass_fraction: None,
509                    e_d_ev: c.e_d_ev,
510                    e_b_ev: c.e_b_ev,
511                    e_s_ev: c.e_s_ev,
512                })
513                .collect(),
514        }
515    }
516}
517
518#[cfg(test)]
519mod tests {
520    use super::*;
521
522    fn rel(a: f64, b: f64) -> f64 {
523        ((a - b) / b).abs()
524    }
525
526    // Hand values: n = rho[kg/m^3] * 1e3 * N_A * (atoms per formula unit)
527    //                  / (sum of atomic weights over the formula unit)
528    // with N_A = 6.02214076e23, weights Si 28.085, O 15.999, Ga 69.723,
529    // As 74.921595, C 12.011, and these densities (g/cm^3):
530    //   Si   2.329  (table value, CRC)
531    //   SiO2 2.648  (alpha-quartz, 20 C; crystallographic density)
532    //   GaAs 5.3176 (zinc blende, 300 K; crystallographic density)
533    //   SiC  3.21   (3C/6H SiC, room temperature)
534    // Results, per m^3 (rho in g/m^3 = 1e6 * the g/cm^3 value):
535    //   Si   4.99397e28 m^-3   (4.994e22 cm^-3)
536    //   SiO2 3 * 2.648e6*N_A/60.083   = 7.9623e28 m^-3
537    //   GaAs 2 * 5.3176e6*N_A/144.644595 = 4.4279e28 m^-3
538    //   SiC  2 * 3.21e6*N_A/40.096    = 9.6424e28 m^-3
539    #[test]
540    fn si_number_density() {
541        let m = Material::from_atom_fractions(&[(14, 1.0)], None).unwrap();
542        assert!(rel(m.atom_number_density(), 4.99397e28) < 1e-5);
543    }
544
545    #[test]
546    fn sio2_number_density() {
547        let m = Material::from_atom_fractions(&[(14, 1.0), (8, 2.0)], Some(2648.0)).unwrap();
548        assert!(rel(m.atom_number_density(), 7.9623e28) < 1e-4);
549        assert!(rel(m.number_density_of(8).unwrap(), 7.9623e28 * 2.0 / 3.0) < 1e-4);
550    }
551
552    #[test]
553    fn gaas_number_density() {
554        let m = Material::from_atom_fractions(&[(31, 1.0), (33, 1.0)], Some(5317.6)).unwrap();
555        assert!(rel(m.atom_number_density(), 4.4279e28) < 1e-4);
556    }
557
558    #[test]
559    fn sic_number_density() {
560        let m = Material::from_atom_fractions(&[(14, 1.0), (6, 1.0)], Some(3210.0)).unwrap();
561        assert!(rel(m.atom_number_density(), 9.6424e28) < 1e-4);
562    }
563
564    #[test]
565    fn mass_atom_fraction_round_trip() {
566        let a = Material::from_atom_fractions(&[(14, 1.0), (8, 2.0)], Some(2648.0)).unwrap();
567        let w = a.mass_fractions();
568        let sum: f64 = w.iter().map(|p| p.1).sum();
569        assert!((sum - 1.0).abs() < 1e-14);
570        // Hand value: w_Si = 28.085 / 60.083 = 0.467435...
571        assert!((w[0].1 - 28.085 / 60.083).abs() < 1e-12);
572        let b = Material::from_mass_fractions(&w, Some(2648.0)).unwrap();
573        for (x, y) in a.atom_fractions().iter().zip(b.atom_fractions()) {
574            assert_eq!(x.0, y.0);
575            assert!((x.1 - y.1).abs() < 1e-14);
576        }
577        assert!(rel(a.atom_number_density(), b.atom_number_density()) < 1e-14);
578    }
579
580    #[test]
581    fn error_paths() {
582        use MaterialError::*;
583        assert_eq!(Material::from_atom_fractions(&[], Some(1.0)), Err(Empty));
584        assert_eq!(
585            Material::from_atom_fractions(&[(0, 1.0)], Some(1.0)),
586            Err(UnknownAtomicNumber(0))
587        );
588        assert_eq!(
589            Material::from_atom_fractions(&[(93, 1.0)], Some(1.0)),
590            Err(UnknownAtomicNumber(93))
591        );
592        assert!(matches!(
593            Material::from_atom_fractions(&[(14, -1.0)], None),
594            Err(InvalidFraction { z: 14, .. })
595        ));
596        assert!(matches!(
597            Material::from_atom_fractions(&[(14, f64::NAN)], None),
598            Err(InvalidFraction { .. })
599        ));
600        assert!(matches!(
601            Material::from_atom_fractions(&[(14, 0.0), (8, 0.0)], Some(1.0)),
602            Err(NotNormalizable(_))
603        ));
604        assert_eq!(
605            Material::from_atom_fractions(&[(14, 1.0), (14, 1.0)], Some(1.0)),
606            Err(DuplicateElement(14))
607        );
608        assert!(matches!(
609            Material::from_atom_fractions(&[(14, 1.0)], Some(-5.0)),
610            Err(InvalidDensity(_))
611        ));
612        assert_eq!(
613            Material::from_atom_fractions(&[(14, 1.0), (8, 2.0)], None),
614            Err(MissingDensity)
615        );
616        // Gas: no tabulated solid density.
617        assert_eq!(
618            Material::from_atom_fractions(&[(18, 1.0)], None),
619            Err(MissingDensity)
620        );
621        let mut m = Material::from_atom_fractions(&[(14, 1.0)], None).unwrap();
622        assert!(matches!(
623            m.set_surface_binding_energy_ev(14, -1.0),
624            Err(InvalidEnergy { .. })
625        ));
626        assert_eq!(
627            m.set_displacement_energy_ev(8, 10.0),
628            Err(ElementNotInMaterial(8))
629        );
630        // Oxygen has no default E_s: must be set by the user.
631        let mut o = Material::from_atom_fractions(&[(14, 1.0), (8, 2.0)], Some(2648.0)).unwrap();
632        assert_eq!(
633            o.surface_binding_energy_ev(8),
634            Err(EnergyNotSet {
635                z: 8,
636                kind: EnergyKind::SurfaceBinding
637            })
638        );
639        assert!(!o.unset_energies().is_empty());
640        o.set_surface_binding_energy_ev(8, 2.0).unwrap();
641        assert_eq!(o.surface_binding_energy_ev(8), Ok(2.0));
642    }
643
644    #[test]
645    fn defaults_and_overrides() {
646        let mut m = Material::from_atom_fractions(&[(14, 1.0)], None).unwrap();
647        assert_eq!(m.surface_binding_energy_ev(14), Ok(4.63));
648        assert_eq!(m.lattice_binding_energy_ev(14), Ok(0.0));
649        m.set_surface_binding_energy_ev(14, 3.0).unwrap();
650        assert_eq!(m.surface_binding_energy_ev(14), Ok(3.0));
651    }
652
653    const TOML_FIXTURE: &str = r#"
654name = "SiO2"
655density_g_cm3 = 2.648
656
657[[elements]]
658symbol = "Si"
659mass_fraction = 28.085
660e_d_ev = 21.0
661
662[[elements]]
663z = 8
664mass_fraction = 31.998
665e_s_ev = 2.0
666"#;
667
668    #[test]
669    fn toml_fixture_and_serde_round_trip() {
670        let m: Material = toml::from_str(TOML_FIXTURE).unwrap();
671        assert_eq!(m.name(), Some("SiO2"));
672        assert!(rel(m.atom_number_density(), 7.9623e28) < 1e-4);
673        assert_eq!(m.displacement_energy_ev(14), Ok(21.0));
674        assert_eq!(m.surface_binding_energy_ev(8), Ok(2.0));
675        let s = toml::to_string(&m).unwrap();
676        let back: Material = toml::from_str(&s).unwrap();
677        assert_eq!(m.components().len(), back.components().len());
678        assert!(rel(m.atom_number_density(), back.atom_number_density()) < 1e-14);
679        assert_eq!(back.surface_binding_energy_ev(8), Ok(2.0));
680        // Deterministic output.
681        assert_eq!(s, toml::to_string(&back).unwrap());
682        let j = serde_json::to_string(&m).unwrap();
683        let _: Material = serde_json::from_str(&j).unwrap();
684    }
685
686    #[test]
687    fn bad_specs_are_rejected() {
688        let mixed = r#"
689[[elements]]
690symbol = "Si"
691atom_fraction = 1.0
692[[elements]]
693symbol = "O"
694mass_fraction = 1.0
695"#;
696        assert!(toml::from_str::<Material>(mixed).is_err());
697        let unknown = "[[elements]]\nsymbol = \"Xx\"\natom_fraction = 1.0\n";
698        assert!(toml::from_str::<Material>(unknown).is_err());
699        assert!(toml::from_str::<Material>("elements = []").is_err());
700    }
701}