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