For those kind of problems, you can choose to either use the method of images in real space or Fourier series. The former is efficient if the potential is short ranged, while the latter works better if the potential is smooth. The issue with the Coulomb potential is that it is neither due to the singularity at zero and the slow decay at infinity. In real space and Fourier space it is:
$$
G(x) = \frac{\Gamma(d/2-1)}{4\pi^{d/2}|x|^{d-2}} \quad \hat G(k) = \frac1{k^2}.
$$
To apply Ewald summation, you will therefore need to carefully split it in a short range part (which will be computed using the method of images) and long range, smooth contribution (which will be computed by Fourier series).
Method 1
You can physically motivate a short range term by adding a shielding length scale:
$$
\hat\phi(k) = \underbrace{\frac1{k^2+\kappa^2}}_{\text{short range: } \phi_s}+\underbrace{\frac{\kappa^2}{k^2(k^2+\kappa^2)}}_{\text{long range: } \phi_l}.
$$
The short range will give you a Yukawa potential:
$$
\Phi = \frac{K_{d/2-1}\left(\kappa|x|\right)}{(2\pi)^d|x|^{d/2-1}}
$$
and in $d=3$:
$$
\Phi = \frac{e^{-\kappa |x|}}{4\pi|x|}.
$$
Using Ewald summation:
$$
\phi(x) = \sum_{y\in\Lambda}\Phi(x-y)+\frac1{L^3}\sum_{k\in\Lambda^*-\{0\}}\frac{\kappa^2}{k^2(k^2+\kappa^2)}e^{ik\cdot x}.
$$
Notice that it only works for $d<4$, the limiting factor being the slow polynomial decay of the long-range term (it is not so smooth).
In practice, you will truncate the sums. If you truncate for $|y|<R$ and $|k|<K$, then the asymptotic ($R\gg L,K\gg1/L$) error will be of the order:
$$
\delta\phi \sim \frac{Re^{-\kappa R}}{\kappa L^3}+\frac{\kappa^2}{K}.
$$
In practice, your resources limit $R,K$, so you can optimise for $\kappa$ to minimise the error. The short range error is minimised when $\kappa\to\infty$ and the long range error is minimised when $\kappa\to0$, so a compromise needs to be found. You can work it out analytically in the case of large $R,K$:
$$
\kappa \sim \frac{\ln(KR^4/L^3)}R
$$
Method 2
The previous technique only works for $d<4$, the limiting factor being the Fourier series which decays only polynomially. An alternative, more efficient method is to use Gaussians instead:
$$
\hat\phi(k) = \underbrace{\frac{1-\exp(-\sigma^2k^2/2)}{k^2}}_{\text{short range: } \phi_s}+\underbrace{\frac{\exp(-\sigma^2k^2/2)}{k^2}}_{\text{long range: } \phi_l}.
$$
with $\sigma$ a characteristic lengthscale that you can tune. Back in real space, the short range potential is:
$$
\Phi(x) = \int \frac{1-\exp(-\sigma^2k^2/2)}{k^2}e^{ik\cdot x}\frac{d^dk}{(2\pi)^d}.
$$
You can compute the integral, or you can solve the Poisson equation:
$$
\frac1{r^{d-1}}\partial_r(r^{d-1}\partial_r\Phi) = \frac{\exp(-r^2/2\sigma^2)}{(2\pi\sigma^2)^{d/2}}
$$
In general, it is (using the upper incomplete gamma function):
$$
\Phi = \frac{\Gamma(d/2,r^2/2\sigma^2)}{2\pi^{d/2}(d-2)r^{d-2}}-\frac{\exp(-r^2/2\sigma^2)}{(2\pi)^{d/2}(d-2)\sigma^{D-2}}\underset{r\to\infty}\sim \frac{\exp(-r^2/2\sigma^2)}{(2\pi)^{d/2}\sigma^{d-2}(r/\sigma)^2}
$$
and in $d=3$:
$$
\Phi = \frac{\operatorname{erfc}(r/\sqrt 2\sigma)}{4\pi r} \underset{r\to\infty}\sim \frac{\exp(-r^2/2\sigma^2)}{(2\pi)^{3/2}\sigma(r/\sigma)^2}
$$
You can now express your periodic potential in terms of quickly, normally converging series (in any dimension):
$$
\phi(x) = \sum_{y\in\Lambda}\Phi(x-y)+\frac1{L^d}\sum_{k\in\Lambda^*-\{0\}}\frac{\exp(-\sigma^2k^2/2)}{k^2}e^{ik\cdot x}.
$$
As in the previous method, depending on the number of terms you can compute, you can similarly adjust $\sigma$ accordingly to minimise the error.
Hope this helps.
