Add extra line to unit tests of lebedev grids

[tovrstra]

@tczorro This is actually more an issue than a PR. It seems there is an error in the datafiles containing the Lebedev grids. I've added a corresponding line to the unit tests, which fails (appropriately). The problem is that some of the Lebedev quadrature weights in the data files are negative while they should be all positive. (Corresponding grids gave rather weird results when I was trying to use them). This is probably a glitch in the script used to generate these data files? Could you add this script to the data directory, such that we can take a closer look? Without it, it is difficult to see what went wrong and how to fix it. Thanks.

(I've used .min() so the failing assert shows the most negative weight.)

[tczorro]

Hi Toon, I checked other Lebedev integral code (quadpy, numgrid). It seems the weight for lebedev scheme with degree(13)- points(74) does have negative weights for several points. Haven't done a complete scan for all degree of freedoms but i guess many of them may encounter negative weights.

[tovrstra]

It seems a bit dangerous to have negative quadrature weights but I also found them in HORTON2 for N = 14, 74, 230 & 266. All the other ones have positive coefficients. I'm wondering if this is actually correct. Looking for other implementations, 14 seems to have positive weights:

• https://github.com/Rufflewind/lebedev_laikov/blob/master/src/Lebedev-Laikov.c#L679
• https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html
• http://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

I think I may have taken a wrong or outdated reference when implementing them in HORTON2.

[tczorro]

It seems a bit dangerous to have negative quadrature weights but I also found them in HORTON2 for N = 14, 74, 230 & 266. All the other ones have positive coefficients. I'm wondering if this is actually correct. Looking for other implementations, 14 seems to have positive weights:

• https://github.com/Rufflewind/lebedev_laikov/blob/master/src/Lebedev-Laikov.c#L679
• https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html
• http://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

I think I may have taken a wrong or outdated reference when implementing them in HORTON2.

I agree it's dangerous, but it probably designed to cancel out something if the grid is taken as a whole.
In the GRID implementation, I found 74, 230, and 266 with some negative weights. I think it may worth giving a warning when user choose these schemes especially partially picking some integral points

[tovrstra]

I see, indeed 14 has positive weights. I also checked the literature reference for N=74 (10.1016/0041-5553(75)90133-0) and the paper has a remark "the quadrature has negative weight" for this one. I guess that is the 1975 equivalent of a Python warning. :)

I'll check if I get systematically worse results for 74, 230, and 266. If so, a warning might indeed be in order. If there is by chance a cusp of a neighbouring atom at a grid point with negative weight, I would expect trouble, because this is not easily compensated for by the other points of the angular grid.

[tczorro]

I'll check if I get systematically worse results for 74, 230, and 266. If so, a warning might indeed be in order. If there is by chance a cusp of a neighbouring atom at a grid point with negative weight, I would expect trouble, because this is not easily compensated for by the other points of the angular grid.

Thanks a lot. Please let me know the results. I will try to make a warning draft just incase. But i feel one point we will need it anyway.

[PaulWAyers]

For multidimensional integration rules, there is no guarantee that the quadrature weights are positive. The best-behaving formulas are positive-weighted (mitigates round-off error), but that's not universal. So it doesn't surprise me that there are some negative weights, as long as they are few and small. It would be good to clearly document them and to print warnings too.

Often the negative weights can be eliminated by using larger (non-minimal) numbers of points, or using a different symmetry (e.g., icosahedral vs. octahedral vs. tetrahedral) symmetry.

The real test of the grids would be to confirm that spherical harmonics up to a specified order are integrated correctly.

[tovrstra]

I'm attaching a script and the resulting plot with errors on simple density integrals over the whole molecular grid, as function of the angular grid size. The radial grid is always the same (see script). Vertical lines correspond to the grids with negative weights. The three cases do not pop out as extremely problematic, so a warning may be too pedantic.

That said, the result for he2_ghost_psi4_1.0.molden is troubling. After some trial and error, this turned out to be related to the fact that the Bragg radius for Helium is set to Nan. This still results in a grid (with all finite weights) without raising an exception nor a warning. I'll open a separate issue for it.

@PaulWAyers Do you think the non-smooth trend of the error as function of grid size is fine? I would think so because one can always get a small error by chance.

lebcheck.zip

[tovrstra]

It turns out I disabled warnings when constructing the grid, hiding the problem. It should raise an exception though, because the result is insensible when Nan radii are used.

[PaulWAyers]

The real solution is to put in a Bragg-Slater radius for Helium or move to a more modern atomic radius set. PySCF uses 1.40 for He, which seems HUGE but probably still is better than nothing. Alternatively we could just use Clementi-Roetti or some other form of radii that are universally defined as the default. PySCF supports several choices.

The up-and-downs is not unexpected; once the error is small you often get an off-trend result when the approximate integral goes from being too low to too high (or vice versa) as one increases grid size, as a "magic grid" that comes close to the crossing point can give a fortuitously small error. Just using my eyes, I think that if you averaged the curves over all the molecules you'd get something with a more consistent (and smoother) decreasing trend.

[FarnazH]

@tczorro can you please mention the negative weights when describing how lebedev points/weights were generated? See: https://github.com/theochem/grid/blob/master/src/grid/lebedev.py#L22

Any other checks or precautions @tovrstra & @PaulWAyers?

[PaulWAyers]

This was really just a discussion/issue. Negative weights occur with Lebedev grids sometimes (not often). We should print a warning and indicate this in the documentation when it arises.

We have raised a new issue related to the atomic radius issue. #158