theochem / grid

Python library for numerical integration, interpolation, and differentiation on (molecular) grids.
https://grid.qcdevs.org/
GNU Lesser General Public License v3.0
43 stars 17 forks source link

BeckeWeights.generate_weights should raise an exception when it tries to use Nan radii #100

Closed tovrstra closed 2 years ago

tovrstra commented 4 years ago

The following example will cause the becke partitioning to use Nan radii. The resulting grid gives large integration errors but not Nans:

import numpy as np
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import GaussChebyshev
from grid.rtransform import BeckeTF
becke = BeckeWeights(order=3)
oned = GaussChebyshev(150)
rgrid = BeckeTF(1e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(np.array([[0, 0, -1], [0, 0, 1]]), np.array([2, 2]), rgrid, 74, becke)
print(grid.weights.min())
print(grid.weights.max())

Output:

/home/toon/miniconda3/envs/horton3/lib/python3.7/site-packages/grid/becke.py:82: RuntimeWarning: invalid value encountered in greater
  alpha[alpha > cutoff] = cutoff
/home/toon/miniconda3/envs/horton3/lib/python3.7/site-packages/grid/becke.py:83: RuntimeWarning: invalid value encountered in less
  alpha[alpha < -cutoff] = -cutoff
-121784734604746.08
109395391879381.38

The resulting grid is not usable, so it would be safer to raise an exception when trying to construct becke weights with atoms for which the radii ar not defined.

PaulWAyers commented 3 years ago

Psi4's radii can be found at: https://github.com/psi4/psi4/blob/b4a272f74b73df8cf0edc4096081a250c495d55e/psi4/src/psi4/libfock/cubature.cc#L111

They just use H's radius for He. The use 3.3 a.u. whenever it isn't available.

tczorro commented 3 years ago

add an error message to nan after the rollback if exists