"""Phase 17a-1 — molecular finite-difference Hessian + harmonic frequencies.
Two-sided finite differences on the analytic atomic gradient give the
3N × 3N nuclear Hessian in Hartree/bohr². Diagonalizing the mass-
weighted Hessian (with translation + rotation modes projected out)
returns the 3N − 6 (or 3N − 5 for linear molecules) harmonic
vibrational frequencies in cm⁻¹.
API
---
::
from vibeqc import compute_hessian_fd, HessianFDOptions
result = compute_hessian_fd(
mol, basis_name="sto-3g", method="RHF",
)
print(result.frequencies_cm1)
For DFT methods, pass ``scf_options`` carrying the functional plus
``grid_options``::
opts = vibeqc.RKSOptions()
opts.functional = "PBE"
result = compute_hessian_fd(
mol, "6-31g*", method="RKS", scf_options=opts,
)
Cost
----
6N gradient evaluations (one ± displacement per Cartesian DOF). Each
evaluation = SCF + analytic gradient — typically a few seconds for
small molecules and order N³ in basis size for the SCF itself. Phase
17b (analytic CPHF) replaces this with a single linear-response solve
and an order-of-magnitude speed-up.
Convention notes
----------------
* Imaginary modes (saddle points / non-stationary geometries) are
reported as **negative** wavenumbers — this matches the convention
used by Gaussian, NWChem, ORCA, and PySCF. Mathematically they're
``i * |ν|``; the negative-real encoding is just convenient for
sorting and plotting.
* Trans / rot projection is "exact zero" — when
``project_trans_rot=True`` (default) the 6 (or 5 for linear) lowest
modes are replaced with literal zeros so the user can rely on
``frequencies_cm1[6:]`` (or ``[5:]``) being the vibrational subset.
* Mass weighting uses standard atomic masses (the same table used by
:func:`vibeqc.center_of_mass`). For isotope effects, edit
``HessianFDOptions.atomic_masses_amu`` to override the table.
* Step size 0.005 bohr (~0.0026 Å) is the default. Smaller steps risk
round-off at SCF tolerance; larger steps include cubic anharmonicity.
Run a step-size convergence check if frequencies look suspicious.
* **Tighten SCF tolerances.** The default ``RHFOptions`` use
``conv_tol_grad = 1e-6``; that is fine for energies but contributes
~1e-6 / step ≈ 2e-4 Ha/bohr² of noise to each Hessian element at
the default step. For µHa-level frequency reproducibility, pass
``scf_options`` with ``conv_tol_energy = 1e-12`` and
``conv_tol_grad = 1e-9`` (the PySCF default for Hessian work).
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Optional, Sequence
import numpy as np
from . import _vibeqc_core as _core
from ._vibeqc_core import (
Atom,
BasisSet,
GridOptions,
Molecule,
RHFOptions,
RKSOptions,
UHFOptions,
UKSOptions,
compute_dipole,
compute_gradient,
compute_gradient_rks,
compute_gradient_uhf,
compute_gradient_uks,
run_rhf,
run_rks,
run_uhf,
run_uks,
)
from .properties import _ATOMIC_MASSES, _total_density, center_of_mass
# ----------------------------------------------------------------------
# Constants
# ----------------------------------------------------------------------
# CODATA 2018: 1 atomic-unit angular frequency in cm⁻¹.
# Derivation: ω in a.u. corresponds to E = ℏω with ℏ = 1, so ω[a.u.]
# numerically equals an energy in Hartree, and 1 Ha = 219474.6313632 cm⁻¹.
_HARTREE_TO_CM_INV: float = 219474.6313632
# 1 atomic mass unit in electron-mass units (so Hessian [Ha/bohr²]
# divided by mass [m_e] gives ω² [a.u.²]). CODATA 2018: m_u/m_e.
_AMU_TO_ELECTRON_MASS: float = 1822.888486209
# IR intensity conversion (Wilson-Decius-Cross formula
# I = (N_A π / (3 c²)) · |∂μ/∂Q|²
# in SI, then rescaled to km/mol).
#
# The widely-quoted 974.864 km/mol per |∂μ/∂Q|²[a.u.] factor (Gaussian,
# Halls-Schlegel JCC 1998) assumes Q is the mass-weighted normal
# coordinate with masses in *amu* (i.e. Q in √amu·bohr). Our
# mass-weighted Hessian uses *electron-mass* units throughout, so Q
# is in √m_e·bohr and |∂μ/∂Q|² differs by a factor m_u/m_e ≈ 1823.
# Multiply through to get the right factor for our convention:
_AU_TO_KM_PER_MOL: float = 974.864 * _AMU_TO_ELECTRON_MASS
# ----------------------------------------------------------------------
# Public dataclasses
# ----------------------------------------------------------------------
[docs]
@dataclass
class HessianFDOptions:
"""Knobs for :func:`compute_hessian_fd`.
step_bohr
Centered-finite-difference displacement (default 0.005 bohr).
Smaller → more sensitive to SCF round-off; larger → more
anharmonicity contamination. 0.001–0.01 bohr is the practical
sweet spot.
sym_strict
Symmetrize the FD Hessian as ``H = (H + H.T) / 2`` before
mass-weighting (default True). Set False to keep the raw FD
matrix — useful for diagnosing translational invariance of
the gradient.
project_trans_rot
Project the 3 (linear: 2) translational + 3 (linear: 2)
rotational zero modes out of the mass-weighted Hessian
before diagonalisation (default True). When True, the
corresponding entries of ``frequencies_cm1`` are exactly
zero so callers can slice ``[6:]`` (or ``[5:]`` for linear
molecules) for the vibrational subset.
atomic_masses_amu
Optional per-atom mass override (length n_atoms, in amu).
Defaults to the standard atomic masses for the molecule's
Z values. Use this for isotope substitution (D for H, ¹³C,
etc.).
"""
step_bohr: float = 0.005
sym_strict: bool = True
project_trans_rot: bool = True
atomic_masses_amu: Optional[Sequence[float]] = None
include_dipole_derivatives: bool = False
"""When True, also evaluate the dipole moment at every displaced
geometry and assemble the dipole-derivative tensor
``∂μ_β/∂R_iα`` (in atomic units, e). Cheap because the SCF is
already done — adds only one ``compute_dipole`` call per
displacement. Required by :func:`ir_intensities`. Default False."""
[docs]
@dataclass
class HessianResult:
"""Output of :func:`compute_hessian_fd`.
hessian
(3N, 3N) matrix in Hartree/bohr². Indexing convention:
row/col ``3*i + α`` is atom ``i``, Cartesian direction
``α ∈ {0, 1, 2}`` for ``{x, y, z}``.
hessian_mw
Mass-weighted Hessian: ``M^(-1/2) H M^(-1/2)`` where ``M``
is the diagonal mass matrix in electron-mass units, with
each atom's mass duplicated 3× (one per Cartesian DOF).
Eigenvalues of this matrix are ``ω²`` in atomic units.
frequencies_cm1
(3N,) array of harmonic frequencies in cm⁻¹, sorted
ascending. Imaginary modes appear as negative numbers.
Trans/rot zero modes appear as exact 0.0 when
``project_trans_rot=True``.
normal_modes
(3N, 3N) orthonormal mass-weighted normal-mode matrix
(each column is a mode). Cartesian displacement pattern
for atom ``i`` in mode ``k`` is
``normal_modes[3i:3i+3, k] / sqrt(M_i)``.
imaginary_count
Number of modes with ω² < 0 after trans/rot projection.
n_displacements
``6 * n_atoms`` for centered FD.
is_linear
True when the molecule's inertia tensor is rank-deficient
(collinear atoms) — 5 zero modes instead of 6.
"""
hessian: np.ndarray
hessian_mw: np.ndarray
frequencies_cm1: np.ndarray
normal_modes: np.ndarray
imaginary_count: int
n_displacements: int
is_linear: bool
dipole_derivatives: Optional[np.ndarray] = None
"""(3N, 3) tensor of ``∂μ_β/∂R_iα`` in atomic units (e), populated
when ``HessianFDOptions.include_dipole_derivatives=True`` (else
None). Row index ``3*i + α``; column ``β ∈ {0,1,2}`` for ``{x,y,z}``.
Used by :func:`ir_intensities` to compute IR band intensities."""
dipole_origin: Optional[np.ndarray] = None
"""(3,) array, the dipole reference origin in bohr — typically the
center of mass of the *reference* (undisplaced) geometry, which
keeps the dipole-derivative tensor origin-consistent across
displacements. None when dipole derivatives weren't requested."""
masses_amu: Optional[np.ndarray] = None
"""(n_atoms,) atomic masses in amu used during mass-weighting.
Stored so :func:`ir_intensities` and downstream thermochemistry
can stay self-contained without re-deriving masses from the
molecule."""
# ----------------------------------------------------------------------
# Internal helpers
# ----------------------------------------------------------------------
def _atom_positions_bohr(mol: Molecule) -> np.ndarray:
"""(n_atoms, 3) array of atom coordinates in bohr."""
return np.asarray(
[[a.xyz[0], a.xyz[1], a.xyz[2]] for a in mol.atoms],
dtype=np.float64,
)
def _atom_z(mol: Molecule) -> np.ndarray:
"""(n_atoms,) int array of atomic numbers."""
return np.asarray([int(a.Z) for a in mol.atoms], dtype=np.int64)
def _displaced_molecule(mol: Molecule, atom_index: int,
cart: int, delta: float) -> Molecule:
"""Return a copy of ``mol`` with one atom shifted by ``delta`` along
Cartesian axis ``cart``. ``delta`` is in bohr and may be negative."""
atoms = []
for j, a in enumerate(mol.atoms):
xyz = [a.xyz[0], a.xyz[1], a.xyz[2]]
if j == atom_index:
xyz[cart] += delta
atoms.append(Atom(int(a.Z), xyz))
return Molecule(atoms, charge=mol.charge, multiplicity=mol.multiplicity)
def _dipole_at(mol: Molecule, basis: BasisSet,
result, origin: np.ndarray) -> np.ndarray:
"""Total dipole moment (electronic + nuclear, in atomic units) of
a converged SCF result, evaluated about ``origin`` (bohr).
Mirrors :func:`vibeqc.dipole_moment` but returns a plain ``ndarray``
instead of the ``DipoleMoment`` dataclass — used inside the FD
Hessian loop where we just need the three numbers.
"""
dip = compute_dipole(basis, [float(x) for x in origin])
Mx = np.asarray(dip.x)
My = np.asarray(dip.y)
Mz = np.asarray(dip.z)
P = _total_density(result)
mu_e = -np.array([
np.einsum("ij,ji->", P, Mx),
np.einsum("ij,ji->", P, My),
np.einsum("ij,ji->", P, Mz),
])
mu_n = np.zeros(3)
for atom in mol.atoms:
z = float(atom.Z)
mu_n[0] += z * (atom.xyz[0] - origin[0])
mu_n[1] += z * (atom.xyz[1] - origin[1])
mu_n[2] += z * (atom.xyz[2] - origin[2])
return mu_e + mu_n
def _scf_and_gradient_at(
mol: Molecule,
basis_name: str,
method: str,
scf_options,
grid_options,
):
"""Run SCF + analytic gradient at ``mol``. Returns
``(basis, result, gradient)`` so the caller can reuse the
converged result (e.g. to compute the dipole moment) without a
redundant SCF."""
basis = BasisSet(mol, basis_name)
if method == "RHF":
result = run_rhf(mol, basis, scf_options)
if not result.converged:
raise RuntimeError(
"FD Hessian: SCF failed to converge at a displaced geometry "
"(method=RHF). Try a tighter scf_options or a smaller step."
)
g = compute_gradient(mol, basis, result)
elif method == "UHF":
result = run_uhf(mol, basis, scf_options)
if not result.converged:
raise RuntimeError(
"FD Hessian: UHF SCF failed to converge at a displaced "
"geometry."
)
g = compute_gradient_uhf(mol, basis, result)
elif method == "RKS":
result = run_rks(mol, basis, scf_options)
if not result.converged:
raise RuntimeError(
"FD Hessian: RKS SCF failed to converge at a displaced "
"geometry."
)
g = compute_gradient_rks(mol, basis, result, grid_options)
elif method == "UKS":
result = run_uks(mol, basis, scf_options)
if not result.converged:
raise RuntimeError(
"FD Hessian: UKS SCF failed to converge at a displaced "
"geometry."
)
g = compute_gradient_uks(mol, basis, result, grid_options)
else:
raise ValueError(
f"FD Hessian: unknown method {method!r}. "
"Expected one of {'RHF', 'UHF', 'RKS', 'UKS'}."
)
return basis, result, np.asarray(g, dtype=np.float64)
def _resolve_masses(mol: Molecule,
override: Optional[Sequence[float]]) -> np.ndarray:
"""Look up atomic masses (in amu) from the standard table, applying
a user-supplied override if given. Returns an (n_atoms,) array."""
n_atoms = len(mol.atoms)
if override is not None:
m = np.asarray(override, dtype=np.float64)
if m.shape != (n_atoms,):
raise ValueError(
f"FD Hessian: atomic_masses_amu has shape {m.shape}, "
f"expected ({n_atoms},)"
)
if np.any(m <= 0):
raise ValueError("FD Hessian: atomic masses must be positive")
return m
masses = np.empty(n_atoms, dtype=np.float64)
for i, atom in enumerate(mol.atoms):
z = int(atom.Z)
if z < len(_ATOMIC_MASSES) and _ATOMIC_MASSES[z] > 0:
masses[i] = _ATOMIC_MASSES[z]
else:
raise ValueError(
f"FD Hessian: no built-in mass for Z={z}. Pass an explicit "
"atomic_masses_amu override."
)
return masses
def _build_trans_rot_modes(positions_bohr: np.ndarray,
masses_e: np.ndarray) -> tuple[np.ndarray, bool]:
"""Build orthonormal mass-weighted translation + rotation vectors,
one column per zero mode. Detects linear molecules (rank-deficient
rotational subspace) and emits 5 modes instead of 6.
Convention: vectors live in mass-weighted Cartesian space, i.e.
if ``v`` is a column then ``v[3i:3i+3]`` is atom ``i``'s
contribution scaled by ``sqrt(m_i)``. Translation/rotation in
*real* coordinates corresponds to no change in potential energy
by Galilean / SO(3) invariance, so these are exact zero modes
of the mass-weighted Hessian.
Returns (modes, is_linear) where modes has shape (3N, n_zero) and
columns are orthonormal.
"""
n_atoms = positions_bohr.shape[0]
n_dof = 3 * n_atoms
sqrt_m = np.sqrt(masses_e) # (n_atoms,)
# --- Translations ---------------------------------------------------
trans = np.zeros((n_dof, 3), dtype=np.float64)
for alpha in range(3):
for i in range(n_atoms):
trans[3 * i + alpha, alpha] = sqrt_m[i]
# --- Rotations (about the center of mass) --------------------------
total_mass = masses_e.sum()
com = (masses_e[:, None] * positions_bohr).sum(axis=0) / total_mass
R = positions_bohr - com # (n_atoms, 3) coords relative to COM
rot = np.zeros((n_dof, 3), dtype=np.float64)
# Generator R_α x: column for rotation about α has components
# (R_i × ê_α) at atom i, scaled by sqrt(m_i).
e = np.eye(3)
for alpha in range(3):
cross = np.cross(R, e[alpha]) # (n_atoms, 3)
for i in range(n_atoms):
rot[3 * i:3 * i + 3, alpha] = sqrt_m[i] * cross[i]
raw = np.concatenate([trans, rot], axis=1) # (n_dof, 6)
# Orthonormalise via QR; drop columns whose norm collapsed below a
# tolerance (linear-molecule case: one rotational mode is null).
Q, R_qr = np.linalg.qr(raw)
diag = np.abs(np.diag(R_qr))
# Tolerance scaled to the largest column so it works for any size
keep = diag > 1e-8 * diag.max()
modes = Q[:, keep]
is_linear = modes.shape[1] == 5
if modes.shape[1] not in (5, 6):
# Single-atom or weirdly degenerate case — fall through with
# whatever survived. Hard to reach in practice for n_atoms ≥ 2.
pass
return modes, is_linear
def _project_out(H_mw: np.ndarray, zero_modes: np.ndarray) -> np.ndarray:
"""Return P · H_mw · P with P = I − Σ_k v_k v_k^T (zero modes
orthonormal)."""
n = H_mw.shape[0]
P = np.eye(n) - zero_modes @ zero_modes.T
return P @ H_mw @ P
def _enforce_projected_zero(eigenvalues: np.ndarray,
eigenvectors: np.ndarray,
n_zero: int) -> tuple[np.ndarray, np.ndarray]:
"""After projection, the lowest ``n_zero`` eigenvalues should be
near-zero (numerical noise of the projector). Replace them with
exact zeros so callers can slice them off cleanly. Eigenvalues
are returned in ascending order — pull the n_zero with smallest
|value| (not most-negative; small noise can be of either sign)."""
abs_idx = np.argsort(np.abs(eigenvalues))
zero_idx = abs_idx[:n_zero]
eigenvalues = eigenvalues.copy()
eigenvalues[zero_idx] = 0.0
return eigenvalues, eigenvectors
def _omega2_to_cm_inv(omega2: np.ndarray) -> np.ndarray:
"""Convert eigenvalues of the mass-weighted Hessian (in atomic units
with masses in m_e) to wavenumbers in cm⁻¹. Imaginary modes
(negative eigenvalues) come out negative."""
out = np.empty_like(omega2)
pos = omega2 >= 0
out[pos] = np.sqrt(omega2[pos]) * _HARTREE_TO_CM_INV
out[~pos] = -np.sqrt(-omega2[~pos]) * _HARTREE_TO_CM_INV
return out
# ----------------------------------------------------------------------
# Public driver
# ----------------------------------------------------------------------
def compute_hessian_fd(
mol: Molecule,
basis_name: str,
method: str = "RHF",
*,
scf_options=None,
grid_options: Optional[GridOptions] = None,
hessian_options: Optional[HessianFDOptions] = None,
) -> HessianResult:
"""Build the 3N × 3N nuclear Hessian by central differences on the
analytic atomic gradient and return mass-weighted normal-mode
frequencies.
See module docstring for API examples and convention notes.
Parameters
----------
mol :
Molecule at the geometry where the Hessian is wanted.
Should be a stationary point for the resulting frequencies
to be physical, but the function doesn't enforce this so
non-equilibrium force constants can be evaluated too.
basis_name :
Basis-set name. Re-built per displacement so basis-center
derivative pieces are absorbed into the FD.
method :
One of ``"RHF"``, ``"UHF"``, ``"RKS"``, ``"UKS"``
(case-insensitive).
scf_options :
Per-method options struct. Default-initialises a fresh
``RHFOptions`` / ``UHFOptions`` / ``RKSOptions`` /
``UKSOptions`` if None. For DFT methods the functional
comes from this struct's ``functional`` field.
grid_options :
DFT integration grid. Forwarded to
:func:`compute_gradient_rks` / ``_uks``. Defaults to a
fresh ``GridOptions()`` for DFT methods (same default as
the SCF call), unused for HF.
hessian_options :
Step size + projection knobs.
"""
if hessian_options is None:
hessian_options = HessianFDOptions()
method_norm = method.upper()
if method_norm not in {"RHF", "UHF", "RKS", "UKS"}:
raise ValueError(
f"FD Hessian: unknown method {method!r}. "
"Expected one of 'RHF', 'UHF', 'RKS', 'UKS'."
)
if scf_options is None:
scf_options = {
"RHF": RHFOptions(),
"UHF": UHFOptions(),
"RKS": RKSOptions(),
"UKS": UKSOptions(),
}[method_norm]
if grid_options is None:
grid_options = GridOptions()
# ---- Build the Hessian via centered FD ---------------------------------
n_atoms = len(mol.atoms)
if n_atoms < 1:
raise ValueError("FD Hessian: empty molecule")
n_dof = 3 * n_atoms
H = np.zeros((n_dof, n_dof), dtype=np.float64)
delta = float(hessian_options.step_bohr)
if delta <= 0:
raise ValueError("FD Hessian: step_bohr must be positive")
# Optional dipole-derivative tensor. Reference origin = COM of the
# undisplaced geometry; using the same origin for all dipole
# evaluations keeps the dipole-derivative tensor free of spurious
# origin-shift contributions (which matter for charged systems and
# leak into IR intensities even for neutral ones if you mix origins).
include_dipole = bool(hessian_options.include_dipole_derivatives)
if include_dipole:
dipole_origin = center_of_mass(mol)
dipole_derivs = np.zeros((n_dof, 3), dtype=np.float64)
else:
dipole_origin = None
dipole_derivs = None
n_disp = 0
for i in range(n_atoms):
for cart in range(3):
mol_plus = _displaced_molecule(mol, i, cart, +delta)
mol_minus = _displaced_molecule(mol, i, cart, -delta)
basis_p, res_p, g_plus = _scf_and_gradient_at(
mol_plus, basis_name, method_norm,
scf_options, grid_options,
)
basis_m, res_m, g_minus = _scf_and_gradient_at(
mol_minus, basis_name, method_norm,
scf_options, grid_options,
)
n_disp += 2
# Column 3i+cart of H = (g(+δ) - g(-δ)) / (2δ), flattened.
col = (g_plus - g_minus).reshape(-1) / (2.0 * delta)
H[:, 3 * i + cart] = col
if include_dipole:
mu_p = _dipole_at(mol_plus, basis_p, res_p, dipole_origin)
mu_m = _dipole_at(mol_minus, basis_m, res_m, dipole_origin)
dipole_derivs[3 * i + cart, :] = (mu_p - mu_m) / (2.0 * delta)
if hessian_options.sym_strict:
H = 0.5 * (H + H.T)
# ---- Mass weight ------------------------------------------------------
masses_amu = _resolve_masses(mol, hessian_options.atomic_masses_amu)
masses_e = masses_amu * _AMU_TO_ELECTRON_MASS
inv_sqrt_m_per_dof = np.repeat(1.0 / np.sqrt(masses_e), 3)
M = inv_sqrt_m_per_dof # diagonal of M^(-1/2) in 3N space
H_mw = (M[:, None] * H) * M[None, :]
H_mw = 0.5 * (H_mw + H_mw.T) # squash residual asymmetry from FD noise
# ---- Project out trans/rot, diagonalize -------------------------------
is_linear = False
if hessian_options.project_trans_rot:
positions = _atom_positions_bohr(mol)
zero_modes, is_linear = _build_trans_rot_modes(positions, masses_e)
H_proj = _project_out(H_mw, zero_modes)
omega2, modes = np.linalg.eigh(0.5 * (H_proj + H_proj.T))
n_zero = zero_modes.shape[1]
omega2, modes = _enforce_projected_zero(omega2, modes, n_zero)
else:
omega2, modes = np.linalg.eigh(H_mw)
# eigh returns ascending eigenvalues already, but the
# "force smallest |value| to zero" step may have reordered them
# via index reuse. Re-sort to keep the contract clean.
order = np.argsort(omega2)
omega2 = omega2[order]
modes = modes[:, order]
freqs = _omega2_to_cm_inv(omega2)
imag_count = int(np.sum(omega2 < -1e-10))
return HessianResult(
hessian=H,
hessian_mw=H_mw,
frequencies_cm1=freqs,
normal_modes=modes,
imaginary_count=imag_count,
n_displacements=n_disp,
is_linear=is_linear,
dipole_derivatives=dipole_derivs,
dipole_origin=(np.asarray(dipole_origin)
if dipole_origin is not None else None),
masses_amu=masses_amu,
)
def ir_intensities(result: HessianResult) -> np.ndarray:
"""Convert a Hessian result's dipole derivatives + normal modes
into IR band intensities in km/mol.
Phase 17a-2.
The harmonic-oscillator IR intensity for the *p*-th normal mode
is the Wilson-Decius-Cross formula
.. math::
I_p = \\frac{N_A \\pi}{3 c^2} \\left| \\frac{\\partial \\mu}{\\partial Q_p} \\right|^2,
where :math:`Q_p` is the mass-weighted normal coordinate of mode
*p* with eigenvector :math:`\\mathbf{L}_p` of the mass-weighted
Hessian. Going from Cartesian to normal-mode dipole derivatives
uses the chain rule
.. math::
\\frac{\\partial \\mu_\\beta}{\\partial Q_p} = \\sum_k
\\frac{\\partial \\mu_\\beta}{\\partial R_k}
\\cdot \\frac{L_{kp}}{\\sqrt{M_k}},
with :math:`M_k` the mass (in electron masses) of the atom
corresponding to flat index *k*. The conversion from atomic units
of :math:`|\\partial \\mu / \\partial Q|^2` to km/mol uses the
standard Gaussian / NWChem / ORCA / PySCF factor of 974.864
(CODATA-derived).
Parameters
----------
result :
:class:`HessianResult` produced with
``HessianFDOptions(include_dipole_derivatives=True)``.
Returns
-------
np.ndarray
``(3N,)`` array of IR intensities in km/mol, indexed by the
same ordering as ``result.frequencies_cm1`` (ascending
frequency). Trans/rot zero modes have intensities ≈ 0 by
construction (the corresponding eigenvector has no overlap
with the dipole-gradient field for translation-invariant
observables).
Raises
------
ValueError
If ``result.dipole_derivatives`` is None (the Hessian was
built without ``include_dipole_derivatives=True``).
"""
if result.dipole_derivatives is None:
raise ValueError(
"ir_intensities: dipole derivatives missing on the "
"HessianResult. Re-run compute_hessian_fd with "
"HessianFDOptions(include_dipole_derivatives=True)."
)
if result.masses_amu is None:
raise ValueError(
"ir_intensities: masses_amu missing on the HessianResult "
"(should be populated by compute_hessian_fd)."
)
B = np.asarray(result.dipole_derivatives, dtype=np.float64) # (3N, 3)
L = np.asarray(result.normal_modes, dtype=np.float64) # (3N, 3N)
masses_e = np.asarray(result.masses_amu, dtype=np.float64) * _AMU_TO_ELECTRON_MASS
inv_sqrt_m_per_dof = np.repeat(1.0 / np.sqrt(masses_e), 3) # (3N,)
# Mass-weighted Cartesian dipole derivative tensor: B_mw[k, β] =
# (∂μ_β/∂R_k) / sqrt(M_k). Then transform to normal-mode basis:
# dμ/dQ[β, p] = Σ_k B_mw[k, β] · L[k, p].
B_mw = B * inv_sqrt_m_per_dof[:, None]
dmu_dQ = B_mw.T @ L # (3, 3N)
intensities = _AU_TO_KM_PER_MOL * np.sum(dmu_dQ ** 2, axis=0) # (3N,)
# Mask out zero-frequency (trans / rot) modes. These modes don't
# oscillate, so the harmonic IR intensity formula doesn't apply to
# them — they're not "vibrations". Numerically, the eigh
# decomposition of the trans/rot-projected Hessian returns an
# arbitrary orthonormal basis of the kernel; that basis can mix
# translations (zero ∂μ/∂Q by sum rule) with rotations (nonzero
# ∂μ/∂Q because rotating a polar molecule changes the dipole
# direction), giving a spurious "intensity" for each kernel
# eigenvector. Standard convention (Gaussian, NWChem, ORCA): drop
# them. Threshold at 1 cm⁻¹ — well below any real vibration but
# well above floating-point noise on the projected eigenvalues.
freqs = np.asarray(result.frequencies_cm1, dtype=np.float64)
intensities = np.where(np.abs(freqs) < 1.0, 0.0, intensities)
return intensities
__all__ = [
"HessianFDOptions",
"HessianResult",
"compute_hessian_fd",
"ir_intensities",
]