lindhard/ion/dynamic/
profiles.rs

1//! Lateral implanted-species and displacement-damage profiles of a 2-D
2//! dynamic run, the film cut-through result and the comparison of a milled
3//! profile with published points.
4//!
5//! # What is tallied
6//!
7//! [`DynamicRun2d::with_profiles`](super::DynamicRun2d::with_profiles) bins the
8//! events of every primary on a [`VoxelTallyGrid`] in the sample frame
9//! (`x` the depth from the original front surface, `y` the lateral
10//! coordinate), reusing [`VoxelTally`] for the event definitions:
11//!
12//! * `implanted`: beam particles that came to rest in the voxel;
13//! * `displacements`: target atoms displaced in the voxel (recoils created
14//!   above the displacement energy), the simulated cascade count. These are
15//!   the definitions of the 1-D [`IonReport`](crate::tally::ion::IonReport)
16//!   and of the [`VoxelTally`]; the NRT estimate is not tallied. Amorphization
17//!   and damage accumulation are not modelled (the count is per ion, summed
18//!   over the run).
19//!
20//! The counts of a step are scaled by the fluence each primary stands for
21//! (`F p_j / count` ions/m² for a primary entering column `j`) and by the
22//! column width `h`, so a voxel holds **atoms per metre of the invariant
23//! axis `z`** (`*_per_m`). Dividing by the voxel cross-section `dx dy` gives
24//! the number density (`*_per_m3`). Counts of the substrate are also kept
25//! apart (`substrate_displacements_per_m`): they are the events whose region
26//! was the substrate layer at the time, which is the damage left in the layer
27//! beneath the film. Rejected adaptive attempts are not accumulated, and the
28//! sums are made in column order, so profiles are bit-identical at any thread
29//! count.
30//!
31//! # Export
32//!
33//! [`ProfileSnapshot::write_csv`] follows the conventions of
34//! [`VoxelReport::write_csv`](crate::tally::voxel::VoxelReport::write_csv): a
35//! `# key = value` header carrying the grid, a column row, then one row per
36//! occupied voxel in flat-index order. The snapshot also implements `serde::Serialize`
37//! (JSON with any serde format crate).
38
39use std::collections::BTreeMap;
40use std::io::{self, Write};
41
42use serde::{Deserialize, Serialize};
43
44use crate::tally::voxel::{VoxelReport, VoxelTallyGrid};
45
46/// Accumulated tallies of one voxel, in atoms per metre of `z`.
47#[derive(Debug, Clone, Copy, PartialEq, Default, Serialize, Deserialize)]
48pub struct ProfileCell {
49    /// Beam atoms at rest here.
50    pub implanted_per_m: f64,
51    /// Target atoms displaced here (all layers).
52    pub displacements_per_m: f64,
53    /// Of those, displaced while in the substrate layer.
54    pub substrate_displacements_per_m: f64,
55}
56
57impl ProfileCell {
58    fn add(&mut self, o: &ProfileCell) {
59        self.implanted_per_m += o.implanted_per_m;
60        self.displacements_per_m += o.displacements_per_m;
61        self.substrate_displacements_per_m += o.substrate_displacements_per_m;
62    }
63}
64
65/// Running profile sums of a run (or of one step before it is accepted).
66#[derive(Debug, Clone, PartialEq)]
67pub(crate) struct ProfileAccumulator {
68    grid: VoxelTallyGrid,
69    cells: BTreeMap<u64, ProfileCell>,
70    outside: ProfileCell,
71}
72
73impl ProfileAccumulator {
74    pub(crate) fn new(grid: VoxelTallyGrid) -> Self {
75        Self {
76            grid,
77            cells: BTreeMap::new(),
78            outside: ProfileCell::default(),
79        }
80    }
81
82    /// Add the counts of all events (`all`) and of the substrate events
83    /// (`substrate`) of one column, each count standing for `scale` atoms/m.
84    pub(crate) fn add_counts(&mut self, all: &VoxelReport, substrate: &VoxelReport, scale: f64) {
85        for v in &all.voxels {
86            let c = self.cells.entry(self.grid.flat_index(v.index)).or_default();
87            c.implanted_per_m += v.cell.implanted as f64 * scale;
88            c.displacements_per_m += v.cell.displacements as f64 * scale;
89        }
90        self.outside.implanted_per_m += all.outside.implanted as f64 * scale;
91        self.outside.displacements_per_m += all.outside.displacements as f64 * scale;
92        for v in &substrate.voxels {
93            let c = self.cells.entry(self.grid.flat_index(v.index)).or_default();
94            c.substrate_displacements_per_m += v.cell.displacements as f64 * scale;
95        }
96        self.outside.substrate_displacements_per_m +=
97            substrate.outside.displacements as f64 * scale;
98    }
99
100    pub(crate) fn merge(&mut self, o: &ProfileAccumulator) {
101        for (&k, c) in &o.cells {
102            self.cells.entry(k).or_default().add(c);
103        }
104        self.outside.add(&o.outside);
105    }
106
107    pub(crate) fn snapshot(&self, beam_z: u8, fluence_m2: f64, ions: u64) -> ProfileSnapshot {
108        ProfileSnapshot {
109            beam_z,
110            peak_fluence_m2: fluence_m2,
111            ions_per_column: ions,
112            grid: self.grid,
113            voxels: self
114                .cells
115                .iter()
116                .map(|(&k, &cell)| (self.grid.unflatten(k), cell))
117                .collect(),
118            outside: self.outside,
119        }
120    }
121}
122
123/// The profiles at one fluence. See the module docs.
124#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
125pub struct ProfileSnapshot {
126    /// Atomic number of the beam species.
127    pub beam_z: u8,
128    /// Peak fluence delivered when the snapshot was taken, ions/m².
129    pub peak_fluence_m2: f64,
130    /// Primaries per illuminated column used.
131    pub ions_per_column: u64,
132    /// The grid the voxel indices refer to.
133    pub grid: VoxelTallyGrid,
134    /// Occupied voxels, `[i, j, k]` with their tallies, in flat-index order.
135    pub voxels: Vec<([usize; 3], ProfileCell)>,
136    /// Events outside the grid.
137    pub outside: ProfileCell,
138}
139
140impl ProfileSnapshot {
141    /// Sum over all voxels and the outside.
142    pub fn grand_total(&self) -> ProfileCell {
143        let mut t = self.outside;
144        for (_, c) in &self.voxels {
145            t.add(c);
146        }
147        t
148    }
149
150    /// Voxel cross-section `dx dy` (m²), the divisor from atoms/m to
151    /// atoms/m³.
152    pub fn voxel_area_m2(&self) -> f64 {
153        let s = self.grid.spacing_m();
154        s[0] * s[1]
155    }
156
157    /// Laterally integrated profile along the depth axis: entry `i` sums the
158    /// voxels with `x` index `i`.
159    pub fn depth_profile(&self) -> Vec<ProfileCell> {
160        let mut out = vec![ProfileCell::default(); self.grid.dims()[0]];
161        for (ijk, c) in &self.voxels {
162            out[ijk[0]].add(c);
163        }
164        out
165    }
166
167    /// Profile along the lateral axis: entry `j` sums the voxels with `y`
168    /// index `j`.
169    pub fn lateral_profile(&self) -> Vec<ProfileCell> {
170        let mut out = vec![ProfileCell::default(); self.grid.dims()[1]];
171        for (ijk, c) in &self.voxels {
172            out[ijk[1]].add(c);
173        }
174        out
175    }
176
177    /// Write the CSV export (see the module docs); densities are in
178    /// atoms/m³.
179    pub fn write_csv<W: Write>(&self, mut w: W) -> io::Result<()> {
180        let g = &self.grid;
181        let j3 = |a: [f64; 3]| format!("{:e},{:e},{:e}", a[0], a[1], a[2]);
182        writeln!(w, "# lindhard 2d profile v1")?;
183        writeln!(w, "# beam_z = {}", self.beam_z)?;
184        writeln!(w, "# peak_fluence_m2 = {:e}", self.peak_fluence_m2)?;
185        writeln!(w, "# ions_per_column = {}", self.ions_per_column)?;
186        writeln!(w, "# origin_m = {}", j3(g.origin_m()))?;
187        writeln!(w, "# spacing_m = {}", j3(g.spacing_m()))?;
188        let d = g.dims();
189        writeln!(w, "# dims = {},{},{}", d[0], d[1], d[2])?;
190        let o = &self.outside;
191        writeln!(
192            w,
193            "# outside_per_m = {:e},{:e},{:e}",
194            o.implanted_per_m, o.displacements_per_m, o.substrate_displacements_per_m
195        )?;
196        writeln!(
197            w,
198            "i,j,k,x_m,y_m,z_m,implanted_per_m,displacements_per_m,\
199             substrate_displacements_per_m,implanted_per_m3,displacements_per_m3,\
200             substrate_displacements_per_m3"
201        )?;
202        let a = self.voxel_area_m2();
203        for (ijk, c) in &self.voxels {
204            let m = g.center_m(*ijk);
205            writeln!(
206                w,
207                "{},{},{},{:e},{:e},{:e},{:e},{:e},{:e},{:e},{:e},{:e}",
208                ijk[0],
209                ijk[1],
210                ijk[2],
211                m[0],
212                m[1],
213                m[2],
214                c.implanted_per_m,
215                c.displacements_per_m,
216                c.substrate_displacements_per_m,
217                c.implanted_per_m / a,
218                c.displacements_per_m / a,
219                c.substrate_displacements_per_m / a
220            )?;
221        }
222        Ok(())
223    }
224
225    /// The CSV export as a string.
226    pub fn to_csv(&self) -> String {
227        let mut buf = Vec::new();
228        self.write_csv(&mut buf).expect("writing to a Vec");
229        String::from_utf8(buf).expect("the CSV is ASCII")
230    }
231}
232
233/// Result of [`DynamicRun2d::run_until_cut_through`](super::DynamicRun2d::run_until_cut_through).
234#[derive(Debug, Clone, PartialEq)]
235pub struct CutThrough {
236    /// Whether an illuminated column became bare substrate within the run.
237    pub reached: bool,
238    /// The first column that did (the lowest index if several in one step).
239    pub column: Option<usize>,
240    /// Peak fluence at the end of the last step before the cut-through,
241    /// ions/m². The true cut-through fluence lies in
242    /// `(fluence_before_m2, fluence_m2]`; the bracket is one step wide, which
243    /// the adaptive policy bounds by the allowed composition change.
244    pub fluence_before_m2: f64,
245    /// Peak fluence at the end of the step in which the cut-through was seen
246    /// (the fluence delivered if it was not reached), ions/m².
247    pub fluence_m2: f64,
248    /// Primaries per column used.
249    pub ions_done: u64,
250    /// The profiles at the stop, if profiles are on. The damage in the layer
251    /// beneath is in `substrate_displacements_per_m`.
252    pub profiles: Option<ProfileSnapshot>,
253}
254
255/// One published point compared with the simulated profile.
256#[derive(Debug, Clone, Copy, PartialEq)]
257pub struct PointDifference {
258    /// Lateral position of the published point (relative to the beam centre),
259    /// m.
260    pub y_m: f64,
261    /// Published depth, m.
262    pub published_m: f64,
263    /// Simulated depth at that position (linear interpolation), m.
264    pub simulated_m: f64,
265}
266
267/// A milled profile compared with published points. See [`compare_profile`].
268#[derive(Debug, Clone, PartialEq)]
269pub struct ProfileComparison {
270    /// Per-point values, in the order given.
271    pub points: Vec<PointDifference>,
272    /// Root-mean-square depth difference, m.
273    pub rms_m: f64,
274    /// Largest absolute depth difference, m.
275    pub max_abs_m: f64,
276    /// `rms_m` over the largest published depth.
277    pub relative_rms: f64,
278    /// The tolerance applied to `relative_rms`.
279    pub tolerance: f64,
280    /// `relative_rms <= tolerance`.
281    pub within_tolerance: bool,
282}
283
284/// Compare a simulated profile with published points.
285///
286/// `sim_y_m` and `sim_depth_m` are the lateral coordinate of the column
287/// centres (relative to the beam centre) and the depth of the surface there
288/// (increasing `y`); `published` holds `(y, depth)` pairs, m, in the same
289/// frame. The simulation is interpolated linearly; a published point outside
290/// the simulated range is an error, never extrapolated. The metric is the
291/// RMS depth difference over the published depth maximum; the tolerance is
292/// the caller's, with nothing fitted: a gap is reported, not removed.
293pub fn compare_profile(
294    sim_y_m: &[f64],
295    sim_depth_m: &[f64],
296    published: &[(f64, f64)],
297    tolerance: f64,
298) -> Result<ProfileComparison, String> {
299    if sim_y_m.len() != sim_depth_m.len() || sim_y_m.len() < 2 {
300        return Err("the simulated profile needs at least two matching points".into());
301    }
302    if published.is_empty() {
303        return Err("no published points".into());
304    }
305    if !(tolerance.is_finite() && tolerance > 0.0) {
306        return Err("the tolerance must be finite and positive".into());
307    }
308    if sim_y_m.windows(2).any(|w| w[1] <= w[0]) {
309        return Err("the simulated lateral coordinates must increase".into());
310    }
311    let mut points = Vec::with_capacity(published.len());
312    for &(y, d) in published {
313        if !(y.is_finite() && d.is_finite()) {
314            return Err("published points must be finite".into());
315        }
316        let (lo, hi) = (sim_y_m[0], sim_y_m[sim_y_m.len() - 1]);
317        if y < lo || y > hi {
318            return Err(format!(
319                "published point y = {y:e} m is outside [{lo:e}, {hi:e}]"
320            ));
321        }
322        let i = sim_y_m
323            .partition_point(|&v| v <= y)
324            .clamp(1, sim_y_m.len() - 1);
325        let t = (y - sim_y_m[i - 1]) / (sim_y_m[i] - sim_y_m[i - 1]);
326        let s = sim_depth_m[i - 1] + t * (sim_depth_m[i] - sim_depth_m[i - 1]);
327        points.push(PointDifference {
328            y_m: y,
329            published_m: d,
330            simulated_m: s,
331        });
332    }
333    let n = points.len() as f64;
334    let rms = (points
335        .iter()
336        .map(|p| (p.simulated_m - p.published_m).powi(2))
337        .sum::<f64>()
338        / n)
339        .sqrt();
340    let max_abs = points
341        .iter()
342        .map(|p| (p.simulated_m - p.published_m).abs())
343        .fold(0.0, f64::max);
344    let scale = points
345        .iter()
346        .map(|p| p.published_m.abs())
347        .fold(0.0, f64::max);
348    if scale == 0.0 {
349        return Err("the largest published depth is zero".into());
350    }
351    let relative_rms = rms / scale;
352    Ok(ProfileComparison {
353        points,
354        rms_m: rms,
355        max_abs_m: max_abs,
356        relative_rms,
357        tolerance,
358        within_tolerance: relative_rms <= tolerance,
359    })
360}
lindhard/ion/dynamic/profiles.rs

lindhard/ion/dynamic/
profiles.rs
