How to Calculate Molecular Electrostatic Potential with Psi4 on a Cartesian Grid

Table of contents

• Scientific introduction
  • Electrostatic potential as a field
  • Nuclear and electronic contributions
  • Relation to reactivity and intermolecular interactions
• Method and basis-set choices
• Coordinate systems and units
• Test progress
• Demo user request
• Demo data
• Reproducibility
• Results and artifacts
• Interpretation and scientific use
• Limitations
• References
• Try this workflow

Molecular electrostatic potential, commonly abbreviated ESP, describes the interaction energy per unit positive test charge at positions around a molecule. It combines the attractive contribution of the nuclei with the repulsive contribution of the electron density, producing a spatial field that helps researchers reason about charge distribution, molecular recognition, hydrogen bonding, and sites that may attract electrophilic or nucleophilic partners.

This validated example uses Psi4 1.11 to calculate the RHF/STO-3G wavefunction of neutral singlet water at a fixed XYZ geometry and then evaluates the ESP at five explicitly supplied Cartesian points. The points are expressed in bohr and the potential is reported in hartree per elementary charge, equivalent to atomic units of electrostatic potential. At the symmetry-related points (2, 0, 0) and (−2, 0, 0) bohr, the potential is exactly −0.01573558401611841 hartree/e in both cases. At (4, 0, 0) bohr its magnitude decreases to 0.0021710524127136477 hartree/e, giving a far-to-near magnitude ratio of 0.13797088245913064.

Those observations establish two inspectable software and science checks: the potential respects the molecular reflection symmetry for the selected pair, and its magnitude decays away from the neutral molecule. They do not establish that RHF/STO-3G provides a quantitatively accurate electrostatic model for every chemical application. The basis is minimal, electron correlation is omitted, only five points were sampled, and no solvent or conformational ensemble was included.

Scientific introduction

Electrostatic potential as a field

For a molecule with nuclei at positions (R_A), nuclear charges (Z_A), and electron density (\rho(r)), the electrostatic potential at a field point (r) can be expressed schematically as the nuclear attraction term minus the electron-density contribution. The sign convention used by a program should always be checked, but the physical idea is simple: positive nuclei and negative electron density contribute oppositely to the potential experienced by a positive probe charge.

ESP is a function of position rather than one scalar property of the molecule. A single value has meaning only together with its coordinate, unit, molecular geometry, electronic method, and basis set. A useful calculation therefore preserves the grid and potential as paired data. Reporting a list of ESP values without coordinates, or coordinates without units, makes the result scientifically ambiguous.

The field can be sampled on arbitrary points, a regular Cartesian lattice, an isodensity surface, or atom-centred shells. Different choices answer different questions. A coarse volume grid can support visualization, while points on a molecular surface may be used for comparing exposed positive and negative regions. Charge-fitting methods select grids and constraints according to their own statistical and chemical assumptions.

Nuclear and electronic contributions

Near a nucleus, the nuclear term becomes strongly positive, while the electron density moderates and reshapes the field. Farther from a neutral molecule, positive and negative total charge cancel at leading order. The remaining long-range behavior is governed by multipoles such as the dipole and quadrupole, so orientation matters.

Water is neutral but polar. Its ESP is therefore not spherically symmetric and does not have the same value at all points with equal distance from the origin. The two selected x-axis points are symmetry-equivalent for the supplied orientation and should agree. The positive and negative z-axis points are not equivalent because the bent molecule has a directed dipole along that axis; the retained values are 0.041488310480275015 and −0.04510275307170053 hartree/e, respectively.

This distinction is important for validation. A test that demanded equality at all opposite points would be scientifically wrong. The correct symmetry assertion follows from the actual molecular plane and orientation. The test records which pair is expected to agree and treats the z-axis difference as meaningful rather than an error.

Relation to reactivity and intermolecular interactions

ESP maps are frequently used as qualitative descriptors of molecular charge distribution. Regions of negative potential may indicate locations favorable for interaction with positive probes, while positive regions may attract negative probes. Such interpretations can help discuss hydrogen-bond acceptors and donors, ligand recognition, or substituent effects.

However, ESP alone is not a complete reactivity theory. Steric accessibility, polarization, charge transfer, solvation, dispersion, orbital interactions, conformational populations, and reaction kinetics may dominate a real system. A visually intense color on an ESP surface is not a binding-energy measurement, reaction barrier, or probability of reaction.

Atomic partial charges can be fitted to reproduce an ESP, but those charges are model parameters rather than directly observable quantum-mechanical operators. Results depend on the fitting grid, restraint scheme, conformation, method, basis, and treatment of buried atoms. This five-point case validates direct potential evaluation; it does not perform a RESP, CHELPG, Merz–Kollman, or other charge fit.

Method and basis-set choices

The electronic state was calculated with restricted Hartree–Fock and the STO-3G basis. Restricted Hartree–Fock is appropriate for this closed-shell singlet fixture because the ten electrons occupy five doubly occupied spatial orbitals. It solves a self-consistent mean-field problem and includes exchange within a single determinant, but it omits dynamical electron correlation.

STO-3G is a minimal basis in which occupied atomic orbitals have very limited flexibility. This makes the calculation small enough for deterministic CPU testing, but minimal bases are generally weak choices for quantitative electrostatic modeling. Polarization functions allow electron density to deform in response to molecular geometry and bonding; diffuse functions can matter for anions and extended density. Larger, carefully selected basis sets and correlated or density-functional methods may substantially alter the ESP.

The total RHF/STO-3G energy was −74.96307837529535 hartree. That energy is retained as a wavefunction consistency check, not as the main ESP result and not as an experimental observable. Matching the expected energy range helps confirm that the supplied water geometry, charge, multiplicity, method, and basis were used before the field was evaluated.

Coordinate systems and units

The molecular geometry is read from data/water.xyz in angstrom:

3
Neutral water geometry in angstrom
O  0.000000  0.000000  0.000000
H  0.000000  0.757160  0.586260
H  0.000000 -0.757160  0.586260

The field points are read independently from data/esp-grid.csv in bohr. Keeping these two unit declarations explicit prevents a common and severe error: interpreting grid coordinates in angstrom when the API expects bohr. One bohr is approximately 0.529177 angstrom, so silently confusing the units changes every probe distance and the resulting potential.

The values are expressed as hartree per elementary charge. In atomic units this is a standard electrostatic-potential unit. Converting for another application requires a documented target convention. A visualization program may rescale values or apply a color range, and a force-field workflow may use different energy and charge units. The raw retained table should remain unchanged so derived views can be audited.

Coordinate orientation must also be preserved. Rotating the molecule and grid together leaves the physics unchanged, but rotating only one of them changes which molecular region is sampled. Translation likewise requires moving both geometry and points consistently. The test uses the original files without optimization or automatic reorientation of the requested grid.

Test progress

The feedback loop rejected two earlier chat attempts. Attempt 2 generated correct numerical values but placed a JSON summary in psi4-output.txt instead of the complete native report. The skill instructions were corrected to direct Psi4 itself to the declared output path. Attempt 3 then produced the genuine report but used the explicit CSV header esp_hartree_per_e, while the initial validator accepted only esp_au. Because these labels describe the same declared atomic potential, the validator was corrected without weakening the numerical checks, and a complete fresh attempt passed.

Demo user request

Load the neutral water geometry from data/water.xyz without changing its coordinates. Use Psi4 and RHF/STO-3G to evaluate the electrostatic potential at the five Cartesian points in data/esp-grid.csv, whose coordinates are in bohr. Save outputs/esp-grid.csv, outputs/esp-summary.json, and the complete native report at outputs/psi4-output.txt. Check the symmetry-equivalent x-axis pair and verify that the potential magnitude decays at the farther x-axis point.

The request is conversational, but it contains the scientific contract: source geometry, electronic method, basis, grid source, grid units, output paths, and expected invariant checks. The agent must execute the installed package rather than substitute generic sequence statistics, an empirical point-charge model, or fabricated values.

Demo data

The linked water.xyz and esp-grid.csv fixtures are small, synthetic, tracked, and free of private or restricted data. The grid has five rows:

Five points are enough to test input loading, units, orientation, symmetry, directional variation, output serialization, and a simple far-field relationship. They are not enough to render a smooth three-dimensional molecular surface or characterize every chemically relevant region. A production ESP map may require thousands or millions of points plus careful surface construction.

Reproducibility

Psi4 1.11 was installed from conda-forge in an isolated skill-owned Python 3.11 environment. The fresh-host solver plan resolved 101 packages totaling 447,057,426 download bytes, below the strict 500,000,000-byte local ceiling. The environment remains under the managed tool root, while downloaded package caches are retained in the skill test dossier for later reuse and separately reviewed offline distribution.

The calculation records the exact geometry and grid files, coordinate units, method, basis, total energy, potential units, program report, package version, operating system, architecture, and CPU-only execution. The CSV preserves full numeric precision rather than only rounded display values.

Reproduction should check file checksums before execution, confirm that Psi4 converged, and reopen every result after execution. The semantic validator asserts five finite values, equality of the x-axis pair within 1e-8, smaller magnitude at 4 bohr than at 2 bohr, a water energy within the expected regression interval, and the RHF, STO-3G, final-energy, and total-energy sections in the native report.

calculator = psi4.core.ESPPropCalc(wavefunction)
values = calculator.compute_esp_over_grid_in_memory(
    psi4.core.Matrix.from_array(points_in_bohr)
)

Users do not need to write this code in the workflow described at the end. It is shown to make the package-level operation transparent and to distinguish direct Psi4 evaluation from a hand-written Coulomb approximation.

Results and artifacts

The x-axis absolute symmetry difference is 0.0 at retained precision. The far-to-near magnitude ratio is 0.13797088245913064, so the 4-bohr value is substantially smaller than the 2-bohr value. These are deterministic fixture checks, not universal thresholds for other molecules or grids.

The first image is a focused application result rather than a file explorer or raw editor. The other two views are derived from retained artifacts. Their source paths, generators, timestamps, focused locator, and SHA-256 digests are recorded in the screenshot and visual-asset manifests.

Interpretation and scientific use

The supported conclusion is that the installed Psi4 workflow correctly evaluated the requested five-point ESP for the supplied water geometry and produced internally consistent artifacts. Symmetry and decay checks behaved as expected. The result also demonstrates that natural-language orchestration preserved separate angstrom and bohr inputs rather than silently merging their units.

To extend the analysis, choose a grid that matches the scientific question. A three-dimensional visualization requires sufficient spatial coverage and resolution. Surface comparison requires a reproducible molecular surface definition. Charge fitting requires a documented fitting algorithm, restraints, conformational strategy, and diagnostics. Comparing ligands may require consistent protonation, tautomer, geometry, method, basis, and orientation choices.

Negative ESP should not be translated directly into “negative atomic charge” at the nearest atom. The potential at a point is generated by the entire molecular charge distribution. Likewise, a local positive region does not prove an electrophilic reaction will occur there. Reactivity and binding involve more than static electrostatics.

Uncertainty is not estimated in this fixture. Numerical precision in the files reflects deterministic computation, not experimental certainty. Model-form error from RHF/STO-3G is likely much larger than the final printed digits. Applied work should compare reasonable electronic methods and basis sets, examine conformational sensitivity, and state whether solvent or environmental polarization is relevant.

Limitations

This test validates one neutral closed-shell molecule, one fixed conformation, one mean-field method, one minimal basis, five grid points, and CPU execution. It does not validate geometry optimization, correlated wavefunctions, density-functional approximations, open-shell systems, excited states, periodic electrostatics, implicit or explicit solvent, protein environments, partial-charge fitting, surface generation, volumetric cube files, or GPU acceleration.

The far-field check compares only two points and does not fit an asymptotic multipole model. Exact x-axis equality follows the chosen water orientation and should not be generalized to arbitrary coordinates. Cross-platform floating-point results may differ slightly and require declared tolerances. macOS, Windows/WSL, and ARM64 remain unverified for this exact end-to-end case.

References

1. Smith et al., “Psi4 1.4: Open-source software for high-throughput quantum chemistry,” Journal of Chemical Physics 152, 184108 (2020).
2. Psi4 official documentation.
3. Psi4 molecular-properties documentation.
4. Murray and Politzer, “The electrostatic potential: an overview,” Wiley Interdisciplinary Reviews: Computational Molecular Science.
5. Bayly et al., “A well-behaved electrostatic potential based method using charge restraints,” Journal of Physical Chemistry.

Try this workflow

MindPlot has built-in support for this scientific skill. You can attach molecular geometry and grid files, describe the desired electrostatic-potential calculation in ordinary language, and receive inspectable JSON, CSV, and native program outputs. Users do not need to write the Python shown above: the MindPlot agent writes and runs the calculation from the installed skill and presents the validated deliverables. Try it at mindplot.ai, or download the desktop version for a better experience and stronger local-data privacy.