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 git+https://github.com/theochem/denspart.git
#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
from urllib.request import urlretrieve
url = "https://media.githubusercontent.com/media/theochem/chemtools/master/chemtools/data/examples/ch2o_q%2B0.fchk"
urlretrieve(url, "ch2o_q+0.fchk")
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 Hirshfeldatomic 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:895: 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 Hirshfeldatomic 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="slater", 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
charges, f_ekl = hirshfeld_i_partitioning(mol, mol_grid, dataset="slater")
formatted_charges = ", ".join(f"{charge:.6f}" for charge in charges)
print(f"The iterative Hirshfeld atomic charges are: {formatted_charges}")
print(f"The KL divergence is: {f_ekl: .6f}")
16.00000757755703
16.00000757755703
Step 1: KL divergence = 0.210637, Charge convergence = 0.085084
Step 2: KL divergence = 0.211514, Charge convergence = 0.041133
Step 3: KL divergence = 0.211633, Charge convergence = 0.026425
Step 4: KL divergence = 0.211254, Charge convergence = 0.017815
Step 5: KL divergence = 0.210808, Charge convergence = 0.012326
Step 6: KL divergence = 0.210434, Charge convergence = 0.008639
Step 7: KL divergence = 0.210154, Charge convergence = 0.006092
Step 8: KL divergence = 0.209952, Charge convergence = 0.004308
Step 9: KL divergence = 0.209809, Charge convergence = 0.003051
Step 10: KL divergence = 0.209707, Charge convergence = 0.002162
Step 11: KL divergence = 0.209636, Charge convergence = 0.001532
Step 12: KL divergence = 0.209586, Charge convergence = 0.001086
Step 13: KL divergence = 0.209550, Charge convergence = 0.000770
The iterative Hirshfeld atomic charges are: -0.375750, 0.317309, 0.029229, 0.029203
The KL divergence is: 0.209550
Minimal Basis Iterative Stockholder (MBIS)#
Minimal Basis Iterative Stockholder (MBIS) analysis is a modern variant of stockholder partitioning. It assigns each atom a minimal basis of exponential functions (1s-type orbitals), then optimizes their effective nuclear charges to minimize the Kullback-Leibler divergence.
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"])
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
MBIS Partitioning
Building local grids
Integral of density: 16.00000757755703
Optimization
#Iter #Call ekld kld -constraint grad.rms cputime (s)
----- ----- ----------- ----------- ----------- ----------- -----------
1 1 0.4285957 0.4286032 -7.5776e-06 1.8799e-01 0.0312500
2 2 0.2924672 0.0568064 2.3566e-01 1.5988e-01 0.0312500
3 3 1.0164087 0.1545821 8.6183e-01 6.8674e-01 0.0312500
4 4 0.1945874 -0.3732477 5.6784e-01 7.4238e-02 0.0468750
5 4 0.1945874 -0.3732477 5.6784e-01 7.4238e-02 0.0468750
6 5 0.1691688 -0.1884295 3.5760e-01 4.6008e-02 0.0468750
7 6 0.1490751 0.1410514 8.0237e-03 1.2698e-02 0.0468750
8 7 101.0193243 90.8367515 1.0183e+01 3.4993e+01 0.0468750
9 8 38.6460663 29.6045556 9.0415e+00 1.9695e+00 0.0468750
10 9 0.1485582 0.1241837 2.4374e-02 1.5509e-02 0.0625000
11 10 0.1482585 0.1366128 1.1646e-02 1.2042e-02 0.0468750
12 11 0.1480900 0.1219824 2.6108e-02 1.0412e-02 0.0468750
13 11 0.1480900 0.1219824 2.6108e-02 1.0412e-02 0.0468750
14 12 0.1471067 0.1434323 3.6744e-03 5.1439e-03 0.0156250
15 13 0.1468437 0.1363995 1.0444e-02 3.2687e-03 0.0312500
16 14 0.1467341 0.1257893 2.0945e-02 3.2799e-03 0.0312500
17 15 1.7832261 -1.3532625 3.1365e+00 6.9541e-01 0.0312500
18 16 0.1467692 0.0582722 8.8497e-02 5.3239e-03 0.0312500
19 17 0.1466437 0.0982640 4.8380e-02 3.3309e-03 0.0312500
20 18 0.1466257 0.0889693 5.7656e-02 3.0806e-03 0.0312500
21 19 0.1466142 0.0836620 6.2952e-02 4.6650e-03 0.0312500
22 20 0.1465833 0.0853964 6.1187e-02 4.2674e-03 0.0312500
23 21 0.1463610 0.1254476 2.0913e-02 2.5517e-03 0.0312500
24 21 0.1463610 0.1254476 2.0913e-02 2.5517e-03 0.0312500
25 22 0.1463134 0.1396872 6.6263e-03 6.9197e-04 0.0312500
26 23 32.7328217 31.4698877 1.2629e+00 2.0671e+01 0.0312500
27 24 0.1884439 0.0278973 1.6055e-01 9.6488e-02 0.0468750
28 25 0.1552530 0.0797075 7.5545e-02 4.7298e-02 0.0312500
29 26 0.1462913 0.1383208 7.9706e-03 1.7764e-03 0.0468750
30 27 0.1462849 0.1383544 7.9305e-03 1.1643e-03 0.0312500
31 28 0.1464835 0.1382201 8.2634e-03 4.3696e-03 0.1093750
32 29 0.1462549 0.1404687 5.7862e-03 1.4548e-03 0.0625000
33 30 0.1462531 0.1406423 5.6108e-03 1.2196e-03 0.0312500
34 31 0.1462490 0.1412579 4.9911e-03 3.9893e-04 0.0312500
35 31 0.1462490 0.1412579 4.9911e-03 3.9893e-04 0.0312500
36 32 1.5891325 4.3380434 -2.7489e+00 7.5917e-01 0.0312500
37 33 0.3055173 0.7013629 -3.9585e-01 1.2228e-01 0.0312500
38 34 0.1773710 0.2123839 -3.5013e-02 5.0770e-02 0.0468750
39 35 0.1462476 0.1430553 3.1923e-03 1.8744e-04 0.0312500
40 36 0.1462468 0.1450331 1.2137e-03 1.6329e-04 0.0312500
41 36 0.1462468 0.1450331 1.2137e-03 1.6329e-04 0.0312500
42 37 0.1462466 0.1455151 7.3153e-04 9.8291e-05 0.0312500
43 38 0.1462466 0.1459874 2.5915e-04 3.4720e-05 0.0156250
44 39 0.1462465 0.1459791 2.6744e-04 2.0650e-05 0.0312500
45 39 0.1462465 0.1459791 2.6744e-04 2.0650e-05 0.0312500
46 40 0.1462465 0.1461684 7.8168e-05 5.5287e-06 0.0156250
47 41 0.1462465 0.1461693 7.7210e-05 7.1008e-06 0.0312500
48 42 0.1462465 0.1461695 7.7039e-05 7.0702e-06 0.0312500
49 43 0.1462465 0.1462118 3.4736e-05 4.0696e-06 0.0312500
50 43 0.1462465 0.1462118 3.4736e-05 4.0696e-06 0.0312500
51 44 0.1462465 0.1462431 3.3808e-06 9.5945e-07 0.0312500
52 44 0.1462465 0.1462431 3.3808e-06 9.5945e-07 0.0312500
53 45 0.1462465 0.1462418 4.7386e-06 3.0263e-07 0.0312500
54 46 0.1462465 0.1462452 1.3464e-06 1.3910e-07 0.0312500
55 46 0.1462465 0.1462452 1.3464e-06 1.3910e-07 0.0312500
56 47 0.1462465 0.1462460 5.0907e-07 3.0463e-08 0.0468750
57 47 0.1462465 0.1462460 5.0907e-07 3.0463e-08 0.0468750
58 48 0.1462465 0.1462466 -9.3156e-08 7.9485e-09 0.0468750
59 48 0.1462465 0.1462466 -9.3156e-08 7.9485e-09 0.0468750
60 49 0.1462465 0.1462466 -1.1231e-07 1.6471e-09 0.0468750
----- ----- ----------- ----------- ----------- ----------- -----------
Optimizer message: "`gtol` termination condition is satisfied."
Total charge: -7.5775570e-06
Sum atomic charges: -7.4652498e-06
Atomic Charges = [-0.41152524 0.33567696 0.03793329 0.03790753]