Constructing Molecular Grids

Molecular grids are particularly useful in quantum chemistry for integrating functions that concentrate and are “spike” around finite set of points, usually the atomic nuclei. This is particularly needed for DFT computation or atom-in-molecule analysis.

Molecular Grid Details

It starts by first decomposing a function into atomic contributions and writing out each atomic integral in spherical coordinates. As proposed by Becke, this is done as

\[\int_{\mathbb{R}^3} f(\vec{\textbf{r}}) d\vec{\textbf{r}} = \sum_{A} \int_{\mathbb{R}^3} w_A(\vec{\textbf{r}}) f(\vec{\textbf{r}}) d\vec{\textbf{r}},\]

where \(w_A\) is known as the nuclear weight function (or atom in molecule weights) such that it has value one close to the center and decay’s continuous over other centers with the condition that \(\sum_{A} w_A(\vec{\textbf{r}}) = 1\) for all \(\vec{\textbf{r}}\). The module Molecular Grids is responsible for computing integrals in this fashion with either Becke or Hirshfeld atom in molecule weights. The process to integrate each individual integral is done by first converting to spherical coordinates \((r, \theta, \phi)\):

\[\int_{\mathbb{R}^3} w_A(\vec{\textbf{r}}) f(\vec{\textbf{r}}) d\vec{\textbf{r}} = \int \int \int w_A(r, \theta, \phi) f(r, \theta, \phi) r^2 \sin(\phi) dr d\theta d\phi\]

Then a one-dimensional grid is chosen over the radial grid (see radial and radial transform module) with weights \(w_i^{rad}\) and an angular grid with weights \(w_j^{ang}\) is chosen to integrate over the angles \((\theta, \phi)\) including the sin factor. The combination of the two is handled by the atomic grid module with weights \(w_{ij} = w_i^{rad} w_j^{ang} r^2\) to achieve the numerical integral

\[\int \int \int w_A(r, \theta, \phi) f(r, \theta, \phi) r^2 \sin(\phi) dr d\theta d\phi \approx \sum_{i=1}^{N_{rad}} \sum_{j=1}^{N_{ang}} w_{ij}(r, \theta, \phi) w_A(r, \theta, \phi) f(r, \theta, \phi).\]

For any general function \(f: \mathbb{R}^3\rightarrow \mathbb{R}\), Grid package offers various grids in Cubic class for constructing hyper-rectangular grids. If \(f\) is periodic, then the periodic module is useful.

Example: Formaldehyde This example computes the integral of the electron density of Formaldehyde obtained from Gaussian calculation. The integral of the electron density should be the sum of the atomic charges (16 for Formaldehyde).

[1]:

# Generate the electron density, coordinates and charges from .fchk file
from iodata import load_one
from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

mol = load_one("ch2o_q+0.fchk")

print(f"Coordinates of the atoms \n {mol.atcoords}\n")
print(f"Atomic numbers of the atom \n {mol.atnums}")

atnums = mol.atnums
atcoords = mol.atcoords

# Construct molecular basis from wave-function information read by IOData
basis = from_iodata(mol)

# Get the RDM to compute electron density
rdm = mol.one_rdms["scf"]

Coordinates of the atoms
 [[ 2.27823914e+00  4.13899085e-07  3.12033662e-07]
 [ 1.01154892e-02  1.09802629e-07 -6.99333116e-07]
 [-1.09577141e+00  1.77311416e+00  1.42544321e-07]
 [-1.09577166e+00 -1.77311468e+00  2.44755133e-07]]

Atomic numbers of the atom
 [8 6 1 1]

Construction of Molecular Grid

Molecular grid offers various methods to construct MolGrid object. Each method requires the atom in molecule weights \(w_A\) to be provided.

MolGrid.__init__: Provide list of atomic grids for each center, and the atomic numbers for each center. This is useful to explicitly place atomic grids at each center.
MolGrid.from_preset(): Generates atomic grid for each center based on the atomic numbers, atomic coordinates, a radial grid to integrate over \([0, \infty)\), and a preset parameter from the list [“coarse”, “medium”, “fine”, “veryfine”, “ultrafine”, “insane”] that controls the degree of the angular grid.
MolGrid.from_size(): Similar to from_preset method but instead of string preset parameter, it has an integer size which controls the size of the angular grid inside the atomic grid class.
MolGrid.from_pruned(): Gives more control to the size of the angular grid based on the radius. This is done by providing r_sectors list and sector_degs that controls what kind of angular grid to place at each radius point. This is useful for placing lower degree angular grids close to the nucleus and larger degree angular grids further away.

See the function documentation for more information. The first step here is to generate an one-dimensional grid to integrate over the radial coordinate \(r \in [0, \infty)\). This is done by constructing a Trapezoidal grid over \([-1, 1]\) and using the Becke Transformation to transform it to \([0, \infty)\).

[2]:

from grid.onedgrid import Trapezoidal
from grid.rtransform import LinearFiniteRTransform

n_rpoints = 400                         # Number of radial points
oned = Trapezoidal(n_rpoints)           # Trapezoidal grid of that size
rmin, R = 0.0, 0.5                      # Parameters for the Becke radial transform
rgrid = LinearFiniteRTransform(rmin, 8).transform_1d_grid(oned)

The Becke weights is constructed as the atom in molecular weights (also called nuclear weight functions)

[3]:

from grid.becke import BeckeWeights

aim_weights = BeckeWeights()

Constructor

The first method of constructing (point 1 above) is using the initializor. This is done by explicitly constructing an atomic grid over each center. The degree of the angular component of the atomic grid should depend on the charge of the atom. It is recommended that the larger the atom, the larger the degree of the angular component. Similarly for the radial grid, the larger the atom, the larger the radial grid component should be.

[4]:

from grid.molgrid import MolGrid
from grid.atomgrid import AtomGrid

# Construct molecular grid
atom_grid_oxygen = AtomGrid(rgrid, degrees=[51], center=atcoords[0])
atom_grid_carbon = AtomGrid(rgrid, degrees=[25], center=atcoords[1])
atom_grid_hydro1 = AtomGrid(rgrid, degrees=[10], center=atcoords[2])
atom_grid_hydro2 = AtomGrid(rgrid, degrees=[10], center=atcoords[3])
atom_grids = [atom_grid_oxygen, atom_grid_carbon, atom_grid_hydro1, atom_grid_hydro2]
molgrid = MolGrid(atnums, atom_grids, aim_weights, store=True)

density = evaluate_density(rdm, basis, molgrid.points)

integral = molgrid.integrate(density)
print(f"The integral of electron density of Formaldehyde is {integral}.")

/mnt/Data/Work/Ayers/QC-Devs/grid/src/grid/atomgrid.py:879: UserWarning: Lebedev weights are negative which can introduce round-off errors.
  sphere_grid = AngularGrid(degree=deg_i, method=method)

The integral of electron density of Formaldehyde is 15.999822947424516.

from_preset method

The disadvantage of the previous method is that it is difficult to choose the size of the angular degrees for each atom, and it can be cumbersome for larger atoms to construct atomic grids for each center. The from_preset method has a preset that handles it for the user that is simpler.

[5]:

from grid.molgrid import MolGrid

molgrid = MolGrid.from_preset(atnums, atcoords, rgrid=rgrid, preset="fine", aim_weights=aim_weights, store=True)

density = evaluate_density(rdm, basis, molgrid.points)
integral = molgrid.integrate(density)
print(f"The integral of electron density of Formaldehyde is {integral}.")

The integral of electron density of Formaldehyde is 15.999839766372306.

from_size method

This method is similar to before except you can explicitly set the size of the angular degrees. This is the same across all atoms.

[6]:

from grid.molgrid import MolGrid

molgrid = MolGrid.from_size(atnums, atcoords, rgrid=rgrid, size=40, aim_weights=aim_weights, store=True)

density = evaluate_density(rdm, basis, molgrid.points)
integral = molgrid.integrate(density)
print(f"The integral of electron density of Formaldehyde is {integral}.")

The integral of electron density of Formaldehyde is 16.000299762134755.

from_pruned method

This method allows explicit control to the angular degree \(L_i\) by breaking up the radial grid \([0, \infty)\) into sectors \([0, R a_1) \cup [R a_1, R a_2) \cup \cdots \cup [R a_Q, \infty)\) where \(R\) is the radius parameter, \(\{a_i\}\) is the r_sectors and each segment is associated to a angular degree \(L_i\), i.e.

\[\begin{split}\begin{align*} &L_1 \text{ when } r < R a_1 \\ &L_2 \text{ when } R a_1 \leq r < R a_2 \\ \vdots \\ &L_{Q+1} \text{ when } R a_{Q} < r. \end{align*}\end{split}\]

[7]:

from grid.molgrid import MolGrid

r_sectors = [
    [0.25, 1.0, 1.5],   # for oxygen
    [0.5, 1.0, 1.5],    # for carbon
    [1.0], [1.0]        # for hydrogen
]
d_sectors = [
    [10, 15, 25, 10],   # for oxygen
    [5, 10, 20, 10],    # for carbon
    [5, 10], [5, 10]    # for hydrogen
]
radius = 1.0
molgrid = MolGrid.from_pruned(atnums, atcoords, radius, r_sectors=r_sectors, d_sectors=d_sectors, aim_weights=aim_weights, rgrid=rgrid)

density = evaluate_density(rdm, basis, molgrid.points)
integral = molgrid.integrate(density)
print(f"The integral of electron density of Formaldehyde is {integral}.")

/mnt/Data/Work/Ayers/QC-Devs/grid/src/grid/atomgrid.py:879: UserWarning: Lebedev weights are negative which can introduce round-off errors.
  sphere_grid = AngularGrid(degree=deg_i, method=method)

The integral of electron density of Formaldehyde is 16.000504626151567.