1use serde::{Deserialize, Serialize};
22
23use super::hist::{Binning, Histogram};
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 escapes: Vec<(FaceAcc, FaceAcc)>,
279 last_site: ReplacementTracker,
281}
282
283impl IonTally {
284 pub fn new(
291 stack: &Stack,
292 species_z: &[u8],
293 config: IonTallyConfig,
294 ) -> Result<Self, MaterialError> {
295 let mut consts = Vec::with_capacity(stack.layers().len());
296 for l in stack.layers() {
297 let m = l.material();
298 let mut elems = Vec::new();
299 let mut mean_z = 0.0;
300 for (z, x) in m.atom_fractions() {
301 elems.push((
302 z,
303 m.displacement_energy_ev(z)?,
304 m.lattice_binding_energy_ev(z)?,
305 ));
306 mean_z += x * f64::from(z);
307 }
308 consts.push(LayerConsts {
309 elems,
310 mean_z,
311 mean_a: m.mean_atomic_weight(),
312 });
313 }
314 let layers = (0..consts.len())
315 .map(|_| LayerAcc {
316 stopped: Moments::new(),
317 depth_hist: Histogram::new(config.depth),
318 nrt: NrtDamage::default(),
319 defects: DefectCounts::default(),
320 })
321 .collect();
322 Ok(Self {
323 config,
324 species_z: species_z.to_vec(),
325 consts,
326 histories: 0,
327 budget: EnergyBudget::default(),
328 depth: Moments::new(),
329 lateral_y: Moments::new(),
330 lateral_z: Moments::new(),
331 radial: Moments::new(),
332 depth_hist: Histogram::new(config.depth),
333 lateral_y_hist: Histogram::new(config.lateral),
334 lateral_z_hist: Histogram::new(config.lateral),
335 radial_hist: Histogram::new(config.radial),
336 layers,
337 displacement_hist: Histogram::new(config.depth),
338 replacement_hist: Histogram::new(config.depth),
339 interstitial_hist: Histogram::new(config.depth),
340 escapes: species_z
341 .iter()
342 .map(|_| (FaceAcc::new(&config), FaceAcc::new(&config)))
343 .collect(),
344 last_site: ReplacementTracker::default(),
345 })
346 }
347
348 pub fn histories(&self) -> u64 {
350 self.histories
351 }
352
353 fn energies(&self, layer: usize, z: u8) -> (f64, f64) {
355 self.consts[layer]
356 .elems
357 .iter()
358 .find(|e| e.0 == z)
359 .map(|e| (e.1, e.2))
360 .expect("recoils are elements of their layer")
361 }
362
363 pub fn report(&self, fit_dual_pearson: bool) -> IonReport {
367 let h = self.histories;
368 let depth = self.depth.summary();
369 let (pearson_iv, pearson_iv_error) = match depth {
370 None => (None, None),
371 Some(s) => match PearsonIv::from_moments(s.mean, s.std_dev, s.skewness, s.kurtosis) {
372 Ok(p) => (Some(p), None),
373 Err(e) => (None, Some(e)),
374 },
375 };
376 let (dual_pearson, dual_pearson_error) = if fit_dual_pearson {
377 match DualPearson::fit(&self.depth_hist) {
378 Ok(f) => (Some(f), None),
379 Err(e) => (None, Some(e)),
380 }
381 } else {
382 (None, None)
383 };
384 let range = RangeReport {
385 stopped: self.depth.n,
386 depth,
387 pearson_iv,
388 pearson_iv_error,
389 dual_pearson,
390 dual_pearson_error,
391 depth_histogram: self.depth_hist.clone(),
392 lateral_y: self.lateral_y.summary(),
393 lateral_z: self.lateral_z.summary(),
394 radial: self.radial.summary(),
395 lateral_y_histogram: self.lateral_y_hist.clone(),
396 lateral_z_histogram: self.lateral_z_hist.clone(),
397 radial_histogram: self.radial_hist.clone(),
398 layers: self
399 .layers
400 .iter()
401 .enumerate()
402 .map(|(i, l)| LayerRange {
403 layer: i,
404 stopped: l.stopped.n,
405 depth: l.stopped.summary(),
406 depth_histogram: l.depth_hist.clone(),
407 })
408 .collect(),
409 };
410 let mut nrt = NrtDamage::default();
411 let mut defects = DefectCounts::default();
412 for l in &self.layers {
413 nrt.merge(&l.nrt);
414 defects.merge(&l.defects);
415 }
416 let mut vacancy_hist = self.displacement_hist.clone();
417 for (v, r) in vacancy_hist
418 .counts
419 .iter_mut()
420 .zip(&self.replacement_hist.counts)
421 {
422 *v -= r;
423 }
424 vacancy_hist.underflow -= self.replacement_hist.underflow;
425 vacancy_hist.overflow -= self.replacement_hist.overflow;
426 let damage = DamageReport {
427 nrt,
428 cascade: defects.report(),
429 layers: self
430 .layers
431 .iter()
432 .enumerate()
433 .map(|(i, l)| LayerDamage {
434 layer: i,
435 nrt: l.nrt,
436 cascade: l.defects.report(),
437 })
438 .collect(),
439 vacancy_histogram: vacancy_hist,
440 interstitial_histogram: self.interstitial_hist.clone(),
441 replacement_histogram: self.replacement_hist.clone(),
442 };
443 let species: Vec<SpeciesEscape> = self
444 .escapes
445 .iter()
446 .enumerate()
447 .map(|(i, (f, b))| SpeciesEscape {
448 species: i,
449 z: self.species_z[i],
450 beam: i == 0,
451 front: f.report(h),
452 back: b.report(h),
453 })
454 .collect();
455 let sputtered: u64 = self.escapes.iter().skip(1).map(|e| e.0.count).sum();
456 let escapes = EscapeReport {
457 backscatter_coefficient: per(self.escapes[0].0.count as f64, h),
458 transmission_coefficient: per(self.escapes[0].1.count as f64, h),
459 energy_reflection_coefficient: if self.budget.incident > 0.0 {
460 self.budget.backscattered / self.budget.incident
461 } else {
462 0.0
463 },
464 sputter_yield: per(sputtered as f64, h),
465 species,
466 };
467 IonReport {
468 histories: h,
469 species_z: self.species_z.clone(),
470 budget: self.budget,
471 range,
472 damage,
473 escapes,
474 }
475 }
476}
477
478impl BcaTally for IonTally {
479 fn begin_history(&mut self, _index: u64) {
480 self.last_site.clear();
481 }
482
483 fn lattice(&mut self, _at: [f64; 3], _layer: usize, kind: LatticeDeposit, _energy_ev: f64) {
484 if kind == LatticeDeposit::Subthreshold {
485 self.last_site.clear();
486 }
487 }
488
489 fn recoil(&mut self, r: &Particle) {
490 let layer = r.layer;
491 let (e_d, e_b) = self.energies(layer, r.z);
492 let acc = &mut self.layers[layer];
493 acc.defects.displacements += 1;
494 self.displacement_hist.fill(r.pos[0]);
495 if r.generation == 1 {
496 let c = &self.consts[layer];
499 let t = r.energy_ev + e_b;
500 let t_dam = damage_energy_ev(t, f64::from(r.z), r.mass_amu, c.mean_z, c.mean_a);
501 acc.nrt.pka_count += 1;
502 acc.nrt.pka_energy_ev += t;
503 acc.nrt.damage_energy_ev += t_dam;
504 acc.nrt.nrt_displacements += nrt_displacements(t_dam, e_d);
505 acc.nrt.kinchin_pease_displacements += kinchin_pease(t_dam, e_d);
506 }
507 self.last_site.recoil(r);
508 }
509
510 fn stopped(&mut self, p: &Particle) {
511 let replaced = self.last_site.stopped(p).is_some();
512 let x = p.pos[0];
513 if replaced {
514 self.layers[p.layer].defects.replacements += 1;
515 self.replacement_hist.fill(x);
516 } else if !p.is_primary() {
517 self.layers[p.layer].defects.interstitials += 1;
518 self.interstitial_hist.fill(x);
519 }
520 if p.is_primary() {
521 let (y, z) = (p.pos[1], p.pos[2]);
522 let r = y.hypot(z);
523 self.depth.push(x);
524 self.lateral_y.push(y);
525 self.lateral_z.push(z);
526 self.radial.push(r);
527 self.depth_hist.fill(x);
528 self.lateral_y_hist.fill(y);
529 self.lateral_z_hist.fill(z);
530 self.radial_hist.fill(r);
531 let l = &mut self.layers[p.layer];
532 l.stopped.push(x);
533 l.depth_hist.fill(x);
534 }
535 }
536
537 fn escaped(&mut self, p: &Particle, face: Face) {
538 self.last_site.clear();
539 let polar = p.dir[0].abs().min(1.0).acos();
540 let pair = &mut self.escapes[p.species];
541 let acc = match face {
542 Face::Front => &mut pair.0,
543 Face::Back => &mut pair.1,
544 Face::Side => return,
547 };
548 acc.count += 1;
549 acc.energy_ev += p.energy_ev;
550 acc.energy_hist.fill(p.energy_ev);
551 acc.polar_hist.fill(polar);
552 }
553
554 fn end_history(&mut self, _index: u64, budget: &EnergyBudget) {
555 self.histories += 1;
556 self.budget.add(budget);
557 self.last_site.clear();
558 }
559
560 fn merge(&mut self, o: Self) {
566 assert!(
567 self.config == o.config && self.species_z == o.species_z && self.consts == o.consts,
568 "merging ion tallies built for different runs"
569 );
570 self.histories += o.histories;
571 self.budget.add(&o.budget);
572 self.depth.merge(&o.depth);
573 self.lateral_y.merge(&o.lateral_y);
574 self.lateral_z.merge(&o.lateral_z);
575 self.radial.merge(&o.radial);
576 self.depth_hist.merge(&o.depth_hist);
577 self.lateral_y_hist.merge(&o.lateral_y_hist);
578 self.lateral_z_hist.merge(&o.lateral_z_hist);
579 self.radial_hist.merge(&o.radial_hist);
580 for (a, b) in self.layers.iter_mut().zip(&o.layers) {
581 a.stopped.merge(&b.stopped);
582 a.depth_hist.merge(&b.depth_hist);
583 a.nrt.merge(&b.nrt);
584 a.defects.merge(&b.defects);
585 }
586 self.displacement_hist.merge(&o.displacement_hist);
587 self.replacement_hist.merge(&o.replacement_hist);
588 self.interstitial_hist.merge(&o.interstitial_hist);
589 for (a, b) in self.escapes.iter_mut().zip(&o.escapes) {
590 a.0.merge(&b.0);
591 a.1.merge(&b.1);
592 }
593 }
594}
595
596#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
598pub struct IonReport {
599 pub histories: u64,
601 pub species_z: Vec<u8>,
603 pub budget: EnergyBudget,
605 pub range: RangeReport,
607 pub damage: DamageReport,
609 pub escapes: EscapeReport,
611}
612
613#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
615pub struct RangeReport {
616 pub stopped: u64,
618 pub depth: Option<MomentSummary>,
621 pub pearson_iv: Option<PearsonIv>,
623 pub pearson_iv_error: Option<PearsonError>,
625 pub dual_pearson: Option<DualPearsonFit>,
627 pub dual_pearson_error: Option<PearsonError>,
629 pub depth_histogram: Histogram,
631 pub lateral_y: Option<MomentSummary>,
633 pub lateral_z: Option<MomentSummary>,
635 pub radial: Option<MomentSummary>,
638 pub lateral_y_histogram: Histogram,
640 pub lateral_z_histogram: Histogram,
642 pub radial_histogram: Histogram,
644 pub layers: Vec<LayerRange>,
646}
647
648#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
650pub struct LayerRange {
651 pub layer: usize,
653 pub stopped: u64,
655 pub depth: Option<MomentSummary>,
657 pub depth_histogram: Histogram,
659}
660
661#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
663pub struct DamageReport {
664 pub nrt: NrtDamage,
666 pub cascade: CascadeDefects,
668 pub layers: Vec<LayerDamage>,
670 pub vacancy_histogram: Histogram,
672 pub interstitial_histogram: Histogram,
674 pub replacement_histogram: Histogram,
676}
677
678#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
680pub struct LayerDamage {
681 pub layer: usize,
683 pub nrt: NrtDamage,
685 pub cascade: CascadeDefects,
687}
688
689#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
691pub struct EscapeReport {
692 pub backscatter_coefficient: f64,
694 pub transmission_coefficient: f64,
696 pub energy_reflection_coefficient: f64,
699 pub sputter_yield: f64,
701 pub species: Vec<SpeciesEscape>,
703}
704
705#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
707pub struct SpeciesEscape {
708 pub species: usize,
710 pub z: u8,
712 pub beam: bool,
715 pub front: FaceEscape,
717 pub back: FaceEscape,
719}
720
721#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
723pub struct FaceEscape {
724 pub count: u64,
726 pub per_ion: f64,
729 pub energy_ev: f64,
731 pub mean_energy_ev: f64,
733 pub energy_histogram: Histogram,
735 pub polar_histogram: Histogram,
737}