"""CPCM SCF driver — macro-iteration coupling of density and ASC (S1b, S1d).
Wraps the molecular SCF entry points (RHF / UHF / RKS / UKS) with a
self-consistent CPCM apparent-surface-charge loop:
* Macro-iter 0: solve the gas-phase SCF for ``D_0``.
* Macro-iter k ≥ 1:
1. Evaluate the molecular electrostatic potential V_i (from
``D_{k−1}`` + nuclei) at every cavity tessellation point.
2. Solve ``A q = −f(ε) V`` for the surface charge.
3. Build the one-electron operator V_q^{μν} =
Σ_i q_i ⟨μ|1/|r − s_i||ν⟩ and add it to ``H_core``.
4. Re-run the inner SCF with the modified ``H_core``; obtain
``D_k``.
* Outer convergence on |E_solv(k) − E_solv(k−1)| ≤ ``tol_e_solv``.
The macro-iteration is the canonical CPCM coupling pattern (see
PySCF's ``solvent.PCM``, ORCA's ``CPCM`` keyword, and Cossi-
Scalmani-Barone 2003 § II.B). Three to five outer cycles typically
converge ΔE_solv to 1e-6 Hartree for neutral closed-shell solutes;
charged or strongly-polar solutes may want 8–10.
This driver builds ``H_core`` from Python via the bundled molecular
integrals and uses the ``run_*_scf_with_jk`` low-level entry points
so the inner SCF — DIIS / level shift / linear-dependence projection
/ canonical orthogonalisation — runs unchanged. ``CavityTessellation``
is geometry-tied and recomputed once per ``run_cpcm_scf`` call.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Any, Callable, Mapping, Optional, Union
import numpy as np
from .cavity import CavityTessellation, build_cavity
from .cpcm import (
CPCMResult,
build_A_matrix,
solve_apparent_charges,
)
from .presets import SOLVENT_PRESETS, canonicalise_solvent_name
@dataclass
class SolventModel:
"""Top-level user-facing solvation spec.
Parameters
----------
epsilon : float
Static dielectric constant of the solvent (dimensionless,
ε_rel). Must be > 1; ``epsilon ≤ 1`` is treated as a gas-phase
no-op by :func:`run_cpcm_scf`.
name : str, optional
Human-readable name for logging (e.g. ``"water"``). Defaults
to ``"custom (ε = ...)"`` if omitted.
variant : {"cpcm", "cosmo"}, default "cpcm"
Dielectric factor formula. ``"cpcm"`` uses
``f = (ε − 1)/ε``; ``"cosmo"`` uses ``f = (ε − 1)/(ε + 0.5)``.
See :func:`vibeqc.solvation.cpcm.dielectric_factor`.
radii : dict[int, float], optional
Per-element vdW-radii overrides (Å). Falls back to the bundled
Bondi table; see :data:`vibeqc.solvation.cavity.BONDI_RADII_ANG`.
radii_scale : float, default 1.20
Multiplier on each Bondi radius (PCM / GEPOL convention).
solvent_probe_radius_ang : float, default 0.0
Adds a probe-radius offset to each atomic sphere. ``0.0``
builds the scaled-vdW (SES) cavity; set to e.g. 1.385 Å for
the water-probe SAS.
n_points_per_sphere : int, default 302
Lebedev order per atomic sphere. See :func:`build_cavity`.
switching_sigma_bohr : float, default 0.5
Scalmani-Frisch switching width (bohr).
max_macro_iter : int, default 30
Outer CPCM macro-iteration cap. Three to five iters is typical;
the cap protects against pathological non-convergence.
tol_e_solv : float, default 1e-6
Outer convergence threshold on |ΔE_solv| (Hartree).
"""
epsilon: float
name: str = ""
variant: str = "cpcm"
radii: Optional[dict[int, float]] = None
radii_scale: float = 1.20
solvent_probe_radius_ang: float = 0.0
n_points_per_sphere: int = 302
switching_sigma_bohr: float = 0.5
max_macro_iter: int = 30
tol_e_solv: float = 1e-6
def __post_init__(self):
if not self.name:
self.name = f"custom (ε = {self.epsilon:.3f})"
@property
def is_gas_phase(self) -> bool:
return self.epsilon <= 1.0
# Accepted user input types for ``solvent=``.
SolventSpec = Union[None, str, float, int, Mapping[str, Any], SolventModel]
def resolve_solvent(spec: SolventSpec) -> Optional[SolventModel]:
"""Coerce a user-supplied ``solvent`` argument to a SolventModel.
Accepts:
* ``None`` — gas phase (returns ``None``).
* ``"water"`` (or any other key/alias in
:data:`vibeqc.solvation.SOLVENT_PRESETS`) — looks up the
preset dielectric.
* a numeric ε — treated as a custom-dielectric solvent.
* a dict — passed to :class:`SolventModel` as kwargs (must include
``epsilon``).
* a :class:`SolventModel` — returned as-is.
The case ``"vacuum"`` / ``"none"`` / ``"gas"`` resolves to ε = 1
and triggers the gas-phase no-op branch in :func:`run_cpcm_scf`.
"""
if spec is None:
return None
if isinstance(spec, SolventModel):
return spec
if isinstance(spec, (int, float)) and not isinstance(spec, bool):
eps = float(spec)
return SolventModel(epsilon=eps)
if isinstance(spec, str):
canon = canonicalise_solvent_name(spec)
eps = SOLVENT_PRESETS[canon]
return SolventModel(epsilon=eps, name=canon)
if isinstance(spec, Mapping):
params = dict(spec)
if "epsilon" not in params:
raise ValueError(
"solvent dict must contain 'epsilon' "
"(see vibeqc.solvation.SolventModel for other keys)."
)
return SolventModel(**params)
raise TypeError(
f"Unrecognised solvent spec: {type(spec).__name__}. "
f"Expected None, str, float, dict, or SolventModel."
)
def _build_point_charge_operators(
basis, # vibeqc.BasisSet
cavity_points: np.ndarray,
) -> np.ndarray:
"""Stack of one-electron operators ⟨μ|1/|r − s_i||ν⟩ per cavity point.
Built once per geometry by calling :func:`vibeqc.compute_nuclear`
on a single-atom fake molecule (Z = 1) at each cavity point. The
libint nuclear-attraction kernel returns
``M_i = −⟨μ|1/r_i|ν⟩``; we negate to expose the positive Coulomb
kernel so callers can write ``F += q[i] * stack[i]`` directly
(matches the standard CPCM Fock-contribution sign convention —
derivation in module docstring of :mod:`vibeqc.solvation.cpcm`).
Shape (n_pts, n_bf, n_bf), dtype float64. Memory is
n_pts · n_bf² doubles — ~5 MB for 240 points × 50 basis fns.
Performance: n_pts independent libint calls; each is microseconds
for a single-atom Z=1 nuclear operator on a typical AO basis.
Total setup cost is comparable to a single full Fock build.
"""
# Local import to avoid a top-level circular import via
# ``vibeqc/__init__.py`` (which imports vibeqc.solvation at
# bottom-of-file scope).
from vibeqc import Atom, Molecule, compute_nuclear
n_pts = int(cavity_points.shape[0])
n_bf = int(basis.nbasis)
stack = np.empty((n_pts, n_bf, n_bf), dtype=np.float64)
for i in range(n_pts):
# Fake single-atom molecule with Z = 1 at the cavity point.
# libint returns V_μν = −Z · ⟨μ|1/|r − s_i||ν⟩ = −⟨μ|1/r_i|ν⟩.
# Negate to surface the positive Coulomb kernel.
# charge=1 + Z=1 → 0 electrons, mult=1 → valid empty system.
# The integral builder only reads Z + position from the atom,
# so this gives M_i = −⟨μ|1/|r − s_i||ν⟩ at unit charge without
# tripping Molecule's electron/multiplicity consistency check.
fake = Molecule([Atom(1, tuple(cavity_points[i]))], 1, 1)
M_i = compute_nuclear(basis, fake)
stack[i] = -np.asarray(M_i)
return stack
def _nuclear_potential_at_cavity(
atom_positions_bohr: np.ndarray,
atomic_numbers: np.ndarray,
cavity_points: np.ndarray,
) -> np.ndarray:
"""V_nuc(s_i) = Σ_A Z_A / |R_A − s_i| at every cavity point.
Vectorised; O(n_atoms · n_pts) flops. Closed-form — no integral
machinery needed since the nuclear charge density is a sum of
delta functions.
"""
R = np.asarray(atom_positions_bohr, dtype=np.float64)
Z = np.asarray(atomic_numbers, dtype=np.float64)
S = np.asarray(cavity_points, dtype=np.float64)
# Pairwise distances (n_atoms, n_pts):
dist = np.sqrt(((R[:, None, :] - S[None, :, :]) ** 2).sum(axis=2))
if np.any(dist == 0.0):
raise ValueError(
"Cavity point lies exactly on a nucleus — this should not "
"happen with the GEPOL switching, but did. Check "
"atom_positions_bohr."
)
return (Z[:, None] / dist).sum(axis=0)
def _density_potential_at_cavity(
density: np.ndarray,
point_charge_operators: np.ndarray,
) -> np.ndarray:
"""V_elec(s_i) = −∫ ρ/|r − s_i| dr = +tr(D · M_i) for each i.
Sign convention: ``point_charge_operators[i]`` is the positive
Coulomb kernel ⟨μ|1/r_i|ν⟩. The electron density carries an
implicit negative charge, so
V_elec(s_i) = −Σ_μν D_μν · ⟨μ|1/r_i|ν⟩
= −tensordot(D, M_i, axes=2).
Returns shape (n_pts,) values typically negative for normal
electron distributions (negative density → negative potential).
"""
return -np.einsum("ab,iab->i", density, point_charge_operators)
def _fock_solvent_contribution(
apparent_charges: np.ndarray,
point_charge_operators: np.ndarray,
) -> np.ndarray:
"""V_q^{μν} = − Σ_i q_i · ⟨μ|1/r_i|ν⟩ — the Fock matrix correction.
Derivation: the one-electron operator an electron (charge −1) at
position **r** feels from external point charges ``q_i`` at
cavity positions ``s_i`` is
v(**r**) = Σ_i (−1) · q_i / |**r** − **s_i**|.
In AO basis ``V_q^{μν} = −Σ_i q_i ⟨μ|1/|r − s_i||ν⟩``, which is
what gets *added* to ``H_core`` when the SCF runs in solvent.
``point_charge_operators[i]`` already exposes the positive
Coulomb kernel ``+⟨μ|1/|r−s_i||ν⟩`` (see
:func:`_build_point_charge_operators`), so the minus sign appears
here.
Equivalent identity: at convergence, ``Tr(D V_q) = +q·V_elec``
(with the sign convention ``V_elec_i = −Σ_{μν} D_{μν} ⟨μ|...|ν⟩``
used in :func:`_density_potential_at_cavity`). That identity is
what the total-energy decomposition in :func:`run_cpcm_scf` relies
on — getting this sign wrong tips the in-solvent SCF energy by
``≈ −2 q·V_elec`` (tens of mHa on polar solutes), so the test
suite asserts ``sol.energy < sol.e_gas`` as a hard guard.
"""
return -np.einsum("i,iab->ab", apparent_charges, point_charge_operators)
[docs]
@dataclass
class SolventResult:
"""Outcome of a CPCM-coupled SCF run.
Carries the original SCF result type (RHFResult / UHFResult /
RKSResult / UKSResult) on ``.scf`` plus the converged solvation
diagnostics.
"""
scf: Any # underlying SCF result
energy: float # total energy in solvent
e_solv: float # solvation energy
e_gas: Optional[float] # gas-phase reference
epsilon: float
solvent_name: str
cavity: CavityTessellation
cpcm: CPCMResult
n_macro_iter: int
macro_history: list[dict[str, float]] = field(default_factory=list)
converged: bool = False
def __repr__(self) -> str:
return (
f"SolventResult(solvent={self.solvent_name!r}, "
f"ε={self.epsilon:.3f}, E_solv={self.e_solv:+.6f} Ha, "
f"E_total={self.energy:.6f} Ha, "
f"macro_iter={self.n_macro_iter}, "
f"converged={self.converged})"
)
def _resolve_method_runners(
method: str,
):
"""Map the high-level method to (gas_runner, jk_inner_runner, needs_grid).
Each runner pair returns ``(SCFResult, density)`` tuples so the
macro-iteration can compose them generically.
"""
from vibeqc import (
run_rhf, run_rks, run_uhf, run_uks,
run_rhf_scf_with_jk, run_rks_scf_with_jk,
run_uhf_scf_with_jk, run_uks_scf_with_jk,
)
m = method.lower()
if m == "rhf":
return ("rhf", run_rhf, run_rhf_scf_with_jk, False)
if m == "uhf":
return ("uhf", run_uhf, run_uhf_scf_with_jk, False)
if m == "rks":
return ("rks", run_rks, run_rks_scf_with_jk, True)
if m == "uks":
return ("uks", run_uks, run_uks_scf_with_jk, True)
raise ValueError(
f"run_cpcm_scf: unsupported method {method!r} "
f"(expected one of: rhf, uhf, rks, uks)."
)
def _density_of(result, *, open_shell: bool) -> np.ndarray:
"""Extract the SCF density (closed-shell or α+β total)."""
if open_shell:
return np.asarray(result.density_alpha) + np.asarray(result.density_beta)
return np.asarray(result.density)
def _ensure_cpcm_aux_basis(options, basis) -> None:
"""If ``options`` asks for density fitting / RIJCOSX but leaves
``aux_basis`` empty, autodetect it from the orbital basis name and
fill it in (in place).
The high-level ``run_rhf`` / ``run_rks`` / ``run_uhf`` / ``run_uks``
drivers — used for the CPCM gas-phase bootstrap SCF — refuse
``density_fit=True`` with an empty ``aux_basis``. Resolving the aux
here keeps the gas SCF, the macro-iteration inner SCFs, and the
JKBuilder all on one consistent auxiliary basis. In-place mutation
matches the established options-handling pattern in
``vibeqc.runner._run_single_point`` (which sets ``opts.functional``
the same way).
"""
density_fit = bool(getattr(options, "density_fit", False))
cosx = bool(getattr(options, "cosx", False))
if not density_fit and not cosx:
return
if (getattr(options, "aux_basis", "") or "").strip():
return
from vibeqc import default_aux_basis_for
orbital_name = (getattr(basis, "name", "") or "").strip()
if not orbital_name:
raise ValueError(
"run_cpcm_scf: density_fit/cosx requested but aux_basis is "
"empty and the orbital BasisSet has no name to autodetect "
"from. Set options.aux_basis explicitly (e.g. "
"'def2-svp-jk')."
)
options.aux_basis = default_aux_basis_for(orbital_name, kind="jk")
def _resolve_cpcm_jk_builder(basis, molecule, options):
"""Pick the JKBuilder the CPCM macro-iteration drives the inner SCF
with, off the method options.
* ``options.cosx`` (with or without ``density_fit``) → RIJCOSX:
RI-J Coulomb + chain-of-spheres exchange
(:func:`vibeqc.make_cosx_jk_builder`).
* ``options.density_fit`` → RI-JK
(:func:`vibeqc.make_df_jk_builder`).
* neither → direct, integral-driven
(:func:`vibeqc.make_direct_jk_builder`) — the v0.9.0 default.
DF / RIJCOSX both need an auxiliary basis. ``options.aux_basis``
is used when set; otherwise it is autodetected from the orbital
basis name via :func:`vibeqc.default_aux_basis_for`. The builder
is constructed once and reused across every CPCM macro-iteration
(the basis doesn't move).
"""
from vibeqc import (
BasisSet,
build_grid,
default_aux_basis_for,
make_cosx_jk_builder,
make_df_jk_builder,
make_direct_jk_builder,
)
density_fit = bool(getattr(options, "density_fit", False))
cosx = bool(getattr(options, "cosx", False))
if not density_fit and not cosx:
return make_direct_jk_builder(basis)
# DF / RIJCOSX path — resolve the auxiliary basis.
aux_name = (getattr(options, "aux_basis", "") or "").strip()
if not aux_name:
orbital_name = getattr(basis, "name", "") or ""
if not orbital_name:
raise ValueError(
"run_cpcm_scf: density_fit/cosx requested but the "
"auxiliary basis could not be autodetected (orbital "
"BasisSet has no name). Set options.aux_basis "
"explicitly (e.g. 'def2-svp-jk')."
)
aux_name = default_aux_basis_for(orbital_name, kind="jk")
aux = BasisSet(molecule, aux_name)
if cosx:
# RIJCOSX — RI-J + seminumerical COSX-K on a quadrature grid.
cosx_grid_opts = getattr(options, "cosx_grid", None)
cosx_grid = (
build_grid(molecule, cosx_grid_opts)
if cosx_grid_opts is not None
else build_grid(molecule)
)
return make_cosx_jk_builder(basis, aux, cosx_grid)
return make_df_jk_builder(basis, aux)
def run_cpcm_scf(
molecule, # vibeqc.Molecule
basis, # vibeqc.BasisSet
*,
method: str = "rhf",
solvent: SolventSpec = "water",
options: Any = None,
xc_grid_options: Any = None,
progress_callback: Optional[Callable[[dict], None]] = None,
output=None,
) -> SolventResult:
"""Self-consistent CPCM SCF for a molecular system.
Parameters
----------
molecule
:class:`vibeqc.Molecule` — atomic structure + charge + spin.
basis
:class:`vibeqc.BasisSet`.
method : {"rhf", "uhf", "rks", "uks"}
Underlying SCF method. KS methods (``rks`` / ``uks``) require
the corresponding ``options`` to carry a ``functional``.
solvent
Anything accepted by :func:`resolve_solvent` — preset name,
ε value, dict, or :class:`SolventModel`. ``None`` /
``"vacuum"`` short-circuits to a plain gas-phase SCF
(still returned as a :class:`SolventResult` with
``e_solv = 0`` for uniform downstream handling).
options
Method-specific options struct (``RHFOptions`` / ``UHFOptions`` /
``RKSOptions`` / ``UKSOptions``). Default uses each method's
bundled default.
xc_grid_options
Optional :class:`GridOptions` for the KS numerical integration
grid. Defaults to the options object's own ``.grid``.
progress_callback
Optional callable invoked once per macro-iteration with a
diagnostic dict (``iter``, ``e_solv``, ``delta_e_solv``,
``total_q``, ``scf_iters``). Useful for live-logging in
:func:`vibeqc.run_job` and notebook progress bars.
Returns
-------
SolventResult
Final SCF result + solvation diagnostics. ``result.energy``
is the *total* in-solvent energy (gas-phase electronic +
nuclear + solvation); ``result.e_solv`` is just the
polarisation contribution.
"""
# ---- Gas-phase short-circuit ----------------------------------------
sm = resolve_solvent(solvent)
if sm is None or sm.is_gas_phase:
return _gas_phase_solvent_result(
molecule, basis, method=method, options=options,
xc_grid_options=xc_grid_options, solvent_model=sm,
)
method_label, gas_runner, jk_runner, needs_grid = _resolve_method_runners(method)
open_shell = method_label in ("uhf", "uks")
# ECP + CPCM is not supported in v0.9.0 — guarded here to fail
# loud rather than return a silently-wrong solvation energy.
#
# The macro-iteration builds Hcore in Python. Correctly coupling
# an ECP system to CPCM needs two pieces, not one: (a) the ECP
# one-electron operator folded into Hcore, and — the subtle part —
# (b) the cavity electrostatic potential V_nuc must use the
# ECP-reduced *effective* nuclear charge (Z − n_core) instead of
# the bare atomic number. With the bare Z the apparent surface
# charge sees a spuriously large net solute charge: e.g. Zn²⁺ with
# the ecp10mdf ECP (n_core = 10) has Z_eff = 20 and 18 valence
# electrons → physical net charge +2, but bare-Z bookkeeping gives
# +12, and since E_solv ∝ q_net² the solvation energy comes out
# (12/2)² = 36× too large. Resolving n_core robustly (it is
# per-element for lanl2dz) is a v0.9.x item.
if options is not None and list(getattr(options, "ecp_centers", []) or []):
raise NotImplementedError(
"run_cpcm_scf: ECP + CPCM is not supported in v0.9.0. The "
"cavity electrostatic potential would use the bare atomic "
"number instead of the ECP-reduced effective nuclear "
"charge (Z − n_core), giving a solvation energy wrong by "
"(Z_bare / Z_eff)². Run the ECP system gas-phase, or use a "
"full all-electron basis with CPCM. Tracked for v0.9.x."
)
# If density-fitting / RIJCOSX was requested, make sure the
# auxiliary basis is resolved *before* the gas-phase SCF — the
# high-level run_rhf/run_rks/... refuse density_fit=True with an
# empty aux_basis. Autodetect from the orbital basis name and fill
# it in so the gas SCF, the macro-iteration inner SCFs, and the
# JKBuilder all see one consistent aux basis.
if options is not None:
_ensure_cpcm_aux_basis(options, basis)
# ---- Step 1: gas-phase SCF ------------------------------------------
if options is None:
gas_result = gas_runner(molecule, basis)
else:
gas_result = gas_runner(molecule, basis, options)
if not gas_result.converged:
raise RuntimeError(
"run_cpcm_scf: gas-phase SCF did not converge — CPCM "
"macro-iteration is meaningless from a non-stationary "
"density. Address the gas-phase convergence problem first "
"(level_shift / damping / better guess)."
)
e_gas = float(gas_result.energy)
D = _density_of(gas_result, open_shell=open_shell)
# ---- Step 2: build cavity + cached operators ------------------------
atoms = list(molecule.atoms)
pos = np.array([list(a.xyz) for a in atoms], dtype=np.float64)
Zs = np.array([int(a.Z) for a in atoms], dtype=int)
cavity = build_cavity(
atom_positions_bohr=pos,
atom_numbers=Zs,
radii=sm.radii,
radii_scale=sm.radii_scale,
solvent_probe_radius_ang=sm.solvent_probe_radius_ang,
n_points_per_sphere=sm.n_points_per_sphere,
switching_sigma_bohr=sm.switching_sigma_bohr,
)
A = build_A_matrix(cavity.points, cavity.weights)
M_stack = _build_point_charge_operators(basis, cavity.points)
V_nuc_cav = _nuclear_potential_at_cavity(pos, Zs.astype(float), cavity.points)
# ---- Step 3: macro-iteration loop -----------------------------------
# Precompute one-electron pieces of Hcore once. The libint Hcore =
# T + V_nuc + (V_ECP, if any) is rebuilt per macro-iter only to
# add the solvent V_q correction — both T and V_nuc are
# geometry-tied and don't change as q varies.
from vibeqc import compute_overlap, compute_kinetic, compute_nuclear
S = np.asarray(compute_overlap(basis))
T = np.asarray(compute_kinetic(basis))
V_nuc_op = np.asarray(compute_nuclear(basis, molecule))
Hcore_gas = T + V_nuc_op
E_nuc = float(molecule.nuclear_repulsion())
n_electrons = int(molecule.n_electrons())
opts_for_jk = options if options is not None else _default_options(method_label)
# JKBuilder reused across macro-iterations (the basis doesn't move).
# Selected off the method options: density_fit → RI-JK, +cosx →
# RIJCOSX, otherwise direct (the v0.9.0 default).
jk_builder = _resolve_cpcm_jk_builder(basis, molecule, opts_for_jk)
# XC grid for KS methods — built once per geometry. RKS/UKS need the
# ``Grid`` C++ object; build via ``vibeqc.build_grid``.
xc_grid = None
if needs_grid:
from vibeqc import build_grid
grid_opts = xc_grid_options
if grid_opts is None and options is not None and hasattr(options, "grid"):
grid_opts = options.grid
xc_grid = build_grid(molecule, grid_opts) if grid_opts is not None \
else build_grid(molecule)
history: list[dict[str, float]] = []
prev_e_solv = 0.0
converged = False
last_inner_result = gas_result
last_cpcm: Optional[CPCMResult] = None
for macro_iter in range(1, sm.max_macro_iter + 1):
# ESP at cavity points from the current density.
V_elec = _density_potential_at_cavity(D, M_stack)
V_tot = V_elec + V_nuc_cav
# Solve CPCM: A q = -f(ε) V.
cpcm = solve_apparent_charges(
A, V_tot, epsilon=sm.epsilon, variant=sm.variant,
)
last_cpcm = cpcm
# V_q matrix and corrected Hcore.
V_q = _fock_solvent_contribution(cpcm.q, M_stack)
Hcore = Hcore_gas + V_q
# Inner SCF on the corrected Hcore. Seed with the current
# density so DIIS bootstraps from where we left off — converges
# in ~2–4 iters typically once the macro-loop is close.
opts = options if options is not None else _default_options(method_label)
if method_label == "rhf":
inner = jk_runner(
basis, n_electrons, S, Hcore, E_nuc, jk_builder, opts, D,
)
elif method_label == "uhf":
n_alpha, n_beta = _alpha_beta_counts(molecule)
D_a = np.asarray(last_inner_result.density_alpha)
D_b = np.asarray(last_inner_result.density_beta)
inner = jk_runner(
basis, n_alpha, n_beta, S, Hcore, E_nuc, jk_builder, opts,
D_a, D_b,
)
elif method_label == "rks":
inner = jk_runner(
basis, n_electrons, S, Hcore, E_nuc, jk_builder, xc_grid,
opts, D,
)
else: # uks
n_alpha, n_beta = _alpha_beta_counts(molecule)
D_a = np.asarray(last_inner_result.density_alpha)
D_b = np.asarray(last_inner_result.density_beta)
inner = jk_runner(
basis, n_alpha, n_beta, S, Hcore, E_nuc, jk_builder, xc_grid,
opts, D_a, D_b,
)
# Pull fresh density for the next macro-iteration.
D = _density_of(inner, open_shell=open_shell)
last_inner_result = inner
delta_e = cpcm.e_solv - prev_e_solv
record = {
"iter": macro_iter,
"e_solv": cpcm.e_solv,
"delta_e_solv": delta_e,
"total_q": float(cpcm.q.sum()),
"scf_iters": int(inner.n_iter),
}
history.append(record)
if progress_callback is not None:
progress_callback(record)
prev_e_solv = cpcm.e_solv
if abs(delta_e) < sm.tol_e_solv:
converged = True
break
assert last_cpcm is not None # macro_iter loop runs at least once
# ---- Standard CPCM total-energy decomposition -----------------------
# The SCF returned ``last_inner_result.energy`` is the SCF energy that
# includes ``V_q`` in the Fock matrix:
# E_SCF^{w/V_q} = E_HF^{gas}[D^{solv}] + Tr(D^{solv} V_q)
# = E_HF^{gas}[D^{solv}] + q^T V_elec
# The standard CPCM total energy (Cossi-Scalmani-Mennucci-Tomasi 2003
# eq. 3, matches PySCF / Gaussian / ORCA convention) is
# E_tot^{PCM} = E_HF^{gas}[D^{solv}] + (1/2) q^T V_tot
# so the "added" piece beyond the in-solvent SCF return value is
# ΔE = (1/2) q^T V_tot − q^T V_elec
# = (1/2) q^T V_nuc − (1/2) q^T V_elec.
# This guards against the common subtle bug where the SCF's already-
# included Tr(D V_q) term is double-counted by adding the full
# ½ q·V on top.
q_dot_Velec = float(np.dot(last_cpcm.q, last_cpcm.V - V_nuc_cav))
e_gas_at_D_solv = float(last_inner_result.energy) - q_dot_Velec
total_energy = e_gas_at_D_solv + last_cpcm.e_solv
# Optional citation siblings — fire the CPCM solvation papers
# (Klamt-Schüürmann 1993 + Cossi-Rega-Scalmani-Barone 2003 +
# Scalmani-Frisch 2010) when the standalone driver was invoked
# with a non-None ``output=`` stem. ``run_job`` already handles
# this for its own pipeline via runner.py.
if output is not None:
try:
from vibeqc.output.citations import emit_citations
_func = getattr(options, "functional", None) if options is not None else None
emit_citations(
output,
method=method,
basis=basis.name,
functional=_func,
uses_cpcm=True,
solvent_variant="cpcm",
)
except Exception:
# Best-effort — never let citation writer failure tank a
# converged CPCM run.
pass
return SolventResult(
scf=last_inner_result,
energy=total_energy,
e_solv=last_cpcm.e_solv,
e_gas=e_gas,
epsilon=sm.epsilon,
solvent_name=sm.name,
cavity=cavity,
cpcm=last_cpcm,
n_macro_iter=len(history),
macro_history=history,
converged=converged,
)
class _SolventAwareSCFResult:
"""Attribute-forwarding wrapper around a pybind11 SCF result.
The C++ result types (``RHFResult`` / ``UHFResult`` / ``RKSResult``
/ ``UKSResult``) are exposed via pybind11 ``def_readonly`` — Python
cannot set attributes on them. To keep
``run_job(..., solvent=...)`` returning something that downstream
code can use exactly like the underlying SCF result *and* exposes
the new solvation diagnostics, wrap the inner result and forward
attribute access via ``__getattr__``.
Exposed extra attributes:
* ``solvent_result`` — the full :class:`SolventResult`.
* ``e_solv`` — polarisation energy ``½ q·V`` (Hartree).
* ``energy_in_solvent`` — standard Cossi-Scalmani total
``E_HF^gas[D^solv] + ½ q·V_tot``.
The ``.energy`` attribute still resolves to the C++ result's
gas-phase electronic-plus-nuclear field (i.e. the
``E_SCF^{w/V_q}`` returned by the inner driver, which is what the
standard ``.out`` writer expects). Use ``energy_in_solvent`` for
the in-solvent total.
"""
__slots__ = ("_inner", "solvent_result", "e_solv", "energy_in_solvent")
def __init__(self, scf_result, solvent_result):
# Use object.__setattr__ to bypass our own __setattr__ guard.
object.__setattr__(self, "_inner", scf_result)
object.__setattr__(self, "solvent_result", solvent_result)
object.__setattr__(self, "e_solv", solvent_result.e_solv)
object.__setattr__(self, "energy_in_solvent", solvent_result.energy)
def __getattr__(self, name):
# __getattr__ is only called when the attribute isn't found
# on this object (or its class); forward to the wrapped result.
return getattr(self._inner, name)
def __setattr__(self, name, value):
# Disallow mutating the wrapped result via the proxy — keeps
# the C++ object's immutability semantics intact.
raise AttributeError(
f"_SolventAwareSCFResult is read-only; cannot set {name!r}"
)
def __repr__(self):
return (
f"_SolventAwareSCFResult("
f"inner={type(self._inner).__name__}, "
f"energy={self._inner.energy:.6f}, "
f"e_solv={self.e_solv:+.6f}, "
f"energy_in_solvent={self.energy_in_solvent:.6f})"
)
def _solvent_aware_scf_result(solvent_result):
"""Public constructor for :class:`_SolventAwareSCFResult`."""
return _SolventAwareSCFResult(solvent_result.scf, solvent_result)
def _alpha_beta_counts(molecule) -> tuple[int, int]:
"""Per-spin electron counts from a Molecule.
For multiplicity m (= 2S + 1) and total n_electrons n:
n_alpha − n_beta = m − 1
n_alpha + n_beta = n
Solving: n_alpha = (n + m − 1) / 2; n_beta = (n − m + 1) / 2.
"""
n_total = int(molecule.n_electrons())
mult = int(molecule.multiplicity)
n_alpha = (n_total + mult - 1) // 2
n_beta = n_total - n_alpha
return n_alpha, n_beta
def _default_options(method: str):
"""Construct each method's bundled default options struct."""
from vibeqc import RHFOptions, RKSOptions, UHFOptions, UKSOptions
m = method.lower()
if m == "rhf":
return RHFOptions()
if m == "uhf":
return UHFOptions()
if m == "rks":
return RKSOptions()
if m == "uks":
return UKSOptions()
raise ValueError(f"_default_options: unknown method {method!r}")
def _gas_phase_solvent_result(
molecule, basis, *, method, options, xc_grid_options, solvent_model,
) -> SolventResult:
"""Gas-phase short-circuit — uniform return type for downstream code.
Called from :func:`run_cpcm_scf` when ``solvent`` resolves to
``None`` or ``ε ≤ 1`` ("vacuum"). Builds an empty cavity (so
``result.cavity`` is still a valid :class:`CavityTessellation`
that callers can introspect) and a zero-charge :class:`CPCMResult`.
"""
method_label, gas_runner, _, _ = _resolve_method_runners(method)
open_shell = method_label in ("uhf", "uks")
if options is None:
gas_result = gas_runner(molecule, basis)
else:
gas_result = gas_runner(molecule, basis, options)
atoms = list(molecule.atoms)
pos = np.array([list(a.xyz) for a in atoms], dtype=np.float64)
Zs = np.array([int(a.Z) for a in atoms], dtype=int)
# Build a minimal (but valid) cavity with the same defaults that
# the in-solvent path would use, so downstream geometry inspection
# is consistent across the gas / solvent branch.
cavity = build_cavity(
atom_positions_bohr=pos, atom_numbers=Zs,
radii_scale=1.20, n_points_per_sphere=110,
)
empty_cpcm = CPCMResult(
q=np.zeros(cavity.n_points),
V=np.zeros(cavity.n_points),
e_solv=0.0,
epsilon=(solvent_model.epsilon if solvent_model is not None else 1.0),
)
name = solvent_model.name if solvent_model is not None else "vacuum"
return SolventResult(
scf=gas_result,
energy=float(gas_result.energy),
e_solv=0.0,
e_gas=float(gas_result.energy),
epsilon=(solvent_model.epsilon if solvent_model is not None else 1.0),
solvent_name=name,
cavity=cavity,
cpcm=empty_cpcm,
n_macro_iter=0,
macro_history=[],
converged=True,
)