Created Nested_Angluar_Grid.ipynb

[izxle]

This PR adds an adaptive nested angular grid using icosphere subdivisions. The method iterates over radial grid points, refining the angular grid until the integral converges within a set threshold.

We welcome any feedback and suggestions to improve this implementation!

[PaulWAyers]

@marco-2023 remind me to discuss this with you in person. It will then be very helpful if you can help review this.

[Ali-Tehrani]

Thank you for the interesting work, I've gone over this, it's rather interesting (because it breaks some of my assumptions) and the results look great.
Some general comments:

• I would recommend to run the code in jupyter notebook (especially the plots) for others to see the results without having to run your code. You can clear output of some cells that isn't needed.

• Your result is rather interesting because it breaks the general assumption that less angular points are needed near the center (but could be explained since the integrand is also highest near the center). (However note that the molecular grids in grid are meant for integrating molecular system rather than atomic).

• I gone over the subdivision scheme of the icosahedron, it makes sense to me. Couple of references for this are 1 and 2. This picture made the most sense to me:
  .

• When integrating any electron-density, I would recommend choosing a reference grid such that it integrates to the correct number of electrons. For example, changing oned_grid = Trapezoidal(npoints=100) to have 250 points or more should give you three-decimal place accuracy of 9.999sh versus 9.0sh.

• The methodology of converging using two consequent angular grids per radial point makes sense to me. However, the goal is to reduce the number of times the integrand is called, so you will need a way to track the new points, and only recompute the integrand for those new points.

• Finally, I would recommend more complicated molecules, e.g. diatomics at different geometries as well!

The important comments I have is the process of choosing the quadrature weights:

• This subdivision methodology is also the nested scheme for spherical quadrature used in the book "Spherical Harmonics and Approximations on the Unit Sphere" Section 5.2. If you choose the weights to be the area of the polygonal component, you get a O(1/N_l) convergence and for each division scheme the error should reduce by four (for twice-continuous function). They also give three other quadrature schemes R1-R3 (Section 5.2.2), that give a amazing: O(1/N^2) O(1/N^3) convergence.
  What's also nice about this subdivision scheme is that the new set of points is also invariant under icosahedron group.

• I suspect that you don't need to do least-square solution on the quadrature weights. In-fact, if you plot the weights, they are all constant and most likely the area of polygon component. I haven't put much thought into this, but I would heavily look into this, and try to derive a formula.
  - This is nice because spherical harmonics of odd degree are odd functions and so a constant weight alongside a symmetric grid points results in the integral of these spherical harmonics to always be zero.

Another comment that is relevant but probably won't help you

• When constructing any sort of Lebedev quadrature on the sphere, you don't just do a least-squares over all spherical harmonics. You remove the spherical-harmonics that are not invariant under some rotational group (e.g. icosahedron group). See Section 2.2.2. in 3 . This produces a over-determined solution whose results is not unique (although in your case this isn't a problem but just important to keep in mind). Also to keep in mind, there are a O(N^3) algorithm for constructing sphere quadrature that are invariant under icosahedral group "Rotationally invariant quadratures for the sphere".

[PaulWAyers]

Thanks @Ali-Tehrani and @izxle !

In-fact, if you plot the weights, they are all constant and most likely the area of polygon component. I haven't put much thought into this, but I would heavily look into this, and try to derive a formula.

This is a really interesting observation. The idea of using the polygonal area is ideal, and does make it easier to also envision (eventually) angular refinement, but I note that the equivalence of R1, R2, and R3 in terms of convergence rate heavily depends on the "cancellation" which isn't going to be true for nonuniform triangulation.

In all these grids, there is a "defect." Specifically, most vertices have 6 neighboring triangles, but this is only 5 for the original vertices of the icosahedron (and 4 for nested octahedral grids and 3 for nested tetrahedral grids). So as one moves away from the vertices, there is "less defect" and I expect the polygonal areas to be more constant. It is possible, I guess, that this all comes out in the wash: the fact that the apical point for example has only five triangular areas may be compensated by other factors.

I am a big perplexed (but I didn't read very hard) that the formulas seem to use the midpoint of the edges (R1) or interior points (R2, R3), which we don't have access to. I would have guessed, then, that the appropriate weights were constant everywhere (unlikely due to the aforementioned defect) or equal to the surface area of the sphere that is closer to this point than any other (which requires a little nontrivial differential geometry, at least to me). It can probably be done (relatively) easily using the (generalized) Girard's theorem. I think the workflow is to construct the angles from the (geodesic) edges using a formula like this. It is nice if we don't have to solve equations for the weights, however.

[PaulWAyers]

@Ali-Tehrani do you have any thoughts on this?

[Ali-Tehrani]

I am still curious about this, and can have an answer next week! I think I'm going to explore it computationally given the code is there. I don't think there is defect in the symmetry of the grid, but need to think/read more into it.