Pruned grids

[leila-pujal]

This PR is related to issue #112 for the implementation of pruned grids SG type and the ones reported in doi:10.1063/1.5049435 (Ochsenfeld). The following code/files have been added:

• method pruned_grid in MolGrid class: The main difference between this and methods from_size and from_preset is that we don’t pass a OneDGrid type object. This is because the number of points of the radial grid is specified by the prune grid and it must be created inside the function. This function accepts any correct combination of radial grid and transformation.

• Folloing npz files for the pruned grids: All have keys = atnum_rad or atnum_npt. atnum_rad for sg_1 correspond to the alpha parameter as they prune the grid based on regions defined by scaling the braag radii by each alpha. prune_grid_sg_1.npz prune_grid_sg_3.npz prune_grid_sg_2.npz prune_grid_sg_0.npz prune_grid_g1.npz prune_grid_g2.npz prune_grid_g3.npz prune_grid_g4.npz prune_grid_g5.npz prune_grid_g6.npz prune_grid_g7.npz

• test_pruned_molgrid in test_molgrid().

Code to test the function

from iodata import load_one
from grid.molgrid import MolGrid
from grid.becke import BeckeWeights

import numpy as np

mol=load_one("H2.fchk")

becke = BeckeWeights(order=3)

pruned_grid = MolGrid.pruned_grid(mol.atnums, mol.atcoords,'UniformInteger',
                                  'ExpRTransform,1e-5, 2e1', 'sg_2', becke, rotate=False, store=True)

test file H2.tar.gz

Obtained values for sg pruned grids for C atom: sg_0 --> 1414 - This is different from what it is reported in here because we don't support the exact number of points for Lebedev grid for 18 (I believe it is the only one missing) sg_1 --> 1928 - This is reported to give around 2000 points for each atom here sg_2 --> 7790 - This is the same as reported here sg_3 --> 17674 - This is the same as reported here

Obtained values for Ochsenfeld pruned grids for C atom: g1 --> 2828 g2 --> 5418 g3 --> 10202 g4 --> 16344 g5 --> 22654 g6 --> 43002 g7 --> 72002

These are not the same as the ones reported in the paper. Explanation here

[FarnazH]

@leila-pujal, see AtomGrid.from_pruned https://github.com/theochem/grid/blob/master/src/grid/atomgrid.py#L151 which should be extended to also include the pruned schemes you are adding (a quick suggestion would be to make radius and sectors_r optional, with default=None, and allow one to assign sectors_r='sg_2')... Adding the MolGrid.from_pruned is on the agenda, which should just use AtomGrid.from_pruned. Hope this helps. I will be going over the problem you reported in reproducing the literature results soon...

[leila-pujal]

Thanks, @FarnazH to answer the pull request. Just to check if I understood your message; so you want me to have the functionality of these pruned grids in AtomGrid class instead of MolGrid class? From your message, I understand it should be extended in the from_pruned method? If this is what you meant in your message the only problem I see is that to use AtomGrid.from_pruned() method you need to pass as an argument a OneDGrid object and for that, you need to specify the radial points. Almost all of the pruned grids in this PR use a specific number of radial points. At this point I don't see clearly how to put the code in the from_pruned method. Probably I am missing somthing but it has been a while since I have used Grid package and you have been working on the code much more. So if you could answer this doubt I will move the code to from_pruned in AtomGrid class. Also if I am moving the code, do you want to keep the function MolGrid.from_pruned or do you want to delete it? I am not sure from your message if you meant to keep it or not.

[FarnazH]

n_rad = get_rgrid_size("sg_2", atnum=6)
oned = GaussChebyshev(n_rad)
rgrid = BeckeTF(1.e-4, 1.5).transform_1d_grid(oned)
g = AtomGrid.from_preset(rgrid, atnum=1, preset="sg_2") 
# check rgrid.size against
# rgrid can be None

[leila-pujal]

I just committed the updates to include the changes we talked about in the past Grid meeting. These are:

• Modify AtomGrid.from_preset to be able to create the new pruned grids
• Add function get_rgrid_size in atomgrid module to obtain the radial points that should be used to create an Atomic/Molgrid of these new types of pruned grids.

Now you can use AtomGrid.from_preset to create a pruned atomic grid of the type sg_X or gX. Also you can use Molgrid.from_preset to create a Molgrid with pruned atomic grids. I did not use from_pruned method in AtomGrid because for that function you don’t pass a string but specify the radial/angular sectors. You can use the code below to test and it should give the same number of points as the ones reported at the beginning of the PR. Now for pruned grids sg_X you can still pass only one instance of atomic grids because for this type all atoms have the same number of points. On the other hand for gX type the user have to create all the atomic grids for all the atoms (see example below in the code).

Test code for get_rgrid_size

n_rad = get_rgrid_size('sg_0', atnums=list(range(1,18)))
print('n_rad)

Test code for AtomGrid.from_preset

n_rad = get_rgrid_size('sg_3', atnums=6)[0]
oned = GaussChebyshev(n_rad)
rgrid = BeckeRTransform(1.e-5, 1.5).transform_1d_grid(oned)
g = AtomGrid.from_preset(rgrid, atnum=6, preset="sg_3")
print(g.points.shape)

Test code for Molgrid.from_preset for sg_2

n_rad = get_rgrid_size('sg_3', atnums=6)[0]
oned = GaussChebyshev(n_rad)
rgrid = BeckeRTransform(1.e-5, 1.5).transform_1d_grid(oned)
becke = BeckeWeights(order=3)
mg = MolGrid.from_preset(mol.atnums, mol.atcoords, rgrid, 'sg_3', becke, store=True)

Test code for Molgrid.from_preset for g2

becke = BeckeWeights(order=3)
oneds = [GaussChebyshev(rad) for rad in get_rgrid_size('g2', atnums=mol.atnums)]
rgrids = [BeckeRTransform(1.e-5, 1.5).transform_1d_grid(oned) for oned in oneds]
mg = MolGrid.from_preset(mol.atnums, mol.atcoords, rgrids, 'g2', becke, store=True)

The last thing is some considerations about sg_1 grid I write below but they can discuss on the next grid meeting:

• Sg_1 pruned grid is the only one of the sg_X grids that is based on distance to divide the points between the angular sectors. In their paper the authors use a Euler-Maclaurin grid radial grid of 50 points and these are divided based on the different distance ranges to obtain the number of points of each angular sector. However they only provide this distance sectors up to Z=18 (Ar) and after that they use n_radial points = 50 and angular =194 (I.e not pruned). I mention this because for this new function get_rgrid_size to obtain the radial points for each of these pruned grids I set it up to give 50 radial points, although you could use any number of points for sg_1 and Z < 18 and it won’t give an error.
• For sg_1 I used GaussChevishev + BeckeRTransformation and now I get around 4000 points and based on this, I should get around 2000. I know we are not using their atomic grid but I just wanted to note it. -To create the npz for sg_1 (I attached the scripts in the next message below) the distance sectors used to divide the radial points are based on bragg radii but using the function from grid.utils import get_cov_radii some of the entries get a nan value (atnum = [2,10,18,36,54,86,87]) so I couldn’t generate an entry in the sg_1 for these.

[leila-pujal]

Here are the scripts I used to generate the npz files for the new pruned grids scripts_generate_npz.tar.gz

[leila-pujal]

Final note: I haven't deleted the original Molgrid.from_pruned. I thought I will do it once we have discussed the new changes.

[PaulWAyers]

I thought we'd fixed the radius issue, but perhaps not. It's discussed here.

@tczorro do you remember what we did on the grid? It would be reasonable to use the PySCF radii, or Psi4, or from some other open-source code. Results should not be that sensitive to it, but it's an (undergrad level) project to see whether or not different radii give better or worse results.

[tczorro]

@PaulWAyers I am not too aware of the question discussed here. But for the radii reference, it happened to the noble gas atom where their radii is not defined/available. How we circumvent was to use the radii of the atom with 1 less atomic number (or 2, 3, 4 less until it's available)

[PaulWAyers]

@tczorro this is a nice pragmatic solution and I think it provides what @leila-pujal needs to solve her NaN problem.

[leila-pujal]

Thanks a lot for your comments @PaulWAyers and @tczorro. I remembered we discussed this some time ago and I thought a solution was suggested but because I still got the nans I was not sure if it was implemented. Now I remember what @tczorro wrote in his comment, that we should get the one with 1 less atomic number. I will use this and update the npz files. Thanks.

[Ali-Tehrani]

I've added a few commits regarding the documentation with additional help from @leila-pujal. Everything looks fine with me, I've added the option to use symmetric spherical t-design, even though the pruned grids weren't intended for them, but don't think it would be a big issue.