Fortran code for Ewald summation

Question

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 ?

Answer 1

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.

Answer 2

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).

Answer 3

I don't think your formula is correct for the force. You can rewrite your energy expression in an easier form, and it makes taking the derivative easier as well.

The way you have it is

\begin{equation} U_{\rm recip} = \frac{1}{2L^3} \sum_{\vec \kappa\ne0} \frac{4 \pi}{\kappa^2} |\rho(\vec\kappa)|^2 \exp^{(-\kappa^2 / 2\alpha)} \end{equation} However, in practice, it is easier to wield if written as \begin{equation} U_{\rm recip} = \frac{1}{2L^3} \sum_{\vec \kappa\ne0} \sum_{i=1} \sum_{j=1} \frac{4 \pi q_i q_j}{\kappa^2} \exp^{[i\vec \kappa \cdot (\vec r_j - \vec r_i)]} \exp^{(-\kappa^2 / 2\alpha)} \end{equation} It is now a lot clearer how to take the gradient with respect to $\vec r_i$

Note that $\frac{d \exp^{[i\vec \kappa \cdot (\vec r_j - \vec r_i)]}}{d \vec r_i} = -i \vec \kappa \exp^{[i\vec \kappa \cdot (\vec r_j - \vec r_i)]} = -\kappa \sin[\vec \kappa \cdot (\vec r_j - \vec r_i)]$

See this vector/matrix cookbook for help with the derivatives math cookbook

So the only part of the reciprocal energy that depends on position $\vec r_i$ is the exponential, which I have shown how to take the derivative of, so we get for $\vec F_i = -\nabla(U)$

\begin{equation} \vec F_i = q_i\frac{4\pi}{L^3} \sum_{j=0} q_j \sum_{\vec \kappa\ne0} \frac{\vec \kappa}{\kappa^2} \exp^{(-\kappa^2/{4\kappa})}\sin[\vec \kappa \cdot \vec r_{ij}] \end{equation}

In some respects, the real contribution is more work since the derivative involves the chain rule, but that also is minor.

You can find more information on the energy terms in Frenkel & Smit and the final answer can also be found (no maths shown) here as well Forces of Ewald.

Answer 4

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