import Base
import ./geom.bend as G
import ./protein.bend as P
import ./topology.bend as T
import ./force.bend as F
import ./rmsd.bend as R
import ./sasa.bend as S
import ./rng.bend as Rng
import ./pdb.bend as Pdb
import ./par.bend as Par
import ./bonded.bend as B
import ./pbc.bend as Pbc
import ./sim.bend as Sim
import ./tpl.bend as Tpl

# Laws: precise specs an AI must prove before merging. Humans write claims here
# and never touch them; PROOF.bend fills each one with a def of the same name.

def empty_atoms() -> List<&2, P.Atom>:
  Nil{}

# Counting contacts in an empty set finds nothing, at any point and cutoff.
law count_below_nil:
  for qp: G.Vec3
  for c2: F32
  {T.count_below(empty_atoms(), qp, c2) == 0n : Nat}

# An empty set has no clash, at any point and threshold.
law has_close_nil:
  for qp: G.Vec3
  for min2: F32
  {T.has_close(empty_atoms(), qp, min2) == False{} : Bool}

def empty_cx() -> P.Complex:
  P.Complex{Nil{}}

# Spread of nothing is zero.
law sum_dist2_nil:
  for qp: G.Vec3
  {T.sum_dist2(empty_atoms(), qp) == 0.0 : F32}

# An empty set links nothing to anything.
law contacts_between_nil:
  for ys: List<&2, P.Atom>
  for c2: F32
  {T.contacts_between(empty_atoms(), ys, c2) == 0n : Nat}

# No closest approach exists from an empty set.
law min_dist2_nil:
  for qp: G.Vec3
  {T.min_dist2_to(empty_atoms(), qp) == None{} : Maybe<&2, F32>}

# An empty complex holds no atoms.
law complex_atoms_nil:
  {P.Complex.num_atoms(empty_cx()) == 0n : Nat}

# LJ energy of nothing is zero, at any parameters.
law lj_nil:
  for eps: F32
  for sig2: F32
  for c2: F32
  {F.lj_total(empty_atoms(), empty_atoms(), eps, sig2, c2) == 0.0 : F32}

# The empty set is at zero RMSD from itself (empty case, by computation).
# (This says nothing about nonempty sets; see rmsd_single. Unequal lengths
# silently answer 0 by zip semantics; see rmsd_checked_some/none.)
law rmsd_nil:
  {R.rmsd(empty_atoms(), empty_atoms()) == 0.0 : F32}

# Nothing has no surface.
law sasa_nil:
  for rad: F32
  for rc2: F32
  for area_pt: F32
  {S.sasa_all(empty_atoms(), S.sasa_offsets(), rad, empty_atoms(), rc2, area_pt) == Nil{} : List<&2, F32>}

# Zero minimization steps change nothing, on any input, not just nothing.
law sd_fuel0:
  for xs: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  for dt: F32
  {F.minimize(0n, xs, eps, sig2, c2, dt) == xs : List<&2, P.Atom>}

# --- Step equations: each fold's cons case, by computation. ---
# Together with the Nil laws above, these pin every recursion's skeleton:
# each function is characterized on all inputs up to its head terms.

# Counting a cons cell is the head hit plus the count of the tail.
law count_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qp: G.Vec3
  for c2: F32
  {T.count_below(h <> t, qp, c2) == Nat.add(T.count_below(t, qp, c2), T.head_hit_nat(h, qp, c2)) : Nat}

# Interface counting peels one head the same way.
law contacts_cons:
  for x: P.Atom
  for xt: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for c2: F32
  {T.contacts_between(x <> xt, ys, c2) == Nat.add(T.contacts_between(xt, ys, c2), T.contacts_of_head(x, ys, c2)) : Nat}

# The fused spread/count pass extends by one head distance.
law sumd2_count_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qp: G.Vec3
  {T.sumd2_count(h <> t, qp) == T.sumd2_extend(T.head_dist2(h, qp), T.sumd2_count(t, qp)) : F32 & Nat}

# Coordination counts peel one head.
law coord_cons:
  for x: P.Atom
  for xt: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for c2: F32
  {T.coord_counts(x <> xt, ys, c2) == T.contacts_of_head(x, ys, c2) <> T.coord_counts(xt, ys, c2) : List<&2, Nat>}

# Closest-approach folds one head into the running best.
law min_dist2_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qp: G.Vec3
  {T.min_dist2_to(h <> t, qp) == T.min_cons(h, T.min_dist2_to(t, qp), qp) : Maybe<&2, F32>}

# Set linkage folds one head linkage into the running best.
law min_link2_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  {T.min_link2(h <> t, ys) == T.min_maybe(T.min_link2(t, ys), T.head_link(h, ys)) : Maybe<&2, F32>}

# An LJ row over a cons cell splits head from tail.
law lj_row_cons:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for h: P.Atom
  for t: List<&2, P.Atom>
  {F.lj_row(pi, eps, sig2, c2, h <> t) == (F.lj_row(pi, eps, sig2, c2, t) + F.lj_head(pi, eps, sig2, c2, h) : F32) : F32}

# An LJ total over a cons cell splits head from tail.
law lj_total_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  {F.lj_total(h <> t, ys, eps, sig2, c2) == (F.lj_total(t, ys, eps, sig2, c2) + F.lj_self_atom(h, ys, eps, sig2, c2) : F32) : F32}

# A Coulomb row splits head charge from tail charges.
law coul_row_cons:
  for pi: G.Vec3
  for qi: F32
  for ke: F32
  for c2: F32
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qh: F32
  for qt: List<&2, F32>
  {F.coul_row(pi, qi, ke, c2, h <> t, qh <> qt) == (F.coul_row(pi, qi, ke, c2, t, qt) + F.coul_head(pi, qi, ke, c2, h, qh) : F32) : F32}

# A Coulomb total splits head charge from tail charges.
law coul_total_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qh: F32
  for qt: List<&2, F32>
  for ys: List<&2, P.Atom>
  for qs_e: List<&2, F32>
  for ke: F32
  for c2: F32
  {F.coul_total(h <> t, qh <> qt, ys, qs_e, ke, c2) == (F.coul_total(t, qt, ys, qs_e, ke, c2) + F.coul_self(h, qh, ys, qs_e, ke, c2) : F32) : F32}

# A sweep over cons cells moves the head and sweeps the tail.
law sd_sweep_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for f: G.Vec3
  for ft: List<&2, G.Vec3>
  for dt: F32
  {F.sd_sweep(h <> t, f <> ft, dt) == F.sd_move(h, f, dt) <> F.sd_sweep(t, ft, dt) : List<&2, P.Atom>}

# A pair-distance fold splits the head pair from the tail pairs.
law pair_sumd2_cons:
  for x: P.Atom
  for xt: List<&2, P.Atom>
  for y: P.Atom
  for yt: List<&2, P.Atom>
  {R.pair_sumd2(x <> xt, y <> yt) == R.pair_head(x, y, R.pair_sumd2(xt, yt)) : F32 & Nat}

# Scoring against a cons ensemble scores the head and the tail.
law rmsd_to_all_cons:
  for q: List<&2, P.Atom>
  for h: List<&2, P.Atom>
  for t: List<&2, List<&2, P.Atom>>
  {R.rmsd_to_all(q, h <> t) == R.rmsd_against(q, h) <> R.rmsd_to_all(q, t) : List<&2, F32>}

# --- Singletons and golden values: the heads compute. ---

# Counting one atom is that atom's hit.
law count_single:
  for h: P.Atom
  for qp: G.Vec3
  for c2: F32
  {T.count_below([h], qp, c2) == T.head_hit_nat(h, qp, c2) : Nat}

# One atom is close exactly when its head test says so.
law has_close_single:
  for h: P.Atom
  for qp: G.Vec3
  for min2: F32
  {T.has_close([h], qp, min2) == T.head_is_close(h, qp, min2) : Bool}

# The closest approach to a singleton is its one distance.
law min_single:
  for s: U32
  for e: U32
  for x: F32
  for y: F32
  for z: F32
  for qx: F32
  for qy: F32
  for qz: F32
  {T.min_dist2_to([P.Atom{s, e, G.V3{x, y, z}}], G.V3{qx, qy, qz}) == Some{G.Vec3.dist2(G.V3{x, y, z}, G.V3{qx, qy, qz})} : Maybe<&2, F32>}

# Spread about a point of nothing is zero through rg_of.
law rg_about_nil:
  for qp: G.Vec3
  {T.rg_about(empty_atoms(), qp) == 0.0 : F32}

# A singleton is at the RMSD its one pair distance dictates.
law rmsd_single:
  for s: U32
  for e: U32
  for x: F32
  for y: F32
  for z: F32
  {R.rmsd([P.Atom{s, e, G.V3{x, y, z}}], [P.Atom{s, e, G.V3{x, y, z}}]) == F32.sqrt(((G.Vec3.dist2(G.V3{x, y, z}, G.V3{x, y, z}) + 0.0 : F32) / F32.from_nat(1n) : F32)) : F32}

# splitmix32 golden values (catches mistyped mixer constants).
law rng_next_0:
  {Rng.rng_next(0) == 1684164658 : U32}

law rng_next_1:
  {Rng.rng_next(1) == 1580013426 : U32}

# One atom with nothing near it exposes all six sample points.
law sasa_single:
  {S.sasa_of([P.Atom{1, 6, G.V3{0.0, 0.0, 0.0}}], 2.0, 16.0, 1.0) == [(1.0 * F32.from_nat(6n) : F32)] : List<&2, F32>}

# --- Truncation pins: mismatched lockstep inputs drop the rest. ---

# Pair folds against an empty list answer the identity.
law pair_trunc_right:
  for x: P.Atom
  for xt: List<&2, P.Atom>
  {R.pair_sumd2(x <> xt, empty_atoms()) == (0.0, 0n) : F32 & Nat}

# Coulomb with no charges answers zero, however many atoms remain.
law coul_trunc_charges:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for qs: List<&2, F32>
  for ke: F32
  for c2: F32
  {F.coul_total(h <> t, Nil{}, ys, qs, ke, c2) == 0.0 : F32}

# A sweep with no forces leaves the atoms in place.
law sd_trunc_forces:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for dt: F32
  {F.sd_sweep(h <> t, Nil{}, dt) == h <> t : List<&2, P.Atom>}

# --- Induction: empty arguments annihilate, by induction on atoms. ---

# Nothing links to anything, however many heads are peeled.
law contacts_nil_right:
  for xs: List<&2, P.Atom>
  for c2: F32
  {T.contacts_between(xs, empty_atoms(), c2) == 0n : Nat}

# No closest approach exists from any set held against nothing.
law min_link2_nil_right:
  for xs: List<&2, P.Atom>
  {T.min_link2(xs, empty_atoms()) == None{} : Maybe<&2, F32>}

# --- Checked entry points: mismatches are values, not silence. ---

# Equal-length singletons check cleanly.
law rmsd_checked_some:
  for x: P.Atom
  for y: P.Atom
  {R.rmsd_checked([x], [y]) == Some{R.rmsd([x], [y])} : Maybe<&2, F32>}

# Unequal lengths answer None{} instead of a silent zero.
law rmsd_checked_none:
  for x: P.Atom
  {R.rmsd_checked([x], empty_atoms()) == None{} : Maybe<&2, F32>}

# --- Structural algebra: the hierarchy agrees with the flat lists. ---

# Appending atom lists associates.
law append_atoms_assoc:
  for xs: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for zs: List<&2, P.Atom>
  {P.append_atoms(P.append_atoms(xs, ys), zs) == P.append_atoms(xs, P.append_atoms(ys, zs)) : List<&2, P.Atom>}

# Flattening distributes over appending residue lists.
law flatten_residues_append:
  for xs: List<&2, P.Residue>
  for ys: List<&2, P.Residue>
  {P.flatten_residues(P.append_residues(xs, ys)) == P.append_atoms(P.flatten_residues(xs), P.flatten_residues(ys)) : List<&2, P.Atom>}

# Adding zero on the right changes nothing (induction seed for comm).
law add_zero_right:
  for b: Nat
  {Nat.add(b, 0n) == b : Nat}

# Adding a successor on the right pulls out front.
law add_succ_right:
  for b: Nat
  for a: Nat
  {Nat.add(b, 1n+a) == 1n+Nat.add(b, a) : Nat}

# Nat addition commutes.
law add_comm:
  for a: Nat
  for b: Nat
  {Nat.add(a, b) == Nat.add(b, a) : Nat}

# Counting an append splits into the counts.
law count_atoms_append:
  for xs: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  {P.count_atoms(P.append_atoms(xs, ys)) == Nat.add(P.count_atoms(xs), P.count_atoms(ys)) : Nat}

# Counting a flattening agrees with counting residues in place.
law count_flatten_eq:
  for rs: List<&2, P.Residue>
  {P.count_atoms(P.flatten_residues(rs)) == P.count_atoms_in(rs) : Nat}

# A complex holds as many atoms as its flattening lists.
law complex_flatten_eq:
  for chains: List<&2, P.Chain>
  {P.Complex.num_atoms(P.Complex{chains}) == P.count_atoms(P.Complex.flatten(P.Complex{chains})) : Nat}

# --- PDB reader: slices, filters, parser goldens, end-to-end counts. ---

# 10^3 is 1000.
law pow10_3:
  {Pdb.pow10(3n) == 1000 : U32}

# Iron reads as 26.
law elem_iron:
  {Pdb.elem_no("FE") == 26 : U32}

# Carbon reads as 6.
law elem_carbon:
  {Pdb.elem_1('C') == 6 : U32}

# ATOM records (6-column tag) pass the filter.
law atom_line_yes:
  {Pdb.is_atom_line("ATOM      1  N   ALA A   1") == True{} : Bool}

# Anything else does not.
law atom_line_no:
  {Pdb.is_atom_line("HETATM 9999  O   HOH A 101") == False{} : Bool}

# No lines in, no lines out.
law atom_lines_nil:
  {Pdb.atom_lines(Nil{}) == Nil{} : List<&2, String>}

# Non-ATOM lines drop, ENDMDL stops the file, the tail never parses.
law atom_lines_keep:
  {Pdb.atom_lines(["ATOM      1  N", "TER", "ENDMDL", "ATOM      2  CA"]) == ["ATOM      1  N"] : List<&2, String>}

# 12.5 parses to 12 + 5/10 (stuck F32, identical both sides).
law parse_f32_one:
  {Pdb.parse_f32("12.5") == Some{(U32.to_f32(12) + (U32.to_f32(5) / U32.to_f32(10) : F32) : F32)} : Maybe<&2, F32>}

# Signs apply outside the magnitude.
law parse_f32_neg:
  {Pdb.parse_f32("-3.25") == Some{(0.0 - (U32.to_f32(3) + (U32.to_f32(25) / U32.to_f32(100) : F32) : F32) : F32)} : Maybe<&2, F32>}

# Non-numeric fields fail instead of answering zero.
law parse_f32_bad:
  {Pdb.parse_f32("abc") == None{} : Maybe<&2, F32>}

# The two-line ALA snippet parses to two atoms ...
law pdb_num_atoms:
  {P.Complex.num_atoms(Pdb.complex_of(Pdb.parse_pdb(Pdb.demo_text()))) == 2n : Nat}

# ... in one residue ...
law pdb_num_residues:
  {P.Complex.num_residues(Pdb.complex_of(Pdb.parse_pdb(Pdb.demo_text()))) == 1n : Nat}

# ... in one chain ...
law pdb_num_chains:
  {P.Complex.num_chains(Pdb.complex_of(Pdb.parse_pdb(Pdb.demo_text()))) == 1n : Nat}

# ... with nothing skipped ...
law pdb_skipped_0:
  {Pdb.skipped_of(Pdb.parse_pdb(Pdb.demo_text())) == 0n : Nat}

# ... while a bad serial skips exactly its line.
law pdb_skipped_1:
  {Pdb.skipped_of(Pdb.parse_pdb(Pdb.demo_bad_text())) == 1n : Nat}

law pdb_bad_atoms:
  {P.Complex.num_atoms(Pdb.complex_of(Pdb.parse_pdb(Pdb.demo_bad_text()))) == 1n : Nat}

# ... and the first atom is nitrogen.
law pdb_first_elem:
  {Pdb.first_elem(Pdb.complex_of(Pdb.parse_pdb(Pdb.demo_text()))) == Some{7} : Maybe<&2, U32>}

# Total area splits head from tail.
law sasa_total_cons:
  for h: F32
  for t: List<&2, F32>
  {S.sasa_total(h <> t) == (h + S.sasa_total(t) : F32) : F32}

# --- Parallel kernels: unrolled steps match, remainders fall back. ---

# Empty input counts nothing, through the fallback arm.
law count_par4_nil:
  for qp: G.Vec3
  for c2: F32
  {Par.count_par4(empty_atoms(), qp, c2) == 0n : Nat}

# A full unrolled step splits four parallel hits plus the tail.
law count_par4_step:
  for h1: P.Atom
  for h2: P.Atom
  for h3: P.Atom
  for h4: P.Atom
  for t: List<&2, P.Atom>
  for qp: G.Vec3
  for c2: F32
  {Par.count_par4(h1 <> h2 <> h3 <> h4 <> t, qp, c2) == Nat.add(Nat.add(T.head_hit_nat(h1, qp, c2), T.head_hit_nat(h2, qp, c2)), Nat.add(T.head_hit_nat(h3, qp, c2), Nat.add(T.head_hit_nat(h4, qp, c2), Par.count_par4(t, qp, c2)))) : Nat}

# Short inputs take the sequential fallback exactly.
law count_par4_small:
  for h: P.Atom
  for qp: G.Vec3
  for c2: F32
  {Par.count_par4([h], qp, c2) == T.count_below([h], qp, c2) : Nat}

# Empty rows sum nothing, through the fallback arm.
law lj_row_par4_nil:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  {Par.lj_row_par4(pi, eps, sig2, c2, empty_atoms()) == 0.0 : F32}

# A full unrolled row splits four parallel heads plus the tail.
law lj_row_par4_step:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for h1: P.Atom
  for h2: P.Atom
  for h3: P.Atom
  for h4: P.Atom
  for t: List<&2, P.Atom>
  {Par.lj_row_par4(pi, eps, sig2, c2, h1 <> h2 <> h3 <> h4 <> t) == ((F.lj_head(pi, eps, sig2, c2, h1) + F.lj_head(pi, eps, sig2, c2, h2) : F32) + ((F.lj_head(pi, eps, sig2, c2, h3) + F.lj_head(pi, eps, sig2, c2, h4) : F32) + Par.lj_row_par4(pi, eps, sig2, c2, t) : F32) : F32) : F32}

# Short rows take the sequential fallback exactly.
law lj_row_par4_small:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for h: P.Atom
  {Par.lj_row_par4(pi, eps, sig2, c2, [h]) == F.lj_row(pi, eps, sig2, c2, [h]) : F32}

# --- Per-element LJ: table goldens and step equations. ---

# Carbon takes illustrative OPLS-style parameters.
law lj_eps_c:
  {F.lj_eps(6) == 0.066 : F32}

law lj_eps_h:
  {F.lj_eps(1) == 0.03 : F32}

law lj_sig_o:
  {F.lj_sig(8) == 2.96 : F32}

# Unknown elements take the generic pair.
law lj_sig_default:
  {F.lj_sig(99) == 3.0 : F32}

# An element row over a cons cell splits head from tail.
law lj_row_elem_cons:
  for si: U32
  for pi: G.Vec3
  for c2: F32
  for h: P.Atom
  for t: List<&2, P.Atom>
  {F.lj_row_elem(si, pi, h <> t, c2) == (F.lj_row_elem(si, pi, t, c2) + F.lj_head_elem(si, pi, c2, h) : F32) : F32}

# An element total over a cons cell splits head from tail.
law lj_total_elem_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for c2: F32
  {F.lj_total_elem(h <> t, ys, c2) == (F.lj_total_elem(t, ys, c2) + F.lj_row_elem(P.Atom.elem(h), P.Atom.pos(h), ys, c2) : F32) : F32}

# --- Bonded terms: empties, steps, lookups, exclusion goldens. ---

# No bonds bind nothing.
law bond_total_nil:
  for k: F32
  {B.bond_total(empty_atoms(), Nil{}, k) == 0.0 : F32}

# A bond total splits head bond from tail bonds.
law bond_total_cons:
  for xs: List<&2, P.Atom>
  for s1: U32
  for s2: U32
  for t: List<&1, U32 & U32>
  for k: F32
  {B.bond_total(xs, (s1, s2) <> t, k) == (B.bond_cov_pair(xs, s1, s2, k) + B.bond_total(xs, t, k) : F32) : F32}

# No angles bend nothing.
law angle_total_nil:
  for k: F32
  for eq: F32
  {B.angle_total(empty_atoms(), Nil{}, k, eq) == 0.0 : F32}

# An angle total splits head triple from tail triples.
law angle_total_cons:
  for xs: List<&2, P.Atom>
  for sa: U32
  for sb: U32
  for sc: U32
  for t: List<&1, U32 & U32 & U32>
  for k: F32
  for eq: F32
  {B.angle_total(xs, (sa, sb, sc) <> t, k, eq) == (B.angle_triple(xs, sa, sb, sc, k, eq) + B.angle_total(xs, t, k, eq) : F32) : F32}

# Nobody is found in nothing.
law pos_of_serial_none:
  for s: U32
  {B.pos_of_serial(empty_atoms(), s) == None{} : Maybe<&2, G.Vec3>}

# A singleton finds its atom's position.
law pos_of_serial_some:
  {B.pos_of_serial([P.Atom{1, 6, G.V3{1.0, 2.0, 3.0}}], 1) == Some{G.V3{1.0, 2.0, 3.0}} : Maybe<&2, G.Vec3>}

# No bonds exclude nothing.
law bonded12_nil:
  for a: U32
  for b: U32
  {B.bonded12_flat(Nil{}, a, b) == False{} : Bool}

# A listed pair (either order) is 1-2 excluded.
law bonded12_hit:
  {B.bonded12_flat([1, 2], 1, 2) == True{} : Bool}

# An unlisted pair is not.
law bonded12_miss:
  {B.bonded12_flat([1, 2], 1, 3) == False{} : Bool}

# Neighbors of 1 in [(1,2),(2,3)] are just [2].
law neighbors_one:
  {B.bond_neighbors(1, [1, 2, 2, 3]) == [2] : List<&2, U32>}

# 1 and 3 share neighbor 2.
law excluded13_one:
  {B.excluded13(1, 3, [1, 2, 2, 3]) == True{} : Bool}

# 2 and 3 are directly bonded (1-2 wins over 1-3 either way).
law excluded_both:
  {B.excluded(2, 3, [1, 2, 2, 3]) == True{} : Bool}

# Pair bonds flatten in order.
law flatten_bonds_one:
  {B.flatten_bonds([(1, 2)]) == [1, 2] : List<&2, U32>}

# Excluded rows sum nothing over nothing.
law lj_row_excl_nil:
  for si: U32
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for bonds: List<&2, U32>
  {B.lj_row_excl(si, pi, eps, sig2, c2, empty_atoms(), bonds) == 0.0 : F32}

law lj_total_excl_nil:
  for ys: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  for bonds: List<&2, U32>
  {B.lj_total_excl(empty_atoms(), ys, eps, sig2, c2, bonds) == 0.0 : F32}

law coul_row_excl_nil:
  for si: U32
  for pi: G.Vec3
  for qi: F32
  for ke: F32
  for c2: F32
  for bonds: List<&2, U32>
  {B.coul_row_excl(si, pi, qi, ke, c2, empty_atoms(), Nil{}, bonds) == 0.0 : F32}

law coul_total_excl_nil:
  for ys: List<&2, P.Atom>
  for qs_e: List<&2, F32>
  for ke: F32
  for c2: F32
  for bonds: List<&2, U32>
  {B.coul_total_excl(empty_atoms(), Nil{}, ys, qs_e, ke, c2, bonds) == 0.0 : F32}

# Element-excluded rows/totals sum nothing over nothing.
law lj_row_elem_excl_nil:
  for si: U32
  for se: U32
  for pi: G.Vec3
  for c2: F32
  for bonds: List<&2, U32>
  {B.lj_row_elem_excl(si, se, pi, empty_atoms(), c2, bonds) == 0.0 : F32}

law lj_total_elem_excl_nil:
  for ys: List<&2, P.Atom>
  for c2: F32
  for bonds: List<&2, U32>
  {B.lj_total_elem_excl(empty_atoms(), ys, c2, bonds) == 0.0 : F32}

# A directly bonded pair contributes exactly 0.0 (wiring pin: the row must
# test serials, not elements; the sum is (0.0 + 0.0) since F32.add is stuck).
law lj_row_elem_excl_hit:
  {B.lj_row_elem_excl(1, 7, G.V3{0.0, 0.0, 0.0}, [P.Atom{2, 6, G.V3{1.0, 0.0, 0.0}}], 144.0, [1, 2]) == (0.0 + 0.0 : F32) : F32}

# Same through a 1-3 link.
law lj_row_elem_excl_link:
  {B.lj_row_elem_excl(1, 7, G.V3{0.0, 0.0, 0.0}, [P.Atom{3, 6, G.V3{2.0, 0.0, 0.0}}], 144.0, [1, 2, 2, 3]) == (0.0 + 0.0 : F32) : F32}

# --- CONECT: filter and parse goldens. ---

# No lines in, no bonds out.
law conect_lines_nil:
  {Pdb.conect_lines(Nil{}) == Nil{} : List<&2, String>}

# CONECT lines survive, everything else drops.
law conect_keep:
  {Pdb.conect_lines(["ATOM      1  N", "CONECT   20  282"]) == ["CONECT   20  282"] : List<&2, String>}

# One record with one partner is one bond.
law parse_conect_one:
  {Pdb.parse_conect_line("CONECT   20  282") == [(20, 282)] : List<&1, U32 & U32>}

# A partnerless record binds nothing.
law parse_conect_none:
  {Pdb.parse_conect_line("CONECT") == Nil{} : List<&1, U32 & U32>}

# --- Masses and Coulomb forces. ---

# Carbon weighs 12.011 amu.
law mass_c:
  {F.mass(6) == 12.011 : F32}

law mass_h:
  {F.mass(1) == 1.008 : F32}

# A Coulomb row over a cons cell splits head from tail.
law coul_frow_cons:
  for pi: G.Vec3
  for qi: F32
  for ke: F32
  for c2: F32
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qh: F32
  for qt: List<&2, F32>
  {F.coul_frow(pi, qi, ke, c2, h <> t, qh <> qt) == G.Vec3.add(F.coul_frow(pi, qi, ke, c2, t, qt), F.coul_fhead(pi, qi, ke, c2, h, qh)) : G.Vec3}

# Coulomb forces over cons cells split head from tail.
law coul_forces_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for qh: F32
  for qt: List<&2, F32>
  for ys: List<&2, P.Atom>
  for qs_e: List<&2, F32>
  for ke: F32
  for c2: F32
  {F.coul_forces(h <> t, qh <> qt, ys, qs_e, ke, c2) == F.coul_fself_atom(h, qh, ys, qs_e, ke, c2) <> F.coul_forces(t, qt, ys, qs_e, ke, c2) : List<&2, G.Vec3>}

# No charges feel nothing.
law coul_frow_nil:
  for pi: G.Vec3
  for qi: F32
  for ke: F32
  for c2: F32
  {F.coul_frow(pi, qi, ke, c2, empty_atoms(), Nil{}) == G.V3{0.0, 0.0, 0.0} : G.Vec3}

# --- PBC: mirror equations and empties. ---

law lj_row_pbc_cons:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  for h: P.Atom
  for t: List<&2, P.Atom>
  {Pbc.lj_row_pbc(pi, eps, sig2, c2, box, h <> t) == (Pbc.lj_row_pbc(pi, eps, sig2, c2, box, t) + Pbc.lj_head_pbc(pi, eps, sig2, c2, box, h) : F32) : F32}

law lj_total_pbc_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  {Pbc.lj_total_pbc(h <> t, ys, eps, sig2, c2, box) == (Pbc.lj_total_pbc(t, ys, eps, sig2, c2, box) + Pbc.lj_self_atom_pbc(h, ys, eps, sig2, c2, box) : F32) : F32}

law lj_row_pbc_nil:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  {Pbc.lj_row_pbc(pi, eps, sig2, c2, box, empty_atoms()) == 0.0 : F32}

law lj_frow_pbc_cons:
  for pi: G.Vec3
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  for h: P.Atom
  for t: List<&2, P.Atom>
  {Pbc.lj_frow_pbc(pi, eps, sig2, c2, box, h <> t) == G.Vec3.add(Pbc.lj_frow_pbc(pi, eps, sig2, c2, box, t), Pbc.lj_fhead_pbc(pi, eps, sig2, c2, box, h)) : G.Vec3}

law lj_forces_pbc_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for ys: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  {Pbc.lj_forces_pbc(h <> t, ys, eps, sig2, c2, box) == Pbc.lj_fself_atom_pbc(h, ys, eps, sig2, c2, box) <> Pbc.lj_forces_pbc(t, ys, eps, sig2, c2, box) : List<&2, G.Vec3>}

law lj_forces_pbc_nil:
  for ys: List<&2, P.Atom>
  for eps: F32
  for sig2: F32
  for c2: F32
  for box: G.Vec3
  {Pbc.lj_forces_pbc(empty_atoms(), ys, eps, sig2, c2, box) == Nil{} : List<&2, G.Vec3>}

# --- Bonded additions: radii, inference, forces. ---

# Carbon's covalent radius is 0.76 A.
law cov_rad_c:
  {B.cov_rad(6) == 0.76 : F32}

# An inference row over nothing finds nothing.
law infer_row_nil:
  for si: U32
  for pi: G.Vec3
  for ei: U32
  {B.infer_row(si, pi, ei, empty_atoms()) == Nil{} : List<&1, U32 & U32>}

# An inference row splits head from tail.
law infer_row_cons:
  for si: U32
  for pi: G.Vec3
  for ei: U32
  for h: P.Atom
  for t: List<&2, P.Atom>
  {B.infer_row(si, pi, ei, h <> t) == B.append_infer(B.infer_one(si, pi, ei, P.Atom.serial(h), P.Atom.elem(h), P.Atom.pos(h)), B.infer_row(si, pi, ei, t)) : List<&1, U32 & U32>}

# Nothing infers from nothing.
law infer_total_nil:
  {B.infer_total(empty_atoms()) == Nil{} : List<&1, U32 & U32>}

# Bond forces sum nothing over nothing.
law bond_forces_nil:
  for bonds: List<&2, U32>
  for k: F32
  {B.bond_forces(empty_atoms(), bonds, k) == Nil{} : List<&2, G.Vec3>}

# Bond forces split head atom from tail atoms.
law bond_forces_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  for bonds: List<&2, U32>
  for k: F32
  {B.bond_forces(h <> t, bonds, k) == B.bond_force_list(P.Atom.serial(h), bonds, h <> t, k) <> B.bond_forces(t, bonds, k) : List<&2, G.Vec3>}

# --- Engine: maps split head from tail; drivers stop on empty fuel. ---

law dyn_of_atoms_nil:
  {Sim.dyn_of_atoms(empty_atoms()) == Nil{} : List<&2, Sim.DynAtom>}

law dyn_of_atoms_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  {Sim.dyn_of_atoms(h <> t) == Sim.dyn_of_atom(h) <> Sim.dyn_of_atoms(t) : List<&2, Sim.DynAtom>}

law dyn_atoms_nil:
  {Sim.dyn_atoms(Nil{}) == Nil{} : List<&2, P.Atom>}

law vv_scale_nil:
  for lam: F32
  {Sim.vv_scale(Nil{}, lam) == Nil{} : List<&2, Sim.DynAtom>}

law vv_scale_cons:
  for h: Sim.DynAtom
  for t: List<&2, Sim.DynAtom>
  for lam: F32
  {Sim.vv_scale(h <> t, lam) == Sim.vv_scale_head(h, lam) <> Sim.vv_scale(t, lam) : List<&2, Sim.DynAtom>}

law vv_half_nil:
  for dt: F32
  {Sim.vv_half(Nil{}, dt) == Nil{} : List<&2, Sim.DynAtom>}

law vv_half_cons:
  for h: Sim.DynAtom
  for t: List<&2, Sim.DynAtom>
  for dt: F32
  {Sim.vv_half(h <> t, dt) == Sim.vv_half_head(h, dt) <> Sim.vv_half(t, dt) : List<&2, Sim.DynAtom>}

law vv_full_nil:
  for dt: F32
  {Sim.vv_full(Nil{}, dt) == Nil{} : List<&2, Sim.DynAtom>}

law vv_full_cons:
  for h: Sim.DynAtom
  for t: List<&2, Sim.DynAtom>
  for dt: F32
  {Sim.vv_full(h <> t, dt) == Sim.vv_full_head(h, dt) <> Sim.vv_full(t, dt) : List<&2, Sim.DynAtom>}

law ke_sim_nil:
  {Sim.ke_sim_of(Nil{}) == 0.0 : F32}

law md_run_zero:
  for ds: List<&2, Sim.DynAtom>
  for qs: List<&2, F32>
  for bonds: List<&2, U32>
  for p: Sim.SimParams
  for dt: F32
  for tau: F32
  for temp0: F32
  for n: Nat
  {Sim.md_run(0n, ds, qs, bonds, p, dt, tau, temp0, n) == ds : List<&2, Sim.DynAtom>}

law min_run_zero:
  for ds: List<&2, Sim.DynAtom>
  for qs: List<&2, F32>
  for bonds: List<&2, U32>
  for p: Sim.SimParams
  for dt: F32
  {Sim.min_run(0n, ds, qs, bonds, p, dt) == ds : List<&2, Sim.DynAtom>}

# Zero charges over nothing is nothing.
law zero_charges_nil:
  {Sim.zero_charges(empty_atoms()) == Nil{} : List<&2, F32>}

# Zero charges split head from tail.
law zero_charges_cons:
  for h: P.Atom
  for t: List<&2, P.Atom>
  {Sim.zero_charges(h <> t) == 0.0 <> Sim.zero_charges(t) : List<&2, F32>}

# Normalized sweeps move nothing over nothing.
law sd_norm_sweep_nil:
  for dt: F32
  {Sim.sd_norm_sweep(Nil{}, dt) == Nil{} : List<&2, Sim.DynAtom>}

# No atoms carry no max force.
law maxforce2_nil:
  {Sim.maxforce2_of(Nil{}) == 0.0 : F32}

# Setting forces over nothing is nothing.
law dyn_set_force_nil:
  for fs: List<&2, G.Vec3>
  {Sim.dyn_set_force(Nil{}, fs) == Nil{} : List<&2, Sim.DynAtom>}

# Adding empties is empty.
law add_forces_nil:
  {Sim.add_forces(Nil{}, Nil{}) == Nil{} : List<&2, G.Vec3>}

# Force addition splits head pair from tail pairs.
law add_forces_cons:
  for x: G.Vec3
  for xt: List<&2, G.Vec3>
  for y: G.Vec3
  for yt: List<&2, G.Vec3>
  {Sim.add_forces(x <> xt, y <> yt) == G.Vec3.add(x, y) <> Sim.add_forces(xt, yt) : List<&2, G.Vec3>}

# No flats pair to nothing.
law pairs_of_flat_nil:
  {Sim.pairs_of_flat(Nil{}) == Nil{} : List<&1, U32 & U32>}

# Flat pairs split head pair from tail pairs.
law pairs_of_flat_cons:
  for x: U32
  for y: U32
  for t: List<&2, U32>
  {Sim.pairs_of_flat(x <> y <> t) == (x, y) <> Sim.pairs_of_flat(t) : List<&1, U32 & U32>}

# Padding and symbols print as written.
law pad_left_golden:
  {Pdb.pad_left("AB", 5n) == "   AB" : String}

law pad_right_golden:
  {Pdb.pad_right("AB", 5n) == "AB   " : String}

law fmt_u32_golden:
  {Pdb.fmt_u32(5n, 42) == "   42" : String}

law elem_sym_golden:
  {Pdb.elem_sym(26) == "FE" : String}

# --- Templates: dispatch goldens, name matching, parse goldens. ---

# ALA links N-CA in the template.
law tpl_has_bond_ala:
  {Tpl.tpl_has_bond("ALA", "N", "CA") == True{} : Bool}

# ... but not N-O.
law tpl_has_bond_miss:
  {Tpl.tpl_has_bond("ALA", "N", "O") == False{} : Bool}

# Unknown residues link nothing.
law tpl_has_bond_unknown:
  {Tpl.tpl_has_bond("ZZZ", "N", "CA") == False{} : Bool}

# HIS variants share heavy topology.
law tpl_has_bond_his:
  {Tpl.tpl_has_bond("HID", "N", "CA") == True{} : Bool}

# Template rows split head from tail.
law tpl_row_nil:
  for a: Tpl.TplAtom
  {Tpl.tpl_row(a, Nil{}) == Nil{} : List<&1, U32 & U32>}

law tpl_row_cons:
  for a: Tpl.TplAtom
  for h: Tpl.TplAtom
  for t: List<&2, Tpl.TplAtom>
  {Tpl.tpl_row(a, h <> t) == Tpl.append_tpl(Tpl.tpl_one(a, h), Tpl.tpl_row(a, t)) : List<&1, U32 & U32>}

law tpl_total_nil:
  {Tpl.tpl_total(Nil{}) == Nil{} : List<&1, U32 & U32>}

# The snippet's first atom parses with identity intact.
law parse_tpl_golden:
  {Tpl.parse_tpl_line("ATOM      1  N   ALA A   1      11.104  13.207   2.100  1.00 13.79           N  ") == Some{Tpl.TA{1, 65, 1n, "ALA", "N"}} : Maybe<&2, Tpl.TplAtom>}

# --- Angle auto-derivation goldens. ---

# Neighbors of 2 excluding 1 in [(1,2),(2,3)] are just [3].
law neighbors_except_golden:
  {B.bond_neighbors_except(2, 1, [1, 2, 2, 3]) == [3] : List<&2, U32>}

# A fan over [3, 4] from (1, 2) lists both triples.
law fan_triples_golden:
  {B.fan_triples(1, 2, [3, 4]) == [(1, 2, 3), (1, 2, 4)] : List<&1, U32 & U32 & U32>}

# A two-bond chain auto-derives both directed angles.
law angle_triples_auto_golden:
  {B.angle_triples_auto([1, 2, 2, 3]) == [(1, 2, 3), (3, 2, 1)] : List<&1, U32 & U32 & U32>}