theochem / grid Goto Github PK

Hi,
I might just be missing something.
However I cannot import this package in my code. It tells me the module does not exist.

"""
ModuleNotFoundError: No module named 'grid'
"""

To install I did the following:

"""
git clone {ssh address}

in initial directory: pip install -e .
"""

I also added an init.py file in the initial directory, to see if that would be a fix, and it was not

I cannot pip install in the grid/src/grid directory because there is no setup.py

When I install it says this:

"""
Running setup.py develop for grid
Successfully installed grid-0.0.0
"""
Is the 0.0.0 an issue?

Here are the imports I have tried:
"""
import grid
from src import grid
from grid.src import grid
"""

All with the same result. The module cannot be found.

Thank you!

Kate

Change functions in onedgrid module to be subclasses of the OneDGrid grid, so that we can add pre_defined(atnum, preset) class method for initializing 1-D grids.

• A function which takes atomic numbers, atomic coordinates, and some description of grid and returns an instance of MolecularGrid class. (Note: here is where tv-*.txt points stored in needed; check #11)

I think the Molecular dipole moment calculation uses the wrong origin.

https://github.com/theochem/grid/blob/71876b912b50c1ab7343e74f86754255bcac93e9/src/grid/utils.py#L851C66-L851C66

The usual convention is to use the center of mass, where the most common isotope is used to define the mass of each element.

This doesn't matter for the case here, because if the monopole moment (charge) is zero the dipole moment is origin-independent.

There should be an easy-to-use functionality to get molecular grids from atomic grids, as well as an overall wrapper to take atomic positions and types and construct atomic grids directly, then molecular grids.

@Ali-Tehrani, based on https://github.com/theochem/grid/blob/master/src/grid/angular.py#L387, shouldn't we add a :math:4\pi to the integral given in https://github.com/theochem/grid/blob/master/src/grid/angular.py#L291? Let me know your thoughts.

From the current code:

• Store Lebedev quadrature (x, y, z) points and weights (on a unit sphere) in a npz file based on l_max.
• Have functions that load points/weights for a given l_max and return UnitSphereGrid.

UnitSphereGrid Class:

1. Initialize with (x, y, z) points and weights

Create a method inside MolGrid called get_atomic_property that takes in a atomic property as an array evaluated on the MolGrid.points and returns another atomic property of that array.

This method uses the aim_weights and MolGrid.indices in-order to get the atomic property.

This is the transformed cubic grid, spun off from Ali's watershed code. It should have the same methods (integration, differentiation, interpolation, etc.) as the Uniform Grid (Issue #7).

Running the tests generates a lot of warning messages, which should be addressed.

I think it may be worth adding spherical designs as an alternative angular quadrature. The rationales are that the fact that Lebedev-Laikov grids sometimes have negative weights gives problems (numerical and practical) and that for very smooth functions, spherical designs have some advantages (they seem competitive with Lebedev quadratures in general; I would speculate that this is because they do a "better job" of approximately integrating higher-degre spherical harmonics (beyond the integration order). Having uniform weights is also convenient, and ensures favorable behavior for, e.g., spherical averages and interpolation.

Other features:

• equal weights, all nonnegative.
• highly symmetric grids.
• competitive in point-count to Lebedev quadrature. Lebedev has $\frac{(L+1)^2}{3}$ points for order L; Spherical designs are $\frac{(L+1)^2}{2}$. See book chapter
• huge numbers of nodes available. (Up to degree 325!) See Wormersley. See also Neil Sloane's list
• It's possible to define nested quadratures.

I'm attaching a script and the resulting plot with errors on simple density integrals over the whole molecular grid, as function of the angular grid size. The radial grid is always the same (see script). Vertical lines correspond to the grids with negative weights. The three cases do not pop out as extremely problematic, so a warning may be too pedantic.

That said, the result for he2_ghost_psi4_1.0.molden is troubling. After some trial and error, this turned out to be related to the fact that the Bragg radius for Helium is set to Nan. This still results in a grid (with all finite weights) without raising an exception nor a warning. I'll open a separate issue for it.

@PaulWAyers Do you think the non-smooth trend of the error as function of grid size is fine? I would think so because one can always get a small error by chance.

lebcheck.zip

Originally posted by @tovrstra in #99 (comment)

The following example will give this error: TypeError: Smallest rgrid.points is negative, got -4.4999
This error is not very helpful, and probably the rgrid should have raised an error not MolGrid.

import numpy as np

from grid.molgrid import MolGrid
from grid.becke import BeckeWeights
from grid.onedgrid import HortonLinear
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 = HortonLinear(100)
rgrid = BeckeTF(1.e-4, 1.5).transform_1d_grid(oned)
grid = MolGrid.horton_molgrid(atcoords, atnums, rgrid, 110, BeckeWeights())

I tried to use MultiExpTF transformation combined with GaussChebyshevType2 radial grid by using the next code:

from grid.onedgrid import GaussChebyshevType2
from grid.rtransform import MultiExpTF

oned = GaussChebyshevType2(30) 
rgrid = MultiExpTF(1e-6, 1.5).transform_1d_grid(oned) 

And got this error:
ValueError Traceback (most recent call last)
in
3
4 oned = GaussChebyshevType2(30)
----> 5 rgrid = MultiExpTF(1e-6, 1.5).transform_1d_grid(oned)

~/grid/grid/src/grid/rtransform.py in transform_1d_grid(self, oned_grid)
99 if new_domain is not None:
100 new_domain = self.transform(np.array(oned_grid.domain))
--> 101 return OneDGrid(new_points, new_weights, new_domain)
102
103 def _convert_inf(self, array, replace_inf=1e16):

~/grid/grid/src/grid/basegrid.py in init(self, points, weights, domain)
256 if len(domain) != 2 or domain[0] > domain[1]:
257 raise ValueError(
--> 258 f"domain should be an ascending tuple of length 2. domain={domain}"
259 )
260 min_p = np.min(points)

ValueError: domain should be an ascending tuple of length 2. domain=[1.e+16 1.e-06]

The path submodule has been deprecated in importlib_resources >= 6.0.

To replace it, one should use the new files API (python/importlib_resources#80) or specify the version importlib_resources < 6.0.

We should have an example of:

• how to use the Local Grid functionality (to integrate local properties in a large molecule) which is based on KDTree
• how to use an atomic grid in a molecule to do local integration
• how to use periodic boundary conditions integrals.

AtomicGridClass:

1. Initialized with RadilGrid & UnitSphereGrid to generate points/weights.
2. Method to integrate.
3. Spherical averaging.

• Store pre-computed tv-*. in json/yml file.
• Have functions to load pre-computed points and return an instance of AtomicGridClass.

It came up yesterday in discussions whether there might be alternatives to spherical designs and Lebedev-style quadrature on the sphere. An alternative, with nonequal but always positive weights, is the maximum determinant stategy.

Up to a degree of 200, the points/weights are available here.

Incidentally, the same web page has one spherical design beyond what we have right now, for degree 325.

This isn't an especially high priority, but is a nice starter issue and something fun to play with.

6-dimension (double-molecular grid) for brute-force integration of Coulomb-like integrals. Depending on structure, promolecular tricks specific to this case may go here or may go somewhere else....

Using the "promolecule tricks" that Ali & Derrick developed for regional and whole-molecule integration of various functions.

At some point we should get rid of all the old HORTON commits. Cloning the repo just takes too long at the moment.

The 18 point degree 5 Lebedev rule is nonstandard.

Procedure:
It should be generated by the 6 points:
(+/- 1, 0, 0)
(0,+/-1,0)
(0,0,+/-1)
(+/-1/sqrt(2),+/-1/sqrt(2),0) (plus the 6 combinations)

Then evaluate the spherical harmonics up to degree 5 and determine the integration weights such that the errors are zero. (Solve linear equations). If the solution is not unique, add degree 6 functions as needed and do least squares.

At the moment the AIM weights are not mutiplied in, which makes the behavior of MolGrid different from its base class.

@Ali-Tehrani from issue #/141:
"""
The biggest issue I have is the radial transformations (LinearRTransform, ExpRTransform, PowerRTransform, HyperbolicRTransform) which were obtained from Horton/grid.py and how the documentation in python grid isn't correct as to what the domain/codomain are and they seem to be very specific to HortonGrid/UniformInteger class. For example, for LinearRTransform it transforms according to r(x) = rmin + x*(rmax - rmin)/(npoint - 1), where npoint is the number of points, rmin is the minimum of grid and rmax is maximum of grid.

The docs claim that it maps to [r_min, r_{max}] but this is only true if you're using HortonLinear/UniformInteger grid. This is also a problem in ExpRTransform and PowerRTransform, where it only maps to [r_{min}, r_{max}] only if you're using HortonLinear grid. One easy way of fixing this, is to replace npoint - 1 to taking the maximum over all grid points i.e. np.max(grid.points), this will hold true for any grid.

For HyperbolicRTransform, the mapping r(x) = ax/ (1 - bx) where a, b are positive scalars. It doesn't map to [0, \infty) but rather it can map everywhere. Additionally, there is the constraint b(npoints - 1) >= 1 that isn't explicit in the docs, I'm assuming it is also specific to HortonGrid where it centers the singularity at npoints - 1 and you only want the positive points on the left-side of singularity, see desmos/grapher, where N=npoints.
"""

Talking with @FarnazH , it would be good to add (to the base class) the ability to dump integration points and weights as a .npz file.

Going further (to give information about sub-grids, atomic grids, etc.) are "expert features" that grid supports, but don't require a high-level API.

@Ali-Tehrani can you add this functionality in short order?

1. Convert to Python
2. Compatible with ChemTools
3. Add Integration methods
4. Differentiation

OneDGrid Class, just a namedtuple with points and weights attributes. Two categories of functions should be present:

1. Generation. Functions to generate an instance of OneDGrid with new points and weights from scratch. These functions should have one mandatory argument, i.e. the first argument: the number of grid points. They could take additional arguments to control other aspects of the integration grids. All these functions should have names starting with generate_.

2. Transformation. Functions that take an existing OneDGrid instance as input (first argument) and optionally some parameters to control the transformation. The return value is a new instance of OneDGrid. All these functions should have names starting with transform_.

Paul's notes and fortran code contain many examples of both types of functions.

No need to make very long function names like generate_onedgrid_* and transform_onedgrid_*, because they already are collected in a module that provides a onedgrid namespace.

See: https://github.com/agronholm/sphinx-autodoc-typehints

The LinearInfiniteRTransform class retains memory about the first grid in which is used. This affects the behaviour of its transform_1d_grid method. When the same instance of the class is used to transform the points of two linear uniform grids with the same number of points but different ranges, the algorithm behaviour depends on the first grid in which is used.

Example:

from grid.onedgrid import UniformInteger    
from grid.rtransform import LinearInfiniteRTransform
import numpy as np

npoints=5
onedgrid_1=UniformInteger(npoints)
onedgrid_2=UniformInteger(npoints)

# create linear infinite R transforms (intemediate)
rtransform_1=LinearInfiniteRTransform(0., 3.)
# change the interval of the second grid to [0, 3]
onedgrid_2=rtransform_1.transform_1d_grid(onedgrid_2)

# create final linear infinite R trasform
rtransform_f=LinearInfiniteRTransform(0., 1.)
rtransform_f_virgin=LinearInfiniteRTransform(0., 1.)

# transform the two grids with the same R transform
rgrid_1=rtransform_f.transform_1d_grid(onedgrid_1)
rgrid_2=rtransform_f.transform_1d_grid(onedgrid_2)
rgrid_2_virgin=rtransform_f_virgin.transform_1d_grid(onedgrid_2)

# print the R transformed grids
print("Printing one-step R transformed grid points:")
print(rgrid_1.points)
print("Printing two-step R transformed grid points:")
print(rgrid_2.points)
print("Printing two-step R transformed grid points (virgin):")
print(rgrid_2_virgin.points)

The output is:

Printing one-step R transformed grid points:
[0.   0.25 0.5  0.75 1.  ]
Printing two-step R transformed grid points:
[0.     0.1875 0.375  0.5625 0.75  ]
Printing two-step R transformed grid points (virgin):
[0.   0.25 0.5  0.75 1.  ]

It can be seen how when the rtransform_f instance is used the second time, the scaling is not properly done. This is due to the method set_maximum_parameter_b of LinearInfiniteRTransform (line 761 of transform module). If a new instance (virgin) of the class is used, the method works.

Is this a feature of a bug?

Among the most efficient ways to integrate of a (hyper)cube are the lattice rules. There are several choices here, but probabilty the embedded lattice rules are the easiest to implement, and they have been tabulated.

For a (simpler) lead reference, I recommend the early paper from Sloan.

It would be good to use lattice rules as an alternative for integration over a cube or (after an affine transformation) an arbitrary parallelepiped. I'm much less sure that higher-dimensions are supported, but it is easy to support arbitrary-dimensional integration on a parallelepiped (beyond $d=3000$, which is more than enough!) with a lattice rule. Interpolation is also relatively straightforward, as it is a generalization of the current structure with specially selected vectors.

Implement SG1 as mentioned in issue #12.

In line 341, the weights equation is not correctly computed. Simply correction fo the value also leads to tests failure. So I think it's better to leave to @rayhe88. Could you take a look at this bug and have it fixed.

At the moment, pre-loading of lebedev grids is done inside the AtomicGrid class, meaning that it gets repeated for every atom, if I understand correctly (should be tested). It can also be implemented in the lebedev.py module, such that preloaded grids stay in memory until the process ends. (The memory overhead is very small.)

I had been trying to use solve_poisson_bvp to compute electronic repulsion between fragments, but apparently, the ode solver didn't converge and printed me the following error:

if I remove the origin it works but becomes extremely slow. I'm testing it using a water dimer and I attached the code with their respective files

ode_issue.zip

This is not urgent, just posting it here to keep track.

The following code can probably be converted to python and with vectorization tricks it should not be too slow: https://github.com/theochem/horton/blob/master/horton/moments.cpp#L59 It computes the real solid spherical harmonics, given Cartesian coordinates. This way one can avoid the series of conversions currently used: Cartesian -> phi, theta -> complex harmonics -> real harmonics.

• Have a pure Python implementation of Becke weights (just for the sake of having it)
• A function which takes all atomic numbers, all atomic coordinates, atomic grid points of one atom and it's an index to generate integration weights for that atom (i.e. something like Hirshfeld weights). These are called aim_weights.

MolecularGrid class:

1. Initialize with aim_weights & list of AtomicGrid instances. (Note: These weights are then multiplied by atomic weights to get the molecular integration weights.). It is probably better to store atomic number and coordinate of each atom in its AtomicGrid object.
2. Method for integration (inherited from base class).

Change default of rotate from false to a given seed value

This is something I intended to carry out shortly after #87 but I got carried away by other things.

The name "subgrid" seems to imply that it will always return a subset of the current grid, which is correct for aperiodic grids. In the case of a periodic grid, the result may contain portions from neighboring images, which could lead to more grid points in the subgrid than the parent grid. The name "localgrid" would make more sense because the algorithm takes a center and a radius as arguments.

Something to reproduce the dot_multi functionality in the old grid code. This is a big win computationally....

@tovrstra has some ideas about how this could be done. There would be a nice functionality in qcgrid for this but probably a bare-bones Python implementation for grid would be useful.

Setuptools version requirement in pyproject.toml is not compatible with Python 3.6:

requires = ["setuptools>=65.0", "setuptools_scm[toml]>=7.1.0"]
build-backend = "setuptools.build_meta"

Fixing this causes further problems when running pip install . and refers grid to as "UNKOWN":

Processing /home/pally/PythonProjects/grid
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
  Installing backend dependencies ... done
    Preparing wheel metadata ... done
Building wheels for collected packages: UNKNOWN
  Building wheel for UNKNOWN (PEP 517) ... done
  Created wheel for UNKNOWN: filename=UNKNOWN-0.0.5a6.post434-py3-none-any.whl size=3902 sha256=06580d1fd849e27839fde35b780de5cfb1b2c5ef1543a97f6aa1e4feccf50550
  Stored in directory: /tmp/pip-ephem-wheel-cache-d4brkxca/wheels/d1/98/24/14bed4afbf15d42cf2596818cc329de4a11f32aaac3c00d1f3
Successfully built UNKNOWN
Installing collected packages: UNKNOWN
  Attempting uninstall: UNKNOWN
    Found existing installation: UNKNOWN 0.0.5a6.post434
    Uninstalling UNKNOWN-0.0.5a6.post434:
      Successfully uninstalled UNKNOWN-0.0.5a6.post434
Successfully installed UNKNOWN-0.0.5a6.post434

The recommended solution is to update setuptools but this didn't work for me. It may be an issue with ubuntu 20.04 and python 3.6. The easiest solution is to require python>=3.7

The current boundary-value-problem for the Poisson solver is pretty robust. The initial value problem is a lot less robust.

It's a good first issue to try to make this work more efficiently. I will post some notes on how to do this.

For moments, we can use as tests a function (e.g. a Gaussian) that is centered at the origin and also one centered off the origin. The change in solid-harmonic coefficients when you shift the origin of a moment expansion is analytically computable.

• There is code in Julia that shifts the origin of a solid harmonic expansion

There are also some analytic formulas available.

• see page 66
• Ruedenberg found a classical result, which is simplified in Eq. (8) of this reference
• There is also a result in this reference
• A tabulation of the coefficients for translation/rotation is available in the supplementary material of this reference

Originally posted by @PaulWAyers in #13 (comment)

I was combining some onedgrid and transformation and I came accross with an error that I am not sure if it is a bug or if this is the intended behaviour of the code.

I combined the following oned grid and rtransform:

oned = Trapezoidal(22)
grid = BeckeTF(1e-4, 1.5).transform_1d_grid(oned) 

and I got the following error:

Traceback (most recent call last):
  File "test_from_predefined.py", line 36, in <module>
    rgrid = MultiExpTF(1e-6, 1.30).transform_1d_grid(oned) #BeckeTF(rmin, rmax)
  File "/home/leila/grid/grid/src/grid/rtransform.py", line 101, in transform_1d_grid
    return OneDGrid(new_points, new_weights, new_domain)
  File "/home/leila/grid/grid/src/grid/basegrid.py", line 272, in __init__
    f"point coordinates should not be above domain! {domain[1] < max_p}"
ValueError: point coordinates should not be above domain! False

I looked a little into this and this is what I understood. Trapezoidal one d grid has a domain=[-1,1], including -1 and 1, meaning the minimum and maximum points generated are -1 and 1. When the transformation is applied, -1 and 1 will be equal to the lower and upper domain values of the transform. In the case of BeckeTF domain is [r_min, inf], thus point with value -1 before transform = rmin and point with value 1 before transform = inf, but this value is trimmed so it will be equal to 1e16 instead of inf, after the transformation is applied. After obtaining this new values of the points, the code creates a new OneDGrid in basegrid.py and it goes through the several checks comparing the new points and the new domain. For the lowest point there's no problem but for the highest, 1e16 in this case, this line is what gives the error:

I am not sure about the purpose of the 1e-7 summation to the domain, but it is clear that by deleting the equal comparison this issue is solved, but as I said, I am not sure if this is a correct error detection. I checked all types of onedgrids and the two that gave this error are Trapezoidal and GaussChebyshevLobatto.

Let me know if something is unclear about the explanation or that if something that I wrote is wrong.

I think it would be helpful to have the functionality to define a uniform grid and export it as a cube file. Doing the reverse (going from a cube to a grid) could be useful too. The idea is that is easier to work with grid than with utilities such as Gaussian's cubman.

Also in this way, we provide an expedited way to visualize results e.g. surfaces that depend on electron densities (create a uniform grid and export to cube).

Get rid of get_fn file, and use importlib.resources to load test files.

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'll expand upon this later. I'm just putting it here not to forget it: during the workshop the emim_bf4 example temporarily needs about 6GB of RAM to construct the grid.

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.