Implicit solvation (CPCM / COSMO)¶

vibe-qc v0.9.0 ships polarisable-continuum implicit solvation in the conductor variant (Cossi-Rega-Scalmani-Barone 2003; “C-PCM”). Aqueous chemistry is most users’ default reaction medium — solvent="water" is now a one-keyword drop-in for run_job.

Quick-start ¶

import vibeqc as vq

mol = vq.Molecule.from_xyz("ethanol.xyz")
result = vq.run_job(
    mol,
    basis="def2-svp",
    method="rks",
    functional="pbe0",
    solvent="water",          # ← the only new ingredient
    output="ethanol_aq",
)

sol = result.solvent_result   # SolventResult bundle
print(f"E_total (in water)   = {sol.energy:.6f} Ha")
print(f"E_total (gas-phase)  = {sol.e_gas:.6f} Ha")
print(f"ΔG_solv              = {sol.e_solv * 627.509:+.3f} kcal/mol")
print(f"Macro-iters          = {sol.n_macro_iter}")

solvent= accepts:

Input form	Meaning
`None` (default)	gas-phase SCF — same as omitting the kwarg.
`"water"`, `"dmso"`, `"acetonitrile"`, …	preset lookup → `vibeqc.SOLVENT_PRESETS`.
`"H2O"`, `"MeOH"`, `"DCM"`	preset aliases — case-insensitive.
`78.39`	numeric ε — custom dielectric.
`{"epsilon": 25.0, "variant": "cosmo"}`	dict → forwarded to :class:`vibeqc.SolventModel`.
`vq.SolventModel(epsilon=25, …)`	direct construction for full cavity control.
`"vacuum"` / `"none"` / `"gas"`	ε = 1 — short-circuits to gas-phase (uniform return).

Solvent preset table ¶

Static dielectric constants ε at 298 K, matching the Gaussian 16 SCRF=Solvent=… table to ±1 %. The full list of preset keys lives in python/vibeqc/solvation/presets.py.

Solvent	ε	Aliases
water	78.39	`h2o`
methanol	32.61	`meoh`
ethanol	24.85	`etoh`
dmso	46.83
dmf	37.22
acetonitrile	35.94	`mecn`, `acn`, `ch3cn`
acetone	20.49
thf	7.43
dichloromethane	8.93	`dcm`, `ch2cl2`
chloroform	4.71	`chcl3`
benzene	2.27
toluene	2.37
n-hexane	1.88	`hexane`
vacuum	1.0	`none`, `gas`

For a solvent outside this table, pass a numeric ε or build a custom :class:vibeqc.SolventModel:

ionic_liquid = vq.SolventModel(
    epsilon=15.6,
    name="1-butyl-3-methylimidazolium-BF4",
    radii_scale=1.25,           # tighter sphere scaling
    n_points_per_sphere=302,    # tighter Lebedev
    variant="cpcm",             # "cpcm" (default) or "cosmo"
)

What CPCM actually computes ¶

The solute sits in a molecular-shaped cavity built as the union of atom-centred spheres (Bondi vdW radii ×1.20). Outside the cavity, the solvent is treated as a homogeneous dielectric. The boundary carries an apparent surface charge \(q\) that satisfies

\[ A \mathbf{q} = -f(\varepsilon)\, \mathbf{V}(R_\text{nuc}, D), \qquad f(\varepsilon) = \frac{\varepsilon - 1}{\varepsilon} \]

where \(A\) is the discrete CPCM matrix (Scalmani-Frisch 2010 diagonal + \(1/|s_i - s_j|\) off-diagonal) and \(\mathbf{V}\) is the molecular electrostatic potential at the cavity tessellation points. The solvation energy is

\[ E_\text{solv} = \tfrac{1}{2} \sum_i q_i V_i \]

and the in-solvent total energy is \(E_\text{tot}^\text{solv} = E_\text{tot}^\text{gas} + E_\text{solv}\).

The full coupled problem \(\{q, D\}\) is solved by a macro-iteration loop (Cossi-Scalmani 2003 § II.B): converge SCF for the current \(V_q\), build a new \(V_q\) from the new density, repeat until \(|\Delta E_\text{solv}| <\) tol_e_solv (default 1 µHa). Three to five outer cycles is typical for small organic molecules in water.

Cavity tessellation ¶

Each atomic sphere is discretised by a Lebedev-Laikov quadrature; points sitting inside a neighbouring atomic sphere are smoothly removed via the Scalmani-Frisch continuous switching function (CSC; \(\sigma = 0.5\,a_0\)). This makes the cavity surface — and hence \(E_\text{solv}(R)\) — \(C^1\) continuous, which is essential for clean geometry optimisation.

Defaults that match Gaussian’s SCRF=PCM cavity:

Parameter	Default	Notes
vdW radii	Bondi 1964	per-element override via `radii=`
`radii_scale`	1.20	PCM / GEPOL convention
`solvent_probe_radius_ang`	0.0	0 = vdW-SES; set 1.385 Å for water SAS
`n_points_per_sphere`	302	Lebedev order 29 (Gaussian default)
`switching_sigma_bohr`	0.5	Scalmani-Frisch CSC width

Supported Lebedev point counts: 6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974. The PySCF / ORCA default is 110; the Gaussian default is 302; Q-Chem defaults to 194.

Inspecting the cavity ¶

from vibeqc.solvation import build_cavity

cav = build_cavity(
    atom_positions_bohr=...,    # (n_atoms, 3) bohr
    atom_numbers=...,           # iterable of Z
    radii_scale=1.20,
    n_points_per_sphere=302,
)

print(f"{cav.n_points} surface points after switching")
print(f"Total area: {cav.total_surface_area_bohr2:.2f} bohr²")

# Dump for visualisation
import numpy as np
np.savetxt("cavity.xyz", np.hstack([cav.points, cav.weights[:, None]]))

Geometry optimisation in solvent ¶

vibe-qc ships a closed-form analytic CPCM solvation gradient, :func:vibeqc.cpcm_gradient. The ASE calculator picks it up automatically when solvent= is set on the calculator:

from ase.optimize import BFGSLineSearch
from vibeqc.ase import VibeQC

atoms.calc = VibeQC(basis="def2-svp", functional="b3lyp",
                    solvent="water")
BFGSLineSearch(atoms).run(fmax=0.05)

The gradient of the in-solvent total energy \(E_\text{tot} = E_\text{HF}^\text{gas}[D^\text{solv}] + \tfrac{1}{2}q^\mathsf{T}V\) decomposes (envelope theorem at SCF + CPCM convergence) into four pieces, each evaluated in closed form:

the standard gas-phase analytic SCF gradient at the in-solvent density,
\(\tfrac{1}{2f}\,q^\mathsf{T}(\partial A/\partial R)\,q\) — the cavity self-interaction-matrix derivative,
\(q^\mathsf{T}(\partial V^\text{nuc}/\partial R)\) — the nuclear-ESP derivative,
\(q^\mathsf{T}(\partial V^\text{elec}/\partial R)\) — the electronic-ESP derivative, evaluated with a single libint nuclear-attraction-gradient pass over the cavity points (compute_external_charge_density_gradient).

No geometry rebuilds, no finite-difference step error. A finite-difference reference gradient, :func:vibeqc.cpcm_gradient_fd (full SCF + CPCM re-converge at \(6\cdot n_\text{atoms}\) displaced geometries), is retained as a cross-check oracle; cpcm_gradient(..., use_fd_electronic=True) likewise swaps just the electronic-ESP piece back to a finite-difference path for verification.

Choosing CPCM vs COSMO ¶

The variant= field selects the dielectric factor:

"cpcm" (default): \(f(\varepsilon) = (\varepsilon - 1)/\varepsilon\). Converges smoothly to the conductor limit at \(\varepsilon \to \infty\).
"cosmo": \(f(\varepsilon) = (\varepsilon - 1)/(\varepsilon + 0.5)\) (Klamt 1993 with \(x = 0.5\)).

The two agree to better than 1 % above \(\varepsilon \approx 30\) — for water and most polar solvents the choice is inconsequential (\(\Delta E_\text{solv} \lesssim 0.1\) kcal/mol). For very low-dielectric solvents (e.g. n-hexane, \(\varepsilon \approx 1.9\)) the choice can shift \(E_\text{solv}\) by ~5 %.

Density fitting and RIJCOSX ¶

The CPCM macro-iteration drives the inner SCF through whatever Fock-build path the method options select — there is nothing solvation-specific about the two-electron build. Set density_fit (and optionally cosx) on the method options and the CPCM SCF picks up RI-JK / RIJCOSX automatically:

import vibeqc as vq

opts = vq.RHFOptions()
opts.density_fit = True            # RI-JK
# opts.aux_basis = "def2-svp-jk"   # optional — autodetected if omitted
# opts.cosx = True                 # RIJCOSX (RI-J + chain-of-spheres K)

result = vq.run_job(mol, basis="def2-tzvp", method="rhf",
                    solvent="water", rhf_options=opts)

When aux_basis is left empty, it is autodetected from the orbital basis name via vibeqc.default_aux_basis_for(orbital, kind="jk") and filled in on the options struct (so the gas-phase bootstrap SCF, the macro-iteration inner SCFs, and the JKBuilder all share one aux basis). Density fitting changes the in-solvent total energy by the usual fitting error — ~1e-4 Ha for def2-svp-jk — and leaves E_solv essentially unchanged (~1e-6 Ha).

The CPCM analytic gradient currently uses conventional (non-DF) two-electron integral derivatives for the gas-phase SCF piece even when the SCF itself was density-fitted; the resulting DF/conventional mismatch is within the DF fitting error and below typical optimisation tolerances. A DF-consistent CPCM gradient is a future refinement.

Known limitations ¶

Molecules only — periodic CPCM is out of scope until v1.x; the conductor model needs adaptation to handle the periodic image problem cleanly.
No non-electrostatic terms — vibe-qc’s CPCM is currently electrostatic-only. The cavitation, dispersion, and repulsion (CDR) terms that a full SMD-style model adds are on the v0.10+ roadmap.

Theoretical references ¶

The implementation follows the standard CPCM literature:

Klamt, A. & Schüürmann, G. J. Chem. Soc. Perkin Trans. 2, 799 (1993) — COSMO conductor model.
Cossi, M., Rega, N., Scalmani, G. & Barone, V. J. Comp. Chem. 24, 669 (2003) — CPCM as the practical numerical variant; cavity + matrix layout used here.
Scalmani, G. & Frisch, M. J. J. Chem. Phys. 132, 114110 (2010) — continuous-surface-charge (CSC) formulation; the Scalmani-Frisch switching function and the \(A_{ii} = 1.0694\sqrt{4\pi / w_i}\) diagonal.
Tomasi, J., Mennucci, B. & Cammi, R. Chem. Rev. 105, 2999 (2005) — comprehensive PCM-family review.
Lange, A. W. & Herbert, J. M. J. Chem. Phys. 133, 244111 (2010) — modern Scalmani-style analytic CPCM gradient.

For citation in published work, see citation discipline.