treams_core/lattice/
real.rs

1//! Real-space terms of the Ewald sum.
2//!
3//! Upstream: `_realsw*` and `_realcw*` of `lattice/_esum.pyx`. The bounded Kambe
4//! integrals at small splits are a treams-rs extension.
5
6use std::f64::consts::PI;
7
8use super::{
9    accuracy::Rounding,
10    wave::{Family, azimuthal, negative_order_sign, solid_jet},
11};
12use crate::{
13    Complex,
14    numerics::Jet,
15    special::{self, SmallSplitKambe, harmonic_normalization},
16};
17
18/// How the real-space terms take their Kambe integrals ([`kambe_jet`]).
19#[derive(Debug)]
20pub(super) enum RealKambe<'a, const N: usize> {
21    /// The integrals of the recurrences, without bounds.
22    Plain,
23    /// Below every automatic split ([`below_automatic`]): the integrals of
24    /// [`SmallSplitKambe`], which keep their closed forms below the rounding `largest`
25    /// of the largest term so far in each component (see [`kambe_floor`]), storing the
26    /// bounds [`cancellation`] reads in `bound`.
27    ///
28    /// [`below_automatic`]: super::accuracy::below_automatic
29    /// [`cancellation`]: super::accuracy::cancellation
30    SmallSplit {
31        kambe: &'a SmallSplitKambe,
32        largest: &'a Rounding<N>,
33        bound: &'a mut Rounding<N>,
34    },
35    /// Only the `eta` derivative, in slot 0, for the split derivative.
36    EtaDerivative,
37}
38
39/// Kambe integral of order `n` at `z` of a real-space term with the factor `factor`, with
40/// its `z` derivative `-z I_(n+2)`, or its `eta` derivative alone.
41fn kambe_jet<const N: usize>(
42    n: i32,
43    z: Jet<N>,
44    eta: Complex,
45    mode: RealKambe<'_, N>,
46    factor: &Jet<N>,
47) -> Jet<N> {
48    let (value, next) = match mode {
49        RealKambe::EtaDerivative => {
50            let derivative =
51                -eta.powi(n) * (0.5 * (1.0 / (eta * eta) - z.value * z.value * eta * eta)).exp();
52            let mut jet = Jet::default();
53            if let Some(slot) = jet.derivative.first_mut() {
54                *slot = derivative;
55            }
56            return jet;
57        }
58        RealKambe::Plain => {
59            let next = (N > 0).then(|| special::kambe(n + 2, z.value, eta));
60            (special::kambe(n, z.value, eta), next)
61        }
62        RealKambe::SmallSplit {
63            kambe,
64            largest,
65            bound,
66        } => small_split_kambe(n, &z, kambe, largest, bound, factor),
67    };
68    let derivative = next.map_or_else(Complex::default, |next| -z.value * next);
69    z.chain(value, derivative)
70}
71
72/// The radial integral shared by all spherical orders of one degree at a point.
73/// Used only above the small-split cancellation regime, where the integral does not
74/// depend on the harmonic's rounding budget.
75pub(super) fn spherical_radial<const N: usize>(
76    l: i32,
77    k: Jet<N>,
78    r: [Jet<N>; 3],
79    eta: Complex,
80) -> Jet<N> {
81    let radius = r.into_iter().map(|r| r * r).sum::<Jet<N>>().sqrt();
82    kambe_jet(2 * l, k * radius, eta, RealKambe::Plain, &Jet::default())
83}
84
85/// A spherical real-space term from its already evaluated radial integral.
86pub(super) fn spherical_from_radial<const N: usize>(
87    l: i32,
88    m: i32,
89    k: Jet<N>,
90    r: [Jet<N>; 3],
91    radial: Jet<N>,
92) -> Jet<N> {
93    let harmonic = solid_jet(l, m, r) * harmonic_normalization(l, m);
94    let factor = -Complex::i() * (2.0 / PI).sqrt() * k.powi(l) * harmonic;
95    factor * radial
96}
97
98/// `I_n` and, for jets, `I_(n+2)` at `z` in [`RealKambe::SmallSplit`], storing their
99/// bounds in the term with the factor `factor` in `bound`.
100///
101/// Kept out of line, like [`Shells::add_compensated`], so the term loop at the other
102/// splits stays small. Calibrated: inlined, it slows 3D jets at the automatic split by
103/// 1% to 2.7%.
104///
105/// [`Shells::add_compensated`]: super::shells::Shells::add_compensated
106#[inline(never)]
107pub(super) fn small_split_kambe<const N: usize>(
108    n: i32,
109    z: &Jet<N>,
110    kambe: &SmallSplitKambe,
111    largest: &Rounding<N>,
112    bound: &mut Rounding<N>,
113    factor: &Jet<N>,
114) -> (Complex, Option<Complex>) {
115    let floor = kambe_floor(largest, factor);
116    let (value, value_bound, next) = if N == 0 {
117        let [(value, value_bound)] = kambe.get([n], z.value, floor);
118        (value, value_bound, None)
119    } else {
120        let [(value, value_bound), next] = kambe.get([n, n + 2], z.value, floor);
121        (value, value_bound, Some(next))
122    };
123    let (next_bound, modulus) = (next.map_or(0.0, |(_, bound)| bound), z.value.l1_norm());
124    *bound = Rounding {
125        value: value_bound,
126        derivative: z.derivative.map(|d| next_bound * modulus * d.l1_norm()),
127    }
128    .times(factor);
129    (value, next.map(|(next, _)| next))
130}
131
132/// The error, in units of epsilon, up to which [`SmallSplitKambe::get`] keeps the closed
133/// forms of the Kambe integral of a real-space term with the factor `factor`: what leaves
134/// less in each component of the term than the rounding `largest` of the largest term so
135/// far, carried by the value of the factor or, where that vanishes (off-axis orders on a
136/// 1D spherical axis, odd `l + m` in a 2D plane), by its derivatives.
137fn kambe_floor<const N: usize>(largest: &Rounding<N>, factor: &Jet<N>) -> f64 {
138    let scale = factor.value.l1_norm();
139    if scale > 0.0 {
140        return largest.value / scale;
141    }
142    largest
143        .derivative
144        .iter()
145        .zip(&factor.derivative)
146        .filter(|&(_, &d)| d != Complex::default())
147        .map(|(largest, d)| largest / d.l1_norm())
148        .fold(f64::INFINITY, f64::min)
149}
150
151/// One real-space summand at the displacement `r = -shift - R`.
152///
153/// # Formula
154///
155/// Spherical: `-i sqrt(2 / pi) k^l |r|^l Y_lm(r / |r|) I_(2l)(k |r|, eta)`, where
156/// `|r|^l Y_lm` is `harmonic_normalization(l, m)` times the solid harmonic. Cylindrical:
157/// `-2 i / pi (-1)^m (k (x ± i y))^|m| I_(2|m| - 1)(k |r|, eta)`, with the sign `(-1)^m`
158/// and `x - i y` for negative `m` only.
159///
160/// [`RealKambe::SmallSplit`] stores first-order bounds on the absolute rounding errors
161/// that its Kambe integral leaves in its value and derivatives (see [`kambe_jet`]).
162pub(super) fn real_term<const N: usize>(
163    wave: Family,
164    k: Jet<N>,
165    r: [Jet<N>; 3],
166    eta: Complex,
167    mode: RealKambe<'_, N>,
168) -> Jet<N> {
169    let radius = r.into_iter().map(|r| r * r).sum::<Jet<N>>().sqrt();
170    let (factor, order) = match wave {
171        Family::Spherical { l, m } => {
172            let harmonic = solid_jet(l, m, r) * harmonic_normalization(l, m);
173            (
174                -Complex::i() * (2.0 / PI).sqrt() * k.powi(l) * harmonic,
175                2 * l,
176            )
177        }
178        Family::Cylindrical { m } => {
179            let (xy, sign) = (azimuthal(r[0], r[1], m), negative_order_sign(m));
180            (
181                -2.0 * Complex::i() / PI * sign * (k * xy).powi(m.abs()),
182                2 * m.abs() - 1,
183            )
184        }
185    };
186    factor * kambe_jet(order, k * radius, eta, mode, &factor)
187}
treams_core/lattice/real.rs

treams_core/lattice/
real.rs
