calculate madelung constant by ewald method

I would like to understand the details of ewald summation. So, I took the calculation of madelung constant (Wikipedia) as an example. I took two method to calculate madelung constant. One is calculate it by expanding cubes method, the other is ewald summation method. A good description of ewald summation can be found here.

After the derivation, the total Coulomb interaction energy can be finally written as,

\begin{array} E E &=& E^S + E^L - E^{self} \\ \end{array}

\begin{array} &=& \frac{1}{4\pi\epsilon_0}\frac{1}{2} \sum_\mathbf{n} \sum_{i=1}^N \sum_{j=1}^N \frac{q_iq_j}{|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}L|} erfc\left( \frac{|\mathbf{r}_i - \mathbf{r}_j + \mathbf{n}L|}{\sqrt{2}\sigma}\right) \end{array}

\begin{array} + + \frac{1}{2V\epsilon_0} \sum_{\mathbf{k}\neq \mathbf{0}} \frac{e^{-\sigma^2k^2/2}}{k^2}|S(\mathbf{k})| -\frac{1}{4\pi\epsilon_0}\frac{1}{\sqrt{2\pi\sigma}} \sum_{i=1}^{N}q_i^2 \end{array}

In the paper I just mentioned, we can find all the details for this formula.

While, we need potential instead of total energy to calculate Madelung constant. We should take derivative over q to calculate the potential. Notice that some term exists which include 1/2. These are used to deal with double counting. When we want to calculate potential, these 1/2 should be erased.

\(test^2\)

madelung_realspace.m madelung_ewald.m