Implement UKS nuclear gradients for Laplacian-dependent meta-GGAs - #229
Open
susilehtola wants to merge 2 commits into
Open
susilehtola wants to merge 2 commits into
susilehtola wants to merge 2 commits into
Conversation
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. One term of its nuclear derivative contracts the BARE
d_x B against 2 P lapl(B); the code contracted (2 P d_x B) against the
bare lapl(B) instead.
The two are equal only after summing over every basis function, so the
total gradient was right and translational invariance held -- but the sum
in exc_grad_local_work_ is restricted to the shells on one atom, so the
force on each individual atom was wrong. Measured against a central
difference of eval_exc over a displaced molecule and displaced basis
(benzene/cc-pVDZ/UltraFine, atom 0, worst of the three components):
before after
MGGA_X_BR89 1.595e-01 1.604e-09
MGGA_X_R2SCANL 1.383e-01 9.365e-08
MGGA_X_SCANL 1.608e-01 --
Three independent Laplacian functionals fail identically beforehand,
including BR89, which has none of SCAN's numerical difficulty; the sign
of the error also differs between them, which no fixed numerical defect
would do. LDA (5.1e-09), GGA (2.1e-09) and tau-only meta-GGA (4.3e-09)
are unaffected and unchanged, as are energies and potentials at every
rung.
The matrix contracted with the Laplacian 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.
RKS only: exc_grad already refuses needs_laplacian with is_uks.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
exc_grad refused the spin-polarized case outright: "UKS Gradients + Laplacian Dependent MGGAs is Not Yet Implemented" Everything below the assembly was already in place -- eval_uvvar_mgga_uks computes lapl for both spins, and lapl/vlapl are allocated per spin -- so only the gradient assembly and the guard were missing. The spin transformation mirrors the LDA, GGA and tau blocks around it: rho_alpha = (rho_s + rho_z)/2, hence the 0.5 (v+ + v-) N + 0.5 (v+ - v-) Z combination, with fac P lapl(B) built for both the scalar and the spin density. One trap: the tau branch immediately above reassigns d2_term_* to the Z hessian term when polarized, and the laplacian block consumed it. Reusing it here would have paired the Z hessian term with the N laplacian contribution. This block now rebuilds hN_* and hZ_* explicitly instead. Verified against a central difference of eval_exc over a displaced molecule and displaced basis, cytosine (doublet) / cc-pVDZ / UltraFine, MGGA_X_R2SCANL+MGGA_C_R2SCANL, atom 0: d/dx analytic -1.5224974586e-01 FD -1.5224974458e-01 rel 8.4e-09 d/dy analytic 3.3558979268e-01 FD 3.3558978846e-01 rel 1.3e-08 d/dz both zero to 1e-10 (planar) Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Builds on #228 (the Laplacian gradient fix); review that first, or this diff will
look like it contains both.
What
eval_exc_gradrefused the spin-polarized case outright:Everything below the assembly was already in place —
eval_uvvar_mgga_ukscomputes
laplfor both spins (lapls,laplz, split intolapl[2i]/lapl[2i+1]), andlapl/vlaplare allocated per spin — so only the gradientassembly and the guard were missing.
How
The spin transformation mirrors the LDA, GGA and tau blocks around it:
rho_alpha = (rho_s + rho_z)/2, hence0.5 (v+ + v-) N + 0.5 (v+ - v-) Z, withfac P lapl(B)built for both thescalar and the spin density (
xLmat,xZLmat; the extra blocks already fit inhost_data.zmat, which is sized5 * spin_dim_scalwhenneeds_laplacian).One trap worth flagging for review: the tau branch immediately above
reassigns
d2_term_*to the Z hessian term when polarized, and thelaplacian block consumed it directly. Reusing it would have paired the Z hessian
term with the N laplacian contribution — a wrong-index bug of exactly the kind
#228 fixes, sitting behind the exception. This block rebuilds
hN_*andhZ_*explicitly instead.
Verification
eval_exc_gradagainst a central difference ofeval_excover a displacedmolecule and displaced basis, so the grid and its partition weights ride on
the atoms and the comparison shares no algebra with the gradient code.
Cytosine (doublet) / cc-pVDZ / UltraFine,
MGGA_X_R2SCANL+MGGA_C_R2SCANL,atom 0, h = 1e-4:
RKS is unchanged and still verifies at 1.6e-09 (BR89) / 9.4e-08 (r2SCANL).
gauxc_test: all 3640979 assertions in 56 test cases pass.Not covered
The device path keeps its own, separate restriction
(
"Device EXC Gradients + Laplacian Dependent MGGAs Not Yet Implemented"),which this does not touch.
🤖 Generated with Claude Code