refeff-io 0.2.0

FEFF file-format readers and writers (feff.inp, .dat/.bin handoffs, PAD encoding) for the refeff FEFF10 port
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142

use ndarray::{Array2, ArrayView3};
use refeff_core::{
    Real, SfconvMomentumSpectralInterpolation, SfconvMomentumSpectralInterpolationInput,
    sfconv_interpolate_momentum_spectral_function,
};

use crate::error::{IoError, Result};

use super::types::{
    SfconvSpecfunctCompatibilityInput, SfconvSpecfunctData, SfconvSpecfunctSpectralRowsInput,
};
use super::validation::{
    validate_compatibility_input, validate_finite_scalar, validate_specfunct_dat,
    validate_specfunct_spectral_rows_input,
};

const SO2CONV_CACHE_SCALAR_TOLERANCE: f64 = 1.0e-12;

pub fn sfconv_specfunct_data_from_spectral_rows(
    input: SfconvSpecfunctSpectralRowsInput<'_>,
) -> Result<SfconvSpecfunctData> {
    validate_specfunct_spectral_rows_input(input)?;
    let data = SfconvSpecfunctData {
        wigner_seitz_radius: input.wigner_seitz_radius,
        core_hole_lifetime: input.core_hole_lifetime,
        asymmetric_phase: input.asymmetric_phase,
        satellite_type: input.satellite_type,
        low_q_mode: input.low_q_mode,
        pole_count: input.pole_count,
        pole_energy: input.pole_energy.to_owned(),
        pole_broadening: input.pole_broadening.to_owned(),
        pole_weight: input.pole_weight.to_owned(),
        spectral_info: input.spectral_info.to_owned(),
        weights: input.weights.to_owned(),
        extrinsic_quasiparticle: spectral_function_row(input.spectral_function, 0),
        extrinsic_satellite: spectral_function_row(input.spectral_function, 1),
        interference_quasiparticle: spectral_function_row(input.spectral_function, 2),
        interference_satellite: spectral_function_row(input.spectral_function, 3),
        intrinsic_satellite: spectral_function_row(input.spectral_function, 4),
        clipped_extrinsic_satellite: spectral_function_row(input.spectral_function, 7),
        energy_grid: input.energy_grid.to_owned(),
    };
    validate_specfunct_dat(&data)?;
    Ok(data)
}

/// Return whether parsed cache data matches the current SO2CONV inputs.
///
/// This mirrors the reuse checks in `SFCONV/so2conv.f90`: material scalars,
/// integer selectors, active pole rows, and the momentum grid must match. The
/// scalar comparison keeps a narrow tolerance because FEFF caches can carry
/// values generated before the rounded output header was written. FEFF compares
/// the momentum grid after converting both values to default `REAL`, so this
/// function compares those entries as `f32`.
pub fn sfconv_specfunct_matches_so2conv_inputs(
    data: &SfconvSpecfunctData,
    input: SfconvSpecfunctCompatibilityInput<'_>,
) -> Result<bool> {
    validate_specfunct_dat(data)?;
    validate_compatibility_input(input)?;

    if !cache_scalar_matches(data.wigner_seitz_radius, input.wigner_seitz_radius)
        || !cache_scalar_matches(data.core_hole_lifetime, input.core_hole_lifetime)
        || data.asymmetric_phase != input.asymmetric_phase
        || data.low_q_mode != input.low_q_mode
        || data.satellite_type != input.satellite_type
        || data.pole_count != input.pole_count
        || data.momentum_count() != input.momentum_grid.len()
    {
        return Ok(false);
    }

    let active_poles_match = (0..data.pole_count).all(|index| {
        cache_scalar_matches(data.pole_energy[index], input.pole_energy[index])
            && cache_scalar_matches(data.pole_broadening[index], input.pole_broadening[index])
            && cache_scalar_matches(data.pole_weight[index], input.pole_weight[index])
    });
    if !active_poles_match {
        return Ok(false);
    }

    let momentum_matches = (0..data.momentum_count()).all(|index| {
        (data.spectral_info[[index, 0]] as f32) == (input.momentum_grid[index] as f32)
    });
    Ok(momentum_matches)
}

fn cache_scalar_matches(left: f64, right: f64) -> bool {
    (left - right).abs() <= SO2CONV_CACHE_SCALAR_TOLERANCE * left.abs().max(right.abs()).max(1.0)
}

/// Build a validated core interpolation view over a `specfunct.dat` cache.
///
/// This maps the FEFF cache layout to `refeff-core`'s momentum spectral
/// interpolation input without copying the cached arrays. The first `sfinfo`
/// column is FEFF `pgrid`; columns 4 through 8 are `se`, `ce`, `width`, `z1`,
/// and `z1i`.
pub fn sfconv_specfunct_momentum_interpolation_input(
    data: &SfconvSpecfunctData,
    photoelectron_momentum: Real,
) -> Result<SfconvMomentumSpectralInterpolationInput<'_>> {
    validate_specfunct_dat(data)?;
    validate_finite_scalar(photoelectron_momentum, "photoelectron momentum")?;

    Ok(SfconvMomentumSpectralInterpolationInput {
        photoelectron_momentum,
        momentum_grid: data.spectral_info.column(0),
        energy_grid: data.energy_grid.view(),
        extrinsic_quasiparticle: data.extrinsic_quasiparticle.view(),
        extrinsic_satellite: data.extrinsic_satellite.view(),
        interference_quasiparticle: data.interference_quasiparticle.view(),
        interference_satellite: data.interference_satellite.view(),
        intrinsic_satellite: data.intrinsic_satellite.view(),
        clipped_extrinsic_satellite: data.clipped_extrinsic_satellite.view(),
        weights: data.weights.view(),
        self_energy_real: data.spectral_info.column(3),
        energy_correction: data.spectral_info.column(4),
        width: data.spectral_info.column(5),
        renormalization_real: data.spectral_info.column(6),
        renormalization_imag: data.spectral_info.column(7),
    })
}

/// Interpolate one cached `specfunct.dat` spectral row to a photoelectron momentum.
///
/// This is the typed handoff from FEFF's binary SO2CONV cache to the core
/// numerical interpolation kernel used by the future full driver.
pub fn sfconv_specfunct_interpolate_momentum(
    data: &SfconvSpecfunctData,
    photoelectron_momentum: Real,
) -> Result<SfconvMomentumSpectralInterpolation> {
    let input = sfconv_specfunct_momentum_interpolation_input(data, photoelectron_momentum)?;
    sfconv_interpolate_momentum_spectral_function(input)
        .map_err(|source| IoError::SpecfunctDatInterpolation { source })
}
fn spectral_function_row(values: ArrayView3<'_, f64>, source_row: usize) -> Array2<f64> {
    let (momentum_count, _, spectral_point_count) = values.dim();
    Array2::from_shape_fn(
        (momentum_count, spectral_point_count),
        |(momentum, point)| values[[momentum, source_row, point]],
    )
}