Molecular dynamics

export const UNITS = {
  reduced: { name: 'reduced', kB: 1, ke: 1, accel: 1, pConv: 1 },
  real: { name: 'real (Angstrom, kcal/mol, amu, ps, e)', kB: 0.0019872041, ke: 332.0637, accel: 418.4, pConv: 68568.4, check: 'not confirmed' },
};

export function erfc(x) { const z = Math.abs(x), t = 1 / (1 + 0.5 * z), r = t * Math.exp(-z * z - 1.26551223 + t * (1.00002368 + t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 + t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277))))))))); return x >= 0 ? r : 2 - r; }

sys.sig = sys.types.map((a) => sys.types.map((b) => (a.sigma + b.sigma) / 2)); sys.eps = sys.types.map((a) => sys.types.map((b) => Math.sqrt(a.eps * b.eps)));
sys.eShift = sys.types.map((_, a) => sys.types.map((__, b) => { if (!sys.shift) return 0; const s = (sys.sig[a][b] / sys.rc) ** 6; return 4 * sys.eps[a][b] * (s * s - s); }));

/** Kinetic energy in the energy unit of the system. */
export function kinetic(s) { let k = 0; for (let i = 0; i < s.N; i++) k += s.mass[i] * (s.v[3 * i] ** 2 + s.v[3 * i + 1] ** 2 + s.v[3 * i + 2] ** 2); return 0.5 * k / s.U.accel; }

export function temperature(s) { return 2 * kinetic(s) / (dof(s) * s.U.kB); }

const pair = (i, j) => {
  let dx = mi(x[3 * i] - x[3 * j], box[0]), dy = mi(x[3 * i + 1] - x[3 * j + 1], box[1]), dz = mi(x[3 * i + 2] - x[3 * j + 2], box[2]);
  const r2 = dx * dx + dy * dy + dz * dz, key = i < j ? i * N + j : j * N + i, excluded = s.excl.size && s.excl.has(key);
  if (r2 > rc2 && !excluded) return;
  let fr = 0; // force / r
  if (excluded) { if (alpha && q[i] && q[j]) { const r = Math.sqrt(r2), er = erf(alpha * r) / r; excl -= ke * q[i] * q[j] * er; fr = ke * q[i] * q[j] * (twoA * Math.exp(-alpha * alpha * r2) - er) / r2; } }
  else {
    const a = s.type[i], b = s.type[j], sg = s.sig[a][b], ep = s.eps[a][b];
    if (ep > 0) { const s2 = sg * sg / r2, s6 = s2 * s2 * s2, e = 4 * ep * (s6 * s6 - s6) - s.eShift[a][b]; lj += e; fr += 24 * ep * (2 * s6 * s6 - s6) / r2; }
    if (alpha && q[i] && q[j]) { const r = Math.sqrt(r2), ec = erfc(alpha * r); coulReal += ke * q[i] * q[j] * ec / r; fr += ke * q[i] * q[j] * (ec / r + twoA * Math.exp(-alpha * alpha * r2)) / r2; }
  }
  if (fr === 0) return;
  W += fr * r2; f[3 * i] += fr * dx; f[3 * i + 1] += fr * dy; f[3 * i + 2] += fr * dz; f[3 * j] -= fr * dx; f[3 * j + 1] -= fr * dy; f[3 * j + 2] -= fr * dz;
};

let eb = 0, ea = 0, ed = 0;
for (const [i, j, k, r0] of s.bonds) { const dx = mi(x[3 * i] - x[3 * j], box[0]), dy = mi(x[3 * i + 1] - x[3 * j + 1], box[1]), dz = mi(x[3 * i + 2] - x[3 * j + 2], box[2]), r = Math.hypot(dx, dy, dz), d = r - r0; eb += 0.5 * k * d * d; const fr = -k * d / (r || 1); W += fr * r * r; f[3 * i] += fr * dx; f[3 * i + 1] += fr * dy; f[3 * i + 2] += fr * dz; f[3 * j] -= fr * dx; f[3 * j + 1] -= fr * dy; f[3 * j + 2] -= fr * dz; }

for (const [i, j, l, k, th0] of s.angles) {   // i-j-l with j the vertex
  const ax = mi(x[3 * i] - x[3 * j], box[0]), ay = mi(x[3 * i + 1] - x[3 * j + 1], box[1]), az = mi(x[3 * i + 2] - x[3 * j + 2], box[2]), bx = mi(x[3 * l] - x[3 * j], box[0]), by = mi(x[3 * l + 1] - x[3 * j + 1], box[1]), bz = mi(x[3 * l + 2] - x[3 * j + 2], box[2]);
  const ra = Math.hypot(ax, ay, az), rb = Math.hypot(bx, by, bz), c = Math.max(-1, Math.min(1, (ax * bx + ay * by + az * bz) / (ra * rb))), th = Math.acos(c), sn = Math.max(Math.sqrt(1 - c * c), 1e-8), dth = th - th0;
  ea += 0.5 * k * dth * dth; const g = k * dth / sn;   // dE/dcos = -k dth / sin
  const fi = [g * (bx / (ra * rb) - c * ax / (ra * ra)), g * (by / (ra * rb) - c * ay / (ra * ra)), g * (bz / (ra * rb) - c * az / (ra * ra))], fl = [g * (ax / (ra * rb) - c * bx / (rb * rb)), g * (ay / (ra * rb) - c * by / (rb * rb)), g * (az / (ra * rb) - c * bz / (rb * rb))];
  for (let d = 0; d < 3; d++) { f[3 * i + d] += fi[d]; f[3 * l + d] += fl[d]; f[3 * j + d] -= fi[d] + fl[d]; }
  W += fi[0] * ax + fi[1] * ay + fi[2] * az + fl[0] * bx + fl[1] * by + fl[2] * bz;
}

for (const [i, j, k, l, kd, n, delta] of s.dihedrals) {
  const P = (a, b) => [mi(x[3 * b] - x[3 * a], box[0]), mi(x[3 * b + 1] - x[3 * a + 1], box[1]), mi(x[3 * b + 2] - x[3 * a + 2], box[2])], cr = (u, v) => [u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]], dt = (u, v) => u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
  const b1 = P(i, j), b2 = P(j, k), b3 = P(k, l), n1 = cr(b1, b2), n2 = cr(b2, b3), lb2 = Math.sqrt(dt(b2, b2)), n1s = dt(n1, n1), n2s = dt(n2, n2);
  if (n1s < 1e-12 || n2s < 1e-12) continue;
  const phi = Math.atan2(lb2 * dt(b1, n2), dt(n1, n2)); ed += kd * (1 + Math.cos(n * phi - delta));
  const dE = -kd * n * Math.sin(n * phi - delta), gi = n1.map((v) => -lb2 / n1s * v), gl = n2.map((v) => lb2 / n2s * v), p = dt(b1, b2) / (lb2 * lb2), r = dt(b3, b2) / (lb2 * lb2);
  const gj = gi.map((v, d) => -v - p * v + r * gl[d]), gk = gl.map((v, d) => -v + p * gi[d] - r * gl[d]);
  for (let d = 0; d < 3; d++) { f[3 * i + d] -= dE * gi[d]; f[3 * j + d] -= dE * gj[d]; f[3 * k + d] -= dE * gk[d]; f[3 * l + d] -= dE * gl[d]; }
}

const km = s.ewald.kmax, V = box[0] * box[1] * box[2], pre = 2 * Math.PI * ke / V, a4 = 1 / (4 * alpha * alpha);
if (s.ewald.pme) {                                             // smooth particle-mesh Ewald (pme.js): grid [Kx, Ky, Kz] powers of two, spline order
  const g = s.ewald.pme.grid || [32, 32, 32], r = pmeRecip({ x: s.x, q, box, alpha, grid: Array.isArray(g) ? g : [g, g, g], order: s.ewald.pme.order || 4, ke }); erec = r.E; for (let i = 0; i < 3 * N; i++) f[i] += r.F[i];
} else for (let nx = -km; nx <= km; nx++) for (let ny = -km; ny <= km; ny++) for (let nz = -km; nz <= km; nz++) {
  if (!nx && !ny && !nz) continue; const n2 = nx * nx + ny * ny + nz * nz; if (n2 > km * km) continue;
  const kx = 2 * Math.PI * nx / box[0], ky = 2 * Math.PI * ny / box[1], kz = 2 * Math.PI * nz / box[2], k2 = kx * kx + ky * ky + kz * kz, A = Math.exp(-k2 * a4) / k2;
  let Sr = 0, Si = 0; const cs = new Float64Array(N), sn = new Float64Array(N);
  for (let i = 0; i < N; i++) { if (!q[i]) continue; const ph = kx * x[3 * i] + ky * x[3 * i + 1] + kz * x[3 * i + 2]; cs[i] = Math.cos(ph); sn[i] = Math.sin(ph); Sr += q[i] * cs[i]; Si += q[i] * sn[i]; }
  erec += pre * A * (Sr * Sr + Si * Si);
  for (let i = 0; i < N; i++) { if (!q[i]) continue; const g = 2 * pre * A * q[i] * (Sr * sn[i] - Si * cs[i]); f[3 * i] += g * kx; f[3 * i + 1] += g * ky; f[3 * i + 2] += g * kz; }
}
let q2 = 0; for (let i = 0; i < N; i++) q2 += q[i] * q[i]; eself = -ke * alpha / Math.sqrt(Math.PI) * q2; W += erec;

/** Pressure in the system's pressure unit (reduced: epsilon/sigma^3; real: atm). */
export function pressure(s) { return (2 * kinetic(s) + s.W) / (3 * volume(s)) * s.U.pConv; }

/** Advance `n` steps of `dt` (ps in real units; reduced time otherwise). opts { thermostat: null | { type: 'langevin', T, gamma } | { type: 'berendsen', T, tau }, barostat: null | { P, tau, kappa }, every, onSample(s) } */
export function run(s, n, dt, { thermostat = null, barostat = null, every = 0, onSample = null } = {}) {
  const { N, v, x, f, mass, U } = s, a = U.accel;
  for (let it = 0; it < n; it++) {
    const xOld = s.constraints.length ? Float64Array.from(x) : null;
    for (let i = 0; i < N; i++) { const c = 0.5 * dt * a / mass[i]; for (let d = 0; d < 3; d++) { v[3 * i + d] += c * f[3 * i + d]; x[3 * i + d] += dt * v[3 * i + d]; } }
    if (xOld) s.shakeInfo = shake(s, s.constraints, dt, xOld);
    if (barostat) { const P = pressure(s), mu = Math.cbrt(Math.max(0.5, 1 - dt / barostat.tau * (barostat.kappa ?? 1) * (barostat.P - P))); for (let d = 0; d < 3; d++) s.box[d] *= mu; for (let i = 0; i < 3 * N; i++) x[i] *= mu; if (s.rc * 2 > Math.min(...s.box)) throw new Error('The box shrank below twice the cutoff.'); }
    computeForces(s);
    for (let i = 0; i < N; i++) { const c = 0.5 * dt * a / mass[i]; for (let d = 0; d < 3; d++) v[3 * i + d] += c * f[3 * i + d]; }
    if (thermostat?.type === 'langevin') { const c1 = Math.exp(-thermostat.gamma * dt), c2 = Math.sqrt(1 - c1 * c1); for (let i = 0; i < N; i++) { const sd = Math.sqrt(U.kB * thermostat.T * a / mass[i]); for (let d = 0; d < 3; d++) v[3 * i + d] = c1 * v[3 * i + d] + c2 * sd * s.rng.n(); } }
    else if (thermostat?.type === 'berendsen') { const T = temperature(s); if (T > 0) { const lam = Math.sqrt(1 + dt / thermostat.tau * (thermostat.T / T - 1)); for (let i = 0; i < 3 * N; i++) v[i] *= lam; } }
    s.step++; s.time += dt;
    if (every && onSample && s.step % every === 0) onSample(s);
  }
  return s;
}

/** Radial distribution accumulator: g = rdf(); rdf.add(sys, typeA?, typeB?); rdf.result(). rmax defaults to half the smallest box side. */
export function rdf(nbins = 100, rmax = null) {
  const hist = new Float64Array(nbins); let frames = 0, rm = rmax, nA = 0, nB = 0, vol = 0, same = true;
  return {
    add(s, ta = null, tb = null) {
      rm = rm ?? Math.min(...s.box) / 2; const dr = rm / nbins; const ia = [], ib = []; for (let i = 0; i < s.N; i++) { if (ta === null || s.types[s.type[i]].name === ta) ia.push(i); if (tb === null || s.types[s.type[i]].name === tb) ib.push(i); }
      same = ta === tb; nA = ia.length; nB = ib.length; vol = volume(s);
      for (const i of ia) for (const j of ib) { if (i === j) continue; const r = Math.hypot(mi(s.x[3 * i] - s.x[3 * j], s.box[0]), mi(s.x[3 * i + 1] - s.x[3 * j + 1], s.box[1]), mi(s.x[3 * i + 2] - s.x[3 * j + 2], s.box[2])); if (r < rm) hist[Math.floor(r / dr)]++; }
      frames++;
    },
    result() { const dr = rm / nbins, r = [], g = []; for (let b = 0; b < nbins; b++) { const shell = 4 / 3 * Math.PI * (((b + 1) * dr) ** 3 - (b * dr) ** 3), ideal = shell * (same ? (nA - 1) : nB) / vol * nA; r.push((b + 0.5) * dr); g.push(hist[b] / (frames * ideal || 1)); } return { r, g, frames }; },
  };
}

/** Mean-square displacement from unwrapped coordinates: m = msd(sys); each sample call m.add(sys, time); fit D = slope/6 over [from, to] of the stored time. */
export function msd(s0) {
  const N = s0.N, ref = Float64Array.from(s0.x), cur = Float64Array.from(s0.x), prev = Float64Array.from(s0.x), pts = [];
  return {
    add(s, t) { for (let i = 0; i < 3 * N; i++) { let d = s.x[i] - prev[i]; const L = s.box[i % 3]; d -= L * Math.round(d / L); cur[i] += d; prev[i] = s.x[i]; } let m = 0; for (let i = 0; i < 3 * N; i++) m += (cur[i] - ref[i]) ** 2; pts.push([t, m / N]); },
    points: pts,
    diffusion(from = 0.3, to = 1) { const a = pts.slice(Math.floor(pts.length * from), Math.ceil(pts.length * to)); if (a.length < 3) return null; const n = a.length, st = a.reduce((x, p) => x + p[0], 0), sm = a.reduce((x, p) => x + p[1], 0), stt = a.reduce((x, p) => x + p[0] * p[0], 0), stm = a.reduce((x, p) => x + p[0] * p[1], 0); return (n * stm - st * sm) / (n * stt - st * st) / 6; },
  };
}

/** Mass density of the system in g/cm3 (real units only). */
export function density(s) { if (s.U.name === 'reduced') throw new Error('Density in g/cm3 needs real units.'); let m = 0; for (let i = 0; i < s.N; i++) m += s.mass[i]; return m / 0.602214076 / volume(s); }   // amu per Angstrom^3 -> g/cm3 : 1.66054 per amu/A^3

/**
 * pme.js — smooth particle-mesh Ewald (SPME) reciprocal sum and SHAKE bond constraints (0.8 Track CH, CH.25; DOM-free, no imports).
 * Clean room from the published papers: Essmann, Perera, Berkowitz, Darden, Lee and Pedersen (1995) J. Chem. Phys. 103, 8577 for SPME (cardinal B-splines,
 * charge spreading on a grid, FFT convolution with the B-spline moduli, forces from the spline derivatives), Ryckaert, Ciccotti and Berendsen (1977) for SHAKE.
 *   pmeRecip({ x, q, box, alpha, grid: [Kx, Ky, Kz] (powers of two), order = 4, ke }) → { E, F (Float64Array 3N, force = -dE/dr), W (virial term, = E as for the plain sum) }
 *   shake(sys, constraints [[i, j, d]], dt, xOld, tol) adjusts positions after an unconstrained step so each |r_i - r_j| = d, and the velocities with them.
 * Not built: RATTLE velocity constraints separately (the velocity correction is the SHAKE displacement over dt), triclinic cells, non-power-of-two grids.
 */
function fft1(re, im, inv) {                      // in place radix-2, unnormalised; inv = +1 for exp(+i), -1 for exp(-i)
  const n = re.length;
  for (let i = 1, j = 0; i < n; i++) { let bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; if (i < j) { let t = re[i]; re[i] = re[j]; re[j] = t; t = im[i]; im[i] = im[j]; im[j] = t; } }
  for (let len = 2; len <= n; len <<= 1) {
    const ang = inv * 2 * Math.PI / len, wr = Math.cos(ang), wi = Math.sin(ang);
    for (let i = 0; i < n; i += len) { let cr = 1, ci = 0; for (let k = 0; k < len / 2; k++) {
      const a = i + k, b = i + k + len / 2, tr = re[b] * cr - im[b] * ci, ti = re[b] * ci + im[b] * cr;
      re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti; const nr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = nr; } }
  }
}

/** Cardinal B-spline M_n(u) (Essmann eq. 4.1). */
export function bspline(n, u) {
  if (u <= 0 || u >= n) return 0;
  if (n === 2) return 1 - Math.abs(u - 1);
  return u / (n - 1) * bspline(n - 1, u) + (n - u) / (n - 1) * bspline(n - 1, u - 1);
}

function modulus2(K, n) {                          // |b(m)|^2 for m = 0..K-1
  const out = new Float64Array(K);
  for (let m = 0; m < K; m++) { let re = 0, im = 0; for (let k = 0; k <= n - 2; k++) { const M = bspline(n, k + 1), a = 2 * Math.PI * m * k / K; re += M * Math.cos(a); im += M * Math.sin(a); } const d = re * re + im * im; out[m] = d > 1e-12 ? 1 / d : 0; }
  // an interior zero (odd order at K/2) is patched by its neighbours, the usual practice
  for (let m = 0; m < K; m++) if (out[m] === 0) out[m] = 0.5 * (out[(m + K - 1) % K] + out[(m + 1) % K]);
  return out;
}

export function pmeRecip({ x, q, box, alpha, grid, order = 4, ke = 1 }) {
  const N = q.length, K = grid, n = order, V = box[0] * box[1] * box[2], tot = K[0] * K[1] * K[2];
  if (K.some((k) => (k & (k - 1)) !== 0 || k < 2 * n)) throw new Error('PME grid sizes must be powers of two, at least twice the spline order');
  const Qr = new Float64Array(tot), Qi = new Float64Array(tot), idx = [], wt = [], dwt = [];
  for (let i = 0; i < N; i++) {
    const I = [], W = [], DW = [];
    for (let d = 0; d < 3; d++) {
      let w = x[3 * i + d] / box[d]; w -= Math.floor(w); const u = w * K[d], fl = Math.floor(u), fr = u - fl, ii = new Int32Array(n), ww = new Float64Array(n), dd = new Float64Array(n);
      for (let j = 0; j < n; j++) { ii[j] = ((fl - j) % K[d] + K[d]) % K[d]; ww[j] = bspline(n, fr + j); dd[j] = dbspline(n, fr + j); }
      I.push(ii); W.push(ww); DW.push(dd);
    }
    idx.push(I); wt.push(W); dwt.push(DW);
    if (!q[i]) continue;
    for (let a = 0; a < n; a++) for (let b = 0; b < n; b++) for (let c = 0; c < n; c++) Qr[(I[0][a] * K[1] + I[1][b]) * K[2] + I[2][c]] += q[i] * wt[i][0][a] * wt[i][1][b] * wt[i][2][c];
  }
  const Q0 = Float64Array.from(Qr);
  fft3(Qr, Qi, K, -1);                              // Q^(m) = sum_k Q(k) exp(-2 pi i m k / K)
  const B = [modulus2(K[0], n), modulus2(K[1], n), modulus2(K[2], n)], pi2a = Math.PI * Math.PI / (alpha * alpha);
  let E = 0;
  for (let a = 0; a < K[0]; a++) { const ma = (a <= K[0] / 2 ? a : a - K[0]) / box[0];
    for (let b = 0; b < K[1]; b++) { const mb = (b <= K[1] / 2 ? b : b - K[1]) / box[1];
      for (let c = 0; c < K[2]; c++) {
        const i = (a * K[1] + b) * K[2] + c, mc = (c <= K[2] / 2 ? c : c - K[2]) / box[2], m2 = ma * ma + mb * mb + mc * mc;
        if (m2 === 0) { Qr[i] = 0; Qi[i] = 0; continue; }
        const G = ke / (Math.PI * V) * Math.exp(-pi2a * m2) / m2 * B[0][a] * B[1][b] * B[2][c];
        E += 0.5 * G * (Qr[i] * Qr[i] + Qi[i] * Qi[i]); Qr[i] *= G; Qi[i] *= G;
      } } }
  fft3(Qr, Qi, K, +1);                              // theta(k) = sum_m G Q^(m) exp(+2 pi i m k / K)
  const F = new Float64Array(3 * N);
  for (let i = 0; i < N; i++) {
    if (!q[i]) continue; const I = idx[i], W = wt[i], D = dwt[i];
    for (let a = 0; a < n; a++) for (let b = 0; b < n; b++) for (let c = 0; c < n; c++) {
      const th = Qr[(I[0][a] * K[1] + I[1][b]) * K[2] + I[2][c]];
      F[3 * i] -= q[i] * th * D[0][a] * W[1][b] * W[2][c] * K[0] / box[0];
      F[3 * i + 1] -= q[i] * th * W[0][a] * D[1][b] * W[2][c] * K[1] / box[1];
      F[3 * i + 2] -= q[i] * th * W[0][a] * W[1][b] * D[2][c] * K[2] / box[2];
    }
  }
  void Q0; return { E, F, W: E };
}

/** SHAKE. sys { N, box, x, v, mass } (typed arrays), constraints [[i, j, d]], xOld = positions before the step (the reference bond directions), dt. Corrects x and v in place. */
export function shake(sys, constraints, dt, xOld, tol = 1e-10, maxIter = 500) {
  const { x, v, mass, box } = sys, mi = (d, L) => d - L * Math.round(d / L); let worst = 0, iters = 0;
  for (let it = 0; it < maxIter; it++) {
    let done = true; iters = it + 1;
    for (const [i, j, d] of constraints) {
      const rx = mi(x[3 * i] - x[3 * j], box[0]), ry = mi(x[3 * i + 1] - x[3 * j + 1], box[1]), rz = mi(x[3 * i + 2] - x[3 * j + 2], box[2]);
      const diff = d * d - (rx * rx + ry * ry + rz * rz); if (Math.abs(diff) <= 2 * tol * d * d) continue; done = false;
      const sx = mi(xOld[3 * i] - xOld[3 * j], box[0]), sy = mi(xOld[3 * i + 1] - xOld[3 * j + 1], box[1]), sz = mi(xOld[3 * i + 2] - xOld[3 * j + 2], box[2]);
      const mu = 1 / mass[i] + 1 / mass[j], g = diff / (2 * mu * (rx * sx + ry * sy + rz * sz));
      const c = [g * sx, g * sy, g * sz];
      for (let k = 0; k < 3; k++) { x[3 * i + k] += c[k] / mass[i]; x[3 * j + k] -= c[k] / mass[j]; v[3 * i + k] += c[k] / mass[i] / dt; v[3 * j + k] -= c[k] / mass[j] / dt; }
    }
    if (done) break;
  }
  for (const [i, j, d] of constraints) { const r = Math.hypot(mi(x[3 * i] - x[3 * j], box[0]), mi(x[3 * i + 1] - x[3 * j + 1], box[1]), mi(x[3 * i + 2] - x[3 * j + 2], box[2])); worst = Math.max(worst, Math.abs(r - d) / d); }
  return { iterations: iters, worst };
}

/** Run statistics from samples: mean and standard error by block averaging in 5 blocks. */
export function stats(values) {
  const n = values.length; if (!n) return { mean: NaN, sem: NaN, n: 0 };
  const mean = values.reduce((a, b) => a + b, 0) / n, nb = Math.min(5, n), bs = Math.floor(n / nb), bm = [];
  for (let b = 0; b < nb; b++) bm.push(values.slice(b * bs, (b + 1) * bs).reduce((a, c) => a + c, 0) / bs);
  const m = bm.reduce((a, c) => a + c, 0) / nb, sem = nb > 1 ? Math.sqrt(bm.reduce((a, c) => a + (c - m) ** 2, 0) / (nb - 1) / nb) : NaN;
  return { mean, sem, n };
}

/** Component fields from an NPT run (real units): densitySamples in g/cm3, D in Angstrom^2/ps. → fields to merge into a component under `md`. */
export function toComponentProps({ T, P, densitySamples = [], D = null, forceField = '', steps = 0, dt = 0 }) {
  const d = stats(densitySamples);
  return { md: { T, P, rhoLiquid: d.n ? d.mean * 1000 : null, rhoSem: d.n ? d.sem * 1000 : null, D: D == null ? null : D * 1e-8, units: { rhoLiquid: 'kg/m3', D: 'm2/s' }, forceField, steps, dt, source: 'molecular dynamics, simulated, NOT confirmed', check: 'not confirmed' } };
}