# Calculation of long-range forces in Molecular Dynamics - 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.

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 ?

EDIT 1: As a note, I have simulated, SPC/E water model using this now. The problem that I am facing right now is potential energy that is being calculated is 5 to 6 times more than the ones reported earlier in journals, but however the energy conservation during the run is intact. I don't understand the reason here !!

I've only ever used the Ewald sum, I've never implemented it myself.

However, you mention that you're not converging as $\kappa$ increases nor are you converging to the correct value. It would seem that regardless of the problem, if your implementation is correct it should converge at some point.

If you do reach convergence wrt $\kappa$; as to the point that you are not getting the reported values. Remember that the size of the simulation box is also important. Are you sure you are using a large enough simulation space with enough water molecules in it?

EDIT:

Well energy should be conserved, :P. If your force field, whatever it is you implemented, is conservative then energy should be conserved (or your motion integrator is wrong).

Which ensemble are you using? The paper? If NVE, are you sure the V is exactly right? remember even 0.1% change in box cell size will cause a huge change in energy for a condensed material. (try compressing water 0.1%!)

Are you simulating fluid water or crystalline water? Either way, are you sure you're not stuck in a local minima? Try pulsating the sim box's volume to get out.

• I am using a simulation box size of 18.63 angstroms with 216 molecules !! I have used the particulars given from an old paper, truncating the ewald sum by the minimum image condition !! How can I possibly debug my code ? Why is the energy conservation still intact ? Commented Sep 19, 2015 at 4:07
• +1 Thank you for the insights. I am using NVE ensemble only, the paper I am referring to is one by M. Prevost et al.. One part I didn't understand about the paper was, he has mentioned using 342 vectors in reciprocal space with a magnitude of 3 in the sum over reciprocal lattice sum. Now I didn't understand which magnitude he's referring to, so I couldn't strictly adhere to the magnitude of the k-vectors, although I kept it below an upper limit !! Commented Sep 22, 2015 at 4:21