lindhard/electron/elastic/corrections.rs
1//! Optional corrections to the static potential: local exchange and
2//! correlation-polarization (issue #91, step 3 of the Mott elastic model).
3//!
4//! Both are **off by default**. They need the atomic electron density
5//! `rho(r)`, which a [`ScreenedPotential`] does not provide, so the density is a
6//! separate, explicit input ([`ElectronDensity`]). For the analytic potentials
7//! of this module the density is the one that generates the potential through
8//! Poisson's equation (Salvat, Martinez, Mayol & Parellada, Phys. Rev. A 36, 467
9//! (1987), doi:10.1103/PhysRevA.36.467, Eqs. (2) and (12)); it is never
10//! inferred from an arbitrary user potential.
11//!
12//! # Sources
13//!
14//! The model is the one of F. Salvat, Phys. Rev. A 68, 012708 (2003),
15//! doi:10.1103/PhysRevA.68.012708, Sec. II (read in the open copy of the
16//! University of Barcelona repository, diposit.ub.edu). That paper is the
17//! published account of the optical-model potential that the ELSEPA paper
18//! (Salvat, Jablonski & Powell, Comput. Phys. Commun. 165, 157 (2005)) uses;
19//! the ELSEPA paper itself is closed access and was not opened, and no ELSEPA
20//! code was read (Tier C).
21//!
22//! * **Exchange**, Furness & McCarthy, J. Phys. B 6, 2280 (1973),
23//! doi:10.1088/0022-3700/6/11/021, as written in Salvat (2003) Eq. (2) (and
24//! identically in Jablonski, Salvat & Powell, J. Phys. Chem. Ref. Data 33,
25//! 409 (2004), Eq. (15)), in hartree atomic units:
26//!
27//! ```text
28//! V_ex(r) = 1/2 [E - V_st(r)] - 1/2 { [E - V_st(r)]^2 + 4 pi rho(r) }^(1/2)
29//! ```
30//!
31//! with `E` the kinetic energy of the projectile and `V_st` the static
32//! potential energy. **The equation number in the 1973 paper itself is not
33//! cited**: that paper is closed access (IOP subscription; Unpaywall lists
34//! only a metadata record), so only the two secondary accounts above were
35//! read. We evaluate the algebraically identical form
36//! `-2 pi rho / (D + sqrt(D^2 + 4 pi rho))`, `D = E - V_st >= 0`, which avoids
37//! the cancellation of the printed form when `4 pi rho << D^2` (high energy,
38//! outer atom). `V_ex <= 0` everywhere and `V_ex = 0` where `rho = 0`.
39//! Salvat (2003) quotes Bransden et al. (1976) that it is accurate above
40//! about 1 hartree.
41//! * **Correlation-polarization**, Salvat (2003) Eqs. (3)-(7) and (9):
42//! - long-range polarization (Buckingham) `V_pol = -alpha_p/(2 (r^2 + d^2)^2)`,
43//! Eq. (3), with `d^4 = (1/2) alpha_p Z^(-1/3) b_pol^2` (bohr units),
44//! Eq. (4), and Seltzer's empirical `b_pol^2 = (E - 50 eV)/(16 eV)`,
45//! Eq. (5) (proposed for `E >= 100 eV`);
46//! - the local-density correlation potential of Perdew & Zunger, Phys. Rev.
47//! B 23, 5048 (1981), in the parametrization printed as Salvat (2003)
48//! Eqs. (6)-(7), with `r_s = (3/(4 pi rho))^(1/3)`:
49//! `V_co = 0.0311 ln r_s - 0.0584 + 0.00133 r_s ln r_s - 0.0084 r_s`
50//! for `r_s < 1` and
51//! `V_co = -beta_0 (1 + 7/6 beta_1 r_s^(1/2) + 4/3 beta_2 r_s)/(1 + beta_1 r_s^(1/2) + beta_2 r_s)^2`
52//! for `r_s >= 1`, `beta_0 = 0.1423`, `beta_1 = 1.0529`, `beta_2 = 0.3334`.
53//! **Sign of Eq. (7a).** Salvat prints Eq. (7a) with an overall minus
54//! sign in front of the bracket. With that sign the two branches jump from
55//! `+0.0668` to `-0.0668` hartree at `r_s = 1` and the high-density
56//! correlation potential would be repulsive. We use the sign that makes
57//! the branches continuous (tested), which is also the sign of the same
58//! Perdew-Zunger expression `A ln r_s + (B - A/3) + (2/3) C r_s ln r_s +
59//! (1/3)(2D - C) r_s` with `A = 0.0311`, `B = -0.048`, `C = 0.002`,
60//! `D = -0.0116` as reproduced (citing Perdew & Zunger) on the Delta
61//! Science Institute DFT course page
62//! <https://www.dsedu.org/courses/dft/xc>, read 2026-10-07. The 1981 paper
63//! itself was not opened (APS bot challenge; Unpaywall's link is the
64//! publisher PDF);
65//! - the join, Eq. (9): `V_cp = max(V_co, V_pol)` for `r < r_cp` and
66//! `V_cp = V_pol` for `r >= r_cp`, where `r_cp` is the outer radius at
67//! which `V_co` and `V_pol` cross.
68//!
69//! # Numerical choices (ours, not from the papers)
70//!
71//! * **Finding `r_cp`.** `V_co - V_pol` is sampled on a geometric grid (ratio
72//! 1.01) from [`JOIN_SCAN_FLOOR`] to the outer radius and the outermost sign
73//! change is refined by bisection. Two crossings closer than one grid ratio
74//! would be missed. If `V_co > V_pol` on the whole grid the crossing lies
75//! below the floor (`V_co` diverges only logarithmically at the nucleus);
76//! then Eq. (9) is `V_pol` on the whole grid and [`r_cp`] is reported as
77//! `None`. If `V_co < V_pol` at the outer radius, the outer crossing lies
78//! beyond it and construction fails.
79//! * **Outer radius.** The `r^-4` tail is integrated out to
80//! [`CorrelationPolarization::outer_radius`] and set to zero beyond it (a
81//! breakpoint and the matching radius). The neglected tail is not covered by
82//! the `|r V| < threshold` rule of [`ScreenedPotential::matching_radius`],
83//! which would put the matching radius at thousands of bohr; its effect is
84//! measured by the convergence tests instead (doubling the outer radius).
85//! The matching radius of the corrected potential is the larger of the outer
86//! radius and the default scan applied to `V_st + V_ex` (both `<= 0`, so no
87//! cancellation can stop the scan early).
88//! * **Start radius.** With the polarization tail, partial waves of `l` up to
89//! about `k` times the outer radius have their turning points where the
90//! potential is not negligible. The step-1 start rule is not accurate there
91//! (measured: errors up to 0.25 rad at `l ~ 1250`, Cu, 10 keV, 100 bohr), so
92//! the solver uses a WKB start criterion instead. Since #131 it does so for
93//! every potential, corrected or not.
94//!
95//! [`r_cp`]: CorrelationPolarizationInfo::join_radius_bohr
96
97use super::{
98 invalid, solve, ElasticError, ElasticResult, ElasticSolver, SalvatDhfs, ScreenedPotential,
99 SolverOptions, Yukawa, HARTREE_EV,
100};
101use std::f64::consts::PI;
102
103/// Model identifier recorded in the metadata when exchange is on.
104pub const EXCHANGE_MODEL: &str =
105 "Furness-McCarthy local exchange, J. Phys. B 6, 2280 (1973), as Salvat, PRA 68, 012708 (2003) Eq. (2)";
106
107/// Model identifier recorded in the metadata when correlation-polarization is on.
108pub const CORRELATION_POLARIZATION_MODEL: &str =
109 "Salvat, PRA 68, 012708 (2003) Eqs. (3)-(7), (9): \
110 Perdew-Zunger LDA correlation joined to a Buckingham polarization tail";
111
112/// Lower end of the grid on which the join radius `r_cp` is searched, bohr.
113pub const JOIN_SCAN_FLOOR: f64 = 1.0e-8;
114
115/// Default outer radius of the polarization tail, bohr (see the module docs).
116pub const DEFAULT_OUTER_RADIUS: f64 = 50.0;
117
118/// Atomic electron density, the input the corrections need beyond the
119/// potential.
120pub trait ElectronDensity: Sync {
121 /// Electron number density at radius `r` bohr, electrons per bohr^3
122 /// (`>= 0`).
123 fn density(&self, r: f64) -> f64;
124
125 /// Where the density comes from (recorded in the result metadata).
126 fn density_source(&self) -> String;
127}
128
129/// Poisson density of the analytic screening function, Salvat et al. (1987)
130/// Eq. (12): `rho(r) = (Z/(4 pi r)) sum_i A_i alpha_i^2 exp(-alpha_i r)`,
131/// summed over the terms in use (two for the asterisked, `A_3 = 0` rows of
132/// Table I).
133impl ElectronDensity for SalvatDhfs {
134 fn density(&self, r: f64) -> f64 {
135 let s: f64 = self
136 .amplitudes()
137 .iter()
138 .zip(self.alphas())
139 .map(|(a, al)| a * al * al * (-al * r).exp())
140 .sum();
141 self.z * s / (4.0 * PI * r)
142 }
143 fn density_source(&self) -> String {
144 format!(
145 "Poisson density of the Salvat et al. (1987) analytic screening function, \
146 PRA 36, 467 Eq. (12), Z={}, A={:?}, alpha={:?} 1/bohr",
147 self.z,
148 self.amplitudes(),
149 self.alphas()
150 )
151 }
152}
153
154/// Poisson density of `V = -Z exp(-r/a)/r` (the one-term case of Salvat et
155/// al. (1987) Eq. (12)): `rho = Z exp(-r/a)/(4 pi a^2 r)`. A fixture density.
156impl ElectronDensity for Yukawa {
157 fn density(&self, r: f64) -> f64 {
158 let a = self.screening_length;
159 self.z * (-r / a).exp() / (4.0 * PI * a * a * r)
160 }
161 fn density_source(&self) -> String {
162 format!(
163 "Poisson density of a Yukawa potential (fixture), Z={}, a={} bohr",
164 self.z, self.screening_length
165 )
166 }
167}
168
169/// How the cutoff `d` of the polarization potential is set (Salvat 2003
170/// Eqs. (4)-(5)).
171#[derive(Debug, Clone, Copy, PartialEq)]
172pub enum PolarizationCutoff {
173 /// Seltzer's recipe `b_pol^2 = (E - 50 eV)/(16 eV)`, Eq. (5); needs
174 /// `E > 50 eV` (proposed for `E >= 100 eV`).
175 Seltzer,
176 /// A caller-chosen `b_pol^2 >= 0` (dimensionless) in Eq. (4).
177 BPolSquared(f64),
178}
179
180/// Inputs of the correlation-polarization correction.
181#[derive(Debug, Clone, PartialEq)]
182pub struct CorrelationPolarization {
183 /// Static dipole polarizability `alpha_p` of the target atom, bohr^3
184 /// (atomic units of polarizability).
185 pub polarizability: f64,
186 /// Citation of `polarizability` (recorded in the metadata; must not be
187 /// empty).
188 pub polarizability_source: String,
189 /// Cutoff parameter of the Buckingham potential.
190 pub cutoff: PolarizationCutoff,
191 /// Radius beyond which the polarization tail is dropped, bohr (a numerical
192 /// setting, see the module docs).
193 pub outer_radius: f64,
194}
195
196impl CorrelationPolarization {
197 /// Seltzer's cutoff and [`DEFAULT_OUTER_RADIUS`].
198 pub fn new(polarizability: f64, polarizability_source: impl Into<String>) -> Self {
199 Self {
200 polarizability,
201 polarizability_source: polarizability_source.into(),
202 cutoff: PolarizationCutoff::Seltzer,
203 outer_radius: DEFAULT_OUTER_RADIUS,
204 }
205 }
206}
207
208/// Which corrections to add to the static potential. The default is none.
209#[derive(Debug, Clone, PartialEq, Default)]
210pub struct Corrections {
211 /// Furness-McCarthy exchange on/off.
212 pub exchange: bool,
213 /// Correlation-polarization, `None` for off.
214 pub correlation_polarization: Option<CorrelationPolarization>,
215}
216
217impl Corrections {
218 /// Both corrections off.
219 pub fn none() -> Self {
220 Self::default()
221 }
222
223 /// True when neither correction is on.
224 pub fn is_none(&self) -> bool {
225 !self.exchange && self.correlation_polarization.is_none()
226 }
227}
228
229/// Metadata of the exchange correction.
230#[derive(Debug, Clone, PartialEq)]
231pub struct ExchangeInfo {
232 /// Model identifier ([`EXCHANGE_MODEL`]).
233 pub model: &'static str,
234}
235
236/// Metadata of the correlation-polarization correction.
237#[derive(Debug, Clone, PartialEq)]
238pub struct CorrelationPolarizationInfo {
239 /// Model identifier ([`CORRELATION_POLARIZATION_MODEL`]).
240 pub model: &'static str,
241 /// Dipole polarizability used, bohr^3.
242 pub polarizability_bohr3: f64,
243 /// Its citation, as supplied.
244 pub polarizability_source: String,
245 /// `b_pol^2` used in Eq. (4), dimensionless.
246 pub b_pol_squared: f64,
247 /// Cutoff `d` of Eq. (3), bohr.
248 pub cutoff_d_bohr: f64,
249 /// Outer crossing `r_cp` of Eq. (9), bohr; `None` when it lies below
250 /// [`JOIN_SCAN_FLOOR`] (then `V_cp = V_pol` on the whole grid).
251 pub join_radius_bohr: Option<f64>,
252 /// Radius beyond which the tail is dropped, bohr.
253 pub outer_radius_bohr: f64,
254}
255
256/// Which corrections a result includes. The default reports both off.
257#[derive(Debug, Clone, PartialEq, Default)]
258pub struct CorrectionMetadata {
259 /// `Some` when exchange is on.
260 pub exchange: Option<ExchangeInfo>,
261 /// `Some` when correlation-polarization is on.
262 pub correlation_polarization: Option<CorrelationPolarizationInfo>,
263 /// [`ElectronDensity::density_source`] of the density used, when any
264 /// correction is on.
265 pub density_source: Option<String>,
266}
267
268#[derive(Debug, Clone, Copy)]
269struct CpState {
270 alpha: f64,
271 d2: f64,
272 r_cp: f64,
273 r_out: f64,
274}
275
276/// The static potential plus the selected corrections, **at one kinetic
277/// energy** (the exchange and the polarization cutoff depend on it). Build the
278/// solver with [`CorrectedPotential::solver`] so the solver energy cannot
279/// differ from the one the potential was built for.
280pub struct CorrectedPotential<'a> {
281 stat: &'a dyn ScreenedPotential,
282 dens: &'a dyn ElectronDensity,
283 energy_ev: f64,
284 e: f64,
285 exchange: bool,
286 cp: Option<CpState>,
287 meta: CorrectionMetadata,
288}
289
290/// Perdew-Zunger correlation potential, Salvat (2003) Eqs. (6)-(7) (sign of
291/// (7a) as discussed in the module docs), hartree, for density `rho`
292/// (bohr^-3). Zero for `rho <= 0` (the `r_s -> infinity` limit).
293pub fn correlation_potential(rho: f64) -> f64 {
294 if rho.is_nan() || rho <= 0.0 {
295 return 0.0;
296 }
297 let rs = (3.0 / (4.0 * PI * rho)).cbrt();
298 if rs < 1.0 {
299 let l = rs.ln();
300 0.0311 * l - 0.0584 + 0.00133 * rs * l - 0.0084 * rs
301 } else {
302 let (b0, b1, b2) = (0.1423, 1.0529, 0.3334);
303 let sq = rs.sqrt();
304 let den = 1.0 + b1 * sq + b2 * rs;
305 -b0 * (1.0 + 7.0 / 6.0 * b1 * sq + 4.0 / 3.0 * b2 * rs) / (den * den)
306 }
307}
308
309/// Furness-McCarthy exchange potential, Salvat (2003) Eq. (2), hartree, for
310/// local kinetic-energy argument `d = E - V_st` (hartree) and density `rho`
311/// (bohr^-3), in the cancellation-free form of the module docs.
312pub fn exchange_potential(d: f64, rho: f64) -> f64 {
313 let rho = rho.max(0.0);
314 let s = (d * d + 4.0 * PI * rho).sqrt();
315 if d >= 0.0 {
316 if rho == 0.0 {
317 return 0.0;
318 }
319 -2.0 * PI * rho / (d + s)
320 } else {
321 0.5 * (d - s)
322 }
323}
324
325// Manual: the borrowed potential and density are trait objects without `Debug`.
326impl std::fmt::Debug for CorrectedPotential<'_> {
327 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
328 f.debug_struct("CorrectedPotential")
329 .field("energy_ev", &self.energy_ev)
330 .field("exchange", &self.exchange)
331 .field("cp", &self.cp)
332 .field("meta", &self.meta)
333 .finish_non_exhaustive()
334 }
335}
336
337impl<'a> CorrectedPotential<'a> {
338 /// Build the corrected potential of `stat` with density `dens` at kinetic
339 /// energy `energy_ev` (eV, > 0).
340 pub fn new(
341 stat: &'a dyn ScreenedPotential,
342 dens: &'a dyn ElectronDensity,
343 energy_ev: f64,
344 corrections: &Corrections,
345 ) -> Result<Self, ElasticError> {
346 if !energy_ev.is_finite() || energy_ev <= 0.0 {
347 return Err(invalid(
348 "energy",
349 format!("{energy_ev} eV (need finite, > 0)"),
350 ));
351 }
352 let e = energy_ev / HARTREE_EV;
353 let mut meta = CorrectionMetadata::default();
354 if corrections.exchange {
355 meta.exchange = Some(ExchangeInfo {
356 model: EXCHANGE_MODEL,
357 });
358 }
359 let mut this = Self {
360 stat,
361 dens,
362 energy_ev,
363 e,
364 exchange: corrections.exchange,
365 cp: None,
366 meta,
367 };
368 if let Some(cp) = &corrections.correlation_polarization {
369 let (state, info) = this.build_cp(cp)?;
370 this.cp = Some(state);
371 this.meta.correlation_polarization = Some(info);
372 }
373 if !corrections.is_none() {
374 this.meta.density_source = Some(dens.density_source());
375 }
376 Ok(this)
377 }
378
379 fn build_cp(
380 &self,
381 cp: &CorrelationPolarization,
382 ) -> Result<(CpState, CorrelationPolarizationInfo), ElasticError> {
383 let alpha = cp.polarizability;
384 if !alpha.is_finite() || alpha <= 0.0 {
385 return Err(invalid(
386 "dipole polarizability",
387 format!("{alpha} bohr^3 (need finite, > 0)"),
388 ));
389 }
390 if cp.polarizability_source.trim().is_empty() {
391 return Err(invalid(
392 "polarizability source",
393 "a citation for the polarizability is required",
394 ));
395 }
396 if !cp.outer_radius.is_finite() || cp.outer_radius <= 0.0 {
397 return Err(invalid(
398 "outer radius",
399 format!("{} bohr (need finite, > 0)", cp.outer_radius),
400 ));
401 }
402 let b2 = match cp.cutoff {
403 PolarizationCutoff::Seltzer => {
404 if self.energy_ev <= 50.0 {
405 return Err(invalid(
406 "polarization cutoff",
407 format!(
408 "Seltzer's b_pol^2 = (E - 50 eV)/16 eV needs E > 50 eV, got {} eV",
409 self.energy_ev
410 ),
411 ));
412 }
413 (self.energy_ev - 50.0) / 16.0
414 }
415 PolarizationCutoff::BPolSquared(b2) => {
416 if !b2.is_finite() || b2 < 0.0 {
417 return Err(invalid("b_pol^2", format!("{b2} (need finite, >= 0)")));
418 }
419 b2
420 }
421 };
422 let z = self.stat.nuclear_charge();
423 if !(z.is_finite() && z > 0.0) {
424 return Err(invalid(
425 "nuclear charge",
426 "the polarization cutoff of Eq. (4) needs the atomic number Z > 0",
427 ));
428 }
429 let d4 = 0.5 * alpha * z.powf(-1.0 / 3.0) * b2;
430 let d2 = d4.sqrt();
431 let r_out = cp.outer_radius;
432 let vpol = |r: f64| -alpha / (2.0 * (r * r + d2).powi(2));
433 let f = |r: f64| correlation_potential(self.dens.density(r)) - vpol(r);
434 if f(r_out) < 0.0 {
435 return Err(invalid(
436 "outer radius",
437 format!(
438 "{r_out} bohr: the correlation potential is still below the polarization \
439 tail there, so the outer crossing r_cp of Eq. (9) lies beyond it"
440 ),
441 ));
442 }
443 // Outermost sign change of f on a geometric grid.
444 let mut grid = Vec::new();
445 let mut r = JOIN_SCAN_FLOOR;
446 while r < r_out {
447 grid.push(r);
448 r *= 1.01;
449 }
450 grid.push(r_out);
451 let vals: Vec<f64> = grid.iter().map(|&r| f(r)).collect();
452 let mut r_cp = None;
453 for i in (1..grid.len()).rev() {
454 if (vals[i - 1] < 0.0) != (vals[i] < 0.0) {
455 let (mut lo, mut hi) = (grid[i - 1], grid[i]);
456 let neg_lo = vals[i - 1] < 0.0;
457 for _ in 0..200 {
458 let mid = 0.5 * (lo + hi);
459 if (f(mid) < 0.0) == neg_lo {
460 lo = mid;
461 } else {
462 hi = mid;
463 }
464 if hi - lo <= 1e-15 * hi {
465 break;
466 }
467 }
468 r_cp = Some(0.5 * (lo + hi));
469 break;
470 }
471 }
472 let state = CpState {
473 alpha,
474 d2,
475 r_cp: r_cp.unwrap_or(JOIN_SCAN_FLOOR),
476 r_out,
477 };
478 let info = CorrelationPolarizationInfo {
479 model: CORRELATION_POLARIZATION_MODEL,
480 polarizability_bohr3: alpha,
481 polarizability_source: cp.polarizability_source.clone(),
482 b_pol_squared: b2,
483 cutoff_d_bohr: d2.sqrt(),
484 join_radius_bohr: r_cp,
485 outer_radius_bohr: r_out,
486 };
487 Ok((state, info))
488 }
489
490 /// Kinetic energy the potential was built for, eV.
491 pub fn energy_ev(&self) -> f64 {
492 self.energy_ev
493 }
494
495 /// What this potential includes.
496 pub fn metadata(&self) -> &CorrectionMetadata {
497 &self.meta
498 }
499
500 /// Exchange potential energy at `r` (0 when exchange is off), hartree.
501 pub fn exchange_energy(&self, r: f64) -> f64 {
502 if !self.exchange {
503 return 0.0;
504 }
505 exchange_potential(self.e - self.stat.energy(r), self.dens.density(r))
506 }
507
508 /// Polarization potential `V_pol(r)`, Eq. (3), without the join or the
509 /// outer cut (0 when the correction is off), hartree.
510 pub fn polarization_energy(&self, r: f64) -> f64 {
511 match self.cp {
512 Some(s) => -s.alpha / (2.0 * (r * r + s.d2).powi(2)),
513 None => 0.0,
514 }
515 }
516
517 /// Correlation-polarization potential `V_cp(r)`, Eq. (9), with the outer
518 /// cut (0 when off), hartree.
519 pub fn correlation_polarization_energy(&self, r: f64) -> f64 {
520 let Some(s) = self.cp else {
521 return 0.0;
522 };
523 if r >= s.r_out {
524 return 0.0;
525 }
526 let vpol = -s.alpha / (2.0 * (r * r + s.d2).powi(2));
527 if r < s.r_cp {
528 vpol.max(correlation_potential(self.dens.density(r)))
529 } else {
530 vpol
531 }
532 }
533
534 /// A solver at the energy this potential was built for.
535 pub fn solver(&self, opts: SolverOptions) -> Result<ElasticSolver<'_>, ElasticError> {
536 ElasticSolver::new(self, self.energy_ev, opts)
537 }
538}
539
540impl ScreenedPotential for CorrectedPotential<'_> {
541 fn energy(&self, r: f64) -> f64 {
542 let v = self.stat.energy(r);
543 if !self.exchange && self.cp.is_none() {
544 return v;
545 }
546 let mut total = v;
547 if self.exchange {
548 total += exchange_potential(self.e - v, self.dens.density(r));
549 }
550 if self.cp.is_some() {
551 total += self.correlation_polarization_energy(r);
552 }
553 total
554 }
555 fn nuclear_charge(&self) -> f64 {
556 self.stat.nuclear_charge()
557 }
558 fn breakpoints(&self) -> Vec<f64> {
559 let mut b = self.stat.breakpoints();
560 if let Some(s) = self.cp {
561 b.push(s.r_out);
562 }
563 b
564 }
565 fn length_scale(&self) -> f64 {
566 self.stat.length_scale()
567 }
568 fn bound_energy_ev(&self) -> Option<f64> {
569 Some(self.energy_ev)
570 }
571 fn correction_metadata(&self) -> CorrectionMetadata {
572 self.meta.clone()
573 }
574 fn matching_radius(&self, threshold: f64) -> f64 {
575 if !self.exchange && self.cp.is_none() {
576 return self.stat.matching_radius(threshold);
577 }
578 // The default rule on V_st + V_ex (both <= 0), beyond every breakpoint
579 // including the outer radius of the polarization tail.
580 let mut r = self.length_scale();
581 let floor = self.breakpoints().into_iter().fold(0.0_f64, f64::max);
582 let short = |r: f64| self.stat.energy(r) + self.exchange_energy(r);
583 while (r < floor || (r * short(r)).abs() >= threshold) && r < 1.0e7 {
584 r *= 1.05;
585 }
586 r
587 }
588}
589
590/// [`solve`] with optional corrections. With both off this **is** [`solve`]
591/// on `pot` (bit-identical results, metadata reporting both off); otherwise
592/// the corrected potential is built at `energy_ev` and solved at the same
593/// energy, and the result's [`ElasticResult::corrections`] records the
594/// choices.
595pub fn solve_corrected(
596 pot: &dyn ScreenedPotential,
597 density: &dyn ElectronDensity,
598 energy_ev: f64,
599 corrections: &Corrections,
600 thetas: &[f64],
601 opts: SolverOptions,
602) -> Result<ElasticResult, ElasticError> {
603 if corrections.is_none() {
604 return solve(pot, energy_ev, thetas, opts);
605 }
606 let corrected = CorrectedPotential::new(pot, density, energy_ev, corrections)?;
607 solve(&corrected, energy_ev, thetas, opts)
608}