Solvation (CPCM/COSMO) handover — v0.9.0 shipped, v0.9.1 in progress (2026-05-20)

Status of the v0.9.0 implicit-solvation headline feature and the v0.9.1 follow-up work. Read this if you’re picking up the solvation chain.

What’s complete on origin/main

The full v0.9.0 CPCM/COSMO stack landed in commits a014833 (feature), cbe713f (tutorial + examples), 4e2abcd (4 rebuild- caught bug fixes), e4a8193 (tutorial numbers corrected to real runs). Verified end-to-end against a live C++ build.

Package: python/vibeqc/solvation/

Module

Role

cavity.py

S1a — Bondi vdW radii × 1.20, Lebedev-Laikov tessellation (6–974 pts/sphere; 302 default), Scalmani-Frisch CSC switching (σ = 0.5 bohr).

cpcm.py

S1b math — build_A_matrix (diagonal 1.0694 √(4π/w_i), off-diagonal `1/

driver.py

S1b coupling — run_cpcm_scf macro-iteration driver (RHF/UHF/RKS/UKS), SolventModel, SolventResult, _SolventAwareSCFResult attribute-forwarding wrapper.

gradient.py

S1c — cpcm_gradient (semi-analytic), cpcm_gradient_fd (full-FD oracle).

presets.py

S1d — 16 solvent presets + aliases.

Wiring

  • run_job(..., solvent="water") — top-level kwarg; attaches SolventResult to the returned object.

  • vibeqc.ase.VibeQC(solvent=...) — ASE calculator path; routes forces through cpcm_gradient.

Energy convention

E_tot = E_HF^gas[D^solv] + ½ q·V_tot (standard Cossi-Scalmani). The in-solvent SCF returns E_SCF^{w/V_q} = E_HF^gas[D] + q·V_elec; the driver decomposes correctly (driver.run_cpcm_scf).

Tests / docs

  • tests/test_solvation_cpcm.py — 25 tests, all green on a live build.

  • docs/user_guide/solvation.md — user-guide chapter.

  • docs/tutorial/39_solvation_water.md — worked tutorial, numbers verified against real runs (CH₂O/PBE0/def2-SVP).

  • examples/molecular/solvation/ — 3 runnable scripts + README.

Verified numbers (CH₂O / PBE0 / def2-SVP, vibe-qc v0.9.0)

  • Single-point hydration: E_solv = −5.011 kcal/mol, ΔG_solv = −4.122 kcal/mol, 7 macro-iters, 1047 cavity points.

  • Geometry opt: r(C=O) 1.1961 Å (gas) → 1.2027 Å (water), Δr = +0.0066 Å.

  • Dielectric scan: textbook (ε−1)/ε saturation; water −5.011, conductor limit reached ~95 % by ε ≈ 20.

v0.9.1 — Item 1 LANDED (this chat, 2026-05-20)

Item 1: fully analytic dV^elec/dR — DONE ✅

The v0.9.0 cpcm_gradient evaluated its electronic-ESP-derivative piece q·∂V^elec/∂R by finite difference. That last FD piece is now closed form. All three milestones complete:

  • Milestone A (C++) ✅ — compute_external_charge_density_gradient in cpp/src/gradient.cpp / cpp/include/vibeqc/gradient.hpp, bound in cpp/src/bindings.cpp (ExternalChargeGradient struct + the function). Drives libint’s nuclear-attraction gradient engine with the cavity points as fractional point charges; returns atom_grad (n_atoms, 3) + point_grad (n_points, 3) of E_ext = Σ_i q_i Tr(D·M_i). Translational invariance verified to machine precision (1e-18).

  • Milestone B (Python) ✅ — gradient._electronic_esp_gradient_analytic consumes the C++ pieces, routes point_grad[i] to cavity.point_atom[i], applies the CPCM sign. cpcm_gradient uses it by default; _electronic_esp_gradient_via_fd reachable via cpcm_gradient(..., use_fd_electronic=True).

  • Milestone C ✅ — 2 new tests in tests/test_solvation_cpcm.py (test_electronic_esp_gradient_analytic_matches_fd, test_cpcm_gradient_analytic_vs_fd_electronic_consistent); 27/27 solvation tests green. docs/user_guide/solvation.md + docs/tutorial/39_solvation_water.md updated; CHANGELOG “### Changed — CPCM solvation gradient is now fully closed-form”.

The CPCM gradient is now closed-form end to end. Verified against the FD oracle on H₂O/STO-3G: electronic-ESP kernel ~1e-4 Ha/bohr, full gradient ~5e-3 Ha/bohr (FD truncation dominates).

Sign accounting (the bug that bit the v0.9.0 ship — get it right):

  • M_i is the positive Coulomb kernel ⟨μ|1/|r−s_i||ν⟩.

  • V^elec_i = −Σ_μν D_μν M_i,μν (electron charge −1).

  • E_ext = Σ_i q_i Tr(D·M_i) = −q·V^elec.

  • The C++ function returns dE_ext/dR; cpcm_gradient needs q·dV^elec/dR = −dE_ext/dR.

  • libint’s Operator::nuclear gives −Σ q_i ⟨1/r⟩; the gradient engine’s contraction with D yields −dE_ext/dR, so the C++ function negates once to return +dE_ext/dR.

Build note for this worktree. third_party/ is symlinked from the parent checkout (ln -s ../../../third_party third_party — already in place; /third_party is gitignored). Per-worktree venv at .venv/. Rebuild after a C++ change with: .venv/bin/pip install --no-build-isolation -e .

v0.9.1 — Item 2 LANDED (this chat, 2026-05-20)

Item 2: DF / RIJCOSX in CPCM — DONE ✅

run_cpcm_scf no longer hard-codes make_direct_jk_builder. The macro-iteration inner SCF now selects the Fock-build path off the method options: density_fit → RI-JK, density_fit + cosx → RIJCOSX, else direct. New helpers _resolve_cpcm_jk_builder + _ensure_cpcm_aux_basis in driver.py (the latter autodetects aux_basis from the orbital basis name and fills it in so the gas-phase bootstrap, the inner SCFs, and the JKBuilder all share one aux). Verified on H₂O/def2-SVP (DF vs direct agree to 5e-5 Ha); tests test_cpcm_density_fit_matches_direct, test_cpcm_rijcosx_converges. Known follow-up: the CPCM analytic gradient’s gas-phase piece still uses conventional integral derivatives even under a DF SCF (mismatch within DF fitting error) — a DF-consistent CPCM gradient is a future refinement.

Item 3: ECP + CPCM — GUARDED, slips to v0.9.x ⚠️

Attempted (2026-05-20) and deliberately not shipped — the combination has a real charge-bookkeeping bug. run_cpcm_scf now raises NotImplementedError when options.ecp_centers is set, rather than returning a silently-wrong solvation energy.

The bug. Folding the ECP one-electron operator into the macro-iteration Hcore is necessary but not sufficient. The cavity electrostatic potential V_nuc must also use the ECP-reduced effective nuclear charge Z − n_core, not the bare atomic number. With bare Z the apparent surface charge sees a spuriously large net solute charge. Measured on Zn²⁺ / ecp10mdf / 6-31g / CPCM-water: the Molecule carries charge = 12 (ECP bookkeeping so n_electrons returns the 18-valence count), V_nuc used bare Z = 30, electrons = 18 → net charge the cavity saw = +12 instead of the physical +2. E_solv ∝ q_net², so the result came out (12/2)² = 36× too large: −22.65 Ha vs a Born estimate of ≈ −0.6 Ha.

What v0.9.x needs. Per-atom effective charge Z_eff = Z − n_core entering both _nuclear_potential_at_cavity and the gradient’s _nuclear_esp_gradient_contribution. n_core is robustly per (element, ECP family) — for ecpNNmdf it is NN, but lanl2dz varies per element, so it must be queried from libecpint rather than parsed from the library name.

Guard + test on main: driver.run_cpcm_scf (early NotImplementedError), test_cpcm_with_ecp_raises_not_implemented.

Item 4+: deferred v0.9.x / v0.10 follow-ups (NOT started)

  • DF-consistent CPCM gradient. The analytic gradient’s gas-phase SCF piece uses conventional 2e integral derivatives even when the SCF was density-fitted. Route _gas_phase_gradient through the DF gradient path when the SCF options carry density_fit. Small.

  • External-code parity validation. No CPCM E_solv has been compared against ORCA / Gaussian / PySCF CPCM. All current tests are internal-consistency (analytic-vs-FD, DF-vs-direct, dielectric monotonicity). A real parity number is needed before claiming benchmark-grade accuracy.

  • Non-electrostatic terms (SMD-style CDR). Cavitation + dispersion + repulsion. v0.10+ roadmap; substantial.

  • Periodic CPCM. Out of scope until v1.x (conductor model needs adaptation for the periodic image problem).

Files to consult

  • python/vibeqc/solvation/gradient.py — the gradient module being extended; _electronic_esp_gradient_via_fd is the function to replace, cpcm_gradient the dispatcher.

  • cpp/src/gradient.cpp::one_electron_gradient_contribution — the nuclear-attraction-gradient block (lines ~189–242) is the direct template for the new C++ function.

  • cpp/src/integrals.cpp::compute_1e_matrix — shows the libint nuclear engine taking a vector<pair<double, array<double,3>>> charge list (fractional charges already supported internally).

  • tests/test_solvation_cpcm.py::test_cpcm_gradient_matches_fd_water_in_water — the analytic-vs-FD regression that must stay green.

Standing rules (CLAUDE.md)

Rebase-then-push to main; never --force; never --no-verify; every commit ships green; new native deps need a roadmap conversation (this work adds none — pure libint). Commit per milestone, rebase before each push.