fix(py): let local_density_of_states choose its broadening

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<f64> = (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

crates/negforge-py/src/lib.rs

@@ -183,23 +183,49 @@ 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<f64>, Vec<Vec<f64>>) {
+    ///
+    /// `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<f64>,
+        d_e: Option<f64>,
+    ) -> PyResult<(Vec<f64>, Vec<Vec<f64>>)> {
+        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
             .iter()
             .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<f64> = (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))
     }
 }
 

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: