Clarification on Weight Distribution in Becke Weights Near Nucleus

[izxle]

Hi QC-devs Team,

We’ve been working with the grid repository as part of a project where we're constructing the promolecular density for a diatomic molecule and integrating various terms using the molecular grids.

While using the Becke weights function, we noticed that the weights near the nucleus are not close to 1, which I had expected based on my understanding. Is this behavior intentional, or might it indicate an issue? If it’s intentional, could you provide some insight into the reasoning behind it or point me toward relevant references?

I thought I’d ask here first to understand the design of the weights better.

Thank you in advance for your help and for providing such valuable resources to the community.

Here below I include an example.

import numpy as np
import matplotlib.pyplot as plt

from atomdb import make_promolecule
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import Trapezoidal
from grid.rtransform import LinearFiniteRTransform


# Define atomic information
atnums = np.array([10, 10])
atcoords = np.array([[0, 0, 0], [0, 0, 10]])

# Create promolecular objects
promol = make_promolecule(atnums, atcoords, dataset='gaussian')

# Generate the radial grid using Trapezoidal method
oned = Trapezoidal(400) # Create a Trapezoidal grid with n_rpoints
rgrid = LinearFiniteRTransform(0, 10).transform_1d_grid(oned) # Apply linear transformation to the grid

# Generate the molecular grid
mol_grid = MolGrid.from_preset(
    atnums=atnumsAB,
    atcoords=atcoordsAB,
    rgrid=rgrid,
    preset="fine",
    aim_weights=BeckeWeights(),
    store=True,
)

# Plotting weights along the z-axis near x=0, y=0 for mol_grid
# Find points where abs(x) < 0.1 and abs(y) < 0.1
th = 1e-1
indices = np.where((np.abs(mol_grid.points[:, 0]) < th) & (np.abs(mol_grid.points[:, 1]) < th))[0]

# Extract z coordinates and weights at those points
z_coords = mol_grid.points[indices, 2]
weights = mol_grid.weights[indices]

# Sort data for plotting
sorted_indices = np.argsort(z_coords)
z_coords = z_coords[sorted_indices]
weights = weights[sorted_indices]

# Plot
plt.plot(z_coords, weights)
plt.xlabel("z-coordinate")
plt.ylabel("Weights")
plt.title("Weights along z-axis (near x<=0.1, y<=0.1)")
plt.show()

[Ali-Tehrani]

Hello @izxle,

The weights you are plotted aren't exactly the Becke-weight, but the molecular grid weights. The molecular grid weights are a multiplication of the Becke-weights of the atom and the radial grid and the angular (spherical) grid weight. I'll recommend the published Grid paper if you want the mathematical details.

You can see the code where it calculates the weights here, specifically Line 103:

grid/src/grid/molgrid.py

Lines 88 to 103 in 8bf58af

	if callable(aim_weights):
	self._aim_weights = aim_weights(self._points, self._atcoords, atnums, self._indices)
	elif isinstance(aim_weights, np.ndarray):
	if aim_weights.size != size:
	raise ValueError(
	"aim_weights is not the same size as grid.\n"
	f"aim_weights.size: {aim_weights.size}, grid.size: {size}."
	)
	self._aim_weights = aim_weights
	else:
	raise TypeError(f"Not supported aim_weights type, got {type(aim_weights)}.")
	# initialize parent class
	super().__init__(self.points, self._atweights * self._aim_weights)

.

What you want is to plot Line 83 to get the Becke-Weights around an atom. Here is a quick example script that shows you get all ones if you only use Becke-Weights near the nucleus

import numpy as np
from grid.becke import BeckeWeights

# Define atomic information
atnums = np.array([10, 10])
atcoords = np.array([[0, 0, 0], [0, 0, 10]])

becke = BeckeWeights()

pts = np.array([
   # Points on Atom A
   [0.0, 0.0, 0.0],
   [0.0, 0.0, 0.05],
   [0.0, 0.05, 0.0],
   [0.05, 0.0, 0.0],

   # Points on Atom B
   [0.0, 0.0, 10.0],
   [0.05, 0.0, 10.0],
   [0.0, 0.05, 10.0]]
)

# Index 0 to 3 corresponds to Atom A in pts
# Index 3 to 6 corresponds to Atom B in pts
indices = np.array([0, 4, 7])
aim_weights = becke(pts, atcoords, atnums, indices=indices)
print(aim_weights)
# Output: [1. 1. 1. 1. 1. 1. 1.]

[PaulWAyers]

@izxle if you want to use promolecules, I'd recommend using the AtomDB package.

The notebook gives a (too simple) example of using the promolecule with the grid.

We'll try to add a better example in the next few days. I have one in hand, just need to get it polished and added to the documentation for AtomDB.