1use serde::{Deserialize, Serialize};
22
23use super::hist::{Binning, Histogram, WeightedHistogram};
24use super::moments::{MomentSummary, Moments};
25use super::pearson::{DualPearson, DualPearsonFit, PearsonError, PearsonIv};
26use crate::geometry::Stack;
27use crate::ion::bca::{BcaTally, EnergyBudget, Face, LatticeDeposit, Particle};
28use crate::ion::damage::{damage_energy_ev, kinchin_pease, nrt_displacements};
29use crate::material::MaterialError;
30
31#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
33pub struct IonTallyConfig {
34 pub depth: Binning,
37 pub lateral: Binning,
39 pub radial: Binning,
41 pub escape_energy: Binning,
43 pub escape_polar: Binning,
46}
47
48#[derive(Debug, Clone, Copy, PartialEq, Default, Serialize, Deserialize)]
59pub struct NrtDamage {
60 pub pka_count: u64,
62 pub pka_energy_ev: f64,
65 pub damage_energy_ev: f64,
67 pub nrt_displacements: f64,
69 pub kinchin_pease_displacements: f64,
71}
72
73impl NrtDamage {
74 fn merge(&mut self, o: &NrtDamage) {
75 self.pka_count += o.pka_count;
76 self.pka_energy_ev += o.pka_energy_ev;
77 self.damage_energy_ev += o.damage_energy_ev;
78 self.nrt_displacements += o.nrt_displacements;
79 self.kinchin_pease_displacements += o.kinchin_pease_displacements;
80 }
81}
82
83#[derive(Debug, Clone, Copy, PartialEq, Eq, Default, Serialize, Deserialize)]
103pub struct CascadeDefects {
104 pub displacements: u64,
106 pub replacements: u64,
108 pub vacancies: u64,
110 pub interstitials: u64,
112}
113
114#[derive(Debug, Clone, Copy, PartialEq, Default)]
115struct DefectCounts {
116 displacements: u64,
117 replacements: u64,
118 interstitials: u64,
119}
120
121impl DefectCounts {
122 fn merge(&mut self, o: &DefectCounts) {
123 self.displacements += o.displacements;
124 self.replacements += o.replacements;
125 self.interstitials += o.interstitials;
126 }
127
128 fn report(&self) -> CascadeDefects {
129 CascadeDefects {
130 displacements: self.displacements,
131 replacements: self.replacements,
132 vacancies: self.displacements - self.replacements,
133 interstitials: self.interstitials,
134 }
135 }
136}
137
138#[derive(Debug, Clone, PartialEq)]
140struct LayerConsts {
141 elems: Vec<(u8, f64, f64)>,
143 mean_z: f64,
145 mean_a: f64,
147}
148
149#[derive(Debug, Clone, PartialEq)]
150struct LayerAcc {
151 stopped: Moments,
152 depth_hist: Histogram,
153 nrt: NrtDamage,
154 defects: DefectCounts,
155}
156
157#[derive(Debug, Clone, PartialEq)]
158struct FaceAcc {
159 count: u64,
160 energy_ev: f64,
161 energy_hist: Histogram,
162 polar_hist: Histogram,
163}
164
165impl FaceAcc {
166 fn new(c: &IonTallyConfig) -> Self {
167 Self {
168 count: 0,
169 energy_ev: 0.0,
170 energy_hist: Histogram::new(c.escape_energy),
171 polar_hist: Histogram::new(c.escape_polar),
172 }
173 }
174
175 fn merge(&mut self, o: &FaceAcc) {
176 self.count += o.count;
177 self.energy_ev += o.energy_ev;
178 self.energy_hist.merge(&o.energy_hist);
179 self.polar_hist.merge(&o.polar_hist);
180 }
181
182 fn report(&self, histories: u64) -> FaceEscape {
183 FaceEscape {
184 count: self.count,
185 per_ion: per(self.count as f64, histories),
186 energy_ev: self.energy_ev,
187 mean_energy_ev: if self.count > 0 {
188 self.energy_ev / self.count as f64
189 } else {
190 0.0
191 },
192 energy_histogram: self.energy_hist.clone(),
193 polar_histogram: self.polar_hist.clone(),
194 }
195 }
196}
197
198fn per(x: f64, histories: u64) -> f64 {
199 if histories == 0 {
200 0.0
201 } else {
202 x / histories as f64
203 }
204}
205
206#[derive(Debug, Clone, Copy, PartialEq)]
208struct Site {
209 pos: [f64; 3],
210 z: u8,
211 species: usize,
212 displacer_generation: u32,
213}
214
215#[derive(Debug, Clone, Copy, PartialEq, Default)]
223pub(crate) struct ReplacementTracker {
224 site: Option<Site>,
225}
226
227impl ReplacementTracker {
228 pub(crate) fn clear(&mut self) {
230 self.site = None;
231 }
232
233 pub(crate) fn recoil(&mut self, r: &Particle) {
235 self.site = Some(Site {
236 pos: r.pos,
237 z: r.z,
238 species: r.species,
239 displacer_generation: r.generation.saturating_sub(1),
240 });
241 }
242
243 pub(crate) fn stopped(&mut self, p: &Particle) -> Option<usize> {
246 let hit = match self.site {
247 Some(s) if s.pos == p.pos && s.z == p.z && s.displacer_generation == p.generation => {
248 Some(s.species)
249 }
250 _ => None,
251 };
252 self.site = None;
253 hit
254 }
255}
256
257#[derive(Debug, Clone, PartialEq)]
259pub struct IonTally {
260 config: IonTallyConfig,
261 species_z: Vec<u8>,
262 consts: Vec<LayerConsts>,
263 histories: u64,
264 budget: EnergyBudget,
265 depth: Moments,
266 lateral_y: Moments,
267 lateral_z: Moments,
268 radial: Moments,
269 depth_hist: Histogram,
270 lateral_y_hist: Histogram,
271 lateral_z_hist: Histogram,
272 radial_hist: Histogram,
273 layers: Vec<LayerAcc>,
274 displacement_hist: Histogram,
275 replacement_hist: Histogram,
276 interstitial_hist: Histogram,
277 nrt_hist: WeightedHistogram,
279 kp_hist: WeightedHistogram,
281 escapes: Vec<(FaceAcc, FaceAcc)>,
283 last_site: ReplacementTracker,
285}
286
287impl IonTally {
288 pub fn new(
295 stack: &Stack,
296 species_z: &[u8],
297 config: IonTallyConfig,
298 ) -> Result<Self, MaterialError> {
299 let mut consts = Vec::with_capacity(stack.layers().len());
300 for l in stack.layers() {
301 let m = l.material();
302 let mut elems = Vec::new();
303 let mut mean_z = 0.0;
304 for (z, x) in m.atom_fractions() {
305 elems.push((
306 z,
307 m.displacement_energy_ev(z)?,
308 m.lattice_binding_energy_ev(z)?,
309 ));
310 mean_z += x * f64::from(z);
311 }
312 consts.push(LayerConsts {
313 elems,
314 mean_z,
315 mean_a: m.mean_atomic_weight(),
316 });
317 }
318 let layers = (0..consts.len())
319 .map(|_| LayerAcc {
320 stopped: Moments::new(),
321 depth_hist: Histogram::new(config.depth),
322 nrt: NrtDamage::default(),
323 defects: DefectCounts::default(),
324 })
325 .collect();
326 Ok(Self {
327 config,
328 species_z: species_z.to_vec(),
329 consts,
330 histories: 0,
331 budget: EnergyBudget::default(),
332 depth: Moments::new(),
333 lateral_y: Moments::new(),
334 lateral_z: Moments::new(),
335 radial: Moments::new(),
336 depth_hist: Histogram::new(config.depth),
337 lateral_y_hist: Histogram::new(config.lateral),
338 lateral_z_hist: Histogram::new(config.lateral),
339 radial_hist: Histogram::new(config.radial),
340 layers,
341 displacement_hist: Histogram::new(config.depth),
342 replacement_hist: Histogram::new(config.depth),
343 interstitial_hist: Histogram::new(config.depth),
344 nrt_hist: WeightedHistogram::new(config.depth),
345 kp_hist: WeightedHistogram::new(config.depth),
346 escapes: species_z
347 .iter()
348 .map(|_| (FaceAcc::new(&config), FaceAcc::new(&config)))
349 .collect(),
350 last_site: ReplacementTracker::default(),
351 })
352 }
353
354 pub fn histories(&self) -> u64 {
356 self.histories
357 }
358
359 fn energies(&self, layer: usize, z: u8) -> (f64, f64) {
361 self.consts[layer]
362 .elems
363 .iter()
364 .find(|e| e.0 == z)
365 .map(|e| (e.1, e.2))
366 .expect("recoils are elements of their layer")
367 }
368
369 pub fn report(&self, fit_dual_pearson: bool) -> IonReport {
373 let h = self.histories;
374 let depth = self.depth.summary();
375 let (pearson_iv, pearson_iv_error) = match depth {
376 None => (None, None),
377 Some(s) => match PearsonIv::from_moments(s.mean, s.std_dev, s.skewness, s.kurtosis) {
378 Ok(p) => (Some(p), None),
379 Err(e) => (None, Some(e)),
380 },
381 };
382 let (dual_pearson, dual_pearson_error) = if fit_dual_pearson {
383 match DualPearson::fit(&self.depth_hist) {
384 Ok(f) => (Some(f), None),
385 Err(e) => (None, Some(e)),
386 }
387 } else {
388 (None, None)
389 };
390 let range = RangeReport {
391 stopped: self.depth.n,
392 depth,
393 pearson_iv,
394 pearson_iv_error,
395 dual_pearson,
396 dual_pearson_error,
397 depth_histogram: self.depth_hist.clone(),
398 lateral_y: self.lateral_y.summary(),
399 lateral_z: self.lateral_z.summary(),
400 radial: self.radial.summary(),
401 lateral_y_histogram: self.lateral_y_hist.clone(),
402 lateral_z_histogram: self.lateral_z_hist.clone(),
403 radial_histogram: self.radial_hist.clone(),
404 layers: self
405 .layers
406 .iter()
407 .enumerate()
408 .map(|(i, l)| LayerRange {
409 layer: i,
410 stopped: l.stopped.n,
411 depth: l.stopped.summary(),
412 depth_histogram: l.depth_hist.clone(),
413 })
414 .collect(),
415 };
416 let mut nrt = NrtDamage::default();
417 let mut defects = DefectCounts::default();
418 for l in &self.layers {
419 nrt.merge(&l.nrt);
420 defects.merge(&l.defects);
421 }
422 let mut vacancy_hist = self.displacement_hist.clone();
423 for (v, r) in vacancy_hist
424 .counts
425 .iter_mut()
426 .zip(&self.replacement_hist.counts)
427 {
428 *v -= r;
429 }
430 vacancy_hist.underflow -= self.replacement_hist.underflow;
431 vacancy_hist.overflow -= self.replacement_hist.overflow;
432 let damage = DamageReport {
433 nrt,
434 cascade: defects.report(),
435 layers: self
436 .layers
437 .iter()
438 .enumerate()
439 .map(|(i, l)| LayerDamage {
440 layer: i,
441 nrt: l.nrt,
442 cascade: l.defects.report(),
443 })
444 .collect(),
445 vacancy_histogram: vacancy_hist,
446 interstitial_histogram: self.interstitial_hist.clone(),
447 replacement_histogram: self.replacement_hist.clone(),
448 nrt_histogram: self.nrt_hist.clone(),
449 kinchin_pease_histogram: self.kp_hist.clone(),
450 };
451 let species: Vec<SpeciesEscape> = self
452 .escapes
453 .iter()
454 .enumerate()
455 .map(|(i, (f, b))| SpeciesEscape {
456 species: i,
457 z: self.species_z[i],
458 beam: i == 0,
459 front: f.report(h),
460 back: b.report(h),
461 })
462 .collect();
463 let sputtered: u64 = self.escapes.iter().skip(1).map(|e| e.0.count).sum();
464 let escapes = EscapeReport {
465 backscatter_coefficient: per(self.escapes[0].0.count as f64, h),
466 transmission_coefficient: per(self.escapes[0].1.count as f64, h),
467 energy_reflection_coefficient: if self.budget.incident > 0.0 {
468 self.budget.backscattered / self.budget.incident
469 } else {
470 0.0
471 },
472 sputter_yield: per(sputtered as f64, h),
473 species,
474 };
475 IonReport {
476 histories: h,
477 species_z: self.species_z.clone(),
478 budget: self.budget,
479 range,
480 damage,
481 escapes,
482 }
483 }
484}
485
486impl BcaTally for IonTally {
487 fn begin_history(&mut self, _index: u64) {
488 self.last_site.clear();
489 }
490
491 fn lattice(&mut self, _at: [f64; 3], _layer: usize, kind: LatticeDeposit, _energy_ev: f64) {
492 if kind == LatticeDeposit::Subthreshold {
493 self.last_site.clear();
494 }
495 }
496
497 fn recoil(&mut self, r: &Particle) {
498 let layer = r.layer;
499 let (e_d, e_b) = self.energies(layer, r.z);
500 let acc = &mut self.layers[layer];
501 acc.defects.displacements += 1;
502 self.displacement_hist.fill(r.pos[0]);
503 if r.generation == 1 {
504 let c = &self.consts[layer];
507 let t = r.energy_ev + e_b;
508 let t_dam = damage_energy_ev(t, f64::from(r.z), r.mass_amu, c.mean_z, c.mean_a);
509 acc.nrt.pka_count += 1;
510 acc.nrt.pka_energy_ev += t;
511 acc.nrt.damage_energy_ev += t_dam;
512 let n_nrt = nrt_displacements(t_dam, e_d);
513 let n_kp = kinchin_pease(t_dam, e_d);
514 acc.nrt.nrt_displacements += n_nrt;
515 acc.nrt.kinchin_pease_displacements += n_kp;
516 self.nrt_hist.fill(r.pos[0], n_nrt);
517 self.kp_hist.fill(r.pos[0], n_kp);
518 }
519 self.last_site.recoil(r);
520 }
521
522 fn stopped(&mut self, p: &Particle) {
523 let replaced = self.last_site.stopped(p).is_some();
524 let x = p.pos[0];
525 if replaced {
526 self.layers[p.layer].defects.replacements += 1;
527 self.replacement_hist.fill(x);
528 } else if !p.is_primary() {
529 self.layers[p.layer].defects.interstitials += 1;
530 self.interstitial_hist.fill(x);
531 }
532 if p.is_primary() {
533 let (y, z) = (p.pos[1], p.pos[2]);
534 let r = y.hypot(z);
535 self.depth.push(x);
536 self.lateral_y.push(y);
537 self.lateral_z.push(z);
538 self.radial.push(r);
539 self.depth_hist.fill(x);
540 self.lateral_y_hist.fill(y);
541 self.lateral_z_hist.fill(z);
542 self.radial_hist.fill(r);
543 let l = &mut self.layers[p.layer];
544 l.stopped.push(x);
545 l.depth_hist.fill(x);
546 }
547 }
548
549 fn escaped(&mut self, p: &Particle, face: Face) {
550 self.last_site.clear();
551 let polar = p.dir[0].abs().min(1.0).acos();
552 let pair = &mut self.escapes[p.species];
553 let acc = match face {
554 Face::Front => &mut pair.0,
555 Face::Back => &mut pair.1,
556 Face::Side => return,
559 };
560 acc.count += 1;
561 acc.energy_ev += p.energy_ev;
562 acc.energy_hist.fill(p.energy_ev);
563 acc.polar_hist.fill(polar);
564 }
565
566 fn end_history(&mut self, _index: u64, budget: &EnergyBudget) {
567 self.histories += 1;
568 self.budget.add(budget);
569 self.last_site.clear();
570 }
571
572 fn merge(&mut self, o: Self) {
578 assert!(
579 self.config == o.config && self.species_z == o.species_z && self.consts == o.consts,
580 "merging ion tallies built for different runs"
581 );
582 self.histories += o.histories;
583 self.budget.add(&o.budget);
584 self.depth.merge(&o.depth);
585 self.lateral_y.merge(&o.lateral_y);
586 self.lateral_z.merge(&o.lateral_z);
587 self.radial.merge(&o.radial);
588 self.depth_hist.merge(&o.depth_hist);
589 self.lateral_y_hist.merge(&o.lateral_y_hist);
590 self.lateral_z_hist.merge(&o.lateral_z_hist);
591 self.radial_hist.merge(&o.radial_hist);
592 for (a, b) in self.layers.iter_mut().zip(&o.layers) {
593 a.stopped.merge(&b.stopped);
594 a.depth_hist.merge(&b.depth_hist);
595 a.nrt.merge(&b.nrt);
596 a.defects.merge(&b.defects);
597 }
598 self.displacement_hist.merge(&o.displacement_hist);
599 self.replacement_hist.merge(&o.replacement_hist);
600 self.interstitial_hist.merge(&o.interstitial_hist);
601 self.nrt_hist.merge(&o.nrt_hist);
602 self.kp_hist.merge(&o.kp_hist);
603 for (a, b) in self.escapes.iter_mut().zip(&o.escapes) {
604 a.0.merge(&b.0);
605 a.1.merge(&b.1);
606 }
607 }
608}
609
610#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
612pub struct IonReport {
613 pub histories: u64,
615 pub species_z: Vec<u8>,
617 pub budget: EnergyBudget,
619 pub range: RangeReport,
621 pub damage: DamageReport,
623 pub escapes: EscapeReport,
625}
626
627#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
629pub struct RangeReport {
630 pub stopped: u64,
632 pub depth: Option<MomentSummary>,
635 pub pearson_iv: Option<PearsonIv>,
637 pub pearson_iv_error: Option<PearsonError>,
639 pub dual_pearson: Option<DualPearsonFit>,
641 pub dual_pearson_error: Option<PearsonError>,
643 pub depth_histogram: Histogram,
645 pub lateral_y: Option<MomentSummary>,
647 pub lateral_z: Option<MomentSummary>,
649 pub radial: Option<MomentSummary>,
652 pub lateral_y_histogram: Histogram,
654 pub lateral_z_histogram: Histogram,
656 pub radial_histogram: Histogram,
658 pub layers: Vec<LayerRange>,
660}
661
662#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
664pub struct LayerRange {
665 pub layer: usize,
667 pub stopped: u64,
669 pub depth: Option<MomentSummary>,
671 pub depth_histogram: Histogram,
673}
674
675#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
677pub struct DamageReport {
678 pub nrt: NrtDamage,
680 pub cascade: CascadeDefects,
682 pub layers: Vec<LayerDamage>,
684 pub vacancy_histogram: Histogram,
686 pub interstitial_histogram: Histogram,
688 pub replacement_histogram: Histogram,
690 pub nrt_histogram: WeightedHistogram,
694 pub kinchin_pease_histogram: WeightedHistogram,
697}
698
699#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
701pub struct LayerDamage {
702 pub layer: usize,
704 pub nrt: NrtDamage,
706 pub cascade: CascadeDefects,
708}
709
710#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
712pub struct EscapeReport {
713 pub backscatter_coefficient: f64,
715 pub transmission_coefficient: f64,
717 pub energy_reflection_coefficient: f64,
720 pub sputter_yield: f64,
722 pub species: Vec<SpeciesEscape>,
724}
725
726#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
728pub struct SpeciesEscape {
729 pub species: usize,
731 pub z: u8,
733 pub beam: bool,
736 pub front: FaceEscape,
738 pub back: FaceEscape,
740}
741
742#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
744pub struct FaceEscape {
745 pub count: u64,
747 pub per_ion: f64,
750 pub energy_ev: f64,
752 pub mean_energy_ev: f64,
754 pub energy_histogram: Histogram,
756 pub polar_histogram: Histogram,
758}