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

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:

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

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 18

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 to Exercise 18

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, turns from electronic structure to the molecule’s shape itself — rotational constants, radial distribution, and surface/volume descriptors.