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

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

Start with the minimal Python example and reproduce the Cartesian-versus-spherical gradient comparison on the stated PySCF and gpu4pyscf versions. Trace the Coulomb 2e gradient contribution identified in the analysis, then verify that Cartesian results agree with PySCF while the spherical comparison remains consistent.