Electrostatic potential

Hi,

I am trying to compute the electrostatic potential of molecules on a grid. Unfortunately, all molecules have already been calculated and their density matrix stored in a database, but no EPS calculation was performed at the time. Thus I am making a fast GPU code that uses basis set information and density matrix to calculate it, assuming this is enough information to do it.

In order to test my code, I started a simple H atom calculation with PSI4, using a completely dull, but very debug-friendly, basis set:

psi4.basis_helper("""
   assign H testset

    [testset]
    spherical
    ****
    H     0
    S   1   1.00
          0.2              1.0000000
    S   1   1.00
          0.3              1.0000000
    ****
""")

and computed the ESP on the H atom, 1 bohr away and 2 bohr away, getting these values:
0 -0.8181848333
1 0.0289214457
2 0.0000480743

The first issue is the finite value in the origin. I always thought that the ESP at point r_p is the sum of electronic part (integral over all space of rho(r)/(r-r_p)) and the nuclear part (simply Z/r_p).
Shouldn’t this diverge in the origin? Is it computed differently in PSI4?

The second issue is that I cannot really get the same numbers.
In my code, the electronic part is computed as:

Sum_i,j DM[i,j] * Integrate[ GTO[i] GTO[j] / (r-r_p) d3r] * k

where DM is the density matrix, GTO are the basis set gaussians, k is a prefactor including the normalisation of angular and radial part of GTOs. The integral has the same mathematical form as a a nuclear attraction integral, where the nuclear position is the grid point where the ESP is being computed, and used those analytical formulas to get it (eq 25 from here pdf). For s orbitals it is fairly simple.
The numbers do no match the output from PSI4, and not just by a scaling factor that I may have forgotten. PSI4 potential decays fast between 1 and 2 bohr distance, while in my calculation it just halves.
So I suspect I am not even using the correct formula to begin with.
How is the EPS calculated in PSI4?

Any help is appreciated!

You may checkout the routine here: psi4/oeprop.cc at master · psi4/psi4 · GitHub

Looks like small deltaRs are screened out:

if (r > 1.0E-8) Vnuc += mol->Z(i) / r;

Thanks for the tip. I had a wrong sign in there, and now my code matches PSI4 output for that simple case.