We've had several questions regarding calculation and plotting of electron density.

The electrostatic potential represents the interaction between a point charge at a given point $\mathbf{r}$ in a system of atoms, given by: $$V(\mathbf{r})=\sum_{A=1}^{N_{\text {atoms }}} \frac{Z_{A}}{\left|\mathbf{R}_{A}-\mathbf{r}\right|}-\int \frac{\rho\left(\mathbf{r}^{\prime}\right)}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|} d \mathbf{r}^{\prime}$$

The first part of this is easy - it's the classical Coulomb formula between the atomic nuclear charges $Z_A$ at positions $R_A$ and the point $\mathbf{r}$.

The second piece involves an integral over the electron density $\rho$

I'd like to evaluate the electrostatic potential at a set of points, either over a grid or vertices of a surface mesh. Ideally, I'd like to see the algorithm as well as an open source implementation.

  • 1
    $\begingroup$ My not-so-hidden motive is to add an implementation to Avogadro after reading in an output file (e.g., Molden, fchk, etc.) $\endgroup$ Jul 6, 2021 at 18:34

2 Answers 2


If you know the Gaussian basis set and the density matrix, the task simplifies to the calculation of nuclear attraction integrals.

The electron density is expanded in the basis set as $$ n({\bf r}) = \sum_{\mu \nu} P_{\mu \nu} \chi_\mu({\bf r}) \chi_\nu({\bf r}).$$ When you substitute this into your equation, you get $$ V({\bf r}) = \sum_{A=1}^{N_{\rm atoms}} \frac {Z_A} {|{\bf r}-{\bf R}_A|} - \sum_{\mu \nu} P_{\mu \nu} \int \frac {\chi_\mu({\bf r}') \chi_\nu({\bf r}')} {{|{\bf r}-{\bf r}'|}} {\rm d}^3r'. $$

The nuclear attraction integral for atom $A$ is $$ V^A_{\mu \nu} = \int \frac {\chi_\mu({\bf r}) \chi_\nu({\bf r})} {|{\bf r}-{\bf R}_A|} {\rm d}^3r $$ so you can see that you just need to compute these integrals for every point on your grid; one-electron integrals are pretty cheap even though the nuclear-attraction integrals are the most expensive to compute.

Nuclear attraction integrals are available e.g. in Sun's libcint library described in J. Comput. Chem. 36, 1664 (2015) or Bin Gao's GEN1INT library described in Int. J. Quantum Chem. 111, 858 (2010).


The open-source program Multiwfn (http://sobereva.com/multiwfn/) contains a module that calculates the electrostatic potentials. According to Tian Lu, it has been quite thoroughly optimized for speed. Depending on your workflow, you may even be able to obtain the electrostatic potentials by calling Multiwfn in a black-box manner through command line, without having to dig into its source code.

  • $\begingroup$ While I like Multiwfn a lot, one of my hopes in asking the question was to have the discussion of "how would someone implement that formula," rather than "there are programs that can calculate it" $\endgroup$ Jul 9, 2021 at 4:30
  • $\begingroup$ My answer also answers that question - you can look into Multiwfn's source code... :) $\endgroup$
    – wzkchem5
    Jul 9, 2021 at 7:40

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .