Add long-range Coulomb Hessians and batched Hessian support

[isayev]

Summary

Enables long-range Coulomb Hessians and lifts the single-structure restriction on Hessian calculations in AIMNet2Calculator. Previously calc(data, hessian=True) raised NotImplementedError for all periodic LR Coulomb methods (dsf, ewald, pme) and supported only a single structure.

• DSF Hessian + training: adds a closed-form, autograd-twice-differentiable pure-PyTorch DSF Coulomb energy (LRCoulomb._coul_dsf_torch) and routes to it for Hessian and force/stress training (the nvalchemiops DSF kernel detaches geometry). Includes the DSF self-energy term; matches the kernel energy to 1e-4, including the periodic path (shifts @ cell).
• Ewald/PME Hessian: provides a Hessian via central finite-difference of the analytic nvalchemiops forces (_coul_nvalchemi_fd_hessian), computed in float64, summed with the NN + short-range autograd Hessian in get_derivatives. The short-range subtraction stays in the autograd graph, so the FD block supplies only the full-periodic term (correct for both subtract_sr settings).
• Batched Hessian: supports 3D batched inputs (stacked (B,N,3,N,3)) and multi-molecule mol_idx inputs (per-molecule list) via a per-structure loop over the existing validated single-structure path. Composes with the periodic paths above.
• Documents the physical distinction for vibrational/IR use: Ewald/PME Hessian is fixed-charge (long-range), DSF is relaxed-charge; the DSF force-shift makes the energy only C¹ at the cutoff.

Test Plan

• DSF Hessian: finite, symmetric; torch energy matches kernel (non-periodic and periodic).
• DSF force/stress training: torch-path forces match kernel inference.
• Ewald/PME Hessian: finite, symmetric, acoustic sum rule; FD-of-forces block matches an independent energy-FD second derivative (float64).
• Batched: 3D batch stacks and matches per-structure; multi-molecule returns a per-molecule list; batched periodic Ewald Hessian stacks.
• tests/test_calculator.py, tests/test_pbc.py, tests/test_lr.py: 227 passed.

Notes

• Requires nvalchemi-toolkit-ops >= 0.3.1 (already declared in pyproject.toml; the hybrid_forces kernel arg the Ewald/PME path uses is absent in 0.3.0).
• Multi-molecule Hessians remain unsupported within a single structure call (existing guard retained; batching dispatches per-structure instead).
• Pre-existing, unrelated: tests/test_pbc.py::TestCalculatorPBC::test_calculator_pbc_both_methods[dsf] fails under certain test orderings with RuntimeError: Cannot access data pointer ... that doesn't have storage (a FakeTensor leak from an earlier torch.compile/vmap test). This reproduces on main and passes in isolation; it is a test-isolation artifact, not a product regression from this PR.

[isayev]

Addressed multi-review findings in 7622349.

Fixed (commit 7622349):

• Precision: the Ewald/PME finite-difference Hessian block is now accumulated in its own float64 precision in get_derivatives instead of being downcast to the model Hessian dtype (was discarding the FD precision).
• State hygiene: _eval_hessian_batched snapshots and restores calculator state (_batch, _max_mol_size, _saved_for_grad, _coulomb_method, external_coulomb.method) around the per-structure loop, so a periodic batch no longer leaves the calculator permanently switched to DSF.
• Batched-split correctness: _split_mol_idx now forwards the pbc key (was dropped for multi-molecule periodic Hessians); _split_batch_dim drops precomputed nbmat/shifts/mol_idx so each subsystem rebuilds them.
• Perf: the FD Hessian loop reuses a single scratch buffer (removes 6N per-iteration full-tensor clones); dropped a redundant q.clone() in the DSF self-energy mask.
• Contract: eval/__call__ are typed dict[str, Any] with a docstring documenting the batched (leading-dim stack) and multi-molecule (per-molecule list) return regimes; the multi-molecule in-call guard message now points users to the supported batched entry.
• Docs: corrected the stale LRCoulomb class docstring (it claimed DSF/Ewald/PME Hessians were rejected) and surfaced the fixed-charge (Ewald/PME) vs relaxed-charge (DSF) Hessian caveat on set_lrcoulomb_method for vibrational/IR users; noted the DSF cutoff C¹/C² discontinuity and FD step/ewald_accuracy coupling.
• Gitignored stray Hessian/optimization run artifacts.

Deferred (not blocking): the larger contract redesign (a separate batched-Hessian entry point vs. the polymorphic dict[str, Any] return) is left as a follow-up — the current change makes the contract explicit via typing + docstring rather than restructuring the public surface.

Full suite green: tests/test_calculator.py + tests/test_pbc.py + tests/test_lr.py = 227 passed.

[isayev]

Added a matrix-free Hessian-vector product (AIMNet2Calculator.hessian_vector_product(data, vectors, *, eps=5e-4)) — folded into this PR (commits 55834ac…7894330).

Why: the dense Hessian is O(N²) memory + 3N backward passes. For negative-curvature saddle checks (Lanczos/LOBPCG) and Hessian-preconditioned optimization (CG-Newton), what's actually needed is H @ v. This provides it without ever forming the dense matrix.

How it mirrors the dense assembly (so it's verifiable against the trusted dense Hessian):

• NN + short-range + simple/dsf Coulomb + DFT-D3: one reverse-mode vjp Hv = -∂(forces·v)/∂coord (reuses the registered conv_sv_2d_sp_bwd_bwd; no triple-backward, no vmap, O(N) memory).
• Ewald/PME long-range: a directional finite difference of the analytic forces along v — 2 force evals per vector instead of the dense block's 2·3N.

Correctness gate: tests/test_hvp.py asserts hvp(v) == H.reshape(3N,3N) @ v against the already-trusted dense Hessian for simple/dsf (rtol 1e-3) and ewald/pme (rtol 5e-2, FD-limited), plus DFT-D3-attached cases, multi-vector batching, shape/validation guards.

Scope & contract: single-structure only (mirrors the dense-Hessian restriction; batched/mol_idx/compiled rejected). Returns a plain Tensor shaped like vectors — does not widen eval's contract (separate method per the architecture review). Periodic block is fixed-charge (same caveat as the dense Ewald/PME Hessian); dtype is float32 for simple/dsf, float64 for ewald/pme.

Multi-agent reviewed (gpu-pytorch / computational-chemist / architect). The review caught and we fixed: DFT-D3 curvature was being dropped from H@v for non-dsf methods (3.8e-3 → 1.35e-5 vs dense), a calculator-state leak (now snapshot/restored like the batched path), and missing species/charge validation (now shared with eval).

Full suite green: test_hvp + test_calculator + test_pbc + test_lr = 239 passed.