/* ──────────────────────────────────────────────────────────────────────────
 * Zeq · Nano-Zeqond Molecular-Dynamics Engine  —  Phase 2 core
 *
 * A precise 2-D molecular-dynamics integrator driven by the same nano-zeqond
 * clock as the field evolver. The HULYAS master equation's NM-family operators
 * (NM19 F=ma, NM30 harmonic bonds, the Lennard-Jones pair force) ARE the field
 * terms here — particles instead of a continuum ϕ. The contract:
 *
 *   • Integrator: velocity-Verlet (kick–drift–kick) — 2nd-order symplectic,
 *     the standard for MD; energy error O(dt²) and BOUNDED → KO42 ≤0.1% drift.
 *   • Forces: shifted-truncated Lennard-Jones (12-6) with minimum-image
 *     periodic boundaries + optional harmonic bonds  ½k(r−r₀)².
 *   • Reduced LJ units (m = ε = σ = k_B = 1). Map to a real species (e.g.
 *     argon: σ=3.405 Å, ε/k_B=119.8 K, τ=2.156 ps) at the UI layer.
 *   • Clock: the metric (and dt) pulses with the 1.287 Hz HulyaPulse
 *     (KO42: dt = dtSim·R(t)); a snapshot PRINTS every 0.777 s zeqond.
 *   • Conserves total energy (drift reported) and total momentum (COM fixed).
 *
 * HONEST throughput: atomic MD resolves femtosecond steps, so this delivers
 * nanoseconds–microseconds of molecular time per wall-second — real, useful MD
 * throughput, not "billions of years" (that framing is for coarse continua).
 * ────────────────────────────────────────────────────────────────────────── */
(function (root, factory) {
  if (typeof module === 'object' && module.exports) module.exports = factory();
  else root.ZeqMdEngine = factory();
}(typeof self !== 'undefined' ? self : this, function () {
  'use strict';

  const F_PULSE = 1.287;     // HulyaPulse frequency (Hz)
  const TAU_Z   = 0.777;     // zeqond (s)
  const ALPHA   = 1.29e-3;   // Zeq proper-time modulation amplitude

  // ── preset systems (real, distinct physics regimes) ──
  const PRESETS = {
    'lj-liquid':  { label: 'Lennard-Jones liquid',  N: 120, rho: 0.78, T0: 1.0,  lattice: false, bonds: false },
    'lj-gas':     { label: 'Lennard-Jones gas',     N: 80,  rho: 0.15, T0: 1.8,  lattice: false, bonds: false },
    'crystal':    { label: 'Crystal → melt',        N: 144, rho: 1.05, T0: 0.35, lattice: true,  bonds: false },
    'polymer':    { label: 'Bonded polymer chains', N: 100, rho: 0.45, T0: 0.9,  lattice: false, bonds: true  },
  };

  class ZeqMdEngine {
    constructor(opts = {}) {
      const presetKey = opts.preset && PRESETS[opts.preset] ? opts.preset : 'lj-liquid';
      const P = PRESETS[presetKey];
      this.presetKey = presetKey;
      this.preset = P;

      this.N    = opts.N    || P.N;
      this.rho  = opts.rho  != null ? opts.rho : P.rho;   // number density (N/area)
      this.T0   = opts.T0   != null ? opts.T0  : P.T0;    // initial temperature (reduced)
      this.eps  = opts.eps  != null ? opts.eps : 1.0;     // LJ ε
      this.sig  = opts.sig  != null ? opts.sig : 1.0;     // LJ σ
      this.rc   = opts.rc   != null ? opts.rc  : 2.5;     // cutoff (·σ)
      this.kBond = opts.kBond != null ? opts.kBond : 30.0; // harmonic bond stiffness
      this.r0Bond = opts.r0Bond != null ? opts.r0Bond : 1.0; // bond rest length

      // square periodic box from density
      this.L = Math.sqrt(this.N / this.rho);

      // state arrays (2 components per particle)
      this.x = new Float64Array(2 * this.N);
      this.v = new Float64Array(2 * this.N);
      this.f = new Float64Array(2 * this.N);
      this.bonds = [];            // [{i,j}]

      // ── clock + resolution (shared model with the field evolver) ──
      this.f_pulse = F_PULSE; this.tauZ = TAU_Z; this.alpha = ALPHA;
      this.phi0 = opts.phi0 != null ? opts.phi0 : 0;
      this.subdiv = opts.subdiv || 1e6;
      this.dtSim = opts.dtSim != null ? opts.dtSim : 0.003;   // MD step (reduced time) — VV-stable
      // simulated time advanced per step. Default = dtSim (1 reduced-time unit
      // per unit), so simTime is honest reduced-time; UI maps to ps/ns.
      this.simDtPerStep = opts.simDtPerStep != null ? opts.simDtPerStep : this.dtSim;

      this.wallT = 0; this.simTime = 0; this.steps = 0;
      this.E0 = null; this.lastZeqond = 0; this.prints = [];

      // precompute LJ energy shift at cutoff (continuity)
      const ir2 = 1 / (this.rc * this.rc), ir6 = ir2 * ir2 * ir2;
      this.ecut = 4 * this.eps * (ir6 * ir6 - ir6);

      this.init();
    }

    /** Initialize positions (lattice or random non-overlapping) + velocities. */
    init() {
      const N = this.N, L = this.L, x = this.x, v = this.v;
      if (this.preset.lattice) {
        // triangular lattice (2-D close packing)
        const cols = Math.ceil(Math.sqrt(N));
        const a = L / cols;
        let k = 0;
        for (let r = 0; r < cols && k < N; r++) {
          for (let c = 0; c < cols && k < N; c++) {
            x[2 * k]     = ((c + 0.5 * (r & 1)) * a) % L;
            x[2 * k + 1] = (r * a * 0.866) % L;
            k++;
          }
        }
      } else if (this.preset.bonds) {
        // build connected chains of 5 beads spaced near the bond rest length,
        // so bonds start near r₀ (no catastrophic initial stretch)
        const CHAIN = 5, step = this.r0Bond;
        let k = 0;
        while (k < N) {
          let px = Math.random() * L, py = Math.random() * L;
          for (let c = 0; c < CHAIN && k < N; c++) {
            x[2 * k] = ((px % L) + L) % L; x[2 * k + 1] = ((py % L) + L) % L;
            if (c < CHAIN - 1) { const a = Math.random() * 2 * Math.PI; px += step * Math.cos(a); py += step * Math.sin(a); }
            if (c < CHAIN - 1 && k + 1 < N) this.bonds.push({ i: k, j: k + 1 });
            k++;
          }
        }
      } else {
        // random placement; overlaps removed by minimization below
        for (let i = 0; i < 2 * N; i++) x[i] = Math.random() * L;
        for (let pass = 0; pass < 20; pass++) this.relaxOverlaps();
      }
      // energy-minimize to a clean, low-force start (kills overlaps, settles
      // bonds & lattice) → NVE run conserves energy to ≤0.1% with no thermal
      // transient injected by violent initial contacts.
      this.minimize(600);
      // Maxwell-Boltzmann velocities at T0, then zero COM momentum
      for (let i = 0; i < 2 * N; i++) v[i] = this.gaussian();
      this.zeroMomentum();
      this.rescaleToT(this.T0);
      this.computeForces();      // prime forces for velocity-Verlet
      this.E0 = this.energy();   // conservation baseline
    }

    /** Capped-displacement steepest descent on the potential energy. */
    minimize(iter) {
      const N = this.N, L = this.L, x = this.x, maxDisp = 0.02;
      for (let it = 0; it < iter; it++) {
        this.computeForces();
        let fmax = 0; for (let i = 0; i < 2 * N; i++) { const a = Math.abs(this.f[i]); if (a > fmax) fmax = a; }
        if (fmax < 1e-4) break;
        const sc = maxDisp / fmax;                  // normalize so largest move = maxDisp
        for (let i = 0; i < 2 * N; i++) x[i] += sc * this.f[i];
        for (let i = 0; i < 2 * N; i++) { let xi = x[i] % L; if (xi < 0) xi += L; x[i] = xi; }
      }
    }

    relaxOverlaps() {
      const N = this.N, L = this.L, x = this.x, min = 0.9 * this.sig, min2 = min * min;
      for (let i = 0; i < N; i++) for (let j = i + 1; j < N; j++) {
        let dx = x[2 * i] - x[2 * j], dy = x[2 * i + 1] - x[2 * j + 1];
        dx -= L * Math.round(dx / L); dy -= L * Math.round(dy / L);
        const r2 = dx * dx + dy * dy;
        if (r2 < min2 && r2 > 1e-9) {
          const r = Math.sqrt(r2), push = 0.5 * (min - r) / r;
          x[2 * i] += push * dx; x[2 * i + 1] += push * dy;
          x[2 * j] -= push * dx; x[2 * j + 1] -= push * dy;
        }
      }
    }

    gaussian() { // Box-Muller
      let u = 0, w = 0; while (u === 0) u = Math.random(); while (w === 0) w = Math.random();
      return Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * w);
    }

    zeroMomentum() {
      const N = this.N, v = this.v; let px = 0, py = 0;
      for (let i = 0; i < N; i++) { px += v[2 * i]; py += v[2 * i + 1]; }
      px /= N; py /= N;
      for (let i = 0; i < N; i++) { v[2 * i] -= px; v[2 * i + 1] -= py; }
    }

    rescaleToT(T) {
      const cur = this.temperature();
      if (cur > 1e-12) { const s = Math.sqrt(T / cur); for (let i = 0; i < 2 * this.N; i++) this.v[i] *= s; }
    }

    /** Compute LJ + bond forces into this.f; accumulate potential energy. */
    computeForces() {
      const N = this.N, L = this.L, x = this.x, f = this.f;
      const rc2 = this.rc * this.rc, eps = this.eps, sig2 = this.sig * this.sig;
      for (let i = 0; i < 2 * N; i++) f[i] = 0;
      let pe = 0;
      // pairwise LJ (all pairs, minimum image)
      for (let i = 0; i < N; i++) {
        const xi = x[2 * i], yi = x[2 * i + 1];
        for (let j = i + 1; j < N; j++) {
          let dx = xi - x[2 * j], dy = yi - x[2 * j + 1];
          dx -= L * Math.round(dx / L); dy -= L * Math.round(dy / L);
          const r2 = dx * dx + dy * dy;
          if (r2 < rc2 && r2 > 1e-12) {
            const inv2 = sig2 / r2, inv6 = inv2 * inv2 * inv2, inv12 = inv6 * inv6;
            pe += 4 * eps * (inv12 - inv6) - this.ecut;
            const fmag = 24 * eps * (2 * inv12 - inv6) / r2;       // (−dU/dr)/r
            f[2 * i] += fmag * dx; f[2 * i + 1] += fmag * dy;
            f[2 * j] -= fmag * dx; f[2 * j + 1] -= fmag * dy;
          }
        }
      }
      // harmonic bonds
      const k = this.kBond, r0 = this.r0Bond;
      for (let b = 0; b < this.bonds.length; b++) {
        const i = this.bonds[b].i, j = this.bonds[b].j;
        let dx = x[2 * i] - x[2 * j], dy = x[2 * i + 1] - x[2 * j + 1];
        dx -= L * Math.round(dx / L); dy -= L * Math.round(dy / L);
        const r = Math.sqrt(dx * dx + dy * dy) || 1e-9;
        pe += 0.5 * k * (r - r0) * (r - r0);
        const fmag = -k * (r - r0) / r;
        f[2 * i] += fmag * dx; f[2 * i + 1] += fmag * dy;
        f[2 * j] -= fmag * dx; f[2 * j + 1] -= fmag * dy;
      }
      this.pe = pe;
    }

    /** One nano-zeqond step: velocity-Verlet with pulse-modulated dt (KO42). */
    step() {
      const N = this.N, x = this.x, v = this.v, f = this.f, L = this.L;
      const Rt = 1 + this.alpha * Math.sin(2 * Math.PI * this.f_pulse * this.wallT + this.phi0);
      const dt = this.dtSim * Rt, half = 0.5 * dt;
      for (let i = 0; i < 2 * N; i++) { v[i] += half * f[i]; x[i] += v[i] * dt; }   // half-kick + drift (m=1)
      // periodic wrap
      for (let i = 0; i < 2 * N; i++) { let xi = x[i] % L; if (xi < 0) xi += L; x[i] = xi; }
      this.computeForces();                                                          // new forces
      for (let i = 0; i < 2 * N; i++) v[i] += half * f[i];                            // second half-kick
      this.steps++;
      this.simTime += this.simDtPerStep;
    }

    block(B) { for (let k = 0; k < B; k++) this.step(); }

    kinetic() { let ke = 0; for (let i = 0; i < 2 * this.N; i++) ke += this.v[i] * this.v[i]; return 0.5 * ke; }
    /** Temperature from equipartition: 2-D, dof = 2N−2 (COM momentum fixed). */
    temperature() { const dof = 2 * this.N - 2; return dof > 0 ? 2 * this.kinetic() / dof : 0; }
    energy() { return this.kinetic() + (this.pe || 0); }      // total (per-system, reduced units)
    momentum() { let px = 0, py = 0; for (let i = 0; i < this.N; i++) { px += this.v[2 * i]; py += this.v[2 * i + 1]; } return Math.sqrt(px * px + py * py); }
    driftPct() { if (this.E0 == null || Math.abs(this.E0) < 1e-12) return 0; return Math.abs((this.energy() - this.E0) / this.E0) * 100; }
    phase() { return (this.wallT % this.tauZ) / this.tauZ; }

    /** Advance wall time, run B steps, print on each zeqond boundary. */
    advance(dtWall, B) {
      this.wallT += dtWall;
      this.block(B);
      const z = Math.floor(this.wallT / this.tauZ);
      if (z > this.lastZeqond) {
        this.lastZeqond = z;
        const p = {
          zeqond: z, wallT: this.wallT, simTime: this.simTime, steps: this.steps,
          KE: this.kinetic(), PE: this.pe || 0, E: this.energy(),
          T: this.temperature(), driftPct: this.driftPct(), momentum: this.momentum(),
          phase: this.phase(),
        };
        this.prints.push(p);
        if (this.prints.length > 512) this.prints.shift();
        return p;
      }
      return null;
    }

    snapshot() {
      return {
        wallT: this.wallT, simTime: this.simTime, steps: this.steps,
        KE: this.kinetic(), PE: this.pe || 0, E: this.energy(), T: this.temperature(),
        driftPct: this.driftPct(), momentum: this.momentum(), phase: this.phase(),
        pulse: 1 + this.alpha * Math.sin(2 * Math.PI * this.f_pulse * this.wallT + this.phi0),
        N: this.N, L: this.L, x: this.x, bonds: this.bonds,
      };
    }
  }

  ZeqMdEngine.PRESETS = PRESETS;
  ZeqMdEngine.F_PULSE = F_PULSE;
  ZeqMdEngine.TAU_Z = TAU_Z;
  return ZeqMdEngine;
}));
