The Porting Process

spice-rs was built by systematically reading ngspice C source code and translating it to Rust. Not reimplementing, not approximating -- translating. The distinction matters.

SPICE has 50 years of accumulated numerical tricks. Every if-statement, every magic constant, every seemingly redundant check exists because someone hit a real circuit that broke without it. The ngspice codebase encodes the hard-won solutions to thousands of convergence edge cases, numerical stability issues, and device physics subtleties. Trying to "improve" or "simplify" this code without understanding every line is how simulator ports fail.

This chapter documents the methodology that achieved 200/224 test circuits passing with 176 producing bit-identical results to ngspice.

Philosophy: Port, Don't Approximate

The core principle of the spice-rs port is a single rule:

Read the actual ngspice C code. Follow the same logic. Don't invent alternative approaches. If your implementation doesn't match ngspice, go back and read more C. The answer is always in the reference source.

Why this matters

SPICE is not a textbook algorithm. It is a textbook algorithm plus 50 years of patches, fixes, workarounds, and numerical tricks accumulated by the Berkeley team, the ngspice maintainers, and the broader SPICE community. The published papers describe the theory. The code describes reality.

Examples of things that exist in the code but not in the papers:

• Voltage limiting in device models. The Shichman-Hodges MOSFET equations are smooth, but Newton-Raphson will overshoot wildly on a junction turn-on without the DEVfetlim / DEVpnjlim voltage limiters.

• The ipass mechanism. When .NODESET is used, ngspice runs the first NR convergence with diagonal elements forcing nodes to their set values (MODEINITFIX). After convergence, it flips to MODEINITFLOAT, removes the forcing, and runs one more iteration. This is not documented anywhere obvious.

• Source stepping fallback. If gmin stepping fails, ngspice switches to source stepping: scale all independent sources from 0 to 1 in numSrcSteps increments.

• Integration order control. The transient engine starts at order 1 (backward Euler) and steps up to order 2 (trapezoidal) after the first accepted step. After a breakpoint, it drops back to order 1.

What this means in practice

Use ngspice variable names. When the C code uses vgs, the Rust code uses vgs. When it uses qgs, we use qgs. This makes side-by-side verification possible.

Preserve control flow. If the C code tests if (mode & MODETRAN) before computing charges, the Rust code tests if mode.is_tran() at the same point.

Keep magic constants. ngspice defines EPSOX = 3.453133e-11 and EPSSI = 1.03594e-10 in the BSIM3 model. These are not the NIST values -- they are the values that the BSIM3 team calibrated against. Using "more accurate" constants breaks parameter extraction.

Comment the C source location. Every function, every significant block, should reference the ngspice file and line number:

// Port of mos1load.c:163-171 -- cutoff region
if vgs <= von {
    gds = 0.0;
    ids = 0.0;
    gm = 0.0;
    gmbs = 0.0;
}

When to deviate

Legitimate reasons to write Rust-idiomatic code:

• Memory safety. Replace raw pointer arithmetic with array indexing.
• Error handling. Replace goto error patterns with Result<T, E>.
• Type safety. Use enums instead of integer flag constants.

But never deviate on numerical behavior. The sequence of floating-point operations must match ngspice.

Investigation-First Method

Every subsystem in spice-rs was ported using the same five-step method.

1. Read the ngspice C code

Read the actual source files -- not the comments, not the man page. For a MOSFET model, this means reading mos1load.c, mos1temp.c, mos1set.c, and mos1defs.h. Read it with a C debugger mindset: follow every pointer, understand every macro expansion, trace every control flow path.

2. Document the algorithm in plain English

Before writing any Rust, write a plain-English description of what the C code does. This catches misunderstandings early. If you can't explain what the C code does in English, you can't translate it to Rust correctly.

3. Identify key data structures and control flow

Map the C data structures to Rust equivalents:

Map the control flow:

4. Translate to Rust

Keep the translation mechanical. Same variable names, same structure, same order of operations. Add comments referencing C file and line numbers.

5. Validate against ngspice output

Run test circuits through both engines and compare. The comparison must be exact for well-conditioned circuits. For complex circuits, they must be within abs=0.01, rel=0.01.

This is constrained translation, not engineering

When something doesn't match:

1. Don't hypothesize. Don't reason about "maybe the NR loop is overshooting."
2. Instrument and run. Add eprintln! tracing to both sides and compare step by step.
3. Read more C. The answer is always in the reference source.

Eval Harness

The spice-eval crate is the validation backbone of the port. It runs test circuits through both spice-rs and ngspice (via FFI), compares results, and reports divergences.

Architecture

spice-eval
|-- src/main.rs          -- CLI, comparison logic, report formatting
+-- eval/
    |-- manifest.toml    -- test circuit registry
    |-- dc/              -- DC operating point circuits
    |-- tran/            -- transient circuits
    +-- ac/              -- AC analysis circuits

CLI modes

Summary mode (default): Runs all 224 circuits and prints a pass/fail table.

cargo run --release --bin spice-eval

Filter mode: Run only circuits matching a substring.

cargo run --release --bin spice-eval -- --filter=mosfet

Diverge mode: Find the first timepoint where divergence exceeds tolerance and show per-device state.

cargo run --release --bin spice-eval -- --diverge="Circuit Name"

Diverge-deep mode: Per-NR-iteration comparison showing RHS vectors, solution vectors, per-device conductances, and stored currents.

cargo run --release --bin spice-eval -- --diverge-deep="Circuit Name"

Parameter check mode: Compare parsed model parameters between engines.

cargo run --release --bin spice-eval -- --check-params

Translate check mode: Compare the TRANSLATE external-to-internal node mapping.

cargo run --release --bin spice-eval -- --check-translate

Case Study: Porting MOSFET Level 1

The Level 1 model illustrates the porting method on a device with moderate complexity: DC I-V equations, voltage limiting, junction diodes, and Meyer charge model capacitances.

spice-rs source: sim/spice-rs/src/device/mosfet1.rs (1111 lines) ngspice source: reference/ngspice/src/spicelib/devices/mos1/mos1load.c, mos1temp.c, mos1set.c

Key challenge: Meyer charge model

The Meyer capacitance model computes gate charges (Qgs, Qgd, Qgb) as functions of the terminal voltages. The capacitances are voltage-dependent and change discontinuously at region boundaries. The Rust translation preserves the exact order of operations from the C code.

Validation results

Level 1 MOSFET circuits match ngspice at machine precision (~1e-14 relative error):

[L3] NMOS Level 1 DC              PASS   7.105e-15   1.017e-14
[L3] NMOS Level 1 Body Effect     PASS   1.421e-14   2.841e-14
[L5] CMOS Inverter                PASS   3.553e-15   2.367e-14

Lessons learned

• Variable naming matters. Keeping the ngspice names (VTO, vto, Von, Vdsat) eliminates ambiguity.
• Mode flags are critical. Getting MODEINITJCT, MODEINITFIX, MODEINITFLOAT, MODETRAN, MODEINITSMSIG wrong produces silent failures.
• Parser validation first. A missing parser case for NSS caused incorrect threshold voltage computation. The --check-params check caught it immediately.

Case Study: Porting BSIM3v3

BSIM3v3 is the most complex device model in spice-rs. The load function alone (b3ld.c) is approximately 5000 lines of C.

Strategy: function-by-function translation

1. Identify logical blocks in b3ld.c (voltage limiting, effective parameters, drain current, output conductance, junction diodes, charge model, integration, matrix stamps).
2. Translate one block at a time.
3. Test after each block.

Constants: use the BSIM3 values

const EPSOX: f64 = 3.453133e-11;    // not 3.45e-11 or eps0 * 3.9
const EPSSI: f64 = 1.03594e-10;     // not eps0 * 11.7
const CHARGE_Q: f64 = 1.60219e-19;  // not 1.602176634e-19
const KB: f64 = 1.3806226e-23;      // not 1.380649e-23

These are the values the Berkeley BSIM team used for parameter extraction.

Lessons learned

• Translate the defaults too. b3set.c contains conditional defaults that must match exactly.
• Watch for sign conventions. BSIM3 internally works with absolute values and applies sign corrections based on NMOS/PMOS type.
• Test with PMOS. NMOS and PMOS exercise different code paths.

Case Study: Porting the Markowitz Sparse Solver

The Markowitz solver is the most structurally complex piece of the port because the C code relies heavily on raw pointers, linked lists, and 1-indexed arrays.

The Rust translation

Arena-based indices instead of pointers. Replace all Element* pointers with u32 arena indices into a Vec<Element>. The sentinel value NONE = u32::MAX replaces NULL.

Preserving 1-indexed convention. Rather than converting to 0-based Rust indexing, sparse-rs preserves the 1-indexed convention. This makes the Rust code line-up with the C code during verification.

The four-level pivot cascade

1. Search for singletons -- rows or columns with exactly one nonzero (zero fill-in).
2. Quick diagonal search -- scan diagonal elements for smallest Markowitz product.
3. Careful diagonal search -- full threshold check against every column element.
4. Full matrix search -- scan every element in the unreduced portion (last resort).

Validation

The Markowitz solver is cross-validated against KLU: both solvers given the same matrix must agree to machine precision. For the spice-rs integration, it is validated indirectly through the full eval harness.