Adaptive Quadrature Methods

[PaulWAyers]

Description

Add functionality to the Grid library for adaptive quadrature. Grid has quite efficient methods for (molecular and periodic) molecular integration. However, in cases where the accuracy of a grid is inadequate, one typically has to recompute, from scratch, with a larger grid. It would be better to have nested grids, where one can (ideally locally) improve the accuracy of the grid, without discarding the points that one has already computed.

📚 Package Description and Impact

Grid is a pure Python library for numerical integration, interpolation and differentiation of interest for the quantum chemistry community.

👷 What will you do?

Use a traditional Becke-style molecular grid, with nested radial (relatively easy) and angular (not too hard; see #258) points. This is effectively a way to adaptively prune a grid.

🏁 Expected Outcomes

1. Implement Issue Nested Angular Grids (for adaptive quadrature) #258 .
2. Implement adaptive radial grids.
3. Use steps 3 and 4 to construct adaptive Becke-style grids.
4. Write comprehensive tests and documentation for all new functionality.
5. Write tutorial Jupyter notebooks that show how to use the new functionality.

Required skills	Python, OOP
Preferred skills	Be comfortable with math and numerical algorithms. Experience with scientific programming can help
Project size	175 hours, Medium
Difficulty	Medium 🤔

🙋 Mentors

Marco Martínez-González	mmg870630_at_gmail_dot_com	@marco-2023
Paul Ayers	ayers_at_mcmaster_dot_ca	@PaulWAyers

[azhar-ikhtiarudin]

Hi @marco-2023 @PaulWAyers

I am Azhar Ikhtiarudin, a recent graduate in Engineering Physics from the Bandung Institute of Technology (ITB), Indonesia. I have experience in machine learning, scientific computation, and quantum computing (particularly for quantum chemistry simulation) through several research internship at ITB, ANU, Indonesian Research Agency, and QOSF. I am eager to learn more about the classical computational aspects of quantum chemistry, and I believe contributing to QC-Devs (hopefully through GSoC 2025) would be a great fit for me to gain this experience.

I am interested in working on this project, as well as another issue provided here, but I am currently unsure which would be the best fit for me and would appreciate further guidance. May I know the best way to communicate further? Can I directly send an email expressing my interest?

Thank you.

Best regards,
Azhar

[PaulWAyers]

You're welcome to e-mail me. I would probably suggest that, if you are interested in quantum chemistry, that this may be a better fit. The other one is more related to model Hamiltonians, and might be a more difficult ramp-up. This is a pretty self-contained project that, nonetheless, would greatly benefit the package.

[azhar-ikhtiarudin]

Thank you for your response! I appreciate your time and insights. I will follow up via email shortly.

[M-Haseeb01]

Dear Mentors @marco-2023 and @PaulWAyers,

I'm a CS student with a solid foundation in Python, OOP, and numerical computing, reinforced by my coursework. I'm very interested in contributing to your project through GSOC 2025 and would appreciate any guidance on where to begin or suggestions for beginner-friendly tasks.

Best regards,
Haseeb

[PaulWAyers]

The "starter project" we're recommending for people is to write tutorial notebooks/documentation for features. Those are "beginner-friendly" tasks.

For the project itself, it's good to look at methods like #258, which is now a pull request.

[Unnati-2061]

Hi @PaulWAyers ,
I've been working on generating nested grids and implemented an icosahedral grid using a subdivision approach. I explored different ways to generate these grids and put together a script that successfully creates subdivisions. While refining the approach, I ran into a few questions:

• The angular grids (48, 192, 768) look reasonable, but (192, 768, 3072) might be too large. Would it be better to adjust the growth rate?
• For integrating spherical harmonics correctly, should we optimize weights at each level manually, or is there a standard approach you'd recommend?
• Regarding radial grids, do we have a preferred way to match them with these angular grids, or should I experiment with different scaling techniques?

I’ll keep refining the weight calculations and test different configurations. Looking forward to your thoughts!

[PaulWAyers]

@SophiaLi20 the grown rate is sort of what we are stuck with; that is the price we pay for nesting. One can slow the growth rate but then the points become either very inhomogeneous or the nesting is compromised. There is slower growth if one seeds with a tetrahedral grid.

The weights do change for each grid. There are (lots of) different ways to choose them. The goal is always to integrate exactly as many spherical harmonics as one can for a given grid size, and then try to reduce the error for spherical harmonics that one cannot compute exactly. Key criteria are minimal norms for the integration weights, nonnegative integration rates, and low errors for higher-order spherical harmonics (that cannot be integrated exactly).

Radial grids, in this approach, are constructed independently. In general, one wishes for the radial-grid-error and angular-grid-error to be similar.

[Unnati-2061]

Hi @PaulWAyers,
Got it, that makes a lot of sense! I see how the growth rate is a trade-off—keeping nesting intact while avoiding inhomogeneous points is tricky. I’ll take a closer look at using a tetrahedral-seeded grid to see if it helps smooth things out.
For the weights, I’ll try out different approaches to balance minimal norms, nonnegative integration rates, and reducing errors for higher-order spherical harmonics. If there are any references or past methods you’d recommend, I’d love to check them out!
On the radial grids, I get that the goal is to keep the radial and angular errors in sync. I’ll run some tests to see how things balance out and whether tweaking the radial scaling improves accuracy. Do you have any specific benchmarks or test cases you’d suggest for checking this?
Looking forward to your thoughts!

[azhar-ikhtiarudin]

Dear @PaulWAyers @marco-2023, I have studied and tested the implementation of the nested angular grid based on the current pull requests and understand the concepts. However, I have some questions regarding the next steps for implementation.

1. Radial Grid Error:

   • In the case of refining the radial grid, should the error estimation involve comparing computed values to a reference true value, such as the total number of electrons obtained using another method? (because I think this approach may not be generalizable).
   • Alternatively, is there a general approach or framework to estimate the radial grid error that doesn't require reference values? For instance, based on my initial research, could self-consistency checks or convergence tests across refinements be sufficient for error analysis?

2. Combined Grids:

   • When constructing a combined adaptive grid, is it preferable to first develop an adaptive radial grid and then optimize the angular grid for each radial point? Or should we explore a more advanced approach, such as simultaneously comparing and balancing errors between the radial and angular grids for improved accuracy and performance?

Thank you so much

[PaulWAyers]

One can estimate the radial grid error based on the difference between tighter and looser radial grids. I think it is sensible to first converge the radial grid for a "reasonable" angular grid and then refine the angular grid. Technically the two should be adapted simultaneously, but that seems less practical to me.

[azhar-ikhtiarudin]

I see, I got the idea. Thank you very much

[PaulWAyers]

This needs cleanup.

[alok-108]

Hello @PaulWAyers, @marco-2023,

I am interested in contributing to the Adaptive Quadrature Methods project. I've been exploring the codebase and seen why adaptive methods are critical--fixed grids can easily miss features in sharp functions unless we use an inefficiently high number of points. I wrote a small test with a sharp Gaussian ($e^{-1000x^2}$) and saw that even 100 Gauss-Legendre points gave ~4.7e-6 error, confirming the need for local refinement.

I have a few clarifying questions regarding the implementation:
(1) Expected Interface: Do you envision a high-level function like grid.adaptive_integrate(f, tol=1e-6) that returns the integral value, or a way to refine the Grid object itself (e.g., grid.refine(f, tol))?
(2) Radial vs. Angular: Should the initial focus be on 1D recursive refinement for the radial grid?
(3) Library Reuse: Is it allowed to use scipy.integrate.quad for error estimation/base integration, or should the logic be implemented from scratch?
(4) Integration Method: Should adaptive refinement be a separate function, or should the integrate method of an existing grid class accept a tol parameter and refine internally?

Looking forward to your guidance!