The canonical electrostatic potential, from the density (#129)

Completes the ESP work: "Density (ESP)" is now offered and computes the
potential a unit positive test charge feels,

    V(r) = sum_A Z_A/|r-R_A| - integral rho(r')/|r-r'| dr'

with the electronic term an integral over the electron density rather
than a sum over fitted point charges.

EspField holds the physics -- density matrix pairs, the electronic term,
the nuclear term -- deliberately outside the rendering code, because the
assembly is where the silent mistakes live: the sign of the electronic
term, the factor of two on off-diagonal density matrix elements, the
basis normalisation.  Any of those wrong gives a smooth, plausible,
wrong surface.

tests/slater checks it against charge neutrality, which is independent
of the shape of the density: far from a molecule the potential must
approach (nuclear charge - electron count)/R.  Neutral H and He reach
zero, a bare proton and He+ reach +1/R, hydride reaches -1/R, and an
H2-like pair reaches (2 - trace(PS))/R -- that last one pins the
off-diagonal factor of two specifically, since counting it once would
be off by exactly the overlap.

buildEspBasis() flattens the basis into one entry per orbital.  It walks
the same nesting as the field evaluation and must stay in step with it,
because what it counts is the column index into the MO coefficients.
Two walks of one structure is the shape this codebase has been bitten by
repeatedly, so it is CHECKED rather than trusted: the density rebuilt
from the flattened basis must reproduce the density the field loop
already computed, at forty points across the grid.  On disagreement it
declines and falls back to the point-charge map, and records by how much
-- a basis this cannot represent is a reason to fall back, not to take
the application down.

Fixed on the way, and this one was live: the density setup tested the
field type by its literal name,

    if (fieldType=="Density" || fieldType=="Spin Density")

so either ESP type skipped it entirely -- no MO range, no density array,
and a null dereference on the first grid point.  Now tested through the
flags.  The same literal comparison appeared three more times in
IsoSurfaceCmd deciding whether to draw a negative lobe; those are now
one boolean asked once, since three copies of a string comparison is how
one of them gets missed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFMUNbw2UpRPmgsqviAowb

Author: andy

include/tdat/EspField.H

@@ -0,0 +1,113 @@
+/**
+ * @file
+ *
+ * The electrostatic potential of a molecule, from its wavefunction.
+ *
+ *     V(r) = sum_A Z_A/|r - R_A|  -  integral rho(r')/|r - r'| dr'
+ *
+ * The electronic term is an integral over the electron density, not a
+ * sum over fitted point charges, which is what makes this the potential
+ * a unit positive test charge actually feels.
+ *
+ * Separate from the rendering code so the assembly -- the sign of the
+ * electronic term, the factor of two on off-diagonal density matrix
+ * elements, the basis normalisation -- can be tested against cases
+ * whose answer is known.  A mistake in any of those produces a smooth,
+ * plausible-looking, wrong surface.
+ */
+
+#ifndef ESPFIELD_HH
+#define ESPFIELD_HH
+
+#include <vector>
+using std::vector;
+
+/**
+ * One atomic orbital, flattened out of the basis set.
+ *
+ * The angular part is a sum of Cartesian monomials, which is what a
+ * spherical function is in this representation (2z^2-x^2-y^2 and so
+ * on), with any shell normalisation already folded into the
+ * coefficients.  Everything is in bohr.
+ */
+struct EspBasisFunction
+{
+   double center[3];
+
+   vector<int>    powerX;            // one entry per angular term
+   vector<int>    powerY;
+   vector<int>    powerZ;
+   vector<double> angularCoef;
+
+   vector<double> exponent;          // one entry per primitive
+   vector<double> contraction;
+
+   EspBasisFunction() { center[0] = center[1] = center[2] = 0.0; }
+
+   bool empty() const { return angularCoef.empty(); }
+
+   /** The value of this basis function at a point. */
+   double value(double x, double y, double z) const;
+};
+
+
+/** A nucleus: a charge at a point. */
+struct EspNucleus
+{
+   double center[3];
+   double charge;
+};
+
+
+class EspField
+{
+  public:
+
+    /**
+     * The significant density matrix pairs, as (mu, nu, weight).
+     *
+     * Built once and reused at every grid point.  Off-diagonal pairs
+     * appear once with their weight doubled, since the density matrix
+     * is symmetric, and pairs below the cutoff are dropped -- most of
+     * them, for anything but a small molecule, and each surviving pair
+     * costs a Coulomb integral at every single grid point.
+     */
+    struct Pairs
+    {
+      vector<size_t> mu;
+      vector<size_t> nu;
+      vector<double> weight;
+
+      size_t size() const { return mu.size(); }
+    };
+
+    /**
+     * Select the pairs worth computing from a density matrix.
+     *
+     * @param density  nbas x nbas, row major
+     * @param cutoff   relative to the largest element
+     */
+    static void selectPairs(const vector<EspBasisFunction>& basis,
+                            const vector<double>& density,
+                            double cutoff, Pairs& pairs);
+
+    /** The electron density at a point, from the selected pairs. */
+    static double density(const vector<EspBasisFunction>& basis,
+                          const Pairs& pairs,
+                          double x, double y, double z);
+
+    /** The potential at a point: nuclei minus electrons. */
+    static double potential(const vector<EspBasisFunction>& basis,
+                            const Pairs& pairs,
+                            const vector<EspNucleus>& nuclei,
+                            const double* point);
+
+    /**
+     * The electronic term alone, which is negative everywhere.
+     */
+    static double electronicPotential(const vector<EspBasisFunction>& basis,
+                                      const Pairs& pairs,
+                                      const double* point);
+};
+
+#endif

include/viz/ComputeMoCmd.H

@@ -10,6 +10,7 @@
 #include "viz/PropCmd.H"
 
 #include "viz/SGContainer.H"
+#include "tdat/EspField.H"
 
 /**
  * Field types naming an electron density surface coloured by the
@@ -39,6 +40,12 @@
 #define ESP_CHARGES_FIELD_TYPE "Density (ESP from atomic charges)"
 
 class TAtm;
+class TGBSAngFunc;
+class PropTable;
+class SGFragment;
+class TGBSConfig;
+class JCode;
+
 
 class TGBSAngFunc;
 class TGBSConfig;
@@ -82,13 +89,48 @@ class  ComputeMoCmd :  public  PropCmd
       float *makeField(int length);
 
       /**
-       * The electrostatic potential on the grid, for colouring the
-       * density surface just computed.  Returns false if the user
-       * interrupted it.
+       * The electrostatic potential from atomic point charges, for
+       * colouring the density surface just computed.  The cheap
+       * approximation.  Returns false if the user interrupted it.
        */
       bool computeEsp(SingleGrid *grid, vector<TAtm*> *atoms,
                       unsigned long gridRes,
                       float xDelta, float yDelta, float zDelta);
+
+      /**
+       * The basis set flattened into one entry per atomic orbital, in
+       * the same order as the MO coefficient columns.
+       *
+       * Returns false if the basis contains something this cannot
+       * represent, in which case no exact potential should be computed
+       * rather than one from a partial basis.
+       */
+      bool buildEspBasis(const SGFragment *sgfrag, TGBSConfig *gbsConfig,
+                         const JCode *code, TGBSAngFunc *angfunc,
+                         int maxShell, const int *length_shell,
+                         vector<EspBasisFunction>& basis);
+
+      /**
+       * The canonical electrostatic potential: nuclei, minus the
+       * integral over the electron density.  Returns false if the user
+       * interrupted it.
+       */
+      bool computeEspExact(SingleGrid *grid, vector<TAtm*> *atoms,
+                           const vector<EspBasisFunction>& basis,
+                           PropTable *moCoefs, const vector<double>& occ,
+                           PropTable *moCoefsBeta,
+                           const vector<double>& occBeta,
+                           const float *densityField,
+                           unsigned long gridRes,
+                           float xDelta, float yDelta, float zDelta);
+
+      /**
+       * How badly the flattened basis disagreed with the field
+       * evaluation, if it did.  Zero when the exact potential was
+       * computed; non-zero means it declined and the caller should say
+       * so rather than showing nothing.
+       */
+      double p_espBasisMismatch;
       void zeroField(float *field,int length);
       bool squareCheckInterrupt(const char *type, int mo, int max);
       bool moCheckInterrupt(const char *type, const char *atom,

src/apps/builder/MoPanel.C

@@ -268,16 +268,19 @@ void MoPanel::updateUIOptions()
       spinOk=true;
    }
 
-   //  The electrostatic potential is mapped onto a density surface, so
-   //  it needs the density AND per-atom charges.  Offered only when both
-   //  are there rather than being offered and then quietly producing a
-   //  plain surface.
-   bool espOk = densityOk && (expt->getProperty("ESPCHARGE") ||
-                              expt->getProperty("MULLIKEN"));
+   //  Both potentials are mapped onto a density surface, so both need
+   //  the density.  The canonical one integrates over that density and
+   //  so needs nothing further; the point-charge one needs charges, and
+   //  is offered only when they exist rather than being offered and
+   //  then quietly producing a plain surface.
+   bool espOk = densityOk;
+   bool espChargesOk = densityOk && (expt->getProperty("ESPCHARGE") ||
+                                     expt->getProperty("MULLIKEN"));
 
    ewxChoice *typewin = (ewxChoice*)FindWindow(ID_CHOICE_MO_TYPE);
    typewin->Clear();
-   if (espOk) typewin->Insert(ESP_CHARGES_FIELD_TYPE,0);
+   if (espChargesOk) typewin->Insert(ESP_CHARGES_FIELD_TYPE,0);
+   if (espOk) typewin->Insert(ESP_FIELD_TYPE,0);
    if (spinOk) typewin->Insert("Spin Density",0);
    if (densityOk) typewin->Insert("Density",0);
    typewin->Insert("MO",0);

src/tdat/chemistry/EspField.C

@@ -0,0 +1,151 @@
+#include <cmath>
+
+#include "tdat/EspField.H"
+#include "tdat/CoulombIntegrals.H"
+
+namespace {
+
+  //  Closer than this to a nucleus the point-charge term diverges while
+  //  the electronic term stays finite.  Clamped rather than skipped: a
+  //  single infinity would take the whole colour range with it, since
+  //  the ramp is scaled to the extreme value.
+  const double NUCLEUS_CLAMP = 1.0e-3;
+
+  /** The Coulomb integral between two basis functions at a point. */
+  double pairIntegral(const EspBasisFunction& a, const EspBasisFunction& b,
+                      const double* point)
+  {
+    double total = 0.0;
+
+    for (size_t ta = 0; ta < a.angularCoef.size(); ta++) {
+      const int la[3] = { a.powerX[ta], a.powerY[ta], a.powerZ[ta] };
+
+      for (size_t tb = 0; tb < b.angularCoef.size(); tb++) {
+        const int lb[3] = { b.powerX[tb], b.powerY[tb], b.powerZ[tb] };
+        const double angular = a.angularCoef[ta]*b.angularCoef[tb];
+        if (angular == 0.0) continue;
+
+        for (size_t ia = 0; ia < a.exponent.size(); ia++) {
+          for (size_t ib = 0; ib < b.exponent.size(); ib++) {
+            const double coef = a.contraction[ia]*b.contraction[ib];
+            if (coef == 0.0) continue;
+
+            total += angular*coef
+                     * CoulombIntegrals::potential(a.center, la,
+                                                   a.exponent[ia],
+                                                   b.center, lb,
+                                                   b.exponent[ib], point);
+          }
+        }
+      }
+    }
+    return total;
+  }
+}
+
+
+double EspBasisFunction::value(double x, double y, double z) const
+{
+  const double dx = x - center[0];
+  const double dy = y - center[1];
+  const double dz = z - center[2];
+  const double rr = dx*dx + dy*dy + dz*dz;
+
+  double radial = 0.0;
+  for (size_t i = 0; i < exponent.size(); i++) {
+    const double ar = exponent[i]*rr;
+    if (ar < 70.0) radial += contraction[i]*exp(-ar);
+  }
+  if (radial == 0.0) return 0.0;
+
+  double angular = 0.0;
+  for (size_t t = 0; t < angularCoef.size(); t++) {
+    angular += angularCoef[t]
+               * pow(dx, powerX[t]) * pow(dy, powerY[t]) * pow(dz, powerZ[t]);
+  }
+  return radial*angular;
+}
+
+
+void EspField::selectPairs(const vector<EspBasisFunction>& basis,
+                           const vector<double>& density,
+                           double cutoff, Pairs& pairs)
+{
+  pairs.mu.clear();
+  pairs.nu.clear();
+  pairs.weight.clear();
+
+  const size_t nbas = basis.size();
+  if (density.size() != nbas*nbas) return;
+
+  double largest = 0.0;
+  for (size_t i = 0; i < density.size(); i++) {
+    const double magnitude = fabs(density[i]);
+    if (magnitude > largest) largest = magnitude;
+  }
+  const double threshold = largest*cutoff;
+
+  for (size_t mu = 0; mu < nbas; mu++) {
+    if (basis[mu].empty()) continue;
+    for (size_t nu = mu; nu < nbas; nu++) {
+      if (basis[nu].empty()) continue;
+
+      const double value = density[mu*nbas + nu];
+      if (fabs(value) <= threshold) continue;
+
+      pairs.mu.push_back(mu);
+      pairs.nu.push_back(nu);
+      //  The density matrix is symmetric, so an off-diagonal pair
+      //  stands for two equal terms and is counted once, doubled.
+      pairs.weight.push_back((mu == nu) ? value : 2.0*value);
+    }
+  }
+}
+
+
+double EspField::density(const vector<EspBasisFunction>& basis,
+                         const Pairs& pairs, double x, double y, double z)
+{
+  double total = 0.0;
+  for (size_t i = 0; i < pairs.size(); i++) {
+    total += pairs.weight[i]
+             * basis[pairs.mu[i]].value(x, y, z)
+             * basis[pairs.nu[i]].value(x, y, z);
+  }
+  return total;
+}
+
+
+double EspField::electronicPotential(const vector<EspBasisFunction>& basis,
+                                     const Pairs& pairs, const double* point)
+{
+  double total = 0.0;
+  for (size_t i = 0; i < pairs.size(); i++) {
+    total += pairs.weight[i]
+             * pairIntegral(basis[pairs.mu[i]], basis[pairs.nu[i]], point);
+  }
+  //  Electrons are negative.
+  return -total;
+}
+
+
+double EspField::potential(const vector<EspBasisFunction>& basis,
+                           const Pairs& pairs,
+                           const vector<EspNucleus>& nuclei,
+                           const double* point)
+{
+  double total = 0.0;
+
+  for (size_t a = 0; a < nuclei.size(); a++) {
+    double r2 = 0.0;
+    for (int k = 0; k < 3; k++) {
+      const double d = point[k] - nuclei[a].center[k];
+      r2 += d*d;
+    }
+    double r = sqrt(r2);
+    if (r < NUCLEUS_CLAMP) r = NUCLEUS_CLAMP;
+    total += nuclei[a].charge/r;
+  }
+
+  return total + electronicPotential(basis, pairs, point);
+}