Becke grid analysis

All numerical integrations are performed by means of Gauss-Chebychev and Lebedev quadrature using the Becke grid [Becke, 1988]. The radial Gauss-Chebychev points describe the distance from an atomic center, while the Lebedev points [Lebedev, 1976] describe directions on the unit sphere. Their tensor product gives an atom-centered grid.

For molecules, these atom-centered grids overlap. Becke weights turn the overlap into a smooth partition of unity: each grid point receives one weight per atom, and the weights at that point sum to one. Integrals over the whole molecule can therefore be evaluated as a weighted sum of atomic-grid contributions.

Atomic fuzzy cells

It is possible to produce a contour plot of the fuzzy cells (molecular weights) on a plane using the pydft.MolecularGrid.calculate_weights_at_points() method. An example is provided below.

from pydft import MolecularGrid
from pyqint import MoleculeBuilder
import numpy as np
import matplotlib.pyplot as plt

# construct molecule
mol = MoleculeBuilder().from_name('benzene')
cgfs, atoms = mol.build_basis('sto3g')

# construct molecular grid
molgrid = MolecularGrid([a for a in mol], cgfs, 
                        nshells = {'C' : 32, 'H' : 32}, 
                        nangpts = {'C' : 194, 'H' : 194}
                        )

# produce grid of sampling points to calculate the atomic
# weight coefficients for
N = 150
sz = 8
x = np.linspace(-sz,sz,N)
xv,yv = np.meshgrid(x,x)
points = np.array([[x,y,0] for x,y in zip(xv.flatten(),yv.flatten())])

# calculate the atomic weights
mweights = molgrid.calculate_weights_at_points(points, k=3)

# plot the atomic weights
plt.figure(dpi=144)
plt.imshow(np.max(mweights,axis=0).reshape((N,N)),
           extent=(-sz,sz,-sz,sz), interpolation='bicubic')
plt.xlabel('x [a.u.]')
plt.ylabel('y [a.u.]')
plt.colorbar()
plt.grid(linestyle='--', color='black', alpha=0.5)

# add the atoms to the plot
r = np.zeros((len(atoms), 3))
for i,at in enumerate(atoms):
    r[i] = at[0]
plt.scatter(r[0:6,0], r[0:6,1], s=50.0, color='grey', edgecolor='black')
plt.scatter(r[6:12,0], r[6:12,1], s=50.0, color='white', edgecolor='black')

plt.tight_layout()

In the script above, we color every point by the maximum value for each of the atomic weights. When this maximum value is one, this implies that the grid point belongs to a single atom. When multiple atoms ‘share’ a grid point, the maximum value among the atomic weights will be lower than one.

Becke fuzzy-cell weights for benzene projected onto the molecular plane.

Note

Producing such a contour plot is only meaningful for planar molecules such as benzene. For more complex molecules such as methane, it is rather difficult to make sense of the fuzzy cells upon projection on a plane.

Grid points

To obtain the set of grid points on which the numerical integration (quadrature), is performed, we can invoke the pydft.MolecularGrid.get_grid_coordinates() method.

from pydft import MolecularGrid
from pyqint import MoleculeBuilder
import matplotlib.pyplot as plt

# construct molecule
mol = MoleculeBuilder().from_name('co')
cgfs, atoms = mol.build_basis('sto3g')

# construct molecular grid
molgrid = MolecularGrid([a for a in mol], cgfs, 
                        nshells = {'C' : 32, 'O' : 32}, 
                        nangpts = {'C' : 50, 'O' : 50},
                        lmax = {'C' : 5, 'O' : 5}
                        )
molgrid.initialize()

# obtain the grid points
gridpoints = molgrid.get_grid_coordinates()

# plot the atomic weights
colors = '#DD0000', '#222222'
fig = plt.figure(dpi=300, figsize=(8,6))
ax = fig.add_subplot(projection='3d')
for i in range(0, len(gridpoints)):
    ax.scatter(gridpoints[i][:,0], gridpoints[i][:,1], gridpoints[i][:,2],
               s = 1.5, alpha=0.5, color=colors[i])

# set axes
ax.set_xlim(-5,5)
ax.set_ylim(-5,5)
ax.set_zlim(-5,5)
ax.set_xlabel('x [a.u.]')
ax.set_ylabel('y [a.u.]')
ax.set_zlabel('z [a.u.]')
ax.set_box_aspect(aspect=None, zoom=0.8)

plt.show()

Atom-centered quadrature points used to perform the molecular integration.