"""Phase 17b-1 — Coupled-Perturbed Hartree-Fock (CPHF) kernel for RHF.
Solves the linear-response equation
.. math::
\\mathbf{A} \\, \\mathbf{U}^x = -\\mathbf{B}^x
where :math:`\\mathbf{U}^x` is the orbital-rotation amplitude in the
occupied-virtual block (response of the MOs to perturbation *x*),
:math:`\\mathbf{A}` is the orbital Hessian (closed-shell, real,
spin-restricted) and :math:`\\mathbf{B}^x` is the "RHS" describing
how perturbation *x* couples to occupied-virtual orbital pairs.
The orbital Hessian acts as
.. math::
[\\mathbf{A} \\, \\mathbf{v}]_{ai} = (\\varepsilon_a - \\varepsilon_i) v_{ai}
+ \\sum_{bj} \\big[ 4(ai|bj) - (ab|ij) - (aj|ib) \\big] v_{bj}
(the factor 4 / -1 / -1 is the closed-shell-RHF combination of the
exchange + Coulomb response). We never form :math:`\\mathbf{A}`
explicitly. Each matrix-vector product builds a Fock-like contribution
in AO basis (:math:`G[D^v]` with :math:`D^v_{\\mu\\nu} = 2 \\sum_{ai}
v_{ai} (C_{\\mu a} C_{\\nu i} + C_{\\mu i} C_{\\nu a})`, factor 2
because the closed-shell amplitudes describe both spins) and
transforms back to the MO occ-vir block. Cost per iteration: one
:math:`O(N^4)` Fock build, identical to one SCF iteration.
The CG solver uses the diagonal :math:`(\\varepsilon_a -
\\varepsilon_i)^{-1}` as a preconditioner — the orbital Hessian is
strongly diagonal-dominant for non-pathological systems, so a few
tens of iterations is typical.
API
---
::
from vibeqc import Atom, Molecule, BasisSet, RHFOptions, run_rhf
from vibeqc import compute_eri, compute_dipole
from vibeqc.cphf import cphf_solve_rhf, dipole_polarizability_rhf
mol = Molecule.from_xyz("water.xyz")
basis = BasisSet(mol, "6-31g*")
result = run_rhf(mol, basis, RHFOptions(conv_tol_grad=1e-9))
alpha = dipole_polarizability_rhf(result, basis, mol)
print(alpha) # (3, 3) tensor in atomic units (e²·a₀²/E_h)
The first arg of :func:`cphf_solve_rhf` is the AO-basis RHS matrix
(or a list of them). Returns the orbital-response amplitudes ``U``
in occ-vir MO basis.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Optional, Sequence
import numpy as np
from ._vibeqc_core import (
BasisSet,
Molecule,
RHFResult,
compute_dipole,
compute_eri,
)
from .properties import center_of_mass
# ----------------------------------------------------------------------
# Public dataclasses
# ----------------------------------------------------------------------
[docs]
@dataclass
class CPHFOptions:
"""Knobs for :func:`cphf_solve_rhf`.
max_iter
Conjugate-gradient iteration cap. Most cases converge in
20–60 iterations on default basis sets; increase if you see
``CPHFConvergenceError``.
tol
Relative residual tolerance ``||r|| < tol·||rhs||``. Default
1e-8 — tight enough for analytic Hessian work, loose enough
that SCF residual noise doesn't dominate.
use_preconditioner
When True (default), use :math:`(\\varepsilon_a -
\\varepsilon_i)^{-1}` as a CG preconditioner. Disable only
for diagnostics — convergence is much slower without it.
"""
max_iter: int = 100
tol: float = 1e-8
use_preconditioner: bool = True
class CPHFConvergenceError(RuntimeError):
"""Raised when the CG solver fails to converge to ``options.tol``
within ``options.max_iter`` iterations."""
# ----------------------------------------------------------------------
# Internal helpers — orbital-Hessian matrix-vector product
# ----------------------------------------------------------------------
def _ao_density_from_ov(v_ia: np.ndarray, C_occ: np.ndarray,
C_vir: np.ndarray) -> np.ndarray:
"""Form the symmetric AO-basis "response density"
D^v_{μν} = Σ_{ia} v_{ia} (C_μa C_νi + C_μi C_νa)
No extra spin factor: the closed-shell orbital-Hessian
coefficients ``4 (ai|jb) − (ab|ij) − (aj|ib)`` already absorb
the doubly-occupied spin algebra. The factor "2" prefactor on
G(D^v) in the orbital-Hessian action gives the correct
combination "2 · (J − ½K)·D^v = 4 J − K" (in MO occ-vir).
Indexing convention: ``v_ia[i, a]`` (occ is row, vir is column).
"""
# T_{μi} = Σ_a v_{ia} C_μa = (C_vir · v_ia^T)_{μi}
T = C_vir @ v_ia.T # (n_basis, n_occ)
# D = T C_occ^T + C_occ T^T (no leading factor of 2)
D = T @ C_occ.T
return D + D.T
def _ao_density_antisym_from_ov(v_ia: np.ndarray, C_occ: np.ndarray,
C_vir: np.ndarray) -> np.ndarray:
"""Form the *antisymmetric* AO-basis "response density"
D^v_{μν} = Σ_{ia} v_{ia} (C_μa C_νi − C_μi C_νa)
This is the density-like matrix the ``(A − B)`` orbital-Hessian
action couples to. Its Coulomb build J[D] is identically zero
(ERI symmetric in λσ, D antisymmetric), so only the exchange
build contributes — see :func:`_orbital_hessian_minus_action`.
Indexing convention: ``v_ia[i, a]`` (occ row, vir column),
matching :func:`_ao_density_from_ov`.
"""
T = C_vir @ v_ia.T # (n_basis, n_occ)
D = T @ C_occ.T
return D - D.T
def _ov_block(M_ao: np.ndarray, C_occ: np.ndarray,
C_vir: np.ndarray) -> np.ndarray:
"""Transform an AO-basis matrix into the occupied-virtual block of
MO basis. Returns shape ``(n_occ, n_vir)`` in ``[i, a]`` indexing
(occ row, vir column) — matches the convention used everywhere
else in this module.
"""
return C_occ.T @ M_ao @ C_vir
def _build_g_ao(D_ao: np.ndarray, eri: np.ndarray) -> np.ndarray:
"""Closed-shell Fock-like contribution from a (possibly non-RHF)
density-like matrix:
G[D] = J[D] − (1/2) K[D]
J[D]_{μν} = Σ_{λσ} (μν|λσ) D_{λσ}
K[D]_{μν} = Σ_{λσ} (μλ|νσ) D_{λσ}
Implemented as numpy einsums on the stored 4-index ERI tensor.
"""
# J: contract λσ on D against the right two indices of (μν|λσ).
J = np.einsum("uvls,ls->uv", eri, D_ao)
# K: contract μ→λ swap. (μλ|νσ) = eri[μ,λ,ν,σ]; contract λσ on D:
# K_{μν} = Σ_{λσ} (μλ|νσ) D_{λσ}
K = np.einsum("ulvs,ls->uv", eri, D_ao)
return J - 0.5 * K
def _build_k_ao(D_ao: np.ndarray, eri: np.ndarray) -> np.ndarray:
"""Plain exchange build (no Coulomb, no 1/2 factor):
K[D]_{μν} = Σ_{λσ} (μλ|νσ) D_{λσ}
Used by the ``(A − B)`` orbital-Hessian action. For an
*antisymmetric* density-like matrix the Coulomb build J[D]
vanishes identically (the ERI is symmetric under λ↔σ while D is
not), so ``(A − B)``'s two-electron part is pure exchange.
"""
return np.einsum("ulvs,ls->uv", eri, D_ao)
def _orbital_hessian_action(v_ia: np.ndarray, C_occ: np.ndarray,
C_vir: np.ndarray, eps_diff_inv: np.ndarray,
eri: np.ndarray) -> np.ndarray:
"""Apply the closed-shell RHF orbital Hessian A to a vector
v_{ia} in MO occ-vir basis. Returns A·v with the same
``(n_occ, n_vir)`` shape and ``[i, a]`` indexing.
A·v decomposes into:
diag: (ε_a − ε_i) v_{ia}
2-electron: 2·G[D^v]^{ov} in the same MO occ-vir basis.
The "2·" prefactor and "+J − (1/2)K" structure together give the
closed-shell spin-restricted "4 J − K − (exchange-mix)" combination
of the textbook orbital Hessian.
"""
# eps_diff_inv has shape (n_occ, n_vir): 1/(ε_a − ε_i)
# Diagonal piece: (ε_a − ε_i) v_{ia} = v_ia / eps_diff_inv.
diag_part = v_ia / eps_diff_inv
# 2-electron piece via Fock-like AO build.
D_v = _ao_density_from_ov(v_ia, C_occ, C_vir)
G_v = _build_g_ao(D_v, eri)
g_part_ia = 2.0 * _ov_block(G_v, C_occ, C_vir)
return diag_part + g_part_ia
def _orbital_hessian_minus_action(v_ia: np.ndarray, C_occ: np.ndarray,
C_vir: np.ndarray,
eps_diff_inv: np.ndarray,
eri: np.ndarray) -> np.ndarray:
"""Apply the closed-shell RHF ``(A − B)`` operator to a vector
``v_{ia}`` in MO occ-vir basis. Returns ``(A − B)·v`` with the
same ``(n_occ, n_vir)`` shape and ``[i, a]`` indexing.
In TD-HF notation the closed-shell singlet super-matrices are
A_{ia,jb} = (ε_a − ε_i) δ_ij δ_ab + 2(ia|jb) − (ij|ab)
B_{ia,jb} = 2(ia|jb) − (ib|ja)
so
(A − B)_{ia,jb} = (ε_a − ε_i) δ_ij δ_ab − (ij|ab) + (ib|ja).
The two-electron part ``(ib|ja) − (ij|ab)`` is exactly the
occ-vir block of the plain exchange build ``K`` evaluated on the
*antisymmetric* response density (Coulomb cancels) — see
:func:`_ao_density_antisym_from_ov` / :func:`_build_k_ao`. No
prefactor: the doubly-occupied spin algebra of the closed-shell
amplitudes is already absorbed, exactly as for ``(A + B)`` in
:func:`_orbital_hessian_action`.
``(A − B)`` is symmetric positive-definite whenever the
closed-shell RHF reference is stable; it is diagonal for a pure
(HF-exchange-free) functional and picks up the exchange terms
above for HF / hybrids.
"""
diag_part = v_ia / eps_diff_inv
D_anti = _ao_density_antisym_from_ov(v_ia, C_occ, C_vir)
K_anti = _build_k_ao(D_anti, eri)
twoel_part = _ov_block(K_anti, C_occ, C_vir)
return diag_part + twoel_part
# ----------------------------------------------------------------------
# Conjugate-gradient solver (preconditioned)
# ----------------------------------------------------------------------
def _pcg(matvec, rhs, precond, tol: float, max_iter: int):
"""Single-RHS preconditioned conjugate gradient. Solves
``matvec(x) = rhs`` with the symmetric matvec implied by matvec.
Returns (x, n_iter, converged).
"""
rhs_norm = np.linalg.norm(rhs)
if rhs_norm == 0:
return np.zeros_like(rhs), 0, True
x = np.zeros_like(rhs)
r = rhs.copy() # r0 = b - A·x0 = b
z = precond(r)
p = z.copy()
rz = float(np.vdot(r.ravel(), z.ravel()).real)
for k in range(1, max_iter + 1):
Ap = matvec(p)
denom = float(np.vdot(p.ravel(), Ap.ravel()).real)
if denom <= 0:
# Non-PD orbital Hessian (instability or near-degeneracy).
# In practice this shouldn't happen for converged RHF on
# closed-shell systems; bail out gracefully.
raise CPHFConvergenceError(
f"CPHF: orbital Hessian non-positive at CG iter {k} "
f"(p^T A p = {denom:.3e}). Likely an RHF instability."
)
alpha = rz / denom
x = x + alpha * p
r = r - alpha * Ap
if np.linalg.norm(r) < tol * rhs_norm:
return x, k, True
z = precond(r)
rz_new = float(np.vdot(r.ravel(), z.ravel()).real)
beta = rz_new / rz
p = z + beta * p
rz = rz_new
return x, max_iter, False
# ----------------------------------------------------------------------
# Public driver: cphf_solve_rhf
# ----------------------------------------------------------------------
def cphf_solve_rhf(
rhs_ao,
rhf_result: RHFResult,
eri: np.ndarray,
*,
options: Optional[CPHFOptions] = None,
) -> np.ndarray:
"""Solve the closed-shell RHF CPHF equations for one or more
AO-basis perturbations.
Each perturbation matrix ``M_ao`` (shape ``(n_basis, n_basis)``)
is transformed to the occupied-virtual block in MO basis to form
the right-hand side ``B_{ai} = (C_vir^T · M_ao · C_occ)_{ai}``,
then the CPHF equations
:math:`\\mathbf{A}\\,\\mathbf{U} = -\\mathbf{B}` are solved by
preconditioned CG.
Parameters
----------
rhs_ao :
Either a single ``(n_basis, n_basis)`` matrix or a sequence
of them (e.g. shape ``(n_rhs, n_basis, n_basis)``). For
polarizability the three components are the dipole-integral
matrices in the three Cartesian directions.
rhf_result :
Converged :class:`vibeqc.RHFResult` providing MO coefficients
and energies.
eri :
AO-basis 4-index ERI tensor (e.g. from
:func:`vibeqc.compute_eri`).
options :
:class:`CPHFOptions` for solver knobs.
Returns
-------
np.ndarray
Orbital-response amplitudes ``U[n_rhs, n_occ, n_vir]`` (or
``U[n_occ, n_vir]`` for a single RHS).
Raises
------
CPHFConvergenceError
If the CG solver fails to converge to ``options.tol`` in
``options.max_iter`` iterations.
"""
if options is None:
options = CPHFOptions()
# --- MO data ---------------------------------------------------------
C = np.asarray(rhf_result.mo_coeffs, dtype=np.float64)
eps = np.asarray(rhf_result.mo_energies, dtype=np.float64)
if not rhf_result.converged:
raise ValueError(
"cphf_solve_rhf: RHFResult is not converged. CPHF requires "
"a stationary HF reference."
)
n_basis = C.shape[0]
# n_occ = (electron count) / 2. Read from the density matrix's trace
# against the overlap, but RHFResult doesn't expose n_occ directly —
# infer from the column shape and the eps-vs-occ structure. The
# cleanest invariant is "number of orbitals with negative occupation
# eigenvalue", but vibe-qc closed-shell results have no fractional
# occupations, so use a plain sign-of-density approach — if the
# density is 2·C_occ·C_occ^T then trace(D·S) = 2·n_occ. We pull n_occ
# from the result.density rather than re-computing it.
D = np.asarray(rhf_result.density, dtype=np.float64)
# trace(D · S) = 2 n_occ. Easier: use the canonical-MO orthonormality —
# for closed-shell RHF with frozen-HF density, the rank of D/2 is n_occ.
# Eigenvalues of D are {2, 2, ..., 2, 0, 0, ..., 0} so n_occ = round(sum(D_diag)/2).
# Actually that's trace(D)/2 only if S is identity — which it is in
# an orthonormal basis but not in AO. Punt: take MO eigenvalues as the
# signal — those below the HOMO-LUMO gap. The cheapest approach is to
# demand the caller knows n_occ from molecule.n_electrons // 2; we do
# that by inspecting how the density was built. Cleaner: we ask the
# user to pass it. For now, re-derive from the eigenvalue spectrum of
# the density matrix in AO basis.
n_occ = _infer_n_occ(D, C)
C_occ = C[:, :n_occ]
C_vir = C[:, n_occ:]
n_vir = C_vir.shape[1]
eps_occ = eps[:n_occ]
eps_vir = eps[n_occ:]
# eps_diff has shape (n_occ, n_vir): ε_a − ε_i with a=virtual, i=occ
eps_diff = eps_vir[None, :] - eps_occ[:, None]
if np.any(eps_diff <= 0):
raise ValueError(
"cphf_solve_rhf: orbital-energy gap non-positive — the RHF "
"reference is unstable (HOMO above LUMO)."
)
eps_diff_inv = 1.0 / eps_diff
eri_arr = np.asarray(eri, dtype=np.float64)
def _matvec(v_ov):
return _orbital_hessian_action(v_ov, C_occ, C_vir, eps_diff_inv,
eri_arr)
if options.use_preconditioner:
def _precond(r):
return r * eps_diff_inv
else:
def _precond(r):
return r
# --- Process inputs --------------------------------------------------
rhs_arr = np.asarray(rhs_ao, dtype=np.float64)
single = rhs_arr.ndim == 2
if single:
rhs_list = [rhs_arr]
else:
if rhs_arr.ndim != 3:
raise ValueError(
f"cphf_solve_rhf: rhs_ao must be 2D ({n_basis}, {n_basis})"
f" or 3D (n_rhs, {n_basis}, {n_basis}); got shape {rhs_arr.shape}"
)
rhs_list = [rhs_arr[k] for k in range(rhs_arr.shape[0])]
# --- Solve per RHS ---------------------------------------------------
U_list = []
for k, M_ao in enumerate(rhs_list):
if M_ao.shape != (n_basis, n_basis):
raise ValueError(
f"cphf_solve_rhf: RHS {k} has shape {M_ao.shape}, "
f"expected ({n_basis}, {n_basis})"
)
# MO occ-vir block of the RHS, sign flipped: A·U = −B.
# Shape (n_occ, n_vir), [i, a] indexing.
B = -_ov_block(M_ao, C_occ, C_vir)
x, n_iter, converged = _pcg(_matvec, B, _precond,
options.tol, options.max_iter)
if not converged:
raise CPHFConvergenceError(
f"cphf_solve_rhf: CG did not converge for RHS {k} "
f"in {options.max_iter} iterations "
f"(residual ||r|| / ||b|| = {np.linalg.norm(_matvec(x) - B) / np.linalg.norm(B):.3e})."
)
U_list.append(x)
if single:
return U_list[0]
return np.stack(U_list, axis=0)
# ----------------------------------------------------------------------
# MO setup shared by the static and dynamic solvers
# ----------------------------------------------------------------------
def _cphf_mo_setup(rhf_result: RHFResult):
"""Return ``(C_occ, C_vir, eps_diff, eps_diff_inv, n_occ, n_vir)``
for a converged closed-shell RHF reference. Raises ValueError on a
non-stationary or HOMO-above-LUMO reference — both solvers need a
stable closed-shell RHF.
"""
if not rhf_result.converged:
raise ValueError(
"cphf: RHFResult is not converged — CPHF / dynamic response "
"require a stationary HF reference."
)
C = np.asarray(rhf_result.mo_coeffs, dtype=np.float64)
eps = np.asarray(rhf_result.mo_energies, dtype=np.float64)
D = np.asarray(rhf_result.density, dtype=np.float64)
n_occ = _infer_n_occ(D, C)
C_occ = C[:, :n_occ]
C_vir = C[:, n_occ:]
n_vir = C_vir.shape[1]
eps_diff = eps[n_occ:][None, :] - eps[:n_occ][:, None] # (n_occ, n_vir)
if np.any(eps_diff <= 0):
raise ValueError(
"cphf: orbital-energy gap non-positive — the RHF reference "
"is unstable (HOMO above LUMO)."
)
return C_occ, C_vir, eps_diff, 1.0 / eps_diff, n_occ, n_vir
# ----------------------------------------------------------------------
# Dynamic (imaginary-frequency) CPHF / TD-HF linear response
# ----------------------------------------------------------------------
def cphf_solve_dynamic_rhf(
rhs_ao,
rhf_result: RHFResult,
eri: np.ndarray,
omega: float,
*,
options: Optional[CPHFOptions] = None,
) -> np.ndarray:
"""Solve the **imaginary-frequency** closed-shell linear-response
(TD-HF) equations for one or more AO-basis perturbations.
The frequency-dependent TD-HF response to a real, symmetric
one-electron perturbation ``P`` is, with ``U₊ = X + Y`` the
symmetric response amplitude,
.. math::
\\big[ (A + B) + \\omega^2 (A - B)^{-1} \\big] \\, U_+ = -2 P
evaluated on the **imaginary** frequency axis ``iω`` (``omega`` is
the real magnitude ``ω ≥ 0``). On the imaginary axis the shifted
operator is real and symmetric positive-definite — there are no
poles — which is exactly why D3/D4 reference dispersion data is
built there.
At ``omega == 0`` this reduces to ``(A + B) U_+ = -2 P``; note the
factor-2 RHS relative to :func:`cphf_solve_rhf` (which solves
``(A + B) U = -P``), so ``U_+ == 2 · U_static`` in the static
limit. :func:`dynamic_polarizability_rhf` absorbs the factor in
its contraction prefactor so ``α(i0)`` equals the static
:func:`dipole_polarizability_rhf` tensor.
Solver: outer preconditioned CG on the ``ω²``-shifted operator;
each application of ``(A - B)^{-1}`` is an inner preconditioned CG
solve against ``(A - B)`` (positive-definite for a stable
closed-shell RHF). For ``omega == 0`` the inner solve is skipped.
Parameters
----------
rhs_ao :
A single ``(n_basis, n_basis)`` perturbation matrix or a
sequence / 3D stack of them.
rhf_result :
Converged closed-shell :class:`RHFResult`.
eri :
AO-basis 4-index ERI tensor.
omega :
Imaginary-frequency magnitude ``ω ≥ 0`` (atomic units). The
physical frequency is ``iω``.
options :
:class:`CPHFOptions` solver knobs (shared by the inner and
outer CG).
Returns
-------
np.ndarray
Symmetric response amplitudes ``U_+`` with shape
``(n_rhs, n_occ, n_vir)`` (or ``(n_occ, n_vir)`` for a single
RHS).
"""
if options is None:
options = CPHFOptions()
omega = float(omega)
if omega < 0:
raise ValueError(
f"cphf_solve_dynamic_rhf: omega must be >= 0 (imaginary-axis "
f"magnitude); got {omega}"
)
C_occ, C_vir, eps_diff, eps_diff_inv, n_occ, n_vir = _cphf_mo_setup(
rhf_result)
n_basis = C_occ.shape[0]
eri_arr = np.asarray(eri, dtype=np.float64)
def _matvec_plus(v):
return _orbital_hessian_action(v, C_occ, C_vir, eps_diff_inv, eri_arr)
def _matvec_minus(v):
return _orbital_hessian_minus_action(v, C_occ, C_vir, eps_diff_inv,
eri_arr)
omega_sq = omega * omega
# Preconditioners. The inner (A−B) solve uses the orbital-energy-
# difference diagonal 1/(ε_a−ε_i) — (A−B) is strongly diagonal-
# dominant. The *outer* operator is (A+B) + ω²(A−B)⁻¹; its
# diagonal is ≈ (ε_a−ε_i) + ω²/(ε_a−ε_i), so the matching
# preconditioner is the inverse of that ω-shifted diagonal.
# Using the plain 1/(ε_a−ε_i) outer preconditioner badly
# under-damps the ω²(A−B)⁻¹ part and stalls the outer CG at high
# ω (the integrand tail of the Casimir-Polder grid). The shifted
# preconditioner reduces to the static one at ω=0.
if options.use_preconditioner:
outer_precond_diag = eps_diff + omega_sq * eps_diff_inv
def _precond_inner(r):
return r * eps_diff_inv
def _precond_outer(r):
return r / outer_precond_diag
else:
def _precond_inner(r):
return r
def _precond_outer(r):
return r
# Inner solve tolerance: tighter than the outer so the shifted
# operator is effectively exact and the outer CG doesn't stall on
# an inconsistent matvec.
inner_tol = max(options.tol * 1.0e-2, 1.0e-12)
def _shifted_op(v):
out = _matvec_plus(v)
if omega_sq != 0.0:
w, _, conv = _pcg(_matvec_minus, v, _precond_inner,
inner_tol, options.max_iter)
if not conv:
raise CPHFConvergenceError(
"cphf_solve_dynamic_rhf: inner (A-B) CG did not "
f"converge in {options.max_iter} iterations."
)
out = out + omega_sq * w
return out
rhs_arr = np.asarray(rhs_ao, dtype=np.float64)
single = rhs_arr.ndim == 2
if single:
rhs_list = [rhs_arr]
else:
if rhs_arr.ndim != 3:
raise ValueError(
f"cphf_solve_dynamic_rhf: rhs_ao must be 2D "
f"({n_basis}, {n_basis}) or 3D (n_rhs, {n_basis}, "
f"{n_basis}); got shape {rhs_arr.shape}"
)
rhs_list = [rhs_arr[k] for k in range(rhs_arr.shape[0])]
U_list = []
for k, M_ao in enumerate(rhs_list):
if M_ao.shape != (n_basis, n_basis):
raise ValueError(
f"cphf_solve_dynamic_rhf: RHS {k} has shape "
f"{M_ao.shape}, expected ({n_basis}, {n_basis})"
)
# RHS is −2 P^{ov} (the factor 2 is the U₊ = X+Y convention).
rhs = -2.0 * _ov_block(M_ao, C_occ, C_vir)
x, _, converged = _pcg(_shifted_op, rhs, _precond_outer,
options.tol, options.max_iter)
if not converged:
raise CPHFConvergenceError(
f"cphf_solve_dynamic_rhf: outer CG did not converge for "
f"RHS {k} at omega={omega:.6g} in {options.max_iter} "
f"iterations (||r||/||b|| = "
f"{np.linalg.norm(_shifted_op(x) - rhs) / np.linalg.norm(rhs):.3e})."
)
U_list.append(x)
if single:
return U_list[0]
return np.stack(U_list, axis=0)
# ----------------------------------------------------------------------
# Application: static dipole polarizability
# ----------------------------------------------------------------------
def dipole_polarizability_rhf(
rhf_result: RHFResult,
basis: BasisSet,
molecule: Molecule,
*,
eri: Optional[np.ndarray] = None,
origin: Optional[Sequence[float]] = None,
options: Optional[CPHFOptions] = None,
) -> np.ndarray:
"""Static (frequency-zero) dipole polarizability tensor
:math:`\\alpha_{\\alpha\\beta}` from CPHF.
.. math::
\\alpha_{\\alpha\\beta} = -2 \\sum_{ai} U^{\\alpha}_{ai} \\, \\mu_{\\beta, ai}
where :math:`\\mu_{\\alpha, ai}` is the AO→MO occ-vir block of
the dipole integral in direction :math:`\\alpha` and
:math:`U^{\\alpha}_{ai}` solves the CPHF equations with the
α-direction dipole as the RHS. The factor 2 is for closed-shell
spin and the negative sign comes from the electronic charge
convention (electrons have ``Q = -1`` in atomic units).
Returns
-------
np.ndarray
Symmetric ``(3, 3)`` polarizability tensor in atomic units
(``e²·a₀²/E_h`` = bohr³ when expressed via dimensional
analysis of the field-induced dipole). Multiply by 0.14818
to get ų.
"""
if origin is None:
origin_vec = center_of_mass(molecule)
else:
origin_vec = np.asarray(origin, dtype=np.float64)
dip = compute_dipole(basis, [float(x) for x in origin_vec])
M = np.stack([np.asarray(dip.x), np.asarray(dip.y),
np.asarray(dip.z)], axis=0) # (3, n_bf, n_bf)
# The dipole-perturbation Hamiltonian for an electron is (-r), but
# the CPHF "RHS" is conventionally μ itself (without the sign);
# the sign convention is reabsorbed into the polarizability formula
# below.
if eri is None:
eri = np.asarray(compute_eri(basis), dtype=np.float64)
U = cphf_solve_rhf(M, rhf_result, eri, options=options) # (3, n_occ, n_vir)
# Build dipole occ-vir block in MO basis for the contraction.
# Shape (n_occ, n_vir), [i, a] indexing.
C = np.asarray(rhf_result.mo_coeffs, dtype=np.float64)
n_occ = _infer_n_occ(np.asarray(rhf_result.density, dtype=np.float64), C)
C_occ = C[:, :n_occ]
C_vir = C[:, n_occ:]
mu_mo = np.zeros((3, n_occ, C_vir.shape[1]))
for axis in range(3):
mu_mo[axis] = C_occ.T @ M[axis] @ C_vir
# α_αβ = −4 Σ_ia U^β_ia · μ_α,ia (closed-shell convention).
# The factor 4 = 2 (closed-shell spin in P = 2 C_occ C_occ^T)
# × 2 (each orbital contributes via both ⟨∂ψ|r|ψ⟩ and ⟨ψ|r|∂ψ⟩
# in the dipole expectation, for real orbitals).
alpha = -4.0 * np.einsum("xia,yia->xy", U, mu_mo)
# Symmetrize (the response is symmetric for a static field; tiny
# residual asymmetry is FD/CG noise).
return 0.5 * (alpha + alpha.T)
# ----------------------------------------------------------------------
# Application: imaginary-frequency dipole polarizability α(iω)
# ----------------------------------------------------------------------
def dynamic_polarizability_rhf(
rhf_result: RHFResult,
basis: BasisSet,
molecule: Molecule,
frequencies,
*,
eri: Optional[np.ndarray] = None,
origin: Optional[Sequence[float]] = None,
options: Optional[CPHFOptions] = None,
) -> np.ndarray:
"""Imaginary-frequency dipole polarizability :math:`\\alpha(i\\omega)`.
Evaluates the dynamic dipole polarizability tensor on the
imaginary frequency axis,
.. math::
\\alpha_{\\alpha\\beta}(i\\omega)
= -2 \\sum_{ia} U_+^{\\beta}(i\\omega)_{ia}\\;\\mu_{\\alpha,ia}
where :math:`U_+^{\\beta}(i\\omega)` solves the imaginary-frequency
TD-HF response equations (see :func:`cphf_solve_dynamic_rhf`).
This is the quantity the D4 dispersion model needs: the reference
:math:`C_6` coefficients come from the Casimir-Polder integral
:math:`C_6 = \\tfrac{3}{\\pi}\\int_0^\\infty \\alpha_A(i\\omega)\\,
\\alpha_B(i\\omega)\\,d\\omega`, and :math:`\\alpha(i\\omega)` is
smooth, real and monotone-decaying there (no poles) — milestone
D4b-1 of ``docs/handover_d4_native.md``.
Parameters
----------
rhf_result, basis, molecule :
Converged closed-shell RHF reference and its system.
frequencies :
A scalar imaginary-frequency magnitude :math:`\\omega \\ge 0`,
or a 1-D sequence of them (atomic units).
eri :
Optional precomputed AO 4-index ERI tensor (built once and
reused across all frequencies when omitted).
origin :
Gauge origin for the dipole integrals; defaults to the centre
of mass (the polarizability is origin-independent for a
neutral closed-shell molecule, so this only matters for ions).
options :
:class:`CPHFOptions` solver knobs.
Returns
-------
np.ndarray
``(3, 3)`` tensor for a scalar ``frequencies``; otherwise an
``(n_freq, 3, 3)`` stack, one tensor per frequency. Atomic
units (bohr³). :math:`\\alpha(i0)` reproduces the static
:func:`dipole_polarizability_rhf` tensor.
"""
scalar_input = np.isscalar(frequencies)
freqs = np.atleast_1d(np.asarray(frequencies, dtype=np.float64))
if np.any(freqs < 0):
raise ValueError(
"dynamic_polarizability_rhf: imaginary-frequency magnitudes "
"must be >= 0."
)
if origin is None:
origin_vec = center_of_mass(molecule)
else:
origin_vec = np.asarray(origin, dtype=np.float64)
dip = compute_dipole(basis, [float(x) for x in origin_vec])
M = np.stack([np.asarray(dip.x), np.asarray(dip.y),
np.asarray(dip.z)], axis=0) # (3, n_bf, n_bf)
if eri is None:
eri = np.asarray(compute_eri(basis), dtype=np.float64)
# Dipole occ-vir MO block, shared across all frequencies.
C = np.asarray(rhf_result.mo_coeffs, dtype=np.float64)
n_occ = _infer_n_occ(np.asarray(rhf_result.density, dtype=np.float64), C)
C_occ = C[:, :n_occ]
C_vir = C[:, n_occ:]
mu_mo = np.zeros((3, n_occ, C_vir.shape[1]))
for axis in range(3):
mu_mo[axis] = C_occ.T @ M[axis] @ C_vir
out = np.zeros((freqs.size, 3, 3))
for f, omega in enumerate(freqs):
U_plus = cphf_solve_dynamic_rhf(M, rhf_result, eri, omega,
options=options) # (3, n_occ, n_vir)
# α_αβ(iω) = −2 Σ U₊^β · μ_α. The U₊ = X+Y convention carries
# the factor 2 relative to the static U; combined with the −2
# here it reproduces the static −4 prefactor at ω = 0.
alpha = -2.0 * np.einsum("xia,yia->xy", U_plus, mu_mo)
out[f] = 0.5 * (alpha + alpha.T)
return out[0] if scalar_input else out
# ----------------------------------------------------------------------
# n_occ inference
# ----------------------------------------------------------------------
def _infer_n_occ(D_ao: np.ndarray, C: np.ndarray) -> int:
"""Recover the number of doubly-occupied orbitals from the closed-
shell density matrix and the MO coefficients.
The MO-basis density matrix is ``C^T · S · D · S · C`` for a
non-orthonormal AO basis, but we only have the AO Fock and
density at hand. A cheaper trick: in the *MO basis*, the closed-
shell RHF density matrix is ``diag(2, 2, ..., 2, 0, ..., 0)`` so
looking at ``(C^T D)`` and counting columns with non-trivial
inner product with C_occ would work, but again needs S.
The simplest robust signal: the density's trace ``Tr(D)`` equals
twice the integrated density (number of electrons) only when AOs
are orthonormal. In general we'd need ``Tr(D · S)`` = N.
For a converged RHF solution though, ``D = 2 · C_occ · C_occ^T``,
so ``D · C = 2 · C_occ · (C_occ^T · C)`` and this projects onto
the occupied subspace. ``np.linalg.matrix_rank`` of ``D`` gives
n_occ since D has rank n_occ exactly (the AO ranks aren't reduced
by C being non-square here; D's rank is set by the rank of the
occupied projection).
Use rank-of-density.
"""
# Eigenvalues of D are {2, 2, ..., 2, 0, ..., 0} in an orthonormal
# basis; in AO basis they're scaled but the rank is preserved.
eigvals = np.linalg.eigvalsh(D_ao)
# Most singular values are 0; n_occ are nonzero positive.
threshold = 1e-8 * eigvals.max() if eigvals.max() > 0 else 1e-10
return int(np.sum(eigvals > threshold))
__all__ = [
"CPHFOptions",
"CPHFConvergenceError",
"cphf_solve_rhf",
"cphf_solve_dynamic_rhf",
"dipole_polarizability_rhf",
"dynamic_polarizability_rhf",
]