Skip to content

Latest commit

 

History

History

perovskite-sim

README.md

perovskite-sim

The Python simulation package, FastAPI backend, and Vite/TypeScript frontend that make up the SolarLab simulator.

Start here: the root README covers installation, physics, equations, UI walkthrough, and shipped presets. This file is a package-level orientation to the perovskite-sim/ subtree. The 2026-08-11 technical manual is distributed separately. Current solver capability contracts and the evidence registry are linked below.


Layout

perovskite-sim/
├── perovskite_sim/   Python library: models, physics, solvers, experiments, 2D
├── backend/          FastAPI HTTP wrapper and SSE job dispatch
├── frontend/         Vite + TypeScript + Plotly workstation
├── configs/          Three shipped research YAMLs; two frontend catalog entries
├── reproducibility/  Frozen baselines, schemas, hashes, benchmarks, P1 gaps
├── docs/             Package-specific model, benchmark, and implementation docs
├── scripts/          CLI, validation, plotting, import, and probe tools
├── tests/            Unit, integration, regression, and validation suites
├── notebooks/        Exploratory benchmarks
├── Dockerfile.backend  Backend container definition
├── pyproject.toml    Package metadata and pytest configuration
├── .gitignore        Package-local generated-output rules
├── CLAUDE.md         Package development and validation guidance
└── README.md         This file

Local runs may create an ignored outputs/ directory. A result is not part of the remote repository or reproducibility registry merely because it exists there.


Quick Install

From this directory, use an activated Python 3.11+ environment and Node.js 22.12+ for the Vite 8 frontend. The core package and HTTP-service dependencies are separate:

python -m pip install -e ".[dev]"
python -m pip install -r backend/requirements.txt
cd frontend
npm ci

Tests

python -m pytest -q tests/reproducibility/test_research_presets.py tests/unit/backend/test_scaps_inline_config.py tests/unit/experiments/test_plot_calado_fig1f.py

These are the research-preset checks; python -m pytest -q runs the default non-slow suite. The current registry contains 55 configurations: 3 shipped and 52 test fixtures. See Research Presets.

Transient solves retain the historical scalar absolute tolerance by default. For the opt-in reference-scaled policy and the required three-level tolerance study, see ComponentwiseTolerancePolicy.md. For compact-support Poole-Frenkel, interface-density, and thermionic-cap sensitivity ladders, see RhsRegularizationPolicy.md. The staged implementation, test matrices, promotion gates, and explicitly parked physics are tracked in the 2026-08-14 physics/numerics hardening roadmap (archived plan, see CLAUDE.md → Archived documents). The corresponding Phase 1 implementation status, immutable certificate IDs, failed/partial lanes, and evidence boundaries are recorded in the Phase 1 implementation and evidence report. The first post-P1 impedance prerequisite, a residual- and conservation-certified mobile-ion DC state, is documented in 1D ion-aware DC certification closure. The opt-in Phase 2 mass/storage and central-finite-difference linearization slice is documented in ion-aware impedance reference engine. The subsequent exact-discrete-Poisson, analytic SG transport, analytic bulk recombination, defect-free local and clamp-inactive cross-node/projected interface SRH, positive-density shared-occupancy and additive two-sided SRH, the residual-resolved interior QSS interface root, selective-contact, analytic field-mobility, and adaptive per-column operator check is documented in ion-aware structured Jacobian comparison.

Evidence labels matter: a passing load_only or internal regression is not an external validation. See reproducibility/README.md for the authoritative status and limitations of every shipped config.


Notebooks

Interactive notebooks under notebooks/:

Notebook Topic
JvHysteresis.ipynb J-V sweep with hysteresis
Impedance.ipynb Impedance spectroscopy (Nyquist plot)
Degradation.ipynb Long-term degradation simulation

Benchmark scripts (.py, run with python):

Script Topic
IonmongerBenchmark.py Exploratory IonMonger paper-informed diagnostic
ComprehensiveBenchmark.py Exploratory multi-physics diagnostic
E2eNotebookVsApi.py Notebook vs API parity check

Physical Model

The simulator models a thin-film solar cell as a stack of semiconductor layers between two metallic contacts. The default solver is 1D along the device stack; the 2D extension extrudes the same stack laterally for microstructure and grain-boundary studies. Light enters from one side, generates electron-hole pairs, and the built-in electric field separates them to produce current.

Electrical coordinate, layer order, and contact-potential sources

Band bending, quasi-Fermi levels, and abrupt-interface closures

Numerical drivers, variable sets, and certification paths

No single driver contains every optional model. transient, steady_state, quasi_fermi, quasi_fermi_frequency, and 2D use different unknowns and certification checks. The experiment selects the driver explicitly; unsupported combinations fail before the numerical solve. Likewise, legacy, fast, and full are feature ceilings rather than accuracy grades.

Current Research Presets

Config Structure Ions Optics
scaps_mirror_v2 Glass / HTL / PVK / ETL, SCAPS partner device No TMM
calado2016_fig1f Calado 2016 Fig. 1f toy device Yes Beer-Lambert
calado2016_ion_sweep Internal dense ion-sweep figure Yes Uniform generation under its declared waveform

Use Fast for the SCAPS reference, Legacy for the original Calado comparison, and Full for the internal ion-sweep preset. All three are API-accessible; the frontend catalog shows the SCAPS and internal ion-sweep entries. Historical configurations discussed below remain under tests/fixtures/configs/ and are not additional shipped presets.

Continuous chi/Eg grading changes electrical transport only by default, so the historical CIGS notch studies must not be interpreted as graded-optics Jsc/PCE optimization. The separate research preset cigs_graded_optics explicitly couples the same composition coordinate to Minoura CIGS n,k TMM slices; it remains an internal constitutive/numerical model rather than a calibrated device.


Initial Conditions and Boundary Conditions

The simulator solves the coupled 1D drift-diffusion + Poisson + mobile-ion system. Below are the conditions applied at the device contacts (boundary conditions) and the strategies used to seed the state vector (initial conditions).

Governing Equations

Equation PDE
Poisson $\partial/\partial x(\varepsilon_0 \varepsilon_r , \partial\varphi/\partial x) = -\rho$
Electron continuity $\partial n/\partial t = (1/q) , \partial J_n/\partial x + G - R$
Hole continuity $\partial p/\partial t = -(1/q) , \partial J_p/\partial x + G - R$
Ion (vacancy) continuity $\partial P/\partial t = -\partial F_P/\partial x$

State vector per grid node: $\mathbf{y} = (n, p, P)$ — electron density, hole density, and positive-ion (vacancy) density. Dual-species mode adds a negative-ion field $P^-$.


Boundary Conditions

Electrostatic Potential (Poisson) — Dirichlet

Contact Value
Left ($x = 0$) $\varphi = 0$ (grounded)
Right ($x = L$) $\varphi = V_{\text{bi}} - V_{\text{app}}$

Forward bias ($V_{\text{app}} > 0$) reduces the built-in field; $V_{\text{app}} \approx V_{\text{oc}}$ yields near-open-circuit conditions.

The Poisson operator uses harmonic-mean face permittivities:

$$\tilde{\varepsilon}_{i+\frac{1}{2}} = \frac{2,\varepsilon_r[i],\varepsilon_r[i+1]}{\varepsilon_r[i] + \varepsilon_r[i+1]}$$

This is the exact series-capacitor result for a sharp dielectric interface and avoids the field concentration artefact of nodal averaging.

Source: perovskite_sim/physics/poisson.py


Electron and Hole Densities — Ohmic or Robin Contacts

By default, both contacts are treated as ideal ohmic contacts. Carrier densities at the boundaries are clamped to the thermal-equilibrium values derived from the doping of the outermost layers:

$$n \cdot p = n_i^2 \quad\text{(mass-action law)}$$

$$n - p = N_D - N_A \quad\text{(charge neutrality)}$$

Solved via the numerically stable two-branch formula (avoids cancellation and overflow):

$$\text{net} = N_D - N_A, \qquad \text{disc} = \sqrt{\text{net}^2 + 4,n_i^2}$$

$$\text{n-type (net} \ge 0\text{):}\quad n = \tfrac{1}{2}(\text{net} + \text{disc}),\quad p = n_i^2 / n$$

$$\text{p-type (net} < 0\text{):}\quad p = \tfrac{1}{2}(-\text{net} + \text{disc}),\quad n = n_i^2 / p$$

These values are computed once per experiment in build_material_arrays() and stored as n_L, p_L, n_R, p_R. With default ohmic contacts, the time derivatives at the contact nodes are set to zero (dn[0] = dn[-1] = dp[0] = dp[-1] = 0) so the Dirichlet values remain constant throughout the transient.

In FULL mode, optional selective-contact coefficients replace the ohmic pin for the configured carrier/side with a Robin-type boundary flux:

$$J_{c,s} = \sigma_{c,s},qS_{c,s}(u_c - u_{c,\mathrm{eq}})$$

where $c \in {n,p}$, $s \in {\text{left},\text{right}}$, $u_n=n$, and $u_p=p$. The YAML schema supports both flat keys (S_n_left, S_p_left, S_n_right, S_p_right) and a nested readable block:

device:
  mode: full
  contacts:
    left:
      S_p: 1.0e3
      S_n: 1.0e-3
    right:
      S_n: 1.0e3
      S_p: 1.0e-3

Missing or null means the default ohmic pin remains active for that carrier/side. A finite value activates the Robin flux; S = 0 is blocking, and large S approaches the ohmic limit. The 2D solver maps left/right onto top/bottom contacts.

Source: perovskite_sim/solver/mol.py, perovskite_sim/physics/contacts.py, perovskite_sim/twod/solver_2d.py


Ion (Vacancy) Densities — Neumann (Zero-Flux)

Ions cannot leave the device. At both contacts the vacancy flux is set to zero:

$$F_P(x = 0) = F_P(x = L) = 0$$

Implemented by padding the internal flux array with zeros at each end before computing the finite-difference divergence. The same zero-flux condition applies to both positive and negative ion species.

A finite-site modified-PNP chemical potential supplies the default crowding term. For one positive species,

$$ F_P=-D_{\text{ion}}\left[ \frac{1}{1-P/P_{\text{lim}}}\frac{\partial P}{\partial x} +\frac{P}{V_t}\frac{\partial\varphi}{\partial x} \right]. $$

The dual-ion default uses shared occupancy $\theta=(P_+ + P_-)/P_{\text{lim}}$. The legacy tier retains the former whole-flux multiplier for frozen benchmark compatibility. P_lim enters the chemical potential; it is not a numerical projection of every trial state.

Source: perovskite_sim/physics/ion_migration.py


2D Extension Boundary Conditions

The 2D transport contract records the carrier-transport assumptions, field-mobility and reabsorption closures, and the scope of the historical uniform-device parity checks.

The Stage A/B 2D solver uses the same physical stack but extrudes it onto a tensor-product grid. The vertical stack direction carries the device-contact boundary conditions: top/bottom carrier rows are ohmic by default, or Robin when FULL-mode selective-contact coefficients are configured. The lateral direction supports periodic boundaries for laterally uniform presets and Neumann symmetry boundaries for finite-domain microstructure studies. The latter can add finite-volume vertical grain boundaries whose physical band width is preserved by sub-cell SRH mixing rather than whole-node lifetime painting; non-empty microstructures currently fail closed on the periodic topology. The default 2D J-V path keeps positive and negative ions as static Poisson background fields and omits the 1D interface-defect/interface-plane channels. That compatibility path is certified only in the lateral-uniform frozen-ion limit and for prescribed lifetime patterns. A separate explicit research path supports one blocking positive ion and a restricted cross-node interface-SRH sheet; it does not turn the default path into a complete 2D perovskite microstructure model.

At solver level, an explicit ion_dynamics="single_mobile" research transient adds one positive-ion state on Neumann-x domains with blocking boundaries and discrete inventory checks. An explicit companion evaluator reconstructs electron, hole, ionic, and instantaneous displacement current from the same semidiscrete RHS and constitutive face laws. A matching explicit jv-2d-execution-protocol-v1 opens the research-only Python/backend J-V path; the workstation and historical carrier-only terminal-current API remain on the frozen compatibility lane. See docs/TwodMobileIonTransient.md and docs/TwodMobileIonCurrent.md, plus docs/TwodJvExecutionProtocol.md.

An independent interface_srh="two_sided_cross_node" research opt-in maps the 1D two-sided interface-defect surface rate onto a conservative Neumann-x 2D sheet sink and composes with the mobile-ion state. It explicitly excludes projection, shared occupancy, QSS/dynamic interface states, and interface charge. It composes with the mobile-ion state in the same strict public research protocol. See docs/TwodTwoSidedInterfaceSrh.md. The exact combined slice has a source-clean 3x3 grid/tolerance certificate; see docs/TwodCombinedNumericalCertificate.md.

Historical 2D presets now live under tests/fixtures/configs/twod/. The Python and backend entry points (kind="jv_2d" and kind="voc_grain_sweep") remain available for explicitly defined configurations within their supported topologies; historical tests still require their original inputs.

Source: perovskite_sim/twod/solver_2d.py, perovskite_sim/twod/mobile_ion_current_2d.py, perovskite_sim/twod/interface_recombination_2d.py, perovskite_sim/twod/experiments/jv_protocol_2d.py, perovskite_sim/twod/experiments/jv_sweep_2d.py


Thermionic Emission at Heterointerfaces

The default density-variable path retains one bidirectional Scharfetter-Gummel face and may apply the empirical thermionic cap above a 0.05 eV band offset. The dimensionally normalized te_physical_norm form remains opt-in.

The supported ion-free QF path can instead use interface_boundary=true. It removes the ordinary SG face and solves reciprocal thermionic transport plus shared-occupancy interface SRH on an exclusive zero-thickness boundary. Unsupported physics stops before Newton. The registered CBO grid envelope passes, but external SCAPS shape agreement does not; this is a development model rather than an externally certified CBO threshold.

Source: perovskite_sim/physics/continuity.py, perovskite_sim/discretization/fe_operators.py


Interface Recombination

At each heterointerface, surface recombination is parameterised by velocities $(v_n, v_p)$ [m/s] carried in DeviceStack.interfaces. The surface SRH rate is converted to a volumetric rate by dividing by the local dual-grid cell width.

These active interface paths are recombination-only. The explicit interface_charge_closure schema defaults to off; its reserved equilibrium_referenced research value is recognized but fails closed in all production experiment routes. The retired iface_state_charge scalar cannot feed the shared-node Poisson scaffold. A dedicated ion-free QF research lane does couple the signed law Delta sigma = -q Nt (f-f_eq) into the outer two-sided Poisson residual and analytic/IFT Jacobian. It is available through the two-call Python API and the fail-closed POST /api/research/interface-charge/steady-state endpoint only; it is not a production trap-electrostatics capability. D4-E2 has current-source internal grid/tolerance certificates for the charge-off, charged, and resolved device-stress lanes; the production gate remains closed pending D4-E3. See docs/InterfaceChargeClosurePolicy.md.


Initial Conditions

Quasi-Neutral Dark Seed (Default)

solve_equilibrium constructs a quasi-neutral carrier seed with a neutral ionic background. It is an initializer, not by itself a residual-certified ion-relaxed equilibrium. At every grid node:

$$n \cdot p = n_i^2(\text{layer}) \quad\text{(mass-action law)}$$

$$n - p = N_D - N_A \quad\text{(charge neutrality, ions treated as neutral background)}$$

The configured vacancy density $P_0$ is treated as a neutral ionic background — it does not appear as net space charge in the initial carrier balance. This avoids the enormous artificial carrier imbalance that arises when $P_0$ is treated as net positive charge.

The ion profile is initialised to the uniform per-layer value P0. Contact nodes are overwritten with the ohmic-contact equilibrium densities.

Source: perovskite_sim/solver/newton.py


Illuminated Carrier Preconditioning

For experiments that begin under illumination (J-V sweep, impedance, degradation), the initial state is obtained by integrating the full MOL system for $t_{\text{settle}} = 1$ ms starting from dark equilibrium, under illumination at the starting voltage:

y_dark  = solve_equilibrium(x, stack)
y_light = run_transient(x, y_dark, [0, t_settle], illuminated=True, V_app)

The default 1 ms integration is a carrier-preconditioning protocol, not a claim of full ion-relaxed equilibrium. Carrier densities typically settle much faster than the ionic profile, so ions remain nearly frozen over this interval. Longer light soaking defines a different device history and can change a history-dependent J-V result. If the transient solve fails or returns an invalid terminal density state, the helper raises; it never substitutes the dark seed as an illuminated result.

Source: perovskite_sim/solver/illuminated_ss.py


Built-in Potential

built_in_potential_mode selects the signed contact potential used by the Poisson boundary:

Mode Source Intended use
semiconductor_work_function Endpoint semiconductor work functions from $\chi$, $E_g$, DOS, doping, grading, and temperature New physical stacks whose outer layers represent the contact reservoirs
metal_work_function Explicit work_function_left_eV - work_function_right_eV Devices with known electrode work functions
legacy_manual Non-negative V_bi_override magnitude with orientation inferred separately Frozen published benchmarks; V_bi remains a deprecated compatibility alias

The fail-closed semiconductor mode requires complete positive band/DOS inputs. Carrier contact kinetics remain independent: Dirichlet pins or finite-rate Robin S_* values do not select the Poisson voltage. With $s=\operatorname{sign}(V_{\text{bi}}^{\text{bc}})$, positive forward bias is mapped as

$$ \varphi(0)=0,\qquad \varphi(L)=V_{\text{bi}}^{\text{bc}}-sV_{\text{app}}. $$

Pre-mode YAML files retain their historical compatibility behavior. New files that omit both manual keys resolve to semiconductor_work_function.

assess_contact_thermodynamics(stack, mat) evaluates the four endpoint Maxwell-Boltzmann quasi-Fermi levels using the exact reservoir densities, band/DOS arrays, device temperature, and Poisson drop consumed by the solver. The fixed internal gate is 5 meV and callers may tighten but not relax it. metal_work_function additionally requires each electrode work function to match its local semiconductor reservoir in absolute energy; matching only the left-right difference is insufficient. Results are labelled certified, inconsistent, compatible_unverified, or not_assessable. Legacy decks are reported rather than silently migrated or rejected; research workflows can use require_contact_thermodynamic_certificate to fail closed. This is an internal boundary-consistency certificate, not validation of a real contact barrier or surface kinetics.

Source: perovskite_sim/physics/contacts.py


J-V Sweep Voltage Range

run_jv_sweep(stack, V_max=None, ...) opens the forward sweep to

$$V_{\text{upper}} = \max\bigl(V_{\text{bi,eff}} \cdot 1.3,; 1.4\ \text{V}\bigr)$$

where $V_{\text{bi,eff}}$ is the magnitude returned by stack.operating_built_in_potential(). Explicit modes use their selected work-function source consistently; pre-mode compatibility stacks retain the historical band-derived operating value. The 1.3× headroom and 1.4 V floor avoid silently clipping a high-$V_{\text{oc}}$ stack before $J=0$. compute_metrics still exposes voc_bracketed=false when the sampled window does not resolve open circuit.

Source: perovskite_sim/experiments/jv_sweep.py::_default_V_max


External Series and Shunt Resistance

Series resistance and shunt leakage are modeled in an independent, area-normalized DC circuit layer rather than as contact-calibration knobs. For photovoltaic output current positive, it maps an intrinsic junction curve as

$$ J_{\mathrm{terminal}}=J_{\mathrm{device}}-\frac{V_j}{R_{\mathrm{sh}}}, \qquad V_{\mathrm{terminal}}=V_j-J_{\mathrm{terminal}}R_s. $$

ExternalCircuitProtocol freezes the topology, sign convention and resistance values; apply_external_circuit retains both junction and terminal arrays, balance residuals, source/protocol/mapping hashes, and adjusted metrics. The zero-coupling default preserves the source curve and 1-sun metrics exactly. Folded terminal curves and uncertified source results fail closed. The separate POST /api/jv/external-circuit route exposes this opt-in result without changing /api/jv or the workstation default. This is a DC post-processing layer, not yet a terminal-voltage-driven transient circuit DAE. See docs/ExternalCircuitDc.md.

The source-clean external-series-shunt-dc-operating-quadrant-v2 lane is internally certified at commit 2392ba3 with 9/9 grid/tolerance cells and certificate SHA-256 a9f6d63a229ec613594d78d73a2ac94e6d0aea756c10e960948dc31123e1bf26. It compares both branches on fixed $V_{terminal}/V_{oc}\in[0,1]$ coordinates, while separately certifying terminal $V_{oc}$, $J_{sc}$, FF, PCE, exact Kirchhoff balances, source/protocol provenance, and the bit-exact zero-coupling limit. The earlier full-junction-trace v1 lane remains partial because its deep-forward-injection tail is not grid converged; that tail is explicitly not covered by the v2 certificate. The resistance values are synthetic stress inputs, so this remains internal numerical evidence rather than measured-device validation.

Source: perovskite_sim/experiments/external_circuit.py


Lumped Thermal Energy Balance

LumpedThermalProtocol defines an explicit area-normalized first-law control volume around the device plus lumped parasitics. solve_lumped_thermal_steady_state and run_lumped_thermal_transient account for absorbed optical power, terminal electrical export, optional internal heat, linear heat rejection, and one-sided far-field radiation. Canonical protocol hashes, immutable cumulative energy ledgers, maximum-temperature envelopes, and recomputed residual checks make invalid or altered evidence fail closed.

The independent transient keeps terminal power constant. The separate opt-in solve_electrothermal_operating_point path closes a steady feedback root by re-running a fresh, strict temperature-dependent J-V experiment, applying the frozen external circuit, and feeding its sampled terminal MPP export into the thermal balance. It retains every temperature trial, full experiment protocol, source/mapping hashes, and first-law residual. This is not yet a joint electrical-thermal transient DAE or a spatial heat equation. Absorbed optical power remains an explicit input and is not inferred from incident irradiance. See docs/ThermalEnergyBalance.md.

The source-clean electrothermal-terminal-mpp-grid-resolved-v3 lane is internally certified at commit 6afb1a6 with 9/9 grid/tolerance cells and certificate SHA-256 67d18979c9ce2d23a0e2d0e848513f50d8e1e4880819c92a5ff4da8719682cf6. Its 60-to-80 terminal-grid differences are 0.7079% for the 300--380 K sampled-MPP power response, 0.00151 K for operating temperature, 0.0160% for terminal MPP current, 0.0141% for power, and 12.37 microvolts for voltage. The 10/15/20 v1 and 20/30/40 v2 lanes remain partial; their preserved failures show that lower grids had not reached the asymptotic regime. This certificate is limited to the frozen synthetic thermal/circuit inputs and fresh-state sampled-MPP protocol, not measured-device temperature validation.

Sources: perovskite_sim/experiments/thermal_balance.py, perovskite_sim/experiments/electrothermal.py


Interface-SRH Identifiability

The opt-in synthetic identifiability workflow evaluates the production interface-SRH rate formula together with the research equilibrium-referenced sheet-charge closure. A frozen protocol binds carrier conditions, parameter bounds and truth, observable noise scales, finite-difference controls, multi-start fitting, Fisher/SVD rank, fixed-parameter profiles, and the fail-closed forward-error policy. Results report analysis_certified and parameters_identifiable separately so a successfully diagnosed rank deficiency cannot be presented as a precise parameter estimate. The CLI scripts/run_interface_srh_identifiability.py and strict POST /api/identifiability/interface-srh-synthetic route expose the same canonical result contract. See docs/InterfaceSrhIdentifiability.md.

The source-clean interface-srh-identifiability-synthetic-v1 lane is internally certified at commit 12fc7cc with 9/9 carrier-grid/finite-difference cells and certificate SHA-256 b5fd5f2d784277c49b0b2720ad7f225b76586945bd2c9d2c772a34fa3dec5643. With trap density, capture scale, and calibration all free, every cell reports rank 2/3 and the expected capture/calibration null direction; no parameter claim is made. Holding capture scale fixed gives rank 2/2 and exact recovery of the synthetic trap-density and calibration truth. This is formula-local, noise-free internal numerical evidence, not full-device, measured-data, or materials-parameter validation.

Source: perovskite_sim/experiments/identifiability.py


Voc Bracket Detection

compute_metrics(V, J) returns a frozen JVMetrics(V_oc, J_sc, FF, PCE, voc_bracketed). The voc_bracketed flag is True when the sweep contains a zero-current crossing (Voc is then computed by linear interpolation between the two adjacent samples whose J flip sign) and False when the sweep window stops before J flips sign — in which case V_oc, FF, and PCE are zeroed sentinels and only J_sc is physically meaningful. The 2D driver run_jv_sweep_2d calls compute_metrics(..., assume_jsc_positive=False) so the 2D Scharfetter–Gummel sign convention is normalised centrally and JV2DResult.metrics matches the 1D semantics bit-for-bit. The workstation 2D J–V pane reads voc_bracketed to render — for Voc / FF / PCE and an inline V_oc not bracketed — increase V_max warning. The 1D and 2D J–V main plots always use the operational range, without a range selector. A display-only y-axis clip $[-0.15,J_{\text{sc}},\ +1.12,J_{\text{sc}}]$ (mA/cm²) and an upper voltage limit of $V_{\text{oc}} + 0.18$ V keep the photovoltaic region readable; the voltage limit never exceeds the largest sample plus 0.05 V. The 1D limits account for both branches. If valid bracketed metrics are unavailable, the y-axis uses autorange and the voltage cap is omitted. All raw V/J samples and computed metrics are preserved.

Source: perovskite_sim/experiments/jv_sweep.py::compute_metrics, frontend/src/workstation/panes/main-plot-pane.ts::renderJV2D


Summary Table

Variable Contact BCs Type Source
$\varphi$ $\varphi(0)=0$, $\varphi(L)=V_{\text{bi}}^{\text{bc}}-sV_{\text{app}}$ Signed Dirichlet models/device.py, solver/mol.py
$n$ default $n(0) = n_L$, $n(L) = n_R$; optional FULL-mode Robin flux for configured sides Dirichlet or Robin solver/mol.py, physics/contacts.py
$p$ default $p(0) = p_L$, $p(L) = p_R$; optional FULL-mode Robin flux for configured sides Dirichlet or Robin solver/mol.py, physics/contacts.py
$P$ (ions) $F(0) = F(L) = 0$ Neumann physics/ion_migration.py
$P^-$ (neg ions) $F(0) = F(L) = 0$ Neumann physics/ion_migration.py


Python-Only Quick Start

from perovskite_sim.models.config_loader import load_device_from_yaml
from perovskite_sim.experiments.jv_sweep import run_jv_sweep

stack = load_device_from_yaml("configs/calado2016_fig1f.yaml")
result = run_jv_sweep(stack, N_grid=100, n_points=61, v_rate=0.04, V_max=1.2)
print(f"PCE: {result.metrics_fwd.PCE*100:.2f} %")

Dark J-V

result = run_jv_sweep(stack, N_grid=80, n_points=40, v_rate=1.0, illuminated=False)
# result.J_fwd is the diode injection-current characteristic

Current Decomposition

from perovskite_sim.experiments.jv_sweep import compute_current_components

cc = compute_current_components(x, y_state, stack, V_app=0.5, mat=mat)
# cc.J_n, cc.J_p, cc.J_ion, cc.J_disp, cc.J_total — all shape (N-1,)

Spatial Profiles

result = run_jv_sweep(stack, N_grid=80, n_points=40, v_rate=1.0, save_snapshots=True)
snap = result.snapshots_fwd[20]  # snapshot at the 20th voltage point
# snap.phi, snap.E, snap.n, snap.p, snap.P, snap.P_neg, snap.rho

Transient Photovoltage (TPV)

from perovskite_sim.experiments.tpv import run_tpv

result = run_tpv(stack, N_grid=80, delta_G_frac=0.05, t_pulse=1e-6, t_decay=50e-6)
print(f"V_oc: {result.V_oc:.4f} V")
if result.tau is not None:
    print(f"Photovoltage relaxation time: {result.tau:.3e} s")
else:
    print(result.fit.status if result.fit is not None else "Fit unavailable")

The reported pulse is result.delta_V, relative to the matched unpulsed result.V_reference. The open-circuit voltage evolves with terminal charge; current continuity and integrated charge are checked throughout the trace. A missing or non-single-exponential decay does not acquire a lifetime. See TpvOpenCircuitContract.md.


Phase 2 Characterisation Experiments

A family of four higher-level experiments that wrap the core drift-diffusion solver and produce characterisation-grade fits — dark-diode ideality, intensity-dependent $V_{\text{oc}}$, wavelength-resolved quantum efficiency, and capacitance-voltage Mott-Schottky analysis.

Each wrapper returns a frozen dataclass that carries both the raw sweep and the derived fit. All four accept an optional progress=lambda stage, k, n, msg: ... callback so they plug into the backend's SSE progress stream unchanged.

Dark J-V with ideality + $J_0$

from perovskite_sim.experiments.dark_jv import run_dark_jv

r = run_dark_jv(stack, V_max=1.2, n_points=60)
# r.V, r.J                      — dark forward sweep (G = 0 everywhere)
# r.n_ideality                  — diode ideality factor from log|J| vs V
# r.J_0                         — saturation current density [A/m²]
# r.V_fit_lo, r.V_fit_hi        — auto-selected exponential-region window

The fit window is chosen as the voltage range where $\log |J(V)|$ is most linear — the selector minimises the mean $|d^2 \log|J| / dV^2|$ on a 6-point sliding window, then grows outward while the neighbour residual stays within 10× of that minimum. This rejects the sub-turn-on leakage tail (where the $-1$ in $J_0 (e^{V/nV_T} - 1)$ dominates) and the high-V series-resistance roll-off.

Source: perovskite_sim/experiments/dark_jv.py


Suns–$V_{\text{oc}}$ with pseudo-JV

from perovskite_sim.experiments.suns_voc import run_suns_voc

r = run_suns_voc(stack, suns_levels=(0.01, 0.1, 1.0, 5.0, 10.0))
# r.suns, r.V_oc, r.J_sc        — per-intensity sweep
# r.J_pseudo_V, r.J_pseudo_J    — Sinton pseudo-JV: V = V_oc(X),
#                                                    J = J_sc_ref − J_sc(X)
# r.pseudo_FF                   — series-resistance-free fill factor

At each illumination level the wrapper scales the cached 1-sun G_optical profile, solves for the illuminated steady state at $V = 0$ (to read $J_{\text{sc}}$), then bisects for $V_{\text{oc}}$. The resulting $V_{\text{oc}}(X)$ vs $\ln X$ slope gives the effective ideality factor; the constructed pseudo-JV curve (Sinton convention) is immune to series resistance because every point is measured at $V_{\text{oc}}$ where $J = 0$.

Source: perovskite_sim/experiments/suns_voc.py


EQE / IPCE

import numpy as np
from perovskite_sim.experiments.eqe import compute_eqe

wavelengths = np.linspace(350.0, 850.0, 20)
r = compute_eqe(stack, wavelengths_nm=wavelengths)
# r.EQE                         — per-λ external quantum efficiency (signed)
# r.J_sc_per_lambda             — J_sc at each probe wavelength [A/m²]
# r.J_sc_integrated             — q · ∫ EQE(λ) · Φ_AM15G(λ) dλ [A/m²]
# r.J_spread_max                — interior face-to-face current spread [A/m²]

At each wavelength, the wrapper builds a single-wavelength TMM generation profile (absorption $A(x; \lambda)$ scaled by Phi_incident), solves the drift-diffusion problem at $V = 0$ under that monochromatic source, and reads $\text{EQE}(\lambda) = \Delta J(\lambda) / (q \cdot \Phi_{\text{inc}})$ from the signed incremental current against the dark baseline. The sign is kept rather than taking a magnitude, so a wrong-way photocurrent (contact polarity, sign convention, or an unconverged dark baseline) surfaces as a negative EQE plus a RuntimeWarning instead of being disguised as a physical curve. A ValueError is raised if no layer carries tabulated optical_material data — Beer-Lambert-only stacks cannot produce a wavelength-resolved EQE. The integrated $J_{\text{sc}}$ cross-checks against a full-spectrum TMM run to within ~25 % on a 15-point wavelength grid. J_spread_max reports how far the settled states were from the uniform $J(x)$ that charge conservation implies at steady state; it is a magnitude diagnostic, not a convergence threshold (see perovskite-sim/CLAUDE.md).

Source: perovskite_sim/experiments/eqe.py


Impedance Protocol and Evidence

run_impedance exposes three deliberately separate engines. The transient_ion_aware path retains mobile ions and lock-in extraction; the qf_frequency_ion_free path retains the QF solver's DC residual and frequency-domain linear-solve diagnostics but rejects mobile ions; and ion_aware_frequency_certified starts from a residual-certified mobile-ion DC state and returns the reference frequency-domain response. All three use a strict perturbation bound below 20 mV.

Every result now carries the exact bias/light/cycle protocol, a DC operating point report, an electrical-grid assessment, and a frequency-window assessment. For an ionic device the assessment reports Debye, blocking-charge, and diffusion frequency estimates. It distinguishes merely bracketing the blocking scale from bracketing the full diffusion/blocking/dielectric envelope and from covering that envelope with one-decade margins and no sampling gap above 0.5 decades. A single marker frequency or two sparse endpoints do not count as branch coverage; recommendations never rewrite requested points.

The certified ion-aware frequency-domain route additionally emits one certificate per frequency instead of only sweep-wide extrema. Each point retains its all-face spread, linear-solve diagnostics, both finite-difference step comparisons, ionic inventory response, current decomposition, and carrier/ion storage components. Its canonical grid-ladder protocol hashes the exact coordinates and DC state on at least three strictly increasing meshes; the finest pair must agree within 2% in impedance magnitude and 1 degree in phase at every frequency. The public result retains the exact DC and frequency protocols, their hashes, the DC-state hash, the complete per-frequency evidence, and separate numerical, frequency-window, grid, and contact-thermodynamic axes. Set require_frequency_window_certificate=True to reject an uncovered ionic frequency range.

The transient engine applies a continuous sinusoidal boundary voltage in one Radau solve per frequency. Edge and midpoint states co-locate centered displacement current with midpoint conduction current; points_per_cycle (default 40) is recorded for time-resolution ladders. Non-finite DC evidence and unsupported dynamic interface-state blocks fail before AC extraction. The historical 1 ms transient preconditioner remains available, but its carrier/ion residuals and contact certificate are reported rather than assumed to prove steady state. Set require_operating_point_certificate=True to fail closed. Candidate residual thresholds still require lane-specific grid/tolerance/amplitude/cycle refinement before a result is internally certified ionic spectroscopy.

Source: perovskite_sim/experiments/impedance.py


Mott-Schottky C-V

import numpy as np
from perovskite_sim.experiments.mott_schottky import run_mott_schottky

r = run_mott_schottky(
    stack,
    V_range=np.linspace(-0.3, 0.4, 8),
    frequency=1e6,       # 1 MHz default; certify a frequency plateau for claims
    # impedance_method="quasi_fermi_frequency",  # audited local ion-free QF only
)
# r.V, r.C, r.one_over_C2       — dark C-V sweep [F/m² and m⁴/F²]
# r.V_bi_fit                    — apparent built-in voltage from 1/C² fit
# r.N_eff_fit                   — ionised-dopant density from slope [m⁻³]
# r.V_fit_lo, r.V_fit_hi        — auto-selected linear window
# r.eps_r_used                  — ε_r taken from the 'absorber'-role layer

A thin wrapper over run_impedance that drives a single AC excitation at frequency at each DC bias and reads capacitance off as $C = \text{Im}(1/Z) / \omega$. A non-positive susceptance is rejected instead of being hidden by an absolute value. Runs dark (illuminated=False) so photogenerated carriers do not screen the depletion capacitance. The linear-fit helper finds the widest contiguous $(V, 1/C^2)$ window whose RMS residual is within 1 % of its ordinate span — rejects the low-bias fully-depleted tail and the high-bias injection tail without a hand-tuned cutoff. On a clean Mott-Schottky curve the p-n-junction fit adds the two-edge thermal correction $2k_BT/q$ to the bare V-axis intercept. Synthetic-data regression tests pin recovery of $V_{\text{bi}}^{\text{app}}$ to $&lt;0.01$ V and $N$ to $&lt;0.02$ decades. The API field remains V_bi_fit for compatibility, but the result is an apparent depletion-model parameter, not an independent measurement of the contact-potential barrier. Flat or positive-slope $1/C^2$ data return an unidentifiable fit rather than a finite but physically meaningless parameter pair. A transient-path physical claim also requires grid, frequency, amplitude, and cycle convergence; a frequency-domain claim instead requires its registered linearization-step and all-face current certificates.

The default impedance_method="transient" retains the general time-domain model. The explicit quasi_fermi_frequency method instead linearizes about a residual-certified dark QF state and requires a nominal perturbation strictly below 20 mV. It currently supports only the audited local ion-free QF subset; mobile ions, selective contacts, thermionic interfaces, and non-local photon recycling fail closed. Its c-Si N=200/300/400 regression recovers the depletion capacitance scale, frequency/grid plateaus, and the independent DC electron and hole inventory derivatives. The corrected finest-grid fit gives $V_{\text{bi}}^{\text{app}}=0.782$ V, N_eff=9.554e21 m^-3, and a 0.111 V gap to the configured 0.893 V contact-potential magnitude. That gap is consistent with the published distributed-carrier p-n intercept range, but no compatible pointwise external curve is frozen. This is an internal numerical certificate, not external C-V validation or a repair of the legacy endpoint-sampled path.

Source: perovskite_sim/experiments/mott_schottky.py


Validation and Model Scope

The current evidence is intentionally split by claim. See the reproducibility registry for commands and the 2026-08-11 manual (distributed separately) for the full traceability matrix.

Registered c-Si QF J-V and C-V grid-ladder observations

This is internal convergence evidence for the restricted local QF driver, not external c-Si device validation and not a transient-driver certificate.

Physical-interface CBO response, grid contraction, and certification gates

The physical-interface CBO campaign passes the registered numerical grid gate but fails the declared SCAPS-shape gate (certified=false).

Registered 1D and 2D parity domain and current model scope

The 1D/2D parity claim covers the registered interface-free, frozen-ion domain. Mobile-ion dynamics and the 1D interface-SRH/physical-QF boundary are not part of that comparison.

Calado 2016 Fig 1e/1f reproduction: forward and reverse J-V at 40 mV/s with and without contact SRH, and figures of merit against the paper

The Calado 2016 Fig 1e/1f lane (configs/calado2016_fig1f.yaml, the paper's SI Table 1 toy stack; scripts/plot_calado_fig1f.py, the paper's −1 → +1.2 V / 3 s hold protocol at 40 mV s−1) is a partial external comparison. The no-contact-SRH control loop closes (HI 0.007 vs 0.00) and the hysteretic reverse branch is quantitative (Jsc 16.0 vs ≈ 16 mA cm−2, Pmax 81 vs ≈ 85 W m−2, Voc 0.78 vs ≈ 0.73 V), but the forward-scan collapse is about 2× too shallow: HI 0.44 against the paper's 1.84 in its Pmax,rev/Pmax,fwd − 1 definition. The cause is open; the preset header records the study.

Scan-rate ladder on the Calado 2016 Fig 1f preset: J-V loops at three rates and the hysteresis-index bell against scan rate

A nine-rate ladder on the same preset (scripts/plot_calado_fig1f_scan_rate.py, same protocol without the +1.2 V hold) gives the scan-rate bell: HI 0.001 at 1 mV s−1, a peak of 39 at 0.3 V s−1, and 0.06 at 100 V s−1, where ions stay frozen in the short-circuit state and both branches collapse together. The paper's 1.84 at 40 mV s−1 falls between the 0.03 (0.28) and 0.1 V s−1 (5.8) rungs, so the residual gap reads as a ≈2× shift of the ionic time scale; it stays open.

Separately, twod-mobile-ion-interface-srh-v1 certifies one synthetic, research-only Neumann-x slice combining a finite-width grain boundary, one blocking positive ion, clamp-inactive cross-node InterfaceDefect SRH, and complete carrier/ion/displacement current. Its source-clean 9-cell certificate is b02bc4f8b3b5d470d599f6dacde746b26c263591aafecd14cc6c890a94b677dd; it is not external solver or measured-device validation.

Phase 4.1 now includes five certified research-only semi-explicit DAE slices. The first three are single-layer, ohmic, no-interface topologies: no ions, one blocking positive ion, and blocking positive/negative ions on a shared finite-site lattice. The fourth is a two-layer, ion-free, charge-off topology with one physical interface and four algebraic Fermi-Richardson trace states. The fifth combines one synthetic blocking positive ion with those four algebraic traces in the same two-layer, charge-off topology. All retain Poisson potential as an algebraic coordinate, construct a residual-certified consistent initial condition, and provide dense reference and sparse analytic Newton paths. The ion slices add bounded coordinates, production blocking flux/tangents, per-species dual-cell inventory, and exact ion-storage derivatives. The interface slice adds exact adjacent bulk/state cross-block tangents for projection, reciprocal exchange, and shared-occupancy SRH, and fails closed if any defining clamp activates. The five registered 9-cell matrices are internally certified at source commits 985a234, 6e9a274, 2d6b32f, 008aef3, and dc540ca (certificates 44807d654d..., 7538fa4ace..., 15a6a4dcf..., 21bb12e465..., and 02f9e0b2bd...). The dual and combined lanes' ion parameters are synthetic protocol inputs. The combined certificate covers a source-clean 9-cell grid/time-step matrix: all cells completed, maximum carrier/ion/interface normalized residuals were 4.97707e-8 / 1.54565e-10 / 2.29723e-14, maximum ion inventory drift was 9.69749e-14, and the sparse path used at most 0.01935 of the dense residual work. It excludes InterfaceDefect, configurable cross-node carrier sampling, dynamic/charged or two-sided interface states, dual ions, selective contacts, field mobility, photon recycling, clamp-active points, experiments, and backend routes. These are internal numerical certificates, not external solver or experimental validation; see the DAE capability boundary.

Phase 4.3 adds a separate, research-only energy-distributed bulk-trap equilibrium closure. A strict standard-SI schema defines integrated trap density, energy distribution, capture cross sections, thermal velocity, and the donor/acceptor neutral-charge reference. Occupancy, SRH recombination, and absolute trap charge share one controlled energy quadrature, and the Poisson Newton solve uses the analytic trap-charge tangent. Default MoL and experiment routes fail closed when this field is active; historical lifetime profiles and SCAPS Gaussian metadata are not reinterpreted. See the P4.3 capability contract.

Phase 4.4 adds an opt-in composition-resolved CIGS optical path. Minoura Tauc-Lorentz reference spectra and critical-point shifts construct Kramers-Kronig-related n,k; Carron's independent absorption expression is a three-composition benchmark. The exact electrical grading coordinate drives the TMM slice compositions, and the shared adapter is used by build-once generation, photon recycling, EQE, EL, and tandem optics. Historical graded CIGS presets remain scalar-optics by default. See the P4.4 CIGS optics contract.