From bdf0c24a1cc7214df968d4b6560ec2ae8e132205 Mon Sep 17 00:00:00 2001 From: Rene Otten Date: Fri, 7 Aug 2026 22:01:22 -0700 Subject: [PATCH] fix(py): let local_density_of_states choose its broadening `local_density_of_states()` always passed `negf::DEFAULT_ETA` (1e-8 eV) to `green_function_sweep`, with no way to override it from Python. That default is ~5 orders of magnitude narrower than the default 1 meV energy grid (`DeviceParams::d_e`), so resonances land between grid points: sampling a 719 x 561 map at the default device leaves only ~0.3% of entries non-negligible. Plotted on a log scale the result looks empty, and the failure is silent -- no error, just a blank figure. It also means that after `solve_self_consistent()` (which defaults to `eta = 0.08`) the LDOS could not be evaluated at the same broadening the loop itself used. Add optional `eta` and `d_e` arguments. Defaults are unchanged, so existing no-argument callers keep the original MATLAB sharp-resonance behaviour. - eta: broadening in eV; rejects non-finite / non-positive values - d_e: per-call energy grid override in eV; defaults to DeviceParams::d_e - document the eta-vs-d_e pitfall and the (nE,) / (nE, n) return shapes Adds a regression test asserting that a larger eta resolves more of the spectrum on a fixed grid (verified to fail if eta is not actually plumbed through). --- crates/negforge-core/src/negf.rs | 54 ++++++++++++++++++++++++++++++++ crates/negforge-py/src/lib.rs | 42 ++++++++++++++++++++----- python/negforge/__init__.py | 25 ++++++++++++--- 3 files changed, 108 insertions(+), 13 deletions(-) diff --git a/crates/negforge-core/src/negf.rs b/crates/negforge-core/src/negf.rs index dbbd7a6..47205bc 100644 --- a/crates/negforge-core/src/negf.rs +++ b/crates/negforge-core/src/negf.rs @@ -370,6 +370,60 @@ mod tests { } } + /// Larger `eta` must produce a smoother, better-resolved spectrum on a + /// fixed energy grid. This is the property the Python `eta` argument on + /// `local_density_of_states()` exists to expose: at the historical + /// `DEFAULT_ETA` (1e-8 eV) the Lorentzian width is ~10^5 times narrower + /// than a typical 1 meV grid step, so resonances fall between grid points + /// and the sampled LDOS is a sparse set of spikes rather than a spectrum. + #[test] + fn larger_eta_resolves_more_of_the_spectrum_on_a_fixed_grid() { + let mut device = Device::new(DeviceParams { + a: 1.0, + l_ch: 10.0, + l_ds: 10.0, + auto_size_contacts: false, + v_g: 0.3, + ..Default::default() + }); + device.calc_potential(); + + let e_min = device.psi_f.iter().cloned().fold(f64::INFINITY, f64::min); + let energies: Vec = (0..400).map(|i| e_min + i as f64 * 0.001).collect(); + + // Fraction of (energy, site) samples carrying non-negligible weight. + let occupancy = |eta: f64| -> f64 { + let result = green_function_sweep(&device, &energies, eta); + let max = result + .g_diag + .iter() + .flat_map(|row| row.iter()) + .fold(0.0f64, |acc, v| acc.max(v.abs())); + assert!(max > 0.0, "sweep produced an all-zero LDOS"); + let total: usize = result.g_diag.iter().map(|row| row.len()).sum(); + let filled = result + .g_diag + .iter() + .flat_map(|row| row.iter()) + .filter(|v| v.abs() > max * 1e-6) + .count(); + filled as f64 / total as f64 + }; + + let sharp = occupancy(DEFAULT_ETA); + let broadened = occupancy(5e-3); + + assert!( + broadened > sharp, + "broadening should resolve more of the spectrum: \ + sharp(eta={DEFAULT_ETA:e})={sharp:.4}, broadened(eta=5e-3)={broadened:.4}" + ); + assert!( + broadened > 0.5, + "eta a few times the grid step should fill most of the map, got {broadened:.4}" + ); + } + /// Documents the corrected sign claim from the module docs' "Contact /// self-energies" section: for the branch of `ka_s`/`ka_d` this module /// actually computes (`ka in (0, pi)`, the reflected branch), the diff --git a/crates/negforge-py/src/lib.rs b/crates/negforge-py/src/lib.rs index 90e660b..cd24678 100644 --- a/crates/negforge-py/src/lib.rs +++ b/crates/negforge-py/src/lib.rs @@ -183,7 +183,38 @@ impl PyDevice { /// Run the NEGF sweep at the device's current potential and return /// `(energies, ldos)` where `ldos[k]` is the local density of states /// row (one value per grid site) at `energies[k]`. - fn local_density_of_states(&self) -> (Vec, Vec>) { + /// + /// `eta` is the imaginary broadening added to the energy argument of the + /// retarded Green's function, in eV. It defaults to + /// [`negforge_core::negf::DEFAULT_ETA`] (`1e-8`), which reproduces the + /// original MATLAB code's sharp-resonance behaviour. That value is far + /// smaller than the energy grid spacing `d_e` (1 meV by default), so + /// resonances are effectively sampled at random and the result is a set of + /// isolated spikes rather than a smooth spectrum — fine for reproducing the + /// original plot, but poor for visualisation or integration. Pass a value + /// of order a few times `d_e` (e.g. `5e-3`) for a smooth, resolvable LDOS + /// map. + /// + /// `d_e` overrides the device's energy grid spacing (eV) for this sweep + /// only; `None` uses `DeviceParams::d_e`. + #[pyo3(signature = (eta=None, d_e=None))] + fn local_density_of_states( + &self, + eta: Option, + d_e: Option, + ) -> PyResult<(Vec, Vec>)> { + let eta = eta.unwrap_or(negforge_core::negf::DEFAULT_ETA); + if !(eta > 0.0) || !eta.is_finite() { + return Err(pyo3::exceptions::PyValueError::new_err( + "eta must be a finite, strictly positive broadening in eV", + )); + } + let d_e = d_e.unwrap_or(self.inner.params.d_e); + if !(d_e > 0.0) || !d_e.is_finite() { + return Err(pyo3::exceptions::PyValueError::new_err( + "d_e must be a finite, strictly positive energy step in eV", + )); + } let e_min = self .inner .psi_f @@ -191,15 +222,10 @@ impl PyDevice { .cloned() .fold(f64::INFINITY, f64::min); let e_max = 0.7 * self.inner.e_max; - let d_e = self.inner.params.d_e; let steps = ((e_max - e_min) / d_e).floor().max(0.0) as usize; let energies: Vec = (0..=steps).map(|k| e_min + k as f64 * d_e).collect(); - let result = negforge_core::negf::green_function_sweep( - &self.inner, - &energies, - negforge_core::negf::DEFAULT_ETA, - ); - (result.energies, result.g_diag) + let result = negforge_core::negf::green_function_sweep(&self.inner, &energies, eta); + Ok((result.energies, result.g_diag)) } } diff --git a/python/negforge/__init__.py b/python/negforge/__init__.py index 41cb5fb..7f140a9 100644 --- a/python/negforge/__init__.py +++ b/python/negforge/__init__.py @@ -132,14 +132,29 @@ def screening_length(self) -> float: """Natural (screening) length lambda, nm.""" return self._inner.screening_length - def local_density_of_states(self): + def local_density_of_states(self, eta: float | None = None, d_e: float | None = None): """NEGF sweep at the device's current potential. - Returns `(energies, ldos)` where `ldos[k]` is the local density of - states across all grid sites at `energies[k]`. Matches - `calc_green()`'s `G_r_diag` in the original code. + Returns `(energies, ldos)` where `energies` has shape `(nE,)` (eV) and + `ldos` has shape `(nE, n)` — one row of grid-site values per energy. + Matches `calc_green()`'s `G_r_diag` in the original code. + + Parameters + ---------- + eta: + Imaginary broadening of the retarded Green's function, in eV. + Defaults to the Rust `negf::DEFAULT_ETA` (`1e-8`), reproducing the + original MATLAB behaviour. Because that is far below the energy + grid spacing `d_e` (1 meV by default), resonances land between grid + points and the result is a handful of isolated spikes — a log-scale + map made from it looks essentially empty. For a smooth, plottable + LDOS use a few times `d_e`, e.g. `eta=5e-3`. The self-consistent + loop separately defaults to `eta=0.08` eV for the same reason. + d_e: + Energy grid spacing for this sweep only, in eV. Defaults to + `DeviceParams::d_e` (1 meV). """ - energies, ldos = self._inner.local_density_of_states() + energies, ldos = self._inner.local_density_of_states(eta, d_e) return np.asarray(energies), np.asarray(ldos) def sweep_v_g(self, v_min: float, v_max: float, step: float, self_consistent: bool = False) -> IVCurve: