(Iterative) Hirshfeld method
In this tutorial, we will introduce how to use the horton_part API method to execute (iterative) Hirshfeld partitioning methods.
First, we import the required libraries.
[1]:
import numpy as np
from gbasis.evals.eval import evaluate_basis
from gbasis.wrappers import from_iodata
from grid import AtomGrid, ExpRTransform, UniformInteger
from iodata import load_one
from setup import prepare_grid_and_dens, print_results
from horton_part import HirshfeldIWPart, HirshfeldWPart, ProAtomDB, ProAtomRecord
np.set_printoptions(precision=3, linewidth=np.inf, suppress=True)
In (iterative) Hirshfeld methods, a database of atomic radial density profiles is required. Therefore, different single-atom DFT or other ab initio calculations should be executed. Here, we use the Gaussian package. The checkpoint files, in fchk format, are generated using Gaussian tools as well.
Next, the radial atomic density is obtained by computing the spherical average of the atomic density. This process can be applied using the gbasis and grid packages.
[2]:
def prepare_record(filename):
"""Prepare molecular grid and density."""
mol = load_one(filename)
# Specify the integration grid
rtf = ExpRTransform(5e-4, 2e1, 120 - 1)
uniform_grid = UniformInteger(120)
rgrid = rtf.transform_1d_grid(uniform_grid)
# Get the spin-summed density matrix
one_rdm = mol.one_rdms.get("post_scf", mol.one_rdms.get("scf"))
basis = from_iodata(mol)
assert len(mol.atnums) == 1
grid = AtomGrid.from_preset(atnum=mol.atnums[0], preset="fine", rgrid=rgrid)
basis_grid = evaluate_basis(basis, grid.points)
rho = np.einsum("ab,bp,ap->p", one_rdm, basis_grid, basis_grid, optimize=True)
spline = grid.spherical_average(rho)
radial_rho = spline(rgrid.points)
record = ProAtomRecord(mol.atnums[0], mol.charge, mol.energy, rgrid, radial_rho)
return record
In this tutorial, we focus on the water molecule and therefore only need data for hydrogen and oxygen atoms. It should be noted that data for both cations and anions are also used.
[3]:
h_1 = prepare_record("./data/atoms/001__h_001_q+00/mult02/atom.fchk")
h_2 = prepare_record("./data/atoms/001__h_002_q-01/mult01/atom.fchk")
o_7 = prepare_record("./data/atoms/008__o_007_q+01/mult04/atom.fchk")
o_8 = prepare_record("./data/atoms/008__o_008_q+00/mult03/atom.fchk")
o_9 = prepare_record("./data/atoms/008__o_009_q-01/mult02/atom.fchk")
records = [h_1, h_2, o_7, o_8, o_9]
db = ProAtomDB(records)
Next, we prepare the molecular grid and densities as described in the Basic Usage section. It should be noted that the radial grids for atoms in the water molecule should match the atomic radial grids used for the corresponding single atoms in the atomic database calculations.
[4]:
def prepare_argument_dict(mol, grid, rho):
"""Prepare basic input arguments for all AIM methods."""
kwargs = {
"coordinates": mol.atcoords,
"numbers": mol.atnums,
"pseudo_numbers": mol.atnums,
"grid": grid,
"moldens": rho,
"lmax": 3,
"proatomdb": db,
}
return kwargs
Hirshfeld method
The Hirshfeld method can be executed using the HirshfeldWPart class.
[5]:
mol, grid, rho = prepare_grid_and_dens("data/h2o.fchk")
kwargs = prepare_argument_dict(mol, grid, rho)
h = HirshfeldWPart(**prepare_argument_dict(mol, grid, rho))
h.do_all()
print_results(h)
The number of electrons: 10.000003764139395
Coordinates of the atoms
[[ 0. 0. 0.224]
[-0. 1.457 -0.896]
[-0. -1.457 -0.896]]
Atomic numbers of the atom
[8 1 1]
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Initialized: HirshfeldWPart
Scheme : Hirshfeld
Proatomic DB : <horton_part.core.proatomdb.ProAtomDB object at 0x107d43190>
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Computing atomic dispersion coefficients.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.305 0.153 0.153]
cartesian multipoles:
[[-0.305 -0. -0. -0.253 -4.416 -0. -0. -3.647 -0. -4.093 -0. -0. -0.637 -0. -0. -0. -0. -0.247 -0. -1.401]
[ 0.153 -0. 0.167 -0.131 -0.758 -0. -0. -0.643 0.035 -0.689 -0. 0.206 -0.2 -0. -0. -0. 0.468 -0.127 0.164 -0.454]
[ 0.153 0. -0.167 -0.131 -0.758 0. 0. -0.643 -0.035 -0.689 -0. -0.206 -0.2 -0. -0. -0. -0.468 -0.127 -0.164 -0.454]]
Iterative Hirshfeld method
The iterative Hirshfeld method can be executed using HirshfeldIWpart class.
[6]:
hi = HirshfeldIWPart(**prepare_argument_dict(mol, grid, rho))
hi.do_all()
print_results(hi)
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Initialized: HirshfeldIWPart
Scheme : Hirshfeld-I
Convergence threshold : 1.0e-06
Maximum iterations : 500
Proatomic DB : <horton_part.core.proatomdb.ProAtomDB object at 0x107d43190>
Iteration Change Entropy
1 4.91894e-02 2.00489e-01
2 3.09349e-02 1.80862e-01
3 1.95877e-02 1.90830e-01
4 1.26979e-02 2.03872e-01
5 8.40308e-03 2.14964e-01
6 5.64913e-03 2.23432e-01
7 3.84155e-03 2.29631e-01
8 2.63387e-03 2.34078e-01
9 1.81637e-03 2.37237e-01
10 1.25775e-03 2.39468e-01
11 8.73435e-04 2.41038e-01
12 6.07777e-04 2.42140e-01
13 4.23518e-04 2.42912e-01
14 2.95413e-04 2.43453e-01
15 2.06199e-04 2.43832e-01
16 1.43998e-04 2.44098e-01
17 1.00594e-04 2.44283e-01
18 7.02895e-05 2.44413e-01
19 4.91226e-05 2.44503e-01
20 3.43338e-05 2.44567e-01
21 2.39993e-05 2.44611e-01
22 1.67764e-05 2.44642e-01
23 1.17278e-05 2.44664e-01
24 8.19875e-06 2.44679e-01
25 5.73174e-06 2.44690e-01
26 4.00710e-06 2.44697e-01
27 2.80143e-06 2.44702e-01
28 1.95853e-06 2.44706e-01
29 1.36925e-06 2.44709e-01
30 9.57279e-07 2.44710e-01
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Computing atomic dispersion coefficients.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.862 0.431 0.431]
cartesian multipoles:
[[-0.862 0. -0. 0.189 -5.147 -0. -0. -4.812 -0. -5.012 -0. -0. -0.105 0. 0. 0. 0. 0.58 -0. 0.174]
[ 0.431 -0. 0.077 -0.04 -0.392 -0. -0. -0.389 0.02 -0.375 -0. 0.088 -0.057 0. 0. -0. 0.264 -0.039 0.083 -0.136]
[ 0.431 0. -0.077 -0.04 -0.392 0. 0. -0.389 -0.02 -0.375 0. -0.088 -0.057 0. -0. 0. -0.264 -0.039 -0.083 -0.136]]
[ ]: