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

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 | 

Start by rerunning the provided Colab example using AtomGrid, the listed grid classes, radial transforms, and promol.density. Compare the integration results with the expected electron count, then determine which grid/transform combinations produce the extreme errors and whether the issue is in numerical integration or parameter use. Done means the discrepancy is explained and appropriate combinations or a corrective change are identified.