# 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 \qquad \qquad &\text{Cartesian moments}\\ m_{lm} &= \int | \textbf{r} - \textbf{R}_c|^l S_{l}^m(\theta, \phi) f(\textbf{r}) d\textbf{r} \qquad \qquad &\text{Spherical moments} \\ m_n &= \int | \textbf{r} - \textbf{R}_c|^{n} f(\textbf{r}) d\textbf{r} \qquad \qquad &\text{Radial moments}\\ m_{nlm} &= \int | \textbf{r} - \textbf{R}_c|^{n+1} S_l^m(\theta, \phi) f(\textbf{r}) d\textbf{r} \qquad \qquad &\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 $$(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 form a wavefunction $$\Psi$$: $$\vec{\mu} = \int \Psi \hat{\mu} \Psi \vec{r}$$ which results in

$\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.

[1]:

from iodata import load_one

print("Dipole Moments ", mol.moments)

Dipole Moments  {(1, 'c'): array([0.16176897, 0.58865022, 0.39905597])}


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

[2]:

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 moment integrals can be calculated.

[3]:

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, type = from_iodata(mol)

# Compute the electron density
rdm = mol.one_rdms["scf"]
electron_density = evaluate_density(rdm, basis, mol_grid.points, coord_type=type)
# electron_density[np.all(mol_grid.points < 1e-10, axis=1)] += 36.0

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 [0.1617419  0.58865474 0.39904899]
Dipole moment true [0.16176897 0.58865022 0.39905597]
Mean error 1.285537519197614e-05 with maximum error 2.7064523889591685e-05