Far Away Grid Points
The following example computes the maximum distance between grid points, and it looks unusually large to me! I noticed it when testing Hirshfeld weights (instead of Becke weights) for integration. They go unnoticed when using Becke weights for integration. I haven't checked whether that is the case for HORTON2, but we need to add tests for grid points coordinates and weights against HORTON2.
import numpy as np
from scipy.spatial.distance import pdist
from grid.molgrid import MolGrid
from grid.becke import BeckeWeights
from grid.onedgrid import GaussChebyshev
from grid.rtransform import BeckeTF
atnums = np.array([8, 1, 1])
atcoords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.79523983], [1.69257101, 0.0, -0.59840635]])
oned = GaussChebyshev(30)
rgrid = BeckeTF(1.e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(atcoords, atnums, rgrid, 110, BeckeWeights())
# compute pairwise distance between grid points
pdist = pdist(grid.points)
print('min & max of pdist = ', np.min(pdist), np.max(pdist))
The output is: min & max of pdist = 0.0002980577 4377.9871590172
I think that the correct way to treat this may be to
- prune out points that have very small weights.
- prune out points at which some easy-to-evaluate benchmark function is sufficiently large (or small). It is difficult for me to conceive of a need to include grid points ~10 Angstrom from any atom.
Add rmax to OneDGrid functions
@FarnazH @PaulWAyers @tovrstra I think I need a little bit clearer implementation idea. So I am thinking add a rmax to OneDGrid default init, so any points after rmax will be trimmed, or a min_wt to truncate any points with less than certain given weight.
This argument can be default to None. Then we need add a rmax or min_wt args (can be **kwargs) to be add to each onedgrid functions. Since these changes are a little ~duplicate~ repetitive, so I wander are we making any change as mentioned in #118 so that each onedgrid function are turned into a class so some common behavior can be shared.
I see this would produce a lot of repetition in the onedgrid module, while the problem with far-away points is only appearing when creating molecular grids (because of the Hirshfeld partitioning). This suggests that it would be fine to support the rmax or min_wt only in higher-level classes such as AtomicGrid. One could even go further and say this is a Hirshfeld-specific issue and should therefore be resolved at that stage, just dropping points that are insensible from the Hirshfeld perspective. That would also make it easier to define a criterion for which point to discard safely.
That makes a lot of sense by putting at AtomicGrid level. This could be easily handled at high end. So without really type any code, I have two implementation ideas.
- when constructing
AtomicGridorMolecularGrid, we set the limit and drop extra radial points - Add a method(function) which return a subset of points from the complete
AtomicGrid. Since we have the pt_ind ( the indices of each shell), this can also be easily implemented. Therefore, we can have more flexibility at controlling dropping points without reallydropthem. (more memory usage to store the complete grid.)
Todo: check whether this problem exits in preset atomic grid (commented by Derrick, Farnaz, Paul)
Add documentation. Warn users that there are (very) far away points and points with zero weights in some cases.
This isn't easy to fix for molecular grids because the indexing refers to the number of points for each atom. So it isn't easy to fix in general.
We should give some examples of how this is dealt with in examples.