Skip to content

PR[9] (combine PR-1-PR-8) feat: combine weeks 1-9 implementations into single branch#245

Merged
marco-2023 merged 35 commits intotheochem:masterfrom
San1357:feat/woc-weeks1-9
Mar 16, 2026
Merged

PR[9] (combine PR-1-PR-8) feat: combine weeks 1-9 implementations into single branch#245
marco-2023 merged 35 commits intotheochem:masterfrom
San1357:feat/woc-weeks1-9

Conversation

@San1357
Copy link
Copy Markdown
Contributor

@San1357 San1357 commented Mar 11, 2026

Summary

  • Combines all week 1-9 implementations into a single branch for review
  • Includes Boys functions, OS+HGP algorithm (VRR+ETR+HRR), Schwarz screening, PySCF validation, and libcint tolerance fix

Changes

Week Branch Feature
1 feat/boys-functions Boys function APIs + attenuated potentials
2 feat/electron-repulsion-improved OS+HGP base rewrite
3 feat/VRR-two-electron-integral VRR recurrence (Eq. 65)
4 feat/ETR-contraction ETR + primitive contraction (Eq. 66)
5 feat/HRR-complete-pipeline HRR complete pipeline (Eq. 67) fix (#216)
6 feat/schwarz-screening Schwarz screening optimization
7 feat/pyscf-tests-tutorial PySCF comparison tests + Tutorial 3 notebook
8 feat/libcint-tolerance-fix libcint ERI tolerance fix (#208)

Test 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!

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 11, 2026

Summary (Final Completion Of WOC #221 )

  • This PR combines the complete implementation of the WoC project "Two-Electron Integrals via OS+HGP Algorithm" completed over Weeks 1–9 into a single branch, as requested.
  • Includes Boys functions, OS+HGP algorithm (VRR+ETR+HRR), Schwarz screening, PySCF validation, and libcint tolerance fix

Boys Functions (Week 1)

  • Implemented Boys function APIs for standard and attenuated potentials (erf, erfc, yukawa, etc.)
  • Added vectorized NumPy implementation with mpmath fallback for high accuracy

OS+HGP Two-Electron Integrals (Weeks 2–5)

  • Implemented complete Obara-Saika + Head-Gordon-Pople algorithm (VRR + ETR + HRR)
  • Supports s, p, d orbitals correctly (fixes [BUG] Electron repulsion energy #216)
  • Validated against PySCF for correctness

Schwarz Screening (Week 6)

  • Implemented Schwarz inequality-based screening to skip negligible shell quartets
  • Added shell-pair prescreening using Gaussian product theorem

PySCF Validation + Tutorial (Week 7)

  • Added comprehensive PySCF comparison tests
  • Added Tutorial 3 notebook demonstrating OS+HGP, Boys functions, and Schwarz screening

libcint Tolerance Fix (Week 8)

Request :

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

@San1357 San1357 requested a review from msricher March 11, 2026 07:09
Copy link
Copy Markdown
Collaborator

@msricher msricher left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for combining it all. Let's await a final review.

Comment on lines +19 to +24
try:
from pyscf import gto
HAS_PYSCF = True
except ImportError:
HAS_PYSCF = False

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@msricher Done! Added pyscf to [dev] dependencies in pyproject.toml and removed the conditional HAS_PYSCF imports from the test file.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Copy link
Copy Markdown
Contributor

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Suggested change
integrals_m = boys_func(orders, weighted_dist[None, :, :, :, :]) * prefactor
integrals_m = boys_func(orders, weighted_dist[None, :, :, :, :], rho=harm_mean) * prefactor

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Comment on lines +133 to +136
# 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")

Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines +251 to +254
# 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"
)
Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Comment on lines +112 to +116
# 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
Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Suggested change
# 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

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

Comment on lines +444 to +451
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,
)
Copy link

Copilot AI Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copilot uses AI. Check for mistakes.
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@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.

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 12, 2026

Thanks for the review @marco-2023! Some really good catches here.

Going through all 6 fixes now:

  1. The missing rho in boys_func is a real bug — erf/erfc would silently fail at runtime. Will pass rho=harm_mean.
  2. integrals_vert over-allocation makes sense to fix, A-side only needs m_max_a.
  3. Same for integrals_etransf, will resize the A-side axes.
  4. Thread-safety issue noted — will refactor Schwarz screener to pass via kwargs instead of class-level state.
  5. Yeah the tautological assertions are a bad miss on my part, will replace with a proper monotonicity check.
  6. Will add screener caching so bounds aren't recomputed on every call.

Will push the fixes shortly!

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 12, 2026

@marco-2023 Fixed @Copilot suggestion 1, 4, 5, and 6! For @copilot suggestion 2 and 3 — the memory over-allocation in integrals_vert and integrals_etransf is a valid point, but resizing the A-side axes to m_max_a would require changes deeper in the VRR/ETR logic. I'd rather not rush that and risk breaking things. Happy to tackle it as a follow-up once the current fixes are reviewed.

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 12, 2026

@msricher, @marco-2023 , @Copilot
Re: Comments 2 & 3 (VRR/ETR A-side over-allocation)

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] + ...
To build C angular momentum from scratch (c=0 → c=angmom_c+angmom_d), the first ETR step reads [a+1 | c=0] — i.e., the VRR output — where a ranges up to angmom_a + angmom_b + angmom_c + angmom_d. Capping the VRR at m_max_a = angmom_a + angmom_b + 1 would leave those higher-a slots as zero, silently corrupting results for any shell quartet with non-zero CD angular momentum (e.g. (ss|pp), (pp|dd)).

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!

Copy link
Copy Markdown
Contributor

@marco-2023 marco-2023 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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,
)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I get this, you run a $O(n^2)$ screening to screen the $O(n^4)$ integrals.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not use the new boys_function_standard?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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,
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 13, 2026

Hi @marco-2023, addressed all your comments:

  • Reverted ElectronRepulsionIntegral to use the original _compute_two_elec_integrals — the unintentional change is fixed.
  • Confirmed — compute_schwarz_bounds runs O(n²) to screen O(n⁴), which is the intended design.
  • Opened follow-up issue Add cache management for _schwarz_screener_cache (unbounded memory growth) #247 for cache management.
  • Updated ElectronRepulsionIntegralImproved.boys_func to use get_boys_function("coulomb") from boys_functions.py.
  • Passing get_boys_function("coulomb") directly to SchwarzScreener instead of going through the class.
  • Replaced inline hyp1f1 in PointChargeIntegral.boys_func with boys_function_standard.

All changes are in the latest push. Let me know if anything needs further revision.

@marco-2023 marco-2023 marked this pull request as ready for review March 14, 2026 01:34
@marco-2023 marco-2023 self-requested a review March 14, 2026 01:34
Copy link
Copy Markdown
Contributor

@marco-2023 marco-2023 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is ready to merge

@San1357
Copy link
Copy Markdown
Contributor Author

San1357 commented Mar 14, 2026

Thanks for the review and approval @marco-2023! Does this need another maintainer's sign-off before merging?

@marco-2023 marco-2023 merged commit 8c8f69c into theochem:master Mar 16, 2026
9 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

WoC: 2-electron integrals: better performance and more flexibility

4 participants