Fix uninitialized read in int3c1e_ip1_r6_origk y component - #132
Fix uninitialized read in int3c1e_ip1_r6_origk y component#132susilehtola wants to merge 2 commits into
Conversation
CINTgout1e_int3c1e_ip1_r6_origk computes the y derivative s[1] with the term 6*g48[ix]*g76[iy]*g3[iz], but g76 was never filled: the G1E_D_I call for it was missing, so the y component of the integral was computed from stale scratch-buffer contents. The x and z components were unaffected. g76 must hold G1E_D_I applied to g12, consistent with the pointer layout in which every bra-derivative array sits 64 slots above its base (g64=D(g0), g67=D(g3), g79=D(g15), g112=D(g48), g124=D(g60), g127=D(g63)). Verified against finite differences of int3c1e_r6_origk with each shell on its own center (-d/dR_i <i|r^6|jk> = <nabla i|r^6|jk>): before the fix, the y component deviated from the numerical derivative by up to ~0.2 for a small s/p/d/f test basis; after the fix all components agree to 1e-10. int3c1e_ip1_r6_origk is used by PySCF (>= 2.13) for the C4 term of GTH pseudopotential nuclear gradients (pyscf/gto/pp_int.py and pyscf/pbc/gto/pseudo/pp_int.py), so gradients of elements whose GTH local potential includes the r^6 term were affected. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TTc5De7stABc2q58NHV22E
…_r*_origk
test_int3c1e.py takes its reference values from the libcgto library of
the installed pyscf, which only provides the integrals of the libcint
release it was built against. The int3c1e_ipip1/ip1ip2/ipvip1/ipip2
integrals were added after the 6.1.3 release, so the test crashed with
an undefined-symbol AttributeError on CI. Skip integrals that the
installed pyscf does not provide.
Also add self-contained finite-difference checks for
int3c1e_ip1_r{2,4,6}_origk: with the bra shell on an atom that carries
neither j nor k, -d/dR_i of the base integral equals the ip1 integral
exactly. This validates int3c1e_ip1_r6_origk without relying on
pyscf, whose bundled libcint <= 6.1.3 computes its y component from an
uninitialized buffer (fixed in the previous commit); the new check
fails on the unfixed code and passes with the fix.
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TTc5De7stABc2q58NHV22E
|
Note on the earlier ubuntu-latest CI failure: it was pre-existing on master, not caused by this fix. The second commit fixes this by skipping integrals the installed pyscf does not provide, and adds self-contained finite-difference checks for |
This PR is an offshoot of sunqm/qcint#24: it turned out that libcint was missing a term that exists in qcint, which is added in this PR. libcint is currently returning incorrect results.
AI summary
CINTgout1e_int3c1e_ip1_r6_origk(src/cint3c1e_a.c) computes the y-components[1]with the term6*g48[ix]*g76[iy]*g3[iz], butg76is never initialized — theG1E_D_Icall for it is missing from the recurrence list. The y component of the integral is therefore computed from stale scratch-buffer contents; the x and z components are unaffected.g76must holdG1E_D_Iapplied tog12, consistent with the pointer layout in which every bra-derivative array sits 64 slots above its base (g64=D(g0),g67=D(g3),g79=D(g15),g112=D(g48),g124=D(g60),g127=D(g63)). This PR adds the missing line:Validation. With one single-primitive shell per center,
-d/dR_i <i|r^6|jk> = <nabla i|r^6|jk>exactly, so central differences ofint3c1e_r6_origkgive an independent reference. Before the fix, the y component deviates from the numerical derivative by up to ~0.2 (on values of order 1) for an s/p/d/f test basis; after the fix, all components of all shell tuples agree to 1e-10 (the finite-difference truncation level). The fixed code also agrees with qcint's machine-generated implementation of the same integral to ~1e-10.Impact. PySCF >= 2.13 uses
int3c1e_ip1_r6_origkfor the C4 term of GTH pseudopotential nuclear gradients (pyscf/gto/pp_int.py:99,pyscf/pbc/gto/pseudo/pp_int.py:186), so molecular and PBC gradients of elements whose GTH local potential includes the r^6 term have had a silently wrong y component.🤖 Generated with Claude Code
https://claude.ai/code/session_01TTc5De7stABc2q58NHV22E