Why a Lagrangian kernel exists.

Why we have fifty physics engines

The engine on the landing page would be three separate codebases anywhere else. Rigid bodies live in PhysX or Bullet. Molecular dynamics lives in LAMMPS or GROMACS. RNA folding lives in ViennaRNA or Rosetta. Fluids live in OpenFOAM or LiquidFun. None of them compose.

The split isn't accidental — it's the structure each tool optimised against. PhysX's contact solver assumes a contact graph, not pairwise potentials over 10⁵ particles. LAMMPS's neighbour list assumes short-cutoff non-bonded forces, not deformable continua. ViennaRNA's dynamic-programming kernel doesn't accept a force field at all. Each engine made a different cut for a different audience, and the cuts don't line up.

So adding fluids to MuJoCo means reimplementing fluids. Adding EM coupling to an RNA scene means rewriting ViennaRNA. Adding deformables to a PhysX game means swapping engines mid-pipeline. The cross-domain tax — porting, validating, re-tuning — usually exceeds the cost of the physics itself.

A Lagrangian doesn't care which subset you pick. HarmonicBond and LennardJones and QcdCoulomb are the same kind of node in the same kind of graph, evaluated by the same kernel. Choosing which terms run for a given scene becomes a scene-level declaration — not a tooling decision. The live 1000-block stack on the landing page runs on the same evaluator the LJ fluid and the pendulum chain run on. That's the whole pitch.

See the live 1000-block stack at index.html#demos — same evaluator, same scene file. uses newtonian.classical_mechanics.avbd; render "rigid_body".

Subsystems as units of compute

.lagra files describe nodes in the graph. Like Lean for proofs, each file is a self-contained unit the compiler can check, optimise, and compose. Coarse-graining is a node choice, not an engine choice.

The scene below is the harmonic-oscillator demo from the landing page — verbatim, no shortened-for-the-website magic. Two particles, one bond, one declared invariant, one solver block:

scene solver_walkthrough_oscillator {
    uses newtonian.classical_mechanics
    description "Two-bead harmonic oscillator. Only the harmonic bond contributes nonzero force."
    wasm true                                // flag for the landing-page bundle

    entity { type: ball_small, position: [-1.35, 0, 0], mass: 1.0, size: 0.35 }
    entity { type: ball_small, position: [ 1.35, 0, 0], mass: 1.0, size: 0.35 }

    bond harmonic [0, 1] { k: 2.0, r0: 2.0, damping: 0.04 }
    @expect momentum_conservation                  // CI gate
    solver { dt: 0.002, steps_per_frame: 12 }
}

import lagra

s = lagra.Scene(
    name="solver_walkthrough_oscillator",
    uses="newtonian.classical_mechanics",
    description="Two-bead harmonic oscillator. Only the harmonic bond contributes nonzero force.",
)

s.add_entity(particle="ball_small", position=[-1.35, 0, 0], mass=1.0, size=0.35)
s.add_entity(particle="ball_small", position=[ 1.35, 0, 0], mass=1.0, size=0.35)

s.add_bond(0, 1, k=2.0, r0=2.0, damping=0.04)
s.expect(["momentum_conservation"])           # CI gate, same validator catalog
s.set_solver(dt=0.002, steps_per_frame=12)
sim = s.build()

The compiler checks before the kernel runs. Conservation invariants are derived from the active graph and inserted as runtime assertions — if the first step would violate them, the scene doesn't compile. Entity catalogs enforce scientific references: entity ball_small {…} in configs/entities/classical.lagra carries the citation alongside the mass, size, and color, and custom entities go through LagraDef validation against those catalogs.

Each node also owns its own optimisation hints. Pairwise terms can request a cell list; bonded terms can request a neighbour list; whole subgraphs can target a GPU compute pass. The planner reads the hints and composes the dispatch — the kernel's shape is fixed but the dispatch is per-scene.