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 ?