DF-RKS returns wrong energies (+59816 Ha, no exception) when run under a non-blocking CuPy stream — v1.8.1, standalone reproducer

Running GPU4PySCF while a non-default cupy.cuda.Stream is current can produce a wrong energy with
no exception raised. On gpu4pyscf-cuda12x==1.8.1 / cupy-cuda12x==14.1.1 / pyscf==2.14.0,
H100 (sm_90), the attached reproducer converges H2O/def2-TZVPP/B3LYP DF-RKS to
−76.4666819831212 Ha on the default stream, and returns +59816.219055041925 Ha —
converged=False, no traceback, exit code 0 — when only the object construction is moved inside
with cupy.cuda.Stream(non_blocking=True):. The same script on an RTX PRO 6000 Blackwell (sm_120)
does the same thing, and additionally returns wrong per-cycle energies when the object is built on
the default stream and only the per-cycle work is moved into a single non-default stream, with
no concurrency at all.

non_blocking=True is load-bearing — please keep it if you paste the one-liner. CuPy's
Stream.__init__ defaults non_blocking=False, and a bare cupy.cuda.Stream() is the blocking
stream, which my runs report as clean on both cards (that is the setup_in_stream_blocking row
below, and the workaround at the end).

I understand that non-default streams may simply not be a supported usage mode: #561 suggests the
library may be intentionally default-stream-only since the pre-allocated streams were removed, and
if that is the answer I am happy to take it. My ask is smaller than stream safety: this outcome
is indistinguishable from a correct calculation except for a convergence flag, so a cheap guard at
SCF entry plus one line in the README saying which streams are supported would close it entirely.

Measured

Two cards, one run each, eight cases per run, each case in a fresh process. Both dispatches ran
one program — blob ce59e8a1c495e309025c812522b230f12ca17cc6 of
repro_stream_unsafety.py — and the copy attached here differs from that blob in comments and
docstrings only, every executable statement byte-identical (checked by comparing the two token
streams with comment and docstring tokens removed, not by eye). Same pins on both runs:
gpu4pyscf-cuda12x 1.8.1, pyscf 2.14.0, cupy-cuda12x 14.1.1, numpy 2.2.6, scipy 1.15.3,
cuda_driver 13000, cuda_runtime 12090, cusolver 11.7.3, devices 1, num_devices 1,
CUDA_VISIBLE_DEVICES=0. What differs is the card and its driver build: NVIDIA H100 80GB HBM3 /
compute_cap 9.0 / driver 580.126.09, and NVIDIA RTX PRO 6000 Blackwell Workstation Edition /
compute_cap 12.0 / driver 580.126.20.

case	what it does	H100 (sm_90)	RTX PRO 6000 Blackwell (sm_120)
control	build + converge entirely on the default stream	converged, −76.4666819831212	converged, −76.46668198312115
setup_in_stream	build inside Stream(non_blocking=True), converge on default; library already warm	wrong number, no exception, converged=False, +59816.219055041925 (|ΔE| 5.989e+04)	wrong number, no exception, converged=False, +113328.8142946389 (|ΔE| 1.134e+05)
setup_in_stream_blocking	same, Stream(non_blocking=False)	clean, converged, |ΔE| 0.0	clean, converged, |ΔE| 0.0
setup_in_stream_cold	same as setup_in_stream, but no warm-up kernel precedes it	clean, converged, |ΔE| 1.279e-13	clean, converged, |ΔE| 2.842e-14
cycle_in_streams (w=4)	per-cycle SCF work inside 4 non-default streams, event-ordered	clean, |ΔE| 0.0 (16 cycles compared)	wrong, |ΔE| 50.005599042155325 (16 compared)
cycle_in_streams_width1	the same with one stream, no concurrency	clean, |ΔE| 0.0 (4 compared)	wrong, |ΔE| 50.005599042155325 (4 compared)
cycle_blocking_stream (w=4)	the same on blocking streams	clean, |ΔE| 0.0 (16 compared)	clean, |ΔE| 8.527e-14 (16 compared)
cycle_patched_width1	experimental, patches one library internal in-process (see "Where I looked")	clean, 0.0 — but the unpatched arm is clean here too, so it discriminates nothing on this card	clean, 0.0, where the unpatched width-1 arm above is wrong by 50.0056 Ha

The comparison tolerance is 1e-8 Ha throughout. On the cycle cases |ΔE| is the largest disagreement
between the in-stream and default-stream per-cycle energies. On the setup cases it is against a
default-stream reference, but not the same reference in every row, and this is why the cold row
does not read exactly 0.0: setup_in_stream and setup_in_stream_blocking compare against a
reference converged in the same process moments earlier (reference from in_process_warmup), while
setup_in_stream_cold must not run a warm-up kernel by definition, so it compares against the
control case's energy passed in from the parent process (reference from parent_supplied). The
1.279e-13 / 2.842e-14 there is cross-process run-to-run noise, not a defect.

One note on the magnitude, because it is not stable: an earlier three-case version of this script,
on the same sm_120 card, returned +321837.2422047062 Ha for a setup case whose stream usage
is code-identical to the current one (the surrounding reference handling changed; the CUDA op
sequence did not). So the failure reproduces reliably and the particular wrong number does not — I
would not read anything into either value beyond "it is not the energy".

Timing, for sizing a re-run: my CI job's work block — which includes install and clone, not just the
calculation — was 86 s on the H100 and 113 s on the RTX PRO 6000.

Three splits that came out of the runs

• Blocking vs non-blocking — holds on both cards. setup_in_stream (non_blocking=True) is
  wrong and setup_in_stream_blocking (non_blocking=False) is exactly right — same process, same
  reference, one flag apart, on sm_90 and on sm_120. Same on the cycle path: the blocking variant
  is clean on both cards, including on the card where the non-blocking one is not. That is a
  workaround users can apply today: if you must run under a stream, make it a blocking one.
• Warm vs cold — holds on both cards. setup_in_stream and setup_in_stream_cold differ only
  in whether a full default-stream kernel() ran earlier in the process. The warm one is wrong;
  the cold one matches to 1.3e-13 (sm_90) and 2.8e-14 (sm_120) Ha. So the corruption needs a
  warmed-up library and a non-synchronising stream.
• Architecture — the cycle path differs. With the object built on the default stream and only
  the per-cycle work moved into a non-default stream, sm_90 is clean and sm_120 is wrong by
  50.005599042155325 Ha. That figure is identical to all printed digits at width 4 and at
  width 1, so it looks like a deterministic offset rather than a race — but I have one run per
  card and I am not claiming a mechanism for it. Two things are worth flagging about the width-1
  case specifically: every producer (s1e/h1e/dm, built on the default stream) is
  event-ordered into the consuming stream before it is read, and there is one stream, so
  there is no concurrency in it at all. Whatever this is, it is not my own cross-stream read.

On the earlier reports, and why I am asking for a contract rather than a fix

#548 ("problem with RHF multi-gpu", Oct 2025) has a similar shape — H2O, an energy of ~−436,937 Ha,
50 unconverged cycles, no exception — and was closed by #561 ("Remove pre-allocated streams",
merged 2025-11-07), whose description says the pre-allocated streams "can break the device
synchronization". Both are in 1.8.1.

I want to be careful not to oversell that parallel, because it may well be a different bug: #548 was
on 2× NVIDIA P100 (sm_60), multi-GPU, with cart=1, and the reporter found that cart=0 made
it go away and suspected _vhfopt/screening. My reproducer pins CUDA_VISIBLE_DEVICES=0, so it
never takes the multi-GPU path at all. Treat #548 as context for why I went looking, not as the
same defect.

What actually motivates the ask is narrower and is visible from outside: the library creates a
non-default io_stream in df/df.py and threads get_current_stream() into a great many kernel
launches, which reads to a caller as "the current stream is honoured" — and I could not find any
statement in the README or docs about which streams are supported. One sentence in the docs, plus a
guard that says so at runtime, would stop a wrong energy from reaching someone's results.

Already fixed on master?

Not as far as I can tell, though I have only checked the stream-relevant axis by inspection and have
not re-run the reproducer against master.

I re-checked the nine files on the exercised path against master
70a48ab6e39b12ab9dd86c19989368350f40175a ("Optimize pp int (#883)", 2026-09-01):
lib/cusolver.py, lib/multi_gpu.py and dft/rks.py are byte-identical to v1.8.1, and
lib/cupy_helper.py, df/df.py, df/df_jk.py, scf/hf.py, scf/jk.py and dft/numint.py have
changed. What I checked in those six is narrower than "nothing changed": no stream, event or
synchronisation statement differs — grepping the six git diff v1.8.1 <master> file diffs for
stream|event|synchroniz|record|wait_event|cusolver|get_current returns zero hits (zero
case-insensitively too), df/df.py's io_stream block is untouched, and scf/jk.py still has the
same six get_current_stream() sites. What did change in those six is the
lr_factor/sr_factor range-separated-Coulomb API threaded through df/df.py, df/df_jk.py,
scf/jk.py and scf/hf.py; a get_jk restructure in scf/hf.py (the module-level _get_jk
folded into SCF.get_jk, plus a get_veff-before-get_hcore reordering in _kernel);
_nearest_power2 moving to lib.utils; direct_scf_tol becoming a deprecated keyword-only
argument of DF.build; and two additions unrelated to this path — a vec_dot kernel in
lib/cupy_helper.py, and is_hybrid_xc/is_nlc plus an inplace kwarg in dft/numint.py.

Since a stream-handling fix has landed in that window, to be explicit about it: #868 ("Fix issues
related to cuda indexing, cutlass stream handling", 630d274, 2026-08-24) is the only commit since
the 1.8.1 release whose subject mentions streams, and it touches
lib/cupy_helper/grouped_dot.cu, lib/cupy_helper/grouped_gemm.cu, gvhf-rys/rys_contract_k.cu,
pbc/rys_contract_k.cu, nac/tdrhf_grad_nacv.py and scf/soscf.py — none of the nine files above,
and nothing on the DF-RKS path this reproducer exercises.

Reproducer

Attached: repro_stream_unsafety.py. It imports only gpu4pyscf, pyscf, cupy and the standard
library, so it runs in a bare gpu4pyscf environment, prints the full version matrix on a
REPRO_START line and one machine-readable REPRO_CASE line per case, and runs each case in a
fresh process (a poisoned CUDA context must not silence later cases). Most of its length is
comments explaining what each case does and does not show.

The setup_in_stream case, minus the harness, is:

def build_mf():
    """A converged-capable H2O/def2-tzvpp DF-RKS object. Every CUDA-touching statement of the
    setup path is inside this function, so a case can put the WHOLE of it inside a stream context
    or none of it."""
    from pyscf import gto
    import gpu4pyscf.dft as gpu_dft

    mol = gto.M(atom=ATOM, basis=BASIS, unit="Angstrom", verbose=0)
    mf = gpu_dft.RKS(mol, xc=XC)
    mf.grids.level = GRID_LEVEL
    mf.grids.prune = None
    mf.conv_tol = CONV_TOL
    mf.max_cycle = MAX_CYCLE
    mf.verbose = 0
    mf = mf.density_fit()
    mf.conv_tol = CONV_TOL
    mf.max_cycle = MAX_CYCLE
    mf.grids.build()
    wdf = getattr(mf, "with_df", None)
    if wdf is not None and hasattr(wdf, "build"):
        wdf.build()
    return mf

# ... and the body of `_setup_case`, after it has converged its default-stream reference:

    stream = cupy.cuda.Stream(non_blocking=bool(non_blocking))
    with stream:
        mf = build_mf()
    stream.synchronize()
    e = _f(mf.kernel())

Both spans are lifted verbatim; what is elided between them is the reference-energy bookkeeping
(which of three reference sources the case uses) and, after them, the comparison against
CORRUPT_ATOL. ATOM/BASIS/XC/GRID_LEVEL/CONV_TOL/MAX_CYCLE are module constants: the
water geometry from gpu4pyscf's own tests/test_rks.py, def2-tzvpp, b3lyp, grid level 3,
conv_tol=1e-10, max_cycle=50. non_blocking is True for setup_in_stream and False for
setup_in_stream_blocking; _f is a device-scalar-to-float helper.

CUDA_VISIBLE_DEVICES=0 python3 repro_stream_unsafety.py          # all eight cases
CUDA_VISIBLE_DEVICES=0 python3 repro_stream_unsafety.py --case setup_in_stream

CUDA_VISIBLE_DEVICES=0 is not cosmetic: with more than one visible device the calculation routes
through lib/multi_gpu.py's ThreadPoolExecutor, which is a different code path, and the script
reports both devices and gpu4pyscf.__config__.num_devices so you can see which one a run took.
Budget a few minutes on a JIT-cold machine; each case has its own 275 s hang ceiling.

Where I looked, and what one experiment showed

I would rather show you the one experiment than argue for it.

Reading the exercised path, gpu4pyscf/lib/cusolver.py calls cusolverDnDsygvd through ctypes
on a per-device cached handle (device.get_cusolver_handle()) and, as far as I can see, never
calls cusolverDnSetStream on it — whereas CuPy's own cuSOLVER wrappers call _setStream(handle)
before every call, on that same per-device handle. So the handle's stream binding is whatever
last bound it, which need not be the stream that is current when gpu4pyscf uses it.

So the reproducer has an explicitly experimental case, cycle_patched_width1: in that child
process only, it binds the handle to the current stream (via
cupy_backends.cuda.libs.cusolver.setStream) immediately before each gpu4pyscf.lib.cusolver.eigh,
and nothing else changes. The patch targets gpu4pyscf.lib.cusolver — that is a statement about
what the script does, not a run observation. What the run reports is a count: the case line says
the patch was applied "on 1 module(s)", which is the number I expected from reading 1.8.1. If you
want the module names for a run of your own, invoke that case on its own —
--case cycle_patched_width1 additionally prints a REPRO_JSON record carrying rebound_modules
(the all-cases invocation above does not emit it). The case also re-checks afterwards that no
module imported later escaped the patch, and reports itself as error rather than as a
measurement if any did. On sm_120, in the same run:

• cycle_in_streams_width1, unpatched: wrong by 50.005599042155325 Ha
• cycle_patched_width1, patched: 0.0

Scope of that result:

• This is one architecture and one run. On the H100 the unpatched cycle path is already clean,
  so the patched case there discriminates nothing and offers no support in either direction.
• It covers the cycle path only. It says nothing about setup_in_stream, which is wrong on
  both cards and is the headline of this report — whatever is happening there, this does not
  explain it.
• The patched case computes its own default-stream reference with the patch installed. So what
  is measured is that under the patch the in-stream and default-stream trajectories agree, where
  unpatched they differ by 50.0056 Ha. That the patched numbers also equal the unpatched control
  is not something this case measures, and I am not asserting it.

If this lead is useful, I am happy to turn it into a PR — and equally happy to be told the handle
is deliberately left unbound and the real answer is elsewhere.

Ask

A guard at SCF entry that raises or warns when a non-default stream is current, and one README line
stating the supported stream contract. If the answer is "default stream only", saying so in code
and in the docs is all I am after.

One caveat on scope, since my own results argue against the broadest version of that guard: a
blanket check on cupy.cuda.get_current_stream() being non-default would also fire on the
blocking stream, which is clean in every case I measured and is the workaround I would
otherwise recommend. Narrowing it to non-blocking streams, or simply accepting the over-warn,
are both fine by me — I do not want to prescribe the mechanism, only to have the contract stated
somewhere.

Happy to test any candidate fix on either card and report back, and happy to open the
cusolverDnSetStream change as a PR if you want to look at it.

repro_stream_unsafety.py.txt

Start by running the attached repro_stream_unsafety.py, especially the setup_in_stream and cycle_in_streams_width1 cases, then inspect stream handling in df/df.py, scf/jk.py, and the other named path files. Check the README and docs for an existing stream-support contract; done means the unsupported usage is clearly documented and the chosen runtime behavior is covered by a regression test or guard.