Awesome
- Changelog
- Licensed under MPL v2.0 (except John Burkardt’s Lebedev code which is redistributed under LGPL v3.0)
Numgrid
Numgrid is a library that produces numerical integration grid for molecules based on atom coordinates, atom types, and basis set information. This library provides Rust and Python bindings.
Who are the people behind this code?
Authors
- Radovan Bast
Contributors
- Roberto Di Remigio (OS X testing, streamlined Travis testing, better C++, error handling)
For a list of all the contributions see https://github.com/dftlibs/numgrid/contributors.
Acknowledgements
- Simon Neville (reporting issues)
- Jaime Axel Rosal Sandberg (reporting issues)
This tool uses SPHERE_LEBEDEV_RULE, a C library written by John Burkardt which computes a Lebedev quadrature rule over the surface of the unit sphere in 3D, see also: http://people.sc.fsu.edu/~jburkardt/c_src/sphere_lebedev_rule/sphere_lebedev_rule.html
This library uses and acknowledges the MolSSI BSE (https://molssi-bse.github.io/basis_set_exchange/), which is a rewrite of the Basis Set Exchange (https://bse.pnl.gov/bse/portal) and is a collaboration between the Molecular Sciences Software Institute (http://www.molssi.org) and the Environmental Molecular Sciences Laboratory (https://www.emsl.pnl.gov).
Citation
If you use this tool in a program or publication, please acknowledge its author(s)
@misc{numgrid,
author = {Bast, Radovan},
title = {Numgrid: Numerical integration grid for molecules},
month = {7},
year = {2024},
publisher = {Zenodo},
version = {v2.1.1},
doi = {10.5281/zenodo.1470276},
url = {https://doi.org/10.5281/zenodo.1470276}
}
@misc{sphere_lebedev_rule,
author = {Burkardt, John},
title = {SPHERE_LEBEDEV_RULE: Quadrature Rules for the Unit Sphere},
year = {2010},
url = {https://people.sc.fsu.edu/~jburkardt/c_src/sphere_lebedev_rule/sphere_lebedev_rule.html}
}
I kindly ask you to also cite the latter since Numgrid is basically a "shell" around SPHERE_LEBEDEV_RULE, with added radial integration and molecular partitioning.
Would you like to contribute?
Yes please! Please follow this excellent guide. We do not require any formal copyright assignment or contributor license agreement. Any contributions intentionally sent upstream are presumed to be offered under terms of the Mozilla Public License Version 2.0.
Requirements
- Python (3.8 - 3.12).
- For the Rust version: A Rust installation.
Installation
Installing via pip
python -m pip install numgrid
Building from sources and testing
Building the code:
cargo build --release
Testing the Rust interface:
cargo test --release
Running also the longer tests:
cargo test --release -- --ignored
Testing the Python layer:
pip install -r requirements.txt # ideally into a virtual environment
maturin develop
pytest tests/test.py
API
The API changed
The API changed (sorry!) for easier maintenance and simpler use:
-
No initialization or deallocation necessary.
-
One-step instead of two steps (since the radial grid generation time is negligible compared to space partitioning, it did not make sense anymore to separate these steps and introduce a state).
-
alpha_min
is given as dictionary which saves an argument and simplifies explaining the API. -
The library now provides Rust and Python bindings. It used to provide C and Fortran bindings. The C/Fortran code lives on the cpp-version branch. I might bring the C interfaces back into the Rust code if there is sufficient interest/need.
Units
Coordinates are in bohr.
Python example
As an example let us generate a grid for the water molecule:
import numgrid
radial_precision = 1.0e-12
min_num_angular_points = 86
max_num_angular_points = 302
proton_charges = [8, 1, 1]
center_coordinates_bohr = [(0.0, 0.0, 0.0), (1.43, 0.0, 1.1), (-1.43, 0.0, 1.1)]
# cc-pVDZ basis
alpha_max = [
11720.0, # O
13.01, # H
13.01, # H
]
alpha_min = [
{0: 0.3023, 1: 0.2753, 2: 1.185}, # O
{0: 0.122, 1: 0.727}, # H
{0: 0.122, 1: 0.727}, # H
]
for center_index in range(len(center_coordinates_bohr)):
# atom grid using explicit basis set parameters
coordinates, weights = numgrid.atom_grid(
alpha_min[center_index],
alpha_max[center_index],
radial_precision,
min_num_angular_points,
max_num_angular_points,
proton_charges,
center_index,
center_coordinates_bohr,
hardness=3,
)
# atom grid using basis set name
# this takes a second or two for the REST API request
coordinates, weights = numgrid.atom_grid_bse(
"cc-pVDZ",
radial_precision,
min_num_angular_points,
max_num_angular_points,
proton_charges,
center_index,
center_coordinates_bohr,
hardness=3,
)
# radial grid (LMG) using explicit basis set parameters
radii, weights = numgrid.radial_grid_lmg(
alpha_min={0: 0.3023, 1: 0.2753, 2: 1.185},
alpha_max=11720.0,
radial_precision=1.0e-12,
proton_charge=8,
)
# radial grid (LMG) using basis set name
radii, weights = numgrid.radial_grid_lmg_bse(
basis_set="cc-pVDZ",
radial_precision=1.0e-12,
proton_charge=8,
)
# radial grid with 100 points using Krack-Koster approach
radii, weights = numgrid.radial_grid_kk(num_points=100)
# angular grid with 14 points
coordinates, weights = numgrid.angular_grid(num_points=14)
Notes and recommendations
-
The smaller the
radial_precision
, the better grid. -
For
min_num_angular_points
andmax_num_angular_points
, see “Angular grid” below. -
alpha_max
is the steepest basis set exponent. -
alpha_min
is a dictionary and holds the smallest exponents for each angular momentum (order does not matter). -
Using
center_index
we tell the code which of the atom centers is the one we have computed the grid for. -
num_angular_grid_points
has to be one of the many supported Lebedev grids (see table on the bottom of this page).
Rust interface
Needs to be documented better but the library exposes functions with the same name as the Python interface and probably the best example on how it can be used are the integration tests.
Saving grid in NumPy format
The current API makes is relatively easy to export the computed grid in NumPy format.
In this example we save the angular grid coordinates and weights to two separate files in NumPy format:
import numgrid
import numpy as np
coordinates, weights = numgrid.angular_grid(14)
np.save("angular_grid_coordinates.npy", coordinates)
np.save("angular_grid_weights.npy", weights)
Parallelization
The Becke partitioning step is parallelized using Rayon. In other words, this step should be able to use all available cores on the computer or computing node. Since grids are currently generated atom by atom, it is also possible to parallelize "outside" by the caller.
If you need to limit the number of cores used by Rayon, you can set (in this case limiting to 4 threads):
export RAYON_NUM_THREADS=4
Space partitioning
The molecular integration grid is generated from atom-centered grids by scaling the grid weights according to the Becke partitioning scheme, JCP 88, 2547 (1988). The default Becke hardness is 3.
Radial grid
Two choices are available:
- Lindh-Malmqvist-Gagliardi (https://dx.doi.org/10.1007/s002140100263)
- Krack-Köster (https://doi.org/10.1063/1.475719)
Advantage of LMG scheme: The range of the radial grid is basis set dependent. The precision can be tuned with one single radial precision parameter. The smaller the radial precision, the better quality grid you obtain. The basis set (more precisely the Gaussian primitives/exponents) are used to generate the atomic radial grid range. This means that a more diffuse basis set generates a more diffuse radial grid.
Advantage of the KK scheme: parameter-free.
Angular grid
The angular grid is generated according to Lebedev and Laikov [A quadrature formula for the sphere of the 131st algebraic order of accuracy, Russian Academy of Sciences Doklady Mathematics, Volume 59, Number 3, 1999, pages 477-481].
The angular grid is pruned. The pruning is a primitive linear interpolation between the minimum number and the maximum number of angular points per radial shell. The maximum number is reached at 0.2 times the Bragg radius of the center.
The higher the values for minimum and maximum number of angular points, the better.
For the minimum and maximum number of angular points the code will use the following table and select the closest number with at least the desired precision:
{6, 14, 26, 38, 50, 74, 86, 110, 146,
170, 194, 230, 266, 302, 350, 434, 590, 770,
974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470,
3890, 4334, 4802, 5294, 5810}
Taking the same number for the minimum and maximum number of angular points switches off pruning.