Molecular Grids

This example on Formaldehyde illustates all the features of using molecular grids.

The format checkpoint file wavefunction is included in “doc/notebooks/ch2o_q+0.fchk” and is read using the IOData package.

import numpy as np
from iodata import load_one

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

print(f"Coordinates of the atoms \n {mol.atcoords}")
print(f"Atomic numbers of the atom \n {mol.atnums}")
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]

The first step to constructing a molecular grid is to define a radial grid to integrate over the radial coordinate. This is done by first defining a GaussLegendre grid over \([-1, 1]\) and transforming it to \([0, \infty)\) using a Becke Transformation.

from grid.onedgrid import GaussLegendre
from grid.rtransform import BeckeRTransform

oned_grid = GaussLegendre(npoints=50)
print(f"One-dimesional points \n {oned_grid.points}")
radial_grid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned_grid)  #BeckeRTransform(0.0, R=2.0).transform_1d_grid(oned_grid)
print(f"Transformed points \n {radial_grid.points}")
One-dimesional points
 [-0.9988664  -0.99403197 -0.98535408 -0.97286439 -0.95661096 -0.93665662
 -0.91307856 -0.88596798 -0.85542977 -0.82158207 -0.78455583 -0.7444943
 -0.70155247 -0.65589647 -0.60770293 -0.5571583  -0.50445814 -0.44980633
 -0.39341431 -0.33550025 -0.27628819 -0.21600724 -0.15489059 -0.0931747
 -0.03109834  0.03109834  0.0931747   0.15489059  0.21600724  0.27628819
  0.33550025  0.39341431  0.44980633  0.50445814  0.5571583   0.60770293
  0.65589647  0.70155247  0.7444943   0.78455583  0.82158207  0.85542977
  0.88596798  0.91307856  0.93665662  0.95661096  0.97286439  0.98535408
  0.99403197  0.9988664 ]
Transformed points
 [8.50678848e-04 4.48941942e-03 1.10654689e-02 2.06316373e-02
 3.32634175e-02 4.90613930e-02 6.81530638e-02 9.06950874e-02
 1.16876073e-01 1.46920031e-01 1.81090580e-01 2.19696072e-01
 2.63095794e-01 3.11707472e-01 3.66016382e-01 4.26586392e-01
 4.94073421e-01 5.69241890e-01 6.52984919e-01 7.46349269e-01
 8.50566286e-01 9.67090581e-01 1.09764867e+00 1.24430061e+00
 1.40951880e+00 1.59628946e+00 1.80824471e+00 2.04983623e+00
 2.32656593e+00 2.64529648e+00 3.01467435e+00 3.44571510e+00
 3.95262549e+00 4.55397903e+00 5.27442985e+00 6.14726583e+00
 7.21830627e+00 8.55201815e+00 1.02414211e+01 1.24247214e+01
 1.53144537e+01 1.92511601e+01 2.48083999e+01 3.30139230e+01
 4.58609073e+01 6.76418771e+01 1.09055814e+02 2.03335260e+02
 5.01178390e+02 2.64494645e+03]

Once a radial grid is defined, there are several methods to constructing a molecular grid. The simplest of such method is used. Here, Becke weights are used as an atom in molecule weights.

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

mol_grid = MolGrid.from_preset(
    atnums=mol.atnums,          # The atomic numbers of Formaldehyde
    atcoords=mol.atcoords,      # The atomic coordinates of Formaldehyde
    rgrid=radial_grid,          # Radial grid used to construct atomic grids over each carbon, and hydrogen.
    preset="fine",              # Controls the angular degree. Can be "coarse", "medium", "fine", "veryfine", "ultrafine", "insane".
    aim_weights=BeckeWeights(), # Atom-in molecular weights: Becke weights
    store=True,                 # Stores the individual atomic grids, used for interpolation.

# Alternative construction includes:
# mol_grid = MolGrid.from_size(
#     atnums=mol.atnums,          # The atomic numbers of Formaldehyde
#     atcoords=mol.atcoords,      # The atomic coordinates of Formaldehyde
#     rgrid=radial_grid,          # Radial grid used to construct atomic grids over each carbon, and hydrogen.
#     size=50,                    # The angular degree of the atomic grid over each carbon, and hydrogen.
#     aim_weights=BeckeWeights(), # Atom-in molecular weights: Becke weights
# )


Example: Electron Density

This example shows how to use molecular grids to integrate the electron density \(\rho(\textbf{r})\). This should equal to the number of electrons 16 of Formaldehyde. This is calculated using the gbasis package

from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

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

# Compute the electron density
rdm = mol.one_rdms["scf"]
electron_density = evaluate_density(rdm, basis, mol_grid.points, coord_type=type)

# Integrate using the molecular grid
print(f"Integration of the electron density: {mol_grid.integrate(electron_density)}.")
Integration of the electron density: 15.999987778747279.


Example: Electron Density

This example shows how to interpolate any function using molecular grids. The electron density is used again here as the example.

# Construct interpolate function
interpolate = mol_grid.interpolate(

# Construct 100 random-set of points between [-5, 5]^3
random_grid = np.random.uniform(-5, 5, size=(100, 3))

# Evalaute using the interpolate
inter_vals = interpolate(random_grid)
actual_vals = evaluate_density(rdm, basis, random_grid, coord_type=type)

# Look at the error-difference
err = np.abs(inter_vals - actual_vals)
print(f"Max Error {np.max(err)},  Mean Error {np.mean(err)},  Std Dev {np.std(err)}")

Max Error 0.001302904658059434,  Mean Error 3.500902876292755e-05,  Std Dev 0.0001483415813999016

Derivative Interpolation

Example: Weizsacker Kinetic Energy Density

The Weizsacker kinetic energy density is defined as:

\[\tau_W (\textbf{r}) = \frac{1}{8} \frac{|\nabla \rho(\textbf{r}) |^2}{\rho(\textbf{r})}.\]

This example will show how to interpolate the Weizsacker kinetic energy density using the interpolation of the electron density \(\rho(\textbf{r})\).

from gbasis.evals.density import evaluate_density_gradient

# Computing the Weizsacker kinetic energy density
actual_deriv = evaluate_density_gradient(rdm, basis, random_grid, coord_type=type)
actual_dens = evaluate_density(rdm, basis, random_grid, coord_type=type)
norm_squared = np.linalg.norm(actual_deriv, axis=1) ** 2.0  # Take the norm-square of the gradient
kinetic_dens = (1.0 / 8.0) * norm_squared / actual_dens

# Evaluate the interpolation of the derivative and the actual derivative.
dens_interp = interpolate(random_grid, deriv=0)
deriv_interp = interpolate(random_grid, deriv=1)
deriv_norm = np.linalg.norm(deriv_interp, axis=1) ** 2.0
kinetic_interp = (1.0 / 8.0) * deriv_norm / dens_interp

err = np.abs(kinetic_dens - kinetic_interp)
print(f"Max Error {np.max(err)},  Mean Error {np.mean(err)},  Std Dev {np.std(err)}")
Max Error 0.0015169498269268777,  Mean Error 7.528647116405785e-05,  Std Dev 0.00023846069080800918

Solving Poisson Equation

Example: Electrostatic Potential

The Poisson equation is defined as:

\[\nabla^2 V(\textbf{r}) = - 4\pi \rho(\textbf{r}),\]

where \(\nabla^2\) is the Laplacian operator, \(V\) is the Couloumb potential and \(\rho\) is the electron density. This example will show how to use molecular grids to solve the Poisson equation for a given electron density. The solution is equivalent to the integral:

\[V(\textbf{r}_1) = \int \frac{\rho(\textbf{r}_2)}{| \textbf{r}_1 - \textbf{r}_2 | } d\textbf{r}_2.\]

This can be used to calculate the electrostatic potential of a molecule with coordinates \(\{R_A\}\) and charge \(\{Z_A\}\):

\[\Phi(\textbf{r}_1) = -\int \frac{\rho(\textbf{r}_2)}{| \textbf{r}_1 - \textbf{r}_2 | } d\textbf{r}_2 + \sum_{A} \frac{Z_A}{|\textbf{r}_1 - R_A|}\]
from grid.rtransform import InverseRTransform, BeckeRTransform
from grid.poisson import solve_poisson_bvp
from gbasis.evals.electrostatic_potential import electrostatic_potential

# Compute the electron density
dens = evaluate_density(rdm, basis, mol_grid.points, coord_type=type)

# Define a transformation from [0, \infty) to finite interval [-1, 1]
inverse_transf = InverseRTransform(BeckeRTransform(1e-30, R=1.5))

potential = solve_poisson_bvp(
    molgrid=mol_grid,                      # Provide the molecular density
    func_vals=dens,                        # The charge distribution
    transform=inverse_transf,              # Transformation to finite-interval
    include_origin=False,                  # Don't include the origin r=0
    remove_large_pts=10.0,                 # Remove points with radius larger than 10

# Calculate electrostatic potential with nuclear charges all set to zero.
desired_potential = -electrostatic_potential(
    basis, rdm, random_grid, nuclear_coords=mol.atcoords, nuclear_charges=np.array([0, 0, 0, 0]), coord_type=type
actual_potential = potential(random_grid)

err = np.abs(desired_potential - actual_potential)
print(f"Max Error {np.max(err)},  Mean Error {np.mean(err)},  Std Dev {np.std(err)}")
Max Error 0.015115750096519776,  Mean Error 0.0038874435708785506,  Std Dev 0.0032290409703695498