1use 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#[derive(Debug, Clone, Copy, PartialEq, Eq)]
35pub enum EnergyKind {
36 Displacement,
38 LatticeBinding,
40 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#[derive(Debug, Clone, PartialEq, thiserror::Error)]
56pub enum MaterialError {
57 #[error("material has no elements")]
59 Empty,
60 #[error("unknown atomic number {0}")]
62 UnknownAtomicNumber(u8),
63 #[error("unknown element symbol {0:?}")]
65 UnknownSymbol(String),
66 #[error("element Z={0} listed more than once")]
68 DuplicateElement(u8),
69 #[error("invalid fraction {value} for element Z={z}: must be finite and non-negative")]
71 InvalidFraction {
72 z: u8,
74 value: f64,
76 },
77 #[error("fractions sum to {0}: cannot be normalised")]
79 NotNormalizable(f64),
80 #[error("invalid density {0} kg/m^3: must be finite and positive")]
82 InvalidDensity(f64),
83 #[error("no density given and none can be inferred for this composition")]
86 MissingDensity,
87 #[error("invalid {kind} {value} eV for element Z={z}: must be finite and non-negative")]
89 InvalidEnergy {
90 z: u8,
92 kind: EnergyKind,
94 value: f64,
96 },
97 #[error("element Z={0} is not in this material")]
99 ElementNotInMaterial(u8),
100 #[error("{kind} for element Z={z} has no default; set it explicitly")]
102 EnergyNotSet {
103 z: u8,
105 kind: EnergyKind,
107 },
108 #[error("invalid material specification: {0}")]
110 InvalidSpec(String),
111}
112
113#[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 pub fn z(&self) -> u8 {
126 self.z
127 }
128 pub fn atom_fraction(&self) -> f64 {
130 self.atom_fraction
131 }
132}
133
134#[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
163fn 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 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 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 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 pub fn with_name(mut self, name: impl Into<String>) -> Self {
244 self.name = Some(name.into());
245 self
246 }
247
248 pub fn name(&self) -> Option<&str> {
250 self.name.as_deref()
251 }
252
253 pub fn components(&self) -> &[Component] {
255 &self.components
256 }
257
258 pub fn mass_density(&self) -> f64 {
260 self.density_kg_m3
261 }
262
263 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 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 pub fn atom_number_density(&self) -> f64 {
282 self.density_kg_m3 * 1.0e3 * AVOGADRO / self.mean_atomic_weight()
284 }
285
286 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 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 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 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 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 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 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 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 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 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#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
400#[serde(deny_unknown_fields)]
401pub struct ElementSpec {
402 #[serde(default, skip_serializing_if = "Option::is_none")]
404 pub symbol: Option<String>,
405 #[serde(default, skip_serializing_if = "Option::is_none")]
407 pub z: Option<u8>,
408 #[serde(default, skip_serializing_if = "Option::is_none")]
410 pub atom_fraction: Option<f64>,
411 #[serde(default, skip_serializing_if = "Option::is_none")]
413 pub mass_fraction: Option<f64>,
414 #[serde(default, skip_serializing_if = "Option::is_none")]
416 pub e_d_ev: Option<f64>,
417 #[serde(default, skip_serializing_if = "Option::is_none")]
419 pub e_b_ev: Option<f64>,
420 #[serde(default, skip_serializing_if = "Option::is_none")]
422 pub e_s_ev: Option<f64>,
423}
424
425#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
427#[serde(deny_unknown_fields)]
428pub struct MaterialSpec {
429 #[serde(default, skip_serializing_if = "Option::is_none")]
431 pub name: Option<String>,
432 #[serde(default, skip_serializing_if = "Option::is_none")]
434 pub density_g_cm3: Option<f64>,
435 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; 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 #[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 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 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 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 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}