From 6f42313a331929db3fd8454f913d12cc357dc906 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Tue, 15 Sep 2026 10:20:36 +0300 Subject: [PATCH 1/2] Fixes issue 225 by explicit symmetrization of density matrices --- ...d_xc_device_integrator_fxc_contraction.hpp | 21 ++++++++++ ...ted_xc_host_integrator_fxc_contraction.hpp | 40 ++++++++++++++++++- 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp index ffc0ca417..f1c1cd70d 100644 --- a/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp +++ b/src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp @@ -13,6 +13,7 @@ #include "incore_replicated_xc_device_integrator.hpp" #include #include +#include namespace GauXC::detail { @@ -66,6 +67,26 @@ namespace GauXC::detail { GAUXC_GENERIC_EXCEPTION("Invalid LDFXCz"); } + // Symmetrize the TRIAL densities -- see the host integrator for why: + // Exc sees only the symmetric part of a density matrix, but the + // gradient channel of the uvvar kernels assumes symmetry rather than + // enforcing it, so an antisymmetric trial direction (a TDDFT + // B-matrix or stability block) produced a nonzero kernel and broke + // the Hermiticity of the Kohn-Sham Hessian for every GGA and + // meta-GGA. https://github.com/wavefunction91/GauXC/issues/225 + std::vector tPs_sym, tPz_sym; + auto symmetrize = [nbf]( const value_type* A, int64_t lda, + std::vector& out ) -> const value_type* { + if( !A ) return nullptr; + out.resize( nbf*nbf ); + for( int64_t j = 0; j < nbf; ++j ) + for( int64_t i = 0; i < nbf; ++i ) + out[i + j*nbf] = 0.5 * ( A[i + j*lda] + A[j + i*lda] ); + return out.data(); + }; + tPs = symmetrize( tPs, ldtps, tPs_sym ); if( tPs ) ldtps = nbf; + tPz = symmetrize( tPz, ldtpz, tPz_sym ); if( tPz ) ldtpz = nbf; + // Get Tasks auto& tasks = this->load_balancer_->get_tasks(); diff --git a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp index 7320039c7..c338acce2 100644 --- a/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp +++ b/src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp @@ -16,6 +16,7 @@ #include "host/local_host_work_driver.hpp" #include "host/blas.hpp" #include +#include namespace GauXC::detail { @@ -57,6 +58,43 @@ void ReferenceReplicatedXCHostIntegrator:: GAUXC_GENERIC_EXCEPTION("Invalid LDFXCZ"); + // Symmetrize the TRIAL densities. + // + // The XC energy is a functional of n = Tr[P B B], which sees only the + // symmetric part of P: the derivative of Exc along an antisymmetric + // direction is identically zero. The density channel gets that for + // free -- Tr[P B B] kills an antisymmetric P by itself -- but the + // GRADIENT channel does not. eval_uvvar_* builds it as + // d_x n = 2 X . d_x B, X = fac * P * B, + // which equals Tr[(P + P^T) B d_x B] only when P is symmetric. That + // holds for the reference density, so exc_vxc and exc_grad are + // unaffected; it does NOT hold for a trial density, which a caller has + // no reason to hand over symmetric -- a TDDFT B-matrix or stability + // block is contracted against the ANTISYMMETRIC combination on + // purpose. The result was an XC kernel that answered nonzero to a + // direction along which Exc does not vary, breaking the Hermiticity + // of the Kohn-Sham Hessian for every GGA and meta-GGA while leaving + // LDA (which never touches the gradient channel) looking correct. + // See https://github.com/wavefunction91/GauXC/issues/225. + // + // Symmetrizing here rather than in the uvvar kernels keeps the fix in + // one place for every backend and every rung, and costs one nbf^2 + // copy against a grid loop. + std::vector tPs_sym, tPz_sym; + auto symmetrize = [nbf]( const value_type* A, int64_t lda, + std::vector& out ) -> const value_type* { + if( !A ) return nullptr; + out.resize( nbf*nbf ); + for( int64_t j = 0; j < nbf; ++j ) + for( int64_t i = 0; i < nbf; ++i ) + out[i + j*nbf] = 0.5 * ( A[i + j*lda] + A[j + i*lda] ); + return out.data(); + }; + const value_type* tPs_use = symmetrize( tPs, ldtps, tPs_sym ); + const value_type* tPz_use = symmetrize( tPz, ldtpz, tPz_sym ); + const int64_t ldtps_use = tPs_use ? nbf : ldtps; + const int64_t ldtpz_use = tPz_use ? nbf : ldtpz; + // Get Tasks auto& tasks = this->load_balancer_->get_tasks(); @@ -66,7 +104,7 @@ void ReferenceReplicatedXCHostIntegrator:: // Compute Local contributions to FXC contraction this->timer_.time_op("XCIntegrator.LocalWork", [&](){ fxc_contraction_local_work_( basis, Ps, ldps, Pz, ldpz, - tPs, ldtps, tPz, ldtpz, + tPs_use, ldtps_use, tPz_use, ldtpz_use, FXCs, ldfxcs, FXCz, ldfxcz, &N_EL, ks_settings, tasks.begin(), tasks.end() ); From 9e17cd21305d56d89a3b1158cb97a10f68bde9c8 Mon Sep 17 00:00:00 2001 From: Susi Lehtola Date: Wed, 16 Sep 2026 11:45:03 +0300 Subject: [PATCH 2/2] Test and document symmetrization of FXC trial densities Address review feedback on #226. Tests: the stored trial densities in the reference data are exactly symmetric, so the existing FXC reference comparison cannot detect how a non-symmetric trial density is handled -- which is why it passed throughout #225. For every functional, spin case and backend exercised by test_fxc_contractioin, build a non-symmetric trial matrix M whose symmetric part is bitwise the stored trial density (strict upper triangle doubled, lower zeroed), and check that * F[M] reproduces the stored FXC reference, and * F[(M - M^T)/2] vanishes. Without the symmetrization fix, all 12 new checks for PBE0, SCAN (RKS) and BLYP, SCAN (UKS) fail by 4e-4..3e-3 against a 1e-10 tolerance, while the LDA sections and all pre-existing checks still pass. With the fix, all 44 assertions pass (host build). Docs: state at the public C++ and C API entry points that the kernel is contracted with the symmetric part of the trial density, and that a linear-response caller must pass P1 = T + T^T; a one-sided transition density yields half that response. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01Sk6Km8XjNxNU8wyRDQ8sxx --- include/gauxc/c/xc_integrator.h | 22 ++++++++++----- include/gauxc/xc_integrator.hpp | 14 ++++++++++ tests/2nd_derivative_test.cxx | 47 +++++++++++++++++++++++++++++++++ 3 files changed, 77 insertions(+), 6 deletions(-) diff --git a/include/gauxc/c/xc_integrator.h b/include/gauxc/c/xc_integrator.h index 76a0095f8..ec15d90da 100644 --- a/include/gauxc/c/xc_integrator.h +++ b/include/gauxc/c/xc_integrator.h @@ -343,10 +343,15 @@ extern void gauxc_integrator_eval_exx_rks( * @param n Number of columns in the density matrix. * @param density_matrix Pointer to the density matrix data. * @param ldp Leading dimension of the density matrix. - * @param t_density_matrix Pointer to the density matrix data. - * @param ldtp Leading dimension of the density matrix. + * @param t_density_matrix Pointer to the trial density matrix data. + * @param ldtp Leading dimension of the trial density matrix. * @param fxc Matrix container to store the FXC contraction result. * @param ldfxc Leading dimension of the FXC contraction matrix. + * @note Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density tP; any antisymmetric part contributes exactly zero. + * For linear response, pass the perturbed density P1 = T + T^T: + * passing a one-sided transition density T yields half that response. */ extern void gauxc_integrator_eval_fxc_contraction_rks( GauXCStatus* status, @@ -371,14 +376,19 @@ extern void gauxc_integrator_eval_fxc_contraction_rks( * @param ldp_s Leading dimension of the total density matrix. * @param density_matrix_z Pointer to the density matrix data for spin density. * @param ldp_z Leading dimension of the spin density matrix. - * @param t_density_matrix_s Pointer to the density matrix data for total density. - * @param ldtp_s Leading dimension of the total density matrix. - * @param t_density_matrix_z Pointer to the density matrix data for spin density. - * @param ldtp_z Leading dimension of the spin density matrix. + * @param t_density_matrix_s Pointer to the trial density matrix data for total density. + * @param ldtp_s Leading dimension of the total trial density matrix. + * @param t_density_matrix_z Pointer to the trial density matrix data for spin density. + * @param ldtp_z Leading dimension of the spin trial density matrix. * @param fxc_s Matrix container to store the FXC contraction result for total density. * @param ldfxc_s Leading dimension of the FXC contraction matrix for total density. * @param fxc_z Matrix container to store the FXC contraction result for spin density. * @param ldfxc_z Leading dimension of the FXC contraction matrix for spin density. + * @note Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density tP; any antisymmetric part contributes exactly zero. + * For linear response, pass the perturbed density P1 = T + T^T: + * passing a one-sided transition density T yields half that response. */ extern void gauxc_integrator_eval_fxc_contraction_uks( GauXCStatus* status, diff --git a/include/gauxc/xc_integrator.hpp b/include/gauxc/xc_integrator.hpp index 03feaf934..8dd8eed65 100644 --- a/include/gauxc/xc_integrator.hpp +++ b/include/gauxc/xc_integrator.hpp @@ -80,6 +80,20 @@ class XCIntegrator { exx_type eval_exx ( const MatrixType&, const IntegratorSettingsEXX& = IntegratorSettingsEXX{} ); + /** + * @brief Contract the XC kernel with a trial density matrix. + * + * RKS: eval_fxc_contraction( P, tP ). + * UKS: eval_fxc_contraction( Ps, Pz, tPs, tPz ), returning ( FXCs, FXCz ). + * + * Exc depends only on the symmetric part of a density matrix, so the + * kernel is contracted with the symmetric part (tP + tP^T)/2 of each + * trial density; any antisymmetric part contributes exactly zero. + * + * In particular, for linear response the trial density should be the + * perturbed density P1 = T + T^T. Passing a one-sided transition density + * T instead yields half of that response. + */ fxc_contraction_type_rks eval_fxc_contraction ( const MatrixType&, const MatrixType&, const IntegratorSettingsXC& = IntegratorSettingsXC{} ); fxc_contraction_type_uks eval_fxc_contraction ( const MatrixType&, const MatrixType&, const MatrixType&, const MatrixType&, diff --git a/tests/2nd_derivative_test.cxx b/tests/2nd_derivative_test.cxx index eedcb27f7..f7bd8ddcc 100644 --- a/tests/2nd_derivative_test.cxx +++ b/tests/2nd_derivative_test.cxx @@ -134,6 +134,53 @@ void test_fxc_contractioin(ExecutionSpace ex, const RuntimeEnvironment& rt, CHECK(FXCz_diff_nrm / basis.nbf() < 1e-10); } + + // The stored trial densities are symmetric, so the comparison above + // cannot see how a non-symmetric trial density (e.g. a one-sided + // transition density) is handled. Exc depends only on the symmetric + // part of a density matrix, so the contraction must act on (M + M^T)/2 + // and vanish identically along (M - M^T)/2, at every rung. The gradient + // channel of GGAs and meta-GGAs used to violate both, see + // https://github.com/wavefunction91/GauXC/issues/225. + // + // M is built so that (M + M^T)/2 equals the stored trial density exactly + // (strict upper triangle doubled, lower triangle zero), so F[M] can be + // checked against the same reference data. + auto make_nonsymmetric = []( const matrix_type& A ) -> matrix_type { + matrix_type M = 2.0 * A.triangularView().toDenseMatrix(); + M.diagonal() = A.diagonal(); + return M; + }; + auto antisymmetric_part = []( const matrix_type& M ) -> matrix_type { + return 0.5 * ( M - M.transpose() ); + }; + + if (rks) { + matrix_type tM = make_nonsymmetric(tP); + matrix_type tA = antisymmetric_part(tM); + REQUIRE( tA.norm() > 0.1 * tP.norm() ); + + auto FXC_M = integrator.eval_fxc_contraction(P, tM); + CHECK((FXC_M - FXC_ref).norm() / basis.nbf() < 1e-10); + + auto FXC_A = integrator.eval_fxc_contraction(P, tA); + CHECK(FXC_A.norm() / basis.nbf() < 1e-10); + } else if (uks) { + matrix_type tMs = make_nonsymmetric(tP); + matrix_type tMz = make_nonsymmetric(tPz); + matrix_type tAs = antisymmetric_part(tMs); + matrix_type tAz = antisymmetric_part(tMz); + REQUIRE( tAs.norm() > 0.1 * tP.norm() ); + REQUIRE( tAz.norm() > 0.1 * tPz.norm() ); + + auto [FXCs_M, FXCz_M] = integrator.eval_fxc_contraction(P, Pz, tMs, tMz); + CHECK((FXCs_M - FXC_ref ).norm() / basis.nbf() < 1e-10); + CHECK((FXCz_M - FXCz_ref).norm() / basis.nbf() < 1e-10); + + auto [FXCs_A, FXCz_A] = integrator.eval_fxc_contraction(P, Pz, tAs, tAz); + CHECK(FXCs_A.norm() / basis.nbf() < 1e-10); + CHECK(FXCz_A.norm() / basis.nbf() < 1e-10); + } } void test_integrator_2nd(std::string reference_file, functional_type& func, PruningScheme pruning_scheme) {