WoC: Arbitrary-order Overlap Integrals (and evaluations enabled thereby)

[PaulWAyers]

Description

Add functionality to the GBasis library for arbitrary-order overlaps. One motivation of this is to support the evaluation of the intracule (the distribution function for the interelectronic distance) and the extracule (the distribution function for the center-of-mass of electron pairs). See these review articles.

📚 Package Description and Impact

GBasis is a pure-Python package for evaluating and analytically integrating Gaussian-type orbitals and their related quantities. The goal is to build a set of tools to the quantum chemistry community that are easily accessible and easy to use as to facilitate future scientific works.

👷 What will you do?

Your main focus will be to add functionality for the overlaps of arbitrarily many Gaussian basis functions. It is important to include screening, as otherwise the evaluations can be quite expensive. An important application, which should be supported with an appropriate API, is evaluation of the intracule and extracule. Another important application is the evaluation of near-arbitrary interparticle repulsions. This is possible because, using tricks popularized by Beylkin and Monzon, almost any interparticle repulsion can be written as a sum of Gaussians. This would include the (many-body) interactions associated with nucleons and other types of cofermions and cobosons. This is a place where some creativity can be exercised, as innovative new algorithms and improved fits are being consistently proposed. That said, these potential applications are beyond the scope of this specific issue, though they are relatively easy extensions thereto.

This functionality also allows us to compute arbitrary positive integer powers of the electron density and related density matrices analytically. (Typical approaches rely, instead, on numerical integration that become increasingly ill-conditioned for high powers.) This is sometimes useful when describing (very) multicenter bonding.

3- and 4-Gaussian overlaps from PySCF can be used for testing. 1-Gaussian overlaps and 2-Gaussian overlaps are built-in already and can be used for testing. An algorithm similar to what is needed here is already included in IOData (for two Gaussians). See also the overlap integrals in GBasis.

🏁 Expected Outcomes

1. Implement the overlap for an arbitrary number of (contracted, Cartesian or spherical) Gaussian basis functions.
2. Add support for screening, with appropriate default parameters.
3. Add functionality for evaluating the intracule and extracule.
4. Utilities for computing positive integer powers of the electron density and relevant density matrices.
5. Support for (near) arbitrary interparticle repulsion integrals. This uses the aforementioned trick of Beylkin and Monzon. A similar trick is used in MADNESS, which can be used for inspiration.
6. Write comprehensive tests and documentation for all new functionality.
7. Write tutorial Jupyter notebooks that show how to use the new functionality.
8. Follow our guidelines for contributors.

Required skills	Python, OOP
Preferred skills	Be comfortable with math and numerical algorithms. Experience with scientific programming and quantum chemistry can help
Project size	350 hours, Large
Difficulty	Medium to Hard 🥵

🙋 Mentors

Marco Martínez-González	mmg870630_at_gmail_dot_com	@marco-2023
Farnaz Heidar-Zadeh	farnaz_dot_heidarzadeh_at_queensu_dot_ca	@FarnazH
Esteban Vöhringer-Martinez	estebanvohringer_at_qcmmlab_dot_com	@evohringer

📝 Notes

An algorithm (two algorithms, in fact) for doing this are included in the attached notes,
which were developed in the context of the intracule and extracule. The attached notes are old and have Fortran-style pseudocode (they are that old). [To be clear, I would not suggest implementing the direct algorithms for the intracule/extracule. It is better to pass through the multi-Gaussian overlaps.]

Note that multi-Gaussian overlaps are extremely sparse, so a sparse tensor structure should be used. The main operation one needs are tensor contractions (usually pairing indices to one or more copies of the 1- or 2-electron density matrix), which will determine the type of structure one uses to store the tensor. One may need to use pytorch functionality for this.

If we had the arbitrary-order overlaps coded, I would argue that we should support some additional evaluations, notably:

1. intracule on a grid.
2. extracule on a grid.
3. integrals of powers of the electron density.
4. (maybe) integrals of moments of powers of the electron density. These integrals show up in some of the moment density functional development work that Parr and Nagy did a couple decades ago. That specific topic is not very active right now, though, so this could definitely be deferred to the (far) future.
5. This supersedes GSoC 2025: Arbitrary-order Overlap Integrals (and evaluations enabled thereby) #183

[MitanshuKumawat]

Hi! I just came across this repository and I am very interested in contributing to this WoC project.

I have experience with Python, numerical methods, and Density Functional Theory from my research work. This project aligns very well with my interest in Gaussian basis methods and electronic structure theory.

Could you please guide me on the steps to get started and what the initial milestones for this project would look like? Also, is there any formal process to participate in WoC (similar to GSoC), or should I simply begin by contributing through issues/PRs?

I would be happy to begin contributing immediately.
Thank you!

[aamrindersingh]

Hi @PaulWAyers

I am very interested in this issue. I have been reading through the GBasis code and the notes for a few days.
Found it very interesting and educational to go deep into the math. Opening a PR in a moment with a basic N-center overlap implementation using Obara-Saika recursion. Looking forward to feedback to start working on further implementation.

Thanks

[kir943]

Hello mentors (@marco-2023 @FarnazH @evohringer),
I’m Kiran Jadhav, a selected contributor for this project under Winter of Code 5.0.
My proposal for QC-DEVS was accepted, and I’m excited to start working on the
“Arbitrary-order Overlap Integrals” task.

I’ve started setting up the GBasis development environment and am going through
the existing 1- and 2-Gaussian overlap implementations to understand the current
code structure.

For Week 1, my plan is to finalize the algorithmic approach for extending overlaps
to arbitrary N-Gaussians and implement a basic working version (without screening),
validated on small test cases.

Please let me know if there’s any part of the codebase or notes you’d recommend
prioritizing first.
Looking forward to working with you. Thank you!

[PaulWAyers]

Hi Kiran! It's great to have you on board. I'll add you to the requisite repositories and teams.

[kir943]

Hi @PaulWAyers,
I’ve been going through the existing 1- and 2-Gaussian overlap implementations and reviewing the background notes shared in the discussions channel.
I’m currently focusing on finalizing the algorithmic structure for extending overlaps to arbitrary N-Gaussians (without screening initially), and identifying which parts of the current overlap code can be reused or generalized.
Next step is to sketch a minimal working version and validate it on small test cases.
Please let me know if there’s a specific file or implementation detail you’d recommend paying special attention to at this stage.

[msricher]

I'm not the foremost expert on the code structure, but I think I would first focus on:

• the abstract base class for contraction-associated arrays, in base.py, since you will need to implement a sparse (2N)-dimensional array structure for an N-contraction overlap array.
• the user-facing function interface, in integrals/overlap.py, since you'll be writing something similar to this.
• The multipole moment code, integrals/_moment_int.py and whatever it depends on, since I think you'll be using this code and/or writing generalizations of it / writing similar code to it.
• the 1- and 2- index screening functions, in screening.py, are also probably worth looking at since you'll need to screen the overlaps for the (2N)-dimensional computation.

[kir943]

Thanks for the pointers, this helps a lot.
I’ll start by focusing on the contraction-associated array base class in base.py, particularly the dimensionality assumptions and how a sparse (2N)-dimensional overlap array could be integrated. I’ll also review the interface in integrals/overlap.py to keep the N-Gaussian extension consistent with the existing API.
I’ll use the multipole moment code in integrals/_moment_int.py as a reference for the algorithmic structure, and revisit screening in screening.py after a minimal working version is in place.
I’ll follow up once I have initial results or concrete questions.

[PaulWAyers]

I think the better reference for algorithmic structure is the pseudocode in the notes and the analogous implementation in IOData.

https://github.com/theochem/iodata/blob/e6a606e985a4a0010d64982c526efa1b95ac09b6/iodata/overlap.py

The compute_overlap_1d functional can perhaps be (re)used. The other moment integrals in gbasis could be coded in this way (and might be better) but I think they are right now done with Obara-Saika.