From 6afe6ff66930a7e1e5c26ca0e00fa5007a312d65 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Tue, 15 Sep 2026 12:18:06 +0300 Subject: [PATCH 1/2] Fix per-atom error in the Laplacian meta-GGA nuclear gradient 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) --- ...replicated_xc_host_integrator_exc_grad.hpp | 28 ++++++++++++++++--- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp index 64c92995..fed74040 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp @@ -203,6 +203,9 @@ void ReferenceReplicatedXCHostIntegrator:: host_data.tau .resize( spin_dim_scal * npts ); host_data.vtau.resize( spin_dim_scal * npts ); if ( needs_laplacian ) { + // One more block holds 2 P lapl(B), which the nuclear derivative + // of lapl(rho) contracts against; see the assembly below. + host_data.zmat .resize( 5 * spin_dim_scal * npts * nbe ); host_data.basis_eval.resize( 24 * npts * nbe ); // 11 + lapl_grad(3) + der3(10) host_data.lapl .resize( spin_dim_scal * npts ); host_data.vlapl.resize( spin_dim_scal * npts ); @@ -222,6 +225,7 @@ void ReferenceReplicatedXCHostIntegrator:: double* xZmat_x = nullptr; double* xZmat_y = nullptr; double* xZmat_z = nullptr; + double* xLmat = nullptr; // 2 P lapl(B), laplacian mGGAs only auto* eps = host_data.eps.data(); auto* gamma = host_data.gamma.data(); @@ -358,6 +362,13 @@ void ReferenceReplicatedXCHostIntegrator:: blas::lacpy( 'A', nbe, npts, d3basis_xxz_eval, nbe, dlgradbasis_z_eval, nbe ); blas::axpy( nbe * npts, 1., d3basis_yyz_eval, 1, dlgradbasis_z_eval, 1); blas::axpy( nbe * npts, 1., d3basis_zzz_eval, 1, dlgradbasis_z_eval, 1); + + // lapl(rho) = 2 X . lapl(B) + 2 sum_c X_c . d_c B, so one term of + // its nuclear derivative contracts the BARE d_x B against + // 2 P lapl(B). That matrix is not otherwise formed, so build it. + 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 ); } if(is_rks) lwd->eval_uvvar_mgga_rks( npts, nbe, basis_eval, dbasis_x_eval, dbasis_y_eval, @@ -557,13 +568,22 @@ void ReferenceReplicatedXCHostIntegrator:: if( needs_laplacian ) { const double vlapl_ipt = weights[ipt] * vlapl[ipt]; - const double lbf = lbasis_eval[mu_i]; const double dlbx = dlgradbasis_x_eval[mu_i]; const double dlby = dlgradbasis_y_eval[mu_i]; const double dlbz = dlgradbasis_z_eval[mu_i]; - d2_term_x = xN * dlbx + xNx * lbf + 2.0*d2_term_x; - d2_term_y = xN * dlby + xNy * lbf + 2.0*d2_term_y; - d2_term_z = xN * dlbz + xNz * lbf + 2.0*d2_term_z; + const double xL = xLmat[mu_i]; + + // d/dA lapl(rho) = -2 [ X d_x lapl(B) + d_x B (2P lapl(B)) + // + 2 sum_c (hess B)_xc X_c ]. + // The middle term pairs the bare gradient with the + // density-matrix-contracted laplacian. Pairing them the + // other way round -- (2P d_x B) against lapl(B) -- is equal + // only after summing over EVERY basis function, so it left + // the total gradient and translational invariance intact + // while putting the wrong force on each individual atom. + d2_term_x = xN * dlbx + dbx * xL + 2.0*d2_term_x; + d2_term_y = xN * dlby + dby * xL + 2.0*d2_term_y; + d2_term_z = xN * dlbz + dbz * xL + 2.0*d2_term_z; g_acc_x += vlapl_ipt * d2_term_x; g_acc_y += vlapl_ipt * d2_term_y; From 3e035f845ff1479ba358d343450677a5ad8a2785 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Tue, 15 Sep 2026 12:24:19 +0300 Subject: [PATCH 2/2] Implement UKS nuclear gradients for Laplacian-dependent meta-GGAs 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) --- ...replicated_xc_host_integrator_exc_grad.hpp | 74 +++++++++++++------ 1 file changed, 52 insertions(+), 22 deletions(-) diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp index fed74040..fa14f80d 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_grad.hpp @@ -118,9 +118,6 @@ void ReferenceReplicatedXCHostIntegrator:: // MGGA constants const bool needs_laplacian = func.needs_laplacian(); - if(needs_laplacian and is_uks) { - GAUXC_GENERIC_EXCEPTION("UKS Gradients + Laplacian Dependent MGGAs is Not Yet Implemented"); - } // Misc KS settings IntegratorSettingsEXC_GRAD exc_grad_settings; @@ -225,7 +222,8 @@ void ReferenceReplicatedXCHostIntegrator:: double* xZmat_x = nullptr; double* xZmat_y = nullptr; double* xZmat_z = nullptr; - double* xLmat = nullptr; // 2 P lapl(B), laplacian mGGAs only + double* xLmat = nullptr; // fac P lapl(B), laplacian mGGAs only + double* xZLmat = nullptr; auto* eps = host_data.eps.data(); auto* gamma = host_data.gamma.data(); @@ -369,6 +367,11 @@ void ReferenceReplicatedXCHostIntegrator:: 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 ); + if( is_uks ) { + xZLmat = xLmat + npts * nbe; + lwd->eval_xmat( npts, nbf, nbe, submat_map, xmat_fac, Pz, ldpz, + lbasis_eval, nbe, xZLmat, nbe, nbe_scr ); + } } if(is_rks) lwd->eval_uvvar_mgga_rks( npts, nbe, basis_eval, dbasis_x_eval, dbasis_y_eval, @@ -567,27 +570,54 @@ void ReferenceReplicatedXCHostIntegrator:: } if( needs_laplacian ) { - const double vlapl_ipt = weights[ipt] * vlapl[ipt]; const double dlbx = dlgradbasis_x_eval[mu_i]; const double dlby = dlgradbasis_y_eval[mu_i]; const double dlbz = dlgradbasis_z_eval[mu_i]; - const double xL = xLmat[mu_i]; - - // d/dA lapl(rho) = -2 [ X d_x lapl(B) + d_x B (2P lapl(B)) - // + 2 sum_c (hess B)_xc X_c ]. - // The middle term pairs the bare gradient with the - // density-matrix-contracted laplacian. Pairing them the - // other way round -- (2P d_x B) against lapl(B) -- is equal - // only after summing over EVERY basis function, so it left - // the total gradient and translational invariance intact - // while putting the wrong force on each individual atom. - d2_term_x = xN * dlbx + dbx * xL + 2.0*d2_term_x; - d2_term_y = xN * dlby + dby * xL + 2.0*d2_term_y; - d2_term_z = xN * dlbz + dbz * xL + 2.0*d2_term_z; - - g_acc_x += vlapl_ipt * d2_term_x; - g_acc_y += vlapl_ipt * d2_term_y; - g_acc_z += vlapl_ipt * d2_term_z; + + // The tau branch above leaves d2_term_* holding the Z term + // when polarized, so rebuild what this block needs rather + // than reuse it. + const double hN_x = d2bxx * xNx + d2bxy * xNy + d2bxz * xNz; + const double hN_y = d2bxy * xNx + d2byy * xNy + d2byz * xNz; + const double hN_z = d2bxz * xNx + d2byz * xNy + d2bzz * xNz; + + // d/dA lapl(rho) = -2 [ X d_x lapl(B) + d_x B (fac P lapl(B)) + // + 2 sum_c (hess B)_xc X_c ] + const double xL = xLmat[mu_i]; + const double lN_x = xN * dlbx + dbx * xL + 2.0*hN_x; + const double lN_y = xN * dlby + dby * xL + 2.0*hN_y; + const double lN_z = xN * dlbz + dbz * xL + 2.0*hN_z; + + if(is_rks) { + const double vlapl_ipt = weights[ipt] * vlapl[ipt]; + g_acc_x += vlapl_ipt * lN_x; + g_acc_y += vlapl_ipt * lN_y; + g_acc_z += vlapl_ipt * lN_z; + } else { + // Same n / m_z transformation the LDA, GGA and tau blocks + // above use: rho_alpha = (rho_s + rho_z)/2. + const double vlaplp_ipt = weights[ipt] * vlapl[spin_dim_scal * ipt + 0]; + const double vlaplm_ipt = weights[ipt] * vlapl[spin_dim_scal * ipt + 1]; + const double vlapln_ipt = vlaplp_ipt + vlaplm_ipt; + const double vlaplz_ipt = vlaplp_ipt - vlaplm_ipt; + + g_acc_x += 0.5 * vlapln_ipt * lN_x; + g_acc_y += 0.5 * vlapln_ipt * lN_y; + g_acc_z += 0.5 * vlapln_ipt * lN_z; + + const double hZ_x = d2bxx * xZx + d2bxy * xZy + d2bxz * xZz; + const double hZ_y = d2bxy * xZx + d2byy * xZy + d2byz * xZz; + const double hZ_z = d2bxz * xZx + d2byz * xZy + d2bzz * xZz; + + const double xZL = xZLmat[mu_i]; + const double lZ_x = xZ * dlbx + dbx * xZL + 2.0*hZ_x; + const double lZ_y = xZ * dlby + dby * xZL + 2.0*hZ_y; + const double lZ_z = xZ * dlbz + dbz * xZL + 2.0*hZ_z; + + g_acc_x += 0.5 * vlaplz_ipt * lZ_x; + g_acc_y += 0.5 * vlaplz_ipt * lZ_y; + g_acc_z += 0.5 * vlaplz_ipt * lZ_z; + } } } }