Fix per-atom error in the Laplacian meta-GGA nuclear gradient

[susilehtola]

Fixes #227.

The defect

lapl(rho) = 2 X . lapl(B) + 2 sum_c X_c . d_c B, with X = fac P B and
X_c = fac P d_c B (host driver, eval_uvvar_mgga_rks). Differentiating with
respect to a nuclear coordinate gives, per basis function on the moving atom,

d/dA_x lapl(rho) = -2 [ X . d_x lapl(B)            (1)
                      + d_x B . (2 P lapl(B))      (2)
                      + 2 sum_c (hess B)_xc . X_c ](3)

(1) and (3) were correct. (2) was not — the code had (2 P d_x B) . lapl(B),
i.e. the two indices swapped:

d2_term_x = xN * dlbx + xNx * lbf + 2.0*d2_term_x;   // xNx * lbf is wrong

The two forms coincide only after summing over every basis function (relabel
mu <-> nu, P symmetric). The sum in exc_grad_local_work_ is restricted to the
shells on one atom, so the total gradient came out right — and translational
invariance held, being imposed explicitly by the parent-atom subtraction — while
the force on each individual atom was wrong by ~15%. A geometry optimisation
drifts; every invariance check passes.

The change

2 P lapl(B) is not otherwise formed, so it is built with one extra
eval_xmat and one extra block of host_data.zmat, allocated only when
needs_laplacian:

xLmat = host_data.zmat.data() + 4 * spin_dim_scal * npts * nbe;
lwd->eval_xmat( npts, nbf, nbe, submat_map, xmat_fac, Ps, ldps,
                lbasis_eval, nbe, xLmat, nbe, nbe_scr );
...
d2_term_x = xN * dlbx + dbx * xL + 2.0*d2_term_x;

+24 / -4, one file, RKS only (exc_grad already refuses needs_laplacian with
is_uks).

Verification

eval_exc_grad against a central difference of eval_exc over a displaced
molecule and displaced basis — the grid and its partition weights ride on
the atoms, so this is a total derivative that shares no algebra with the
gradient code. Benzene / cc-pVDZ / UltraFine, atom 0, h = 1e-4, worst relative
deviation over the three Cartesian components:

functional	Laplacian?	before	after
SVWN5 (LDA)	no	5.1e-09	5.1e-09
PBE (GGA)	no	2.1e-09	2.1e-09
LAK (mGGA, tau)	no	6.2e-09	4.3e-09
BR89	yes	1.595e-01	1.604e-09
r2SCANL	yes	1.383e-01	9.365e-08
SCANL	yes	1.608e-01	(same path)

Three independent Laplacian functionals fail identically beforehand, including
BR89, which is plain Becke-Roussel exchange with none of SCAN's numerical
difficulty. The sign of the error differs between BR89 and r2SCANL/SCANL,
which no fixed numerical defect would do. The analytic value moves onto the
finite difference that was already there:

BR89  d/dx   analytic  -5.4963948423e-02  ->  -4.7401223822e-02
               FD      -4.7401223640e-02      -4.7401223746e-02

Energies and potentials were never affected and still are not: VXC against a
central difference of Exc in the density matrix gives the ratio 2 (the RKS
convention) to 8 digits at every rung, Laplacian included.

ctest / gauxc_test: all 3640979 assertions in 56 test cases pass,
unchanged. No test covers this path — tests/ref_data holds two Laplacian
meta-GGA reference files and neither contains a gradient, which is how the error
survived.

Not covered here

• The device path needs no matching change: the device integrator refuses
  Laplacian functionals outright
  ("Device EXC Gradients + Laplacian Dependent MGGAs Not Yet Implemented"),
  so there is no device Laplacian gradient to correct.
• UKS + Laplacian gradients remain unimplemented (the existing GAUXC_GENERIC_EXCEPTION).
  Happy to follow up with that as a separate PR.

🤖 Generated with Claude Code

[NatMatLogica]

I reproduced #227 and validated the exact change in this PR on base commit 162e4562552323a871af17ae4acd73b71071bd24.

I used GauXC's benzene/cc-pVDZ fixture and density with builtin MGGA_X_BR89, the host reference integrator, and an UltraFine unpruned Mura–Knowles grid. For atom 0 x, the base implementation gives:

eval_exc_grad                  -7.231719000000e-02
central FD of eval_exc         -6.339166098712e-02   (h = 1e-4 bohr)
relative error                  14.079973%

The FD rebuilds the molecule, displaced basis, grid, load balance and molecular weights for both geometries.

With PR commit 6afe6ff66930a7e1e5c26ca0e00fa5007a312d65, the analytic value is -6.339166144985e-02, an absolute difference of 4.63e-10 from the FD.

As an unaffected control, SVWN5 agrees with its molecular FD to 3.87e-11 on the same build.

I also instrumented the implicated host loop to export both contraction forms from that same molecular run. An AADC-NG call-site trace over the exported 12-atom state attributes atom 0 x as:

fixed gradient                         -0.0633916614498507
old xNx*lbf middle contribution        -0.0024569418587278274
correct dbx*xL middle contribution     +0.0064685916051371081
reconstructed pre-fix gradient         -0.0723171949137156
direct pre-fix gradient                -0.0723171900000000

The traced correction is -0.00892553346386493; the reconstructed old value matches the direct old run to 4.91e-9.

Both all-atom sums remain near zero (~1.4e-13), which independently confirms why a global translational check can pass while individual atomic components are wrong.

This supports the contraction diagnosis and patch already documented here. If useful, I can contribute the compact independent regression fixture and the live-state trace summary.

For disclosure: the source-attributed replay used MatLogica AADC-NG. The GauXC molecular reproducer and central-difference oracle are independent and do not require proprietary tooling. The trace replays values exported inside GauXC's real host loop; it does not operator-overload the whole GauXC integrator.