lindhard/input/electron/

resolve.rs

//! Validation and resolution: turns an [`ElectronInput`] into a
//! [`ResolvedElectron`].

use std::collections::BTreeMap;
use std::path::Path;

use sha2::{Digest, Sha256};

use crate::electron::boundary::{BandModel, BandStructure};
use crate::electron::data::{OpticalElf, Subshell, SubshellBindingTable};
use crate::electron::elastic::corrections::{CorrelationPolarization, PolarizationCutoff};
use crate::electron::elastic::table::DEFAULT_MAX_ENERGY_EV;
use crate::electron::inelastic::{FitWeighting, MerminFitOptions, PennAlgorithm, MAX_OSCILLATORS};
use crate::electron::phonon::{FrohlichPhonon, InsulatorChannels, PolaronTrapping, Sio2LoMode};
use crate::electron::secondary::SecondaryModel;
use crate::electron::transport::{
    BeamProfile, BeamSpread, BoundaryModel, CutoffReference, EscapeRule, Primary, TransportConfig,
};
use crate::elements::element_by_symbol;
use crate::geometry::Stack;
use crate::input::{divergence_from_spec, invalid, resolve_material, InputError, ResolvedLayer};
use crate::ion::crystal::Divergence;
use crate::material::Material;
use crate::tally::{
    BeamProfile as ExposureBeamProfile, Binning, CartesianGrid, CylindricalGrid, DepthDoseConfig,
    DoseConfig, ElectronTallyConfig, ExposureConfig, ExposureVolume, InterfaceRefine,
    LogRadialBinning, PsfConfig, PsfModel, PsfNormalization, VolumeShape, VoxelTallyGrid,
};

use super::resolved::*;
use super::schema::*;

const NM: f64 = 1e-9;

/// XORed into the run seed for the exposure tally's landing-offset stream.
const EXPOSURE_SEED_SALT: u64 = 0x6578_706f_7375_7265;

fn finite_pos(v: f64) -> bool {
    v.is_finite() && v > 0.0
}

/// Validate the landing point, profile and divergence of `[electron.beam]`.
fn resolve_beam_spread(b: &ElectronBeamSpec) -> Result<BeamSpread, InputError> {
    if !b.landing_nm.iter().all(|v| v.is_finite()) {
        return Err(invalid("electron.beam.landing_nm", "must be finite"));
    }
    let size = |field: &str, v: f64| {
        if finite_pos(v) {
            Ok(v * NM)
        } else {
            Err(invalid(
                format!("electron.beam.profile.{field}"),
                "must be finite and positive",
            ))
        }
    };
    let profile = match b.profile {
        None | Some(BeamProfileSpec::Point) => BeamProfile::Point,
        Some(BeamProfileSpec::Gaussian { sigma_nm }) => BeamProfile::Gaussian {
            sigma_m: size("sigma_nm", sigma_nm)?,
        },
        Some(BeamProfileSpec::Disc { radius_nm }) => BeamProfile::Disc {
            radius_m: size("radius_nm", radius_nm)?,
        },
        Some(BeamProfileSpec::Rectangle {
            width_y_nm,
            width_z_nm,
        }) => BeamProfile::Rectangle {
            width_y_m: size("width_y_nm", width_y_nm)?,
            width_z_m: size("width_z_nm", width_z_nm)?,
        },
    };
    let divergence = match &b.divergence {
        None => Divergence::None,
        Some(d) => divergence_from_spec(d, "electron.beam.divergence")?,
    };
    Ok(BeamSpread {
        landing_m: [b.landing_nm[0] * NM, b.landing_nm[1] * NM],
        profile,
        divergence,
    })
}

/// `band` with the optional valence binding of a metal applied.
fn with_valence_binding(
    band: BandStructure,
    binding_ev: Option<f64>,
) -> Result<BandStructure, crate::electron::boundary::BandError> {
    match binding_ev {
        None => Ok(band),
        Some(b) => band.with_valence_binding_ev(b),
    }
}

fn read_data_file(
    base_dir: &Path,
    field: &str,
    path: &str,
) -> Result<(String, DataFile), InputError> {
    let joined = base_dir.join(path);
    let bytes = std::fs::read(&joined)
        .map_err(|e| invalid(field, format!("cannot read {}: {e}", joined.display())))?;
    let text = String::from_utf8(bytes.clone()).map_err(|_| {
        invalid(
            field,
            format!("{} is not valid UTF-8 text", joined.display()),
        )
    })?;
    let sha256 = Sha256::digest(&bytes)
        .iter()
        .map(|b| format!("{b:02x}"))
        .collect();
    Ok((
        text,
        DataFile {
            path: path.to_string(),
            resolved_path: std::fs::canonicalize(&joined).unwrap_or(joined),
            sha256,
        },
    ))
}

fn binning_nm(field: &str, b: &BinsNm) -> Result<Binning, InputError> {
    Binning::new(b.lo_nm * NM, b.hi_nm * NM, b.bins)
        .map_err(|e| invalid(field, format!("{e} (edges in nm)")))
}

/// The smallest stopping threshold over the layers (the transport's rule:
/// the cutoff, or the inner potential plus the cutoff from the vacuum level).
fn lowest_threshold(
    materials: &[ResolvedElectronMaterial],
    layer_material: &[usize],
    reference: CutoffReference,
    cutoff_ev: f64,
) -> f64 {
    layer_material
        .iter()
        .map(|&i| match reference {
            CutoffReference::BandBottom => cutoff_ev,
            CutoffReference::VacuumLevel => {
                materials[i]
                    .band
                    .as_ref()
                    .map_or(0.0, BandStructure::inner_potential_ev)
                    + cutoff_ev
            }
        })
        .fold(f64::INFINITY, f64::min)
}

/// The tally pieces `resolve_in` hands to the resolved input.
struct ResolvedTally {
    tally: ElectronTallyConfig,
    voxel_grid: Option<VoxelTallyGrid>,
    psf_fits: Vec<PsfModel>,
    psf_normalization: PsfNormalization,
}

/// Refuse the options in `fields` when they are set but the choice they
/// belong to (`applies_with`, e.g. `secondaries = "kieft-bosch"`) is not
/// selected.
fn reject_unless_selected(
    fields: [(&str, Option<bool>); 2],
    applies_with: &str,
) -> Result<(), InputError> {
    for (k, v) in fields {
        if v.is_some() {
            return Err(invalid(
                format!("electron.transport.{k}"),
                format!("only applies with {applies_with}"),
            ));
        }
    }
    Ok(())
}

/// Validate the energy and direction of `[electron.beam]`.
fn resolve_primary(b: &ElectronBeamSpec) -> Result<Primary, InputError> {
    if !finite_pos(b.energy_ev) {
        return Err(invalid(
            "electron.beam.energy_ev",
            "must be finite and positive",
        ));
    }
    if !(b.tilt_deg.is_finite() && (0.0..90.0).contains(&b.tilt_deg)) {
        return Err(invalid(
            "electron.beam.tilt_deg",
            "must be in [0, 90) degrees",
        ));
    }
    if !b.azimuth_deg.is_finite() {
        return Err(invalid("electron.beam.azimuth_deg", "must be finite"));
    }
    let (s, c) = b.tilt_deg.to_radians().sin_cos();
    let (sa, ca) = b.azimuth_deg.to_radians().sin_cos();
    Ok(Primary {
        energy_ev: b.energy_ev,
        direction: [c, s * ca, s * sa],
    })
}

/// Validate `[electron.transport]`, fill in the echoed defaults and build
/// the transport configuration. Also returns whether the choices need band
/// structures.
fn resolve_transport(
    t: &TransportSpec,
    echo: &mut ElectronInput,
) -> Result<(TransportConfig, bool), InputError> {
    if !finite_pos(t.cutoff_ev) {
        return Err(invalid(
            "electron.transport.cutoff_ev",
            "must be finite and positive",
        ));
    }
    if t.max_events == 0 {
        return Err(invalid(
            "electron.transport.max_events",
            "must be at least 1",
        ));
    }
    let secondaries = match t.secondaries {
        SecondaryChoice::Off => {
            reject_unless_selected(
                [
                    ("instantaneous_momentum", t.instantaneous_momentum),
                    ("momentum_conservation", t.momentum_conservation),
                ],
                "secondaries = \"kieft-bosch\"",
            )?;
            SecondaryModel::Off
        }
        SecondaryChoice::KieftBosch => {
            let im = t.instantaneous_momentum.unwrap_or(true);
            let mc = t.momentum_conservation.unwrap_or(true);
            echo.electron.transport.instantaneous_momentum = Some(im);
            echo.electron.transport.momentum_conservation = Some(mc);
            SecondaryModel::KieftBosch {
                instantaneous_momentum: im,
                momentum_conservation: mc,
            }
        }
    };
    let boundary = match t.boundary {
        BoundaryChoice::Transparent => {
            reject_unless_selected(
                [
                    ("quantum_transmission", t.quantum_transmission),
                    ("refraction", t.refraction),
                ],
                "boundary = \"step-barrier\"",
            )?;
            BoundaryModel::Transparent
        }
        BoundaryChoice::StepBarrier => {
            let qt = t.quantum_transmission.unwrap_or(true);
            let rf = t.refraction.unwrap_or(true);
            echo.electron.transport.quantum_transmission = Some(qt);
            echo.electron.transport.refraction = Some(rf);
            BoundaryModel::StepBarrier {
                quantum_transmission: qt,
                refraction: rf,
            }
        }
    };
    let cutoff_reference = match t.cutoff_reference {
        CutoffReferenceChoice::BandBottom => CutoffReference::BandBottom,
        CutoffReferenceChoice::VacuumLevel => CutoffReference::VacuumLevel,
    };
    let config = TransportConfig {
        cutoff_ev: t.cutoff_ev,
        escape_rule: match t.escape_rule {
            EscapeRuleChoice::BothFaces => EscapeRule::BothFaces,
            EscapeRuleChoice::FrontOnly => EscapeRule::FrontOnly,
        },
        max_events: t.max_events,
        secondaries,
        boundary,
        cutoff_reference,
    };
    let needs_bands = secondaries != SecondaryModel::Off
        || boundary != BoundaryModel::Transparent
        || cutoff_reference == CutoffReference::VacuumLevel;
    Ok((config, needs_bands))
}

/// Validate `[electron.elastic]`, including the correlation-polarization
/// inputs against the elements of the target `layers`.
fn resolve_elastic(
    el: &ElasticSpec,
    layers: &[ResolvedLayer],
) -> Result<ElasticChoice, InputError> {
    let mut cp = BTreeMap::new();
    if let Some(c) = &el.correlation_polarization {
        let f = "electron.elastic.correlation_polarization";
        let cutoff = match c.b_pol_squared {
            None => PolarizationCutoff::Seltzer,
            Some(b) if b.is_finite() && b >= 0.0 => PolarizationCutoff::BPolSquared(b),
            Some(b) => {
                return Err(invalid(
                    format!("{f}.b_pol_squared"),
                    format!("{b} must be finite and non-negative"),
                ))
            }
        };
        if !finite_pos(c.outer_radius_bohr) {
            return Err(invalid(
                format!("{f}.outer_radius_bohr"),
                "must be finite and positive",
            ));
        }
        let mut zs: Vec<u8> = Vec::new();
        for l in layers {
            for comp in l.material.components() {
                if comp.atom_fraction() > 0.0 && !zs.contains(&comp.z()) {
                    zs.push(comp.z());
                }
            }
        }
        for (sym, p) in &c.polarizability {
            let pf = format!("{f}.polarizability.{sym}");
            let elm = element_by_symbol(sym)
                .ok_or_else(|| invalid(&pf, format!("unknown element symbol {sym:?}")))?;
            if !zs.contains(&elm.z) {
                return Err(invalid(&pf, format!("{sym} is not in any target layer")));
            }
            if !finite_pos(p.bohr3) {
                return Err(invalid(
                    format!("{pf}.bohr3"),
                    "must be finite and positive",
                ));
            }
            if p.source.trim().is_empty() {
                return Err(invalid(
                    format!("{pf}.source"),
                    "a polarizability needs a source (a citation of the value)",
                ));
            }
            cp.insert(
                elm.z,
                CorrelationPolarization {
                    polarizability: p.bohr3,
                    polarizability_source: p.source.trim().to_string(),
                    cutoff,
                    outer_radius: c.outer_radius_bohr,
                },
            );
        }
        for z in zs {
            if !cp.contains_key(&z) {
                let sym = crate::elements::element(z).expect("valid").symbol;
                return Err(invalid(
                    format!("{f}.polarizability"),
                    format!("no polarizability for {sym}; give every target element one"),
                ));
            }
        }
    }
    Ok(ElasticChoice {
        model: el.model,
        potential: el.potential,
        exchange: el.exchange,
        correlation_polarization: cp,
    })
}

/// Validate `[electron.inelastic]` against the resolved `materials`; fills
/// in the echoed Mermin fit floors.
fn resolve_inelastic(
    e: &ElectronSpec,
    materials: &[ResolvedElectronMaterial],
    echo: &mut ElectronInput,
) -> Result<PennAlgorithm, InputError> {
    let inelastic: PennAlgorithm =
        e.inelastic
            .model
            .parse()
            .map_err(|err: crate::electron::data::ElectronDataError| {
                invalid("electron.inelastic.model", err.to_string())
            })?;
    if !(e.inelastic.fermi_energy_ev.is_finite() && e.inelastic.fermi_energy_ev >= 0.0) {
        return Err(invalid(
            "electron.inelastic.fermi_energy_ev",
            "must be finite and non-negative",
        ));
    }
    // The fit options only shape the `mermin-melf` table; anywhere else
    // they would be silently ignored.
    for (name, spec) in &e.materials {
        if spec.mermin_fit.is_some() && inelastic != PennAlgorithm::Mermin {
            return Err(invalid(
                format!("electron.materials.{name}.mermin_fit"),
                format!(
                    "only applies with electron.inelastic.model = \"{}\" (the model is \"{}\")",
                    PennAlgorithm::Mermin.label(),
                    inelastic.label()
                ),
            ));
        }
    }
    // Shell channels are tabulated with the single-pole model only
    // (`build_shell_channel_tables`).
    for (name, spec) in &e.materials {
        if spec.inner_shells.is_some() && inelastic != PennAlgorithm::SinglePole {
            return Err(invalid(
                format!("electron.materials.{name}.inner_shells"),
                format!(
                    "only applies with electron.inelastic.model = \"{}\" (the model is \"{}\")",
                    PennAlgorithm::SinglePole.label(),
                    inelastic.label()
                ),
            ));
        }
    }
    for m in materials {
        if let FitWeighting::Relative { floor } = m.mermin_fit.weighting {
            if let Some(f) = echo
                .electron
                .materials
                .get_mut(&m.name)
                .and_then(|s| s.mermin_fit.as_mut())
            {
                f.relative_floor = Some(floor);
            }
        }
    }
    // The inelastic table of a material with a band takes its Fermi
    // energy from the band (#241); a second one would count it twice.
    if e.inelastic.fermi_energy_ev != 0.0 {
        if let Some(m) = materials.iter().find(|m| m.band.is_some()) {
            return Err(invalid(
                "electron.inelastic.fermi_energy_ev",
                format!(
                    "must be 0 when a material has a band (material {} does): the \
                     inelastic table of a material with a band takes its Fermi energy \
                     from the band",
                    m.name
                ),
            ));
        }
    }
    Ok(inelastic)
}

/// Validate `[electron.tables]` and build the shared energy grid (eV);
/// fills in the echoed upper edge and pushes coverage warnings.
fn resolve_table_grid(
    e: &ElectronSpec,
    config: &TransportConfig,
    materials: &[ResolvedElectronMaterial],
    layer_material: &[usize],
    echo: &mut ElectronInput,
    warnings: &mut Vec<String>,
) -> Result<Vec<f64>, InputError> {
    let g = &e.tables;
    let max_inner = materials
        .iter()
        .filter_map(|m| m.band.as_ref().map(BandStructure::inner_potential_ev))
        .fold(0.0, f64::max);
    let top = match config.boundary {
        BoundaryModel::StepBarrier { .. } => e.beam.energy_ev + max_inner,
        BoundaryModel::Transparent => e.beam.energy_ev,
    };
    let max_energy = g.max_energy_ev.unwrap_or(top);
    echo.electron.tables.max_energy_ev = Some(max_energy);
    if !(finite_pos(g.min_energy_ev) && max_energy.is_finite() && max_energy > g.min_energy_ev) {
        return Err(invalid(
            "electron.tables",
            format!(
                "need 0 < min_energy_ev < max_energy_ev, got {} and {max_energy} eV",
                g.min_energy_ev
            ),
        ));
    }
    let table_energy_ev = crate::electron::elastic::table::log_energy_grid(
        g.min_energy_ev,
        max_energy,
        g.points_per_decade,
    )
    .map_err(|err| invalid("electron.tables", err.to_string()))?;
    if max_energy < top {
        warnings.push(format!(
            "electron.tables.max_energy_ev: {max_energy} eV is below the largest energy \
             an electron can have ({top} eV); rates above it are held at the last row \
             (results.table_coverage counts how often)"
        ));
    }
    // The lowest energy an electron is followed at: the cutoff, or
    // `U + cutoff` from the vacuum level (the transport's thresholds).
    let lowest = lowest_threshold(
        materials,
        layer_material,
        config.cutoff_reference,
        e.transport.cutoff_ev,
    );
    if g.min_energy_ev > lowest {
        warnings.push(format!(
            "electron.tables.min_energy_ev: {} eV is above the lowest stopping threshold \
             ({lowest} eV); rates below it are held at the first row \
             (results.table_coverage counts how often)",
            g.min_energy_ev
        ));
    }
    if max_energy > DEFAULT_MAX_ENERGY_EV {
        warnings.push(format!(
            "electron.tables.max_energy_ev: {max_energy} eV is above the {DEFAULT_MAX_ENERGY_EV} eV \
             the elastic tables are tested to; the inelastic models are nonrelativistic"
        ));
    }
    Ok(table_energy_ev)
}

impl ElectronInput {
    /// Whether a TOML document is an electron input: it has a top-level
    /// `electron` table.
    pub fn is_electron_toml(text: &str) -> bool {
        text.parse::<toml::Table>()
            .map(|t| t.contains_key("electron"))
            .unwrap_or(false)
    }

    /// Parse a TOML document (schema errors name the key and line).
    pub fn from_toml_str(text: &str) -> Result<Self, InputError> {
        if text
            .parse::<toml::Table>()
            .is_ok_and(|t| t.contains_key("crystal"))
        {
            return Err(invalid(
                "crystal",
                "not supported with an [electron] run (crystal transport is for ion runs)",
            ));
        }
        toml::from_str(text).map_err(|e| InputError::Parse(e.to_string()))
    }

    /// The input as a TOML document in the same schema.
    pub fn to_toml_string(&self) -> Result<String, InputError> {
        toml::to_string(self).map_err(|e| InputError::Parse(e.to_string()))
    }

    /// Validate everything with relative data paths taken from the current
    /// directory.
    pub fn resolve(&self) -> Result<ResolvedElectron, InputError> {
        self.resolve_in(Path::new("."))
    }

    /// Validate everything, read the data files (paths relative to
    /// `base_dir`, the input file's directory) and build the engine's
    /// inputs. Errors name the offending key.
    pub fn resolve_in(&self, base_dir: &Path) -> Result<ResolvedElectron, InputError> {
        let mut warnings = Vec::new();
        let e = &self.electron;
        let mut echo = self.clone();
        echo.run.threads = None;

        self.resolve_run()?;

        let primary = resolve_primary(&e.beam)?;
        let beam = resolve_beam_spread(&e.beam)?;

        let (config, needs_bands) = resolve_transport(&e.transport, &mut echo)?;

        // Materials (all of them, so a broken unused entry is still reported).
        for (name, spec) in &self.materials {
            Material::try_from(spec.clone())
                .map_err(|err| invalid(format!("materials.{name}"), err.to_string()))?;
        }

        let (stack, layers, fields) = self.resolve_target()?;

        let (materials, layer_material) =
            self.resolve_layer_materials(base_dir, &layers, &fields, needs_bands, &mut warnings)?;

        let elastic = resolve_elastic(&e.elastic, &layers)?;

        let inelastic = resolve_inelastic(e, &materials, &mut echo)?;

        let table_energy_ev = resolve_table_grid(
            e,
            &config,
            &materials,
            &layer_material,
            &mut echo,
            &mut warnings,
        )?;

        let ResolvedTally {
            tally,
            voxel_grid,
            psf_fits,
            psf_normalization,
        } = self.resolve_tally(&mut echo)?;

        Ok(ResolvedElectron {
            input: echo,
            stack,
            layers,
            layer_material,
            materials,
            config,
            primary,
            beam,
            elastic,
            inelastic,
            inelastic_fermi_ev: e.inelastic.fermi_energy_ev,
            table_energy_ev,
            tally,
            voxel_grid,
            psf_fits,
            psf_normalization,
            warnings,
        })
    }

    /// Validate `[run]`.
    fn resolve_run(&self) -> Result<(), InputError> {
        if self.run.histories == 0 {
            return Err(invalid("run.histories", "must be at least 1"));
        }
        if self.run.threads == Some(0) {
            return Err(invalid(
                "run.threads",
                "must be at least 1 (omit for all cores)",
            ));
        }
        Ok(())
    }

    /// Validate `[target]` and build the stack. Returns the stack, the
    /// layers in stack order (the substrate last) and the input key of each
    /// layer's material, for error messages.
    fn resolve_target(&self) -> Result<(Stack, Vec<ResolvedLayer>, Vec<String>), InputError> {
        let tg = &self.target;
        if tg.layers.is_empty() && tg.substrate.is_none() {
            return Err(invalid(
                "target",
                "needs at least one [[target.layers]] entry or a substrate",
            ));
        }
        let mut layers = Vec::new();
        let mut fields = Vec::new();
        let mut thick = Vec::new();
        for (i, l) in tg.layers.iter().enumerate() {
            let field = format!("target.layers[{i}]");
            if !finite_pos(l.thickness_nm) {
                return Err(invalid(
                    format!("{field}.thickness_nm"),
                    format!("{} nm must be finite and positive", l.thickness_nm),
                ));
            }
            let f = format!("{field}.material");
            layers.push(resolve_material(&self.materials, &f, &l.material)?);
            fields.push(f);
            thick.push(l.thickness_nm * NM);
        }
        if let Some(sub) = &tg.substrate {
            layers.push(resolve_material(&self.materials, "target.substrate", sub)?);
            fields.push("target.substrate".to_string());
        }
        let substrate = tg.substrate.as_ref().map(|_| {
            layers
                .last()
                .expect("substrate pushed last")
                .material
                .clone()
        });
        let finite: Vec<(Material, f64)> = layers
            .iter()
            .zip(&thick)
            .map(|(l, &t)| (l.material.clone(), t))
            .collect();
        let stack =
            Stack::new(finite, substrate).map_err(|err| invalid("target", err.to_string()))?;
        Ok((stack, layers, fields))
    }

    /// Electron data, one entry per material name the target uses. Returns
    /// the materials and, for each layer, the index of its material.
    fn resolve_layer_materials(
        &self,
        base_dir: &Path,
        layers: &[ResolvedLayer],
        fields: &[String],
        needs_bands: bool,
        warnings: &mut Vec<String>,
    ) -> Result<(Vec<ResolvedElectronMaterial>, Vec<usize>), InputError> {
        let e = &self.electron;
        let mut materials: Vec<ResolvedElectronMaterial> = Vec::new();
        let mut layer_material = Vec::with_capacity(layers.len());
        for (l, field) in layers.iter().zip(fields) {
            if l.source == "inline" {
                return Err(invalid(
                    field,
                    "an electron run needs a named material (a [materials] key or an \
                     element symbol), so that [electron.materials.<name>] can give its \
                     electron data",
                ));
            }
            if let Some(i) = materials.iter().position(|m| m.name == l.source) {
                layer_material.push(i);
                continue;
            }
            let name = &l.source;
            let mfield = format!("electron.materials.{name}");
            let spec = e.materials.get(name).ok_or_else(|| {
                invalid(
                    field,
                    format!(
                        "material {name:?} has no electron data; add [{mfield}] with at \
                         least optical_elf"
                    ),
                )
            })?;
            materials.push(self.resolve_electron_material(
                base_dir,
                &mfield,
                name,
                &l.material,
                spec,
                needs_bands,
            )?);
            layer_material.push(materials.len() - 1);
        }
        for name in e.materials.keys() {
            if !materials.iter().any(|m| &m.name == name) {
                warnings.push(format!(
                    "electron.materials.{name}: not used by any target layer"
                ));
            }
        }
        Ok((materials, layer_material))
    }

    /// Validate `[electron.tally]` and build the tally configuration; fills
    /// in the echoed defaults.
    fn resolve_tally(&self, echo: &mut ElectronInput) -> Result<ResolvedTally, InputError> {
        let e = &self.electron;
        let b = &e.beam;
        let ts = &e.tally;
        let f = "electron.tally";
        if !(ts.se_bse_split_ev.is_finite() && ts.se_bse_split_ev >= 0.0) {
            return Err(invalid(
                format!("{f}.se_bse_split_ev"),
                "must be finite and non-negative",
            ));
        }
        let emax = ts.escape_energy_max_ev.unwrap_or(b.energy_ev);
        echo.electron.tally.escape_energy_max_ev = Some(emax);
        let escape_energy = Binning::new(0.0, emax, ts.escape_energy_bins).map_err(|err| {
            invalid(
                format!("{f}.escape_energy_max_ev / escape_energy_bins"),
                err.to_string(),
            )
        })?;
        if !(finite_pos(ts.escape_polar_max_deg) && ts.escape_polar_max_deg <= 180.0) {
            return Err(invalid(
                format!("{f}.escape_polar_max_deg"),
                "must be in (0, 180] degrees",
            ));
        }
        let escape_polar = Binning::new(
            0.0,
            ts.escape_polar_max_deg
                .to_radians()
                .min(std::f64::consts::PI),
            ts.escape_polar_bins,
        )
        .map_err(|err| invalid(format!("{f}.escape_polar_bins"), err.to_string()))?;
        let mut tally = ElectronTallyConfig::new(escape_energy, escape_polar);
        tally.se_bse_split_ev = ts.se_bse_split_ev;
        if let Some(c) = &ts.cartesian {
            tally.cartesian = Some(CartesianGrid {
                x: binning_nm(&format!("{f}.cartesian.x"), &c.x)?,
                y: binning_nm(&format!("{f}.cartesian.y"), &c.y)?,
                z: binning_nm(&format!("{f}.cartesian.z"), &c.z)?,
            });
        }
        if let Some(c) = &ts.cylindrical {
            if c.r.lo_nm < 0.0 {
                return Err(invalid(format!("{f}.cylindrical.r.lo_nm"), "must be >= 0"));
            }
            tally.cylindrical = Some(CylindricalGrid {
                r: binning_nm(&format!("{f}.cylindrical.r"), &c.r)?,
                depth: binning_nm(&format!("{f}.cylindrical.depth"), &c.depth)?,
            });
        }
        if let Some(d) = &ts.dose {
            let depth = match d.depth_bins_per_layer {
                None => {
                    if d.refine_width_nm.is_some()
                        || d.refine_bins.is_some()
                        || d.substrate_depth_nm.is_some()
                    {
                        return Err(invalid(
                            format!("{f}.dose.depth_bins_per_layer"),
                            "required when a depth-profile key is given",
                        ));
                    }
                    None
                }
                Some(bins_per_layer) => {
                    if bins_per_layer == 0 {
                        return Err(invalid(
                            format!("{f}.dose.depth_bins_per_layer"),
                            "must be at least 1",
                        ));
                    }
                    let refine = match d.refine_width_nm {
                        None => None,
                        Some(w) => {
                            if !(w.is_finite() && w >= 0.0) {
                                return Err(invalid(
                                    format!("{f}.dose.refine_width_nm"),
                                    "must be finite and non-negative",
                                ));
                            }
                            let bins = d.refine_bins.unwrap_or(10);
                            if bins == 0 {
                                return Err(invalid(
                                    format!("{f}.dose.refine_bins"),
                                    "must be at least 1",
                                ));
                            }
                            Some(InterfaceRefine {
                                width_m: w * NM,
                                bins,
                            })
                        }
                    };
                    Some(DepthDoseConfig {
                        bins_per_layer,
                        refine,
                        substrate_depth_m: d.substrate_depth_nm.unwrap_or(0.0) * NM,
                    })
                }
            };
            tally.dose = Some(DoseConfig {
                depth,
                transmission: d.transmission,
            });
        }
        let mut voxel_grid = None;
        if let Some(v) = &ts.voxel {
            for a in 0..3 {
                if !v.origin_nm[a].is_finite() {
                    return Err(invalid(format!("{f}.voxel.origin_nm"), "must be finite"));
                }
                if !finite_pos(v.spacing_nm[a]) {
                    return Err(invalid(
                        format!("{f}.voxel.spacing_nm"),
                        "must be finite and positive",
                    ));
                }
                if v.dims[a] == 0 {
                    return Err(invalid(format!("{f}.voxel.dims"), "must be at least 1"));
                }
            }
            voxel_grid = Some(
                VoxelTallyGrid::new(
                    v.origin_nm.map(|x| x * NM),
                    v.spacing_nm.map(|x| x * NM),
                    v.dims,
                )
                .map_err(|e| invalid(format!("{f}.voxel"), e.to_string()))?,
            );
        }
        let mut psf_fits: Vec<PsfModel> = Vec::new();
        let mut psf_normalization = PsfNormalization::default();
        if let Some(p) = &ts.psf {
            let nm = |key: &str, v: f64| -> Result<f64, InputError> {
                if v.is_finite() {
                    Ok(v * NM)
                } else {
                    Err(invalid(format!("{f}.psf.{key}"), "must be finite"))
                }
            };
            let depth_lo = nm("depth_lo_nm", p.depth_lo_nm)?;
            let depth_hi = nm("depth_hi_nm", p.depth_hi_nm)?;
            if depth_hi <= depth_lo {
                return Err(invalid(
                    format!("{f}.psf.depth_hi_nm"),
                    "must be greater than depth_lo_nm",
                ));
            }
            let r_min = nm("r_min_nm", p.r_min_nm)?;
            let r_max = nm("r_max_nm", p.r_max_nm)?;
            if r_min <= 0.0 {
                return Err(invalid(format!("{f}.psf.r_min_nm"), "must be positive"));
            }
            if r_max <= r_min {
                return Err(invalid(
                    format!("{f}.psf.r_max_nm"),
                    "must be greater than r_min_nm",
                ));
            }
            if p.bins == 0 {
                return Err(invalid(format!("{f}.psf.bins"), "must be at least 1"));
            }
            let radial = LogRadialBinning::new(r_min, r_max, p.bins)
                .map_err(|e| invalid(format!("{f}.psf.bins"), e.to_string()))?;
            tally.psf = Some(
                PsfConfig::new(radial, depth_lo, depth_hi)
                    .map_err(|e| invalid(format!("{f}.psf"), e.to_string()))?,
            );
            for m in p.fits.iter().copied().map(PsfModel::from) {
                if !psf_fits.contains(&m) {
                    psf_fits.push(m);
                }
            }
            psf_normalization = p.normalization.into();
        }

        if let Some(x) = &ts.exposure {
            tally.exposure = Some(self.resolve_exposure(x, echo)?);
        }
        Ok(ResolvedTally {
            tally,
            voxel_grid,
            psf_fits,
            psf_normalization,
        })
    }

    /// Validate `[electron.tally.exposure]` and build its configuration.
    fn resolve_exposure(
        &self,
        x: &ExposureSpec,
        echo: &mut ElectronInput,
    ) -> Result<ExposureConfig, InputError> {
        let k = "electron.tally.exposure".to_string();
        let nm = |key: &str, v: f64| -> Result<f64, InputError> {
            if v.is_finite() {
                Ok(v * NM)
            } else {
                Err(invalid(format!("{k}.{key}"), "must be finite"))
            }
        };
        let beam = match x.beam {
            ExposureBeamSpec::Point => ExposureBeamProfile::Point,
            ExposureBeamSpec::Gaussian { sigma_nm } => ExposureBeamProfile::Gaussian {
                sigma_m: nm("beam.sigma_nm", sigma_nm)?,
            },
            ExposureBeamSpec::UniformDisc { radius_nm } => ExposureBeamProfile::UniformDisc {
                radius_m: nm("beam.radius_nm", radius_nm)?,
            },
        };
        let mut volumes = Vec::new();
        for (i, v) in x.volumes.iter().enumerate() {
            let key = format!("volumes[{i}]");
            if v.name.is_empty() || x.volumes[..i].iter().any(|o| o.name == v.name) {
                return Err(invalid(
                    format!("{k}.{key}.name"),
                    "must be non-empty and unique",
                ));
            }
            let need = |name: &str, v: Option<f64>| -> Result<f64, InputError> {
                let v = v.ok_or_else(|| {
                    invalid(format!("{k}.{key}.{name}"), "required for this shape")
                })?;
                nm(&format!("{key}.{name}"), v)
            };
            let shape = match v.shape {
                ExposureShapeSpec::Cylinder => VolumeShape::Cylinder {
                    radius_m: need("radius_nm", v.radius_nm)?,
                },
                ExposureShapeSpec::Box => VolumeShape::Box {
                    half_y_m: need("half_y_nm", v.half_y_nm)?,
                    half_z_m: need("half_z_nm", v.half_z_nm)?,
                },
            };
            volumes.push(ExposureVolume {
                name: v.name.clone(),
                shape,
                depth_lo_m: nm(&format!("{key}.depth_lo_nm"), v.depth_lo_nm)?,
                depth_hi_m: nm(&format!("{key}.depth_hi_nm"), v.depth_hi_nm)?,
                offset_y_m: nm(&format!("{key}.offset_y_nm"), v.offset_y_nm)?,
                offset_z_m: nm(&format!("{key}.offset_z_nm"), v.offset_z_nm)?,
            });
        }
        let x_max_ev = x.x_max_ev.unwrap_or(self.electron.beam.energy_ev);
        echo.electron.tally.exposure.as_mut().unwrap().x_max_ev = Some(x_max_ev);
        let cfg = ExposureConfig {
            volumes,
            beam,
            // A stream independent of the transport's (same run seed).
            seed: self.run.seed ^ EXPOSURE_SEED_SALT,
            bins: x.bins,
            x_max_ev,
            lambda: x.lambda,
            quantiles: x.quantiles.clone(),
            below_ev: x.below_ev.clone(),
            above_ev: x.above_ev.clone(),
            tail_tolerance: x.tail_tolerance,
        };
        cfg.validate()
            .map_err(|e| invalid(k.clone(), e.to_string()))?;
        Ok(cfg)
    }

    fn resolve_electron_material(
        &self,
        base_dir: &Path,
        field: &str,
        name: &str,
        material: &Material,
        spec: &ElectronMaterialSpec,
        needs_bands: bool,
    ) -> Result<ResolvedElectronMaterial, InputError> {
        let ef = format!("{field}.optical_elf");
        let (text, file) = read_data_file(base_dir, &ef, &spec.optical_elf)?;
        // The loader refuses a missing or blank provenance.
        let optical_elf = OpticalElf::from_toml_str(&text)
            .map_err(|err| invalid(&ef, format!("{}: {err}", spec.optical_elf)))?;

        let band = match &spec.band {
            None if needs_bands => {
                return Err(invalid(
                    format!("{field}.band"),
                    "required with secondaries, the step barrier or the vacuum-level \
                     cutoff (they need the material's Fermi energy and inner potential)",
                ))
            }
            None => None,
            Some(bs) => {
                let r = match bs {
                    BandSpec::Metal {
                        fermi_ev,
                        work_function_ev,
                        valence_binding_ev,
                        provenance,
                    } => BandStructure::new(
                        BandModel::Metal {
                            fermi_ev: *fermi_ev,
                            work_function_ev: *work_function_ev,
                        },
                        provenance.trim(),
                    )
                    .and_then(|b| with_valence_binding(b, *valence_binding_ev)),
                    BandSpec::Insulator {
                        valence_band_width_ev,
                        band_gap_ev,
                        affinity_ev,
                        provenance,
                    } => BandStructure::new(
                        BandModel::Insulator {
                            valence_band_width_ev: *valence_band_width_ev,
                            band_gap_ev: *band_gap_ev,
                            affinity_ev: *affinity_ev,
                        },
                        provenance.trim(),
                    ),
                    BandSpec::FreeElectronMetal {
                        valence_electrons_per_atom,
                        work_function_ev,
                        valence_binding_ev,
                        provenance,
                    } => BandStructure::free_electron_metal(
                        material,
                        *valence_electrons_per_atom,
                        *work_function_ev,
                        provenance.trim(),
                    )
                    .and_then(|b| with_valence_binding(b, *valence_binding_ev)),
                };
                Some(r.map_err(|err| invalid(format!("{field}.band"), err.to_string()))?)
            }
        };

        let mut channels = InsulatorChannels::none();
        if let Some(p) = &spec.phonon {
            let pf = format!("{field}.phonon");
            let ph = match p.preset.as_deref() {
                Some(preset) => {
                    if p.hbar_omega_ev.is_some()
                        || p.eps_static.is_some()
                        || p.eps_high_frequency.is_some()
                        || p.provenance.is_some()
                    {
                        return Err(invalid(
                            &pf,
                            "a preset fixes the parameters and their source; give only \
                             preset and temperature_k",
                        ));
                    }
                    let mode = match preset {
                        "sio2-63mev" => Sio2LoMode::Mev63,
                        "sio2-153mev" => Sio2LoMode::Mev153,
                        other => {
                            return Err(invalid(
                                format!("{pf}.preset"),
                                format!(
                                    "unknown preset {other:?}; expected \"sio2-63mev\" or \
                                     \"sio2-153mev\""
                                ),
                            ))
                        }
                    };
                    FrohlichPhonon::sio2(mode, p.temperature_k)
                }
                None => {
                    let need = |k: &str, v: Option<f64>| {
                        v.ok_or_else(|| invalid(format!("{pf}.{k}"), "required without a preset"))
                    };
                    FrohlichPhonon::new(
                        need("hbar_omega_ev", p.hbar_omega_ev)?,
                        need("eps_static", p.eps_static)?,
                        need("eps_high_frequency", p.eps_high_frequency)?,
                        p.temperature_k,
                        p.provenance.as_deref().unwrap_or("").trim(),
                    )
                }
            };
            channels.phonon = Some(ph.map_err(|err| invalid(&pf, err.to_string()))?);
        }
        if let Some(p) = &spec.polaron {
            channels.polaron = Some(
                PolaronTrapping::new(p.c_per_nm / NM, p.gamma_per_ev, p.provenance.trim())
                    .map_err(|err| invalid(format!("{field}.polaron"), err.to_string()))?,
            );
        }

        let mermin_fit = match &spec.mermin_fit {
            None => MerminFitOptions::default(),
            Some(f) => mermin_fit_options(&format!("{field}.mermin_fit"), f)?,
        };

        let inner_shells = match &spec.inner_shells {
            None => None,
            Some(i) => Some(resolve_inner_shells(
                base_dir,
                &format!("{field}.inner_shells"),
                material,
                i,
            )?),
        };

        let e_pair = match &spec.e_pair {
            None => None,
            Some(p) => {
                if !(p.ev.is_finite() && p.ev > 0.0) {
                    return Err(invalid(
                        format!("{field}.e_pair.ev"),
                        format!("must be finite and positive, got {}", p.ev),
                    ));
                }
                if p.source.trim().is_empty() {
                    return Err(invalid(
                        format!("{field}.e_pair.source"),
                        "a citation is required: E_pair has no default",
                    ));
                }
                Some(ResolvedEPair {
                    ev: p.ev,
                    source: p.source.trim().to_string(),
                })
            }
        };

        Ok(ResolvedElectronMaterial {
            name: name.to_string(),
            material: material.clone(),
            optical_elf,
            optical_elf_file: file,
            band,
            channels,
            mermin_fit,
            inner_shells,
            e_pair,
        })
    }
}

/// The subshells and binding table of an `inner_shells` table; errors name
/// `field` and the key. The table must hold every requested subshell of every
/// element of `material`.
fn resolve_inner_shells(
    base_dir: &Path,
    field: &str,
    material: &Material,
    spec: &InnerShellsSpec,
) -> Result<ResolvedInnerShells, InputError> {
    let sf = format!("{field}.subshells");
    if spec.subshells.is_empty() {
        return Err(invalid(&sf, "must list at least one subshell"));
    }
    let mut subshells: Vec<Subshell> = Vec::new();
    for label in &spec.subshells {
        let s = match Subshell::from_label(label) {
            Some(s) if SUPPORTED_SUBSHELLS.contains(&s.label()) => s,
            _ => {
                return Err(invalid(
                    &sf,
                    format!(
                        "unknown or unsupported subshell {label:?}; expected one of \"K\", \
                         \"L1\", \"L2\", \"L3\" (no hydrogenic formula is implemented for the M \
                         shell and above)"
                    ),
                ))
            }
        };
        if subshells.contains(&s) {
            return Err(invalid(&sf, format!("{label:?} is listed twice")));
        }
        subshells.push(s);
    }
    let (binding, binding_file) = match &spec.binding {
        None => (SubshellBindingTable::eadl2017(), None),
        Some(path) => {
            let bf = format!("{field}.binding");
            let (text, file) = read_data_file(base_dir, &bf, path)?;
            // The loader refuses a missing or blank provenance.
            let table = SubshellBindingTable::from_toml_str(&text)
                .map_err(|err| invalid(&bf, format!("{path}: {err}")))?;
            (table, Some(file))
        }
    };
    for c in material.components() {
        if c.atom_fraction() <= 0.0 {
            continue;
        }
        let z = c.z();
        let sym = crate::elements::element(z).map_or("?", |e| e.symbol);
        if binding.atom(z).is_none() {
            return Err(invalid(
                format!("{field}.binding"),
                format!("the table has no entry for {sym} (Z = {z}), an element of the material"),
            ));
        }
        for s in &subshells {
            if binding.atom(z).and_then(|a| a.shell(*s)).is_none() {
                return Err(invalid(
                    &sf,
                    format!(
                        "{} is not in the binding table for {sym} (Z = {z}), an element of \
                         the material",
                        s.label()
                    ),
                ));
            }
        }
    }
    Ok(ResolvedInnerShells {
        subshells,
        binding,
        binding_file,
    })
}

/// The subshells with a hydrogenic formula (`electron::inelastic::shell_elf`).
const SUPPORTED_SUBSHELLS: [&str; 4] = ["K", "L1", "L2", "L3"];

/// The fit options of a `mermin_fit` table; errors name `field` and the key.
fn mermin_fit_options(field: &str, f: &MerminFitSpec) -> Result<MerminFitOptions, InputError> {
    if !(1..=MAX_OSCILLATORS).contains(&f.oscillators) {
        return Err(invalid(
            format!("{field}.oscillators"),
            format!("must be in 1..={MAX_OSCILLATORS}, got {}", f.oscillators),
        ));
    }
    let weighting = match f.weighting {
        MerminWeightingChoice::Uniform => {
            if f.relative_floor.is_some() {
                return Err(invalid(
                    format!("{field}.relative_floor"),
                    "only applies with weighting = \"relative\"",
                ));
            }
            FitWeighting::Uniform
        }
        MerminWeightingChoice::Relative => {
            let floor = f.relative_floor.unwrap_or_else(default_relative_floor);
            if !(floor.is_finite() && floor > 0.0 && floor <= 1.0) {
                return Err(invalid(
                    format!("{field}.relative_floor"),
                    format!("must be in (0, 1], got {floor}"),
                ));
            }
            FitWeighting::Relative { floor }
        }
    };
    Ok(MerminFitOptions {
        n_oscillators: f.oscillators,
        weighting,
        ..MerminFitOptions::default()
    })
}