"""CPCM (Conductor-like Polarisable Continuum Model) — core math (S1b).
The conductor-screening CPCM equations on a discrete cavity
tessellation are
A q = − f(ε) · V with f(ε) = (ε − 1) / ε
where ``A`` is the (n_pts, n_pts) cavity self-interaction matrix,
``V`` is the gas-phase molecular electrostatic potential evaluated at
the cavity surface points, and ``q`` is the apparent surface charge.
The total solvation energy is
E_solv = (1/2) Σ_i q_i V_i.
This module assembles ``A`` and solves for ``q``; the SCF wiring
lives in :mod:`vibeqc.solvation.driver`.
References
----------
* Klamt-Schüürmann 1993 — original COSMO.
* Cossi-Rega-Scalmani-Barone 2003 — CPCM matrix layout used here.
* Scalmani-Frisch 2010 — diagonal element from continuous-surface-
charge (CSC) formulation.
"""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
# Scalmani-Frisch (2010) diagonal constant.
#
# A_ii = (CPCM_DIAG_ALPHA × √(4π)) / √(w_i)
# = (1.0694 × √(4π)) / √(w_i)
#
# where ``w_i`` is the actual surface area element at the cavity point
# (bohr², so the dimensional units work out: 1 / √(bohr²) = bohr⁻¹).
# The numerical factor 1.0694 is the self-energy of a uniformly charged
# disk relative to the equivalent Gaussian — see Scalmani-Frisch 2010
# eq. (11) and Cossi-Scalmani-Mennucci-Tomasi 2003 eq. (18). The
# baked-in √(4π) makes the diagonal dominant over typical nearest-
# neighbour 1/r off-diagonals on Lebedev-discretised spheres, which
# the bare 1.0694/√w_i form fails for moderately dense tessellations
# (110+ points per sphere).
#
# Matches PySCF (pyscf/solvent/pcm.py::PCM.get_F0) and Q-Chem (CPCM
# default; see Stratmann-Scuseria-Frisch JCP 109 8218).
import math as _math
CPCM_DIAG_ALPHA_BARE = 1.0694
CPCM_DIAG_ALPHA = CPCM_DIAG_ALPHA_BARE * _math.sqrt(4.0 * _math.pi)
[docs]
@dataclass(frozen=True)
class CPCMResult:
"""Outcome of one apparent-surface-charge solve.
Attributes
----------
q : ndarray, shape (n_pts,)
Apparent surface charges, atomic units.
V : ndarray, shape (n_pts,)
Total molecular electrostatic potential (electron +
nuclear) at the cavity points, atomic units.
e_solv : float, Hartree
Solvation energy E_solv = ½ Σ_i q_i V_i.
epsilon : float
Solvent dielectric constant used.
"""
q: np.ndarray
V: np.ndarray
e_solv: float
epsilon: float
@property
def total_charge(self) -> float:
return float(self.q.sum())
def build_A_matrix(cavity_points: np.ndarray,
cavity_weights: np.ndarray) -> np.ndarray:
"""Assemble the CPCM cavity self-interaction matrix A.
The diagonal uses Scalmani-Frisch's disk self-energy
A_ii = α/√w_i (eq. 11 of their 2010 paper). The off-diagonal is
the Coulomb kernel A_ij = 1/|s_i − s_j|.
All inputs in atomic units (bohr for distances, bohr² for
weights). The returned matrix is symmetric, positive-definite —
suitable for ``numpy.linalg.solve`` or a Cholesky.
Memory is O(n_pts²) doubles. For 240 cavity points that's
~440 KB; for 5000 cavity points that's ~190 MB — typical of
"tight Gaussian-default cavity on a 50-heavy-atom molecule".
Use ``solve_apparent_charges`` to factor + solve in one shot.
"""
pts = np.asarray(cavity_points, dtype=np.float64)
w = np.asarray(cavity_weights, dtype=np.float64)
if pts.ndim != 2 or pts.shape[1] != 3:
raise ValueError(
f"cavity_points must be (n_pts, 3); got {pts.shape}"
)
if w.shape != (pts.shape[0],):
raise ValueError(
f"cavity_weights must be (n_pts,); got {w.shape}"
)
# Pairwise Euclidean distances (bohr) — full N×N construction is
# the same cost as the matvec and avoids a second loop later.
diff = pts[:, None, :] - pts[None, :, :]
dist = np.sqrt(np.einsum("ijk,ijk->ij", diff, diff))
# Off-diagonal: 1/r_ij. Diagonal will overwrite shortly.
with np.errstate(divide="ignore"):
A = 1.0 / dist
# Scalmani-Frisch disk self-energy on the diagonal: A_ii = α/√w_i.
# This stays finite for the smallest-weight (most heavily-switched)
# points so A is positive-definite end-to-end.
np.fill_diagonal(A, CPCM_DIAG_ALPHA / np.sqrt(w))
return A
def dielectric_factor(epsilon: float, *, variant: str = "cpcm") -> float:
"""Conductor screening factor f(ε).
``variant = "cpcm"`` (default) uses the Cossi-Scalmani-Barone
factor ``f = (ε − 1)/ε`` — converges smoothly to the conductor
limit at high ε.
``variant = "cosmo"`` uses the Klamt-Schüürmann factor
``f = (ε − 1)/(ε + x)`` with ``x = 0.5`` (the COSMO X parameter;
Klamt 1993). The two differ by O(1/ε) and agree to better than 1%
above ε ≈ 30 — for water and most polar solvents the choice
moves E_solv by ≲ 0.1 kcal/mol.
"""
if epsilon <= 1.0:
raise ValueError(
f"dielectric_factor: epsilon must be > 1; got {epsilon} "
f"(use variant='gas' or set solvent=None to skip CPCM)."
)
v = variant.lower()
if v == "cpcm":
return (epsilon - 1.0) / epsilon
if v == "cosmo":
x = 0.5
return (epsilon - 1.0) / (epsilon + x)
raise ValueError(
f"dielectric_factor: unknown variant {variant!r} "
f"(use 'cpcm' or 'cosmo')."
)
def solve_apparent_charges(
A: np.ndarray,
V_at_cavity: np.ndarray,
*,
epsilon: float,
variant: str = "cpcm",
) -> CPCMResult:
"""Solve A q = −f(ε) V for the apparent surface charges.
Returns a :class:`CPCMResult` carrying ``q``, ``V``, the solvation
energy ``½ q·V``, and the dielectric used. The matrix ``A`` is
factorised by ``numpy.linalg.solve`` — for repeated solves with
the same cavity, prefer caching a Cholesky outside (the cavity is
fixed for the duration of one SCF macro-iteration set).
"""
V = np.asarray(V_at_cavity, dtype=np.float64).reshape(-1)
if V.shape[0] != A.shape[0]:
raise ValueError(
f"solve_apparent_charges: V length {V.shape[0]} does "
f"not match A shape {A.shape}"
)
f = dielectric_factor(epsilon, variant=variant)
rhs = -f * V
q = np.linalg.solve(A, rhs)
e_solv = 0.5 * float(np.dot(q, V))
return CPCMResult(q=q, V=V, e_solv=e_solv, epsilon=float(epsilon))