This set of notes shows how to do Hirshfeld(-style) partitioning using AtomDB, Grid, GBasis, and IOdata.
14. Hirshfeld Partitioning#
Given a molecule, build the promolecule using neutral atoms, $\( \rho^0(\mathbf{r}) = \sum_A \rho_A^0(\mathbf{r} - \mathbf{R}_A) \)\( The atomic weights and charges are \)\( w_A(\mathbf{r}) = \frac{\rho_A^0(\mathbf{r})}{\rho^0(\mathbf{r})} \)\( \)\( q_A = Z_A - \int w_A(\mathbf{r}) \rho(\mathbf{r}) d\mathbf{r} \)\( The KL information loss is \)\( \mathcal{F} = \int \rho^0(\mathbf{r}) \log \left( \frac{\rho^0(\mathbf{r})}{\rho(\mathbf{r})} \right) d\mathbf{r} \)$
# If running on Google Colab or any environment where the following packages are not installed, uncomment the following lines
#!pip install qc-gbasis
#!pip install qc-iodata
#!pip install qc-atomDB
#!pip install qc-grid
#!pip install qc-denspart
#Key import statements
from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density
import numpy as np
import atomdb
from iodata import load_one
from grid.onedgrid import GaussLegendre
from grid.rtransform import BeckeRTransform
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[2], line 7
3 from gbasis.evals.density import evaluate_density
4 import numpy as np
5 import atomdb
6 from iodata import load_one
----> 7 from grid.onedgrid import GaussLegendre
8 from grid.rtransform import BeckeRTransform
9 from grid.becke import BeckeWeights
10 from grid.molgrid import MolGrid
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/__init__.py:25
21 """Grid Module."""
22 # flake8: noqa
---> 25 from grid.angular import *
26 from grid.atomgrid import *
27 from grid.basegrid import *
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/angular.py:75
72 import numpy as np
73 from importlib_resources import files
---> 75 from grid.basegrid import Grid
77 # Lebedev dictionary for grid's number of points (keys) and degrees (values)
78 LEBEDEV_NPOINTS = {
79 6: 3,
80 18: 5,
(...) 110 5810: 131,
111 }
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/basegrid.py:25
22 import numpy as np
23 from scipy.spatial import cKDTree
---> 25 from grid.utils import (
26 convert_cart_to_sph,
27 generate_orders_horton_order,
28 solid_harmonics,
29 )
32 class Grid:
33 """Basic Grid class for grid information storage."""
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/utils.py:23
20 """Utility Module."""
22 import numpy as np
---> 23 from scipy.special import sph_harm
25 _bragg = np.array(
26 [
27 np.nan, # index 0, place holder
(...) 114 ]
115 )
117 _cambridge = np.array(
118 [
119 np.nan, # index 0, place holder
(...) 206 ]
207 )
ImportError: cannot import name 'sph_harm' from 'scipy.special' (/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/special/__init__.py)
def make_grid(mol, rad_points=100, preset="ultrafine"):
'''Make a grid object from a molecule
Args:
mol (IOData): molecule object
rad_points (int): number of radial points
preset (str): preset grid type
"coarse", "medium", "fine", "veryfine", "ultrafine", "insane"
Returns:
mol_grid (Grid): grid object
'''
# create a radial grid with rad_points points
oned_grid = GaussLegendre(rad_points)
radial_grid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned_grid)
# create a pruned molecular grid with the preset angular convention
mol_grid = MolGrid.from_preset(mol.atnums,mol.atcoords,preset,rgrid=radial_grid,store=True)
return mol_grid
# Evaluate the molecular density on a grid
def eval_density(mol, mol_grid, spin="a+b"):
'''Evaluate the molecular density on a grid
Args:
mol (IOData): molecule object
mol_grid (Grid): grid object
spin (str): spin type
"a+b" (total), "a", "b", "a-b" (magnetization)
Returns:
dens (np.ndarray): molecular density
'''
# Construct molecular basis from wave-function information read by IOData
basis = from_iodata(mol)
# Compute the 1DM and then the electron density
if spin == "a+b":
rdm = mol.one_rdms["scf"]
dens = evaluate_density(rdm, basis, mol_grid.points)
else:
rdm = mol.one_rdms["scf_spin"]
if spin == "a":
dens = evaluate_density(rdm[0], basis, mol_grid.points)
elif spin == "b":
dens = evaluate_density(rdm[1], basis, mol_grid.points)
elif spin == "a-b":
dens = evaluate_density(rdm[0] - rdm[1], basis, mol_grid.points)
else:
raise ValueError(f"Invalid spin type: {spin}")
return dens
# KL Divergence
def ext_kl_divergence(rho0, rho, mol_grid):
'''Compute the extended KL divergence between two densities
Args:
rho0 (np.ndarray): reference density
rho (np.ndarray): density
mol_grid (Grid): grid object
Returns:
F (float): KL divergence
'''
# We only evaluate the integrand in locations where the molecular
# and promolecular density are sufficiently large.
mask = (rho > np.finfo(float).eps) & (rho0 > np.finfo(float).eps)
# Initialize the integrand with zeros
integrand = np.zeros_like(rho)
# Compute the integrand only where the mask is true
integrand[mask] = rho[mask] * np.log(rho[mask] / rho0[mask]) + rho0[mask] - rho[mask]
return mol_grid.integrate(integrand)
def hirshfeld_partitioning(mol, mol_grid, dataset="gaussian"):
'''Hirshfeld partitioning of a IOData molecule
Args:
mol (IOData): molecule object
mol_grid (Grid): grid object
dataset (AtomDB): atomic database
"gaussian", "slater", "hci", "numeric"
Returns:
q (np.ndarray): atomic atomic charges
f_ekl (float): extended KL information loss
'''
# evaluate the molecular density on the grid
dens = eval_density(mol, mol_grid)
print(mol_grid.integrate(dens))
# evaluate the promolecular density on the grid
promol = atomdb.make_promolecule(mol.atnums, mol.atcoords, dataset=dataset)
promol_dens = promol.density(mol_grid.points, log=True)
# Compute the KL divergence
f_ekl = ext_kl_divergence(promol_dens, dens, mol_grid)
# compute the atomic weights, populations, and
at_charges = np.zeros(len(mol.atnums))
weight = np.zeros(len(dens))
for i in range(len(mol.atnums)):
proatom = atomdb.make_promolecule(mol.atnums[i,np.newaxis], mol.atcoords[i, np.newaxis], dataset=dataset)
proatom_dens = proatom.density(mol_grid.points, log=True)
# Do not attempt to compute the weight if the pro-atom density is zero
# Set it equal to zero directly.
mask = proatom_dens > np.finfo(float).eps
weight = np.zeros_like(proatom_dens)
weight[mask] = proatom_dens[mask] / promol_dens[mask]
population = mol_grid.integrate(dens * weight)
at_charges[i] = mol.atnums[i] - population
return at_charges, f_ekl
# Example of Hirshfeld partitioning of a molecule
mol = load_one("ch2o_q+0.fchk")
mol_grid = make_grid(mol, rad_points=100, preset="ultrafine")
promol = atomdb.make_promolecule(mol.atnums, mol.atcoords, dataset="gaussian")
print(promol.coeffs)
promol_dens = promol.density(mol_grid.points, log=True)
mol_grid.integrate(promol_dens)
charges, f_ekl = hirshfeld_partitioning(mol, mol_grid, dataset="gaussian")
formatted_charges = ", ".join(f"{charge:.6f}" for charge in charges)
print(f"The atomic charges are: {formatted_charges}")
print(f"The KL divergence is: {f_ekl: .6f}")
c:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:896: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
[1.0, 1.0, 1.0, 1.0]
16.00000757755703
The atomic charges are: -0.216955, 0.124497, 0.046231, 0.046218
The KL divergence is: 0.225480
14. Iterative Hirshfeld Partitioning#
Iterative Hirshfeld partitioning is the same, but one uses charged proatoms. The pro-atom charge is taken from the previous iteration.
I.e., the atomic weights are: $\( w_{A,k}(\mathbf{r}) = \frac{\rho_{A,q_{k-1}}^0(\mathbf{r})}{\rho_{k-1}^0(\mathbf{r})} \)\( where the atomic densities are chosen to have the charge from the previous iteration (initialized to the Hirshfeld densities), \)\( q_{A;k-1} = Z_A - \int w_{A;k-1}(\mathbf{r}) \rho(\mathbf{r}) d\mathbf{r} \)$ The charged atomic densities are defined by linear interpolation, as is “correct” in DFT.
def hirshfeld_i_partitioning(mol, mol_grid, dataset="gaussian", max_iter=100, tol=1e-3):
'''Iterative Hirshfeld partitioning of a IOData molecule
Args:
mol (IOData): molecule object
mol_grid (Grid): grid object
dataset (AtomDB): atomic database
"gaussian", "slater", "hci", "numeric"
max_iter (int): maximum number of iterations
tol (float): convergence tolerance
Returns:
q (np.ndarray): atomic atomic charges
f_ekl (float): extended KL information loss
'''
# Initialize with Hirshfeld partitioning
at_charges, f_ekl = hirshfeld_partitioning(mol, mol_grid, dataset)
# evaluate the molecular density on the grid
dens = eval_density(mol, mol_grid)
print(mol_grid.integrate(dens))
# Iterate until convergence or maximum number of iterations is exceeded
for i in range(max_iter):
#Store old charges
at_charges_old = at_charges.copy()
# evaluate the promolecular density on the grid
promol = atomdb.make_promolecule(mol.atnums, mol.atcoords, charges=at_charges, dataset=dataset)
promol_dens = promol.density(mol_grid.points, log=True)
# Compute the KL divergence
f_ekl = ext_kl_divergence(promol_dens, dens, mol_grid)
# compute the atomic weights, populations, and
weight = np.zeros(len(dens))
for j_atom in range(len(mol.atnums)):
mol.atnums[j_atom,np.newaxis], mol.atcoords[j_atom, np.newaxis],
proatom = atomdb.make_promolecule(mol.atnums[j_atom,np.newaxis], mol.atcoords[j_atom,np.newaxis],
charges=at_charges[j_atom,np.newaxis],dataset=dataset)
proatom_dens = proatom.density(mol_grid.points, log=True)
# Do not attempt to compute the weight if the pro-atom density is zero
# Set it equal to zero directly.
mask = proatom_dens > np.finfo(float).eps
weight = np.zeros_like(proatom_dens)
weight[mask] = proatom_dens[mask] / promol_dens[mask]
population = mol_grid.integrate(dens * weight)
at_charges[j_atom] = mol.atnums[j_atom] - population
# Check to see if the maximum change in the charges is
# less than the tolerance
at_charge_convergence = np.max(np.abs(at_charges - at_charges_old))
print(f"Step {i+1}: KL divergence = {f_ekl: .6f}, Charge convergence = {at_charge_convergence: .6f}")
if at_charge_convergence < tol:
break
if i == max_iter - 1:
print("Maximum number of iterations exceeded")
return at_charges, f_ekl
Minimal Basis Iterative Stockholder (MBIS)#
from denspart.mbis import MBISProModel
from denspart.vh import optimize_reduce_pro_model
from denspart.properties import compute_radial_moments, compute_multipole_moments, safe_ratio
mol = load_one("ch2o_q+0.fchk")
mol_grid = make_grid(mol, rad_points=100, preset="ultrafine")
density = eval_density(mol, mol_grid, spin="a+b")
print("MBIS Partitioning")
pro_model_init = MBISProModel.from_geometry(mol.atnums, mol.atcoords)
pro_model, localgrids = optimize_reduce_pro_model(pro_model_init,mol_grid,density)
results = pro_model.to_dict()
results.update({"charges": pro_model.charges})
print("Atomic Charges = ", results["charges"])