fastesp

Exact electrostatic potential from a Gaussian-basis wavefunction — evaluated analytically, and solved on the GPU.

fastesp computes

V(C) = sum_A Z_A/|C - R_A|  -  sum_{mu nu} P_{mu nu} < mu | 1/|r - C| | nu >

at an arbitrary set of points, without any real-space quadrature. The one-electron Coulomb integrals are evaluated in closed form (Hermite / McMurchie–Davidson expansion of the Gaussian product, Boys functions for the 1/r operator), so the result is exact to machine precision. It is a drop-in faster and more accurate replacement for the "electrostatic potential at a set of points" function of Multiwfn (menu 5 → 12 → 100).

Benchmark

54-atom protein fragment, ORCA wB97M-D3BJ/def2-TZVPPD (1450 spherical / 1634 Cartesian functions), points on the union of the Bondi vdW spheres. Multiwfn 8 CPU threads vs one RTX 5090, measured with benchmarks/benchmark.py under SLURM on a dedicated node.

* extrapolated; Multiwfn is linear in the point count at 86.2–87.6 ms/point. The part that replaces Multiwfn's integration is 0.030 ms/point vs 87.6 — a 2900x speedup — because the table build is point-independent.

Accuracy at the identical 300 points:

The 1.75e-7 a.u. gap is Multiwfn's quadrature error, and this is pinned down by three measurements: (1) the electron density — a pure basis evaluation on Multiwfn's side — agrees to 9e-10 relative, so both programs read the wavefunction identically; (2) Multiwfn's answer is bit-for-bit unchanged when its Becke grid is swept from 75×434 (1.8e6 points) to 200×2354 (2.5e7 points); (3) its three internal ESP engines (iESPcode 1/2/3) agree exactly with each other. See benchmarks/check_density_mwfn.py and benchmarks/mwfn_grid_convergence.py.

Install

pip install fastesp              # NumPy reference (exact, no toolchain needed)
fastesp build-native             # optional: compile the C++/CUDA backends
fastesp selftest                 # validate the analytic integrals

The native backends are compiled into ~/.cache/fastesp (override with FASTESP_CACHE) and need g++ (with OpenMP) and, for the GPU path, nvcc. pip install alone already gives exact results — just slower.

Use

# ESP at the points in pts.txt (Bohr, one "x y z" per line)
fastesp esp mol.molden.input --points pts.txt -o esp.txt
fastesp esp mol.molden.input --points pts.txt --backend numpy   # force reference

  backend      : cuda
  basis        : 1634 Cartesian functions, 2080 primitives
  pair table   :   0.2600 s  (point-independent)
  ESP kernel   :   0.0089 s
  ESP range    : -2.2055 .. +1.9854 V

As a library:

import numpy as np
from fastesp import read_molden, Basis, build_table

wfn    = read_molden("mol.molden.input")
basis  = Basis(wfn)
dens   = basis.density_matrix_cartesian()

table  = build_table(wfn, basis, dens, points_bohr)   # point-independent tables
esp    = table.esp()                                  # Hartree/e

examples/vdw_surface_esp.py shows the full production workflow: sample the Bondi vdW surface of a molecule and evaluate the ESP on it.

How it works

For a pair of primitive Cartesian Gaussians, the product is one Gaussian times a polynomial. Expanding about the Gaussian product centre P in Hermite Gaussians L_t gives the exact closed form

g_A^a g_B^b = K sum_t E_t L_t(r - P)
int L_t(r')/|R - r'| dr' = (-1)^|t| (d/dR)^t F(R),   F(R) = (2 pi / p) F_0(p R^2)

and (d/dR)^t F is expanded exactly by the one-dimensional recurrence a[t+1][j] = -a[t][j] + (t-2j+2) a[t][j-1].

The structural step that makes this fast: everything except R = C - P is independent of the sampling point. Contracting the E-coefficients, the (-1)^|t| signs and the (2p)^m factors over the Cartesian components and over the primitive pairs — weighted by the density matrix — yields, once per primitive pair, a sparse table G[alpha, m] such that

I_pair(C) = K * sum_{alpha,m} G[alpha,m] * R^alpha * F_m(p R^2)

Per point that is one Boys evaluation, a short monomial schedule and one dot product — a shape that maps directly onto GPU threads. The O(N_pairs) table build is independent of the point set (and of how many points there are), so it is done once and reused; only the O(N_pairs x N_points) evaluation runs per point set.

Screening uses the rigorous bound |I_pair| <= K * int|rho_pair| / d, with int|rho_pair| obtained from the |He_n| moment products computed at build time. The builder also prints the sum rule Tr(P S) = N_elec as a sanity check on the basis it just read.

MOLDEN conventions

The MOLDEN format does not pin these down. All three were determined from the file itself and are verified by Tr(P S) = N_elec and max|C^T S C - I|:

1. Primitive normalisation. ORCA writes GTO coefficients that already include N_s = (2 alpha/pi)^{3/4} (4 alpha)^{l/2} — the radial prefactor without the angular (2l-1)!! factor of a fully normalised Cartesian Gaussian. Verified against every single-primitive p/d/f shell in the test file to 8 digits.
2. Angular ordering. The spherical functions are the real solid harmonics in MOLDEN order m = 0, +1, -1, +2, -2, ... (fastesp/basis.py, _POLY).
3. Cartesian ↔ spherical transform. Rows are S-orthonormalised so that T S T^T = I.

If you have a wavefunction whose conventions differ, Basis + overlap_matrix will show it: Tr(P S) will not reproduce the MOLDEN occupancy sum.

Repository layout

Limitations

• Only the vacuum wavefunction ESP is implemented. Adding external point charges (e.g. a CPCM reaction field) is a one-line addition to the nuclear term.
• l = 4 (g) solid harmonics are tabulated but unverified — the test file has no g functions.
• The GPU solver is fp64 throughout; there is headroom for a mixed-precision build if needed.
• Tested against ORCA-generated MOLDEN files; other writers may use different conventions (see above).

Citation

If you use this in published work, please cite the repository (https://github.com/ZidongS/fastesp).

License

MIT — see LICENSE.