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
from iodata import load_one
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}")
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")
radial_grid = BeckeRTransform(0.0, R=1.5).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.
[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_deriv = evaluate_density_gradient(rdm, basis, random_grid)
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