I believe that you can make a small improvement, but not a dramatic one. Assume that your force loop consists of an outer loop over cells, and an inner one over nearest-neighbour cells (and then, more loops over the atoms within those cells). For the pair distribution function calculation, the inner loop can go over not just the nearest-neighbour cells, but over all cells covering a roughly spherical region out to the desired distance, relative to the cell being considered in the outer loop. You can construct a list of these cells (with indices relative to the central cell) at the start of the program. How you do this depends on how you handle the cell indexing, as you mentioned in your earlier question.
If the maximum distance of interest is not as large as half the box length $L$, this can make a reasonable saving of time over the all-pairs method. However, if you want to calculate out to half the box length, the maximum saving is the ratio of the spherical volume you are looking at, $\frac{4}{3}\pi (L/2)^3$ to the volume $L^3$ that you would look at with all pairs: around 50%. Essentially, you are just avoiding looking in the "corners" of the periodic box surrounding a given cell. For a system of size $N=1000$, I would think that it's not worth it, and you would do best to stick to all-pairs (discarding any for which the separation is larger than $L/2$).
The silver lining here is that the statistics of the calculation improve along with the extra work done. In a sense, you get what you pay for: you look at more pairs, for the higher values of $r$, so the cost goes up, but the precision of the calculation is better for large $r$. If you wish, you can offset this by doing the calculation less often (certainly not every step!).
There is another approach, which should give you the long-range part of $g(r)$ more cheaply, but again it is probably only worth doing if you have much larger systems. You can create a 3D histogram of particle density, $\rho(\mathbf{r})$ from a snapshot of the configuration, by counting the atoms in a cubic grid of cells, but on a smaller length scale than for the force calculation. The cell size should be a fraction of the particle diameter, typically. Then perform a 3D discrete Fourier transform, using an FFT algorithm, giving the instantaneous $\hat{\rho}(\mathbf{k})$ and take the square modulus. This function of wave-vector $|\hat{\rho}(\mathbf{k})|^2$ can be averaged over snapshots taken at intervals, and eventually gives the structure factor $S(\mathbf{k})$. Inverse-Fourier-transforming this would give you $g(r)$. Again, however, I think this would be too much trouble for $N=1000$ atoms.
