The Cartesian basis DF-J gradient disagrees from the PySCF on CPU

[vmitq]

For a Cartesian AO/auxiliary basis, the default gpu4pyscf DF-RHF nuclear gradient differs from PySCF’s CPU-based one, evaluated on the same converged wavefunction. In the case of H2O/cc-pvdz evaluated with cart=True the difference is 2.4e-7 a.u., while for the spherical basis it is 4.5e-12 a.u. For larger molecules like taxol, the discrepancy increases significantly, exceeding 1e-4. A detailed analysis (computing gradient contributions from integrals and finite differences checks) shows that the source of the error is the Coulomb 2e gradient contribution in gpu4pyscf.

Environment

• gpu4pyscf-cuda12x 1.7.1 (latest PyPI release, 2026-05-28), pyscf 2.13.1
• CUDA 12.x, NVIDIA A100 80GB
• Basis cc-pvdz, auxiliary cc-pvdz-jkfit

Minimal example

import numpy as np
import cupy
from pyscf import gto, scf

H2O = "O 0.000 0.000 0.117; H 0.000 0.757 -0.467; H 0.000 -0.757 -0.467"

for cart in (True, False):
    mol = gto.M(atom=H2O, unit="angstrom", basis="cc-pvdz", cart=cart, verbose=0)
    mf = scf.RHF(mol).density_fit(auxbasis="cc-pvdz-jkfit").run()
    de_cpu = mf.nuc_grad_method().kernel()
    de_gpu = cupy.asnumpy(mf.to_gpu().nuc_grad_method().kernel())
    print(f"{'CARTESIAN' if cart else 'SPHERICAL'}: "
          f"max|pyscf - gpu4pyscf| = {np.abs(de_cpu - de_gpu).max():.2e}")

Output:

CARTESIAN: max|pyscf - gpu4pyscf| = 2.42e-07
SPHERICAL: max|pyscf - gpu4pyscf| = 4.79e-12

[sunqm]

This difference is due to the different DF metric solver used for Cartesian case. gpu4pyscf uses an eigen-decomposition. pyscf uses np.linalg.solve. the relevant code is here

gpu4pyscf/gpu4pyscf/df/grad/rhf.py

Line 124 in 3ae9805

metric = aux_coeff.dot(_gen_metric_solver(j2c, 'ED')(aux_coeff.T))

If you force gpu4pyscf to use linalg.solve in this place, gpu4pyscf can produce results closer.