theochem / grid Goto Github PK
View Code? Open in Web Editor NEWPython library for numerical integration, interpolation, and differentiation on (molecular) grids.
Home Page: https://grid.qcdevs.org/
License: GNU Lesser General Public License v3.0
Python library for numerical integration, interpolation, and differentiation on (molecular) grids.
Home Page: https://grid.qcdevs.org/
License: GNU Lesser General Public License v3.0
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:
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.
We should have an example of:
Change default of rotate
from false
to a given seed value
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.
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())
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....
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:
Line 266 in d3d0590
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.
@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
.
"""
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
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).
Using the "promolecule tricks" that Ali & Derrick developed for regional and whole-molecule integration of various functions.
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 are also some analytic formulas available.
Originally posted by @PaulWAyers in #13 (comment)
aim_weights
.MolecularGrid
class:
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.MolecularGrid
class. (Note: here is where tv-*.txt
points stored in needed; check #11)Running the tests generates a lot of warning messages, which should be addressed.
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
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?
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.
AtomicGridClass
:
RadilGrid
& UnitSphereGrid
to generate points/weights.tv-*.
in json/yml
file.AtomicGridClass
.From the current code:
npz
file based on l_max
.l_max
and return UnitSphereGrid
.UnitSphereGrid
Class:
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?
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.
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
OneDGrid
Class, just a namedtuple with points
and weights
attributes. Two categories of functions should be present:
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_
.
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.
At the moment the AIM weights are not mutiplied in, which makes the behavior of MolGrid different from its base class.
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
.
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.
Originally posted by @tovrstra in #99 (comment)
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.)
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.
At some point we should get rid of all the old HORTON commits. Cloning the repo just takes too long at the moment.
Implement SG1 as mentioned in issue #12.
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
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.
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.
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.
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
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).
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.
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.
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]
@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.
Get rid of get_fn
file, and use importlib.resources
to load test files.
I think the Molecular dipole moment calculation uses the wrong origin.
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.
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.