Electrostatic potential on atoms lead to divide by zero

[dspoel]

It seems that the following code could lead to an overflow (1/0) if the grid points coincide with the atoms.

void compute_grid_nucpot(double* coordinates, double* charges, long natom,
                         double* points, double* output, long npoint) {
    for (long ipoint=npoint-1; ipoint >= 0; ipoint--) {
        double tmp = 0.0;
        for (long iatom=natom-1; iatom >= 0; iatom--) {
            double delta = points[0]-coordinates[3*iatom];
            double dsq = delta*delta;
            delta = points[1]-coordinates[3*iatom+1];
            dsq += delta*delta;
            delta = points[2]-coordinates[3*iatom+2];
            dsq += delta*delta;
            tmp += charges[iatom]/sqrt(dsq);
        }
        *output += tmp;

        points += 3;
        output++;
    }
}

[PaulWAyers]

Yeah, ouch, even if it is the "right answer" it's nicer to handle the exception.

We do not actively maintain horton anymore, now that it is mostly being fragmented into bits and pieces. I'm not sure if gbasis (what we currently use for evaluating the electrostatic potential) has the same problem, but we should check that @marco-2023 and @FarnazH .

[dspoel]

I may be mistaken, but isn't this code used in Psi4 which is rather popular? That is how I found this out anyway.

[tovrstra]

It would be surprising if this code also exists in Psi4 as such. (A quick search did not reveal a match.) It can be found in Horton here:

https://github.com/theochem/horton/blob/master/horton/nucpot.cpp

I suppose the same feature is currently also implemented elsewhere with vectorized NumPy calls? It would be easier to deal with exception handling in such an implementation.