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

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 !!

Answer 1

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:

Based on your comments below:

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.