Density Integration Inaccuracy with AtomDB Densities and Specific Grid/Transform Combinations

[izxle]

When using electron densities from the AtomDB repository within the Grid library, the integrated density does not correctly sum to the expected number of electrons for certain combinations of grids and radial transforms. This discrepancy suggests a potential issue with the numerical integration methods or their interaction with AtomDB density functions.

Any insights into why specific grid/transform combinations might lead to discrepancies would be helpful. Additionally, guidance on choosing appropriate grid and radial transformations parameters to improve accuracy would be greatly appreciated.

Thank you for looking into this issue!

Expected Behavior

The total integrated electron density should match the expected number of electrons for the given atomic system.

Observed Behavior

The integrated density deviates from the expected value depending on the grid and transformation settings.

Steps to Reproduce

1. Load an atomic density from AtomDB.

2. Choose different grid configurations and radial transformations in Grid.

3. Perform numerical integration of the density.

4. Compare the integrated value to the expected electron count.

Example Code

The code was executed on a Google Colab notebook.

# !pip install -q qc-AtomDB git+https://github.com/theochem/grid.git

import numpy as np
import pandas as pd

from grid import AtomGrid
from grid.onedgrid import GaussLegendre, GaussChebyshevType2, GaussChebyshevLobatto, RectangleRuleSineEndPoints, TanhSinh, Simpson, MidPoint, ClenshawCurtis, FejerFirst, FejerSecond, Trapezoidal
from grid.rtransform import BeckeRTransform, LinearFiniteRTransform, MultiExpRTransform, KnowlesRTransform, HandyRTransform, HandyModRTransform
from atomdb import make_promolecule


n_electrons = 10
n_points = 151
r_min = 1e-3
r_max = 4
R = 1.5

promol = make_promolecule(atnums=[n_electrons], coords=[0, 0, 0], dataset="gaussian")

grid_types = {
    "Gauss-Chebyshev Type 2": GaussChebyshevType2,
    "GaussLegendre": GaussLegendre,
    "Gauss-Chebyshev-Lobatto": GaussChebyshevLobatto,
    "Rectangle-rule sine": RectangleRuleSineEndPoints,
    "tanh-sinh": TanhSinh,
    "Simpson": Simpson,
    "Midpoint": MidPoint,
    "Clenshaw-Curtis": ClenshawCurtis,
    "Fejer First": FejerFirst,
    "Fejer second": FejerSecond,
    "Trapezoidal": Trapezoidal
}

rtransforms = {
    "Becke": BeckeRTransform,
    "Linear Finite": LinearFiniteRTransform,
    "Multi Exp": MultiExpRTransform,
    "Knowles": KnowlesRTransform,
    "Handy": HandyRTransform,
    "Handy Mod": HandyModRTransform
}

results = []
for radial_name, grid_class in grid_types.items():
    for rtransform_name, rtransform_class in rtransforms.items():
        print(f"Performing integration for grid {radial_name} and {rtransform_name}")
        oned = grid_class(n_points)
        
        if rtransform_name == "Linear Finite":
            rgrid = rtransform_class(r_min, r_max).transform_1d_grid(oned)
        elif rtransform_name in ["Handy", "Handy Mod"]:
            rgrid = rtransform_class(r_min, r_max, 1).transform_1d_grid(oned)
        elif rtransform_name == "Knowles":
            rgrid = rtransform_class(r_min, R, 1).transform_1d_grid(oned)
        else:
            rgrid = rtransform_class(r_min, R=R).transform_1d_grid(oned)

        grid = AtomGrid(rgrid)

        dens = promol.density(grid.points)
        n_integral = grid.integrate(dens)

        results.append({
            "Radial Grid": radial_name,
            "RTransform": rtransform_name,
            "N integral": n_integral,
            "Error": abs(n_integral - n_electrons),
        })

df = pd.DataFrame(results).sort_values('Error')

# This part is for ipython display
n_figures = 12
formatted_display = df.style.format({
    'N integral': f'{{:15.8e}}',
    'Error': f'{{:10.2e}}'
})
display(formatted_display)

Here's the result in the form of a table:

| Radial Grid             | RTransform    | N integral      | Error           |
+-------------------------+---------------+-----------------+-----------------+
| Midpoint                | Knowles       |  9.99992458e+00 |  7.54192219e-05 | 
| Rectangle-rule sine     | Knowles       |  1.00001142e+01 |  1.14233352e-04 | 
| Rectangle-rule sine     | Becke         |  1.00001200e+01 |  1.20011617e-04 | 
| Gauss-Chebyshev-Lobatto | Knowles       |  9.99986593e+00 |  1.34067443e-04 | 
| Fejer second            | Knowles       |  9.99986579e+00 |  1.34212830e-04 | 
| Gauss-Chebyshev Type 2  | Knowles       |  9.99986562e+00 |  1.34383769e-04 | 
| GaussLegendre           | Knowles       |  9.99986538e+00 |  1.34618943e-04 | 
| Fejer First             | Knowles       |  9.99986524e+00 |  1.34763446e-04 | 
| Rectangle-rule sine     | Linear Finite |  1.00004079e+01 |  4.07924476e-04 | 
| Midpoint                | Becke         |  1.00006409e+01 |  6.40887330e-04 | 
| Simpson                 | Linear Finite |  9.99800847e+00 |  1.99153299e-03 | 
| Midpoint                | Linear Finite |  9.99788443e+00 |  2.11556703e-03 | 
| Midpoint                | Handy Mod     |  9.99660407e+00 |  3.39592724e-03 | 
| tanh-sinh               | Linear Finite |  9.99660400e+00 |  3.39600319e-03 | 
| Simpson                 | Handy Mod     |  9.99656963e+00 |  3.43036560e-03 | 
| Rectangle-rule sine     | Handy Mod     |  9.99656786e+00 |  3.43213554e-03 | 
| Fejer second            | Linear Finite |  9.99656081e+00 |  3.43918714e-03 | 
| Gauss-Chebyshev Type 2  | Linear Finite |  9.99655997e+00 |  3.44003093e-03 | 
| Clenshaw-Curtis         | Linear Finite |  9.99655690e+00 |  3.44309656e-03 | 
| Fejer First             | Handy Mod     |  9.99655611e+00 |  3.44388530e-03 | 
| Gauss-Chebyshev-Lobatto | Linear Finite |  9.99655606e+00 |  3.44394323e-03 | 
| GaussLegendre           | Handy Mod     |  9.99655571e+00 |  3.44429136e-03 | 
| Clenshaw-Curtis         | Handy Mod     |  9.99655561e+00 |  3.44438795e-03 | 
| Fejer second            | Handy Mod     |  9.99655514e+00 |  3.44486073e-03 | 
| Gauss-Chebyshev-Lobatto | Handy Mod     |  9.99655433e+00 |  3.44567404e-03 | 
| Gauss-Chebyshev Type 2  | Handy Mod     |  9.99655391e+00 |  3.44608711e-03 | 
| GaussLegendre           | Linear Finite |  9.99655222e+00 |  3.44777624e-03 | 
| Fejer First             | Linear Finite |  9.99654769e+00 |  3.45230902e-03 | 
| tanh-sinh               | Handy Mod     |  9.99654454e+00 |  3.45546377e-03 | 
| Trapezoidal             | Handy Mod     |  9.99645736e+00 |  3.54264319e-03 | 
| Trapezoidal             | Linear Finite |  9.99476213e+00 |  5.23786575e-03 | 
| Rectangle-rule sine     | Handy         |  1.00063460e+01 |  6.34599704e-03 | 
| Midpoint                | Handy         |  1.02793950e+01 |  2.79394987e-01 | 
| Fejer First             | Multi Exp     | -9.99986524e+00 |  1.99998652e+01 | 
| GaussLegendre           | Multi Exp     | -9.99986538e+00 |  1.99998654e+01 | 
| Gauss-Chebyshev Type 2  | Multi Exp     | -9.99986562e+00 |  1.99998656e+01 | 
| Fejer second            | Multi Exp     | -9.99986579e+00 |  1.99998658e+01 | 
| Midpoint                | Multi Exp     | -9.99992458e+00 |  1.99999246e+01 | 
| Rectangle-rule sine     | Multi Exp     | -1.00001142e+01 |  2.00001142e+01 | 
| Gauss-Chebyshev-Lobatto | Becke         |  6.21826194e+05 |  6.21816194e+05 | 
| Gauss-Chebyshev Type 2  | Becke         |  7.29031985e+05 |  7.29021985e+05 | 
| Fejer second            | Becke         |  8.59371365e+05 |  8.59361365e+05 | 
| GaussLegendre           | Becke         |  2.22528357e+07 |  2.22528257e+07 | 
| Gauss-Chebyshev-Lobatto | Handy         |  2.24174169e+08 |  2.24174159e+08 | 
| Gauss-Chebyshev Type 2  | Handy         |  2.62805886e+08 |  2.62805876e+08 | 
| Fejer second            | Handy         |  3.09791995e+08 |  3.09791985e+08 | 
| Fejer First             | Becke         |  4.83005375e+09 |  4.83005374e+09 | 
| GaussLegendre           | Handy         |  8.01374189e+09 |  8.01374188e+09 | 
| Fejer First             | Handy         |  1.73795354e+12 |  1.73795354e+12 | 
| tanh-sinh               | Knowles       |  7.34918175e+60 |  7.34918175e+60 | 
| Clenshaw-Curtis         | Knowles       |  2.11908098e+72 |  2.11908098e+72 | 
| Clenshaw-Curtis         | Handy         |  2.11908098e+72 |  2.11908098e+72 | 
| Clenshaw-Curtis         | Becke         |  2.11908098e+72 |  2.11908098e+72 | 
| tanh-sinh               | Becke         |  8.55393215e+72 |  8.55393215e+72 | 
| Simpson                 | Handy         |  2.11898680e+74 |  2.11898680e+74 | 
| Simpson                 | Becke         |  2.11898680e+74 |  2.11898680e+74 | 
| Simpson                 | Knowles       |  2.11898680e+74 |  2.11898680e+74 | 
| Trapezoidal             | Handy         |  3.17848020e+74 |  3.17848020e+74 | 
| Trapezoidal             | Becke         |  3.17848020e+74 |  3.17848020e+74 | 
| Trapezoidal             | Knowles       |  3.17848020e+74 |  3.17848020e+74 | 
| tanh-sinh               | Handy         |  3.07594237e+75 |  3.07594237e+75 | 
| Simpson                 | Multi Exp     |            -inf |             inf | 
| Clenshaw-Curtis         | Multi Exp     |            -inf |             inf | 
| Trapezoidal             | Multi Exp     |            -inf |             inf | 
| Gauss-Chebyshev-Lobatto | Multi Exp     |             nan |             nan | 
| tanh-sinh               | Multi Exp     |             nan |             nan | 

[Ali-Tehrani]

This is very cool and interesting. Thank you for submitting this issue. Just couple of quick thoughts I have, that may or may not help with your question

• One way to determine if a quadrature is accurate should depend on what kinds of integrand the quadrature it's accurate for! Unfortunately, we don't have a list of this.
  Reference [1] has a well-written methodology and reasoning of choosing some of these transforms.

• One method is to plot the output of the RTransform from using the OneDGrid and see how the points are distributed.

• I would recommend placing a smaller rmin from 1e-3 to 1e-6/1e-8 a.u. So that it captures near the nucleus better.

• One important thing to note is that your promolecular example is symmetric, so certain grids may be advantageous over others (due to cancellation of odd-terms within its Taylor Series). I'm guessing changing your AtomicGrid to use angular grid with "symmetric spherical t-design" would give lower errors because the weights are all positive. This can be done by setting `AtomGrid(method="spherical"). I assume spherical t-design have mirror symmetry for this to work.

• I would recommend PowerRTransform based on prior work done by analyzing diatomics integration accuracy from QC-Devs, you can use it naturally with AtomGrid.from_preset([n_electrons]).

• It's not surprising that Knowles gives the highest accuracy based on literature [1, 2].

• Multi-Exp should have the best level of accuracy (since it's meant for multi exponential integrand) [1]. But it's definitely not working and giving opposite, negative values is very surprising, indicating there is definitely a bug in the grid. Having a quick look shows there is a bug in calculating its derivatives

[1] Gill, Peter MW, and Siu‐Hung Chien. "Radial quadrature for multiexponential integrands." Journal of computational chemistry 24.6 (2003): 732-740.
[2] Kakhiani, Khatuna, Kakha Tsereteli, and Paata Tsereteli. "A program to generate a basis set adaptive radial quadrature grid for density functional theory." Computer Physics Communications 180.2 (2009): 256-268.

[izxle]

Thank you for the detailed feedback! I’ve implemented some of the suggestions and have new results to share.

Changes made based on your recommendations:

• Adjusted rmin to 1e-8 to better capture density near the nucleus.
• Set the angular grid method to "spherical" to improve symmetry handling.
• Used PowerRTransform via AtomGrid.from_preset(n_electrons, preset='fine'). I added it to the table as "Fine preset, Power".

I appreciate the insight regarding quadrature accuracy and the recommendation to visualize point distribution using OneDGrid with RTransform. I’ll explore that further to understand how different grid choices affect accuracy.

Additionally, the issue with Multi-Exp returning negative values is interesting. If there's anything I can do to help resolve the bug, let me know!

Thanks again for the helpful references and suggestions. Let me know if there are any other adjustments you'd recommend.

Updated results:

| Radial Grid             | RTransform    | N integral      | Error           |
+-------------------------+---------------+-----------------+-----------------+
| Gauss-Chebyshev-Lobatto | Knowles       |  9.99986873e+00 |  1.31268658e-04 | 
| Gauss-Chebyshev Type 2  | Knowles       |  9.99986783e+00 |  1.32171682e-04 | 
| Fejer second            | Knowles       |  9.99986782e+00 |  1.32180905e-04 | 
| Fejer First             | Knowles       |  9.99986772e+00 |  1.32276244e-04 | 
| GaussLegendre           | Knowles       |  9.99986758e+00 |  1.32420527e-04 | 
| Midpoint                | Knowles       |  9.99986705e+00 |  1.32952274e-04 | 
| Rectangle-rule sine     | Knowles       |  1.00001467e+01 |  1.46698469e-04 | 
| Rectangle-rule sine     | Becke         |  1.00001527e+01 |  1.52724442e-04 | 
| Fine preset             | Power         |  9.99982368e+00 |  1.76324840e-04 | 
| Midpoint                | Becke         |  1.00005830e+01 |  5.82956102e-04 | 
| Rectangle-rule sine     | Linear Finite |  1.00007957e+01 |  7.95743170e-04 | 
| Simpson                 | Linear Finite |  9.99809418e+00 |  1.90582147e-03 | 
| Midpoint                | Linear Finite |  9.99742737e+00 |  2.57262562e-03 | 
| tanh-sinh               | Linear Finite |  9.99660646e+00 |  3.39354204e-03 | 
| Rectangle-rule sine     | Handy Mod     |  9.99659508e+00 |  3.40491579e-03 | 
| Clenshaw-Curtis         | Linear Finite |  9.99656135e+00 |  3.43864648e-03 | 
| Fejer second            | Linear Finite |  9.99656105e+00 |  3.43895252e-03 | 
| Gauss-Chebyshev-Lobatto | Linear Finite |  9.99656098e+00 |  3.43902465e-03 | 
| Gauss-Chebyshev Type 2  | Linear Finite |  9.99656066e+00 |  3.43934198e-03 | 
| Simpson                 | Handy Mod     |  9.99656001e+00 |  3.43999317e-03 | 
| Midpoint                | Handy Mod     |  9.99655879e+00 |  3.44121375e-03 | 
| Fejer First             | Handy Mod     |  9.99655857e+00 |  3.44142889e-03 | 
| GaussLegendre           | Handy Mod     |  9.99655832e+00 |  3.44168268e-03 | 
| Fejer second            | Handy Mod     |  9.99655772e+00 |  3.44227564e-03 | 
| Clenshaw-Curtis         | Handy Mod     |  9.99655768e+00 |  3.44232239e-03 | 
| Gauss-Chebyshev Type 2  | Handy Mod     |  9.99655666e+00 |  3.44334032e-03 | 
| Gauss-Chebyshev-Lobatto | Handy Mod     |  9.99655656e+00 |  3.44344287e-03 | 
| Trapezoidal             | Handy Mod     |  9.99655429e+00 |  3.44571416e-03 | 
| GaussLegendre           | Linear Finite |  9.99655268e+00 |  3.44732334e-03 | 
| tanh-sinh               | Handy Mod     |  9.99655087e+00 |  3.44913268e-03 | 
| Fejer First             | Linear Finite |  9.99655006e+00 |  3.44993500e-03 | 
| Trapezoidal             | Linear Finite |  9.99566479e+00 |  4.33521467e-03 | 
| Rectangle-rule sine     | Handy         |  1.00067207e+01 |  6.72069887e-03 | 
| Midpoint                | Handy         |  1.02789326e+01 |  2.78932566e-01 | 
| Midpoint                | Multi Exp     | -9.99986705e+00 |  1.99998670e+01 | 
| GaussLegendre           | Multi Exp     | -9.99986758e+00 |  1.99998676e+01 | 
| Fejer First             | Multi Exp     | -9.99986772e+00 |  1.99998677e+01 | 
| Fejer second            | Multi Exp     | -9.99986782e+00 |  1.99998678e+01 | 
| Gauss-Chebyshev Type 2  | Multi Exp     | -9.99986783e+00 |  1.99998678e+01 | 
| Rectangle-rule sine     | Multi Exp     | -1.00001467e+01 |  2.00001467e+01 | 
| Gauss-Chebyshev-Lobatto | Becke         |  6.21825966e+05 |  6.21815966e+05 | 
| Gauss-Chebyshev Type 2  | Becke         |  7.29031725e+05 |  7.29021725e+05 | 
| Fejer second            | Becke         |  8.59371059e+05 |  8.59361059e+05 | 
| GaussLegendre           | Becke         |  2.22528310e+07 |  2.22528210e+07 | 
| Gauss-Chebyshev-Lobatto | Handy         |  2.24174139e+08 |  2.24174129e+08 | 
| Gauss-Chebyshev Type 2  | Handy         |  2.62805851e+08 |  2.62805841e+08 | 
| Fejer second            | Handy         |  3.09791954e+08 |  3.09791944e+08 | 
| Fejer First             | Becke         |  4.83005332e+09 |  4.83005331e+09 | 
| GaussLegendre           | Handy         |  8.01374126e+09 |  8.01374125e+09 | 
| Fejer First             | Handy         |  1.73795348e+12 |  1.73795348e+12 | 
| tanh-sinh               | Knowles       |  7.34918175e+60 |  7.34918175e+60 | 
| Clenshaw-Curtis         | Knowles       |  2.11908098e+72 |  2.11908098e+72 | 
| Clenshaw-Curtis         | Handy         |  2.11908098e+72 |  2.11908098e+72 | 
| Clenshaw-Curtis         | Becke         |  2.11908098e+72 |  2.11908098e+72 | 
| tanh-sinh               | Becke         |  8.55393215e+72 |  8.55393215e+72 | 
| Simpson                 | Handy         |  2.11898680e+74 |  2.11898680e+74 | 
| Simpson                 | Becke         |  2.11898680e+74 |  2.11898680e+74 | 
| Simpson                 | Knowles       |  2.11898680e+74 |  2.11898680e+74 | 
| Trapezoidal             | Becke         |  3.17848020e+74 |  3.17848020e+74 | 
| Trapezoidal             | Knowles       |  3.17848020e+74 |  3.17848020e+74 | 
| Trapezoidal             | Handy         |  3.17848020e+74 |  3.17848020e+74 | 
| tanh-sinh               | Handy         |  3.07594237e+75 |  3.07594237e+75 | 
| Simpson                 | Multi Exp     |            -inf |             inf | 
| Trapezoidal             | Multi Exp     |            -inf |             inf | 
| Clenshaw-Curtis         | Multi Exp     |            -inf |             inf | 
| Gauss-Chebyshev-Lobatto | Multi Exp     |             nan |             nan | 
| tanh-sinh               | Multi Exp     |             nan |             nan |