# Fortran code for Ewald summation

I am trying to write a code to calculate the potential and forces, for the same using ewald summation.For this purpose, the formula for potential and force I have used is :

$$U = U^{(r)} + U^{(k)} + U^{(bc)} + U^{self}$$

where the k-space contribution of potential is given by $$U^{(k)} = \frac{1}{2\pi L^{3}}\sum_{\textbf{k}\ne0} \frac{4\pi^2}{k^2}\text {e}^{-\frac{k^2}{4\kappa^2}}|S(\textbf k)|^2 \qquad S(\textbf k) = \sum_{i=1}^N z_i\text{e}^{i\textbf{k}.\text{r}_i}$$ $$U^{(r)} = \sum_{j<i} \sum_{\textbf n=0}^{\infty}z_iz_j\frac{\text{erfc}(\kappa|\textbf{r}_{ij}+\textbf n|)}{|\textbf{r}_{ij}+\textbf n|}$$ $$U^{(bc)} = \frac{2\pi}{3L^3}\bigg|\sum_{i=1}^Nz_i\textbf r_i \bigg|^2 \qquad U^{(self)} = \frac{\kappa}{\sqrt{\pi}}\sum_{i=1}^Nz_i^2$$

and the force equations - $$\textbf F_i = \textbf F_i^{(k)} +\textbf F_i^{(r)}+\textbf F_i^{(bc)}$$

$$\textbf F_i^{(k)} = \frac{4\pi z_i}{2\pi L^3}\sum_{\textbf k\ne 0}\frac{\textbf k \text e^{\frac{-k^2}{4\kappa^2}}}{k^2}\bigg(\sin(\textbf k_i. \textbf r_i) \text {Re}(S(\textbf k)) + \cos(\textbf k_i. \textbf r_i) \text {Im}(S(\textbf k))\bigg)$$ and two more equations, which am tired of writing but kind of sure that they are correct !

I am currently using this method to solve for water molecules molecular dynamics simulation. The problem I am facing is the calculation of potential is highly dependent on $\kappa$, which it can be, but there is vast change in potential if I change $\kappa$ slightly. Secondly shouldn't the change with respect diminish beyond certain value, I don't see that also happening. Here is a segment of code used in k-space summation, I have my doubts on this part only :

!calculates the structure factor for each k-vector
do i=0,kreq
call struct_fact(kx(i),ky(i),kz(i),rx,ry,rz,q,nm,ns,sk_r(i),sk_i(i),cs,sn)
end do

!summing over all k-vectors
k_sum = 0.d0
const = 1/(4.d0*kappa*kappa)
do i=0,kreq
k_sum = k_sum + exp(-ksq(i)*const)*( sk_r(i)**2 + sk_i(i)**2 )/ksq(i)
end do

!calculation of force-vector for all sites in k-space
ew_kfsx = 0.d0
ew_kfsy = 0.d0
ew_kfsz = 0.d0
do k=0,kreq
do i=1,nm
do j=1,ns
temp = sn(i,j)*sk_r(k) + cs(i,j)*sk_i(k)
!print *, ( q(i,j)*temp*kx(i)*exp(-ksq(i)*const) )/ksq(i)
ew_kfsx(i,j) = ( q(i,j)*temp*kx(k)*exp(-ksq(k)*const) )/ksq(k)
ew_kfsy(i,j) = ( q(i,j)*temp*ky(k)*exp(-ksq(k)*const) )/ksq(k)
ew_kfsz(i,j) = ( q(i,j)*temp*kz(k)*exp(-ksq(k)*const) )/ksq(k)
end do
end do
end do

k_sum = (twopi/boxL**3)*k_sum
ew_kfsx = (twopi/boxL**3) * ew_kfsx
ew_kfsy = (twopi/boxL**3) * ew_kfsy
ew_kfsz = (twopi/boxL**3) * ew_kfsz


I here wish to know, if there is a precise way to find out if my code based on the above formulas are working correctly ?

Notice that the summation in $U^{(r)}$ is incorrect. You want to sum over all the copies of the atoms in the lattice of periodic boxes, not just those whose indices satisfy $i > j$. In the original box, of course you want to avoid self interactions $i = j$ but just in that one box. In other words, do not discard the electrostatic interactions between atom $i$ in the original box with atom $i$ in the neighboring periodic boxes. If you are using the minimum image convention and the value of $\kappa$ is right (I have not checked it, though, you can do so yourself simply by plotting ${\rm erfc}(\kappa \, r)$ for $r \ge 0$ and seeing if it decays fast enough), then you are already summing correctly.

The summation in the $U^{(k)}$ term goes over all the triplets of integers (positive and negative) except the one that is zero in all components (i.e., $\mathbf{k} \in \mathbb{Z}^3 \setminus \{ (0, 0, 0) \}$). It is not clear from your Fortran snippet that this is being taken into account.

There is a test that is computationally expensive but will make sure your Ewald summation code is working properly. Just compute the electrostatic energy,

$$E=\frac{1}{2}\sum_{i=1}^{N}\sum_{\stackrel{j=1}{i\neq j}}^{N} \frac{z_i z_j}{|\mathbf{r}_{ij} |}+\frac{1}{2}\sum_{\stackrel{\mathbf{n}\in\mathbb{Z}^{3% }}{\mathbf{n}\neq{0}}}\sum_{i=1}^{N}\sum_{j=1}^{N} \frac{z_i z_j}{|\mathbf{r}_{ij}+\mathbf{n}|},$$

over many cubic shells (as many as you need for $E$ to converge up to a desired accuracy), as specified here. The result of your Ewald summation must converge to $E$.

Changing $\kappa$ determines how much the work is distributed between the direct, $U^{(r)}$, and the Fourier space, $U^{(k)}$, parts of the computation. In principle, you should get the same value for the potential energy even if the values of $\kappa$ change.

I wrote a tutorial on this topic with an accompanying implementation of the Ewald summation and Particle-Mesh Ewald methods for dispersive (not Coulomb) potentials that you may find helpful. You can find it here.

• But following Allen and Tildesley, I used minimum image condition in the real space space summation, that is to say only $\textbf n = 0$ and with it used the value of $\kappa =5/L$. I thought since it was found in the book of Allen and Tildesley, it was a standard one. Sep 27 '15 at 3:15
• By many cubic shells, can you be specific about the number !! Sep 27 '15 at 3:18
• Using the minimum image convention is indeed equivalent for certain ranges of $\kappa$ to what is written in Allen and Tildesley but not to what was written in your question (notice the apostrophe in their summation symbol, etc.). Sep 27 '15 at 3:23
• I've added a remark following the "many cubic shells" sentence to clarify it. Sep 27 '15 at 3:24
• I am sorry, I am at loss to understand the ranges of $\kappa$ you're referring to, since I am using the one given from the same book itself. As for the apostraphe, doesn't it mean to say that its not including the i=j term in the $\textbf n=0$ term. Sep 27 '15 at 3:27

The National Institute of Standards and Technology (NIST) proposes different inputs and outputs to test your implementation. Have a look there:

https://www.nist.gov/mml/csd/chemical-informatics-research-group/spce-water-reference-calculations-10%C3%A5-cutoff

See paragraph 6. Reference Calculations of Intermolecular Energy for SPC/E that includes the usual contributions $E_{real}$, $E_{fourier}$ etc. with everything detailed (cutoffs, numerical precision, etc).

I was having a similar problem as the one you had.

The k-space, real part, and self-interaction are fairly dependent on $\kappa$, once you sum them up, the effects of $\kappa$ will cancel out.

You can check my (Python) implementation at this link: https://github.com/chemlab/chemlab/blob/180287718380ae635aa6401cbd183d93e3bdc81c/chemlab/md/ewald.py