Fixes issue 225 by explicit symmetrization of density matrices

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,

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&,

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 <gauxc/util/misc.hpp>
 #include <gauxc/util/unused.hpp>
+#include <vector>
 
 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<value_type> tPs_sym, tPz_sym;
+    auto symmetrize = [nbf]( const value_type* A, int64_t lda,
+                             std::vector<value_type>& 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();
 

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 <stdexcept>
+#include <vector>
 
 namespace GauXC::detail {
 
@@ -57,6 +58,43 @@ void ReferenceReplicatedXCHostIntegrator<ValueType>::
     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<value_type> tPs_sym, tPz_sym;
+  auto symmetrize = [nbf]( const value_type* A, int64_t lda,
+                           std::vector<value_type>& 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<ValueType>::
   // 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() );

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<Eigen::StrictlyUpper>().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) {