"""Phase 17a-1 -- molecular finite-difference Hessian + harmonic frequencies.
Two-sided finite differences on the atomic gradient give the
3N x 3N nuclear Hessian in Hartree/bohr^2. 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). OM2/OM3
use the existing checked native energy differences inside this outer stencil.
The initial inner step is 0.001 bohr; the native guard may refine individual
components by factors of 0.1 up to three times. Actual accepted inner steps
are not exposed by the gradient API. Their cost is quadratic in Cartesian degrees of
freedom times the SCF cost; IR intensities are unavailable. These Hessians
describe the implemented OMx approximations, not published-method parity.
For Gaussian methods, each
evaluation = SCF + analytic gradient -- typically a few seconds for
small molecules and order N^3 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 inserts exact zeros when
``project_trans_rot=True`` (default). Negative vibrational frequencies
sort before these rigid-body zeros; do not discard the first six
frequencies to select vibrations at a saddle point.
* 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^2 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).
* **ECPs and the JK backend follow the SCF options.** The gradient that
is differentiated is built from ``scf_options`` unless an explicit
``gradient_options`` is passed: density fitting / COSX and the ECP
fields of both routes (libecpint XML-library ``ecp_centers`` and
inline ``ecp_primitive_blocks``) are mirrored, and the ECP centres --
absolute coordinates on the options -- are moved with each displaced
atom (#576). Pre-fix the default gradient was the bare-Z all-electron
one and the centres stayed behind, so an ECP-bearing atom's SCF
diverged at the first displacement while the all-electron atoms'
columns were silently wrong. The analytic Hessian kernels
(``compute_hessian_*_analytic``) still take no ECP options; use this
FD path for ECP systems.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Optional, Sequence
import numpy as np
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_uks,
)
from .ecp_metadata import (
attach_inline_ecp_options_from_basis_sidecar,
ecp_centre_atom_indices,
ecp_centres_follow_atoms,
effective_nuclear_charges,
options_carry_ecp,
)
from .density_fitting import resolve_default_jk_aux
from .gradient_options import gradient_options_from_scf
from .properties import _ATOMIC_MASSES, _total_density, center_of_mass
# ----------------------------------------------------------------------
# Constants
# ----------------------------------------------------------------------
# CODATA 2018: 1 atomic-unit angular frequency in cm⁻¹.
# Derivation: w in a.u. corresponds to E = ℏw with ℏ = 1, so w[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^2]
# divided by mass [m_e] gives w^2 [a.u.^2]). CODATA 2018: m_u/m_e.
_AMU_TO_ELECTRON_MASS: float = 1822.888486209
# IR intensity conversion (Wilson-Decius-Cross formula
# I = (N_A pi / (3 c^2)) . |dmu/dQ|^2
# in SI, then rescaled to km/mol).
#
# The widely-quoted 974.864 km/mol per |dmu/dQ|^2[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 |dmu/dQ|^2 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, ¹^3C,
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
``dmu_b/dR_ia`` (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."""
frozen_indices: Optional[Sequence[int]] = None
"""Atom indices to hold fixed (not displaced). When set, the result
is the **partial-Hessian** sub-block over the *unfrozen* atoms only --
the ``(3M, 3M)`` block of the full Hessian for ``M`` unfrozen atoms --
and its ``3M`` frequencies are all vibrational: the frozen atoms
anchor the system, so there are no translational/rotational zero
modes and ``project_trans_rot`` is forced off. This is the
partial-Hessian vibrational analysis (PHVA; Li & Jensen, *Theor.
Chem. Acc.* 2002) used to characterise a transition state on a slab
while holding the bulk-like substrate layers fixed -- it cuts the
displacement count from ``6N`` to ``6M``. ``None`` => every atom is
displaced (the full Hessian). Downstream
:func:`vibeqc.compute_thermochemistry` treats a partial-Hessian
result as a vibrational-only (anchored) system: no gas-phase
translational / rotational partition functions."""
[docs]
@dataclass
class HessianResult:
"""Output of :func:`compute_hessian_fd`.
hessian
(3N, 3N) matrix in Hartree/bohr^2. Indexing convention:
row/col ``3*i + a`` is atom ``i``, Cartesian direction
``a in {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 3x (one per Cartesian DOF).
Eigenvalues of this matrix are ``w^2`` 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 w^2 < -1e-10 (atomic units), the existing
summary threshold. The full signed frequencies retain smaller
negative values too; the count is not a stability certificate.
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 ``dmu_b/dR_ia`` in atomic units (e), populated
when ``HessianFDOptions.include_dipole_derivatives=True`` (else
None). Row index ``3*i + a``; column ``b in {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 -- for
a partial Hessian (frozen atoms) this is the ``M`` unfrozen masses,
matching the ``3M`` rows/columns. Stored so :func:`ir_intensities`
and downstream thermochemistry can stay self-contained without
re-deriving masses from the molecule."""
solvent: Optional[str] = None
"""``None`` for a gas-phase Hessian. For a Hessian of the in-solvent
surface, the label of the medium (for example ``"CPCM water (eps = 78.39)"``):
the displaced gradients were the solvated analytic gradients, with the
apparent surface charges re-converged at every displacement, so the
frequencies are those of the *equilibrium*-solvation surface."""
unfrozen_indices: Optional[np.ndarray] = None
"""``None`` for a full Hessian. When ``HessianFDOptions.frozen_indices``
was set, the ``(M,)`` array of *unfrozen* atom indices this partial
Hessian's ``3M`` rows/columns correspond to (row/col ``3*k + a`` is
unfrozen atom ``unfrozen_indices[k]``, Cartesian ``a``). Lets callers
map normal modes back to atoms and signals
:func:`vibeqc.compute_thermochemistry` to use the vibrational-only
(anchored-system) partition function."""
# ----------------------------------------------------------------------
# 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,
nuclear_charges: 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.
``nuclear_charges`` is the per-atom charge the valence-only Hamiltonian
sees (``Z - n_core`` on ECP atoms, bare ``Z`` otherwise;
:func:`vibeqc.ecp_metadata.effective_nuclear_charges`). It is the same
vector for every displaced geometry, so the caller resolves it once on
the reference molecule. With bare ``Z`` on an ECP system each ECP
atom's dipole-derivative row picked up a spurious ``n_core`` per
Cartesian direction and leaked into every IR intensity of a mode that
moves that atom (GitLab #642).
"""
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 z, atom in zip(nuclear_charges, mol.atoms):
z = float(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,
gradient_options=None,
*,
reference_mol: Optional[Molecule] = None,
energy_only: bool = False,
solvent=None,
):
"""Run SCF + analytic gradient at ``mol``. Returns
``(basis, result, gradient)`` (gradient is None for energy_only) so the caller can reuse the
converged result (e.g. to compute the dipole moment) without a
redundant SCF.
The gradient differentiates the Hamiltonian the SCF at ``mol``
solves (#576):
* ``gradient_options=None`` derives the ``GradientOptions`` from
``scf_options`` -- JK backend (density_fit / cosx / aux_basis) and
ECP fields, both routes -- through
:func:`vibeqc.gradient_options.gradient_options_from_scf`. An
explicit object is used as given, except that its ECP fields are
synchronised to the SCF options for the call.
* When ``reference_mol`` is given, ``scf_options`` describes that
geometry and ``mol`` is a displaced copy of it: the ECP centres
(``ecp_centers`` / ``ecp_primitive_centers``, which are absolute
coordinates) are moved onto the displaced atoms for the duration
of the call and restored afterwards
(:func:`vibeqc.ecp_metadata.ecp_centres_follow_atoms`). Without
this the displaced atom silently reverts to bare Z while V_ECP
keeps acting and the SCF diverges.
"""
with ecp_centres_follow_atoms(
scf_options, reference_mol, mol, gradient_options
):
if gradient_options is None:
# The displaced SCFs below call the bare drivers, so give a
# fitted request its default JK aux first (#362); the gradient
# options derived next then mirror the resolved name.
resolve_default_jk_aux(scf_options, basis_name)
gradient_options = gradient_options_from_scf(scf_options)
# Preloaded READ matrices belong to the reference geometry. Project
# from that fixed basis at every displacement and restore them after
# the call, so projections do not accumulate along the FD stencil.
restored = {}
if reference_mol is not None and scf_options is not None:
from ._initial_guess import coerce_initial_guess
from ._vibeqc_core import InitialGuess
if coerce_initial_guess(scf_options.initial_guess) == InitialGuess.READ:
from .guess_read import _to_current_basis
source_basis = BasisSet(reference_mol, basis_name)
target_basis = BasisSet(mol, basis_name)
for name in ("read_density", "read_density_alpha", "read_density_beta"):
value = getattr(scf_options, name, None)
if value is not None and np.asarray(value).size:
restored[name] = np.asarray(value).copy()
setattr(scf_options, name, _to_current_basis(
value, source_basis, target_basis,
))
try:
return _scf_and_gradient_at_impl(
mol, basis_name, method, scf_options, grid_options,
gradient_options, energy_only=energy_only, solvent=solvent,
)
finally:
for name, value in restored.items():
setattr(scf_options, name, value)
def _solvated_scf_and_gradient_at(
mol, basis, basis_name, method, scf_options, grid_options,
gradient_options, solvent, energy_only,
):
"""The in-solvent counterpart of the gas-phase body of
:func:`_scf_and_gradient_at_impl`: ``(basis, result, gradient)`` on the
surface the solvated job reports.
A Hessian that finite-differences the gas-phase gradient after a solvated
SCF describes the wrong surface and says nothing: every number is smooth
and plausible (#359 was the same defect for optimizations). So each
displaced geometry gets its own :func:`run_cpcm_scf`, whose converged charges
the matching :func:`cpcm_gradient` differentiates.
"""
from .solvation import cpcm_gradient, run_cpcm_scf
from .solvation.driver import _solvent_aware_scf_result
lower = method.lower()
if lower not in ("rhf", "uhf", "rks", "uks"):
raise NotImplementedError(
f"FD Hessian in solvent: method {method!r} has no solvated analytic "
"gradient (RHF, UHF, RKS and UKS do)."
)
sol = run_cpcm_scf(mol, basis, method=lower, solvent=solvent, options=scf_options)
if not sol.converged:
raise RuntimeError(
f"FD Hessian: the solvated SCF ({method}, {sol.solvent_name}) failed "
"to converge at a displaced geometry. Try tighter scf_options or a "
"smaller step."
)
result = _solvent_aware_scf_result(sol)
if not np.isfinite(float(sol.energy)):
raise ValueError("FD Hessian: non-finite in-solvent energy at a displacement.")
if energy_only:
return basis, result, None
if lower in ("rks", "uks"):
from .gradient_terms import functional_gradient_terms_missing
if functional_gradient_terms_missing(
scf_options.functional, spin=1 if lower == "rks" else 2
):
g = _full_energy_ks_gradient(
mol, basis_name, method, scf_options, grid_options,
gradient_options, solvent=solvent,
)
else:
g = cpcm_gradient(
sol.scf, mol, basis, sol, method=lower,
options=gradient_options, grid_options=grid_options,
)
else:
g = cpcm_gradient(
sol.scf, mol, basis, sol, method=lower, options=gradient_options,
)
return basis, result, np.asarray(g, dtype=np.float64)
def _scf_and_gradient_at_impl(
mol: Molecule,
basis_name: str,
method: str,
scf_options,
grid_options,
gradient_options,
*,
energy_only: bool = False,
solvent=None,
):
"""Body of :func:`_scf_and_gradient_at`; ``gradient_options`` is a
concrete ``GradientOptions`` already matched to ``scf_options``.
A resolved ``solvent`` (:class:`vibeqc.SolventModel`) puts the evaluation on
the in-solvent surface: the SCF is :func:`~vibeqc.solvation.run_cpcm_scf`
and the gradient :func:`~vibeqc.solvation.cpcm_gradient`, both at this
displaced geometry, with the apparent charges re-converged there.
"""
# Public wrappers apply selector validation and project file-backed READ
# sources onto each displaced geometry before entering the native driver.
from . import run_rhf, run_uhf, run_rks, run_uks
basis = BasisSet(mol, basis_name)
if solvent is not None:
return _solvated_scf_and_gradient_at(
mol, basis, basis_name, method, scf_options, grid_options,
gradient_options, solvent, energy_only,
)
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, gradient_options)
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, gradient_options)
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."
)
if energy_only:
if not np.isfinite(result.energy):
raise ValueError("FD Hessian: non-finite full-energy displacement.")
return basis, result, None
from .gradient_terms import functional_gradient_terms_missing
if functional_gradient_terms_missing(scf_options.functional, spin=1):
g = _full_energy_ks_gradient(
mol, basis_name, method, scf_options, grid_options, gradient_options,
)
else:
g = compute_gradient_rks(mol, basis, result, grid_options,
gradient_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."
)
if energy_only:
if not np.isfinite(result.energy):
raise ValueError("FD Hessian: non-finite full-energy displacement.")
return basis, result, None
from .gradient_terms import functional_gradient_terms_missing
if functional_gradient_terms_missing(scf_options.functional, spin=2):
g = _full_energy_ks_gradient(
mol, basis_name, method, scf_options, grid_options, gradient_options,
)
else:
g = compute_gradient_uks(mol, basis, result, grid_options,
gradient_options)
elif method == "ROHF":
# Restricted open-shell HF: analytic gradient with the ROHF
# energy-weighted density (vibeqc.rohf.compute_rohf_gradient).
from .rohf import ROHFOptions, compute_rohf_gradient, run_rohf
result = run_rohf(
mol, basis, scf_options if scf_options is not None else ROHFOptions()
)
if not result.converged:
raise RuntimeError(
"FD Hessian: ROHF SCF failed to converge at a displaced "
"geometry."
)
g = compute_rohf_gradient(
mol, basis, result, gradient_options=gradient_options
)
else:
raise ValueError(
f"FD Hessian: unknown method {method!r}. "
"Expected one of {'RHF', 'UHF', 'RKS', 'UKS', 'ROHF'}."
)
return basis, result, np.asarray(g, dtype=np.float64)
# Inner central-energy step follows the reported nested-FD Hessian convention.
_KS_ENERGY_GRADIENT_STEP_BOHR = 1e-3
def _full_energy_ks_gradient(mol, basis_name, method, scf_options, grid_options,
gradient_options, solvent=None):
"""Differentiate the complete SCF energy, preserving ECP/READ transport.
Larsen et al., J. Phys.: Condens. Matter 29, 273002 (2017), sec.8.1:
the outer Hessian differentiates gradients; here each missing analytic
gradient is itself a two-sided energy difference. No incomplete kernel is
called, and every inner SCF must converge on the same requested surface.
"""
step = _KS_ENERGY_GRADIENT_STEP_BOHR
gradient = np.zeros((len(mol.atoms), 3), dtype=np.float64)
for atom in range(len(mol.atoms)):
for axis in range(3):
energies = []
for sign in (-1, 1):
moved = _displaced_molecule(mol, atom, axis, sign * step)
# ``solvent`` is passed only when there is one, so the gas-phase
# call is exactly the one it was (tests pin its keywords).
_, result, _ = _scf_and_gradient_at(
moved, basis_name, method, scf_options, grid_options,
gradient_options, reference_mol=mol, energy_only=True,
**({} if solvent is None else {"solvent": solvent}),
)
energies.append(float(result.energy))
gradient[atom, axis] = (energies[1] - energies[0]) / (2 * step)
return gradient
def _omx_and_gradient_at(mol: Molecule, method: str):
"""Evaluate the implemented OM2/OM3 surface, without a Gaussian basis.
The dispatcher retains the native energy/parameter snapshot used by
its checked numerical gradient. Inner displaced SCF failures propagate;
this is a nested finite difference, not an analytic OMx derivative.
"""
from .semiempirical.runner import run_semiempirical
result = run_semiempirical(method.lower(), mol)
if not bool(result.converged):
raise RuntimeError(
f"FD Hessian: {method} SCF failed at a displaced geometry."
)
if not np.isfinite(result.energy):
raise ValueError("FD Hessian: non-finite OMx energy.")
gradient = np.asarray(result.gradient(), dtype=np.float64)
if gradient.shape != (len(mol.atoms), 3):
raise ValueError("FD Hessian: OMx gradient must have shape (n_atoms, 3).")
if not np.all(np.isfinite(gradient)):
raise ValueError("FD Hessian: non-finite OMx gradient.")
return None, result, gradient
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_a x: column for rotation about a has components
# (R_i x ê_a) 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 - S_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,
grid_level: str = "orca-defgrid3",
gradient_options=None,
hessian_options: Optional[HessianFDOptions] = None,
dispersion_params=None,
source=None,
solvent=None,
) -> HessianResult:
"""Build the 3N x 3N nuclear Hessian by central differences on the
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. Use an empty string
for basis-free OM2/OM3.
method :
One of ``"RHF"``, ``"UHF"``, ``"ROHF"``, ``"RKS"``, ``"UKS"``
or molecular ``"OM2"`` / ``"OM3"`` (case-insensitive). OMx uses
the run_job dispatcher defaults and retained-result numerical
gradients; Gaussian SCF/grid/gradient options and dipole derivatives
are refused. Charge, multiplicity and element gates are unchanged.
scf_options :
Per-method options struct. Default-initialises a fresh
``RHFOptions`` / ``UHFOptions`` / ``ROHFOptions`` /
``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
the resolved SCF grid for DFT methods, unused for HF. An explicit
grid is used by both the displaced SCFs and their gradients.
grid_level :
Preset for an untouched SCF grid (default ``"orca-defgrid3"``).
Pass ``"legacy"`` for construction defaults; custom grid fields win.
gradient_options :
``GradientOptions`` for the analytic gradient that is
differentiated. ``None`` (the default) derives it from
``scf_options`` -- JK backend (density_fit / aux_basis / cosx)
and ECP fields, both the XML-library and the inline-primitive
route -- so the differentiated gradient describes the same
Hamiltonian as the SCF energy (#576). Pass an explicit object
only to override the JK backend; its ECP fields are always
re-synchronised to ``scf_options`` at every displaced geometry.
On ECP systems the ECP centres carry absolute coordinates and
are moved with the displaced atom; a centre that coincides with
no atom of ``mol`` is refused up front (``ValueError``).
hessian_options :
Step size + projection knobs.
dispersion_params :
A :class:`~vibeqc.dispersion.D3BJParams` or
:class:`~vibeqc.dispersion.D4DispersionSpec`. Its analytic gradient is
recomputed at every displaced geometry and added to the SCF gradient,
so the Hessian differentiates the complete dispersion-corrected
gradient (#548). ``None`` (default) differentiates the bare SCF
surface. Refused for OM2/OM3.
source :
An energy/gradient source such as
:class:`vibeqc.manta_geometry.MolecularHFGeometrySource`: a callable
``source(molecule) -> (energy, gradient)`` exposing ``last_result``
and ``last_basis`` for the geometry it just evaluated, and
``options`` (its SCF options). The displaced gradients come from
the source instead of the legacy SCF/gradient drivers, so a
selected integral engine's frequencies differentiate the same
analytic gradient its optimizer walks; the source owns the ECP
centre following and any dispersion fold, so ``gradient_options``,
``grid_options`` and ``dispersion_params`` are refused with it.
``scf_options`` defaults to ``source.options``.
solvent :
Anything :func:`vibeqc.solvation.resolve_solvent` accepts (a preset
name, a dielectric constant, a :class:`vibeqc.SolventModel`). ``None``
or the gas phase is the ordinary gas-phase Hessian. Otherwise the
displaced gradients are the in-solvent analytic gradients
(:func:`~vibeqc.solvation.cpcm_gradient`), each at its own re-converged
apparent charges, so the frequencies are those of the
*equilibrium*-solvation surface: the solvent follows every displacement
(no non-equilibrium or dynamic response). The geometry should be a
stationary point of that same surface, so optimize in the solvent
first. Supported for ``RHF``, ``UHF``, ``RKS`` and ``UKS``; a
Direct COSMO-RS model is refused, because its energy has a cusp at a
symmetric geometry that a Hessian from gradient differences would
straddle (see the user guide, "Symmetric solutes"). Dipole derivatives,
when requested, are the solute's dipole at the in-solvent density,
without the solvent's polarization response or a local-field
correction.
"""
if hessian_options is None:
hessian_options = HessianFDOptions()
method_norm = method.upper()
from .output.formats.vibrational import FD_HESSIAN_SURFACES
if method_norm.lower() not in FD_HESSIAN_SURFACES:
raise ValueError(
f"FD Hessian: unknown method {method!r}. "
"Expected one of 'RHF', 'UHF', 'ROHF', 'RKS', 'UKS', 'OM2', 'OM3'."
)
is_omx = method_norm in {"OM2", "OM3"}
solvent_model = None
if solvent is not None:
from .solvation.driver import resolve_solvent
solvent_model = resolve_solvent(solvent)
if solvent_model is not None and solvent_model.is_gas_phase:
solvent_model = None # a vacuum "solvent" is the gas-phase Hessian
if solvent_model is not None:
if method_norm not in {"RHF", "UHF", "RKS", "UKS"}:
raise NotImplementedError(
f"FD Hessian in solvent: method {method!r} has no solvated "
"analytic gradient; RHF, UHF, RKS and UKS do. Computing it on "
"the gas-phase surface would describe a different molecule "
"than the solvated job reports."
)
if source is not None:
raise ValueError(
"FD Hessian: a gradient source (selected integral engine) and "
"a solvent do not compose; the solvated gradient is evaluated "
"by the solvation driver."
)
if solvent_model.variant == "dcosmo-rs":
raise NotImplementedError(
"FD Hessian in solvent: Direct COSMO-RS (variant='dcosmo-rs') "
"is refused. Its energy has a cusp at a symmetric geometry (the "
"COSMO FINE cavity's molecule-fixed frame ties there), so a "
"Hessian from gradient differences straddles it: a measured "
"slope change of 1.4e-3 Ha/bohr across a 0.01 bohr stencil is "
"an error of the order of 0.14 Ha/bohr^2 (1.4e-3/0.01) on a "
"curvature of about 0.4. Use variant='cpcm' or 'cosmo' for "
"frequencies."
)
_warn_if_frame_tie(solvent_model, mol, hessian_options.step_bohr)
if source is not None:
if is_omx:
raise ValueError("FD Hessian: a gradient source applies to Gaussian-basis SCF surfaces, not to OM2/OM3.")
if gradient_options is not None or grid_options is not None:
raise ValueError("FD Hessian: the gradient source owns its gradient; gradient_options and grid_options do not apply.")
if dispersion_params is not None:
raise ValueError("FD Hessian: fold dispersion into the gradient source (dispersion_params=) instead of passing it here.")
if scf_options is None:
scf_options = getattr(source, "options", None)
if is_omx and dispersion_params is not None:
raise ValueError(
"FD Hessian: dispersion_params applies to Gaussian-basis SCF "
"surfaces, not to OM2/OM3."
)
if is_omx:
if basis_name or any(
value is not None for value in (scf_options, grid_options, gradient_options)
):
raise ValueError(
"FD Hessian: OM2/OM3 require an empty basis name and use the "
"molecular dispatcher defaults, without Gaussian options."
)
if hessian_options.include_dipole_derivatives:
raise NotImplementedError("FD Hessian: OM2/OM3 dipole derivatives are unavailable.")
else:
if scf_options is None:
if method_norm == "ROHF":
from .rohf import ROHFOptions
scf_options = ROHFOptions()
else:
scf_options = {
"RHF": RHFOptions(),
"UHF": UHFOptions(),
"RKS": RKSOptions(),
"UKS": UKSOptions(),
}[method_norm]
if method_norm in {"RKS", "UKS"}:
from .runner import apply_ks_grid_default
if grid_options is not None:
scf_options.grid = grid_options
else:
apply_ks_grid_default(scf_options, grid_level)
grid_options = scf_options.grid
elif grid_options is None:
grid_options = GridOptions()
# A direct caller commonly supplies only an ECP-paired basis name. Attach
# that sidecar before deriving GradientOptions or entering the displacement
# loop; otherwise the first displaced SCF would discover the ECP too late
# and its gradient would still describe the all-electron defaults.
attach_inline_ecp_options_from_basis_sidecar(
scf_options, mol, BasisSet(mol, basis_name)
)
# ECP centres are absolute coordinates on the SCF options. Refuse a
# set that does not sit on the atoms of ``mol`` before spending any
# SCF on it: such an SCF is already inconsistent (bare Z on the atom,
# V_ECP still acting), and the displaced-geometry loop below relies on
# this mapping to move each centre with its atom (#576).
ecp_centre_atom_indices(scf_options, mol)
# Nuclear charges for the dipole-derivative rows: Z - n_core on ECP
# atoms (both routes), bare Z otherwise. Resolved once here -- the
# displaced geometries carry the same ECP assignment (#642).
_dipole_nuclear_charges = effective_nuclear_charges(mol, scf_options)
if source is None and method_norm in {"RKS", "UKS"}:
from .gradient_terms import functional_gradient_terms_missing
missing = functional_gradient_terms_missing(
scf_options.functional, spin=1 if method_norm == "RKS" else 2,
)
if missing:
from .output import write
write(
f" Note: FD Hessian for {method_norm} functional "
f"{scf_options.functional!r} uses full-energy central-difference "
f"gradients (inner step {_KS_ENERGY_GRADIENT_STEP_BOHR:g} bohr); "
f"the analytic gradient omits {'; '.join(missing)}.\n"
)
# ---- 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 not np.isfinite(delta) or delta <= 0:
raise ValueError("FD Hessian: step_bohr must be finite and positive")
# Partial-Hessian (frozen atoms): displace only the unfrozen atoms;
# the result is the (3M, 3M) sub-block over them (see HessianFDOptions
# .frozen_indices). Frozen columns are never filled and frozen rows
# are dropped at the end.
frozen_set = {int(i) for i in (hessian_options.frozen_indices or ())}
bad = sorted(i for i in frozen_set if i < 0 or i >= n_atoms)
if bad:
raise ValueError(
f"FD Hessian: frozen_indices {bad} out of range [0, {n_atoms})"
)
is_partial = len(frozen_set) > 0
unfrozen = [i for i in range(n_atoms) if i not in frozen_set]
if not unfrozen:
raise ValueError(
"FD Hessian: every atom is frozen -- nothing to displace"
)
unfrozen_dof = [3 * i + c for i in unfrozen for c in range(3)]
# 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
omx_parameter_snapshot = None
for i in range(n_atoms):
if i in frozen_set:
continue # frozen atom: not displaced (partial Hessian)
for cart in range(3):
mol_plus = _displaced_molecule(mol, i, cart, +delta)
mol_minus = _displaced_molecule(mol, i, cart, -delta)
if is_omx:
basis_p, res_p, g_plus = _omx_and_gradient_at(mol_plus, method_norm)
basis_m, res_m, g_minus = _omx_and_gradient_at(mol_minus, method_norm)
for displaced_result in (res_p, res_m):
snapshot = (
displaced_result.parameter_identity,
displaced_result.parameter_sha256,
)
if not all(snapshot):
raise RuntimeError("FD Hessian: missing OMx parameter identity.")
if omx_parameter_snapshot is None:
omx_parameter_snapshot = snapshot
elif snapshot != omx_parameter_snapshot:
raise RuntimeError(
"FD Hessian: OMx parameters changed across displacements."
)
elif source is not None:
# The source evaluates its own SCF and analytic gradient at
# each displaced geometry (ECP centres and dispersion included).
_, g_plus = source(mol_plus)
basis_p, res_p = source.last_basis, source.last_result
_, g_minus = source(mol_minus)
basis_m, res_m = source.last_basis, source.last_result
g_plus = np.asarray(g_plus, dtype=float).reshape(n_atoms, 3)
g_minus = np.asarray(g_minus, dtype=float).reshape(n_atoms, 3)
else:
_solvent_kw = {} if solvent_model is None else {"solvent": solvent_model}
basis_p, res_p, g_plus = _scf_and_gradient_at(
mol_plus, basis_name, method_norm,
scf_options, grid_options, gradient_options,
reference_mol=mol, **_solvent_kw,
)
basis_m, res_m, g_minus = _scf_and_gradient_at(
mol_minus, basis_name, method_norm,
scf_options, grid_options, gradient_options,
reference_mol=mol, **_solvent_kw,
)
if dispersion_params is not None:
# The dispersion gradient at each displaced geometry, so the
# Hessian differentiates the complete gradient (#548).
from .dispersion import compute_dispersion
g_plus = g_plus + np.asarray(compute_dispersion(
mol_plus, dispersion_params, with_gradient=True,
).gradient, dtype=float).reshape(np.shape(g_plus))
g_minus = g_minus + np.asarray(compute_dispersion(
mol_minus, dispersion_params, with_gradient=True,
).gradient, dtype=float).reshape(np.shape(g_minus))
n_disp += 2
# Larsen et al., J. Phys.: Condens. Matter 29, 273002 (2017),
# sec. 8.1, eq. (3): H = -dF/dR = dg/dR (g is not force).
# Column 3i+cart of H = (g(+d) - g(-d)) / (2d), 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, _dipole_nuclear_charges
)
mu_m = _dipole_at(
mol_minus, basis_m, res_m, dipole_origin, _dipole_nuclear_charges
)
dipole_derivs[3 * i + cart, :] = (mu_p - mu_m) / (2.0 * delta)
# Partial Hessian: keep only the unfrozen x unfrozen sub-block (frozen
# columns were never filled; the frozen rows are dropped here too).
if is_partial:
H = H[np.ix_(unfrozen_dof, unfrozen_dof)]
if include_dipole:
dipole_derivs = dipole_derivs[unfrozen_dof, :]
if hessian_options.sym_strict:
H = 0.5 * (H + H.T)
# ---- Mass weight ------------------------------------------------------
# For a partial Hessian, mass-weight with the unfrozen masses only so
# the (3M, 3M) block matches its 3M rows/columns.
masses_amu_all = _resolve_masses(mol, hessian_options.atomic_masses_amu)
masses_amu = masses_amu_all[unfrozen] if is_partial else masses_amu_all
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)
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 -------------------------------
# A partial Hessian is anchored by its frozen atoms -- it has no
# translational/rotational zero modes, so projection is skipped and all
# 3M modes are vibrational.
is_linear = False
if hessian_options.project_trans_rot and not is_partial:
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,
unfrozen_indices=(np.asarray(unfrozen, dtype=int)
if is_partial else None),
solvent=(solvent_label(solvent_model) if solvent_model is not None else None),
)
def _warn_if_frame_tie(model, mol, step_bohr: float) -> bool:
"""Warn when a solvated Hessian would straddle the COSMO FINE cavity's frame tie.
In the molecule-fixed frame a symmetric solute's frame picks one of several
symmetry-equivalent atoms, so the energy has a cusp at the symmetric
geometry (#535). A Hessian differences gradients at displaced geometries,
which sit on opposite branches of that cusp whenever the solute is symmetric
to within the step, and the antisymmetric modes inherit the jump. Measured on
symmetric water, RHF/STO-3G, fine cavity at 0.35 A: the asymmetric stretch is
4381 cm^-1 against 4426 on the default Lebedev cavity and 4424 on the fine
cavity in the lab frame, an error of 45 cm^-1 that nothing else announces.
The gradient projection of #535 cannot help: the displaced geometries are far
outside its tolerance.
Returns whether it warned.
"""
import warnings
if getattr(model, "cavity", "lebedev") != "fine" or getattr(model, "fine_frame", "") != "molecular":
return False
from .solvation.point_group import symmetry_operations
atoms = list(mol.atoms)
ops = symmetry_operations(
[int(a.Z) for a in atoms],
np.array([list(a.xyz) for a in atoms], dtype=np.float64),
max(2.0 * float(step_bohr), 1e-3),
)
# An operation that moves no atom (a planar solute's own plane) cannot tie the
# frame's selections; only one that permutes atoms can.
n_atoms = len(atoms)
moving = [op for op in ops if not np.array_equal(op[1], np.arange(n_atoms))]
if not moving:
return False
warnings.warn(
"FD Hessian in solvent on the COSMO FINE cavity in the molecule-fixed "
"frame, at a geometry with point-group symmetry (%d operations that "
"permute atoms, within the displacement): the frame ties here, the "
"energy has a cusp, and the "
"displaced gradients straddle it, so modes that break the symmetry are "
"wrong (measured: 45 cm^-1 on water's asymmetric stretch). Use "
"SolventModel(cavity='fine', fine_frame='lab') or the default Lebedev "
"cavity for frequencies." % len(moving),
UserWarning,
stacklevel=3,
)
return True
def solvent_label(model) -> str:
"""A short description of a resolved solvent model for output and manifests,
for example ``"CPCM water (eps = 78.39)"``."""
family = {"cpcm": "CPCM", "cosmo": "COSMO", "dcosmo-rs": "Direct COSMO-RS"}.get(
str(model.variant), str(model.variant)
)
eps = float(model.epsilon)
eps_text = "inf" if eps == float("inf") else f"{eps:g}"
cavity = "" if str(getattr(model, "cavity", "lebedev")) == "lebedev" else ", fine cavity"
return f"{family} {model.name} (eps = {eps_text}{cavity})"
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, b] =
# (dmu_b/dR_k) / sqrt(M_k). Then transform to normal-mode basis:
# dmu/dQ[b, p] = S_k B_mw[k, b] . 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 dmu/dQ by sum rule) with rotations (nonzero
# dmu/dQ 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",
]