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.

from iodata import load_one

mol = load_one("./ch2o_q+0.fchk")

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.

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

# Construct a radial grid
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.

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