Fixes issue 225 by explicit symmetrization of density matrices

[susilehtola]

Any asymmetry in the density matrix leaks into the GGA Hessian, this PR fixes the problem.

Closes #225

AI Summary

eval_fxc_contraction mishandles a non-symmetric trial density matrix in
its density-gradient channel. Since a transition density matrix is not
symmetric, every GGA and meta-GGA matrix–vector product picks up a spurious
contribution from the trial matrix's antisymmetric part, which destroys the
Hermiticity of the Kohn–Sham Hessian. LDA is immune because it never reads the
density gradient — which is exactly why SVWN5 looked healthy and PBE did not.

The reference density is symmetric, so eval_exc_vxc and eval_exc_grad are
not affected. The defect is confined to the fxc path, and it is present
identically in the host and CUDA backends.

Root cause

ReferenceLocalHostWorkDriver::eval_uvvar_gga_rks (and its uks / mgga
siblings; CUDA uvvars_gga.hpp is line-for-line the same) builds the density
and its gradient from X = fac · P · B:

den_eval[i]    =      blas::dot( inbe, basis_eval + ioff, 1, X_i, 1 );
const auto dx  = 2. * blas::dot( inbe, dbasis_x_eval + ioff, 1, X_i, 1 );

Write the two channels out in terms of the matrix that was fed in:

channel	what the code computes	what it should be
n	Tr[P χχ]	Tr[P χχ] ✅ — already sees only the symmetric part
∇n	2·Tr[P χ ∂χ]	Tr[(P + Pᵀ) χ ∂χ]

The two agree only when P = Pᵀ. That holds for the ground-state density,
so the ground-state code paths are correct and have always been correct. It
does not hold for a trial density, and eval_fxc_contraction reuses the very
same routine on one.

Consequently the kernel returns a non-zero result for a direction along which
Exc does not vary at all:

                          F[(M − Mᵀ)/2]     operator asymmetry
                          (must be 0)       Tr(X·F[Y]) − Tr(Y·F[X]), rel.
  SVWN5  (LDA)             8.85e-18  ✅        9.7e-16   ✅
  PBE0   (GGA)             1.31e-01  ❌        1.79e-01  ❌
  LAK    (meta-GGA)        2.91e+00  ❌        1.872     ❌

Benzene / cc-pVDZ / UltraFine, two random non-symmetric trial matrices. Every
rung above LDA is affected. The meta-GGA case is not marginal: an asymmetry of
1.872 relative means Tr(X·F[Y]) and Tr(Y·F[X]) came out with opposite
signs (−15.08 against the correct +7.36).

This is a structural defect, not a numerical one — an antisymmetric input must
map to exactly zero whatever the functional is — so it is worth demonstrating on
a well-behaved meta-GGA rather than on one with difficult numerics. LAK is used
here for that reason; it is reachable through ExchCXX's libxc-name constructor,
XCKernel(libxc_name_string("MGGA_X_LAK"), spin), without an ExchCXX enum entry.

(The tau channel is not implicated: τ = ½Tr[P ∇χ·∇χ] is a symmetric
bilinear form in the basis pair, so like the density it annihilates an
antisymmetric input by itself. Meta-GGAs are affected because they carry a
gradient channel, not because of τ.)

Why the reproducer triggers it

demonstrate_nonhermitian_hessian.cpp passes a symmetric ground-state density
d — max|d − dᵀ| = 0 — but the ten trial matrices are AO transition
density matrices, and those are not symmetric. Each is rank 1, i.e. an outer
product c_occ ⊗ c_virt, whose antisymmetric part is as large as its
symmetric part:

  tdm      |sym part|   |antisym part|   rank
   [0]      3.914e-01      3.914e-01      1
   [2]      5.592e-01      5.241e-01      1
   [5]      5.876e-01      4.604e-01      1
   ...      (all ten: rank 1, antisym ~ sym)

A perturbed (response) density matrix P¹ = T + Tᵀ would be symmetric by
construction and would not have exposed this. A one-sided transition density
matrix does.

Fix

Symmetrize the trial density once, at both integrator entry points:

• src/xc_integrator/replicated/host/reference_replicated_xc_host_integrator_fxc_contraction.hpp
• src/xc_integrator/replicated/device/incore_replicated_xc_device_integrator_fxc_contraction.hpp

out[i + j*nbf] = 0.5 * ( A[i + j*lda] + A[j + i*lda] );

One nbf² copy against a grid loop, one place per backend, covering every rung
and both spin cases. Doing it here rather than inside the uvvar kernels keeps
a single site per backend and leaves the ground-state paths untouched.

The fix cannot change any LDA result. Symmetrizing preserves the density
channel identically — Tr[M χχ] = Tr[Mˢʸᵐ χχ] — and repairs the gradient
channel, since Tr[(M + Mᵀ) χ ∂χ] = 2·Tr[Mˢʸᵐ χ ∂χ] is precisely what the
existing code computes from Mˢʸᵐ. Verified: the SVWN5 numbers are bit-identical
before and after.

Verification

1. The reporter's own demonstrator, built against the tree with and without
the fix (git revert --no-commit):

	SVWN5	PBE
before	0 asymmetries	12 asymmetries, largest 9.19e-04
after	0	0

Asymmetry found at (2,5) vs (5,2):        ← before
  <v_2, A v_5> = 0.0830374
  <v_5, A v_2> = 0.0821179
  Difference   = 0.000919462

2. Operator symmetry and the antisymmetric null direction, every rung,
before and after:

	F[(M−Mᵀ)/2], must be 0		Tr(X·F[Y]) − Tr(Y·F[X]), rel.
	before	after	before	after
SVWN5 (LDA)	8.85e-18	0.0	9.7e-16	1.5e-15
PBE0 (GGA)	1.31e-01	0.0	1.79e-01	3.9e-15
LAK (meta-GGA)	2.91e+00	0.0	1.872	0.0
SCAN (meta-GGA)	1.28e+00	0.0	2.10e-01	1.0e-14

and F[M] vs F[(M+Mᵀ)/2] for PBE0: 9.7e-02 → 1.6e-16.

3. Against a reference GauXC shares no algebra with. The shipped
tests/2nd_derivative_test.cxx compares FXC to stored values that GauXC
generated itself, so it proves self-consistency only and passed throughout.
Comparing instead against a central difference of eval_exc_vxc — which is what
the contraction is defined to be — for a symmetric trial matrix:

  SVWN5   max|F − FD| / scale = 1.00e-07,  best-fit F = 1.0000000 · FD
  PBE0    max|F − FD| / scale = 1.03e-07,  best-fit F = 1.0000000 · FD

unchanged before and after the fix (the residual is finite-difference
truncation).

4. GauXC's own suite, host build: all 3,640,979 assertions in 56 test
cases pass.

Note for callers

After the fix, eval_fxc_contraction(P, M) returns the kernel applied to
(M + Mᵀ)/2 — the unique answer consistent with Exc depending only on the
symmetric part of a density matrix. A TDDFT driver that wants the response to
the perturbed density P¹ = T + Tᵀ must therefore pass T + Tᵀ, not the
one-sided T; passing T yields half of that. This contract is worth stating
in the API docs, since nothing enforces it.

Suggested regression guard

The existing test cannot see this class of error. Three lines that can:

// Exc does not vary along an antisymmetric direction, so the kernel
// contracted with one must vanish identically -- at every rung.
matrix_type A = 0.5*(tP - tP.transpose());
auto FA = integrator.eval_fxc_contraction(P, A);
CHECK( FA.cwiseAbs().maxCoeff() < 1e-12 );

Out of scope but noticed

• The CUDA half of the fix is a textual mirror of the host change; it has not
  been compiled here (no GPU available).
• The HIP backend cannot link at all, independently of this issue.
  scheme1_base.cxx is compiled for any device build and calls
  eval_tmat_{lda,gga,mgga}, eval_vvars_*_trial,
  eval_kern_vxc_fxc_{lda,gga,mgga}, eval_kern_exc_vxc_mgga and
  exx_ek_screening_bfn_stats, none of which HIP defines. hip/hipify.sh
  converts a hard-coded file list that never gained zmat_fxc.cu,
  uvvars_{lda,gga,mgga}.hpp, increment_exc_grad.cu or
  exx_ek_screening_bfn_stats.cu; hip/kernels/uvvars.hip still carries
  pre-integrator_ks_scheme signatures, so HIP predates UKS/GKS, meta-GGA and
  the fxc contraction entirely. Worth its own issue.

[Copilot]

🟡 Changes recommended

The new API semantics need public documentation and targeted regression coverage.

Get a fresh assessment by requesting another Copilot review.

Pull request overview

Projects trial density matrices onto their symmetric components to preserve Hermiticity in GGA/meta-GGA FXC contractions.

Changes:

• Symmetrizes RKS/UKS trial densities before contraction.
• Applies equivalent handling across host and device backends.

File summaries

File	Description
reference_replicated_xc_host_integrator_fxc_contraction.hpp	Adds host-side trial-density symmetrization.
incore_replicated_xc_device_integrator_fxc_contraction.hpp	Mirrors symmetrization for device execution.

Review details

• Files reviewed: 2/2 changed files
• Comments generated: 2
• Review effort level: Balanced

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[susilehtola]

Addressed in 9e17cd2.

Worth noting why the existing test was blind to this: every stored trial density in the reference data is exactly symmetric (max|tP - tPᵀ| = 0.0), so antisymmetrizing the stored matrices would also have given a vacuous test.

Instead, test_fxc_contractioin — the shared helper, so it runs for Host and every Device section — now builds a non-symmetric M whose symmetric part is bitwise the stored trial density (strict upper triangle doubled, lower zeroed), and checks, for RKS and both UKS channels:

• F[M] reproduces the stored FXC reference, i.e. F[M] == F[(M+Mᵀ)/2] against independent data;
• F[(M−Mᵀ)/2] vanishes.

Verified on a host build that it discriminates:

	new checks	pre-existing checks
without the fix	12 fail — PBE0, SCAN (RKS), BLYP, SCAN (UKS), all 4e-4…3e-3 vs. 1e-10 tolerance	pass
with the fix	pass	pass

LDA sections pass either way, as expected. 44/44 assertions with the fix. Device path not run (no GPU here).

[susilehtola]

Addressed in 9e17cd2. Both XCIntegrator::eval_fxc_contraction overloads (C++) and gauxc_integrator_eval_fxc_contraction_{rks,uks} (C) now document that the kernel is contracted with the symmetric part (tP + tPᵀ)/2 of each trial density, that an antisymmetric part contributes exactly zero, and that a linear-response caller must pass P1 = T + Tᵀ — a one-sided transition density T yields half that response. The C @param entries for the trial densities, which were copy-pasted from the reference density, now say "trial density" too.