# 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](../20-guide/properties/esp-surfaces.md)
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](../20-guide/properties/esp-surfaces.md) 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](linear-response-cphf.md) 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:

```python
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}
:label: ex-pcm-theory

1. 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.
2. 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.
3. 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](linear-response-cphf.md)'s own
   criterion for when an external potential needs a response contribution, explain this difference in one
   sentence.
:::

:::{solution} ex-pcm-theory
:class: dropdown

1. 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$.
2. 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.
3. 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](dispersion-theory.md) 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.