Whenever one is doing grid-based calculations for particles interacting via the Coulomb potential, the singularity of the Coulomb potential $w(\mathbf r_1,\mathbf r_2) = \frac1{|\mathbf r_1 - \mathbf r_2|}$complicates the whole thing tremendously. Particularly for three-dimensional cartesian grids, whenever two gridpoint-triple $(x_i,y_j,z_k)$ coincide one gets infinity and that spoils any numerics.
I know of the following methods to avoid this unpleasant behavior (plus I add my experiences):
Cutoff parameters: introduce a suitable number $\alpha$ and replace the Coulomb potential by $\frac1{|\mathbf r_1 - \mathbf r_2 + \alpha^2|}$. This is very simple but it's not easy to find the "right" cutoff parameter $\alpha$ (meaning that the singularity is removed while the result is the same as for an exact method).
Applying the Poisson equation: the idea here is to evaluate the electrostatic potential of one particle in the given representation, then project on the other. The singularity somehow disappears in this process. This approach is generic, but rather costly (at least O(gridpoints))
Use some alternative representation for the Coulomb potential, such as Fourier or Laplace transform, other coordinate systems, projection on other basis functions etc. I have the impression these are often some kind of workarounds where the singularity is usually hidden in an infinite sum which is then truncated, or the like.
I'd like to do some large scale quantum many-body calculations on extensive grids, too large to apply accurate evaluation methods such as 2. or 3., still at the same time I want correct results. Is there any way to accomplish that?