# Electrostatic potential surfaces

The **electrostatic potential (ESP)** is the closest thing to a "what a probe charge feels" map of a
molecule — it is what predicts where a nucleophile is welcomed and an electrophile is repelled, before any
orbital interaction is even considered. This chapter covers the ESP as a real-space field, its extrema, and
the surfaces mapped with it.

## Theory: the potential a point charge would feel

The molecular electrostatic potential at a point $\mathbf r$ is the classical Coulomb potential from the
**nuclei** (point charges) plus the **electron density** (a continuous charge distribution):

$$
V(\mathbf r) = \underbrace{\sum_A \frac{Z_A}{|\mathbf r-\mathbf R_A|}}_{\text{nuclear}}
\;-\;
\underbrace{\int \frac{\rho(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf r'}_{\text{electronic}} .
$$

The nuclear term is a trivial point-charge sum; the electronic term is exactly the same $1/r$ integral that
builds the nuclear-attraction matrix in the SCF itself (`int1e_rinv`, evaluated with the operator centred at
each *evaluation point* $\mathbf r$ rather than at a nucleus), contracted with the converged density — so
computing $V(\mathbf r)$ at a grid of points costs one such integral per point. $V(\mathbf r)>0$ means a
positive test charge is repelled (an electron-poor region — where a **nucleophile** is welcomed); $V(\mathbf
r)<0$ means it is attracted (electron-rich — where an **electrophile** attacks).

### Extrema and $\sigma$-holes

The chemically informative points are the ESP's **local extrema** on a surface (usually the van der Waals
isodensity surface, $\rho=0.001$ a.u.), found by evaluating $V$ on a dense triangulated mesh and applying
spatial non-maximum suppression (a point counts as an extremum only if it is the most extreme value within a
neighbourhood radius — this avoids reporting hundreds of mesh-noise "extrema"). A **minimum** marks a lone
pair, a $\pi$-system, or an anion's excess density; a **maximum** marks an electron-poor site — including the
famous **$\sigma$-hole**: an anisotropic *positive* region directly opposite a polarizable substituent's bond
(a halogen, say), responsible for halogen bonding.

## Usage

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r").run()

ext = qc.prop.esp.extrema(m)
ext["v_min"], ext["v_max"]     # global extrema, kcal/mol
ext["maxima"], ext["minima"]   # full lists: {value, position, atom, symbol}
```

For water, this finds **two O–H maxima** (≈ +44.8 and +44.6 kcal/mol, one behind each O–H bond — where a
nucleophile would approach) and **one O lone-pair minimum** (≈ −42.0 kcal/mol) — exactly the electrophilic
and nucleophilic sites organic-chemistry intuition already expects, now with numbers attached. A molecule with
a $\sigma$-hole (e.g. a C–Cl bond) shows a *positive* maximum on the extension of that bond axis, even though
chlorine is usually thought of as electron-rich — the ESP surface is what makes that counterintuitive
positive region visible and quantifiable.

### Surface summaries

`qc.prop.esp.surface(mychk)` gives whole-surface statistics rather than individual extrema:

```python
surf = qc.prop.esp.surface(m)
surf["v_max"], surf["v_min"], surf["v_avg"]   # the overall spread and average
```

`v_avg` and the spread of $V$ over the surface are used in some solvation and reactivity correlations
(they summarize "how polar the accessible surface is" in one or two numbers) — a coarser but cheaper
signal than the full extrema list.

## Worked example

```python
import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r").run()

ext = qc.prop.esp.extrema(m)
print(f"global: v_min={ext['v_min']:.2f}, v_max={ext['v_max']:.2f} kcal/mol")
for pt in ext["maxima"]:
    print(f"  max {pt['value']:.2f} kcal/mol near {pt['symbol']}{pt['atom']}")
for pt in ext["minima"]:
    print(f"  min {pt['value']:.2f} kcal/mol near {pt['symbol']}{pt['atom']}")
# global: v_min=-42.03, v_max=44.82 kcal/mol
#   max 44.82 kcal/mol near H1
#   max 44.61 kcal/mol near H2
#   min -42.03 kcal/mol near O0
```

:::{exercise}
:label: ex-esp

1. Why is the nuclear part of $V(\mathbf r)$ a simple point-charge sum while the electronic part needs an
   integral over $\rho$?
2. A $\text{CF}_3\text{Cl}$ molecule shows a *positive* ESP maximum directly opposite the C–Cl bond, even
   though chlorine is more electronegative than carbon. What is this feature called, and why does it not
   contradict chlorine being "electron-rich" overall?
3. You want one or two numbers that summarize "how polar is this molecule's accessible surface" without
   listing every individual extremum. Which `qc.prop.esp` call gives you that?
:::

:::{solution} ex-esp
:class: dropdown

1. The nuclei are point charges, so their contribution to $V$ is an exact closed-form Coulomb sum. The
   electrons are described by a continuous density $\rho(\mathbf r')$, so their contribution is a Coulomb
   *integral* over all space — the same $1/r$ kernel, just smeared over a distribution instead of points.
2. It is a **$\sigma$-hole** — an anisotropic region of positive potential on the extension of a polarizable
   atom's bond axis, caused by the anisotropic (not spherical) distribution of that atom's electron density.
   Chlorine's *overall* charge is still negative (as ordinary charge schemes report); the $\sigma$-hole is a
   *directional* depletion at one specific point on its surface, not a contradiction of its net charge.
3. **`qc.prop.esp.surface(mychk)`** — its `v_min`/`v_max`/`v_avg` summarize the whole surface in a few
   numbers, rather than `esp.extrema`'s full per-point list.
:::

The last properties chapter, [Geometric analysis](geometric-analysis.md), turns from electronic structure to
the molecule's shape itself — rotational constants, radial distribution, and surface/volume descriptors.