π density field engine

The page π electron density and pKa maps asks one question: when a proton leaves, where does its electron pair go? It answers with maps of the π electron density, its interference part, the shift between acid and base, the electric potential, and the individual π orbitals. Since version 16 it works for any molecule you type, not only the 20 reference molecules.

#Two sources, one format

Both come from the same code, the orbital layer, and use the same payload. The browser never knows which one it is drawing.

Routes: /pka-density (canonical), /pi-map and /learn/pka-density (301 redirects). Deep link to a molecule: /pka-density?smiles=Oc1ccc(cc1)C#N. Diagnostics: GET /api/pi-density-health.

#The payload of one form

Each acid or base form is a set of atoms and one entry per conjugated π component:

{
  "atoms": [{"el": "C", "x": 0.0, "y": 1.39, "z": 0.0, "Q": 0.062, "pi": true, "formal": 0}],
  "site": 6,
  "probe": [-0.914, 3.054, 0.0],
  "components": [{
    "atoms": [0, 1, 2, 3, 4, 5, 6],
    "orbitals": [{"atom": 0, "axis": [0, 0, 1], "type": "C"}],
    "basis": {"positions": [[0, 1.39, 0]], "radial": [[[2, 1.1438, 0.6976], [2, 2.3736, 0.3722]]], "M": [[1.15]]},
    "x": [2.1, 1.0, 1.0], "C": [[0.41]], "occ": [2, 2, 2],
    "homo": 3, "lumo": 4, "n_electrons": 8,
    "bond_orders": [[0, 1, 0.66]]
  }]
}

• axis is each p orbital's own direction from the conformer, so twisted and non-planar molecules are handled.
• radial gives each orbital's radial function as Slater terms [N, ζ, c]: the element's free-atom Hartree–Fock valence orbital, the same function the orbital layer and the 3D viewer use (the orbital layer, section 3). Payloads from 16.0 with zeta, n_star are still read as one Slater term.
• M \(= S^{-1/2}\) is the Löwdin matrix of that component.
• x, C, occ are the Hückel roots (\(E = \alpha + x\beta\)), coefficients and occupations.
• Q is the net atomic charge used for the electric potential; probe is where the proton sits (or sat).

#What the engine computes

static/js/pi-field-engine.js has no DOM dependency and exposes:

QOEngine.makeFrame     QOEngine.prepareForm   QOEngine.makeGrid
QOEngine.columnDensity QOEngine.sliceBasis    QOEngine.moMap
QOEngine.potentialMap  QOEngine.potentialAt   QOEngine.dipole
QOEngine.contour       QOEngine.linfit
QOEngine.volumeBox     QOEngine.volume        QOEngine.sampleVolume     (3D, since 16.1)

For an orthogonalised basis the density matrix is \(P_{rs} = \sum_k n_k c_{rk} c_{sk}\). With the Löwdin matrix the density in space is

and each orbital is \(\psi_k(\mathbf r) = \boldsymbol\phi(\mathbf r)^T M\,\mathbf c_k\). The atomic part keeps only the diagonal of P in the Löwdin basis; the interference part is the remainder,

The integral is zero because Löwdin orbitals are orthonormal: interference moves charge into the bonds without creating any (Ruedenberg 1962). Under every density map the page prints both integrals, so you can watch the sum rule hold.

#Any molecule, any orientation

The p functions use their own axes,

and the map is drawn in a view frame chosen from the molecule:

• flat π systems (π atoms in one plane, p axes along its normal): the plane itself. The column density through it is integrated over one half-space with 20 Gauss–Legendre nodes and doubled by symmetry.
• everything else: the best plane through the π atoms (their principal axes). The column density is integrated through that plane in four panels (−9…−3, −3…0, 0…3, 3…9 bohr, 48 nodes), because a tilted p lobe is no longer symmetric about it. The page says when this frame is in use.

Checks (16.1, Hartree–Fock radial functions): the column density integrates to the π electron count (phenol 8.000, 4-nitrophenol 12.000), the interference part to zero within 10⁻⁴, and a density value agrees with the Python layer to seven digits.

#Electric potential

from the atomic net charges, softened near the nuclei for display. V at the proton position is the descriptor V_probe of the pKa forecast.

#The same fields in 3D

QOEngine.volume(form, {kind}) evaluates π density, interference, one molecular orbital or the electric potential on a 3D box in model coordinates (at most about 110 000 points, step ≥ 0.16 Å). The 3D viewer uses it with the pi_field of /api/model, which comes from the same conformer as the atoms it draws, and turns the volumes into isosurfaces: π density at a fixed e/bohr³, interference as red (charge gained) and blue (charge taken) surfaces, orbitals in amber (ψ > 0) and teal (ψ < 0), and the potential as colour on the π density surface. Checks: the 3D density of phenol integrates to 7.995 of 8 π electrons inside the box, and a point value agrees with the Python layer to seven digits.

#Why the browser does the fields

A 168 × 168 map of a 10-atom π system needs millions of basis-function evaluations per view (28 224 pixels × 10 orbitals × 20–48 points through the plane). Done on the server, that cost would grow with every visitor and every slider move. In the browser it is the visitor's CPU, and the server's cost per molecule stays one orbital-layer run. The payload is small: the S and M matrices grow with the square of the π atoms in one component, not with the grid.

#Scope

The maps show a π-electron model: Hückel theory with heteroatom parameters on one conformer, without self-consistency, σ polarisation or solvent. The forecasts were fitted on a curated set (12 O–H acids, 5 amines) and come with their leave-one-out errors. They rank substituents; they do not replace a measurement. The page lists these limits next to the results.