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