Multipole Moments

Every grid class has the ability to compute the multipole moment integral of a function over various centers. It can compute the following types:

\begin{split}\begin{align*} m_{n_x, n_y, n_z} &= \int (x - X_c)^{n_x} (y - Y_c)^{n_y} (z - Z_c)^{n_z} f(r) dr \quad &\text{Cartesian moments}\\ m_{lm} &= \int | \textbf{r} - \textbf{R}_c|^l S_{l}^m(\theta, \phi) f(\textbf{r}) d\textbf{r} \quad &\text{Spherical moments} \\ m_n &= \int | \textbf{r} - \textbf{R}_c|^{n} f(\textbf{r}) d\textbf{r} \quad &\text{Radial moments}\\ m_{nlm} &= \int | \textbf{r} - \textbf{R}_c|^{n+1} S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r} \quad &\text{Radial with spherical moments} \end{align*}\end{split}

for some function $$f : \mathbb{R}^3\rightarrow \mathbb{R}$$, where $$S_l^m$$ is the regular, real solid harmonics, $$(n_x, n_y, n_z)$$ are the Cartesian orders over some center $$\textbf{R}_c = (X_c, Y_c, Z_c)$$ and $$(l, m)$$ are the angular order and degree.

This example illustrates how to compute the dipole moment of water. This is defined as the observable acting on a wavefunction $$\Psi$$: $$\vec{\mu} = \int \Psi \hat{\mu} \Psi \vec{r}$$ which results in the calculation of the dipole moment as

$\vec{\mu} = \sum_{i=1}^{N_{atoms}} Z_i (\vec{R_i} - \vec{R_c}) - \int (\vec{r} - \vec{R_c}) \rho(\vec{r}) dr,$

where $$N_{atoms}$$ is the number of atoms, $$Z_i$$ is the atomic charge of the ith atom, $$\vec{R_i}$$ is the ith coordinate of the atom, $$\vec{R_c}$$ is the center of the molecule and $$\rho$$ is the electron density of the molecule.

IOData is used to first read the wavefunction information of Formaldehyde.

[4]:

from iodata import load_one

print("Dipole Moments ", mol.moments)

Dipole Moments  {(1, 'c'): array([-9.39793529e-01,  2.44832724e-08, -2.02253053e-07])}


In order to compute the moment integral, a molecular grid class is constructed.

[5]:

from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from grid.onedgrid import GaussLegendre
from grid.rtransform import BeckeRTransform

oned_grid = GaussLegendre(npoints=150)
radial_grid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned_grid)  #BeckeRTransform(0.0, R=2.0).transform_1d_grid(oned_grid)

# Construct Molecular grid with angular degree of 50 for each atom.
mol_grid = MolGrid.from_size(
atnums=mol.atnums,          # The atomic numbers of Formaldehyde
atcoords=mol.atcoords,      # The atomic coordinates of Formaldehyde
rgrid=radial_grid,          # Radial grid used to construct atomic grids over each carbon, and hydrogen.
size=130,                    # The angular degree of the atomic grid over each carbon, and hydrogen.
aim_weights=BeckeWeights(), # Atom-in molecular weights: Becke weights,
)


The dipole moment can then be calculated.

[6]:

import numpy as np
from gbasis.wrappers import from_iodata
from gbasis.evals.density import evaluate_density

from grid.utils import dipole_moment_of_molecule

# Construct molecular basis from wave-function information read by IOData
basis = from_iodata(mol)

# Compute the electron density
rdm = mol.one_rdms["scf"]
electron_density = evaluate_density(rdm, basis, mol_grid.points)

true = dipole_moment_of_molecule(mol_grid, electron_density, mol.atcoords, mol.atnums)
desired = mol.moments[(1, "c")]
print(f"Dipole moment calculated {true}")
print(f"Dipole moment true {desired}")

err = np.abs(true - desired)
print(f"Mean error {np.mean(err)} with maximum error {np.max(err)}")

Dipole moment calculated [-9.39751887e-01  6.42314269e-05  4.02806905e-05]
Dipole moment true [-9.39793529e-01  2.44832724e-08 -2.02253053e-07]
Mean error 4.877715321240814e-05 with maximum error 6.42069436543974e-05