"""Post-SCF molecular properties: atomic charges, bond orders, dipole.
The standard "sanity-check" output you expect from every QC program,
built on top of a converged SCF result.
Public API
----------
.. autofunction:: mulliken_charges
.. autofunction:: loewdin_charges
.. autofunction:: mayer_bond_orders
.. autofunction:: dipole_moment
.. autofunction:: center_of_mass
All functions accept the ``result`` object returned by
:func:`vibeqc.run_rhf`, :func:`vibeqc.run_uhf`, :func:`vibeqc.run_rks`,
or :func:`vibeqc.run_uks`; the matching :class:`Molecule` and
:class:`BasisSet` used to run the SCF; and for dipole moments an
optional origin.
Implementation notes
--------------------
Mulliken and Löwdin population analyses are AO-basis-dependent in well-
known ways — Mulliken in particular is very sensitive to diffuse
functions. The output is useful for trend-watching across a series
(e.g. charge transfer along a reaction coordinate) but atomic charges
beyond the leading digit should never be taken too seriously. Mayer
bond orders are rotation-invariant and less basis-sensitive.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import TYPE_CHECKING, Iterable, Optional, Sequence
import numpy as np
from ._vibeqc_core import (
BasisSet,
Molecule,
compute_dipole,
compute_overlap,
)
if TYPE_CHECKING: # pragma: no cover — import cycle avoidance
from ._vibeqc_core import RHFResult # noqa: F401
__all__ = [
"DipoleMoment",
"mulliken_charges",
"loewdin_charges",
"mayer_bond_orders",
"dipole_moment",
"center_of_mass",
]
# Standard atomic masses (u), index 0 unused. Dalton approximate values —
# good to ~1e−3 u — plenty for centre-of-mass computation. We ship a
# subset sufficient for H–Kr; for heavier elements callers pass an
# explicit origin.
_ATOMIC_MASSES: tuple[float, ...] = (
0.0,
1.008, 4.003, # H He
6.94, 9.012, 10.81, 12.011, 14.007, 15.999, 18.998, 20.180, # Li–Ne
22.990, 24.305, 26.982, 28.085, 30.974, 32.06, 35.45, 39.948, # Na–Ar
39.098, 40.078, 44.956, 47.867, 50.942, 51.996, 54.938, 55.845, # K–Fe
58.933, 58.693, 63.546, 65.38, 69.723, 72.630, 74.922, 78.971, # Co–Se
79.904, 83.798, # Br Kr
)
# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------
def _total_density(result) -> np.ndarray:
"""Return the total (closed-shell + open-shell) density matrix.
RHF/RKS store the already-combined density on ``.density``; UHF/UKS
expose ``density_alpha`` and ``density_beta`` separately."""
if hasattr(result, "density_alpha"):
return np.asarray(result.density_alpha) + np.asarray(result.density_beta)
return np.asarray(result.density)
def _shell_to_atom(basis: BasisSet) -> np.ndarray:
"""1-D int array, length ``nbasis``, mapping each AO to the 0-based
atom index it lives on. Computed from ``basis.shells()`` (which is
public C++ API exposed for basis-set I/O)."""
shells = basis.shells()
per_ao: list[int] = []
for shell in shells:
# 2L+1 spherical AOs per shell (vibe-qc forces set_pure(true)
# at basis construction; Cartesian shells would have a
# different multiplicity, but that is not our configuration).
n = 2 * shell.l + 1
per_ao.extend([int(shell.atom_index)] * n)
return np.asarray(per_ao, dtype=np.int64)
def _per_atom_sum(ao_values: np.ndarray, ao_to_atom: np.ndarray,
n_atoms: int) -> np.ndarray:
"""Sum an nbasis-length array into per-atom totals."""
out = np.zeros(n_atoms, dtype=np.float64)
for a, v in zip(ao_to_atom, ao_values):
out[a] += v
return out
[docs]
def center_of_mass(molecule: Molecule) -> np.ndarray:
"""Centre of mass (bohr). Atomic masses from a built-in table up to
Z = 36; callers with heavier elements should specify the origin to
dipole_moment directly."""
atoms = list(molecule.atoms)
if not atoms:
return np.zeros(3)
total_mass = 0.0
com = np.zeros(3)
for atom in atoms:
z = int(atom.Z)
if z < len(_ATOMIC_MASSES):
m = _ATOMIC_MASSES[z]
else:
# Crude fall-back: 2 u per nucleon ≈ A, and A ≈ 2 Z on average
# for light elements — good enough not to throw. Heavy-element
# users should pass an explicit origin.
m = 2.0 * z
pos = np.array([atom.xyz[0], atom.xyz[1], atom.xyz[2]])
com += m * pos
total_mass += m
if total_mass == 0.0:
return np.zeros(3)
return com / total_mass
# ---------------------------------------------------------------------------
# Mulliken population analysis
# ---------------------------------------------------------------------------
[docs]
def mulliken_charges(result, basis: BasisSet, molecule: Molecule) -> np.ndarray:
"""Mulliken atomic partial charges q_A = Z_A − Σ_{μ ∈ A} (P·S)_μμ.
Returns a 1-D array of length ``n_atoms``. Charges sum to the total
molecular charge (``molecule.charge``) to machine precision. The
overall partition is AO-basis-dependent; Mulliken is sensitive to
diffuse functions and should be used for trend analysis rather than
for quantitative charge assignment.
"""
P = _total_density(result)
S = np.asarray(compute_overlap(basis))
PS_diag = np.einsum("ij,ji->i", P, S) # diag(P · S)
ao_to_atom = _shell_to_atom(basis)
n_atoms = len(molecule.atoms)
electron_pop = _per_atom_sum(PS_diag, ao_to_atom, n_atoms)
Z = np.array([atom.Z for atom in molecule.atoms], dtype=np.float64)
return Z - electron_pop
# ---------------------------------------------------------------------------
# Löwdin population analysis
# ---------------------------------------------------------------------------
[docs]
def loewdin_charges(result, basis: BasisSet, molecule: Molecule) -> np.ndarray:
"""Löwdin (symmetric orthogonalisation) atomic partial charges:
q_A = Z_A − Σ_{μ ∈ A} (S^{1/2} · P · S^{1/2})_μμ.
Compared to Mulliken this partitions the off-diagonal overlap
equally between the two AOs involved — less basis-sensitive in
practice.
"""
P = _total_density(result)
S = np.asarray(compute_overlap(basis))
# S = U · diag(s) · U.T → S^{1/2} = U · diag(√s) · U.T
w, V = np.linalg.eigh(S)
if np.min(w) < 1e-10:
# Defensive: raise the way run_rhf's symmetric-orthogonalisation
# does. In practice vibe-qc's basis sets never trigger this on
# molecules; periodic near-linear-dependence surfaces earlier.
raise ValueError(
f"loewdin_charges: overlap matrix is near-singular "
f"(min eigenvalue {np.min(w):.2e})"
)
S_half = V @ np.diag(np.sqrt(w)) @ V.T
SPS = S_half @ P @ S_half
diag = np.diag(SPS)
ao_to_atom = _shell_to_atom(basis)
n_atoms = len(molecule.atoms)
electron_pop = _per_atom_sum(diag, ao_to_atom, n_atoms)
Z = np.array([atom.Z for atom in molecule.atoms], dtype=np.float64)
return Z - electron_pop
# ---------------------------------------------------------------------------
# Mayer bond orders
# ---------------------------------------------------------------------------
[docs]
def mayer_bond_orders(result, basis: BasisSet,
molecule: Molecule) -> np.ndarray:
"""Mayer bond-order matrix, shape (n_atoms, n_atoms).
For closed-shell systems:
B_AB = Σ_{μ ∈ A} Σ_{ν ∈ B} (P·S)_μν · (P·S)_νμ
For UHF / UKS we replace (P·S)(P·S) with the spin-resolved
definition 2·[(Pα·S)(Pα·S) + (Pβ·S)(Pβ·S)] which reduces to the
closed-shell formula when Pα = Pβ = P/2.
Diagonal entries are the "free valence" of each atom. The off-
diagonals are the chemically-meaningful bond orders.
"""
S = np.asarray(compute_overlap(basis))
ao_to_atom = _shell_to_atom(basis)
n_atoms = len(molecule.atoms)
if hasattr(result, "density_alpha"):
Pa = np.asarray(result.density_alpha)
Pb = np.asarray(result.density_beta)
PS_a = Pa @ S
PS_b = Pb @ S
# Element-wise Mayer: M_μν = 2·[(PS_a)_μν (PS_a)_νμ
# + (PS_b)_μν (PS_b)_νμ]
M = 2.0 * (PS_a * PS_a.T + PS_b * PS_b.T)
else:
P = np.asarray(result.density)
PS = P @ S
M = PS * PS.T # broadcasting element-wise; equivalent to
# M_μν = (PS)_μν · (PS)_νμ since the matrix is
# real.
bond_orders = np.zeros((n_atoms, n_atoms), dtype=np.float64)
for mu in range(M.shape[0]):
a = ao_to_atom[mu]
for nu in range(M.shape[1]):
b = ao_to_atom[nu]
if a != b:
bond_orders[a, b] += M[mu, nu]
# Symmetrise numerically — analytical B_AB = B_BA for real AOs.
bond_orders = 0.5 * (bond_orders + bond_orders.T)
return bond_orders
def prominent_bonds(
bond_orders: np.ndarray,
molecule: Molecule,
*,
threshold: float = 0.10,
) -> list[tuple[int, int, float]]:
"""Return ``[(i, j, B_ij)]`` pairs with ``i < j`` and
``B_ij >= threshold`` — convenience for formatting a compact
bond-order table in the log output."""
n = bond_orders.shape[0]
out: list[tuple[int, int, float]] = []
for i in range(n):
for j in range(i + 1, n):
if bond_orders[i, j] >= threshold:
out.append((i, j, float(bond_orders[i, j])))
# Sort by descending bond order so the covalent bonds float to the top.
out.sort(key=lambda t: -t[2])
return out
# ---------------------------------------------------------------------------
# Dipole moment
# ---------------------------------------------------------------------------
_BOHR_TO_DEBYE = 2.541746473 # 1 e·bohr = 2.541746473 Debye
[docs]
@dataclass
class DipoleMoment:
"""Dipole moment components in atomic units (e·bohr), plus Debye."""
x: float
y: float
z: float
origin: tuple[float, float, float]
@property
def total(self) -> float:
r"""\|mu\| in atomic units (e\*bohr)."""
return float(np.sqrt(self.x ** 2 + self.y ** 2 + self.z ** 2))
@property
def total_debye(self) -> float:
return self.total * _BOHR_TO_DEBYE
[docs]
def components_debye(self) -> tuple[float, float, float]:
return (self.x * _BOHR_TO_DEBYE,
self.y * _BOHR_TO_DEBYE,
self.z * _BOHR_TO_DEBYE)
[docs]
def dipole_moment(
result,
basis: BasisSet,
molecule: Molecule,
*,
origin: Optional[Sequence[float]] = None,
) -> DipoleMoment:
"""Electric dipole moment of a converged SCF calculation.
``origin`` (bohr) defaults to the molecular centre of mass, which
makes the dipole origin-independent for neutral systems (the
convention every standard QC code uses). For charged systems the
dipole depends on origin; pass an explicit vector if you need a
particular reference.
"""
if origin is None:
origin_vec = center_of_mass(molecule)
else:
origin_vec = np.asarray(origin, dtype=np.float64)
if origin_vec.shape != (3,):
raise ValueError("dipole_moment: origin must be a 3-vector (bohr)")
dip = compute_dipole(basis, [float(x) for x in origin_vec])
Mx = np.asarray(dip.x)
My = np.asarray(dip.y)
Mz = np.asarray(dip.z)
P = _total_density(result)
# Electronic contribution (electrons are negative): -tr(P · M_c).
mu_e_x = -np.einsum("ij,ji->", P, Mx)
mu_e_y = -np.einsum("ij,ji->", P, My)
mu_e_z = -np.einsum("ij,ji->", P, Mz)
# Nuclear contribution (with origin shift): Σ_A Z_A (R_A − O).
mu_n_x = 0.0
mu_n_y = 0.0
mu_n_z = 0.0
for atom in molecule.atoms:
z = float(atom.Z)
mu_n_x += z * (atom.xyz[0] - origin_vec[0])
mu_n_y += z * (atom.xyz[1] - origin_vec[1])
mu_n_z += z * (atom.xyz[2] - origin_vec[2])
return DipoleMoment(
x=float(mu_e_x + mu_n_x),
y=float(mu_e_y + mu_n_y),
z=float(mu_e_z + mu_n_z),
origin=(float(origin_vec[0]),
float(origin_vec[1]),
float(origin_vec[2])),
)