"""CPCM (Conductor-like Polarisable Continuum Model) -- core math (S1b).
The conductor-screening CPCM equations on a discrete cavity
tessellation are
A q = - f(e) . V with f(e) = (e - 1) / e
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) S_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 x √(4pi)) / √(w_i)
# = (1.0694 x √(4pi)) / √(w_i)
#
# where ``w_i`` is the actual surface area element at the cavity point
# (bohr^2, so the dimensional units work out: 1 / √(bohr^2) = 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 √(4pi) 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)
def _native_cpcm():
try:
from vibeqc import _vibeqc_core as _core
except ImportError:
return None
required = (
"cpcm_build_A_matrix",
"cpcm_build_capped_A_matrix",
"cpcm_dielectric_factor",
"cpcm_solve_apparent_charges",
)
if not all(hasattr(_core, name) for name in required):
return None
return _core
[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 = 1/2 S_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 = a/√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^2 for
weights). The returned matrix is symmetric, positive-definite --
suitable for ``numpy.linalg.solve`` or a Cholesky.
Memory is O(n_pts^2) 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}"
)
native = _native_cpcm()
if native is not None:
return np.asarray(native.cpcm_build_A_matrix(pts, w), dtype=np.float64)
# Pairwise Euclidean distances (bohr) -- full NxN 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 = a/√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(e).
``variant = "cpcm"`` (default) uses the Cossi-Scalmani-Barone
factor ``f = (e - 1)/e`` -- converges smoothly to the conductor
limit at high e.
``variant = "cosmo"`` uses the Klamt-Schüürmann factor
``f = (e - 1)/(e + x)`` with ``x = 0.5`` (the COSMO X parameter;
Klamt 1993). The two differ by O(1/e) and agree to better than 1%
above e ≈ 30 -- for water and most polar solvents the choice
moves E_solv by ≲ 0.1 kcal/mol.
"""
native = _native_cpcm()
if native is not None:
try:
return float(
native.cpcm_dielectric_factor(float(epsilon), variant.lower())
)
except RuntimeError as exc:
raise ValueError(str(exc)) from exc
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(e) V for the apparent surface charges.
Returns a :class:`CPCMResult` carrying ``q``, ``V``, the solvation
energy ``1/2 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}"
)
native = _native_cpcm()
if native is not None:
try:
res = native.cpcm_solve_apparent_charges(
np.asarray(A, dtype=np.float64),
V,
float(epsilon),
variant.lower(),
)
except RuntimeError as exc:
raise ValueError(str(exc)) from exc
return CPCMResult(
q=np.asarray(res.q, dtype=np.float64),
V=np.asarray(res.V, dtype=np.float64),
e_solv=float(res.e_solv),
epsilon=float(res.epsilon),
)
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))