gto/ecp: fix premature convergence in radial quadrature at large exponents#3250
Merged
sunqm merged 2 commits intoJun 13, 2026
Merged
Conversation
…quadrature
The adaptive Gauss-Chebyshev radial quadrature in ECPtype1_cart,
ECPtype2_cart, and ECPtype_so_cart declared per-primitive-pair
convergence the first time CLOSE_ENOUGH(plast, prad) held between two
successive doubled grids. For sharply-peaked integrands (combined AO+ECP
exponent (2a+g) large, n >= 2, especially at higher AO angular momentum)
two coarse rules can agree to 1e-12 relative while both still
under-sample the peak at r ~ (2a+g)^-1/2, freezing in a wrong answer.
The n=1 channel happened to stay exact because r * exp(-c r^2) is far
smoother under the log-quadratic change of variables used.
Promote the per-pair flag to a counter and require two consecutive
CLOSE_ENOUGH matches before the pair is removed from refinement.
Also rewrite CLOSE_ENOUGH as a clean relative test (max of |x|, |y|
instead of |y| in the denominator, plus the absolute fallback). The
previous 1e-12 absolute floor falsely matched high-l / high-exponent
integrals whose true magnitude was below it; the new form catches the
0 == 0 case naturally and has no magic absolute threshold.
For a Kr atom with a 48-term large-exponent local ECP this removes a
constant ~4e-5 Ha absolute-energy error. Worst-case relative error on a
full alpha x g x n x l sweep (local + semilocal channels, alpha in
{1, 1e2, 1e4, 3e5}, g in {1e1, 1e3, 1e5, 1e7}, n in {1, 2}, l in
{0, 1, 2}) drops from 5.5e-5 to 4.7e-13. LANL2DZ benchmark (Cu + 6H)
shows no measurable slowdown.
Closes pyscf#3249.
…ntegrals
For a same-center single-primitive AO of angular momentum l with a
one-term local or semilocal ECP channel c * r^(n-2) * exp(-g r^2), the
matrix element factorises so that the ratio I(g1) / I(g2) at fixed
alpha is
((2 alpha + g2) / (2 alpha + g1)) ** ((n + 2 l + 1) / 2)
independent of the AO normalisation. This is a stringent and cheap
closed-form check on the radial quadrature.
Sweep n in {1, 2}, l in {0, 1, 2}, alpha in {1, 1e2, 1e4, 3e5},
g in {1e1, 1e3, 1e5, 1e7}, local and semilocal channels.
Before the convergence-check fix worst-case rel. err was ~5.5e-5
(local) and ~1e-1 (semilocal at l = 2, large g); now 4.7e-13.
Refs pyscf#3249.
sunqm
approved these changes
Jun 13, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #3249.
Summary
ECPtype1_cart,ECPtype2_cart, andECPtype_so_cartdeclared per-primitive-pairconvergence the first time
CLOSE_ENOUGH(plast, prad)held. Forsharply-peaked integrands at large combined AO + ECP exponent,
two coarse rules could agree to 1e-12 while both still
under-sampled the peak — freezing in a wrong answer. Promote the
per-pair flag to a counter and require two consecutive matches.
CLOSE_ENOUGHas a clean relative test|x − y| ≤ 1e-12 · max(|x|, |y|). The previous1e-12absolutefallback falsely matched high-
l/ high-exponent integrals whosetrue magnitude was below it; the new form catches the
0 == 0case naturally and has no magic absolute threshold.
semilocal channels across n ∈ {1, 2}, l ∈ {0, 1, 2}, α ∈ {1, 1e2,
1e4, 3e5}, g ∈ {1e1, 1e3, 1e5, 1e7}.
Before / after
Same-center s-shell with one-term local ECP
n · g · 1.0, α = AOexponent (existing PySCF radial integrator vs. closed form
<s_α|r^(n-2) e^(-g r²)|s_α>):Closed-form ratio test (
I(g₁)/I(g₂)against((2α+g₂)/(2α+g₁))^((n+2l+1)/2))across the full sweep — worst-case relative error:
Test plan
test_large_exponent_ecp_closed_formpasses (worst rel. err4.7e-13, target 1e-10).
pyscf/gto/test/test_ecp.py(14 tests) continues to pass.pyscf/gto/test/test_basis_parser.py,pyscf/gto/test/test_moleintor.py,pyscf/grad/test/test_rhf.pyECP-touching tests continue to pass.
ECPscalarbuild) shows nomeasurable slowdown — 5.0 vs 5.1 ms, within run-to-run noise.
Co-Authored-By: Claude Opus 4.7 noreply@anthropic.com