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.

[Ali-Tehrani]

Okay, I've gotten into a rabbit hole about this, sorry for the delay. I doubled checked, the initial icosahedron and first subdivision have equal weights, and they're precisely the spherical design. The next subdivision has four positive, unique weights, then the subdivision afterwards had negative weights, then subdivision afterwards had positive weights.

These are the degrees that they correspond to, starting with the icosahedron, respectively, 2, 5, 11, 24, and 49. The number of points (12, 42, 162, 642, 2562) is usually double of Lebedev and Spherical design grids.

I was wrong about the weights being the area, and also that the area of each triangle is equal (so there is a defect). I've implemented it using Girards thoerem, and played around with it. That led me to do a deep delve into literature. Particularly, this book chapter was a nice read

To develop a subdivision scheme for spherical quadrature. Any spherical quadrature scheme is evaluated based on:

1. Number of points
2. Maximum degree of spherical harmonic it can integrate accurately
3. Minimal $\sum |w_i|$ ( See Theorem 10 in the book chapter)
4. Positive Weights (if true, then $\sum |w_i| = \sum w_i = 4\pi$ and minimizes the Theorem 10 Error)
5. Symmetry of the grid
6. Whether you want good, fast interpolation (e.g. max determinant grids) or nested properties.

The nice thing about this subdivision is that the weights are usually positive, which is known to be a problem when enforcing polynomial accuracy (Section 4). The (area of triangle) centroid method I've talked about earlier always has positive weights, but they don't enforce polynomial accuracy.

Other Nested Grids: Another nested grid is based on having equal-area points (Section 6.1), also known as HEALPix grids. They are very easy to implement. One could enforce polynomial accuracy, but it is a question whether the weights are positive and the number of points.
Recent math paper that shows the existence of nested spherical design grids.

Ideas: I think it would be nice to implement a general angular grid, where the user gives a random set of points (Section 5, scattered points), performs Delaunay triangulation, calculates the weights using triangle area, or it implements the calculated weights using this procedure, and tells the user if the weights are negative or absolute sum is large.
Also I wonder if it is better instead of adding points on the edges, to add points on the center of the triangles, this would reduce the number of points, and consequently smaller angular degree.

[PaulWAyers]

I think if you add a point at the center of the triangle you end up with very non-equilateral triangulations. That's probably not ideal. Nesting always compromises accuracy.

For our applications it is really important that they are nested. So I think just the weights that Ali mentions for the grids are fine; we could just store them and use them. It's unfortunate that there are negative weights in some cases, but that was true for some Lebedev grids too....probably it is inevitable.

I also like the HEALPix grid which is, as @Ali-Tehrani notes, very easy to implement see equation two in the reference Ali gave and beautifully nested when $N_{side}$ is incremented by factors of 2 (or any other integer, I think, but 2 is the obvious choice). The growth rate goes up by factors of 4, which is relatively slow. It's equi-area so one can just use equal weights, as in a spherical design.

Perhaps we just store the weights for the icosahedral division up to 2562 points and have HEALPix as a secondary option with slower growth? That would give us two options.

I'm not averse to doing something with random weights too, though without nesting it's not directly useful for the applications I have in mind.

One nice thing about the Theorem 10 that @Ali-Tehrani cites is that it gives a very good way to control what to do when one can satisfy the "interpolatory" property up to $L$ but this does not fully determine the equation: one uses SVD to minimize the norm of the weights (subject to the constraint that the rule is interpolatory). This theorem justifies that this choice, which was intuitive to me, is rigorous. Note that one can imagine doing an "interpolatory" HEALpix grid too.

Perhaps the right workflow is to write generic code to construct the "best" interpolatory quadrature formula for a given set of points. Perhaps @Ali-Tehrani already has this from the study on the icosahedron?

[Ali-Tehrani]

Sorry never mind, unfortunately, HEALPix isn't nested, I thought it was from naively looking at the pictures. But they are very good it seems like at adaptive grid refinement. I guess the only other nested grid that comes in-mind are doing Clenshaw-Curtis tensor product grids. These do cluster points at the poles, but I think if you rotate them to align with the bond axis, that might actually be a good thing.

Also forgot to mention, the other benefit of positive weights is when the integrands are noisy. This is my guess, as to why the geoscience community cares so much. Also, only 6 points out of 642 that has negative weights, so it isn't too bad.

I think solving for the nodes and weights, storing them as npz, and adding a method=icosahedral flag to AngularGrid would be an option. It might be a little weird with how AngularGrid takes in degree, making it difficult to specify the next subdivision. Also the points should be ordered, so that the new points are appended at the end from the previous subdivision, i.e. from index= $N_{i-1}, \cdots, N_i - 1$, where $N_i$ is the number of points in the $i$ th subdivision.

The scatter points/random points I was thinking of having applications outside quantum chemistry. Or, more broadly towards geosciences, i.e. any of the methods within the two book chapters.

I would suggest for the generic code to get quadrature on scatter points with enforced polynomial accuracy is to follow the same approach as Oscar but:

1. Do bounded Linear Least-squares scipy.optimize.lsq_linear(spherical_harmonics, b, bounds=(1e-5, 4 * np.pi)).
   • The alternative option to do non-negative least-squares but this gives way too many zeros.
2. Check if the cost function is very low < 1e-12, if not, accept you'll have negative weights and do normal linear least-squares.
3. Report if any of the weights are zero, negative or $\sum |w_i|$ is large.

[PaulWAyers]

OK, so we have the icosahedral grid, which is awesome. @Ali-Tehrani 's idea for structuring the points makes sense.

For the random points, I understand @Ali-Tehrani 's protocol. I would only add that the lower bound should decrease as the number of grid points increases, so maybe something like $1e-3*N_{grid}^{-1}$ is a better lower bound.

When one can't do a full fit with degree $L$, one has the choice of adding negative weights or reducing to grid $L-1$ and then trying to minimize the L2 norm of the weights. I can see that being better when there are substantial negative weights and/or the weights get large.

P.S. @Ali-Tehrani I also fell for the HEALPix based on the description. It looks so nested (even the formulas) but I didn't confirm. But from the paper you have above, it seems like a nesting scheme is possible and the code from that paper may implement it. Perhaps an interesting direction, as it gives a nice way to do adaptive grids in specific directions, which is (generically) hard.