Implicit solvation theory: PCM

Implicit solvation theory: PCM#

Modeling a solute surrounded by an explicit shell of solvent molecules is expensive and adds substantial configurational sampling to an already-expensive calculation. The Polarizable Continuum Model (PCM) instead replaces the solvent with a structureless dielectric continuum outside a molecule-shaped cavity, and solves for the surface charge that continuum develops in response to the solute’s own electrostatic potential — a reaction field that in turn polarizes the solute itself, coupling into the SCF self-consistency loop. This chapter derives the boundary-integral electrostatics behind IEF-PCM and C-PCM, the GePol cavity that discretizes the molecular surface, and the SMD extension that adds non-electrostatic solvation physics — grounded in qc-rs’s pcm-core crate (a pure-Rust, PCMSolver-faithful port) and AGENTS.md’s pcm-core rules.

The boundary-integral picture#

PCM’s electrostatics reduces the full 3-D dielectric continuum problem to a 2-D one: instead of solving Poisson’s equation everywhere outside the cavity, it is enough to find the apparent surface charge \(q\) — a set of point charges on a discretized cavity surface — that reproduces the correct reaction potential the continuum would generate. The cavity surface is divided into small polygonal patches called tesserae, each carrying an area \(a_k\), a center, and an outward normal; the unknown is one apparent charge per tessera. Two classical boundary-integral operators act between tesserae \(i,j\) built from the vacuum Green’s function \(G(\mathbf r,\mathbf r')=1/|\mathbf r-\mathbf r'|\):

\[ S_{ij} = G(\mathbf r_i,\mathbf r_j), \qquad D_{ij} = \frac{\partial G}{\partial\mathbf n_i}\bigg|_{\mathbf r_j} = \hat{\mathbf n}_i\cdot\frac{\mathbf r_j-\mathbf r_i}{|\mathbf r_j-\mathbf r_i|^3} \ \ (i\ne j), \]

the single-layer (\(S\), an ordinary Coulomb-like kernel between tesserae) and double-layer (\(D\), its normal derivative — the field a point charge at tessera \(j\) projects along tessera \(i\)’s outward normal) operators — PCMSolver’s sign convention, which qc-rs matches exactly: \(\hat{\mathbf n}_{\text{probe}}\cdot (\text{source}-\text{probe})/|\text{source}-\text{probe}|^3\). The diagonal (\(i=j\)) elements need special treatment since the kernel is singular there; the collocation method qc-rs implements uses closed-form diagonal formulas instead of a naive (divergent) point evaluation,

\[ S_{ii} = c\sqrt{\frac{4\pi}{a_i}}, \qquad D_{ii} = -c\sqrt{\frac{\pi}{a_i}}\cdot\frac{1}{R_i}, \]

with \(c\approx1.07\) (the PCMSolver default collocation factor) and \(R_i\) the radius of the sphere tessera \(i\) belongs to — these diagonal formulas are what actually make collocation a robust, cheap default rather than requiring a genuine singular-integral quadrature at every tessera.

IEF-PCM: the general boundary-integral solve#

The Integral Equation Formalism (IEF) is the general PCM formulation, valid for both isotropic dielectrics and more exotic outside environments (ionic liquids, anisotropic media). For an isotropic dielectric of permittivity \(\varepsilon\) outside the cavity and vacuum (\(\varepsilon=1\)) inside, with \(A\) the diagonal matrix of tessera areas, the system matrices are

\[ \mathbf T(\varepsilon) = \Bigl[2\pi f(\varepsilon)\,\mathbf I - \mathbf D_i\mathbf A\Bigr]\,\mathbf S_i, \qquad f(\varepsilon) = \frac{\varepsilon+1}{\varepsilon-1}, \qquad \mathbf R_\infty = 2\pi\,\mathbf I - \mathbf D_i\mathbf A, \]

built from the inside (vacuum) Green’s-function operators \(S_i,D_i\), and the apparent surface charge solves the linear system

\[ \mathbf T(\varepsilon)\,\mathbf q = -\mathbf R_\infty\,\mathbf v, \]

where \(\mathbf v\) is the solute’s electrostatic potential sampled at every tessera (the molecular electrostatic potential from both nuclei and electron density, exactly the ESP this manual derives elsewhere, just evaluated at the cavity surface rather than a general grid). The factor \(f(\varepsilon)\) is what encodes how strongly the surrounding medium screens the solute’s field — as \(\varepsilon\to\infty\) (a perfect conductor), \(f(\varepsilon)\to1\) and the IEF formulation reduces smoothly toward the conductor-like limit derived in the next section, C-PCM. For a genuinely non-uniform outside environment (an ionic liquid, an anisotropic dielectric), the anisotropic generalization replaces the single Green’s function with separate inside/outside operators \(S_o,D_o\) (a strict superset of the isotropic case, not a different formulation):

\[ \mathbf T = (2\pi\mathbf I - \mathbf D_o\mathbf A)\mathbf S_i + \mathbf S_o(2\pi\mathbf I + \mathbf A\mathbf D_i^{\mathsf T}), \qquad \mathbf R = (2\pi\mathbf I - \mathbf D_o\mathbf A) - \mathbf S_o\mathbf S_i^{-1}(2\pi\mathbf I - \mathbf D_i\mathbf A), \]

with the same charge equation \(\mathbf T\mathbf q=-\mathbf R\mathbf v\). qc-rs restricts non-uniform outside Green’s functions to IEF-PCM specifically — C-PCM’s simpler derivation (below) assumes a uniform outside dielectric from the start and has no analogous anisotropic generalization.

C-PCM: the conductor-like limit#

C-PCM (also called COSMO) starts from a physically simpler picture: pretend the solvent is a perfect conductor rather than a finite-\(\varepsilon\) dielectric, solve the much simpler conductor boundary condition, then apply an empirical scaling correction to recover finite-\(\varepsilon\) behavior. On an ideal conductor surface the reaction potential exactly cancels the solute’s own potential, giving the unscaled conductor charge equation \(\mathbf S\mathbf q=-\mathbf v\); the empirical fix scales the single-layer matrix by a factor depending on \(\varepsilon\):

\[ \mathbf S_{\text{scaled}} = \frac{\mathbf S}{f_\varepsilon}, \qquad f_\varepsilon = \frac{\varepsilon-1}{\varepsilon+k}, \qquad \mathbf S_{\text{scaled}}\,\mathbf q = -\mathbf v, \]

with \(k\) a small empirical correction (qc-rs follows the PySCF convention: \(k=0\) for plain C-PCM, \(k=0.5\) for the closely related COSMO variant — the same underlying solver, just a different scaling constant). C-PCM’s derivation never needed the double-layer operator \(D\) at all — the conductor boundary condition is purely a single-layer (\(S\)) relation — which is both why it is markedly cheaper to assemble than IEF-PCM and why it has no natural extension to a non-uniform outside environment (there is no \(D_o/S_o\) split to generalize in the first place). For most solvents at ordinary dielectric constants, C-PCM and IEF-PCM agree closely (both are approximating the same physical reaction field, just via different closed-form routes), but they are not identical formulations and can diverge more for very small \(\varepsilon\) or unusual geometries — this is a real difference in solved equations, not merely two numerical routes to the same answer.

The GePol cavity: from atomic spheres to a tessellated surface#

Both solvers need a discretized cavity surface as input — a set of tesserae with positions, normals, and areas. qc-rs’s GePol cavity (a Rust port of the PEDRA algorithm from the PCMSolver project) builds this from a union of atom-centered spheres (radii from a standard table — Bondi or similar van der Waals radii, usually scaled by a factor \(\sim1.2\)) using a specific, deliberate construction: the algorithm always starts from the \(D_{2h}\) tessellation of a single sphere — it tessellates one-eighth of a sphere (one octant) directly, then uses the three reflection planes (\(Oyz\), \(Oxz\), \(Oxy\)) to replicate that one patch to the full sphere. This exploits Abelian point-group structure directly in the cavity construction algorithm itself, independent of whether the molecule as a whole happens to have any symmetry — it is a computational convenience for generating one sphere’s tesserae efficiently, not a statement that only symmetric molecules get a cavity. Overlapping spheres (from atoms close enough together that their van der Waals spheres intersect) need their intersection regions excluded so the final cavity encloses the whole molecule with no interior gaps or spurious re-entrant surface — a standard, well-established geometric step in GePol-family cavity construction that the tessellation and reflection procedure above must be combined with to produce a physically sensible molecular surface.

Coupling into the SCF: the reaction-field Fock term#

The apparent surface charge \(\mathbf q\) is not a static, one-time correction — it depends on the solute’s electron density (through \(\mathbf v\)), and the reaction field it produces polarizes that same density in turn, so PCM must participate in the SCF loop as a genuinely density-dependent Fock term, qc-scf’s PcmTerm. Each SCF cycle: build \(\mathbf v\) (the electrostatic potential at every tessera) from the current density and nuclei, solve for \(\mathbf q\) via whichever solver (\(\mathbf T\mathbf q=-\mathbf R_\infty\mathbf v\) for IEF-PCM, \(\mathbf S_{\text{scaled}}\mathbf q=-\mathbf v\) for C-PCM), then add the resulting reaction-field potential back into the Fock matrix via the same int1e_rinv-type tessera potential integrals the ESP chapter uses for evaluating a potential at an arbitrary point (here, evaluated at every tessera instead of a general grid). This is structurally the same seam pattern the linear-response chapter describes for any genuinely density-dependent environment term: unlike a fixed ECP operator, PCM’s reaction field responds to any density — including a transition density — so it must also contribute its own term to the response kernel \(\sigma\)-build whenever a stability analysis, CPHF solve, or (in principle) TDDFT calculation probes a density-dependent perturbation in a solvated system, not merely to the ground-state Fock build.

SMD: adding non-electrostatic solvation physics#

IEF-PCM and C-PCM capture only the electrostatic part of solvation — the reaction field from charge polarization. Real solvation free energies also have a substantial non-electrostatic contribution: the work needed to form a cavity in the solvent (cavitation), dispersion interactions between solute and solvent, and structural/exchange-repulsion effects at the solute-solvent interface. SMD (Solvation Model based on Density) captures this with a single additional term built from geometry alone:

\[ E_{\text{solv}} = \underbrace{E_{\text{electrostatic}}(\text{IEF-PCM}, R_{\text{SMD}})}_{\text{same IEF-PCM machinery, with SMD's own intrinsic-Coulomb radii}} + \underbrace{G_{\text{CDS}}}_{\text{cavity-dispersion- solvent-structure}}, \]

where \(G_{\text{CDS}}\) is a geometry-only scalar — a sum of atomic and molecular surface-tension-like parameters multiplied by a solvent-accessible surface area (SASA), with no additional Fock-matrix term and no additional SCF self-consistency beyond the electrostatic IEF-PCM piece it wraps. The electrostatic part genuinely reuses the SWIG IEF-PCM solver path unchanged — SMD’s distinctive contribution is entirely the radii it feeds into that shared machinery (its own intrinsic-Coulomb radii, distinct from the plain van der Waals radii IEF-PCM/C-PCM use by default) plus the purely geometric \(G_{\text{CDS}}\) correction added afterward.

Verified example — water, in vacuum and in three solvation treatments:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

m_gas = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r").run()
m_gas.scf.energy   # -74.96294665653868

m_ief = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r",
                    pcm={"solvent": "water", "model": "IEF-PCM"}).run()
m_ief.scf.energy, m_ief.scf.energy - m_gas.scf.energy
# (-74.9690349700267, -0.006088313488021413)

m_cpcm = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r",
                     pcm={"solvent": "water", "model": "C-PCM"}).run()
m_cpcm.scf.energy, m_cpcm.scf.energy - m_gas.scf.energy
# (-74.96907206979104, -0.006125413252362932)

m_smd = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r",
                    pcm={"solvent": "water", "model": "smd"}).run()
m_smd.scf.energy, m_smd.scf.energy - m_gas.scf.energy
# (-74.97017942182167, -0.007232765282992659)

All three treatments stabilize the energy relative to vacuum, as physically expected for a polar solute in a polar solvent. IEF-PCM and C-PCM agree closely but not identically (a \(\sim4\times10^{-5}\ E_h\) difference here) — exactly the “closely related but genuinely different equations” relationship derived above, not a numerical-precision artifact. SMD’s larger stabilization reflects the additional non-electrostatic \(G_{\text{CDS}}\) contribution on top of an electrostatic part built with its own, different radii — not merely a rescaled IEF-PCM number.

Exercise 14

C-PCM’s derivation never introduces the double-layer operator \(D\), only \(S\). Explain in one sentence why the conductor boundary condition (the physical starting point C-PCM’s derivation begins from) makes \(D\) unnecessary, whereas IEF-PCM’s finite-\(\varepsilon\) boundary condition cannot avoid it.
The GePol cavity construction always tessellates a single sphere’s \(D_{2h}\) (one octant, reflected three times), regardless of whether the actual molecule has any symmetry at all. Explain why this is not a contradiction — why building a symmetric single-sphere tessellation is compatible with producing a final cavity for a molecule with no point-group symmetry whatsoever.
PCM must contribute a term to the linear-response \(\sigma\)-build (not just the ground-state Fock build), while a fixed ECP operator does not. Using the linear-response chapter’s own criterion for when an external potential needs a response contribution, explain this difference in one sentence.

Solution to Exercise 14

An ideal conductor’s surface is, by definition, an equipotential — the reaction field the induced surface charge produces exactly and locally cancels the solute’s potential at every point of the surface, giving a boundary condition that involves only potentials (i.e. only the single-layer operator \(S\)). A finite-\(\varepsilon\) dielectric has no such simple equipotential condition — the correct boundary condition there involves matching both the potential and its normal derivative (the normal component of the displacement field) across the interface, which is exactly what pulls the double-layer operator \(D\) (a normal-derivative kernel) into the IEF-PCM system matrices \(T(\varepsilon)\) and \(R_\infty\).
The \(D_{2h}\)-of-one-sphere tessellation is a computational device for efficiently generating the tesserae of a single, individually symmetric object (a sphere is always symmetric regardless of what molecule it happens to be attached to) — it says nothing about how multiple such spheres, placed at the actual (possibly completely asymmetric) atomic positions of a real molecule, combine. The final cavity is the union of many individually-generated, individually-symmetric-sphere tessellations positioned at whatever atomic geometry the molecule actually has, including overlap trimming between neighboring spheres — the molecule’s own point group (or lack of one) is a property of where those spheres sit relative to each other, not a property of how any single sphere’s own surface happens to be tessellated.
The linear-response chapter’s criterion is that only density-dependent terms need a response contribution — the response kernel is a sum of “the response of each density-dependent term” to a change in density. A fixed ECP operator does not depend on the density at all (it is built once from the converged geometry and stays fixed across SCF iterations), so its response to any density change is exactly zero. PCM’s reaction field, by contrast, is defined by the solute’s electrostatic potential \(\mathbf v\), which is itself a density-dependent quantity — the reaction field responds to any density, including a transition density used in a stability check or a CPHF right-hand side, so it genuinely needs its own contribution wherever the response kernel is evaluated.

Solvation and the ECP are qc-rs’s two Hamiltonian-modification chapters in Part II; dispersion theory covers a third, geometry-only correction that composes with either (or neither) of them, needed because neither Hartree-Fock nor most density functionals capture long-range van der Waals attraction on their own.