Ewald Sums

In 1921 Ewald (79) proposed an efficient way to recast the summation 4.9 in two rapidly converging series. Here in order to derive in a systematic and controlled way the final result we consider a Yukawa potential $v(r)= e^{ - \epsilon \vert r\vert }/\vert r\vert$ and take the limit $\vert\epsilon\vert \to 0$ only in the final expression. Following the Ewald's idea we split the potential in two parts:

$$v(\vert r\vert) = v_{long} (\vert r\vert) + v_{short} (\vert r\vert)$$ (4.10)

where

$$v_{short} (\vert r\vert) = v(\vert r\vert) erfc( \sqrt{\alpha} \vert r\vert)$$ (4.11)

$$v_{long} (\vert r\vert) = v(\vert r\vert) - v_{short} (\vert r\vert)=erf(\sqrt{\alpha} \vert r\vert) v(\vert r\vert),$$ (4.12)

$erf$ is the error function and $erfc$ the complementary one. Notice that the long range part has several important properties:

$$\lim\limits_{r \to 0} v_{long} (r) = 2 \sqrt{ \alpha/\pi}$$ (4.13)

$$v_{long}(k) = 4 \pi/k^2 e^{ -k^2/(4 \alpha) } { \rm ~ for~ } \epsilon \to 0$$ (4.14)

$$v_{long} (k=0) = 4 \pi/\epsilon^2 {\rm ~~ finite~ only~for~~} \epsilon >0$$ (4.15)

On the other hand the short range potential decays very fast in real space and the sum converges very quickly. Since $U$ in Eq.(4.9) depends linearly on the potential, we can easily decompose two contributions: a short-range and a long range one. Then the latter can be more easily evaluated in Fourier space:

$$U = U_{short-range} + U_{long-range}$$

$$U_{short-range} = \frac{1}{2} \sum\limits_{\xi_i\ne \xi_j+R_s,R_s} q_i q_j v_{short} ( \vert\xi_i -\xi_j + R_s\vert )$$

$$U_{long-range} = \frac{1}{2} \sum_{\xi_i,\xi_j,R_s} q_i q_j v_{long} (\vert\xi_i -\xi_j+R_s\vert) - \frac{1}{2} \sum_{i} v_{long} (r \rightarrow 0) q_i^2$$

$$= \frac{1}{2 V} \sum_{ \vec k \ne 0} \sum\limits_{i,j} q_i q_j e^{ i \vec k (\xi_i-\xi_j) } v_{long} (k) - \sqrt{\frac{\alpha}{\pi}} \sum_{i} q_i^2$$

where $V$ is the volume of the unit cell and the sum over the momenta are on the discrete $\vec k$ values allowed by the periodicity $\vec k \cdot \vec R_s= 2 \pi n$ . In the latter expression we have used Eq.(4.15) and the fact that the charge neutrality $\sum_i q_i=0$ implies that the $\vec k=0$ term can be omitted in the sum for any $\epsilon>0$ . In this way the limit $\epsilon \to 0$ can be found consistently also for long range potentials by replacing expression (4.14) in the corresponding Fourier transform for $v_{long}$ . For a non neutral system instead the Ewald sum is divergent as expected.

For Coulomb interaction the potential energy becomes:

$$U = \frac{1}{2} \sum_{i,j}^N \sum_{R_s} \frac{q_i q_j}{\vec r_{ij} + \vec R_s} erfc \left( \sqrt{\alpha}(\vec r_{ij} +\vec R_s ) \right)$$

$$+ \frac{1}{2} \sum_{\vec k \neq 0} \sum_{i,j}^N \frac{4 \pi q_i q_j... ...\vec k \right\vert^2 /4\alpha} - \sqrt{ \frac{\alpha}{\pi}} \sum_{i=1}^N q^2_i.$$ (4.16)

In the potential energy 4.16 the parameters $\alpha$ determines the convergence speed in the real and Fourier space series. For a given choice of $\alpha$ we have chosen a real-space cutoff distance $r_c$ and a $k_c$ cutoff in the Fourier space. The cutoff $k_c$ determines the total number of Fourier components, $(4 \pi/3) n_c^3$ , where $n_c$ is a positive integer. This parameter has been choosen in such a way that the error on the Ewald summation is much smaller than the Quantum Monte Carlo statistical one.
A careful choice of the parameter $\alpha$ can minimize the error in the summation (see Ref. (80)). In our simulation we have chosen $\alpha = L/5$ , where $L$ is the size of the simulation box. With this cutoff it is sufficient to sum the short range part in the eq. 4.16 only on the first image of each particle. Notice that during each VMC or DMC simulation the ionic coordinates are fixed throughout the calculation. Therefore the contribution of the ion-ion Coulomb interaction in the short-range part can be evaluated only at the beginning of the simulation. As an electron $i$ is moved during a VMC calculation the sum of the short range part of the eq. 4.16 is easily updated subtracting the old contribution electron-electron and electron-ion due to the electron $i$ , and adding the new one.
The sum in Fourier space can be written as:

$$U_k = \sum_{ \vec k \neq 0} { \frac{4 \pi }{V \vert\vec k\vert^2}... ...r_i) \right)^2 + \left ( \sum_i^N \cos ( \vec k \vec r_i) \right)^2 \right ] },$$ (4.17)

then for each $\vec k$ vector all $sin$ and $cos$ are stored in such a way that when an electron moves, the sum can be easily updated without calculating all the elements from scratch.
It is easy to understand that the Ewald summation scales as $O(N^2)$ . In fact the updating of the eq. 4.16 costs $N$ times the number of Fourier's components. Then the number of Fourier component goes as $(1/\alpha)^3$ where $\alpha$ is proportional $L$ and for a given density scales as $N$ . The Ewald sums are faster than the QMC sweep and so, even if nowadays other faster techniques exist, as for instance particle-mesh-based one, it was not necessary to adopt other more complicated methods in our calculation.