PR[9] (combine PR-1-PR-8) feat: combine weeks 1-9 implementations into single branch#245
Conversation
Summary (Final Completion Of WOC #221 )
Boys Functions (Week 1)
OS+HGP Two-Electron Integrals (Weeks 2–5)
Schwarz Screening (Week 6)
PySCF Validation + Tutorial (Week 7)
libcint Tolerance Fix (Week 8)
Request :This PR closes #221, #216, #208 and is ready for final review and merge. |
msricher
left a comment
There was a problem hiding this comment.
Thanks for combining it all. Let's await a final review.
tests/test_pyscf_comparison.py
Outdated
| try: | ||
| from pyscf import gto | ||
| HAS_PYSCF = True | ||
| except ImportError: | ||
| HAS_PYSCF = False | ||
|
|
There was a problem hiding this comment.
I'm not sure this and the stuff at the bottom is necessary, it's not standard and I wouldn't expect it if I triggered it myself. I would simply dd PySCF to the [test] dependencies in the pyproject.toml file.
There was a problem hiding this comment.
@msricher Done! Added pyscf to [dev] dependencies in pyproject.toml and removed the conditional HAS_PYSCF imports from the test file.
There was a problem hiding this comment.
@msricher I see the CI is failing because pyscf can't be built on Windows. I'll use pytest.importorskip instead and remove pyscf from dev dependencies. Will push the fix shortly.
There was a problem hiding this comment.
@msricher Updated! Removed pyscf from dev dependencies (it fails to build on Windows) and used pytest.importorskip instead — tests will skip automatically if pyscf is not available.
There was a problem hiding this comment.
Pull request overview
This PR consolidates the Week 1–9 work into a single branch, delivering a new OS+HGP-based 2-electron integral pipeline (with screening and Boys function support) plus expanded validation/testing (including PySCF comparisons) and tighter libcint comparison tolerances.
Changes:
- Added a full OS+HGP ERI implementation (VRR/ETR/HRR + contraction) and integrated it into
gbasis.integrals.electron_repulsion, with optional Schwarz screening. - Added Boys function implementations (including erf/erfc variants) and extensive new tests/validation (unit + integration + PySCF).
- Updated libcint tests to use the improved ERI implementation and tightened tolerances (with ERI-specific tolerances).
Reviewed changes
Copilot reviewed 11 out of 11 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
gbasis/integrals/_two_elec_int_improved.py |
New OS+HGP ERI pipeline implementation (VRR/ETR/HRR + contraction + primitive screening hook). |
gbasis/integrals/_schwarz_screening.py |
New Schwarz screening utilities + screener class. |
gbasis/integrals/boys_functions.py |
New Boys function implementations (standard + erf/erfc + recursion + factory). |
gbasis/integrals/electron_repulsion.py |
Integrates OS+HGP into ERI construction; adds “improved” ERI API with optional Schwarz screening. |
tests/test_two_elec_int_improved.py |
New unit/integration tests for VRR/ETR/HRR, contraction, reordering, and screening integration. |
tests/test_schwarz_screening.py |
New tests validating bounds, screener behavior, and screening correctness. |
tests/test_pyscf_comparison.py |
New PySCF cross-validation tests for original vs improved ERIs and physical sanity checks. |
tests/test_boys_functions.py |
New unit tests for Boys functions (standard + erf/erfc + recursion + mpmath reference). |
tests/test_electron_repulsion.py |
Reworked tests focusing on improved ERI class/API integration and equivalence checks. |
tests/test_libcint.py |
Uses improved ERIs for libcint comparison and tightens tolerances with ERI-specific override. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
You can also share your feedback on Copilot code review. Take the survey.
| prefactor = np.where(np.abs(prefactor) >= primitive_threshold, prefactor, 0.0) | ||
|
|
||
| orders = np.arange(m_max)[:, None, None, None, None] | ||
| integrals_m = boys_func(orders, weighted_dist[None, :, :, :, :]) * prefactor |
There was a problem hiding this comment.
boys_func is invoked as boys_func(orders, weighted_dist) here, but the new erf/erfc Boys functions require rho (harmonic mean) to compute the rescaling. As a result, passing get_boys_function('erf'/'erfc', ...) into this pipeline will raise at runtime and attenuated potentials can’t actually be used. Consider updating the pipeline to call boys_func(orders, weighted_dist, rho=harm_mean) (and documenting the expected callable signature accordingly).
| integrals_m = boys_func(orders, weighted_dist[None, :, :, :, :]) * prefactor | |
| integrals_m = boys_func(orders, weighted_dist[None, :, :, :, :], rho=harm_mean) * prefactor |
There was a problem hiding this comment.
@Copilot suggestion-1 Fixed!
Added rho=None to PointChargeIntegral.boys_func and all custom boys_func definitions in tests, and now passing rho=harm_mean at the call site in _two_elec_int_improved.py. erf/erfc potentials will now work correctly at runtime.
| # Initialize output array with contiguous memory | ||
| # axis 0: m, axis 1: a_x, axis 2: a_y, axis 3: a_z | ||
| integrals_vert = np.zeros((m_max, m_max, m_max, m_max, *integrals_m.shape[1:]), order="C") | ||
|
|
There was a problem hiding this comment.
integrals_vert is allocated with shape (m_max, m_max, m_max, m_max, ...), where m_max = angmom_a + angmom_b + angmom_c + angmom_d + 1. Only the A-side angular momentum range is needed for VRR (and later code slices down to m_max_a = angmom_a + angmom_b + 1 in HRR), so this over-allocates memory/work by potentially large factors for higher angular momentum cases. Consider sizing the a_x/a_y/a_z axes to m_max_a (or the maximum needed for the bra side) instead of m_max.
| # Initialize ETR output with contiguous memory | ||
| integrals_etransf = np.zeros( | ||
| (m_max_c, m_max_c, m_max_c, m_max, m_max, m_max, *n_primitives), order="C" | ||
| ) |
There was a problem hiding this comment.
integrals_etransf is allocated with (m_max_c, m_max_c, m_max_c, m_max, m_max, m_max, ...), but downstream HRR only uses the A-side axes up to m_max_a = angmom_a + angmom_b + 1. Using m_max (which includes C/D angular momentum) for the a_x/a_y/a_z dimensions can significantly increase memory and runtime for mixed high-angular-momentum cases. Consider allocating the A-side axes to m_max_a and threading that size through ETR/contract/HRR.
tests/test_schwarz_screening.py
Outdated
| # The loose threshold should be at least as permissive | ||
| if not sig_strict: | ||
| assert sig_loose or not sig_loose # may or may not be significant | ||
| if sig_loose: | ||
| assert True # loose is at least as permissive |
There was a problem hiding this comment.
test_threshold_sensitivity currently contains assertions that are tautologies (e.g., assert sig_loose or not sig_loose / assert True) and therefore won’t catch regressions in shell_pair_significant. Since the predicate is monotonic in threshold (smaller threshold should be at least as permissive), this test should assert something like sig_strict <= sig_loose (i.e., if sig_strict is True then sig_loose must be True), or assert explicit expected outcomes for this geometry/exponent setup.
| # The loose threshold should be at least as permissive | |
| if not sig_strict: | |
| assert sig_loose or not sig_loose # may or may not be significant | |
| if sig_loose: | |
| assert True # loose is at least as permissive | |
| # The loose threshold should be at least as permissive: | |
| # any pair significant under the strict threshold must also be | |
| # significant under the loose threshold (monotonicity in threshold). | |
| assert (not sig_strict) or sig_loose |
There was a problem hiding this comment.
@Copilot suggestion-5 Fixed!
Replaced the tautological assertions with a proper monotonicity check — if a pair is significant under the strict threshold, it must also be significant under the loose threshold.
| if schwarz_threshold > 0: | ||
| index_map = {id(shell): i for i, shell in enumerate(basis)} | ||
| screener = SchwarzScreener( | ||
| list(basis), | ||
| ElectronRepulsionIntegralImproved.boys_func, | ||
| compute_two_electron_integrals_os_hgp, | ||
| schwarz_threshold, | ||
| ) |
There was a problem hiding this comment.
When schwarz_threshold > 0, a new SchwarzScreener (and full bounds matrix) is recomputed on every call. Since compute_schwarz_bounds() evaluates (ab|ab) via the full OS+HGP ERI routine for each shell pair, this setup cost can be substantial and may dominate for repeated ERI builds with the same basis. Consider caching the screener/bounds keyed by the basis (or exposing a way to reuse a precomputed screener) to avoid repeated O(n_shell^2) bound evaluations.
There was a problem hiding this comment.
@Copilot suggestion-6 Fixed!
Fixed! Added a module-level cache _schwarz_screener_cache keyed by (basis_ids, threshold). Now the SchwarzScreener and its bounds matrix are only computed once per unique basis+threshold combination — subsequent calls reuse the cached screener directly.
|
Thanks for the review @marco-2023! Some really good catches here. Going through all 6 fixes now:
Will push the fixes shortly! |
…threshold_sensitivity
|
@marco-2023 Fixed @Copilot suggestion 1, 4, 5, and 6! For @copilot suggestion 2 and 3 — the memory over-allocation in |
|
@msricher, @marco-2023 , @Copilot Thanks for the review! After careful analysis, the A-side axes in both integrals_vert and integrals_etransf must be sized m_max, not m_max_a. The ETR uses the recursion: [a|c+1] = ... − (ζ/η)[a+1|c] + ... This is a known characteristic of the OS/HGP algorithm: the VRR must "over-build" A-side angular momentum to provide the ETR with enough input. A triangular indexing scheme could reduce memory usage here, but that's a more involved restructuring. Happy to discuss further if I've misread the intent of the comment! |
marco-2023
left a comment
There was a problem hiding this comment.
Overall, the code is clean, nice work. I added some comments that I would like to get feedback on before this is merged.
| cont_four.angmom_components_cart, | ||
| cont_four.exps, | ||
| cont_four.coeffs, | ||
| ) |
There was a problem hiding this comment.
I thought the idea was adding the new class ElectronRepulsionIntegralImproved, which is going to ultimately substitute this one. Why the modifications to the original class?
There was a problem hiding this comment.
You're right @marco-2023 — this was an unintentional change. The original ElectronRepulsionIntegral class got modified in a commit that was meant to be a style/formatting fix, but it accidentally replaced the old _compute_two_elec_integrals implementation with compute_two_electron_integrals_os_hgp.
The original class should remain as-is (using the old algorithm as reference), and only ElectronRepulsionIntegralImproved should use the new OS+HGP pipeline. As a side effect, this also means the comparison tests in test_electron_repulsion.py (result_old vs result_new) will be meaningful again — right now both classes call the same function so the tests can't catch regressions.
Will revert this in the next push.
| self.n_computed = 0 | ||
|
|
||
| # Precompute Schwarz bounds | ||
| self.bounds = compute_schwarz_bounds(contractions, boys_func, compute_integral_func) |
There was a problem hiding this comment.
Here, compute_schwarz_bounds calls compute_schwarz_bound_shell_pair for each pair of shells. Then compute_schwarz_bound_shell_pair computes the integrals for each shell pair. Isn't this the same as computing all the integrals unscreened?
There was a problem hiding this comment.
I think I get this, you run a
There was a problem hiding this comment.
Exactly, @marco-2023 — compute_schwarz_bounds runs in O(n²) over shell pairs, computing (ab|ab) type integrals. This upfront cost then screens out negligible shell quartets from the O(n⁴) ERI computation, which is where the real savings come from for larger basis sets.
|
|
||
| """ | ||
|
|
||
| boys_func = PointChargeIntegral.boys_func |
There was a problem hiding this comment.
Why not use the new boys_function_standard?
There was a problem hiding this comment.
Makes sense. Updated ElectronRepulsionIntegralImproved.boys_func to use get_boys_function("coulomb") from boys_functions.py — no reason to borrow from PointChargeIntegral when we have a dedicated module for this.
| if cache_key not in _schwarz_screener_cache: | ||
| _schwarz_screener_cache[cache_key] = SchwarzScreener( | ||
| list(basis), | ||
| ElectronRepulsionIntegralImproved.boys_func, |
There was a problem hiding this comment.
I would import and use the boys' function from the new module. Maybe I am wrong, but adding the function as a member of a class to then pass the method unbounded as an argument seems overcomplicated. Your thoughts @San1357
There was a problem hiding this comment.
passing get_boys_function("coulomb") directly to SchwarzScreener is cleaner. Going through ElectronRepulsionIntegralImproved.boys_func was an unnecessary indirection, especially since get_boys_function is already imported. Fixed in the latest commit.
…lsionIntegralImproved
…ralImproved.boys_func
|
Hi @marco-2023, addressed all your comments:
All changes are in the latest push. Let me know if anything needs further revision. |
marco-2023
left a comment
There was a problem hiding this comment.
I think it is ready to merge
|
Thanks for the review and approval @marco-2023! Does this need another maintainer's sign-off before merging? |
Summary
Changes
feat/boys-functionsfeat/electron-repulsion-improvedfeat/VRR-two-electron-integralfeat/ETR-contractionfeat/HRR-complete-pipelinefeat/schwarz-screeningfeat/pyscf-tests-tutorialfeat/libcint-tolerance-fixTest Results
Related PRs
This PR closes #221, #216, #208 and is ready for final review and merge.
Please let me know (@PaulWAyers , @marco-2023 , @msricher ) if any further changes or refinements are needed. Thank you!