Help with deconvolving photon reads distributed across neighboring pixels
There is some work on resolving point sources from their smoothed signal called super-resolution, along with some elegant methods to solve it using convex-programming. The idea is that if you know that your signal consists of spikes (single points of energy), you can resolve them from the smoothed signal, even when the smoothing causes interference between spikes.

I used this method in my PhD to find the amplitude of 1D spikes from low-pass filtered signals, using the research published here (not mine). I didn't look at 2D signals like yours, but I recall that some work has been done on it, so you might be able to find some code.

For your method, I'm guessing the negative pixels occur because the smearing isn't exactly modelled by your PSF

What about something simple:


Assume that the smear occurs within a 3 x 3 radius.
Convolve the image with a 3 x 3 filter of all ones. This will add all the values in that area.
Find the points of local maxima in a 3 x 3 radius in the original image.
Replace the value of those points with those of the convolved image and zero the rest.