`DF.reset()` keeps a stale `auxmol` after in-place geometry change — silently wrong SCF energies (diverges from PySCF CPU behavior)

[Isitea]

Summary

After changing the geometry of mf.mol in place (mol.set_geom_(...)) and
calling mf.reset(), a density-fitted calculation in gpu4pyscf rebuilds
_cderi with the auxiliary basis still sitting at the old coordinates.
The SCF then converges normally and returns a silently wrong energy — for a
1e-3 Bohr displacement of one atom in water/def2-SVP we get an error of
~2e-2 Ha, with no warning or error raised.

The same sequence is handled correctly by CPU PySCF, so code that follows the
documented scanner/reset pattern gives correct results on CPU and corrupted
results on GPU. We hit this in production: our finite-difference Hessians
(displace → reset() → kernel() loop) were entirely corrupted, producing
non-physical frequencies in the 10^4–10^5 cm⁻¹ range before we tracked it down.

Minimal reproduction

from pyscf import gto, dft

WATER = "O 0 0 0; H 0.9572 0 0; H -0.2400 0.9266 0"

def make_mol(displaced=False):
    mol = gto.M(atom=WATER, basis="def2-svp", unit="Angstrom", verbose=0)
    if displaced:
        c = mol.atom_coords(unit="Bohr").copy()
        c[0, 0] += 1e-3
        mol.set_geom_(c, unit="Bohr")
        mol.build(dump_input=False)
    return mol

# Reference: fresh mf at the displaced geometry
e_fresh = dft.RKS(make_mol(True), xc="wb97x").density_fit().to_gpu().kernel()

# Reuse: kernel at equilibrium, displace in place, reset, kernel again
mf = dft.RKS(make_mol(), xc="wb97x").density_fit().to_gpu()
mf.kernel()
c = mf.mol.atom_coords(unit="Bohr").copy()
c[0, 0] += 1e-3
mf.mol.set_geom_(c, unit="Bohr")
mf.mol.build(dump_input=False)
mf.reset()
e_reused = mf.kernel()

print(f"fresh  = {float(e_fresh):.10f}")
print(f"reused = {float(e_reused):.10f}")
print(f"|diff| = {abs(float(e_reused) - float(e_fresh)):.2e} Ha")
# stale auxiliary basis: max deviation equals the displacement
import numpy as np
aux = mf.with_df.auxmol.atom_coords(unit="Bohr")
print(f"auxmol max deviation from current mol geometry: "
      f"{np.abs(aux - mf.mol.atom_coords(unit='Bohr')).max():.2e} Bohr")

Observed (gpu4pyscf 1.7.4, cupy 14.1.1, pyscf 2.13.1, RTX 3070, CUDA 12.4/WSL2):

fresh  = -76.3496849152
reused = -76.3700807832
|diff| = 2.03e-02 Ha
auxmol max deviation from current mol geometry: 1.00e-03 Bohr

Expected: |diff| at SCF-convergence level (the same script with
.density_fit() removed, or run on CPU PySCF, gives ~1e-11 Ha).

Controls we ran (same displace+reset+kernel sequence):

variant	result
GPU DF + SMD	wrong (2.04e-2 Ha)
GPU DF only	wrong (2.03e-2 Ha)
GPU SMD only (no DF)	correct
GPU bare RKS	correct
CPU DF (+/- SMD)	correct

So the issue is specific to the density-fitting path on GPU; solvent models
are unaffected on their own.

Root cause

CPU PySCF df.DF.build regenerates the auxiliary Mole unconditionally
from the current self.mol on every build
(pyscf/df/df.py,
auxmol = self.auxmol = addons.make_auxmol(self.mol, self.auxbasis)),
so a stale auxmol can never survive a rebuild.

gpu4pyscf df.DF.build only regenerates it when it is missing:

if auxmol is None:
    self.auxmol = auxmol = addons.make_auxmol(mol, self.auxbasis)

and DF.reset(mol=None) only clears auxmol when a new mol object is
passed:

def reset(self, mol=None):
    if mol is not None:
        self.mol = mol
        self.auxmol = None
        ...
    self._cderi = None
    ...

mf.reset() (no argument) therefore keeps the old-geometry auxmol (and on
current master also intopt/j_engine), while _cderi is rebuilt against
it — producing wrong 3c2e integrals for the new geometry with no diagnostic.

Suggested fix

Either

• clear self.auxmol (and intopt/j_engine) unconditionally in
  DF.reset(), or
• rebuild auxmol from self.mol unconditionally in DF.build(), matching
  CPU PySCF semantics.

If keeping the cache is intentional for scanner performance, a geometry
consistency check between mol and auxmol at build time would at least
turn silent corruption into a loud error.

Environment

Minimal reproduction above was run on:

• gpu4pyscf 1.7.4 (code path unchanged on current master as of 2026-07-23)
• pyscf 2.13.1, cupy 14.1.1
• NVIDIA GeForce RTX 3070, CUDA 12.4, WSL2 (Linux 6.6)

The original production corruption (finite-difference Hessians) was also
observed on RTX 4070 Ti and A100 SXM systems with CUDA 13.x, so the issue is
not specific to a particular GPU architecture or CUDA major version.

Possibly related

#552 reported numerical vs analytical Hessian disagreement with DF on GPU.
If that workflow reused a live mf across displacements (in-place
set_geom_ + reset()), the stale-auxmol behavior described here could
explain it.

[Isitea]

Follow-up while preparing a fix: the bug turned out to have a second layer beyond the conditional cache invalidation in DF.reset().

Unconditionally discarding auxmol on reset is not sufficient: make_auxmol re-parses mol._atom (the input record), and after to_gpu() mf.mol and mf.with_df.mol are distinct objects sharing _env. An in-place mf.mol.set_geom_() therefore updates the runtime coordinates but not with_df.mol._atom, so even a freshly rebuilt auxmol sits at the old coordinates. (The AO side is unaffected — SortedMole.from_mol works from mol._env.)

Fix for both layers, with regression tests and A/B validation (2.03e-2 Ha → 4.3e-11 Ha on the reproduction above): #828

[sunqm]

This is a misuse of mol.set_geom_(). Although this function supports modifying mol instance inplace, the resulting mol should still be treated as a new Mole instance. You should call mf.reset(mol) to update the state of mf and its attributes.

Silently modifying the contents of mol and expecting all dependent objects to update automatically is not guaranteed.

[Isitea]

Following up on the discussion in #828 — the API contract makes sense, and we have since updated our code to follow the mf.reset(mol) convention. It is working well. Thank you again for the guidance.

If it's not too much, I'd like to share one small reflection from our side, only in case it is ever useful. The reason we didn't notice our mistake for so long was that the same usage happened to give correct results with CPU PySCF, so when our workflow moved to GPU4PySCF we simply assumed the behavior would carry over. That was our own assumption to check, of course — I only mention it because other users migrating CPU-based scripts might walk the same path we did.

I imagine the GPU-side caching is there for good performance reasons, so please take this as nothing more than a gentle thought: if someday there happens to be a low-cost way to make the two backends behave alike in this corner (or to make the difference visible when it matters), it might quietly save someone a debugging journey like ours. If not, that is completely understandable, and we are perfectly happy with the documented convention as clarified.

As for this issue itself, please feel free to close it whenever you see fit — from our side everything is resolved, and I'll leave the decision entirely in your hands.

Thanks again for your time and patience with this issue.