When attempting to integrate GauXC into Psi4, I tested using the KS electronic Hessian's matrix-vector products delivered by GauXC's eval_fxc_contraction to drive linear-response TDDFT. While I saw results in excellent agreement with those using Psi4-native matrix-vector products for SVWN5 with a STO-6G water molecule, I observed that when the functional was switched to PBE, the Davidson algorithm had a large residual, preventing convergence with default settings. Further testing revealed a much more specific problem: the PBE KS Hessian fails to even be Hermitian. This is mathematically impossible by equality of mixed partials and prevents the use of standard TDDFT algorithms.
I've attached a reproducer file, which is my Psi4 example translated to use only GauXC and ExchCXX. Unfortunately, there is a large amount of boilerplate. The key lines begin at l. 122. The current implementation uses the SVWN5 functional and computes the hessian in the basis of the 10 transition density matrix. The hessian is hermitian to 1e-10, as expected. But when the functional is switched to PBE, we don't even have hermiticity to within 1e-3.
[ -0.133002 6.72205e-18 -3.22212e-18 3.13957e-18 -5.8534e-19 3.15194e-19 2.43347e-17 7.26674e-19 2.25932e-18 2.35862e-18]
[-8.67362e-18 -0.108432 1.22851e-18 -4.98166e-18 9.48157e-18 3.64521e-18 -9.60879e-19 -8.43604e-18 2.26738e-18 9.63097e-18]
[-1.90947e-18 -3.07472e-18 -0.159646 7.71841e-13 -8.44221e-13 0.0830374 0.0316462 -7.24844e-13 -0.0598633 5.47641e-15]
[-4.53788e-18 2.71539e-18 8.86596e-13 -0.160716 0.0951338 -9.15934e-16 3.81639e-17 0.000560716 2.83866e-15 -0.0636568]
[-1.54859e-18 8.16588e-19 -1.01986e-12 0.0899572 -0.183399 6.66134e-16 -3.1225e-16 -0.133125 -3.98347e-15 -0.0577205]
[ 1.00931e-17 4.54297e-18 0.0821179 -6.80012e-16 6.52256e-16 -0.280076 -0.117921 -8.60423e-16 -0.0435206 -4.63908e-15]
[ 7.40146e-19 1.21981e-18 0.0235284 4.19803e-16 -6.73073e-16 -0.114792 -0.274316 -8.1532e-16 -0.250539 -1.90554e-15]
[-2.28289e-19 1.55441e-18 -7.7645e-13 -0.00116171 -0.126807 -1.17267e-15 -8.25728e-16 -0.138134 -1.91137e-15 -0.0953222]
[ 1.63569e-18 7.59115e-18 -0.0689541 3.66742e-15 -5.24017e-15 -0.0423594 -0.2421 -2.50922e-15 -0.349724 -9.54589e-17]
[-1.54643e-18 1.73319e-18 1.41269e-14 -0.0623118 -0.0462046 -6.02936e-15 -2.53595e-15 -0.0935492 -9.6743e-17 -0.139341]
Asymmetry found at (2,5) vs (5,2):
<v_2, A v_5> = 0.0830374
<v_5, A v_2> = 0.0821179
Difference = 0.000919462
Asymmetry found at (2,6) vs (6,2):
<v_2, A v_6> = 0.0316462
<v_6, A v_2> = 0.0235284
Difference = 0.00811782
Asymmetry found at (2,8) vs (8,2):
<v_2, A v_8> = -0.0598633
<v_8, A v_2> = -0.0689541
Difference = 0.00909089
Asymmetry found at (3,4) vs (4,3):
<v_3, A v_4> = 0.0951338
<v_4, A v_3> = 0.0899572
Difference = 0.00517658
Asymmetry found at (3,7) vs (7,3):
<v_3, A v_7> = 0.000560716
<v_7, A v_3> = -0.00116171
Difference = 0.00172242
Asymmetry found at (3,9) vs (9,3):
<v_3, A v_9> = -0.0636568
<v_9, A v_3> = -0.0623118
Difference = 0.00134503
Asymmetry found at (4,7) vs (7,4):
<v_4, A v_7> = -0.133125
<v_7, A v_4> = -0.126807
Difference = 0.00631806
Asymmetry found at (4,9) vs (9,4):
<v_4, A v_9> = -0.0577205
<v_9, A v_4> = -0.0462046
Difference = 0.0115158
Asymmetry found at (5,6) vs (6,5):
<v_5, A v_6> = -0.117921
<v_6, A v_5> = -0.114792
Difference = 0.00312888
Asymmetry found at (5,8) vs (8,5):
<v_5, A v_8> = -0.0435206
<v_8, A v_5> = -0.0423594
Difference = 0.00116119
Asymmetry found at (6,8) vs (8,6):
<v_6, A v_8> = -0.250539
<v_8, A v_6> = -0.2421
Difference = 0.0084395
Asymmetry found at (7,9) vs (9,7):
<v_7, A v_9> = -0.0953222
<v_9, A v_7> = -0.0935492
Difference = 0.00177304
When attempting to integrate GauXC into Psi4, I tested using the KS electronic Hessian's matrix-vector products delivered by GauXC's
eval_fxc_contractionto drive linear-response TDDFT. While I saw results in excellent agreement with those using Psi4-native matrix-vector products for SVWN5 with a STO-6G water molecule, I observed that when the functional was switched to PBE, the Davidson algorithm had a large residual, preventing convergence with default settings. Further testing revealed a much more specific problem: the PBE KS Hessian fails to even be Hermitian. This is mathematically impossible by equality of mixed partials and prevents the use of standard TDDFT algorithms.I've attached a reproducer file, which is my Psi4 example translated to use only GauXC and ExchCXX. Unfortunately, there is a large amount of boilerplate. The key lines begin at l. 122. The current implementation uses the SVWN5 functional and computes the hessian in the basis of the 10 transition density matrix. The hessian is hermitian to 1e-10, as expected. But when the functional is switched to PBE, we don't even have hermiticity to within 1e-3.
demonstrate_nonhermitian_hessian.cpp