From 44e66aa8d1abb2a53f1a1c8c35dfbf49c8a8905b Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Sat, 19 Sep 2026 23:57:53 +0300 Subject: [PATCH] Add the XC nuclear Hessian (RKS; LDA, GGA, meta-GGA) eval_exc_hess(P, settings) returns the 3N x 3N second derivative of the XC energy at fixed density matrix, on the host replicated integrator. IntegratorSettingsEXC_HESS mirrors IntegratorSettingsEXC_GRAD: - include_weight_derivatives = true (default): the basis functions ride their atoms, the grid points ride their parent atoms, and the partition weights change. A task's energy depends on R_A - R_C only (C its parent atom), so each task is accumulated in coordinates excluding C and C's rows and columns are restored by translational invariance. The weight class adds w'' e + w' e'^T + e' w'^T per point, with the Becke or SSF weight derivatives in log form dw/w = a_C - g, g = sum_D pi_D a_D d2w/w = (a_C-g)(a_C-g)^T + B_C - sum_D pi_D [(a_D-g)(a_D-g)^T + B_D] (pi_D = P_D/Z, a_D = d ln P_D), finite as P_C -> 0. LKO throws. - include_weight_derivatives = false: the basis class alone (fixed grid and weights), the analogue of the gradient's Hellmann-Feynman option. The basis class sums generated per-function rows over the shells of each atom, contracts the two row sets as an outer product, and adds the Pulay and same-atom terms. The Pulay term is exactly bilinear in the displaced rows, alpha(u)^T W(g) beta(v), so it is contracted with GEMMs. The per-point kernels (host/gauxc_hess_kernel.hpp) and their call sites (host/gauxc_hess_call_*.inc) are generated by xckernel and should be regenerated, not edited: python -m xckernel.emitters.gauxcwriter --emit-dir src/xc_integrator/local_work_driver/host Call-site arguments are bound by name from each kernel's signature: GauXC's names (vgamma, ...) sort differently from Libxc's (vsigma, ...), and a hand-written meta-GGA call site had mis-bound five of seven arguments. The generator also proves symbolically that the assembly recipe equals the unfactorised Hessian, and generates the confocal coordinate mu_DE with its derivatives and the Becke/SSF cell functions with t = s'/s, u = s''/s from the factored ln s (stable as mu -> 1). Validation, benzene/cc-pVDZ, LDA, PBE, TPSS and LAK, against central differences (Richardson) of eval_exc_grad with weight derivatives, displacing atom + basis and rebuilding the grid: diagonal and off-diagonal blocks agree to ~1e-12 (Becke) and ~1e-9 (SSF) on the PySCF0 grid; on UltraFine Becke agrees to 8e-12 and SSF to 2e-7, where the finite-difference reference is itself asymmetric by the same amount. H is symmetric to 3e-16 and obeys the translational sum rule to 3e-15. The basis class alone agrees with finite differences over the basis function origins at the truncation level of the reference. Not included: UKS; Laplacian-dependent meta-GGAs (they need fourth collocation derivatives). The Becke weight term costs O(N^3) per grid point; SSF skips inactive atom pairs. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_0135evJ9zgNL1y8U6T9bQ3UT --- include/gauxc/xc_integrator.hpp | 4 + include/gauxc/xc_integrator/impl.hpp | 7 + .../gauxc/xc_integrator/replicated/impl.hpp | 15 + .../replicated_xc_integrator_impl.hpp | 10 + .../replicated_xc_integrator.hpp | 2 + .../xc_integrator/xc_integrator_impl.hpp | 16 + include/gauxc/xc_integrator_settings.hpp | 4 + .../host/gauxc_hess_call_cell.inc | 9 + .../host/gauxc_hess_call_egrad.inc | 11 + .../host/gauxc_hess_call_mu.inc | 5 + .../host/gauxc_hess_call_pair.inc | 11 + .../host/gauxc_hess_call_pulayW.inc | 11 + .../host/gauxc_hess_call_rows.inc | 26 + .../host/gauxc_hess_call_same.inc | 11 + .../host/gauxc_hess_kernel.hpp | 411 ++++++++++++ ...eference_replicated_xc_host_integrator.cxx | 1 + ...eference_replicated_xc_host_integrator.hpp | 8 + ...replicated_xc_host_integrator_exc_hess.hpp | 595 ++++++++++++++++++ 18 files changed, 1157 insertions(+) create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_cell.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_egrad.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_mu.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_pair.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_pulayW.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_rows.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_call_same.inc create mode 100644 src/xc_integrator/local_work_driver/host/gauxc_hess_kernel.hpp create mode 100644 src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_hess.hpp diff --git a/include/gauxc/xc_integrator.hpp b/include/gauxc/xc_integrator.hpp index 03feaf934..30dfe2fbf 100644 --- a/include/gauxc/xc_integrator.hpp +++ b/include/gauxc/xc_integrator.hpp @@ -39,6 +39,8 @@ class XCIntegrator { using exc_vxc_type_uks = std::tuple< value_type, matrix_type, matrix_type >; using exc_vxc_type_gks = std::tuple< value_type, matrix_type, matrix_type, matrix_type, matrix_type >; using exc_grad_type = std::vector< value_type >; + //! Nuclear Hessian, row-major (3*natoms) x (3*natoms). + using exc_hess_type = std::vector< value_type >; using exx_type = matrix_type; using fxc_contraction_type_rks = matrix_type; using fxc_contraction_type_uks = std::tuple< matrix_type, matrix_type >; @@ -77,6 +79,8 @@ class XCIntegrator { exc_grad_type eval_exc_grad( const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); exc_grad_type eval_exc_grad( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); + exc_hess_type eval_exc_hess( const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); + exx_type eval_exx ( const MatrixType&, const IntegratorSettingsEXX& = IntegratorSettingsEXX{} ); diff --git a/include/gauxc/xc_integrator/impl.hpp b/include/gauxc/xc_integrator/impl.hpp index 400afb7c7..2bd35a2ce 100644 --- a/include/gauxc/xc_integrator/impl.hpp +++ b/include/gauxc/xc_integrator/impl.hpp @@ -92,6 +92,13 @@ typename XCIntegrator::exc_grad_type return pimpl_->eval_exc_grad(Ps, Pz, ks_settings); }; +template +typename XCIntegrator::exc_hess_type + XCIntegrator::eval_exc_hess( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) { + if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED(); + return pimpl_->eval_exc_hess(P, ks_settings); +}; + template typename XCIntegrator::exx_type XCIntegrator::eval_exx( const MatrixType& P, diff --git a/include/gauxc/xc_integrator/replicated/impl.hpp b/include/gauxc/xc_integrator/replicated/impl.hpp index bfc95fc88..0d82e4ea9 100644 --- a/include/gauxc/xc_integrator/replicated/impl.hpp +++ b/include/gauxc/xc_integrator/replicated/impl.hpp @@ -175,6 +175,21 @@ typename ReplicatedXCIntegrator::exc_grad_type } +template +typename ReplicatedXCIntegrator::exc_hess_type + ReplicatedXCIntegrator::eval_exc_hess_( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) { + + if( not pimpl_ ) GAUXC_PIMPL_NOT_INITIALIZED(); + + const size_t n3 = 3*pimpl_->load_balancer().molecule().natoms(); + std::vector EXC_HESS( n3*n3 ); + pimpl_->eval_exc_hess( P.rows(), P.cols(), P.data(), P.rows(), + EXC_HESS.data(), ks_settings ); + + return EXC_HESS; + +} + template typename ReplicatedXCIntegrator::exc_grad_type ReplicatedXCIntegrator::eval_exc_grad_( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) { diff --git a/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp b/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp index 457315122..f020c564c 100644 --- a/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp +++ b/include/gauxc/xc_integrator/replicated/replicated_xc_integrator_impl.hpp @@ -80,6 +80,11 @@ class ReplicatedXCIntegratorImpl { virtual void eval_exc_grad_( int64_t m, int64_t n, const value_type* P, int64_t ldp, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ) = 0; + virtual void eval_exc_hess_( int64_t m, int64_t n, const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ) { + GAUXC_GENERIC_EXCEPTION("EXC Hessian Not Implemented For This Integrator"); + } + virtual void eval_exc_grad_( int64_t m, int64_t n, const value_type* P, int64_t ldps, const value_type* Pz, int64_t lpdz, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ) = 0; virtual void eval_exx_( int64_t m, int64_t n, const value_type* P, @@ -155,6 +160,11 @@ class ReplicatedXCIntegratorImpl { void eval_exc_grad( int64_t m, int64_t n, const value_type* P, int64_t ldp, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ); + void eval_exc_hess( int64_t m, int64_t n, const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ) { + eval_exc_hess_(m, n, P, ldp, EXC_HESS, ks_settings); + } + void eval_exc_grad( int64_t m, int64_t n, const value_type* Ps, int64_t ldps, const value_type* Pz, int64_t ldpz, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ); diff --git a/include/gauxc/xc_integrator/replicated_xc_integrator.hpp b/include/gauxc/xc_integrator/replicated_xc_integrator.hpp index 1ca53f917..074c41d2c 100644 --- a/include/gauxc/xc_integrator/replicated_xc_integrator.hpp +++ b/include/gauxc/xc_integrator/replicated_xc_integrator.hpp @@ -36,6 +36,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl { using exc_vxc_type_uks = typename XCIntegratorImpl::exc_vxc_type_uks; using exc_vxc_type_gks = typename XCIntegratorImpl::exc_vxc_type_gks; using exc_grad_type = typename XCIntegratorImpl::exc_grad_type; + using exc_hess_type = typename XCIntegratorImpl::exc_hess_type; using exx_type = typename XCIntegratorImpl::exx_type; using fxc_contraction_type_rks = typename XCIntegratorImpl::fxc_contraction_type_rks; using fxc_contraction_type_uks = typename XCIntegratorImpl::fxc_contraction_type_uks; @@ -56,6 +57,7 @@ class ReplicatedXCIntegrator : public XCIntegratorImpl { exc_vxc_type_gks eval_exc_vxc_ ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override; exc_grad_type eval_exc_grad_( const MatrixType&, const IntegratorSettingsXC& ) override; exc_grad_type eval_exc_grad_( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override; + exc_hess_type eval_exc_hess_( const MatrixType&, const IntegratorSettingsXC& ) override; exx_type eval_exx_ ( const MatrixType&, const IntegratorSettingsEXX& ) override; fxc_contraction_type_rks eval_fxc_contraction_ ( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& ) override; fxc_contraction_type_uks eval_fxc_contraction_ ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, const IntegratorSettingsXC&) override; diff --git a/include/gauxc/xc_integrator/xc_integrator_impl.hpp b/include/gauxc/xc_integrator/xc_integrator_impl.hpp index ba7bebebb..027207374 100644 --- a/include/gauxc/xc_integrator/xc_integrator_impl.hpp +++ b/include/gauxc/xc_integrator/xc_integrator_impl.hpp @@ -28,6 +28,7 @@ class XCIntegratorImpl { using exc_vxc_type_uks = typename XCIntegrator::exc_vxc_type_uks; using exc_vxc_type_gks = typename XCIntegrator::exc_vxc_type_gks; using exc_grad_type = typename XCIntegrator::exc_grad_type; + using exc_hess_type = typename XCIntegrator::exc_hess_type; using exx_type = typename XCIntegrator::exx_type; using fxc_contraction_type_rks = typename XCIntegrator::fxc_contraction_type_rks; using fxc_contraction_type_uks = typename XCIntegrator::fxc_contraction_type_uks; @@ -48,6 +49,12 @@ class XCIntegratorImpl { const IntegratorSettingsXC& ks_settings ) = 0; virtual exc_grad_type eval_exc_grad_( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) = 0; virtual exc_grad_type eval_exc_grad_( const MatrixType& Ps, const MatrixType& Pz, const IntegratorSettingsXC& ks_settings ) = 0; + /** Nuclear Hessian. Defaulted to a refusal rather than pure virtual: + * only the reference host integrator implements it so far, and every + * other backend would otherwise need a stub. */ + virtual exc_hess_type eval_exc_hess_( const MatrixType&, const IntegratorSettingsXC& ) { + GAUXC_GENERIC_EXCEPTION("EXC Hessian Not Implemented For This Integrator"); + } virtual exx_type eval_exx_ ( const MatrixType& P, const IntegratorSettingsEXX& settings ) = 0; virtual fxc_contraction_type_rks eval_fxc_contraction_ ( const MatrixType& P, @@ -142,6 +149,15 @@ class XCIntegratorImpl { return eval_exc_grad_(Ps, Pz, ks_settings); } + /** Integrate the EXC nuclear Hessian for RKS + * + * @param[in] P The alpha density matrix + * @returns EXC Hessian, row-major (3*natoms) x (3*natoms) + */ + exc_hess_type eval_exc_hess( const MatrixType& P, const IntegratorSettingsXC& ks_settings ) { + return eval_exc_hess_(P, ks_settings); + } + /** Integrate Exact Exchange for RHF * * @param[in] P The alpha density matrix diff --git a/include/gauxc/xc_integrator_settings.hpp b/include/gauxc/xc_integrator_settings.hpp index 1ec26d0e6..2dd5dcf21 100644 --- a/include/gauxc/xc_integrator_settings.hpp +++ b/include/gauxc/xc_integrator_settings.hpp @@ -29,4 +29,8 @@ struct IntegratorSettingsEXC_GRAD : public IntegratorSettingsKS { bool include_weight_derivatives= true; // whether to include grid weight contribution and employ translational invariance, or just use Hellmann-Feynman gradient }; +struct IntegratorSettingsEXC_HESS : public IntegratorSettingsKS { + bool include_weight_derivatives = true; // as for EXC_GRAD: moving grid (translational invariance) plus the partition-weight first and second derivatives; false gives the fixed-grid, fixed-weight (basis-only) Hessian +}; + } diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_cell.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_cell.inc new file mode 100644 index 000000000..8f00d1c23 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_cell.inc @@ -0,0 +1,9 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: cell-function call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_becke ) { + xckernel::xck_gauxc_weight_cell_becke( mu, s, t, u ); +} else { + xckernel::xck_gauxc_weight_cell_ssf( mu, s, t, u ); +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_egrad.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_egrad.inc new file mode 100644 index 000000000..6a709b2b8 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_egrad.inc @@ -0,0 +1,11 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: egrad call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_mgga ) { + xckernel::xck_gauxc_hess_egrad_mgga_tau( ROW(a,d,0,ip), ROW(a,d,1,ip), ROW(a,d,2,ip), vrho[ip], vgamma[ip], vtau[ip], de ); +} else if( is_gga ) { + xckernel::xck_gauxc_hess_egrad_gga( ROW(a,d,0,ip), ROW(a,d,1,ip), vrho[ip], vgamma[ip], de ); +} else { + xckernel::xck_gauxc_hess_egrad_lda( ROW(a,d,0,ip), vrho[ip], de ); +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_mu.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_mu.inc new file mode 100644 index 000000000..e0ee43c68 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_mu.inc @@ -0,0 +1,5 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: mu call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +xckernel::xck_gauxc_weight_mu( RD[0], RD[1], RD[2], RE[0], RE[1], RE[2], rg[0], rg[1], rg[2], mu, dmu[0], dmu[1], dmu[2], dmu[3], dmu[4], dmu[5], d2mu[0], d2mu[1], d2mu[2], d2mu[3], d2mu[4], d2mu[5], d2mu[6], d2mu[7], d2mu[8], d2mu[9], d2mu[10], d2mu[11], d2mu[12], d2mu[13], d2mu[14], d2mu[15], d2mu[16], d2mu[17], d2mu[18], d2mu[19], d2mu[20], d2mu[21], d2mu[22], d2mu[23], d2mu[24], d2mu[25], d2mu[26], d2mu[27], d2mu[28], d2mu[29], d2mu[30], d2mu[31], d2mu[32], d2mu[33], d2mu[34], d2mu[35] ); +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pair.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pair.inc new file mode 100644 index 000000000..9a7a8a909 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pair.inc @@ -0,0 +1,11 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: pair call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_mgga ) { + xckernel::xck_gauxc_hess_pair_mgga_tau( ROW(a,dx,0,ip), ROW(b,dy,0,ip), ROW(a,dx,1,ip), ROW(b,dy,1,ip), ROW(a,dx,2,ip), ROW(b,dy,2,ip), ROW(a,dx,3,ip), ROW(a,dx,4,ip), ROW(a,dx,5,ip), ROW(b,dy,3,ip), ROW(b,dy,4,ip), ROW(b,dy,5,ip), v2rho2[ip], v2rhogamma[ip], v2rhotau[ip], v2gamma2[ip], v2gammatau[ip], v2tau2[ip], vgamma[ip], h ); +} else if( is_gga ) { + xckernel::xck_gauxc_hess_pair_gga( ROW(a,dx,0,ip), ROW(b,dy,0,ip), ROW(a,dx,1,ip), ROW(b,dy,1,ip), ROW(a,dx,2,ip), ROW(a,dx,3,ip), ROW(a,dx,4,ip), ROW(b,dy,2,ip), ROW(b,dy,3,ip), ROW(b,dy,4,ip), v2rho2[ip], v2rhogamma[ip], v2gamma2[ip], vgamma[ip], h ); +} else { + xckernel::xck_gauxc_hess_pair_lda( ROW(a,dx,0,ip), ROW(b,dy,0,ip), v2rho2[ip], h ); +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pulayW.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pulayW.inc new file mode 100644 index 000000000..d34bfc2ac --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_pulayW.inc @@ -0,0 +1,11 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: pulayW call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_mgga ) { + xckernel::xck_gauxc_hess_pulayW_mgga_tau( dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], vtau[ip], weights[ip], W[0], W[1], W[2], W[3], W[4], W[5], W[8], W[10], W[12], W[15] ); +} else if( is_gga ) { + xckernel::xck_gauxc_hess_pulayW_gga( dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], weights[ip], W[0], W[1], W[2], W[3], W[4], W[8], W[12] ); +} else { + xckernel::xck_gauxc_hess_pulayW_lda( vrho[ip], weights[ip], W[0] ); +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_rows.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_rows.inc new file mode 100644 index 000000000..b82438fd4 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_rows.inc @@ -0,0 +1,26 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: rows call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_mgga ) { + xckernel::xck_gauxc_hess_rows_mgga_tau( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], dden_x[ip], dden_y[ip], dden_z[ip], F_rho, F_sigma, F_tau, Gx, Gy, Gz ); + ROW(a,d,0,ip) += F_rho; + ROW(a,d,1,ip) += F_sigma; + ROW(a,d,2,ip) += F_tau; + ROW(a,d,3,ip) += Gx; + ROW(a,d,4,ip) += Gy; + ROW(a,d,5,ip) += Gz; +} else if( is_gga ) { + xckernel::xck_gauxc_hess_rows_gga( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], dden_x[ip], dden_y[ip], dden_z[ip], F_rho, F_sigma, Gx, Gy, Gz ); + ROW(a,d,0,ip) += F_rho; + ROW(a,d,1,ip) += F_sigma; + ROW(a,d,2,ip) += Gx; + ROW(a,d,3,ip) += Gy; + ROW(a,d,4,ip) += Gz; +} else { + xckernel::xck_gauxc_hess_rows_lda( U0, U1, U2, U3, dchi, ddchi[0], ddchi[1], ddchi[2], F_rho, Gx, Gy, Gz ); + ROW(a,d,0,ip) += F_rho; + ROW(a,d,1,ip) += Gx; + ROW(a,d,2,ip) += Gy; + ROW(a,d,3,ip) += Gz; +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_call_same.inc b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_same.inc new file mode 100644 index 000000000..8f8dcaab6 --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_call_same.inc @@ -0,0 +1,11 @@ +// ==> BEGIN GENERATED CODE [xckernel gauxcwriter: same call site] <== +// Arguments are bound BY NAME from each kernel's own signature; +// regenerate rather than edit: python -m xckernel.emitters.gauxcwriter --emit-dir +if( is_mgga ) { + xckernel::xck_gauxc_hess_same_mgga_tau( xmat[k], xmat_x[k], xmat_y[k], xmat_z[k], d2c, d3c[0], d3c[1], d3c[2], dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], vtau[ip], weights[ip], sv ); +} else if( is_gga ) { + xckernel::xck_gauxc_hess_same_gga( xmat[k], xmat_x[k], xmat_y[k], xmat_z[k], d2c, d3c[0], d3c[1], d3c[2], dden_x[ip], dden_y[ip], dden_z[ip], vrho[ip], vgamma[ip], weights[ip], sv ); +} else { + xckernel::xck_gauxc_hess_same_lda( xmat[k], d2c, vrho[ip], weights[ip], sv ); +} +// ==> END GENERATED CODE <== diff --git a/src/xc_integrator/local_work_driver/host/gauxc_hess_kernel.hpp b/src/xc_integrator/local_work_driver/host/gauxc_hess_kernel.hpp new file mode 100644 index 000000000..a62d1419d --- /dev/null +++ b/src/xc_integrator/local_work_driver/host/gauxc_hess_kernel.hpp @@ -0,0 +1,411 @@ +// Machine-generated by xckernel; do not edit. +// Reproduce with: python -m xckernel.emitters.gauxcwriter --emit +// Copyright (c) 2026 Susi Lehtola. +#pragma once +#include + +namespace GauXC { +namespace xckernel { +/* machine-generated by xckernel; do not edit. */ +// lda: per-function nuclear-displacement field rows. +// Sum over the basis functions of one atom to get the perturbed +// fields rho^(A,x), grad rho^(A,x) and tau^(A,x) at this point. +// dchi/ddchi carry the -d/dr sign and the atom mask. +static inline void xck_gauxc_hess_rows_lda(double U0_u, double U1_u, double U2_u, double U3_u, double dchi_gA_u, double ddchi_gA_u_x, double ddchi_gA_u_y, double ddchi_gA_u_z, double & F_rho, double & G_x, double & G_y, double & G_z) { + F_rho = 2*U0_u*dchi_gA_u; + G_x = 2*(U0_u*ddchi_gA_u_x + U1_u*dchi_gA_u); + G_y = 2*(U0_u*ddchi_gA_u_y + U2_u*dchi_gA_u); + G_z = 2*(U0_u*ddchi_gA_u_z + U3_u*dchi_gA_u); +} + +/* machine-generated by xckernel; do not edit. */ +// lda: non-Pulay pair term, per grid point. +// The fxc quadratic form on the two perturbed field sets, PLUS +// the sigma-curvature term 2 vsigma G^A . G^B -- sigma is +// quadratic in grad rho, so d2 sigma / dA dB survives for A /= B. +// Multiply by the quadrature weight at the call site. +static inline void xck_gauxc_hess_pair_lda(double F_rho_A, double F_rho_B, double v2rho2, double & h) { + h = F_rho_A*F_rho_B*v2rho2; +} + +/* machine-generated by xckernel; do not edit. */ +// lda: Pulay term per basis-function PAIR, with the +// density-matrix factor D_uv DIVIDED OUT -- multiply by it, +// and by the quadrature weight, at the call site. +// Cannot factorise into an outer product; this is the only +// piece that needs the function pair. +static inline void xck_gauxc_hess_pulay_lda(double dchi_gA_u, double dchi_gB_v, double vrho, double w, double & s) { + s = 2*dchi_gA_u*dchi_gB_v*vrho*w; +} + +/* machine-generated by xckernel; do not edit. */ +// lda: same-atom (delta_AB) term, per basis function. +// Complete over Cartesian components; multiply by the +// quadrature weight at the call site. +static inline void xck_gauxc_hess_same_lda(double U0_u, double d2chi_g2_u, double vrho, double w, double & s) { + s = 2*U0_u*d2chi_g2_u*vrho*w; +} + +/* machine-generated by xckernel; do not edit. */ +// lda: Pulay point-weight matrix W_ab, weight included. +// Slots: 0 = displaced function, 1-3 = its displaced gradient. +// pulay_uv = sum_ab alpha_a(u) W_ab beta_b(v): contract the rows +// as matrix products rather than per function pair. +static inline void xck_gauxc_hess_pulayW_lda(double vrho, double w, double & W00) { + W00 = 2*vrho*w; +} + +/* machine-generated by xckernel; do not edit. */ +// lda: d e(r_g) / d R_A, e = the unweighted energy density, +// basis class only (points fixed). No quadrature weight. +static inline void xck_gauxc_hess_egrad_lda(double F_rho_A, double vrho, double & de) { + de = F_rho_A*vrho; +} + +/* machine-generated by xckernel; do not edit. */ +// gga: per-function nuclear-displacement field rows. +// Sum over the basis functions of one atom to get the perturbed +// fields rho^(A,x), grad rho^(A,x) and tau^(A,x) at this point. +// dchi/ddchi carry the -d/dr sign and the atom mask. +static inline void xck_gauxc_hess_rows_gga(double U0_u, double U1_u, double U2_u, double U3_u, double dchi_gA_u, double ddchi_gA_u_x, double ddchi_gA_u_y, double ddchi_gA_u_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double & F_rho, double & F_sigma, double & G_x, double & G_y, double & G_z) { + const double t0 = U0_u*ddchi_gA_u_x + U1_u*dchi_gA_u; + const double t1 = U0_u*ddchi_gA_u_y + U2_u*dchi_gA_u; + const double t2 = U0_u*ddchi_gA_u_z + U3_u*dchi_gA_u; + F_rho = 2*U0_u*dchi_gA_u; + F_sigma = 4*grad_rho_x*t0 + 4*grad_rho_y*t1 + 4*grad_rho_z*t2; + G_x = 2*t0; + G_y = 2*t1; + G_z = 2*t2; +} + +/* machine-generated by xckernel; do not edit. */ +// gga: non-Pulay pair term, per grid point. +// The fxc quadratic form on the two perturbed field sets, PLUS +// the sigma-curvature term 2 vsigma G^A . G^B -- sigma is +// quadratic in grad rho, so d2 sigma / dA dB survives for A /= B. +// Multiply by the quadrature weight at the call site. +static inline void xck_gauxc_hess_pair_gga(double F_rho_A, double F_rho_B, double F_sigma_A, double F_sigma_B, double G_A_x, double G_A_y, double G_A_z, double G_B_x, double G_B_y, double G_B_z, double v2rho2, double v2rhosigma, double v2sigma2, double vsigma, double & h) { + h = F_rho_A*F_rho_B*v2rho2 + F_sigma_A*F_sigma_B*v2sigma2 + v2rhosigma*(F_rho_A*F_sigma_B + F_rho_B*F_sigma_A) + 2*vsigma*(G_A_x*G_B_x + G_A_y*G_B_y + G_A_z*G_B_z); +} + +/* machine-generated by xckernel; do not edit. */ +// gga: Pulay term per basis-function PAIR, with the +// density-matrix factor D_uv DIVIDED OUT -- multiply by it, +// and by the quadrature weight, at the call site. +// Cannot factorise into an outer product; this is the only +// piece that needs the function pair. +static inline void xck_gauxc_hess_pulay_gga(double dchi_gA_u, double dchi_gB_v, double ddchi_gA_u_x, double ddchi_gA_u_y, double ddchi_gA_u_z, double ddchi_gB_v_x, double ddchi_gB_v_y, double ddchi_gB_v_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double w, double & s) { + const double t0 = 2*vsigma; + const double t1 = dchi_gA_u*t0; + const double t2 = dchi_gB_v*t0; + s = 2*w*(dchi_gA_u*dchi_gB_v*vrho + ddchi_gA_u_x*grad_rho_x*t2 + ddchi_gA_u_y*grad_rho_y*t2 + ddchi_gA_u_z*grad_rho_z*t2 + ddchi_gB_v_x*grad_rho_x*t1 + ddchi_gB_v_y*grad_rho_y*t1 + ddchi_gB_v_z*grad_rho_z*t1); +} + +/* machine-generated by xckernel; do not edit. */ +// gga: same-atom (delta_AB) term, per basis function. +// Complete over Cartesian components; multiply by the +// quadrature weight at the call site. +static inline void xck_gauxc_hess_same_gga(double U0_u, double U1_u, double U2_u, double U3_u, double d2chi_g2_u, double d3chi_g2_u_x, double d3chi_g2_u_y, double d3chi_g2_u_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double w, double & s) { + const double t0 = 2*vsigma; + const double t1 = U0_u*t0; + const double t2 = d2chi_g2_u*t0; + s = 2*w*(U0_u*d2chi_g2_u*vrho + U1_u*grad_rho_x*t2 + U2_u*grad_rho_y*t2 + U3_u*grad_rho_z*t2 + d3chi_g2_u_x*grad_rho_x*t1 + d3chi_g2_u_y*grad_rho_y*t1 + d3chi_g2_u_z*grad_rho_z*t1); +} + +/* machine-generated by xckernel; do not edit. */ +// gga: Pulay point-weight matrix W_ab, weight included. +// Slots: 0 = displaced function, 1-3 = its displaced gradient. +// pulay_uv = sum_ab alpha_a(u) W_ab beta_b(v): contract the rows +// as matrix products rather than per function pair. +static inline void xck_gauxc_hess_pulayW_gga(double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double w, double & W00, double & W01, double & W02, double & W03, double & W10, double & W20, double & W30) { + const double t0 = 4*vsigma*w; + const double t1 = grad_rho_x*t0; + const double t2 = grad_rho_y*t0; + const double t3 = grad_rho_z*t0; + W00 = 2*vrho*w; + W01 = t1; + W02 = t2; + W03 = t3; + W10 = t1; + W20 = t2; + W30 = t3; +} + +/* machine-generated by xckernel; do not edit. */ +// gga: d e(r_g) / d R_A, e = the unweighted energy density, +// basis class only (points fixed). No quadrature weight. +static inline void xck_gauxc_hess_egrad_gga(double F_rho_A, double F_sigma_A, double vrho, double vsigma, double & de) { + de = F_rho_A*vrho + F_sigma_A*vsigma; +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: per-function nuclear-displacement field rows. +// Sum over the basis functions of one atom to get the perturbed +// fields rho^(A,x), grad rho^(A,x) and tau^(A,x) at this point. +// dchi/ddchi carry the -d/dr sign and the atom mask. +static inline void xck_gauxc_hess_rows_mgga_tau(double U0_u, double U1_u, double U2_u, double U3_u, double dchi_gA_u, double ddchi_gA_u_x, double ddchi_gA_u_y, double ddchi_gA_u_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double & F_rho, double & F_sigma, double & F_tau, double & G_x, double & G_y, double & G_z) { + const double t0 = U0_u*ddchi_gA_u_x + U1_u*dchi_gA_u; + const double t1 = U0_u*ddchi_gA_u_y + U2_u*dchi_gA_u; + const double t2 = U0_u*ddchi_gA_u_z + U3_u*dchi_gA_u; + F_rho = 2*U0_u*dchi_gA_u; + F_sigma = 4*grad_rho_x*t0 + 4*grad_rho_y*t1 + 4*grad_rho_z*t2; + F_tau = U1_u*ddchi_gA_u_x + U2_u*ddchi_gA_u_y + U3_u*ddchi_gA_u_z; + G_x = 2*t0; + G_y = 2*t1; + G_z = 2*t2; +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: non-Pulay pair term, per grid point. +// The fxc quadratic form on the two perturbed field sets, PLUS +// the sigma-curvature term 2 vsigma G^A . G^B -- sigma is +// quadratic in grad rho, so d2 sigma / dA dB survives for A /= B. +// Multiply by the quadrature weight at the call site. +static inline void xck_gauxc_hess_pair_mgga_tau(double F_rho_A, double F_rho_B, double F_sigma_A, double F_sigma_B, double F_tau_A, double F_tau_B, double G_A_x, double G_A_y, double G_A_z, double G_B_x, double G_B_y, double G_B_z, double v2rho2, double v2rhosigma, double v2rhotau, double v2sigma2, double v2sigmatau, double v2tau2, double vsigma, double & h) { + h = F_rho_A*F_rho_B*v2rho2 + F_sigma_A*F_sigma_B*v2sigma2 + F_tau_A*F_tau_B*v2tau2 + v2rhosigma*(F_rho_A*F_sigma_B + F_rho_B*F_sigma_A) + v2rhotau*(F_rho_A*F_tau_B + F_rho_B*F_tau_A) + v2sigmatau*(F_sigma_A*F_tau_B + F_sigma_B*F_tau_A) + 2*vsigma*(G_A_x*G_B_x + G_A_y*G_B_y + G_A_z*G_B_z); +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: Pulay term per basis-function PAIR, with the +// density-matrix factor D_uv DIVIDED OUT -- multiply by it, +// and by the quadrature weight, at the call site. +// Cannot factorise into an outer product; this is the only +// piece that needs the function pair. +static inline void xck_gauxc_hess_pulay_mgga_tau(double dchi_gA_u, double dchi_gB_v, double ddchi_gA_u_x, double ddchi_gA_u_y, double ddchi_gA_u_z, double ddchi_gB_v_x, double ddchi_gB_v_y, double ddchi_gB_v_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double vtau, double w, double & s) { + const double t0 = 4*vsigma; + const double t1 = dchi_gA_u*t0; + const double t2 = dchi_gB_v*t0; + s = w*(2*dchi_gA_u*dchi_gB_v*vrho + ddchi_gA_u_x*ddchi_gB_v_x*vtau + ddchi_gA_u_x*grad_rho_x*t2 + ddchi_gA_u_y*ddchi_gB_v_y*vtau + ddchi_gA_u_y*grad_rho_y*t2 + ddchi_gA_u_z*ddchi_gB_v_z*vtau + ddchi_gA_u_z*grad_rho_z*t2 + ddchi_gB_v_x*grad_rho_x*t1 + ddchi_gB_v_y*grad_rho_y*t1 + ddchi_gB_v_z*grad_rho_z*t1); +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: same-atom (delta_AB) term, per basis function. +// Complete over Cartesian components; multiply by the +// quadrature weight at the call site. +static inline void xck_gauxc_hess_same_mgga_tau(double U0_u, double U1_u, double U2_u, double U3_u, double d2chi_g2_u, double d3chi_g2_u_x, double d3chi_g2_u_y, double d3chi_g2_u_z, double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double vtau, double w, double & s) { + const double t0 = 4*vsigma; + const double t1 = U0_u*t0; + const double t2 = d2chi_g2_u*t0; + s = w*(2*U0_u*d2chi_g2_u*vrho + U1_u*d3chi_g2_u_x*vtau + U1_u*grad_rho_x*t2 + U2_u*d3chi_g2_u_y*vtau + U2_u*grad_rho_y*t2 + U3_u*d3chi_g2_u_z*vtau + U3_u*grad_rho_z*t2 + d3chi_g2_u_x*grad_rho_x*t1 + d3chi_g2_u_y*grad_rho_y*t1 + d3chi_g2_u_z*grad_rho_z*t1); +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: Pulay point-weight matrix W_ab, weight included. +// Slots: 0 = displaced function, 1-3 = its displaced gradient. +// pulay_uv = sum_ab alpha_a(u) W_ab beta_b(v): contract the rows +// as matrix products rather than per function pair. +static inline void xck_gauxc_hess_pulayW_mgga_tau(double grad_rho_x, double grad_rho_y, double grad_rho_z, double vrho, double vsigma, double vtau, double w, double & W00, double & W01, double & W02, double & W03, double & W10, double & W11, double & W20, double & W22, double & W30, double & W33) { + const double t0 = 4*vsigma*w; + const double t1 = grad_rho_x*t0; + const double t2 = grad_rho_y*t0; + const double t3 = grad_rho_z*t0; + const double t4 = vtau*w; + W00 = 2*vrho*w; + W01 = t1; + W02 = t2; + W03 = t3; + W10 = t1; + W11 = t4; + W20 = t2; + W22 = t4; + W30 = t3; + W33 = t4; +} + +/* machine-generated by xckernel; do not edit. */ +// mgga_tau: d e(r_g) / d R_A, e = the unweighted energy density, +// basis class only (points fixed). No quadrature weight. +static inline void xck_gauxc_hess_egrad_mgga_tau(double F_rho_A, double F_sigma_A, double F_tau_A, double vrho, double vsigma, double vtau, double & de) { + de = F_rho_A*vrho + F_sigma_A*vsigma + F_tau_A*vtau; +} + +/* machine-generated by xckernel; do not edit. */ +// Becke/SSF confocal coordinate mu_DE = (|r-R_D| - |r-R_E|)/|R_D-R_E|, +// its gradient and full Hessian in (R_D, R_E) at fixed r. +static inline void xck_gauxc_weight_mu(double R_D_x, double R_D_y, double R_D_z, double R_E_x, double R_E_y, double R_E_z, double r_x, double r_y, double r_z, double & mu, double & dmu_0, double & dmu_1, double & dmu_2, double & dmu_3, double & dmu_4, double & dmu_5, double & d2mu_0_0, double & d2mu_0_1, double & d2mu_0_2, double & d2mu_0_3, double & d2mu_0_4, double & d2mu_0_5, double & d2mu_1_0, double & d2mu_1_1, double & d2mu_1_2, double & d2mu_1_3, double & d2mu_1_4, double & d2mu_1_5, double & d2mu_2_0, double & d2mu_2_1, double & d2mu_2_2, double & d2mu_2_3, double & d2mu_2_4, double & d2mu_2_5, double & d2mu_3_0, double & d2mu_3_1, double & d2mu_3_2, double & d2mu_3_3, double & d2mu_3_4, double & d2mu_3_5, double & d2mu_4_0, double & d2mu_4_1, double & d2mu_4_2, double & d2mu_4_3, double & d2mu_4_4, double & d2mu_4_5, double & d2mu_5_0, double & d2mu_5_1, double & d2mu_5_2, double & d2mu_5_3, double & d2mu_5_4, double & d2mu_5_5) { + const double t0 = R_D_x - R_E_x; + const double t1 = (t0*t0); + const double t2 = R_D_y - R_E_y; + const double t3 = (t2*t2); + const double t4 = R_D_z - R_E_z; + const double t5 = (t4*t4); + const double t6 = t1 + t3 + t5; + const double t7 = pow(t6, -1.0/2.0); + const double t8 = -r_x; + const double t9 = R_D_x + t8; + const double t10 = (t9*t9); + const double t11 = -r_y; + const double t12 = R_D_y + t11; + const double t13 = (t12*t12); + const double t14 = -r_z; + const double t15 = R_D_z + t14; + const double t16 = (t15*t15); + const double t17 = t10 + t13 + t16; + const double t18 = sqrt(t17); + const double t19 = R_E_x + t8; + const double t20 = (t19*t19); + const double t21 = R_E_y + t11; + const double t22 = (t21*t21); + const double t23 = R_E_z + t14; + const double t24 = (t23*t23); + const double t25 = t20 + t22 + t24; + const double t26 = sqrt(t25); + const double t27 = t18 - t26; + const double t28 = 1.0/t18; + const double t29 = 1.0/t6; + const double t30 = t27*t29; + const double t31 = t0*t30; + const double t32 = t2*t30; + const double t33 = t30*t4; + const double t34 = 1.0/t26; + const double t35 = t19*t34; + const double t36 = t21*t34; + const double t37 = t23*t34; + const double t38 = 1.0/t17; + const double t39 = t28*t9; + const double t40 = t0*t39; + const double t41 = 2*t29; + const double t42 = 3*t29; + const double t43 = t30*(t1*t42 - 1); + const double t44 = pow(t17, -3.0/2.0); + const double t45 = t44*t9; + const double t46 = t12*t28; + const double t47 = t0*t46; + const double t48 = t2*t39; + const double t49 = 3*t2; + const double t50 = pow(t6, -2); + const double t51 = t0*t27*t50; + const double t52 = t49*t51; + const double t53 = -t7*(t12*t45 + t29*t47 + t29*t48 - t52); + const double t54 = t15*t28; + const double t55 = t0*t54; + const double t56 = t39*t4; + const double t57 = 3*t4; + const double t58 = t51*t57; + const double t59 = -t7*(t15*t45 + t29*t55 + t29*t56 - t58); + const double t60 = pow(t6, -3.0/2.0); + const double t61 = t0*t35; + const double t62 = 3*t30; + const double t63 = t60*(-t1*t62 + t27 + t40 + t61); + const double t64 = t0*t36; + const double t65 = -t31*t49; + const double t66 = t60*(t48 + t64 + t65); + const double t67 = t0*t37; + const double t68 = -t31*t57; + const double t69 = t60*(t56 + t67 + t68); + const double t70 = t2*t46; + const double t71 = t30*(t3*t42 - 1); + const double t72 = t2*t54; + const double t73 = t4*t46; + const double t74 = t27*t4*t49*t50; + const double t75 = -t7*(t12*t15*t44 + t29*t72 + t29*t73 - t74); + const double t76 = t2*t35; + const double t77 = t60*(t47 + t65 + t76); + const double t78 = t2*t36; + const double t79 = t60*(t27 - t3*t62 + t70 + t78); + const double t80 = t2*t37; + const double t81 = -t32*t57; + const double t82 = t60*(t73 + t80 + t81); + const double t83 = t4*t54; + const double t84 = t30*(t42*t5 - 1); + const double t85 = t35*t4; + const double t86 = t60*(t55 + t68 + t85); + const double t87 = t36*t4; + const double t88 = t60*(t72 + t81 + t87); + const double t89 = t37*t4; + const double t90 = t60*(t27 - t5*t62 + t83 + t89); + const double t91 = 1.0/t25; + const double t92 = pow(t25, -3.0/2.0); + const double t93 = t19*t92; + const double t94 = t7*(t21*t93 - t29*t64 - t29*t76 + t52); + const double t95 = t7*(t23*t93 - t29*t67 - t29*t85 + t58); + const double t96 = t7*(t21*t23*t92 - t29*t80 - t29*t87 + t74); + mu = t27*t7; + dmu_0 = t7*(t28*t9 - t31); + dmu_1 = t7*(t12*t28 - t32); + dmu_2 = t7*(t15*t28 - t33); + dmu_3 = t7*(t31 - t35); + dmu_4 = t7*(t32 - t36); + dmu_5 = t7*(t33 - t37); + d2mu_0_0 = -t7*(t28*(t10*t38 - 1) + t40*t41 - t43); + d2mu_0_1 = t53; + d2mu_0_2 = t59; + d2mu_0_3 = t63; + d2mu_0_4 = t66; + d2mu_0_5 = t69; + d2mu_1_0 = t53; + d2mu_1_1 = -t7*(t28*(t13*t38 - 1) + t41*t70 - t71); + d2mu_1_2 = t75; + d2mu_1_3 = t77; + d2mu_1_4 = t79; + d2mu_1_5 = t82; + d2mu_2_0 = t59; + d2mu_2_1 = t75; + d2mu_2_2 = -t7*(t28*(t16*t38 - 1) + t41*t83 - t84); + d2mu_2_3 = t86; + d2mu_2_4 = t88; + d2mu_2_5 = t90; + d2mu_3_0 = t63; + d2mu_3_1 = t77; + d2mu_3_2 = t86; + d2mu_3_3 = t7*(t34*(t20*t91 - 1) - t41*t61 + t43); + d2mu_3_4 = t94; + d2mu_3_5 = t95; + d2mu_4_0 = t66; + d2mu_4_1 = t79; + d2mu_4_2 = t88; + d2mu_4_3 = t94; + d2mu_4_4 = t7*(t34*(t22*t91 - 1) - t41*t78 + t71); + d2mu_4_5 = t96; + d2mu_5_0 = t69; + d2mu_5_1 = t82; + d2mu_5_2 = t90; + d2mu_5_3 = t95; + d2mu_5_4 = t96; + d2mu_5_5 = t7*(t34*(t24*t91 - 1) - t41*t89 + t84); +} + +/* machine-generated by xckernel; do not edit. */ +// becke cell function s(mu) with t = s'/s and u = s''/s, from +// the factored ln s. Valid strictly inside the switching region +// (|mu| < 1 Becke, |mu| < a SSF); the host handles the outside. +static inline void xck_gauxc_weight_cell_becke(double mu, double & s, double & t, double & u) { + const double t0 = mu - 1; + const double t1 = mu + 2; + const double t2 = (mu*mu*mu); + const double t3 = (mu*mu); + const double t4 = t3 - 3; + const double t5 = t2*(t4*t4*t4); + const double t6 = t3 - 1; + const double t7 = t3*(t4*t4); + const double t8 = -4*t3 + t6*t7 + 4; + const double t9 = 36*mu; + const double t10 = -12*t2 + t5 + t9 + 32; + const double t11 = 3*mu - t2 + 4; + const double t12 = 1.0/t11; + const double t13 = -6*t12*t6 + 4/t1 + 8/t0; + const double t14 = (t6*t6); + const double t15 = -t4; + s = -1.0/256.0*pow(t0, 8)*(t1*t1*t1*t1)*(((3.0/2.0)*mu - 1.0/2.0*t2 + 2)*((3.0/2.0)*mu - 1.0/2.0*t2 + 2))*(-9.0/4.0*mu + (3.0/4.0)*t2 - 1.0/16.0*t5 - 2); + t = t13 + 9*t8/t10; + u = -12*mu*t12 + 18*mu*(3*t14*t4 + t7 - 4)/t10 + ((t13 + 9*(-(t15*t15)*t3*t6 + 4*t3 - 4)/((t15*t15*t15)*t2 + 12*t2 - t9 - 32))*(t13 + 9*(-(t15*t15)*t3*t6 + 4*t3 - 4)/((t15*t15*t15)*t2 + 12*t2 - t9 - 32))) - 18*t14/(t11*t11) - 81*(t8*t8)/(t10*t10) - 4/(t1*t1) - 8/(t0*t0); +} + +/* machine-generated by xckernel; do not edit. */ +// ssf cell function s(mu) with t = s'/s and u = s''/s, from +// the factored ln s. Valid strictly inside the switching region +// (|mu| < 1 Becke, |mu| < a SSF); the host handles the outside. +static inline void xck_gauxc_weight_cell_ssf(double mu, double & s, double & t, double & u) { + const double t0 = (mu*mu); + const double t1 = (mu*mu*mu); + const double t2 = 25*mu - 16; + const double t3 = 16000*mu + 9375*t0 + 7424; + const double t4 = 185600*mu + 200000*t0 + 78125*t1 + 65536; + const double t5 = 1.0/t4; + const double t6 = t3*t5 + 4/t2; + s = (1.0/32.0)*(((25.0/16.0)*mu - 1)*((25.0/16.0)*mu - 1)*((25.0/16.0)*mu - 1)*((25.0/16.0)*mu - 1))*((725.0/16.0)*mu + (3125.0/64.0)*t0 + (78125.0/4096.0)*t1 + 16); + t = 25*t6; + u = -625*(t3*t3)/(t4*t4) + 6250*t5*(75*mu + 64) + 625*(t6*t6) - 2500/(t2*t2); +} +} // namespace xckernel +} // namespace GauXC diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx index 6695d9121..008840366 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.cxx @@ -13,6 +13,7 @@ #include "reference_replicated_xc_host_integrator_exc.hpp" #include "reference_replicated_xc_host_integrator_exc_vxc.hpp" #include "reference_replicated_xc_host_integrator_exc_grad.hpp" +#include "reference_replicated_xc_host_integrator_exc_hess.hpp" #include "reference_replicated_xc_host_integrator_exx.hpp" #include "reference_replicated_xc_host_integrator_fxc_contraction.hpp" #include "reference_replicated_xc_host_integrator_dd_psi.hpp" diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp index a32748eb5..7bc216eb3 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator.hpp @@ -82,6 +82,10 @@ class ReferenceReplicatedXCHostIntegrator : void eval_exc_grad_( int64_t m, int64_t n, const value_type* Ps, int64_t ldps, const value_type* Pz, int64_t lpdz, value_type* EXC_GRAD, const IntegratorSettingsXC& settings ) override; + /// RKS EXC nuclear Hessian + void eval_exc_hess_( int64_t m, int64_t n, const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& settings ) override; + /// sn-LinK void eval_exx_( int64_t m, int64_t n, const value_type* P, int64_t ldp, value_type* K, int64_t ldk, @@ -131,6 +135,10 @@ class ReferenceReplicatedXCHostIntegrator : void exc_grad_local_work_( const value_type* Ps, int64_t ldps, const value_type* Pz, int64_t ldpz, value_type* EXC_GRAD, const IntegratorSettingsXC& ks_settings ); + // Implementation details of exc_hess + void exc_hess_local_work_( const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ); + // Implementation details of sn-LinK void exx_local_work_( const value_type* P, int64_t ldp, value_type* K, int64_t ldk, const IntegratorSettingsEXX& settings ); diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_hess.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_hess.hpp new file mode 100644 index 000000000..84df9cb13 --- /dev/null +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_exc_hess.hpp @@ -0,0 +1,595 @@ +/** + * GauXC Copyright (c) 2020-2024, The Regents of the University of California, + * through Lawrence Berkeley National Laboratory (subject to receipt of + * any required approvals from the U.S. Dept. of Energy). + * + * (c) 2024-2025, Microsoft Corporation + * + * All rights reserved. + * + * See LICENSE.txt for details + */ +#pragma once + +#include "reference_replicated_xc_host_integrator.hpp" +#include "integrator_util/integrator_common.hpp" +#include "host/local_host_work_driver.hpp" +#include "host/blas.hpp" +#include "host/gauxc_hess_kernel.hpp" +#include "host/util.hpp" +#include "common/integrator_constants.hpp" +#include +#include +#include +#include + +namespace GauXC::detail { + +/** + * Log-derivatives of one grid point's partition weight w = q P_C / Z, in + * REDUCED coordinates: every atom but the parent C, the point held fixed. + * The parent's rows and columns follow by translational invariance. + * + * dlw = (d w) / w = a_C - g, g = sum_D pi_D a_D + * d2lw = (d2 w) / w = (a_C-g)(a_C-g)^T + B_C + * - sum_D pi_D [ (a_D-g)(a_D-g)^T + B_D ] + * + * with pi_D = P_D / Z, a_D = d ln P_D = sum_E t_DE dmu_DE and + * B_D = sum_E (u_DE - t_DE^2) dmu_DE dmu_DE^T + t_DE d2mu_DE, where + * t = s'/s and u = s''/s of the cell function. Every term stays finite + * as P_C -> 0. mu_DE with its derivatives, and s, t, u, are GENERATED + * (xckernel gauxcwriter) -- t and u from the factored ln s, which does + * not cancel as mu -> 1. + * + * Returns false when the derivatives vanish identically: an SSF point + * inside the parent's cutoff sphere, or a zero partition. + */ +struct WeightDerivScratch { + std::vector rd, P, S, T, U, a, g; +}; + +inline bool partition_weight_log_derivs( bool is_becke, const Molecule& mol, + const MolMeta& meta, int32_t iC, double dist_nearest, const double* rg, + WeightDerivScratch& w, double* dlw, double* d2lw ) { + + const int32_t N = static_cast(mol.natoms()); + const int32_t n3 = 3*N; + const auto& RAB = meta.rab(); + const double a_ssf = integrator::magic_ssf_factor<>; + w.rd.resize(N); w.P.resize(N); + w.S.resize(size_t(N)*N); w.T.resize(size_t(N)*N); w.U.resize(size_t(N)*N); + w.a.resize(size_t(N)*n3); w.g.resize(n3); + + for( int32_t D = 0; D < N; ++D ) { + const double dx = rg[0]-mol[D].x, dy = rg[1]-mol[D].y, dz = rg[2]-mol[D].z; + w.rd[D] = std::sqrt( dx*dx + dy*dy + dz*dz ); + } + // SSF: inside this sphere the partition is exactly 1, flat in every atom + if( not is_becke and w.rd[iC] < 0.5*(1.-a_ssf)*dist_nearest ) return false; + + for( int32_t D = 0; D < N; ++D ) { + w.P[D] = 1.; + for( int32_t E = 0; E < N; ++E ) { + if( E == D ) continue; + const double mu = (w.rd[D] - w.rd[E]) / RAB[E + D*N]; + double s = 1., t = 0., u = 0.; + if( is_becke ? (1. - mu < 1e-10) : (mu >= a_ssf) ) { + s = 0.; // flat zero (Becke: the point sits on E) + } else if( is_becke or mu > -a_ssf ) { + #include "host/gauxc_hess_call_cell.inc" + } + w.S[size_t(D)*N+E] = s; w.T[size_t(D)*N+E] = t; w.U[size_t(D)*N+E] = u; + w.P[D] *= s; + } + } + double Z = 0.; + for( int32_t D = 0; D < N; ++D ) Z += w.P[D]; + if( Z <= 0. or w.P[iC] <= 0. ) return false; + + std::fill( w.g.begin(), w.g.end(), 0. ); + std::fill( d2lw, d2lw + size_t(n3)*n3, 0. ); + auto active = [&]( int32_t D ) { return D == iC or w.P[D] > 1e-16*Z; }; + + for( int32_t D = 0; D < N; ++D ) { + if( not active(D) ) continue; + double* aD = w.a.data() + size_t(D)*n3; + std::fill( aD, aD + n3, 0. ); + const double piD = w.P[D] / Z; + const double bcoef = (D == iC ? 1. : 0.) - piD; + for( int32_t E = 0; E < N; ++E ) { + if( E == D ) continue; + const double t = w.T[size_t(D)*N+E], u = w.U[size_t(D)*N+E]; + if( t == 0. and u == 0. ) continue; + const double RD[3] = { mol[D].x, mol[D].y, mol[D].z }; + const double RE[3] = { mol[E].x, mol[E].y, mol[E].z }; + double mu, dmu[6], d2mu[36]; + #include "host/gauxc_hess_call_mu.inc" + (void)mu; + int idx[6]; + for( int c = 0; c < 3; ++c ) { + idx[c] = (D == iC) ? -1 : 3*D + c; // the parent is not a free + idx[3+c] = (E == iC) ? -1 : 3*E + c; // coordinate here + } + for( int i = 0; i < 6; ++i ) if( idx[i] >= 0 ) aD[idx[i]] += t * dmu[i]; + if( bcoef != 0. ) + for( int i = 0; i < 6; ++i ) if( idx[i] >= 0 ) + for( int j = 0; j < 6; ++j ) if( idx[j] >= 0 ) + d2lw[ size_t(idx[i])*n3 + idx[j] ] += + bcoef * ( (u - t*t)*dmu[i]*dmu[j] + t*d2mu[6*i+j] ); + } + for( int32_t k = 0; k < n3; ++k ) w.g[k] += piD * aD[k]; + } + + const double* aC = w.a.data() + size_t(iC)*n3; + for( int32_t k = 0; k < n3; ++k ) dlw[k] = aC[k] - w.g[k]; + for( int32_t k = 0; k < n3; ++k ) + for( int32_t l = 0; l < n3; ++l ) d2lw[ size_t(k)*n3 + l ] += dlw[k]*dlw[l]; + for( int32_t D = 0; D < N; ++D ) { + if( not active(D) ) continue; + const double piD = w.P[D] / Z; + const double* aD = w.a.data() + size_t(D)*n3; + for( int32_t k = 0; k < n3; ++k ) { + const double vk = piD * (aD[k] - w.g[k]); + if( vk == 0. ) continue; + for( int32_t l = 0; l < n3; ++l ) + d2lw[ size_t(k)*n3 + l ] -= vk * (aD[l] - w.g[l]); + } + } + return true; +} + +/** + * RKS nuclear Hessian of the XC energy at fixed density matrix. + * + * The XC energy depends on the nuclear coordinates three ways: the basis + * functions ride their atoms, the grid points ride their parent atoms, and + * the partition weights change. IntegratorSettingsEXC_HESS mirrors + * exc_grad: + * + * - include_weight_derivatives = false: the basis class alone -- fixed + * grid, fixed weights, the analogue of exc_grad's Hellmann-Feynman + * option. + * - include_weight_derivatives = true (default): all three. A task's + * points ride its parent C, so its energy depends on R_A - R_C only; + * each task is accumulated in coordinates that exclude C and C's rows + * and columns are restored by translational invariance. The weight + * class adds w'' e + w' e'^T + e' w'^T per point, with the Becke or + * SSF weight derivatives from partition_weight_log_derivs (LKO is not + * implemented). + * + * The per-point kernels are GENERATED (xckernel/emitters/gauxcwriter.py from + * engine/geometric.geometric_hessian), so the contraction cannot drift from + * the expressions that tool validates against finite differences. The + * assembly recipe here -- sum the per-function rows over the shells of one + * atom to get the nuclear-perturbed fields, contract the two row sets as an + * OUTER PRODUCT, then add the Pulay and same-atom seeds -- is itself checked + * against the unfactorised expression by that tool, because every kernel can + * be individually correct while the recipe is a wrong reading of them. + */ +template +void ReferenceReplicatedXCHostIntegrator:: + eval_exc_hess_( int64_t m, int64_t n, const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& ks_settings ) { + + const auto& basis = this->load_balancer_->basis(); + const int64_t nbf = basis.nbf(); + if( m != n ) GAUXC_GENERIC_EXCEPTION("P Must Be Square"); + if( m != nbf ) GAUXC_GENERIC_EXCEPTION("P Must Have Same Dimension as Basis"); + if( ldp < nbf ) GAUXC_GENERIC_EXCEPTION("Invalid LDP"); + + this->load_balancer_->get_tasks(); + + this->timer_.time_op("XCIntegrator.LocalWork", [&](){ + exc_hess_local_work_( P, ldp, EXC_HESS, ks_settings ); + }); + + this->timer_.time_op("XCIntegrator.Allreduce", [&](){ + if( not this->reduction_driver_->takes_host_memory() ) + GAUXC_GENERIC_EXCEPTION("This Module Only Works With Host Reductions"); + const int n3 = 3*static_cast(this->load_balancer_->molecule().natoms()); + this->reduction_driver_->allreduce_inplace( EXC_HESS, n3*n3, ReductionOp::Sum ); + }); +} + + +template +void ReferenceReplicatedXCHostIntegrator:: + exc_hess_local_work_( const value_type* P, int64_t ldp, + value_type* EXC_HESS, const IntegratorSettingsXC& settings ) { + + auto* lwd = dynamic_cast(this->local_work_driver_.get()); + const auto& basis = this->load_balancer_->basis(); + const auto& mol = this->load_balancer_->molecule(); + const auto& func = *this->func_; + + if( func.needs_laplacian() ) + GAUXC_GENERIC_EXCEPTION("EXC Hessian Not Implemented For Laplacian-Dependent MGGAs"); + if( func.is_polarized() ) + GAUXC_GENERIC_EXCEPTION("EXC Hessian Only Implemented For RKS"); + + BasisSetMap basis_map(basis, mol); + const int32_t nbf = basis.nbf(); + const int32_t natoms = static_cast(mol.natoms()); + const int32_t n3 = 3*natoms; + + auto& tasks = this->load_balancer_->get_tasks(); + + for( int32_t i = 0; i < n3*n3; ++i ) EXC_HESS[i] = 0.; + + const bool is_gga = func.is_gga() or func.is_mgga(); + const bool is_mgga = func.is_mgga(); + const int nfield = is_mgga ? 3 : (is_gga ? 2 : 1); // rho[, sigma[, tau]] + + IntegratorSettingsEXC_HESS hess_settings; + if( auto* tmp = dynamic_cast(&settings) ) + hess_settings = *tmp; + const bool full = hess_settings.include_weight_derivatives; + + const auto& molmeta = this->load_balancer_->molmeta(); + const auto& lb_state = this->load_balancer_->state(); + if( full and not lb_state.modified_weights_are_stored ) + GAUXC_GENERIC_EXCEPTION("Weights Have Not Been Modified"); + const XCWeightAlg weight_alg = lb_state.weight_alg; + if( full and weight_alg == XCWeightAlg::LKO ) + GAUXC_GENERIC_EXCEPTION("EXC Hessian Weight Derivatives Not Implemented For LKO"); + const bool is_becke = weight_alg == XCWeightAlg::Becke; + const bool partitioned = weight_alg != XCWeightAlg::NOTPARTITIONED; + + const size_t ntasks = tasks.size(); + #pragma omp parallel + { + XCHostData host_data; + std::vector hess_local(size_t(n3)*n3, 0.); + // this task's contribution, before the translational-invariance fold + std::vector HR(size_t(n3)*n3); + std::vector dlw(n3), d2lw( full ? size_t(n3)*n3 : 0 ); + WeightDerivScratch wscr; + + #pragma omp for schedule(dynamic) + for( size_t iT = 0; iT < ntasks; ++iT ) { + auto& task = tasks[iT]; + std::fill( HR.begin(), HR.end(), 0. ); + + const int32_t npts = static_cast(task.points.size()); + const int32_t nbe = task.bfn_screening.nbe; + const int32_t nshells = static_cast(task.bfn_screening.shell_list.size()); + const auto* points = task.points.data()->data(); + const auto* weights = task.weights.data(); + const int32_t* shell_list = task.bfn_screening.shell_list.data(); + + // ---- scratch ------------------------------------------------------ + // basis + grad(3) + hess(6) + der3(10): the Hessian needs one more + // derivative than the gradient does, because d2 chi / dA dB of the + // GRADIENT ingredient is a third derivative of the basis function. + // Even an LDA needs the basis HESSIAN: the delta_AB term carries a + // double displacement of one function. A GGA needs der3 on top, + // because d2/dAdB of the GRADIENT ingredient is a third derivative. + const int ncomp = is_gga ? 20 : 10; + host_data.basis_eval.resize( size_t(ncomp) * npts * nbe ); + host_data.nbe_scr .resize( size_t(nbe) * nbe ); + host_data.zmat .resize( size_t(4) * npts * nbe ); + host_data.den_scr .resize( size_t(4) * npts ); + host_data.eps .resize( npts ); + host_data.vrho .resize( npts ); + host_data.v2rho2 .resize( npts ); + if( is_gga ) { + host_data.gamma .resize( npts ); + host_data.vgamma .resize( npts ); + host_data.v2rhogamma.resize( npts ); + host_data.v2gamma2 .resize( npts ); + } + if( is_mgga ) { + host_data.tau .resize( npts ); + host_data.vtau .resize( npts ); + host_data.v2rhotau .resize( npts ); + host_data.v2gammatau.resize( npts ); + host_data.v2tau2 .resize( npts ); + } + + auto* basis_eval = host_data.basis_eval.data(); + auto* nbe_scr = host_data.nbe_scr.data(); + auto* xmat = host_data.zmat.data(); + auto* den_eval = host_data.den_scr.data(); + + auto* dbf_x = basis_eval + size_t(1)*npts*nbe; + auto* dbf_y = basis_eval + size_t(2)*npts*nbe; + auto* dbf_z = basis_eval + size_t(3)*npts*nbe; + value_type* d2[6] = {nullptr,nullptr,nullptr,nullptr,nullptr,nullptr}; + value_type* d3[10]; + for( int i = 0; i < 10; ++i ) d3[i] = nullptr; + for( int i = 0; i < 6; ++i ) d2[i] = basis_eval + size_t(4+i)*npts*nbe; + if( is_gga ) + for( int i = 0; i < 10; ++i ) d3[i] = basis_eval + size_t(10+i)*npts*nbe; + + auto [submat_map, foo] = + gen_compressed_submat_map( basis_map, task.bfn_screening.shell_list, nbf, nbf ); + + if( is_gga ) + lwd->eval_collocation_der3( npts, nshells, nbe, points, basis, shell_list, + basis_eval, dbf_x, dbf_y, dbf_z, d2[0],d2[1],d2[2],d2[3],d2[4],d2[5], + d3[0],d3[1],d3[2],d3[3],d3[4],d3[5],d3[6],d3[7],d3[8],d3[9] ); + else + lwd->eval_collocation_hessian( npts, nshells, nbe, points, basis, shell_list, + basis_eval, dbf_x, dbf_y, dbf_z, d2[0],d2[1],d2[2],d2[3],d2[4],d2[5] ); + + // X = 2 P B and its gradient rows; the 4 blocks are contiguous, matching + // the contiguous basis+gradient layout eval_xmat expects. + const double xmat_fac = 2.0; + lwd->eval_xmat( 4*npts, nbf, nbe, submat_map, xmat_fac, P, ldp, basis_eval, nbe, + xmat, nbe, nbe_scr ); + auto* xmat_x = xmat + size_t(1)*npts*nbe; + auto* xmat_y = xmat + size_t(2)*npts*nbe; + auto* xmat_z = xmat + size_t(3)*npts*nbe; + + auto* dden_x = den_eval + size_t(1)*npts; + auto* dden_y = den_eval + size_t(2)*npts; + auto* dden_z = den_eval + size_t(3)*npts; + auto* gamma = host_data.gamma.data(); + auto* tau = host_data.tau.data(); + + if( is_mgga ) + lwd->eval_uvvar_mgga_rks( npts, nbe, basis_eval, dbf_x, dbf_y, dbf_z, nullptr, + xmat, nbe, xmat_x, xmat_y, xmat_z, nbe, den_eval, dden_x, dden_y, dden_z, + gamma, tau, nullptr ); + else if( is_gga ) + lwd->eval_uvvar_gga_rks( npts, nbe, basis_eval, dbf_x, dbf_y, dbf_z, + xmat, nbe, den_eval, dden_x, dden_y, dden_z, gamma ); + else + lwd->eval_uvvar_lda_rks( npts, nbe, basis_eval, xmat, nbe, den_eval ); + + auto* vrho = host_data.vrho.data(); + auto* vgamma = host_data.vgamma.data(); + auto* vtau = host_data.vtau.data(); + auto* v2rho2 = host_data.v2rho2.data(); + auto* v2rhogamma = host_data.v2rhogamma.data(); + auto* v2gamma2 = host_data.v2gamma2.data(); + auto* v2rhotau = host_data.v2rhotau.data(); + auto* v2gammatau = host_data.v2gammatau.data(); + auto* v2tau2 = host_data.v2tau2.data(); + + if( is_mgga ) + func.eval_vxc_fxc( npts, den_eval, gamma, nullptr, tau, vrho, vgamma, + nullptr, vtau, v2rho2, v2rhogamma, nullptr, v2rhotau, v2gamma2, + nullptr, v2gammatau, nullptr, nullptr, v2tau2 ); + else if( is_gga ) + func.eval_vxc_fxc( npts, den_eval, gamma, vrho, vgamma, + v2rho2, v2rhogamma, v2gamma2 ); + else + func.eval_vxc_fxc( npts, den_eval, vrho, v2rho2 ); + + // the energy density itself, for the weight-class terms + auto* eps = host_data.eps.data(); + if( full and partitioned ) { + if( is_mgga ) func.eval_exc( npts, den_eval, gamma, nullptr, tau, eps ); + else if( is_gga ) func.eval_exc( npts, den_eval, gamma, eps ); + else func.eval_exc( npts, den_eval, eps ); + } + + // ---- atoms touched by this task ----------------------------------- + std::vector atoms; // local index -> global atom + std::vector sh_atom(nshells); // shell -> LOCAL atom index + for( int32_t ish = 0; ish < nshells; ++ish ) { + const int iAt = basis_map.shell_to_center( shell_list[ish] ); + auto it = std::find( atoms.begin(), atoms.end(), iAt ); + if( it == atoms.end() ) { sh_atom[ish] = static_cast(atoms.size()); + atoms.push_back(iAt); } + else { sh_atom[ish] = static_cast(it - atoms.begin()); } + } + const int32_t nat_loc = static_cast(atoms.size()); + + // rows[(atom,dir)][field][pt]: the nuclear-perturbed fields. Summing the + // generated per-function rows over the shells of one atom is what turns + // the non-Pulay pair term into a rank update. + const int nrow = nfield + 3; // F_* plus G_x,G_y,G_z + std::vector rows( size_t(nat_loc)*3*nrow*npts, 0. ); + auto ROW = [&]( int a, int d, int k, int ip ) -> double& { + return rows[ (((size_t(a)*3 + d)*nrow + k)*npts) + ip ]; + }; + + size_t bf_off = 0; + for( int32_t ish = 0; ish < nshells; ++ish ) { + const int sh_sz = basis[shell_list[ish]].size(); + const int a = sh_atom[ish]; + for( int ibf = 0, mu = static_cast(bf_off); ibf < sh_sz; ++ibf, ++mu ) { + for( int32_t ip = 0; ip < npts; ++ip ) { + const size_t k = size_t(mu) + size_t(ip)*nbe; + const double U0 = xmat[k]; + const double U1 = is_gga ? xmat_x[k] : 0.0; + const double U2 = is_gga ? xmat_y[k] : 0.0; + const double U3 = is_gga ? xmat_z[k] : 0.0; + const double db[3] = { dbf_x[k], dbf_y[k], dbf_z[k] }; + // d2[] packing is xx,xy,xz,yy,yz,zz + const int h_idx[3][3] = { {0,1,2}, {1,3,4}, {2,4,5} }; + for( int d = 0; d < 3; ++d ) { + // the -d/dr sign of a displaced basis function, folded in + const double dchi = -db[d]; + double ddchi[3] = {0.,0.,0.}; + if( is_gga ) + for( int i = 0; i < 3; ++i ) ddchi[i] = -d2[ h_idx[d][i] ][k]; + double F_rho=0., F_sigma=0., F_tau=0., Gx=0., Gy=0., Gz=0.; + // generated: kernel call AND the row writes, so the row + // layout here and in the pair reads below come from one table + #include "host/gauxc_hess_call_rows.inc" + } + } + } + bf_off += sh_sz; + } + + // ---- outer-product (non-Pulay) pair term -------------------------- + for( int32_t a = 0; a < nat_loc; ++a ) + for( int32_t b = 0; b < nat_loc; ++b ) + for( int dx = 0; dx < 3; ++dx ) + for( int dy = 0; dy < 3; ++dy ) { + double acc = 0.; + for( int32_t ip = 0; ip < npts; ++ip ) { + double h = 0.; + #include "host/gauxc_hess_call_pair.inc" + acc += weights[ip] * h; + } + HR[ size_t(3*atoms[a]+dx)*n3 + (3*atoms[b]+dy) ] += acc; + } + + const int h_idx[3][3] = { {0,1,2}, {1,3,4}, {2,4,5} }; + + // ---- same-atom (delta_AB) term ------------------------------------- + // Both displacements hit the SAME function: no pair, no factorisation, + // but only nbe * 9 * npts work. + { + // d3 packing from eval_collocation_der3: xxx,xxy,xxz,xyy,xyz,xzz, + // yyy,yyz,yzz,zzz + static const int d3_idx[3][3][3] = {{{0,1,2},{1,3,4},{2,4,5}}, + {{1,3,4},{3,6,7},{4,7,8}}, + {{2,4,5},{4,7,8},{5,8,9}}}; + size_t off_u = 0; + for( int32_t ish = 0; ish < nshells; ++ish ) { + const int sh_u = basis[shell_list[ish]].size(); + const int a = sh_atom[ish]; + for( int iu = 0, mu = static_cast(off_u); iu < sh_u; ++iu, ++mu ) + for( int dx = 0; dx < 3; ++dx ) + for( int dy = 0; dy < 3; ++dy ) { + double acc = 0.; + for( int32_t ip = 0; ip < npts; ++ip ) { + const size_t k = size_t(mu) + size_t(ip)*nbe; + // two displacement signs cancel: this carries +d2/dxdy + const double d2c = d2[ h_idx[dx][dy] ][k]; + double d3c[3] = {0.,0.,0.}; + if( is_gga ) + for( int i = 0; i < 3; ++i ) d3c[i] = d3[ d3_idx[dx][dy][i] ][k]; + double sv = 0.; + #include "host/gauxc_hess_call_same.inc" + acc += sv; + } + HR[ size_t(3*atoms[a]+dx)*n3 + (3*atoms[a]+dy) ] += acc; + } + off_u += sh_u; + } + } + + // ---- Pulay term, as matrix products -------------------------------- + // The generated weights W_st(g) make the Pulay term bilinear in the + // displaced rows of the two functions, + // pulay_uv = sum_g sum_st Phi^dx_s(u,g) W_st(g) Phi^dy_t(v,g), + // slot 0 the displaced function, slots 1-3 its displaced gradient. So + // M^{dx,dy} = Phi^dx (W Phi^dy)^T is one GEMM, and the Hessian block + // is its D_uv-weighted sum over the functions of each atom pair. + { + std::vector Ploc( size_t(nbe)*nbe ); + detail::submat_set( nbf, nbf, nbe, nbe, P, ldp, Ploc.data(), nbe, submat_map ); + + std::vector bf_atom( nbe ); + for( int32_t ish = 0, mu = 0; ish < nshells; ++ish ) + for( int i = 0; i < basis[shell_list[ish]].size(); ++i ) + bf_atom[mu++] = atoms[ sh_atom[ish] ]; + + // per-point weights, Wbuf[(4*s+t)*npts + ip] + std::vector Wbuf( size_t(16)*npts, 0. ); + for( int32_t ip = 0; ip < npts; ++ip ) { + double W[16] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + #include "host/gauxc_hess_call_pulayW.inc" + for( int i = 0; i < 16; ++i ) Wbuf[ size_t(i)*npts + ip ] = W[i]; + } + + // displaced rows, the -d/dr sign folded in: nbe x (ns*npts) per direction + const int ns = is_gga ? 4 : 1; + const size_t blk = size_t(npts)*nbe; + std::vector Phi( 3*ns*blk ), Z( ns*blk ), M( size_t(nbe)*nbe ); + auto PHI = [&]( int d, int s ) { return Phi.data() + (size_t(d)*ns + s)*blk; }; + value_type* dbf[3] = { dbf_x, dbf_y, dbf_z }; + for( int d = 0; d < 3; ++d ) + for( size_t k = 0; k < blk; ++k ) { + PHI(d,0)[k] = -dbf[d][k]; + for( int i = 0; i < ns-1; ++i ) PHI(d,i+1)[k] = -d2[ h_idx[d][i] ][k]; + } + + for( int dy = 0; dy < 3; ++dy ) { + // Z_s = sum_t W_st Phi^dy_t, column (point) scaling + std::fill( Z.begin(), Z.end(), 0. ); + for( int s = 0; s < ns; ++s ) + for( int t = 0; t < ns; ++t ) { + const double* w = Wbuf.data() + size_t(4*s+t)*npts; + const double* ph = PHI(dy,t); + double* z = Z.data() + s*blk; + for( int32_t ip = 0; ip < npts; ++ip ) { + if( w[ip] == 0. ) continue; + for( int32_t mu = 0; mu < nbe; ++mu ) + z[ mu + size_t(ip)*nbe ] += w[ip] * ph[ mu + size_t(ip)*nbe ]; + } + } + for( int dx = 0; dx < 3; ++dx ) { + blas::gemm( 'N', 'T', nbe, nbe, ns*npts, 1., PHI(dx,0), nbe, + Z.data(), nbe, 0., M.data(), nbe ); + for( int32_t nu = 0; nu < nbe; ++nu ) + for( int32_t mu = 0; mu < nbe; ++mu ) + HR[ size_t(3*bf_atom[mu]+dx)*n3 + (3*bf_atom[nu]+dy) ] += + xmat_fac * Ploc[ size_t(mu) + size_t(nu)*nbe ] * M[ size_t(mu) + size_t(nu)*nbe ]; + } + } + } + + // ---- weight class ------------------------------------------------ + // d2(w e) = w'' e + w' e'^T + e' w'^T + w e'' + // in reduced coordinates; w e'' is the basis class above. e' is the + // basis-class energy-density gradient, from the same per-atom rows. + if( full and partitioned ) { + const int32_t iC = task.iParent; + std::vector de_loc( size_t(nat_loc)*3*npts ); + for( int32_t a = 0; a < nat_loc; ++a ) + for( int d = 0; d < 3; ++d ) + for( int32_t ip = 0; ip < npts; ++ip ) { + double de = 0.; + #include "host/gauxc_hess_call_egrad.inc" + de_loc[ (size_t(a)*3 + d)*npts + ip ] = de; + } + for( int32_t ip = 0; ip < npts; ++ip ) { + if( weights[ip] == 0. ) continue; + if( not partition_weight_log_derivs( is_becke, mol, molmeta, iC, + task.dist_nearest, points + 3*size_t(ip), wscr, dlw.data(), d2lw.data() ) ) + continue; + const double we = weights[ip] * den_eval[ip] * eps[ip]; + for( size_t k = 0; k < size_t(n3)*n3; ++k ) HR[k] += we * d2lw[k]; + for( int32_t a = 0; a < nat_loc; ++a ) { + if( atoms[a] == iC ) continue; // C's functions ride the grid + for( int d = 0; d < 3; ++d ) { + const double v = weights[ip] * de_loc[ (size_t(a)*3 + d)*npts + ip ]; + const size_t K = 3*size_t(atoms[a]) + d; + for( int32_t J = 0; J < n3; ++J ) { + HR[ size_t(J)*n3 + K ] += dlw[J] * v; + HR[ K*n3 + J ] += v * dlw[J]; + } + } + } + } + } + + // ---- fold into the Hessian ----------------------------------------- + // Basis-only: as is. Full: the task's points ride the parent C, so its + // energy depends on R_A - R_C only. Drop C's rows and columns (C's own + // functions ride the grid) and restore them by translational invariance. + if( not full ) { + for( size_t k = 0; k < size_t(n3)*n3; ++k ) hess_local[k] += HR[k]; + } else { + const int32_t iC = task.iParent; + for( int32_t A = 0; A < natoms; ++A ) if( A != iC ) + for( int dx = 0; dx < 3; ++dx ) + for( int32_t B = 0; B < natoms; ++B ) if( B != iC ) + for( int dy = 0; dy < 3; ++dy ) { + const double v = HR[ size_t(3*A+dx)*n3 + (3*B+dy) ]; + if( v == 0. ) continue; + hess_local[ size_t(3*A +dx)*n3 + (3*B +dy) ] += v; + hess_local[ size_t(3*A +dx)*n3 + (3*iC+dy) ] -= v; + hess_local[ size_t(3*iC+dx)*n3 + (3*B +dy) ] -= v; + hess_local[ size_t(3*iC+dx)*n3 + (3*iC+dy) ] += v; + } + } + + } // tasks + + #pragma omp critical + { + for( int32_t i = 0; i < n3*n3; ++i ) EXC_HESS[i] += hess_local[i]; + } + } // omp parallel +} + +} // namespace GauXC::detail