# 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.

[1]:

import numpy as np

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

[2]:

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

oned_grid = GaussLegendre(npoints=50)
print(f"One-dimesional points \n {oned_grid.points}\n")


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.

[3]:

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
# )


## Integration

### 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

[4]:

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

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

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

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

Integration of the electron density:  15.999988.


## Local Properties

This example shows how to extract atomic grids from a molecular grid and integrate the electron density to get atomic charges.

[5]:

# Extract the atomic grid of the Oxygen atom
atom_index = 1
o_grid = mol_grid.get_atomic_grid(atom_index)
print(f"{len(mol_grid.points)} points in molecular grid")
print(f"{len(o_grid.points)} points in atomic grid of atom {atom_index}\n")

# partition electron density into atomic electron densities (each point is assigned to an atom)
# 1. Compute weighted electron density by Becke integration weights
w_dens = dens * mol_grid.aim_weights
# 2. Get weighted electron density of points assigned to Oxygen atom
o_idcs = np.arange(mol_grid.indices[atom_index], mol_grid.indices[atom_index + 1])
o_dens = w_dens[o_idcs]
print(f"Shape of molecular electron density : {dens.shape}")
print(f"Shape of atomic density of Oxygen   : {o_dens.shape}\n")

# Integrate using the atomic grid
o_charge = o_grid.integrate(o_dens)

print(f"Integration of the Oxygen atom electron density on the Oxygen grid: {o_charge: .6f}")
print(f"Atomic charge of the Oxygen atom (q = Z - N): {mol.atnums[1] - o_charge: .6f}")

16448 points in molecular grid
4980 points in atomic grid of atom 1

Shape of molecular electron density : (16448,)
Shape of atomic density of Oxygen   : (4980,)

Integration of the Oxygen atom electron density on the Oxygen grid:  6.213951
Atomic charge of the Oxygen atom (q = Z - N): -0.213951


## Interpolation

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

[6]:

# Construct interpolate function
interpolate = mol_grid.interpolate(func_vals=dens)

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

# Evaluate using the interpolate
inter_vals = interpolate(random_grid)
actual_vals = evaluate_density(rdm, basis, random_grid)

# Compute the interpolation error
err = np.abs(inter_vals - actual_vals)
print(
f"Max Error={np.max(err): .6f}; Mean Error={np.mean(err): .6f}; Std Dev Error={np.std(err): .6f}"
)

Max Error= 0.000438; Mean Error= 0.000014; Std Dev Error= 0.000054


## Derivative Interpolation

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})$$.

[7]:

from gbasis.evals.density import evaluate_density_gradient

# Computing the Weizsacker kinetic energy density
actual_dens = evaluate_density(rdm, basis, random_grid)
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

# Compute the interpolation error
err = np.abs(kinetic_dens - kinetic_interp)
print(
f"Max Error={np.max(err): .6f}; Mean Error={np.mean(err): .6f}; Std Dev Error={np.std(err): .6f}"
)

Max Error= 0.003902; Mean Error= 0.000068; Std Dev Error= 0.000395


## Solving Poisson Equation

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|}$
[8]:

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)

# 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]),
)
actual_potential = potential(random_grid)

# Compute the interpolation error
err = np.abs(desired_potential - actual_potential)
print(
f"Max Error={np.max(err): .6f}; Mean Error={np.mean(err): .6f}; Std Dev Error={np.std(err): .6f}"
)

Max Error= 0.021631; Mean Error= 0.003719; Std Dev Error= 0.003231