diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..3220cb9
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,21 @@
+/target/
+
+# Python
+__pycache__/
+*.pyc
+.venv/
+venv/
+*.egg-info/
+build/
+dist/
+*.whl
+*.so
+*.pyd
+*.dylib
+
+# Jupyter
+.ipynb_checkpoints/
+
+# Editors / OS
+.DS_Store
+*.swp
diff --git a/Cargo.lock b/Cargo.lock
new file mode 100644
index 0000000..72a6fda
--- /dev/null
+++ b/Cargo.lock
@@ -0,0 +1,237 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "approx"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "autocfg"
+version = "1.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53"
+
+[[package]]
+name = "cfg-if"
+version = "1.0.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801"
+
+[[package]]
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
+
+[[package]]
+name = "indoc"
+version = "2.0.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "79cf5c93f93228cf8efb3ba362535fb11199ac548a09ce117c9b1adc3030d706"
+dependencies = [
+ "rustversion",
+]
+
+[[package]]
+name = "libc"
+version = "0.2.189"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2"
+
+[[package]]
+name = "memoffset"
+version = "0.9.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "488016bfae457b036d996092f6cb448677611ce4449e970ceaf42695203f218a"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "negforge-core"
+version = "0.1.0"
+dependencies = [
+ "approx",
+ "num-complex",
+ "thiserror",
+]
+
+[[package]]
+name = "negforge-py"
+version = "0.1.0"
+dependencies = [
+ "negforge-core",
+ "pyo3",
+]
+
+[[package]]
+name = "num-complex"
+version = "0.4.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495"
+dependencies = [
+ "num-traits",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "once_cell"
+version = "1.21.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50"
+
+[[package]]
+name = "portable-atomic"
+version = "1.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3d20d5497ef88037a52ff98267d066e7f11fcc5e99bbfbd58a42336193aacec3"
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.107"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "pyo3"
+version = "0.22.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f402062616ab18202ae8319da13fa4279883a2b8a9d9f83f20dbade813ce1884"
+dependencies = [
+ "cfg-if",
+ "indoc",
+ "libc",
+ "memoffset",
+ "once_cell",
+ "portable-atomic",
+ "pyo3-build-config",
+ "pyo3-ffi",
+ "pyo3-macros",
+ "unindent",
+]
+
+[[package]]
+name = "pyo3-build-config"
+version = "0.22.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b14b5775b5ff446dd1056212d778012cbe8a0fbffd368029fd9e25b514479c38"
+dependencies = [
+ "once_cell",
+ "target-lexicon",
+]
+
+[[package]]
+name = "pyo3-ffi"
+version = "0.22.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9ab5bcf04a2cdcbb50c7d6105de943f543f9ed92af55818fd17b660390fc8636"
+dependencies = [
+ "libc",
+ "pyo3-build-config",
+]
+
+[[package]]
+name = "pyo3-macros"
+version = "0.22.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0fd24d897903a9e6d80b968368a34e1525aeb719d568dba8b3d4bfa5dc67d453"
+dependencies = [
+ "proc-macro2",
+ "pyo3-macros-backend",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "pyo3-macros-backend"
+version = "0.22.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "36c011a03ba1e50152b4b394b479826cad97e7a21eb52df179cd91ac411cbfbe"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "pyo3-build-config",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "quote"
+version = "1.0.47"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "rustversion"
+version = "1.0.23"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cf54715a573b99ac80df0bc206da022bcd442c974952c7b9720069370852e21f"
+
+[[package]]
+name = "syn"
+version = "2.0.119"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "target-lexicon"
+version = "0.12.16"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "61c41af27dd6d1e27b1b16b489db798443478cef1f06a660c96db617ba5de3b1"
+
+[[package]]
+name = "thiserror"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b6aaf5339b578ea85b50e080feb250a3e8ae8cfcdff9a461c9ec2904bc923f52"
+dependencies = [
+ "thiserror-impl",
+]
+
+[[package]]
+name = "thiserror-impl"
+version = "1.0.69"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4fee6c4efc90059e10f81e6d42c60a18f76588c3d74cb83a0b242a2b6c7504c1"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.24"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e6e4313cd5fcd3dad5cafa179702e2b244f760991f45397d14d4ebf38247da75"
+
+[[package]]
+name = "unindent"
+version = "0.2.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7264e107f553ccae879d21fbea1d6724ac785e8c3bfc762137959b5802826ef3"
diff --git a/Cargo.toml b/Cargo.toml
new file mode 100644
index 0000000..723f0bb
--- /dev/null
+++ b/Cargo.toml
@@ -0,0 +1,12 @@
+[workspace]
+resolver = "2"
+members = ["crates/negforge-core", "crates/negforge-py"]
+
+[workspace.package]
+version = "0.1.0"
+edition = "2021"
+license = "MIT"
+
+[profile.release]
+opt-level = 3
+lto = true
diff --git a/README.md b/README.md
index 56b84c5..8b4f451 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,186 @@
-# quantumsim
+# NEGForge
+
+A 1D ballistic-MOSFET electrostatics + NEGF (Non-Equilibrium Green's
+Function) transport simulator. The numerical engine is Rust; the frontend
+is a Jupyter notebook, talking to the engine through a native PyO3
+extension module.
+
+This is a rewrite of a MATLAB nanoelectronics coursework project (see
+[`legacy_matlab/`](legacy_matlab/)) that implemented a self-consistent-style
+1D electrostatic solve, a ballistic (Landauer) current calculation, and an
+NEGF local-density-of-states / charge-density calculation — but never
+actually closed the loop between electrostatics and charge. This rewrite
+does.
+
+## Architecture
+
+```
+crates/negforge-core/   Pure-Rust physics engine (no Python dependency)
+crates/negforge-py/     PyO3 bindings, compiled to the native module
+                         negforge._negforge
+python/negforge/        Pythonic wrapper package (numpy arrays, keyword
+                         device construction, IVCurve helper)
+notebooks/               Jupyter notebook frontend
+legacy_matlab/           Original MATLAB code, kept for reference
+```
+
+`negforge-core` has no PyO3/Python dependency at all — it's a normal Rust
+library with its own unit and integration tests, usable from any Rust
+program. `negforge-py` is a thin binding layer on top of it.
+
+### `negforge-core` modules
+
+- `constants` — physical constants (elementary charge, k_B, h, hbar,
+  electron mass), replacing the undefined `util.const` package the
+  original MATLAB code depended on but never included.
+- `tridiag` — real and complex Thomas-algorithm tridiagonal solvers, used
+  by both the electrostatics and NEGF modules.
+- `device` — device geometry/parameters and the electrostatic
+  (`calc_potential`) and ballistic-current (`calc_current`) calculations.
+- `negf` — the NEGF retarded Green's-function calculation. Rather than the
+  original's dense O(N^3) matrix inversion (`inv()` in MATLAB, fine for a
+  one-off plot but too slow to call repeatedly in a self-consistent loop
+  or a bias sweep), this uses a left-connected recursive Green's-function
+  sweep (O(N)) for the diagonal (local density of states) and a direct
+  tridiagonal solve (O(N)) for each of the two contact-column quantities
+  needed for the charge density — both cross-checked against a dense
+  reference solver in the test suite.
+- `charge` — NEGF-derived electron density, with two bug fixes relative to
+  the original `calc_n()` (see "Deviations from the original" below).
+- `selfconsistent` — **new**: the Poisson↔NEGF self-consistency loop.
+- `sweep` — bias sweeps (`sweep_v_g`, `sweep_v_ds`) and subthreshold-swing
+  extraction, mirroring `plot_Vg_I` / `plot_Vds_I` / `plot_S`.
+
+## Building and running
+
+### Rust engine only
+
+```bash
+cargo build --workspace
+cargo test --workspace
+```
+
+### Python package + notebook
+
+```bash
+python3 -m venv .venv && source .venv/bin/activate
+pip install maturin
+maturin develop --release          # builds negforge-py and installs it editable
+pip install -e ".[notebook]"       # jupyter, matplotlib, ipywidgets
+jupyter notebook notebooks/negforge_demo.ipynb
+```
+
+`maturin develop` compiles `crates/negforge-py` and installs the
+resulting extension inside the `python/negforge/` package (as
+`negforge._negforge`); `python/negforge/__init__.py` wraps it in a
+friendlier, numpy-returning API (`negforge.Device`, `negforge.IVCurve`).
+
+## Physics model
+
+Source-channel-drain chain, 1D, with source/drain treated as long, heavily
+screened contact regions:
+
+- **Electrostatics**: not a full 2D/3D Poisson solve, but a 1D "natural
+  length" (`lambda`) approximation common in compact models for
+  thin-body/double-gate MOSFETs. `lambda = sqrt(k_si/k_ox * d_ch * d_ox /
+  geo)`. The electrostatic operator is tridiagonal with a reflective
+  (Neumann-like) boundary condition at both ends, solved via the Thomas
+  algorithm.
+- **Ballistic current**: the Landauer formula, `I = 2e/h * integral[f_s(E)
+  - f_d(E)] dE`, with Fermi-Dirac occupation of the source/drain
+  reservoirs.
+- **NEGF**: a tight-binding discretization (on-site energy `2t + Psi_f(x)`,
+  hopping `-t`) with open-boundary self-energies at the two contacts,
+  giving the retarded Green's function, local density of states, and
+  (with the self-consistent extension) the charge density.
+
+Units follow the original code throughout: lengths in nm, energies and
+potentials in eV, temperature in Kelvin.
+
+## Deviations from the original MATLAB code
+
+The user requested a faithful port plus closing the missing
+self-consistency loop, not a from-scratch physics redesign. Everything
+below is a deliberate, documented decision, not an accident:
+
+1. **Self-consistent Poisson↔NEGF loop** (`selfconsistent.rs`). The
+   original computed the electrostatic potential once with `rho = 0` and,
+   separately, an NEGF charge density — but never fed one into the other.
+   This is genuinely new: solve electrostatics, compute the NEGF electron
+   density, convert it to a charge density, damp/mix it into `rho`,
+   re-solve, iterate to convergence (or report a `NotConverged` error
+   rather than silently returning garbage).
+
+2. **Charge-density unit fix** (needed for #1 to be numerically meaningful
+   at all). `calc_potential`'s RHS divides `rho + N_dot` by `EPS_0 *
+   k_si`, with `EPS_0` in SI units (F/m); for that division to land in the
+   same eV/nm^2 ballpark as the model's other terms, `rho` must be a
+   genuine volume charge density in C/m^3. The original never exercised
+   this since `rho` was always zero. The self-consistent loop converts the
+   NEGF electron density (computed in nm^-1, the model's native length
+   unit) to C/m^3 by treating the 1D chain as having an implicit unit (1
+   m^2) cross-section and multiplying by `-e`. Skipping this and feeding a
+   raw, differently-scaled density into `rho` was tried first and made the
+   fixed-point iteration diverge by many orders of magnitude — see the
+   module docs in `charge.rs` and `selfconsistent.rs` for the full
+   derivation.
+
+3. **NEGF broadening (`eta`) inside the self-consistent loop.** The
+   original hardcoded `eta = 1e-8` for a one-off local-density-of-states
+   plot. That's fine for a static plot but not for a feedback loop: a
+   bound-state resonance whose energy happens to land within `eta` of an
+   energy-grid point produces a `|G|^2` spike orders of magnitude larger
+   than neighboring grid points, and feeding that grid-alignment-dependent
+   spike back into the electrostatic solve makes the iteration diverge or
+   oscillate. `SelfConsistentOptions::eta` defaults to `0.08` eV (tuned
+   empirically against the default device geometry — see `negf.rs` module
+   docs), large enough to resolve resonances smoothly across iterations.
+   The original's small `eta = 1e-8` is preserved as `negf::DEFAULT_ETA`
+   for the standalone, non-self-consistent `local_density_of_states` /
+   LDOS-plot path, where the original's sharp-peak behavior is what you
+   want to see.
+
+4. **Two `calc_n()` bug fixes** (`charge.rs`), needed for the
+   self-consistent charge density to respond to the actual physics rather
+   than being a constant:
+   - The original multiplies the whole energy sum by a single scalar
+     `f(E_fs)` / `f(E_fd)` (evaluated once, outside the sum) instead of the
+     energy-dependent Fermi occupation `f_s(E)` / `f_d(E)` used everywhere
+     else in the model (e.g. `calc_current`). Fixed to use the proper
+     per-energy weight.
+   - The original reuses one `mask = E > Psi_f(1)` (the *source* band
+     edge) for both the source and drain contact terms. Each contact's
+     contribution is now masked by its own band edge, consistent with how
+     `calc_green` already conditions each contact's self-energy on its own
+     band edge.
+
+5. **O(N) NEGF instead of O(N^3).** See "Architecture" above. Purely a
+   performance change (validated against a dense reference solver in
+   tests); the physics is unchanged.
+
+Everything else — the electrostatic operator, the ballistic current
+formula, the contact self-energy sign convention, the general unit
+handling (nm/eV/K) — is a direct, unmodified port. In particular, the
+model's overall dimensional consistency is inherited as-is from the
+original teaching code (e.g. the electrostatic equation isn't a fully
+rigorous SI-unit Poisson equation); this rewrite does not attempt to
+re-derive the model's physics from first principles, only to make it run,
+close its one clearly-missing feedback loop, and fix the bugs that stood
+in the way of that loop actually doing something.
+
+## Testing
+
+- `crates/negforge-core/src/*.rs` — unit tests per module, including
+  cross-checks of the O(N) tridiagonal/recursive-Green's-function solvers
+  against dense (Gaussian-elimination / full-matrix-inversion) reference
+  implementations on small systems.
+- `crates/negforge-core/tests/self_consistent_realistic_device.rs` — an
+  integration test at the model's default (non-toy) device scale,
+  confirming the self-consistent loop converges and reproduces the
+  expected ballistic-MOSFET trend (current increasing with gate bias).
+- `notebooks/negforge_demo.ipynb` has been executed end-to-end
+  (`jupyter nbconvert --execute`) to confirm the full frontend path works;
+  outputs are cleared before committing since they go stale the moment the
+  engine changes.
+
+Run everything with `cargo test --workspace`.
diff --git a/crates/negforge-core/Cargo.toml b/crates/negforge-core/Cargo.toml
new file mode 100644
index 0000000..625b2fb
--- /dev/null
+++ b/crates/negforge-core/Cargo.toml
@@ -0,0 +1,13 @@
+[package]
+name = "negforge-core"
+version.workspace = true
+edition.workspace = true
+license.workspace = true
+description = "1D ballistic MOSFET electrostatics + NEGF transport engine"
+
+[dependencies]
+num-complex = "0.4"
+thiserror = "1"
+
+[dev-dependencies]
+approx = "0.5"
diff --git a/crates/negforge-core/src/charge.rs b/crates/negforge-core/src/charge.rs
new file mode 100644
index 0000000..22d8529
--- /dev/null
+++ b/crates/negforge-core/src/charge.rs
@@ -0,0 +1,83 @@
+//! NEGF-derived charge density, feeding the self-consistent electrostatic
+//! loop.
+//!
+//! This is a port of `calc_n()` in `legacy_matlab/quantumsim.m`, with two
+//! deliberate bug fixes documented here and in the top-level README:
+//!
+//! 1. The original multiplies the whole energy sum by a single scalar
+//!    `f(E_fs)` / `f(E_fd)` (the Fermi function evaluated once, outside the
+//!    sum) instead of the energy-dependent occupation `f_s(E)` / `f_d(E)`
+//!    used everywhere else in the model (e.g. in `calc_current`). That
+//!    collapses the energy dependence of the injected charge to a single
+//!    number and made the original `calc_n` effectively decorative — it was
+//!    never fed back into `calc_potential`. Here the proper per-energy
+//!    Fermi weight is used, matching `calc_current`.
+//! 2. The original reuses one `mask = E > Psi_f(1)` (the *source* band
+//!    edge) for both the source and drain contact terms. Each contact's
+//!    contribution is now masked by its own band edge, consistent with how
+//!    `calc_green` already conditions each contact's self-energy on its own
+//!    band edge.
+//!
+//! **Third deviation, required for the self-consistent loop to be
+//! numerically meaningful at all:** the original additionally multiplied by
+//! `1e9` to convert the `1/a`-scaled sum from nm^-1 to m^-1 — sensible for
+//! plotting a density on a real-world axis, but `a` and `lambda` (the two
+//! length scales that actually drive `calc_potential`) are expressed in nm
+//! *everywhere else* in this model. Feeding a value that's ~1e9x too large
+//! back into `rho` (which shares scale with `Psi_g/lambda^2` etc., all
+//! still in nm) overwhelms the electrostatic solve and the fixed-point
+//! iteration diverges. [`electron_density`] therefore returns the density
+//! in nm^-1, matching the model's own native length unit. The caller
+//! (`selfconsistent.rs`) is responsible for converting to a proper volume
+//! charge density (m^-1 -> m^-3 assuming unit cross-section, times `-e`)
+//! before assigning it to `rho` — see that module's docs.
+
+use crate::constants::{E, K_B};
+use crate::negf::GreenFunctionResult;
+
+/// Compute the electron number density `n(x)` from a Green's-function
+/// sweep, using the same convention as `legacy_matlab/quantumsim.m`
+/// (see module docs for the two bug fixes applied).
+#[allow(clippy::too_many_arguments)]
+pub fn electron_density(
+    green: &GreenFunctionResult,
+    t_hop: f64,
+    a: f64,
+    psi_f_source: f64,
+    psi_f_drain: f64,
+    e_fs: f64,
+    e_fd: f64,
+    temperature: f64,
+    d_e: f64,
+) -> Vec<f64> {
+    let n_sites = green.g_col_source[0].len();
+    let mut n = vec![0.0; n_sites];
+
+    let f_s = |energy: f64| 1.0 / (((energy - e_fs) * E / (K_B * temperature)).exp() + 1.0);
+    let f_d = |energy: f64| 1.0 / (((energy - e_fd) * E / (K_B * temperature)).exp() + 1.0);
+
+    for (k, &energy) in green.energies.iter().enumerate() {
+        let gamma_s = if energy > psi_f_source {
+            t_hop * green.k_sa[k].re.sin()
+        } else {
+            0.0
+        };
+        let gamma_d = if energy > psi_f_drain {
+            t_hop * green.k_da[k].re.sin()
+        } else {
+            0.0
+        };
+
+        let weight_s = gamma_s * f_s(energy);
+        let weight_d = gamma_d * f_d(energy);
+
+        for (i, n_i) in n.iter_mut().enumerate().take(n_sites) {
+            *n_i += weight_s * green.g_col_source[k][i] + weight_d * green.g_col_drain[k][i];
+        }
+    }
+
+    for v in n.iter_mut() {
+        *v = *v / std::f64::consts::PI * d_e / a;
+    }
+    n
+}
diff --git a/crates/negforge-core/src/constants.rs b/crates/negforge-core/src/constants.rs
new file mode 100644
index 0000000..0fd2982
--- /dev/null
+++ b/crates/negforge-core/src/constants.rs
@@ -0,0 +1,18 @@
+//! Physical constants (SI unless noted). Values are CODATA 2018 recommended
+//! values. The original MATLAB code referenced an undefined `util.const`
+//! package with these same names (`e`, `k_b`, `h`, `h_bar`, `m_e`) — this
+//! module is the concrete replacement for it.
+
+/// Elementary charge, C.
+pub const E: f64 = 1.602_176_634e-19;
+/// Boltzmann constant, J/K.
+pub const K_B: f64 = 1.380_649e-23;
+/// Planck constant, J*s.
+pub const H: f64 = 6.626_070_15e-34;
+/// Reduced Planck constant, J*s.
+pub const H_BAR: f64 = 1.054_571_817e-34;
+/// Electron rest mass, kg.
+pub const M_E: f64 = 9.109_383_701_5e-31;
+/// Vacuum permittivity, F/m. The legacy MATLAB code hardcoded `8.85e-12`
+/// inline instead of using a constants module; we use the precise value.
+pub const EPS_0: f64 = 8.854_187_812_8e-12;
diff --git a/crates/negforge-core/src/device.rs b/crates/negforge-core/src/device.rs
new file mode 100644
index 0000000..e948385
--- /dev/null
+++ b/crates/negforge-core/src/device.rs
@@ -0,0 +1,343 @@
+//! Device geometry, electrostatics and ballistic-current model.
+//!
+//! This is a port of `legacy_matlab/quantumsim.m`. It models a 1D
+//! source-channel-drain chain (source and drain treated as long, heavily
+//! screened contact regions) with:
+//!
+//! - an electrostatic potential solved from a 1D "natural length" (`lambda`)
+//!   approximation rather than a full 2D/3D Poisson equation — the same
+//!   simplification the original coursework code used, common in compact
+//!   models for thin-body/double-gate MOSFETs;
+//! - a ballistic (Landauer) drain current from Fermi-Dirac occupation of the
+//!   source and drain reservoirs.
+//!
+//! Units follow the original code: lengths in nm, energies/potentials in eV,
+//! temperature in Kelvin. `Psi_g`/`Psi_bi`/`Psi_f` store *potential energy*
+//! (already in eV), not electrostatic potential in volts — consistent with
+//! `E_f`/`E_g` also being in eV.
+
+use crate::constants::{E, EPS_0, H_BAR, K_B, M_E};
+use crate::tridiag;
+
+/// Configuration for a [`Device`]. Field defaults match the original
+/// `quantumsim.m` `properties` block.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct DeviceParams {
+    /// Grid spacing, nm.
+    pub a: f64,
+    /// Fermi energy, eV.
+    pub e_f: f64,
+    /// Band gap, eV.
+    pub e_g: f64,
+    /// Drain-source voltage, V.
+    pub v_ds: f64,
+    /// Gate voltage (Psi_g = -e*V_g), V.
+    pub v_g: f64,
+    /// Oxide thickness, nm.
+    pub d_ox: f64,
+    /// Channel (body) thickness, nm.
+    pub d_ch: f64,
+    /// Relative permittivity of the channel (silicon by default).
+    pub k_si: f64,
+    /// Relative permittivity of the oxide.
+    pub k_ox: f64,
+    /// Gate geometry factor (number of gates).
+    pub geo: f64,
+    /// Channel length, nm.
+    pub l_ch: f64,
+    /// If `true` (default, matching the original constructor), the
+    /// source/drain contact length is auto-sized to `floor(lambda) * 15`
+    /// instead of using `l_ds` directly.
+    pub auto_size_contacts: bool,
+    /// Source/drain contact region length, nm. Only used verbatim when
+    /// `auto_size_contacts` is `false`.
+    pub l_ds: f64,
+    /// Fixed dopant charge term added to the mobile charge density in the
+    /// electrostatic solve.
+    pub n_dot: f64,
+    /// Fermi-function tolerance used to bound the ballistic energy window.
+    pub epsilon: f64,
+    /// Source Fermi level, eV.
+    pub e_fs: f64,
+    /// Temperature, K.
+    pub t: f64,
+    /// Energy integration step, eV.
+    pub d_e: f64,
+    /// Effective mass, kg.
+    pub m_eff: f64,
+}
+
+impl Default for DeviceParams {
+    fn default() -> Self {
+        Self {
+            a: 0.5,
+            e_f: 0.15,
+            e_g: 1.0,
+            v_ds: 0.0,
+            v_g: 0.0,
+            d_ox: 5.0,
+            d_ch: 5.0,
+            k_si: 11.2,
+            k_ox: 3.9,
+            geo: 1.0,
+            l_ch: 40.0,
+            auto_size_contacts: true,
+            l_ds: 40.0,
+            n_dot: 0.0,
+            epsilon: 10e-15,
+            e_fs: 0.05,
+            t: 300.0,
+            d_e: 0.001,
+            m_eff: 0.9 * M_E,
+        }
+    }
+}
+
+/// A 1D ballistic-MOSFET device: geometry, electrostatics and derived
+/// transport quantities.
+#[derive(Debug, Clone)]
+pub struct Device {
+    pub params: DeviceParams,
+
+    /// Natural (screening) length, nm.
+    pub lambda: f64,
+    /// Number of grid points along the channel.
+    pub n: usize,
+    /// Last grid index (0-based) belonging to the source region.
+    pub n_left: usize,
+    /// First grid index (0-based) belonging to the drain region.
+    pub n_right: usize,
+
+    /// Gate potential energy profile, eV.
+    pub psi_g: Vec<f64>,
+    /// Built-in potential energy profile, eV.
+    pub psi_bi: Vec<f64>,
+    /// Mobile charge density term entering the electrostatic solve. Zero
+    /// until a self-consistent (NEGF-fed) solve updates it.
+    pub rho: Vec<f64>,
+    /// Solved electrostatic potential energy profile, eV.
+    pub psi_f: Vec<f64>,
+    /// `max(psi_f)`, used as the lower bound of the ballistic energy window.
+    pub psi_0: f64,
+
+    /// Drain Fermi level, eV (`-V_ds + 0.05`, matching the original code).
+    pub e_fd: f64,
+    /// Upper bound of the ballistic energy integration window, eV.
+    pub e_max: f64,
+    /// Nearest-neighbor hopping parameter, eV.
+    pub t_hop: f64,
+}
+
+impl Device {
+    /// Build a new device, replicating the `quantumsim()` MATLAB
+    /// constructor: `lambda` is computed first, then (if
+    /// `auto_size_contacts`) `l_ds` is overwritten from it before the grid
+    /// is sized.
+    pub fn new(mut params: DeviceParams) -> Self {
+        let lambda = Self::compute_lambda(&params);
+        if params.auto_size_contacts {
+            params.l_ds = (lambda.floor()) * 15.0;
+        }
+
+        let (n, n_left, n_right) = Self::compute_grid(&params);
+        let e_fd = -params.v_ds + 0.05;
+        let e_max = params.e_fs - K_B * params.t * params.epsilon.ln() / E;
+        let t_hop = H_BAR.powi(2) / (2.0 * params.m_eff * (params.a * 1e-9).powi(2) * E);
+
+        let mut device = Self {
+            params,
+            lambda,
+            n,
+            n_left,
+            n_right,
+            psi_g: vec![0.0; n],
+            psi_bi: vec![0.0; n],
+            rho: vec![0.0; n],
+            psi_f: vec![0.0; n],
+            psi_0: 0.0,
+            e_fd,
+            e_max,
+            t_hop,
+        };
+        device.init_vectors();
+        device
+    }
+
+    fn compute_lambda(params: &DeviceParams) -> f64 {
+        (params.k_si / params.k_ox * params.d_ch * params.d_ox / params.geo).sqrt()
+    }
+
+    fn compute_grid(params: &DeviceParams) -> (usize, usize, usize) {
+        let l_g = 2.0 * params.l_ds + params.l_ch;
+        let n = (l_g / params.a).floor() as usize + 1;
+        let n_left = (params.l_ds / params.a).floor() as usize;
+        let n_right = n - n_left;
+        (n, n_left, n_right)
+    }
+
+    /// Rebuild `psi_g`, `psi_bi` and reset `rho` from the current
+    /// parameters and region boundaries. Mirrors `init_vectors()`.
+    pub fn init_vectors(&mut self) {
+        self.psi_g = vec![0.0; self.n];
+        self.psi_bi = vec![0.0; self.n];
+        self.rho = vec![0.0; self.n];
+
+        // Channel: indices n_left..n_right-1 (0-based, matching MATLAB's
+        // 1-based `n_left+1 : n_right-1`).
+        for i in self.n_left..self.n_right.saturating_sub(1) {
+            self.psi_g[i] = -self.params.v_g;
+            self.psi_bi[i] = self.params.e_f + self.params.e_g / 2.0;
+        }
+        // Drain: indices n_right-1..N (0-based, matching `n_right:end`).
+        for i in (self.n_right.saturating_sub(1))..self.n {
+            self.psi_g[i] = 0.0;
+            self.psi_bi[i] = -self.params.v_ds;
+        }
+        // Source region (0..n_left) stays at zero, as in the original.
+    }
+
+    /// Build the tridiagonal second-difference operator with the
+    /// reflective (Neumann-like) boundary doubling used for the
+    /// electrostatic solve, scaled by `1/a^2`. Returns `(sub, diag, sup)`.
+    fn laplacian(&self) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
+        let n = self.n;
+        let inv_a2 = 1.0 / self.params.a.powi(2);
+        let diag = vec![-2.0 * inv_a2; n];
+        let mut sup = vec![inv_a2; n - 1];
+        let mut sub = vec![inv_a2; n - 1];
+        sup[0] = 2.0 * inv_a2;
+        sub[n - 2] = 2.0 * inv_a2;
+        (sub, diag, sup)
+    }
+
+    /// Solve for the self-consistent-length electrostatic potential energy
+    /// profile given the current `rho`, storing the result in `psi_f` and
+    /// updating `psi_0 = max(psi_f)`. Mirrors `calc_potential()`.
+    pub fn calc_potential(&mut self) {
+        let (sub, mut diag, sup) = self.laplacian();
+        let inv_lambda2 = 1.0 / self.lambda.powi(2);
+        for d in diag.iter_mut() {
+            *d -= inv_lambda2;
+        }
+
+        let rhs: Vec<f64> = (0..self.n)
+            .map(|i| {
+                (self.rho[i] + self.params.n_dot) / (EPS_0 * self.params.k_si)
+                    - inv_lambda2 * (self.psi_g[i] + self.psi_bi[i])
+            })
+            .collect();
+
+        self.psi_f = tridiag::solve_real(&sub, &diag, &sup, &rhs);
+        self.psi_0 = self.psi_f.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+    }
+
+    /// Ballistic Landauer current for the current potential profile.
+    /// Mirrors `calc_current()`. Units match the original scaling
+    /// convention (see `legacy_matlab/README.md`).
+    pub fn calc_current(&self) -> f64 {
+        let f_s = |e: f64| 1.0 / (((e - self.params.e_fs) * E / (K_B * self.params.t)).exp() + 1.0);
+        let f_d = |e: f64| 1.0 / (((e - self.e_fd) * E / (K_B * self.params.t)).exp() + 1.0);
+
+        let n_steps = ((self.e_max - self.psi_0) / self.params.d_e)
+            .floor()
+            .max(0.0) as usize;
+        let mut sum = 0.0;
+        for k in 0..=n_steps {
+            let energy = self.psi_0 + k as f64 * self.params.d_e;
+            sum += f_s(energy) - f_d(energy);
+        }
+        2.0 * E / crate::constants::H * sum * self.params.d_e * E * 1e-3
+    }
+
+    pub fn set_v_ds(&mut self, v: f64) {
+        self.params.v_ds = v;
+        self.e_fd = -v + 0.05;
+        self.init_vectors();
+    }
+
+    pub fn set_v_g(&mut self, v: f64) {
+        self.params.v_g = v;
+        self.init_vectors();
+    }
+
+    pub fn set_l_ch(&mut self, l: f64) {
+        self.params.l_ch = l;
+        let (n, n_left, n_right) = Self::compute_grid(&self.params);
+        self.n = n;
+        self.n_left = n_left;
+        self.n_right = n_right;
+        self.init_vectors();
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn lambda_matches_closed_form() {
+        let params = DeviceParams::default();
+        let dev = Device::new(params);
+        let expected = (params.k_si / params.k_ox * params.d_ch * params.d_ox / params.geo).sqrt();
+        assert_relative_eq!(dev.lambda, expected, epsilon = 1e-12);
+    }
+
+    #[test]
+    fn zero_bias_zero_gate_gives_symmetric_flat_contacts() {
+        // With V_g = V_ds = 0 and no charge, psi_bi is zero everywhere
+        // outside the channel and the channel barrier is E_f + E_g/2 — the
+        // solved potential should be flat (zero) far into the contacts,
+        // since there is no driving term there.
+        let dev = Device::new(DeviceParams::default());
+        assert_relative_eq!(dev.psi_g[0], 0.0);
+        assert_relative_eq!(dev.psi_bi[0], 0.0);
+        assert_relative_eq!(dev.psi_bi[dev.n - 1], 0.0);
+    }
+
+    #[test]
+    fn calc_potential_is_finite_and_bounded_by_source_terms() {
+        let mut dev = Device::new(DeviceParams::default());
+        dev.calc_potential();
+        assert_eq!(dev.psi_f.len(), dev.n);
+        assert!(dev.psi_f.iter().all(|v| v.is_finite()));
+        // The channel barrier term is E_f + E_g/2; the solved potential
+        // (a screened, smoothed version of the driving terms) should not
+        // wildly overshoot it.
+        let max_drive = dev.params.e_f + dev.params.e_g / 2.0;
+        assert!(dev.psi_0 <= max_drive * 1.01);
+    }
+
+    #[test]
+    fn increasing_v_ds_increases_current() {
+        let mut dev = Device::new(DeviceParams::default());
+        dev.calc_potential();
+        let i0 = dev.calc_current();
+
+        dev.set_v_ds(0.3);
+        dev.calc_potential();
+        let i1 = dev.calc_current();
+
+        assert!(
+            i1 > i0,
+            "current should increase with drain bias: {i0} -> {i1}"
+        );
+    }
+
+    #[test]
+    fn increasing_v_g_increases_current_ballistic_mosfet() {
+        let mut dev = Device::new(DeviceParams::default());
+        dev.set_v_ds(0.3);
+        dev.calc_potential();
+        let i0 = dev.calc_current();
+
+        dev.set_v_g(0.3);
+        dev.calc_potential();
+        let i1 = dev.calc_current();
+
+        assert!(
+            i1 > i0,
+            "current should increase with gate bias: {i0} -> {i1}"
+        );
+    }
+}
diff --git a/crates/negforge-core/src/error.rs b/crates/negforge-core/src/error.rs
new file mode 100644
index 0000000..7ded803
--- /dev/null
+++ b/crates/negforge-core/src/error.rs
@@ -0,0 +1,15 @@
+use thiserror::Error;
+
+#[derive(Debug, Error)]
+pub enum NegForgeError {
+    #[error("self-consistent loop did not converge after {iterations} iterations (residual {residual:.3e}, tolerance {tolerance:.3e})")]
+    NotConverged {
+        iterations: usize,
+        residual: f64,
+        tolerance: f64,
+    },
+    #[error("invalid parameter: {0}")]
+    InvalidParameter(String),
+}
+
+pub type Result<T> = std::result::Result<T, NegForgeError>;
diff --git a/crates/negforge-core/src/lib.rs b/crates/negforge-core/src/lib.rs
new file mode 100644
index 0000000..c97b713
--- /dev/null
+++ b/crates/negforge-core/src/lib.rs
@@ -0,0 +1,21 @@
+//! NEGForge: a Rust port of a 1D ballistic-MOSFET electrostatics + NEGF
+//! transport model (originally MATLAB, see `legacy_matlab/`).
+//!
+//! See the crate modules for the physics of each piece, and the top-level
+//! repository README for the overall architecture, build instructions, and
+//! a list of deliberate deviations from the original code.
+
+pub mod charge;
+pub mod constants;
+pub mod device;
+pub mod error;
+pub mod negf;
+pub mod selfconsistent;
+pub mod sweep;
+mod tridiag;
+
+pub use device::{Device, DeviceParams};
+pub use error::{NegForgeError, Result};
+pub use negf::GreenFunctionResult;
+pub use selfconsistent::{SelfConsistentOptions, SelfConsistentResult};
+pub use sweep::IvPoint;
diff --git a/crates/negforge-core/src/negf.rs b/crates/negforge-core/src/negf.rs
new file mode 100644
index 0000000..59e447b
--- /dev/null
+++ b/crates/negforge-core/src/negf.rs
@@ -0,0 +1,282 @@
+//! Non-Equilibrium Green's Function (NEGF) transport calculation.
+//!
+//! Port of `calc_green` / `calc_n` in `legacy_matlab/quantumsim.m`: builds
+//! the tight-binding Hamiltonian for the potential profile `psi_f`, attaches
+//! open-boundary self-energies at the source/drain contacts, and computes
+//! the retarded Green's function `G^r(E) = [(E + i*eta) I - H]^-1` for a
+//! grid of energies.
+//!
+//! Only the diagonal (local density of states) and the first/last columns
+//! (needed for the injected charge density) are needed, so rather than a
+//! dense O(N^3) matrix inversion (what the original MATLAB `inv()` call
+//! did) this uses:
+//!
+//! - a left-connected recursive Green's function sweep (O(N)) for the full
+//!   diagonal, and
+//! - a direct tridiagonal solve (O(N), via [`tridiag::solve_complex`]) for
+//!   each of the two boundary columns, sidestepping the (easy to get wrong)
+//!   off-diagonal recursive Green's function formulas entirely.
+//!
+//! Both pieces are cross-checked against a dense reference solver in the
+//! test module below.
+
+use num_complex::Complex64;
+
+use crate::device::Device;
+use crate::tridiag;
+
+/// Default imaginary broadening added to the energy, `E + i*eta`, standing
+/// in for an infinitesimal escape rate. Matches `eta = 1i*1e-8` in the
+/// original code; suitable for a one-off local-density-of-states plot at a
+/// fixed potential.
+///
+/// This value is deliberately *not* used inside the self-consistent loop
+/// (see `selfconsistent.rs`): with such a small broadening, a bound-state
+/// resonance whose energy happens to land within `eta` of an energy-grid
+/// point produces a `|G|^2` spike many orders of magnitude larger than
+/// neighboring grid points (the finite chain's true poles are only
+/// infinitesimally broadened). Feeding that spiky, grid-alignment-dependent
+/// charge estimate back into the electrostatic solve makes the fixed-point
+/// iteration diverge. The self-consistent loop instead uses a broadening
+/// comparable to the energy-grid spacing, so resonances are resolved over
+/// several grid points and the charge estimate varies smoothly between
+/// iterations. This is standard NEGF numerical practice (`eta` is already
+/// an artificial regularization in the model, not a physical dephasing
+/// rate) rather than a physics change.
+pub const DEFAULT_ETA: f64 = 1e-8;
+
+/// Result of a [`green_function_sweep`] call: retarded-Green's-function
+/// quantities on an energy grid, one entry per energy point.
+pub struct GreenFunctionResult {
+    /// Energy grid, eV.
+    pub energies: Vec<f64>,
+    /// `Im(G_ii(E)) / a` for every site `i`, indexed `[energy][site]`.
+    /// Proportional to the local density of states (matches `G_r_diag`).
+    pub g_diag: Vec<Vec<f64>>,
+    /// `|G_{i,0}(E)|^2` for every site `i`, indexed `[energy][site]`
+    /// (matches `G_r_1N(:,:,1)`, the source-injected column).
+    pub g_col_source: Vec<Vec<f64>>,
+    /// `|G_{i,N-1}(E)|^2` for every site `i`, indexed `[energy][site]`
+    /// (matches `G_r_1N(:,:,2)`, the drain-injected column).
+    pub g_col_drain: Vec<Vec<f64>>,
+    /// Source contact wavevector at each energy (complex: evanescent below
+    /// the contact band edge). Matches `k_sa`.
+    pub k_sa: Vec<Complex64>,
+    /// Drain contact wavevector at each energy. Matches `k_da`.
+    pub k_da: Vec<Complex64>,
+}
+
+/// Run the NEGF sweep over `energies` for the device's current `psi_f`,
+/// with imaginary broadening `eta` (see [`DEFAULT_ETA`] for guidance on
+/// choosing it). Mirrors `calc_green()`. The caller chooses the energy grid
+/// (the original hardcoded `min(Psi_f) : dE : 0.7*E_max`).
+pub fn green_function_sweep(device: &Device, energies: &[f64], eta: f64) -> GreenFunctionResult {
+    let n = device.n;
+    let t = device.t_hop;
+    let a = device.params.a;
+    let psi_0s = device.psi_f[0];
+    let psi_0d = device.psi_f[n - 1];
+
+    // Bulk on-site energies (before contact self-energy is applied).
+    let bulk_diag: Vec<f64> = device.psi_f.iter().map(|p| 2.0 * t + p).collect();
+    let c1 = bulk_diag[0];
+    let c_n = bulk_diag[n - 1];
+    let off_diag = vec![t; n - 1]; // constant real off-diagonal of (E+i eta)I - H
+
+    let mut g_diag = Vec::with_capacity(energies.len());
+    let mut g_col_source = Vec::with_capacity(energies.len());
+    let mut g_col_drain = Vec::with_capacity(energies.len());
+    let mut k_sa = Vec::with_capacity(energies.len());
+    let mut k_da = Vec::with_capacity(energies.len());
+
+    for &e in energies {
+        let ka_s = Complex64::new(-(2.0 * t - e + psi_0s) / (2.0 * t), 0.0).acos();
+        let ka_d = Complex64::new(-(2.0 * t - e + psi_0d) / (2.0 * t), 0.0).acos();
+        k_sa.push(ka_s);
+        k_da.push(ka_d);
+
+        let mut diag: Vec<Complex64> = bulk_diag
+            .iter()
+            .map(|&d| Complex64::new(e, eta) - Complex64::new(d, 0.0))
+            .collect();
+
+        if e >= psi_0s {
+            let sigma_l = t * (Complex64::i() * ka_s).exp();
+            diag[0] = Complex64::new(e, eta) - (Complex64::new(c1, 0.0) + sigma_l);
+        }
+        if e >= psi_0d {
+            let sigma_r = t * (Complex64::i() * ka_d).exp();
+            diag[n - 1] = Complex64::new(e, eta) - (Complex64::new(c_n, 0.0) + sigma_r);
+        }
+
+        let full_diag = full_diagonal(&off_diag, &diag);
+
+        let mut e_source = vec![Complex64::new(0.0, 0.0); n];
+        e_source[0] = Complex64::new(1.0, 0.0);
+        let col_source = tridiag::solve_complex(&off_diag, &diag, &off_diag, &e_source);
+
+        let mut e_drain = vec![Complex64::new(0.0, 0.0); n];
+        e_drain[n - 1] = Complex64::new(1.0, 0.0);
+        let col_drain = tridiag::solve_complex(&off_diag, &diag, &off_diag, &e_drain);
+
+        g_diag.push(full_diag.iter().map(|g| g.im / a).collect());
+        g_col_source.push(col_source.iter().map(|g| g.norm_sqr()).collect());
+        g_col_drain.push(col_drain.iter().map(|g| g.norm_sqr()).collect());
+    }
+
+    GreenFunctionResult {
+        energies: energies.to_vec(),
+        g_diag,
+        g_col_source,
+        g_col_drain,
+        k_sa,
+        k_da,
+    }
+}
+
+/// Full diagonal of `A^-1` for a tridiagonal `A` with constant real
+/// off-diagonal `t` (i.e. `A[i][i+1] = A[i+1][i] = t` for all `i`) and
+/// complex diagonal `diag`, via the standard left-connected recursive
+/// Green's function sweep.
+fn full_diagonal(off_diag: &[f64], diag: &[Complex64]) -> Vec<Complex64> {
+    let n = diag.len();
+    let mut g_left = vec![Complex64::new(0.0, 0.0); n];
+    g_left[0] = Complex64::new(1.0, 0.0) / diag[0];
+    for i in 1..n {
+        let t = off_diag[i - 1];
+        g_left[i] = Complex64::new(1.0, 0.0) / (diag[i] - g_left[i - 1] * (t * t));
+    }
+
+    let mut full = vec![Complex64::new(0.0, 0.0); n];
+    full[n - 1] = g_left[n - 1];
+    for i in (0..n - 1).rev() {
+        let t = off_diag[i];
+        full[i] = g_left[i] + g_left[i] * g_left[i] * (t * t) * full[i + 1];
+    }
+    full
+}
+
+#[cfg(test)]
+#[allow(clippy::needless_range_loop)]
+mod tests {
+    use super::*;
+    use crate::device::{Device, DeviceParams};
+    use approx::assert_relative_eq;
+
+    /// Dense reference: build the full complex tridiagonal matrix and
+    /// invert it via Gauss-Jordan elimination, for cross-checking the O(N)
+    /// recursive/solve-based implementation above on small systems.
+    fn dense_inverse(off_diag: &[f64], diag: &[Complex64]) -> Vec<Vec<Complex64>> {
+        let n = diag.len();
+        let mut a = vec![vec![Complex64::new(0.0, 0.0); n]; n];
+        for i in 0..n {
+            a[i][i] = diag[i];
+            if i > 0 {
+                a[i][i - 1] = Complex64::new(off_diag[i - 1], 0.0);
+            }
+            if i < n - 1 {
+                a[i][i + 1] = Complex64::new(off_diag[i], 0.0);
+            }
+        }
+        let mut inv = vec![vec![Complex64::new(0.0, 0.0); n]; n];
+        for i in 0..n {
+            inv[i][i] = Complex64::new(1.0, 0.0);
+        }
+        for col in 0..n {
+            let mut pivot = col;
+            for row in (col + 1)..n {
+                if a[row][col].norm() > a[pivot][col].norm() {
+                    pivot = row;
+                }
+            }
+            a.swap(col, pivot);
+            inv.swap(col, pivot);
+            let pivot_val = a[col][col];
+            for k in 0..n {
+                a[col][k] /= pivot_val;
+                inv[col][k] /= pivot_val;
+            }
+            let pivot_row_a = a[col].clone();
+            let pivot_row_inv = inv[col].clone();
+            for row in 0..n {
+                if row == col {
+                    continue;
+                }
+                let factor = a[row][col];
+                for k in 0..n {
+                    a[row][k] -= factor * pivot_row_a[k];
+                    inv[row][k] -= factor * pivot_row_inv[k];
+                }
+            }
+        }
+        inv
+    }
+
+    #[test]
+    fn recursive_diagonal_matches_dense_inverse() {
+        let n = 9;
+        let off_diag: Vec<f64> = vec![0.37; n - 1];
+        let diag: Vec<Complex64> = (0..n)
+            .map(|i| Complex64::new(-1.5 - 0.3 * i as f64, 1e-6))
+            .collect();
+
+        let recursive = full_diagonal(&off_diag, &diag);
+        let dense = dense_inverse(&off_diag, &diag);
+
+        for i in 0..n {
+            assert_relative_eq!(recursive[i].re, dense[i][i].re, epsilon = 1e-8);
+            assert_relative_eq!(recursive[i].im, dense[i][i].im, epsilon = 1e-8);
+        }
+    }
+
+    #[test]
+    fn boundary_columns_match_dense_inverse() {
+        let n = 9;
+        let off_diag: Vec<f64> = vec![0.37; n - 1];
+        let diag: Vec<Complex64> = (0..n)
+            .map(|i| Complex64::new(-1.5 - 0.3 * i as f64, 1e-6))
+            .collect();
+
+        let mut e0 = vec![Complex64::new(0.0, 0.0); n];
+        e0[0] = Complex64::new(1.0, 0.0);
+        let col0 = tridiag::solve_complex(&off_diag, &diag, &off_diag, &e0);
+
+        let mut e_last = vec![Complex64::new(0.0, 0.0); n];
+        e_last[n - 1] = Complex64::new(1.0, 0.0);
+        let col_last = tridiag::solve_complex(&off_diag, &diag, &off_diag, &e_last);
+
+        let dense = dense_inverse(&off_diag, &diag);
+
+        for i in 0..n {
+            assert_relative_eq!(col0[i].re, dense[i][0].re, epsilon = 1e-8);
+            assert_relative_eq!(col0[i].im, dense[i][0].im, epsilon = 1e-8);
+            assert_relative_eq!(col_last[i].re, dense[i][n - 1].re, epsilon = 1e-8);
+            assert_relative_eq!(col_last[i].im, dense[i][n - 1].im, epsilon = 1e-8);
+        }
+    }
+
+    #[test]
+    fn green_function_sweep_runs_and_is_finite() {
+        let mut device = Device::new(DeviceParams {
+            a: 1.0,
+            l_ch: 10.0,
+            l_ds: 10.0,
+            auto_size_contacts: false,
+            ..Default::default()
+        });
+        device.calc_potential();
+
+        let e_min = device.psi_f.iter().cloned().fold(f64::INFINITY, f64::min);
+        let energies: Vec<f64> = (0..20).map(|i| e_min + i as f64 * 0.01).collect();
+        let result = green_function_sweep(&device, &energies, DEFAULT_ETA);
+
+        assert_eq!(result.g_diag.len(), energies.len());
+        for row in &result.g_diag {
+            assert_eq!(row.len(), device.n);
+            assert!(row.iter().all(|v| v.is_finite()));
+        }
+        for row in &result.g_col_source {
+            assert!(row.iter().all(|v| v.is_finite() && *v >= 0.0));
+        }
+    }
+}
diff --git a/crates/negforge-core/src/selfconsistent.rs b/crates/negforge-core/src/selfconsistent.rs
new file mode 100644
index 0000000..970b862
--- /dev/null
+++ b/crates/negforge-core/src/selfconsistent.rs
@@ -0,0 +1,203 @@
+//! Self-consistent Poisson <-> NEGF loop.
+//!
+//! **This module has no equivalent in the original MATLAB code.** The
+//! legacy `quantumsim.m` computed the electrostatic potential once with
+//! `rho = 0` (`calc_potential`) and, separately, an NEGF charge density
+//! (`calc_n`) from that potential — but never fed the computed charge back
+//! into another electrostatic solve. This module closes that loop: solve
+//! electrostatics, compute the NEGF electron density from the result, turn
+//! it into a charge density, damp/mix it into `rho`, and re-solve, until
+//! the potential profile stops changing (or a caller-supplied iteration
+//! budget is exhausted, which is reported as an error rather than silently
+//! returning an unconverged result).
+//!
+//! ## Charge-density units
+//!
+//! `calc_potential`'s RHS divides `rho + N_dot` by `EPS_0 * k_si`, with
+//! `EPS_0` in SI units (F/m). For that division to land back in the same
+//! eV/nm^2 ballpark as the other RHS term (`(Psi_g + Psi_bi) / lambda^2`),
+//! `rho` has to actually be a charge density in C/m^3 — the original code
+//! never exercised this (`rho` was always zero) so it was never validated.
+//! This loop therefore converts the NEGF electron density to a proper
+//! volume charge density before assigning it to `rho`:
+//!
+//! 1. [`charge::electron_density`] returns a linear density in nm^-1
+//!    (natural units for this model's nm-scaled grid).
+//! 2. Convert to m^-1 (`* 1e9`) and treat the 1D chain as having an
+//!    implicit unit (1 m^2) cross-section, so the linear density doubles as
+//!    a volume density.
+//! 3. Multiply by `-e` (electrons are negatively charged) to get `rho` in
+//!    C/m^3: `rho(x) = -e * n_electron(x) * 1e9`.
+//!
+//! `N_dot` (fixed dopant charge) is left as-is (default zero); if it is
+//! used with a nonzero value it should be supplied in the same C/m^3 units
+//! for consistency.
+
+use crate::charge;
+use crate::constants::E as ELEMENTARY_CHARGE;
+use crate::device::Device;
+use crate::error::{NegForgeError, Result};
+use crate::negf::{self, GreenFunctionResult};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct SelfConsistentOptions {
+    /// Maximum number of Poisson<->NEGF iterations before giving up.
+    pub max_iterations: usize,
+    /// Converged when the largest change in `psi_f` between iterations
+    /// (eV) drops below this.
+    pub tolerance: f64,
+    /// Linear mixing factor in `(0, 1]` applied to the charge-density
+    /// update (simple damping to stabilize the fixed-point iteration).
+    pub mixing: f64,
+    /// Upper bound of the NEGF energy grid, as a fraction of `device.e_max`
+    /// (matches the `0.7 * E_max` hardcoded in the original `calc_green`).
+    pub green_energy_fraction: f64,
+    /// Imaginary broadening used for the NEGF evaluations inside the loop.
+    /// Should be several times the energy-grid step (`device.params.d_e`)
+    /// so that resonances are resolved smoothly across iterations instead
+    /// of spiking whenever a grid point happens to land near a pole — see
+    /// the module docs and [`crate::negf::DEFAULT_ETA`]. Increase this if
+    /// the loop fails to converge; decrease for sharper resonance
+    /// resolution (at the cost of a noisier, potentially non-converging
+    /// iteration).
+    pub eta: f64,
+}
+
+impl Default for SelfConsistentOptions {
+    fn default() -> Self {
+        Self {
+            max_iterations: 50,
+            tolerance: 1e-6,
+            mixing: 0.3,
+            green_energy_fraction: 0.7,
+            eta: 0.08,
+        }
+    }
+}
+
+pub struct SelfConsistentResult {
+    pub iterations: usize,
+    pub residual: f64,
+    pub green: GreenFunctionResult,
+    pub electron_density: Vec<f64>,
+}
+
+fn negf_energy_grid(device: &Device, fraction: f64) -> Vec<f64> {
+    let e_min = device.psi_f.iter().cloned().fold(f64::INFINITY, f64::min);
+    let e_max = fraction * device.e_max;
+    let d_e = device.params.d_e;
+    let steps = ((e_max - e_min) / d_e).floor().max(0.0) as usize;
+    (0..=steps).map(|k| e_min + k as f64 * d_e).collect()
+}
+
+/// Run the self-consistent Poisson<->NEGF loop on `device`, mutating its
+/// `rho` and `psi_f` in place. Returns diagnostics plus the final NEGF
+/// sweep and electron density on success.
+pub fn solve_self_consistent(
+    device: &mut Device,
+    opts: &SelfConsistentOptions,
+) -> Result<SelfConsistentResult> {
+    if !(0.0 < opts.mixing && opts.mixing <= 1.0) {
+        return Err(NegForgeError::InvalidParameter(format!(
+            "mixing must be in (0, 1], got {}",
+            opts.mixing
+        )));
+    }
+
+    device.calc_potential();
+
+    let mut last_residual = f64::INFINITY;
+    for iteration in 1..=opts.max_iterations {
+        let energies = negf_energy_grid(device, opts.green_energy_fraction);
+        let green = negf::green_function_sweep(device, &energies, opts.eta);
+        let n_electron = charge::electron_density(
+            &green,
+            device.t_hop,
+            device.params.a,
+            device.psi_f[0],
+            device.psi_f[device.n - 1],
+            device.params.e_fs,
+            device.e_fd,
+            device.params.t,
+            device.params.d_e,
+        );
+
+        let mixing = opts.mixing;
+        for (rho_i, &n_i) in device.rho.iter_mut().zip(n_electron.iter()) {
+            let target_rho = -ELEMENTARY_CHARGE * n_i * 1e9;
+            *rho_i += mixing * (target_rho - *rho_i);
+        }
+
+        let psi_before = device.psi_f.clone();
+        device.calc_potential();
+
+        let residual = psi_before
+            .iter()
+            .zip(device.psi_f.iter())
+            .map(|(a, b)| (a - b).abs())
+            .fold(0.0_f64, f64::max);
+        last_residual = residual;
+
+        if residual < opts.tolerance {
+            return Ok(SelfConsistentResult {
+                iterations: iteration,
+                residual,
+                green,
+                electron_density: n_electron,
+            });
+        }
+    }
+
+    Err(NegForgeError::NotConverged {
+        iterations: opts.max_iterations,
+        residual: last_residual,
+        tolerance: opts.tolerance,
+    })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::device::DeviceParams;
+
+    #[test]
+    fn converges_for_a_small_device() {
+        let mut device = Device::new(DeviceParams {
+            a: 1.0,
+            l_ch: 10.0,
+            l_ds: 10.0,
+            auto_size_contacts: false,
+            d_e: 0.01,
+            ..Default::default()
+        });
+
+        let opts = SelfConsistentOptions {
+            max_iterations: 100,
+            tolerance: 1e-5,
+            mixing: 0.3,
+            green_energy_fraction: 0.7,
+            eta: 0.04,
+        };
+        let result = solve_self_consistent(&mut device, &opts).expect("should converge");
+        assert!(result.residual < opts.tolerance);
+        assert_eq!(result.electron_density.len(), device.n);
+        assert!(result.electron_density.iter().all(|v| v.is_finite()));
+        assert!(device.psi_f.iter().all(|v| v.is_finite()));
+    }
+
+    #[test]
+    fn rejects_invalid_mixing() {
+        let mut device = Device::new(DeviceParams {
+            a: 1.0,
+            l_ch: 10.0,
+            l_ds: 10.0,
+            auto_size_contacts: false,
+            ..Default::default()
+        });
+        let opts = SelfConsistentOptions {
+            mixing: 0.0,
+            ..Default::default()
+        };
+        assert!(solve_self_consistent(&mut device, &opts).is_err());
+    }
+}
diff --git a/crates/negforge-core/src/sweep.rs b/crates/negforge-core/src/sweep.rs
new file mode 100644
index 0000000..4688103
--- /dev/null
+++ b/crates/negforge-core/src/sweep.rs
@@ -0,0 +1,144 @@
+//! Bias sweeps and derived figures of merit (subthreshold swing), mirroring
+//! `plot_Vg_I`, `plot_Vds_I` and `plot_S` in `legacy_matlab/quantumsim.m`.
+
+use crate::device::Device;
+use crate::error::Result;
+use crate::selfconsistent::{self, SelfConsistentOptions};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct IvPoint {
+    pub voltage: f64,
+    pub current: f64,
+}
+
+/// If `Some`, each sweep point solves the full self-consistent Poisson<->NEGF
+/// loop; if `None`, each point uses the decoupled electrostatic solve only
+/// (matching the original scripts' `calc_potential` + `calc_current`).
+pub type SelfConsistency<'a> = Option<&'a SelfConsistentOptions>;
+
+fn inclusive_range(min: f64, max: f64, step: f64) -> Vec<f64> {
+    let steps = ((max - min) / step).round().max(0.0) as usize;
+    (0..=steps).map(|i| min + i as f64 * step).collect()
+}
+
+/// Gate-voltage sweep at fixed drain bias. Mirrors `plot_Vg_I`.
+pub fn sweep_v_g(
+    device: &mut Device,
+    v_min: f64,
+    v_max: f64,
+    step: f64,
+    self_consistency: SelfConsistency,
+) -> Result<Vec<IvPoint>> {
+    let mut points = Vec::new();
+    for v in inclusive_range(v_min, v_max, step) {
+        device.set_v_g(v);
+        match self_consistency {
+            Some(opts) => {
+                selfconsistent::solve_self_consistent(device, opts)?;
+            }
+            None => device.calc_potential(),
+        }
+        points.push(IvPoint {
+            voltage: v,
+            current: device.calc_current(),
+        });
+    }
+    Ok(points)
+}
+
+/// Drain-voltage sweep at fixed gate bias. Mirrors `plot_Vds_I`.
+pub fn sweep_v_ds(
+    device: &mut Device,
+    v_min: f64,
+    v_max: f64,
+    step: f64,
+    self_consistency: SelfConsistency,
+) -> Result<Vec<IvPoint>> {
+    let mut points = Vec::new();
+    for v in inclusive_range(v_min, v_max, step) {
+        device.set_v_ds(v);
+        match self_consistency {
+            Some(opts) => {
+                selfconsistent::solve_self_consistent(device, opts)?;
+            }
+            None => device.calc_potential(),
+        }
+        points.push(IvPoint {
+            voltage: v,
+            current: device.calc_current(),
+        });
+    }
+    Ok(points)
+}
+
+/// Ordinary least-squares fit `y = slope * x + intercept`.
+fn linear_fit(x: &[f64], y: &[f64]) -> (f64, f64) {
+    let n = x.len() as f64;
+    let mean_x = x.iter().sum::<f64>() / n;
+    let mean_y = y.iter().sum::<f64>() / n;
+    let mut num = 0.0;
+    let mut den = 0.0;
+    for (&xi, &yi) in x.iter().zip(y.iter()) {
+        num += (xi - mean_x) * (yi - mean_y);
+        den += (xi - mean_x).powi(2);
+    }
+    let slope = num / den;
+    let intercept = mean_y - slope * mean_x;
+    (slope, intercept)
+}
+
+/// Subthreshold swing (mV/decade equivalent, in the original's V/decade
+/// units) from a gate-voltage sweep, computed as `1 / slope` of a
+/// `log10(I)` vs `V_g` linear fit over `[fit_v_min, fit_v_max]`. Mirrors
+/// the `polyfit` step inside `plot_Vg_I`.
+pub fn subthreshold_swing(points: &[IvPoint], fit_v_min: f64, fit_v_max: f64) -> f64 {
+    let (xs, ys): (Vec<f64>, Vec<f64>) = points
+        .iter()
+        .filter(|p| p.voltage >= fit_v_min && p.voltage <= fit_v_max)
+        .map(|p| (p.voltage, p.current.log10()))
+        .unzip();
+    let (slope, _intercept) = linear_fit(&xs, &ys);
+    1.0 / slope
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::device::DeviceParams;
+
+    #[test]
+    fn v_g_sweep_produces_monotonic_current_for_ballistic_mosfet() {
+        let mut device = Device::new(DeviceParams {
+            v_ds: 0.3,
+            ..Default::default()
+        });
+        let points = sweep_v_g(&mut device, 0.0, 0.4, 0.05, None).unwrap();
+        assert_eq!(points.len(), 9);
+        for w in points.windows(2) {
+            assert!(
+                w[1].current >= w[0].current,
+                "current should be monotonic in V_g"
+            );
+        }
+    }
+
+    #[test]
+    fn subthreshold_swing_is_positive_for_a_reasonable_device() {
+        let mut device = Device::new(DeviceParams {
+            v_ds: 0.3,
+            ..Default::default()
+        });
+        let points = sweep_v_g(&mut device, 0.0, 0.4, 0.02, None).unwrap();
+        let s = subthreshold_swing(&points, 0.0, 0.4);
+        assert!(s > 0.0 && s.is_finite());
+    }
+
+    #[test]
+    fn linear_fit_recovers_known_line() {
+        let x = vec![0.0, 1.0, 2.0, 3.0];
+        let y = vec![1.0, 3.0, 5.0, 7.0];
+        let (slope, intercept) = linear_fit(&x, &y);
+        assert!((slope - 2.0).abs() < 1e-9);
+        assert!((intercept - 1.0).abs() < 1e-9);
+    }
+}
diff --git a/crates/negforge-core/src/tridiag.rs b/crates/negforge-core/src/tridiag.rs
new file mode 100644
index 0000000..5f55caf
--- /dev/null
+++ b/crates/negforge-core/src/tridiag.rs
@@ -0,0 +1,160 @@
+//! Tridiagonal linear system solvers (Thomas algorithm), real and complex.
+//!
+//! Both the electrostatic Poisson-like solve and the NEGF Green's-function
+//! recursion operate on tridiagonal matrices, so this is the one place that
+//! implements the O(N) forward-elimination/back-substitution algorithm.
+
+use num_complex::Complex64;
+
+/// Solve `A x = rhs` for a real tridiagonal `A` given as three diagonals.
+///
+/// `sub[i]` is `A[i+1][i]` (length N-1), `diag[i]` is `A[i][i]` (length N),
+/// `sup[i]` is `A[i][i+1]` (length N-1).
+pub fn solve_real(sub: &[f64], diag: &[f64], sup: &[f64], rhs: &[f64]) -> Vec<f64> {
+    let n = diag.len();
+    assert_eq!(sub.len(), n - 1);
+    assert_eq!(sup.len(), n - 1);
+    assert_eq!(rhs.len(), n);
+
+    let mut c_prime = vec![0.0; n - 1];
+    let mut d_prime = vec![0.0; n];
+
+    c_prime[0] = sup[0] / diag[0];
+    d_prime[0] = rhs[0] / diag[0];
+
+    for i in 1..n {
+        let denom = diag[i]
+            - sub[i - 1]
+                * if i - 1 < c_prime.len() {
+                    c_prime[i - 1]
+                } else {
+                    0.0
+                };
+        if i < n - 1 {
+            c_prime[i] = sup[i] / denom;
+        }
+        d_prime[i] = (rhs[i] - sub[i - 1] * d_prime[i - 1]) / denom;
+    }
+
+    let mut x = vec![0.0; n];
+    x[n - 1] = d_prime[n - 1];
+    for i in (0..n - 1).rev() {
+        x[i] = d_prime[i] - c_prime[i] * x[i + 1];
+    }
+    x
+}
+
+/// Complex counterpart of [`solve_real`], used by the NEGF module where the
+/// matrix `(E + i*eta) I - H` has a real tridiagonal skeleton but complex
+/// diagonal entries (contact self-energies) and the right-hand side may be
+/// complex (unit vectors when extracting a single column of the Green's
+/// function).
+pub fn solve_complex(
+    sub: &[f64],
+    diag: &[Complex64],
+    sup: &[f64],
+    rhs: &[Complex64],
+) -> Vec<Complex64> {
+    let n = diag.len();
+    assert_eq!(sub.len(), n - 1);
+    assert_eq!(sup.len(), n - 1);
+    assert_eq!(rhs.len(), n);
+
+    let mut c_prime = vec![Complex64::new(0.0, 0.0); n - 1];
+    let mut d_prime = vec![Complex64::new(0.0, 0.0); n];
+
+    c_prime[0] = Complex64::new(sup[0], 0.0) / diag[0];
+    d_prime[0] = rhs[0] / diag[0];
+
+    for i in 1..n {
+        let denom = diag[i] - c_prime[i - 1] * sub[i - 1];
+        if i < n - 1 {
+            c_prime[i] = Complex64::new(sup[i], 0.0) / denom;
+        }
+        d_prime[i] = (rhs[i] - d_prime[i - 1] * sub[i - 1]) / denom;
+    }
+
+    let mut x = vec![Complex64::new(0.0, 0.0); n];
+    x[n - 1] = d_prime[n - 1];
+    for i in (0..n - 1).rev() {
+        x[i] = d_prime[i] - c_prime[i] * x[i + 1];
+    }
+    x
+}
+
+#[cfg(test)]
+#[allow(clippy::needless_range_loop)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    /// Dense reference solver (Gaussian elimination with partial pivoting)
+    /// used only in tests to validate the O(N) Thomas-algorithm solvers.
+    fn dense_solve_real(sub: &[f64], diag: &[f64], sup: &[f64], rhs: &[f64]) -> Vec<f64> {
+        let n = diag.len();
+        let mut a = vec![vec![0.0; n]; n];
+        for i in 0..n {
+            a[i][i] = diag[i];
+            if i > 0 {
+                a[i][i - 1] = sub[i - 1];
+            }
+            if i < n - 1 {
+                a[i][i + 1] = sup[i];
+            }
+        }
+        let mut b = rhs.to_vec();
+        for col in 0..n {
+            let mut pivot = col;
+            for row in (col + 1)..n {
+                if a[row][col].abs() > a[pivot][col].abs() {
+                    pivot = row;
+                }
+            }
+            a.swap(col, pivot);
+            b.swap(col, pivot);
+            for row in (col + 1)..n {
+                let factor = a[row][col] / a[col][col];
+                for k in col..n {
+                    a[row][k] -= factor * a[col][k];
+                }
+                b[row] -= factor * b[col];
+            }
+        }
+        let mut x = vec![0.0; n];
+        for i in (0..n).rev() {
+            let mut s = b[i];
+            for j in (i + 1)..n {
+                s -= a[i][j] * x[j];
+            }
+            x[i] = s / a[i][i];
+        }
+        x
+    }
+
+    #[test]
+    fn thomas_matches_dense_gaussian_elimination() {
+        let n = 12;
+        let sub: Vec<f64> = (0..n - 1).map(|i| 1.0 + 0.1 * i as f64).collect();
+        let sup: Vec<f64> = (0..n - 1).map(|i| 0.7 + 0.05 * i as f64).collect();
+        let diag: Vec<f64> = (0..n).map(|i| -3.0 - 0.2 * i as f64).collect();
+        let rhs: Vec<f64> = (0..n).map(|i| (i as f64).sin() + 1.0).collect();
+
+        let fast = solve_real(&sub, &diag, &sup, &rhs);
+        let reference = dense_solve_real(&sub, &diag, &sup, &rhs);
+
+        for (a, b) in fast.iter().zip(reference.iter()) {
+            assert_relative_eq!(a, b, epsilon = 1e-9);
+        }
+    }
+
+    #[test]
+    fn identity_system_returns_rhs() {
+        let n = 5;
+        let sub = vec![0.0; n - 1];
+        let sup = vec![0.0; n - 1];
+        let diag = vec![1.0; n];
+        let rhs = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        let x = solve_real(&sub, &diag, &sup, &rhs);
+        assert_eq!(x, rhs);
+    }
+}
diff --git a/crates/negforge-core/tests/self_consistent_realistic_device.rs b/crates/negforge-core/tests/self_consistent_realistic_device.rs
new file mode 100644
index 0000000..3dd9eb8
--- /dev/null
+++ b/crates/negforge-core/tests/self_consistent_realistic_device.rs
@@ -0,0 +1,37 @@
+//! Integration test at the model's default (non-toy) device scale — the
+//! same geometry a notebook user gets from `Device::new(DeviceParams
+//!::default())`. Confirms the self-consistent loop converges reliably and
+//! produces the expected ballistic-MOSFET trend (current increasing with
+//! gate bias) rather than just checking the tiny unit-test-scale device in
+//! `selfconsistent.rs`.
+
+use negforge_core::selfconsistent::solve_self_consistent;
+use negforge_core::{Device, DeviceParams, SelfConsistentOptions};
+
+#[test]
+fn self_consistent_current_increases_with_gate_bias_at_default_scale() {
+    let opts = SelfConsistentOptions::default();
+
+    let mut low_vg = Device::new(DeviceParams {
+        v_ds: 0.3,
+        v_g: 0.0,
+        ..Default::default()
+    });
+    let r_low = solve_self_consistent(&mut low_vg, &opts).expect("should converge at V_g=0");
+    assert!(r_low.residual < opts.tolerance);
+    let i_low = low_vg.calc_current();
+
+    let mut high_vg = Device::new(DeviceParams {
+        v_ds: 0.3,
+        v_g: 0.3,
+        ..Default::default()
+    });
+    let r_high = solve_self_consistent(&mut high_vg, &opts).expect("should converge at V_g=0.3");
+    assert!(r_high.residual < opts.tolerance);
+    let i_high = high_vg.calc_current();
+
+    assert!(
+        i_high > i_low,
+        "self-consistent current should increase with gate bias: {i_low:e} -> {i_high:e}"
+    );
+}
diff --git a/crates/negforge-py/Cargo.toml b/crates/negforge-py/Cargo.toml
new file mode 100644
index 0000000..4269e45
--- /dev/null
+++ b/crates/negforge-py/Cargo.toml
@@ -0,0 +1,14 @@
+[package]
+name = "negforge-py"
+version.workspace = true
+edition.workspace = true
+license.workspace = true
+description = "PyO3 bindings for negforge-core"
+
+[lib]
+name = "_negforge"
+crate-type = ["cdylib"]
+
+[dependencies]
+negforge-core = { path = "../negforge-core" }
+pyo3 = { version = "0.22", features = ["extension-module"] }
diff --git a/crates/negforge-py/src/lib.rs b/crates/negforge-py/src/lib.rs
new file mode 100644
index 0000000..90e660b
--- /dev/null
+++ b/crates/negforge-py/src/lib.rs
@@ -0,0 +1,210 @@
+// PyO3's #[pymethods] macro expands `?` in ways that trip this lint on
+// PyResult-returning methods; it's a known false positive, not a real
+// no-op conversion in our code.
+#![allow(clippy::useless_conversion)]
+
+use pyo3::exceptions::PyRuntimeError;
+use pyo3::prelude::*;
+
+use negforge_core::{Device, DeviceParams, SelfConsistentOptions};
+
+fn to_py_err(e: negforge_core::NegForgeError) -> PyErr {
+    PyRuntimeError::new_err(e.to_string())
+}
+
+/// Python-facing wrapper around [`negforge_core::Device`].
+#[pyclass(name = "Device")]
+struct PyDevice {
+    inner: Device,
+}
+
+#[pymethods]
+impl PyDevice {
+    #[new]
+    #[pyo3(signature = (
+        a=0.5, e_f=0.15, e_g=1.0, v_ds=0.0, v_g=0.0, d_ox=5.0, d_ch=5.0,
+        k_si=11.2, k_ox=3.9, geo=1.0, l_ch=40.0, auto_size_contacts=true,
+        l_ds=40.0, n_dot=0.0, epsilon=10e-15, e_fs=0.05, t=300.0, d_e=0.001
+    ))]
+    #[allow(clippy::too_many_arguments)]
+    fn new(
+        a: f64,
+        e_f: f64,
+        e_g: f64,
+        v_ds: f64,
+        v_g: f64,
+        d_ox: f64,
+        d_ch: f64,
+        k_si: f64,
+        k_ox: f64,
+        geo: f64,
+        l_ch: f64,
+        auto_size_contacts: bool,
+        l_ds: f64,
+        n_dot: f64,
+        epsilon: f64,
+        e_fs: f64,
+        t: f64,
+        d_e: f64,
+    ) -> Self {
+        let params = DeviceParams {
+            a,
+            e_f,
+            e_g,
+            v_ds,
+            v_g,
+            d_ox,
+            d_ch,
+            k_si,
+            k_ox,
+            geo,
+            l_ch,
+            auto_size_contacts,
+            l_ds,
+            n_dot,
+            epsilon,
+            e_fs,
+            t,
+            d_e,
+            m_eff: 0.9 * negforge_core::constants::M_E,
+        };
+        Self {
+            inner: Device::new(params),
+        }
+    }
+
+    fn calc_potential(&mut self) {
+        self.inner.calc_potential();
+    }
+
+    fn calc_current(&self) -> f64 {
+        self.inner.calc_current()
+    }
+
+    fn set_v_ds(&mut self, v: f64) {
+        self.inner.set_v_ds(v);
+    }
+
+    fn set_v_g(&mut self, v: f64) {
+        self.inner.set_v_g(v);
+    }
+
+    fn set_l_ch(&mut self, l: f64) {
+        self.inner.set_l_ch(l);
+    }
+
+    #[pyo3(signature = (max_iterations=50, tolerance=1e-6, mixing=0.3, eta=0.08))]
+    fn solve_self_consistent(
+        &mut self,
+        max_iterations: usize,
+        tolerance: f64,
+        mixing: f64,
+        eta: f64,
+    ) -> PyResult<(usize, f64)> {
+        let opts = SelfConsistentOptions {
+            max_iterations,
+            tolerance,
+            mixing,
+            green_energy_fraction: 0.7,
+            eta,
+        };
+        let result = negforge_core::selfconsistent::solve_self_consistent(&mut self.inner, &opts)
+            .map_err(to_py_err)?;
+        Ok((result.iterations, result.residual))
+    }
+
+    #[getter]
+    fn psi_f(&self) -> Vec<f64> {
+        self.inner.psi_f.clone()
+    }
+
+    #[getter]
+    fn rho(&self) -> Vec<f64> {
+        self.inner.rho.clone()
+    }
+
+    #[getter]
+    fn n(&self) -> usize {
+        self.inner.n
+    }
+
+    #[getter]
+    fn a(&self) -> f64 {
+        self.inner.params.a
+    }
+
+    /// Natural (screening) length, nm. Named `screening_length` rather than
+    /// `lambda` on the Python side since `lambda` is a reserved keyword.
+    #[getter]
+    fn screening_length(&self) -> f64 {
+        self.inner.lambda
+    }
+
+    fn positions_nm(&self) -> Vec<f64> {
+        (0..self.inner.n)
+            .map(|i| i as f64 * self.inner.params.a)
+            .collect()
+    }
+
+    fn sweep_v_g(
+        &mut self,
+        v_min: f64,
+        v_max: f64,
+        step: f64,
+        self_consistent: bool,
+    ) -> PyResult<(Vec<f64>, Vec<f64>)> {
+        let opts = SelfConsistentOptions::default();
+        let sc = if self_consistent { Some(&opts) } else { None };
+        let points = negforge_core::sweep::sweep_v_g(&mut self.inner, v_min, v_max, step, sc)
+            .map_err(to_py_err)?;
+        Ok((
+            points.iter().map(|p| p.voltage).collect(),
+            points.iter().map(|p| p.current).collect(),
+        ))
+    }
+
+    fn sweep_v_ds(
+        &mut self,
+        v_min: f64,
+        v_max: f64,
+        step: f64,
+        self_consistent: bool,
+    ) -> PyResult<(Vec<f64>, Vec<f64>)> {
+        let opts = SelfConsistentOptions::default();
+        let sc = if self_consistent { Some(&opts) } else { None };
+        let points = negforge_core::sweep::sweep_v_ds(&mut self.inner, v_min, v_max, step, sc)
+            .map_err(to_py_err)?;
+        Ok((
+            points.iter().map(|p| p.voltage).collect(),
+            points.iter().map(|p| p.current).collect(),
+        ))
+    }
+
+    /// 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>>) {
+        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)
+    }
+}
+
+#[pymodule]
+fn _negforge(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_class::<PyDevice>()?;
+    Ok(())
+}
diff --git a/legacy_matlab/README.md b/legacy_matlab/README.md
new file mode 100644
index 0000000..21b124f
--- /dev/null
+++ b/legacy_matlab/README.md
@@ -0,0 +1,22 @@
+# Legacy MATLAB implementation
+
+This directory preserves the original MATLAB coursework code (1D ballistic
+MOSFET / NEGF transport simulator) that the Rust engine in `crates/` and the
+notebook in `notebooks/` are based on. It is kept for reference and is not
+part of the build.
+
+Notes on the original code, carried over into the rewrite's design notes
+(see the top-level README "Physics model & known limitations" section):
+
+- `quantumsim.m` is the most complete version (electrostatics + ballistic
+  current + NEGF local density of states / charge density).
+- `simulation.m`, `simulation_v2.m`, `simulation_v3.m` are earlier
+  procedural-script duplicates of the same model with minor sweep variants.
+- `aufgabe1e.m` is a standalone, simpler Poisson-only exercise.
+- All of them reference an undefined `util.const` package (elementary
+  charge, Boltzmann constant, Planck constant, electron mass) that is not
+  present in this repository, so none of this code runs as-is.
+- The electrostatic solve and the NEGF charge calculation were never wired
+  together into a feedback loop in the original code (`calc_potential`
+  always used `rho = 0`). The Rust rewrite adds that self-consistent loop;
+  see the top-level README for details and caveats.
diff --git a/aufgabe1e.m b/legacy_matlab/aufgabe1e.m
similarity index 100%
rename from aufgabe1e.m
rename to legacy_matlab/aufgabe1e.m
diff --git a/legacy_matlab/quantumsim.m b/legacy_matlab/quantumsim.m
new file mode 100644
index 0000000..af64af5
--- /dev/null
+++ b/legacy_matlab/quantumsim.m
@@ -0,0 +1,337 @@
+classdef quantumsim < handle
+    %UNTITLED2 Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+        %constants, lenth in nm, energy in eV
+        % # of grid points
+        a=0.5;
+        E_f=0.15;       % fermi energy in eV
+        E_g=1;        % band gap in eV
+        V_ds=0.;       % drain-source voltage in V
+        V_g=0;          % gate potential Psi_g=-e*V_g in eV
+        d_ox=5;         % oxide thickness in nm
+        d_ch=5;        % channel thickness in nm
+        e=util.const.e; % elementary charge
+        k_0=8.85e-12;   % dielectric constant
+        k_Si=11.2;     % dielectric constant for Silicon
+        k_ox=3.9;       % dielectric constant oxide
+        geo=1;          % Geometriefaktor # of gates, 'w' wrapgate
+        l_ch=40;        % channel length
+        l_ds=40;        % length of drain and source regions
+        N_dot=0;        % # Dopands
+        lambda
+        L_sparse
+        n_left
+        n_right
+        N
+        Psi_g
+        Psi_bi
+        rho
+        Psi_f
+        Psi_0
+        epsilon=10e-15;         % tolarance for fermi function
+        E_fs=0.05;  %?????
+        E_fd
+        T=300;  %temperature in Kelvin
+        dE=0.001; % energy step
+        DOS
+        
+        E_max
+        E
+        % From here: (Eigenenergies)
+        H
+        G_r_diag
+        G_r_1N
+        t
+        k_sa
+        k_da
+        m=0.9*util.const.m_e;
+    end
+    
+    properties (Dependent)
+        
+    end
+    
+    methods
+        function self =  quantumsim() %constructor
+            self.get_lambda();
+            self.l_ds=floor(self.lambda)*15.;
+            self.get_N();
+            self.get_laplacian();
+            self.E_fd=-self.V_ds+0.05;
+            self.E_max=self.E_fs-util.const.k_b*self.T*log(self.epsilon)/util.const.e;
+            self.E=self.Psi_0:self.dE:self.E_max;
+            self.t = util.const.h_bar^2/(2*self.m*(self.a*10^(-9))^2*util.const.e);
+        end
+        
+        function get_lambda(self)
+            self.lambda=sqrt(self.k_Si/self.k_ox*self.d_ch*self.d_ox/self.geo);
+        end
+        
+        function get_laplacian(self)
+            %tic;
+            super=zeros(self.N,1);
+            super(2:end)=1;
+            super(2)=2;
+            sub=zeros(self.N,1);
+            sub(1:end-1)=1;
+            sub(end-1)=2;
+            middle=zeros(self.N,1);
+            middle(:)=-2;
+            self.L_sparse=spdiags([super,middle,sub],[1,0,-1],self.N,self.N)/self.a.^2;
+            %toc;
+        end
+        
+        function get_N(self)
+            l_g=2*self.l_ds+self.l_ch;
+            self.N=floor(l_g/(self.a))+1;
+            
+            self.n_left=floor(self.l_ds/self.a);      % source region 1 to n-left
+            self.n_right = self.N - self.n_left;       % drain region n_right-end
+            
+        end
+        
+        function init_vectors(self)
+            %create empty arrays
+            self.Psi_g=zeros(self.N,1);
+            self.Psi_bi=zeros(self.N,1);
+            
+            % soure both 0
+            
+            % channel
+            self.Psi_g(self.n_left+1:self.n_right-1) =-self.V_g;
+            self.Psi_bi(self.n_left+1:self.n_right-1)=self.E_f+self.E_g/2.;
+            
+            %drain
+            self.Psi_g(self.n_right:end) =0;
+            self.Psi_bi(self.n_right:end)=-self.V_ds;
+            
+            % create charge density
+            self.rho=zeros(self.N,1);
+        end
+        
+        function calc_potential(self)
+            %tic;
+            self.Psi_f=(self.L_sparse-spdiags(zeros(self.N,1)+1/self.lambda^2,0,self.N,self.N))\...
+                ((self.rho+self.N_dot)/self.k_0/self.k_Si-1/self.lambda^2*(self.Psi_g+self.Psi_bi));
+            %toc;
+            self.Psi_0=max(self.Psi_f);
+            
+        end
+        
+        function I = calc_current(self)
+            self.E=self.Psi_0:self.dE:self.E_max;
+            f_s = @(E) 1./(exp((E-self.E_fs).*util.const.e./util.const.k_b./self.T)+1);
+            f_d = @(E) 1./(exp((E-self.E_fd).*util.const.e./util.const.k_b./self.T)+1);
+            I=2*util.const.e/util.const.h*sum(f_s(self.E)-f_d(self.E))*self.dE*util.const.e*1e-3;
+        end
+        
+        function plot_potential(self)
+            figure, plot((0:self.N-1).*self.a,self.Psi_f);
+        end
+        
+        function set_V_ds(self,V)
+            self.V_ds=V;
+            self.E_fd=-self.V_ds+0.05;
+            self.init_vectors();
+        end
+        
+        function set_V_g(self,V)
+            self.V_g=V;
+            self.init_vectors();
+        end
+        
+        function set_l_ch(self,l)
+            self.l_ch=l;
+            self.get_N();
+            self.get_laplacian();
+            self.init_vectors();
+            
+        end
+        
+        function S = plot_Vg_I(self,min,max,step,plotme)
+            I=[];
+            for V=min:step:max
+                self.set_V_g(V);
+                self.calc_potential();
+                I(end+1)=self.calc_current();
+            end
+            x=min:step:max;
+            y=log10(I);
+            index = (x >= 0) & (x <= 0.4);
+            p = polyfit(x(index),y(index),1);  %# Fit polynomial coefficients for line
+            yfit = p(2)+x.*p(1);                %# Compute the best-fit line
+            
+            if nargin > 4
+            figure, plot(x,y);                      %# Plot the data
+            hold on;              %# Add to the plot
+            plot(x,yfit,'r');     %# Plot the best-fit line           
+            end
+            S=1/p(1);
+        end
+        
+        function plot_Vds_I(self,min,max,step)
+            I=[];
+            for V=min:step:max
+                self.set_V_ds(V);
+                self.calc_potential();
+                I(end+1)=self.calc_current();
+            end
+            x=min:step:max;
+            y=I;                    
+            figure, plot(x,y);                                      
+        end
+        
+        function plot_S(self)
+            S=[];
+            for l=15:5:100
+                set_l_ch(self,l);
+                S(end+1)=self.plot_Vg_I(0,1,0.01);
+            end
+            figure,plot(15:5:100,S);
+        end
+        
+        %task 3
+        function task3(self)
+            
+            self.set_l_ch(40);
+            %check other stuff!
+            Vg_min=0;
+            Vg_max=1;
+            Vg_step=0.1;
+            
+            Vd_min=0;
+            Vd_max=1;
+            Vd_step=0.01;
+            
+            figure();
+            hold on;
+            
+            for Vd=Vg_min:Vg_step:Vg_max
+            self.set_V_g(Vd);
+            I=[];
+            for V=Vd_min:Vd_step:Vd_max
+                self.set_V_ds(V);
+                self.calc_potential();
+                I(end+1)=self.calc_current();
+            end
+            x=Vd_min:Vd_step:Vd_max;
+            y=I;                          
+            plot(x,y);                      %# Plot the data
+                         %# Add to the plot                       
+            end
+            hold off;
+            
+            self.plot_S();
+            
+        end
+        
+
+        
+        function calc_green(self, verbose)
+            self.E=min(self.Psi_f):self.dE:0.7*self.E_max;
+            
+            
+            %self.set_l_ch(100);
+            super=zeros(self.N,1);
+            super(2:end)=1;
+            sub=zeros(self.N,1);
+            sub(1:end-1)=1;
+            middle=zeros(self.N,1);
+            middle(:)=-2;
+
+            
+                  
+            %pot = zeros(self.N,1);
+            pot = self.Psi_f;
+            eta = (1i*1e-8);
+           
+            Phi_0s=self.Psi_f(1);
+            Phi_0d=self.Psi_f(end);
+            
+            faktor_rechts=1;
+            faktor_links=1;
+            
+            self.t = util.const.h_bar^2/(2*self.m*(self.a*10^(-9))^2*util.const.e);
+            self.H=-self.t.*spdiags([super,middle,sub],[1,0,-1],self.N,self.N)+spdiags(pot,0,self.N,self.N);
+            self.G_r_diag = zeros(length(self.E),self.N);
+            self.G_r_1N = zeros(length(self.E),self.N,2);
+            
+            c1=self.H(1);
+            c2=self.H(end);
+            
+            self.k_da=acos(-(2*self.t-self.E+Phi_0d)/2./self.t);
+            self.k_sa=acos(-(2*self.t-self.E+Phi_0s)/2./self.t);
+            
+            for k=1:length(self.E)
+                
+                %linker Kontakt
+                if self.E(k)>=self.Psi_f(1)
+                    %self.k_sa=acos(-(2*self.t-self.E(k)+Phi_0s)/2./self.t);
+                    self.H(1)=c1+faktor_links*self.t*exp(1i*self.k_sa(k)); 
+                else
+                    self.H(1)=c1;
+                end
+                
+                %rechter Kontakt
+                if self.E(k)>=self.Psi_f(end)                   
+                    %self.k_da=acos(-(2*self.t-self.E(k)+Phi_0d)/2./self.t);          
+                    self.H(end)=c2+faktor_rechts*self.t*exp(1i*self.k_da(k));
+                else
+                    self.H(end)=c2;
+                end
+                
+                temp=inv(spdiags(ones(self.N,1)*...
+                    self.E(k)+eta,0,self.N,self.N)-self.H);
+                
+                self.G_r_diag(k,:) = imag(diag(temp/self.a)); 
+                self.G_r_1N(k,:,1) = abs(temp(:,1).^2);
+                self.G_r_1N(k,:,2) = abs(temp(:,end).^2);
+                
+            end
+            
+            %self.DOS=1/pi*self.G_r;
+            
+            if (nargin>1)
+                disp(verbose);
+                figure();
+                hold on
+                imagesc(1:self.N,self.E,1/pi*self.G_r_diag);
+                plot(1:self.N,self.Psi_f,'r','LineWidth',2);
+                xlim(gca,[1 self.N]);
+                ylim(gca,[self.E(1) self.E(end)]);
+                hold off
+                
+                colorbar();
+                set(gca,'Ydir','Normal');
+            end
+          
+            
+            %Ladungstr?gerdichte
+           
+            
+            %figure, plot(n)
+
+        end
+        
+        function calc_n(self)
+            f = @(E) 1./(exp((E).*util.const.e./util.const.k_b./self.T)+1);
+           
+            mask=self.E>self.Psi_f(1);
+            
+            n1 = (1/pi*self.G_r_diag')*f(self.E-self.E_f)'*self.dE/self.a*1e9; %mal dE
+            disp('stop');
+            n=1/pi*(self.t*sin(self.k_sa.*mask)*self.G_r_1N(:,:,1)*f(self.E_fs)+...
+                self.t*sin(self.k_da.*mask)*self.G_r_1N(:,:,2)*f(self.E_fd))*self.dE/self.a*1e9;
+            figure, plot(1:self.N,(n));
+            
+        end
+        
+        %set a, N functionen 
+        
+        
+    end
+    
+end
+
diff --git a/simulation.m b/legacy_matlab/simulation.m
similarity index 52%
rename from simulation.m
rename to legacy_matlab/simulation.m
index 980cd32..422cdf9 100644
--- a/simulation.m
+++ b/legacy_matlab/simulation.m
@@ -1,25 +1,34 @@
 %% constants, lenth in nm, energy in eV
-N=500;          % # of grid points
-
-E_f=0.01;       % fermi energy in eV
-E_g=0.1;        % band gap in eV
-V_ds=0.1;       % drain-source voltage in V
-V_g=0;          % gate potential Psi_g=-e*V_g in eV
-d_ox=1;         % oxide thickness in nm
-d_ch=10;        % channel thickness in nm
+N=171;          % # of grid points
+
+E_f=0.05;       % fermi energy in eV
+E_g=1;        % band gap in eV
+V_ds=0.5;       % drain-source voltage in V
+V_g=0.7;          % gate potential Psi_g=-e*V_g in eV
+d_ox=5;         % oxide thickness in nm
+d_ch=5;        % channel thickness in nm
 e=util.const.e; % elementary charge 
 k_0=8.85e-12;   % dielectric constant
-k_Si=11.68;     % dielectric constant for Silicon
+k_Si=11.2;     % dielectric constant for Silicon
 k_ox=3.9;       % dielectric constant oxide
 geo=1;          % Geometriefaktor # of gates, 'w' wrapgate
-l_ch=10;        % channel length
-l_ds=10;        % length of drain and source regions
+l_ch=30;        % channel length
+l_ds=40;        % length of drain and source regions
 N_dot=0;        % # Dopands
 
 
 lambda=sqrt(k_Si/k_ox*d_ch*d_ox/geo); 
 
-%% create sparse laplacian
+%% calculate source, drain and channel regions
+% calc overall length
+l_g=2*l_ds+l_ch;
+% calculate spacing
+a=l_g/(N-1);
+% calculate indices
+n_left=floor(l_ds/a);   % source region 1 to n-left
+n_right=N-n_left;       % drain region n_right-end
+
+% create sparse laplacian
 tic;
 super=zeros(N,1);
 super(2:end)=1;
@@ -29,19 +38,11 @@
 sub(end-1)=2;
 middle=zeros(N,1);
 middle(:)=-2;
-L_sparse=spdiags([super,middle,sub],[1,0,-1],N,N);
+L_sparse=spdiags([super,middle,sub],[1,0,-1],N,N)/a.^2;
 toc;
 
-%% calculate source, drain and channel regions
-% calc overall length
-l_g=2*l_ds+l_ch;
-% calculate spacing
-a=l_g/(N-1);
-% calculate indices
-n_left=floor(l_ds/a);   % source region 1 to n-left
-n_right=N-n_left;       % drain region n_right-end
 
-%% assign Psi_g and Psi_b for all regions
+% assign Psi_g and Psi_b for all regions
 %create empty arrays
 Psi_g=zeros(N,1);
 Psi_bi=zeros(N,1);
@@ -59,12 +60,28 @@
 % create charge density
 rho=zeros(N,1);
 
-%% solve eq. for Psi_f
+% solve eq. for Psi_f
 % solve sparse (much quicker!)
 tic;
 Psi_f=(L_sparse-spdiags(zeros(N,1)+1/lambda^2,0,N,N))\((rho+N_dot)/k_0/k_Si-1/lambda^2*(Psi_g+Psi_bi));
 toc;
-%% Plot
+% Plot
 figure, plot((0:N-1).*a,Psi_f);
+Psi_0=max(Psi_f);
+
+% 2nd Lecture
+epsilon=10e-15;         % tolarance for fermi function
+E_fs=0.05;  %?????
+E_fd=-V_ds+0.05;
+T=300;  %temperature in Kelvin
+dE=0.001; % energy step
+
+E_max=E_fs-util.const.k_b*T*log(epsilon)/util.const.e;
+E=Psi_0:dE:E_max;
+
+f_s = @(E) 1./(exp((E-E_fs).*util.const.e./util.const.k_b./T)+1);
+f_d = @(E) 1./(exp((E-E_fd).*util.const.e./util.const.k_b./T)+1);
+
+I=2*util.const.e/util.const.h*sum(f_s(E)-f_d(E))*dE*e*1e-3;
 
 
diff --git a/legacy_matlab/simulation_v2.m b/legacy_matlab/simulation_v2.m
new file mode 100644
index 0000000..d865348
--- /dev/null
+++ b/legacy_matlab/simulation_v2.m
@@ -0,0 +1,104 @@
+%% constants, lenth in nm, energy in eV
+N=171;          % # of grid points
+
+E_f=0.05;       % fermi energy in eV
+E_g=1;        % band gap in eV
+V_ds=0.5;       % drain-source voltage in V
+V_g=0;          % gate potential Psi_g=-e*V_g in eV
+d_ox=5;         % oxide thickness in nm
+d_ch=5;        % channel thickness in nm
+e=util.const.e; % elementary charge 
+k_0=8.85e-12;   % dielectric constant
+k_Si=11.2;     % dielectric constant for Silicon
+k_ox=3.9;       % dielectric constant oxide
+geo=1;          % Geometriefaktor # of gates, 'w' wrapgate
+l_ch=30;        % channel length
+l_ds=40;        % length of drain and source regions
+N_dot=0;        % # Dopands
+
+
+lambda=sqrt(k_Si/k_ox*d_ch*d_ox/geo); 
+
+%% calculate source, drain and channel regions
+
+S=[];
+for l_ch=5:50
+% calc overall length
+l_g=2*l_ds+l_ch;
+% calculate spacing
+a=l_g/(N-1);
+% calculate indices
+n_left=floor(l_ds/a);   % source region 1 to n-left
+n_right=N-n_left;       % drain region n_right-end
+
+% create sparse laplacian
+tic;
+super=zeros(N,1);
+super(2:end)=1;
+super(2)=2;
+sub=zeros(N,1);
+sub(1:end-1)=1;
+sub(end-1)=2;
+middle=zeros(N,1);
+middle(:)=-2;
+L_sparse=spdiags([super,middle,sub],[1,0,-1],N,N)/a.^2;
+toc;
+
+
+I=[];
+for V_g=0:0.01:1
+% assign Psi_g and Psi_b for all regions
+%create empty arrays
+Psi_g=zeros(N,1);
+Psi_bi=zeros(N,1);
+
+% soure both 0
+
+% channel
+Psi_g(n_left+1:n_right-1) =-V_g;
+Psi_bi(n_left+1:n_right-1)=E_f+E_g/2.;
+
+%drain
+Psi_g(n_right:end) =0;
+Psi_bi(n_right:end)=-V_ds;
+
+% create charge density
+rho=zeros(N,1);
+
+% solve eq. for Psi_f
+% solve sparse (much quicker!)
+tic;
+Psi_f=(L_sparse-spdiags(zeros(N,1)+1/lambda^2,0,N,N))\((rho+N_dot)/k_0/k_Si-1/lambda^2*(Psi_g+Psi_bi));
+toc;
+% Plot
+%figure, plot((0:N-1).*a,Psi_f);
+Psi_0=max(Psi_f);
+
+% 2nd Lecture
+epsilon=10e-15;         % tolarance for fermi function
+E_fs=0.05;  %?????
+E_fd=-V_ds+0.05;
+T=300;  %temperature in Kelvin
+dE=0.001; % energy step
+
+E_max=E_fs-util.const.k_b*T*log(epsilon)/util.const.e;
+E=Psi_0:dE:E_max;
+
+f_s = @(E) 1./(exp((E-E_fs).*util.const.e./util.const.k_b./T)+1);
+f_d = @(E) 1./(exp((E-E_fd).*util.const.e./util.const.k_b./T)+1);
+
+I(end+1)=2*util.const.e/util.const.h*sum(f_s(E)-f_d(E))*dE*e*1e-3;
+end
+
+x=0:0.01:1;
+y=log10(I);
+index = (x >= 0) & (x <= 0.4);
+p = polyfit(x(index),y(index),1);  %# Fit polynomial coefficients for line
+yfit = p(2)+x.*p(1);  %# Compute the best-fit line
+%figure, plot(x,y);            %# Plot the data
+%hold on;              %# Add to the plot
+%close all
+%plot(x,yfit,'r');     %# Plot the best-fit line
+S(end+1)=1/p(1);
+end
+figure,plot(5:50,S);
diff --git a/legacy_matlab/simulation_v3.m b/legacy_matlab/simulation_v3.m
new file mode 100644
index 0000000..2fed0af
--- /dev/null
+++ b/legacy_matlab/simulation_v3.m
@@ -0,0 +1,101 @@
+%% constants, lenth in nm, energy in eV
+N=170;          % # of grid points
+
+E_f=0.05;       % fermi energy in eV
+E_g=1;        % band gap in eV
+V_ds=0.5;       % drain-source voltage in V
+V_g=0;          % gate potential Psi_g=-e*V_g in eV
+d_ox=1;         % oxide thickness in nm
+d_ch=5;        % channel thickness in nm
+e=util.const.e; % elementary charge 
+k_0=8.85e-12;   % dielectric constant
+k_Si=11.2;     % dielectric constant for Silicon
+k_ox=3.9;       % dielectric constant oxide
+geo=1;          % Geometriefaktor # of gates, 'w' wrapgate
+l_ch=5;        % channel length
+l_ds=40;        % length of drain and source regions
+N_dot=0;        % # Dopands
+
+
+lambda=sqrt(k_Si/k_ox*d_ch*d_ox/geo); 
+
+%% calculate source, drain and channel regions
+
+
+% calc overall length
+l_g=2*l_ds+l_ch;
+% calculate spacing
+a=l_g/(N-1);
+% calculate indices
+n_left=floor(l_ds/a);   % source region 1 to n-left
+n_right=N-n_left;       % drain region n_right-end
+
+% create sparse laplacian
+tic;
+super=zeros(N,1);
+super(2:end)=1;
+super(2)=2;
+sub=zeros(N,1);
+sub(1:end-1)=1;
+sub(end-1)=2;
+middle=zeros(N,1);
+middle(:)=-2;
+L_sparse=spdiags([super,middle,sub],[1,0,-1],N,N)/a.^2;
+toc;
+
+
+I=[];
+for V_g=0:0.01:4
+% assign Psi_g and Psi_b for all regions
+%create empty arrays
+Psi_g=zeros(N,1);
+Psi_bi=zeros(N,1);
+
+% soure both 0
+
+% channel
+Psi_g(n_left+1:n_right-1) =-V_g;
+Psi_bi(n_left+1:n_right-1)=E_f+E_g/2.;
+
+%drain
+Psi_g(n_right:end) =0;
+Psi_bi(n_right:end)=-V_ds;
+
+% create charge density
+rho=zeros(N,1);
+
+% solve eq. for Psi_f
+% solve sparse (much quicker!)
+tic;
+Psi_f=(L_sparse-spdiags(zeros(N,1)+1/lambda^2,0,N,N))\((rho+N_dot)/k_0/k_Si-1/lambda^2*(Psi_g+Psi_bi));
+toc;
+% Plot
+%figure, plot((0:N-1).*a,Psi_f);
+Psi_0=max(Psi_f);
+
+% 2nd Lecture
+epsilon=10e-15;         % tolarance for fermi function
+E_fs=0.05;  %?????
+E_fd=-V_ds+0.05;
+T=300;  %temperature in Kelvin
+dE=0.001; % energy step
+
+E_max=E_fs-util.const.k_b*T*log(epsilon)/util.const.e;
+E=Psi_0:dE:E_max;
+
+f_s = @(E) 1./(exp((E-E_fs).*util.const.e./util.const.k_b./T)+1);
+f_d = @(E) 1./(exp((E-E_fd).*util.const.e./util.const.k_b./T)+1);
+
+I(end+1)=2*util.const.e/util.const.h*sum(f_s(E)-f_d(E))*dE*e*1e-3;
+end
+
+x=0:0.01:4;
+y=log10(I);
+index = (x >= 0) & (x <= 0.4);
+p = polyfit(x(index),y(index),1);  %# Fit polynomial coefficients for line
+yfit = p(2)+x.*p(1);  %# Compute the best-fit line
+figure, plot(x,y);            %# Plot the data
+hold on;              %# Add to the plot
+%close all
+plot(x,yfit,'r');     %# Plot the best-fit line
+S=1/p(1);
diff --git a/legacy_matlab/untitled3.m b/legacy_matlab/untitled3.m
new file mode 100644
index 0000000..239a2da
--- /dev/null
+++ b/legacy_matlab/untitled3.m
@@ -0,0 +1,7 @@
+
+sim=quantumsim();
+sim.init_vectors();
+sim.calc_potential();
+sim.plot_potential();
+sim.calc_green();
+sim.calc_n();
\ No newline at end of file
diff --git a/notebooks/negforge_demo.ipynb b/notebooks/negforge_demo.ipynb
new file mode 100644
index 0000000..df4122a
--- /dev/null
+++ b/notebooks/negforge_demo.ipynb
@@ -0,0 +1,257 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "fa7c50df",
+   "metadata": {},
+   "source": [
+    "# negforge: 1D ballistic MOSFET simulator\n",
+    "\n",
+    "This notebook is the Python frontend for the `negforge` Rust engine (`crates/negforge-core`, exposed here via PyO3 as the `negforge` package). It reproduces every plot from the original MATLAB coursework code in `legacy_matlab/` (`plot_potential`, `plot_Vg_I`, `plot_Vds_I`, the NEGF local-density-of-states heatmap) plus a self-consistent Poisson&harr;NEGF solve that the original code never had.\n",
+    "\n",
+    "**Before running:** build and install the extension once from the repository root:\n",
+    "```bash\n",
+    "pip install maturin\n",
+    "maturin develop --release\n",
+    "```\n",
+    "See the top-level `README.md` for the physics background, unit conventions, and a list of deliberate deviations from the original MATLAB model (most importantly: the charge-feedback loop and two bug fixes in the NEGF charge-density formula)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2e8ff408",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import negforge\n",
+    "\n",
+    "plt.rcParams[\"figure.figsize\"] = (7, 4.5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d849bc3b",
+   "metadata": {},
+   "source": [
+    "## 1. Electrostatic potential profile\n",
+    "\n",
+    "A device with a gate and drain bias applied. `calc_potential()` is the **decoupled** solve (matches the original `calc_potential()`: charge density `rho` fixed at zero)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f2989dc5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dev = negforge.Device(v_ds=0.3, v_g=0.2)\n",
+    "dev.calc_potential()\n",
+    "\n",
+    "plt.plot(dev.positions_nm, dev.psi_f)\n",
+    "plt.xlabel(\"position (nm)\")\n",
+    "plt.ylabel(r\"$\\Psi_f$ (eV)\")\n",
+    "plt.title(\"Decoupled electrostatic potential (V$_g$=0.2V, V$_{ds}$=0.3V)\")\n",
+    "plt.grid(alpha=0.3)\n",
+    "plt.show()\n",
+    "\n",
+    "print(dev)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ab9f3af6",
+   "metadata": {},
+   "source": [
+    "## 2. Self-consistent Poisson&harr;NEGF solve\n",
+    "\n",
+    "This closes the feedback loop the original MATLAB code never implemented: the NEGF-computed electron charge is fed back into the electrostatic solve and iterated to convergence. Compare against the decoupled potential above &mdash; the self-consistent charge screens (flattens) the barrier somewhat."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9e0ec08e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dev_decoupled = negforge.Device(v_ds=0.3, v_g=0.2).calc_potential()\n",
+    "psi_decoupled = dev_decoupled.psi_f.copy()\n",
+    "\n",
+    "dev_sc = negforge.Device(v_ds=0.3, v_g=0.2)\n",
+    "iterations, residual = dev_sc.solve_self_consistent()\n",
+    "print(f\"converged in {iterations} iterations, residual={residual:.3e} eV\")\n",
+    "\n",
+    "plt.plot(dev_decoupled.positions_nm, psi_decoupled, label=\"decoupled (rho=0)\")\n",
+    "plt.plot(dev_sc.positions_nm, dev_sc.psi_f, label=\"self-consistent\")\n",
+    "plt.xlabel(\"position (nm)\")\n",
+    "plt.ylabel(r\"$\\Psi_f$ (eV)\")\n",
+    "plt.title(\"Decoupled vs. self-consistent potential\")\n",
+    "plt.legend()\n",
+    "plt.grid(alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4d529976",
+   "metadata": {},
+   "source": [
+    "## 3. Transfer characteristic (I&ndash;V$_g$) and subthreshold swing\n",
+    "\n",
+    "Mirrors `plot_Vg_I` / the `polyfit`-based swing extraction in the original code. The self-consistent sweep is slower (a full Poisson&ndash;NEGF solve per point) so it uses fewer points than the decoupled one."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0c36f653",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dev = negforge.Device(v_ds=0.3)\n",
+    "curve_decoupled = dev.sweep_v_g(0.0, 0.6, 0.01, self_consistent=False)\n",
+    "\n",
+    "dev_sc = negforge.Device(v_ds=0.3)\n",
+    "curve_sc = dev_sc.sweep_v_g(0.0, 0.4, 0.05, self_consistent=True)\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "ax.semilogy(curve_decoupled.voltage, curve_decoupled.current, label=\"decoupled\")\n",
+    "ax.semilogy(curve_sc.voltage, curve_sc.current, \"o-\", label=\"self-consistent\")\n",
+    "ax.set_xlabel(r\"$V_g$ (V)\")\n",
+    "ax.set_ylabel(\"current (a.u., see README for units)\")\n",
+    "ax.set_title(\"Transfer characteristic\")\n",
+    "ax.legend()\n",
+    "ax.grid(alpha=0.3, which=\"both\")\n",
+    "plt.show()\n",
+    "\n",
+    "s_decoupled = curve_decoupled.subthreshold_swing(0.0, 0.4)\n",
+    "s_sc = curve_sc.subthreshold_swing(0.0, 0.4)\n",
+    "print(f\"subthreshold swing (decoupled):        {s_decoupled * 1000:.1f} mV/decade\")\n",
+    "print(f\"subthreshold swing (self-consistent):  {s_sc * 1000:.1f} mV/decade\")\n",
+    "print(\"(ideal thermal limit at 300 K is ln(10)*kT/e = 59.6 mV/decade)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "572cc359",
+   "metadata": {},
+   "source": [
+    "## 4. Output characteristic (I&ndash;V$_{ds}$)\n",
+    "\n",
+    "Mirrors `plot_Vds_I`, swept at a few gate voltages."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "40921152",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots()\n",
+    "for v_g in [0.0, 0.1, 0.2, 0.3]:\n",
+    "    dev = negforge.Device(v_g=v_g)\n",
+    "    curve = dev.sweep_v_ds(0.0, 0.6, 0.02, self_consistent=False)\n",
+    "    ax.plot(curve.voltage, curve.current, label=f\"$V_g$={v_g} V\")\n",
+    "ax.set_xlabel(r\"$V_{ds}$ (V)\")\n",
+    "ax.set_ylabel(\"current (a.u., see README for units)\")\n",
+    "ax.set_title(\"Output characteristic\")\n",
+    "ax.legend()\n",
+    "ax.grid(alpha=0.3)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a0cc69e1",
+   "metadata": {},
+   "source": [
+    "## 5. Local density of states (NEGF)\n",
+    "\n",
+    "Mirrors `calc_green()`'s `imagesc` plot: the local density of states as a function of position and energy, with the electrostatic potential profile overlaid."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47ab13a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dev = negforge.Device(v_ds=0.3, v_g=0.2)\n",
+    "dev.solve_self_consistent()\n",
+    "energies, ldos = dev.local_density_of_states()\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(8, 5))\n",
+    "mesh = ax.pcolormesh(\n",
+    "    dev.positions_nm, energies, ldos, shading=\"auto\", cmap=\"inferno\", vmax=np.percentile(ldos, 99)\n",
+    ")\n",
+    "ax.plot(dev.positions_nm, dev.psi_f, color=\"cyan\", linewidth=1.5, label=r\"$\\Psi_f$\")\n",
+    "ax.set_xlabel(\"position (nm)\")\n",
+    "ax.set_ylabel(\"energy (eV)\")\n",
+    "ax.set_title(\"Local density of states\")\n",
+    "ax.legend(loc=\"upper right\")\n",
+    "fig.colorbar(mesh, ax=ax, label=\"LDOS (a.u.)\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8ba69efc",
+   "metadata": {},
+   "source": [
+    "## 6. Interactive exploration\n",
+    "\n",
+    "Drag the sliders to see how gate and drain bias reshape the (decoupled) potential profile in real time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9627ce67",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ipywidgets import interact, FloatSlider\n",
+    "\n",
+    "\n",
+    "def plot_potential(v_g=0.2, v_ds=0.3):\n",
+    "    dev = negforge.Device(v_g=v_g, v_ds=v_ds)\n",
+    "    dev.calc_potential()\n",
+    "    plt.figure(figsize=(7, 4.5))\n",
+    "    plt.plot(dev.positions_nm, dev.psi_f)\n",
+    "    plt.xlabel(\"position (nm)\")\n",
+    "    plt.ylabel(r\"$\\Psi_f$ (eV)\")\n",
+    "    plt.ylim(-1.0, 1.0)\n",
+    "    plt.grid(alpha=0.3)\n",
+    "    plt.title(f\"V$_g$={v_g:.2f} V, V$_{{ds}}$={v_ds:.2f} V\")\n",
+    "    plt.show()\n",
+    "\n",
+    "\n",
+    "interact(\n",
+    "    plot_potential,\n",
+    "    v_g=FloatSlider(min=-0.5, max=1.0, step=0.05, value=0.2),\n",
+    "    v_ds=FloatSlider(min=0.0, max=1.0, step=0.05, value=0.3),\n",
+    ");"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..5dd36da
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,21 @@
+[build-system]
+requires = ["maturin>=1.5,<2.0"]
+build-backend = "maturin"
+
+[project]
+name = "negforge"
+version = "0.1.0"
+description = "1D ballistic-MOSFET electrostatics + NEGF transport simulator (Rust engine, Python frontend)"
+requires-python = ">=3.9"
+readme = "README.md"
+license = { text = "MIT" }
+dependencies = ["numpy>=1.22"]
+
+[project.optional-dependencies]
+notebook = ["jupyter", "matplotlib", "ipywidgets"]
+
+[tool.maturin]
+manifest-path = "crates/negforge-py/Cargo.toml"
+python-source = "python"
+module-name = "negforge._negforge"
+features = ["pyo3/extension-module"]
diff --git a/python/negforge/__init__.py b/python/negforge/__init__.py
new file mode 100644
index 0000000..41cb5fb
--- /dev/null
+++ b/python/negforge/__init__.py
@@ -0,0 +1,155 @@
+"""NEGForge: 1D ballistic-MOSFET electrostatics + NEGF transport simulator.
+
+The numerical engine (self-consistent Poisson/NEGF solve, tridiagonal linear
+algebra, recursive Green's-function evaluation) lives in Rust
+(``negforge-core``, exposed here via PyO3 as ``negforge._negforge``).
+This package is a thin, notebook-friendly Python wrapper around it: numpy
+arrays instead of raw tuples, keyword-argument device construction, and a
+small `IVCurve` helper for bias sweeps. See the top-level repository README
+for the physics background, unit conventions, and known limitations.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ._negforge import Device as _RustDevice
+
+__all__ = ["Device", "IVCurve"]
+
+
+class IVCurve:
+    """Voltage/current arrays from a bias sweep ([`Device.sweep_v_g`] /
+    [`Device.sweep_v_ds`])."""
+
+    def __init__(self, voltage, current):
+        self.voltage = np.asarray(voltage)
+        self.current = np.asarray(current)
+
+    def __iter__(self):
+        return iter((self.voltage, self.current))
+
+    def __repr__(self):
+        return f"IVCurve({len(self.voltage)} points)"
+
+    def subthreshold_swing(self, v_min: float = 0.0, v_max: float = 0.4) -> float:
+        """Subthreshold swing from a `log10(I)` vs voltage linear fit over
+        `[v_min, v_max]`, matching `plot_Vg_I`'s `polyfit` step in the
+        original MATLAB code."""
+        mask = (self.voltage >= v_min) & (self.voltage <= v_max)
+        if mask.sum() < 2:
+            raise ValueError("not enough points in [v_min, v_max] to fit a swing")
+        log_i = np.log10(self.current[mask])
+        slope, _ = np.polyfit(self.voltage[mask], log_i, 1)
+        return 1.0 / slope
+
+
+class Device:
+    """A 1D ballistic MOSFET device.
+
+    Thin, Pythonic wrapper around the Rust engine: numpy arrays for
+    profiles, and dataclass-like results instead of raw tuples. Keyword
+    arguments match the Rust `DeviceParams` field names (see the top-level
+    README for units and defaults), e.g.::
+
+        dev = negforge.Device(v_ds=0.3, v_g=0.2, l_ch=40.0)
+        dev.solve_self_consistent()
+        current = dev.calc_current()
+    """
+
+    def __init__(self, **kwargs):
+        self._inner = _RustDevice(**kwargs)
+        self._inner.calc_potential()
+
+    def __repr__(self):
+        return f"Device(n={self.n}, a={self.a} nm, screening_length={self.screening_length:.3f} nm)"
+
+    # -- electrostatics ----------------------------------------------------
+
+    def calc_potential(self) -> "Device":
+        """Decoupled electrostatic solve (`rho` unchanged). Matches the
+        original `calc_potential()`."""
+        self._inner.calc_potential()
+        return self
+
+    def solve_self_consistent(
+        self,
+        max_iterations: int = 50,
+        tolerance: float = 1e-6,
+        mixing: float = 0.3,
+        eta: float = 0.08,
+    ) -> tuple[int, float]:
+        """Run the self-consistent Poisson<->NEGF loop (not present in the
+        original MATLAB code — see the top-level README). Returns
+        `(iterations, residual)` on success; raises `RuntimeError` if it
+        does not converge within `max_iterations`.
+        """
+        return self._inner.solve_self_consistent(max_iterations, tolerance, mixing, eta)
+
+    # -- bias control --------------------------------------------------------
+
+    def set_v_ds(self, v: float) -> "Device":
+        self._inner.set_v_ds(v)
+        return self
+
+    def set_v_g(self, v: float) -> "Device":
+        self._inner.set_v_g(v)
+        return self
+
+    def set_l_ch(self, l: float) -> "Device":
+        self._inner.set_l_ch(l)
+        return self
+
+    # -- observables -----------------------------------------------------------
+
+    def calc_current(self) -> float:
+        return self._inner.calc_current()
+
+    @property
+    def psi_f(self) -> np.ndarray:
+        """Solved electrostatic potential energy profile, eV."""
+        return np.asarray(self._inner.psi_f)
+
+    @property
+    def rho(self) -> np.ndarray:
+        """Charge density term entering the electrostatic solve."""
+        return np.asarray(self._inner.rho)
+
+    @property
+    def positions_nm(self) -> np.ndarray:
+        return np.asarray(self._inner.positions_nm())
+
+    @property
+    def n(self) -> int:
+        return self._inner.n
+
+    @property
+    def a(self) -> float:
+        return self._inner.a
+
+    @property
+    def screening_length(self) -> float:
+        """Natural (screening) length lambda, nm."""
+        return self._inner.screening_length
+
+    def local_density_of_states(self):
+        """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.
+        """
+        energies, ldos = self._inner.local_density_of_states()
+        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:
+        """Gate-voltage sweep at the device's current drain bias. Matches
+        `plot_Vg_I`."""
+        voltage, current = self._inner.sweep_v_g(v_min, v_max, step, self_consistent)
+        return IVCurve(voltage, current)
+
+    def sweep_v_ds(self, v_min: float, v_max: float, step: float, self_consistent: bool = False) -> IVCurve:
+        """Drain-voltage sweep at the device's current gate bias. Matches
+        `plot_Vds_I`."""
+        voltage, current = self._inner.sweep_v_ds(v_min, v_max, step, self_consistent)
+        return IVCurve(voltage, current)