Force Error

The Force error is the error incurred when we cut the potential interaction after a certain distance. Following the works of [Dharuman et al., 2017, Kolafa and Perram, 1992, Stern and Calkins, 2008] we define the total force error for our P3M algorithm as

$$\mathrm{\Delta }{F}_{\text{tot}}=\sqrt{\mathrm{\Delta }{F}_{\mathcal{R}}^2+\mathrm{\Delta }{F}_{\mathcal{F}}^2}$$

where $\mathrm{\Delta }{F}_{\mathcal{R}}$ is the error obtained in the PP part of the force calculation and $\mathrm{\Delta }{F}_{\mathcal{F}}$ is the error obtained in the PM part, the subscripts $\mathcal{R}\mathcal{,}\mathcal{F}$ stand for real space and Fourier space respectively. $\mathrm{\Delta }{F}_{\mathcal{R}}$ is calculated as follows

$$\mathrm{\Delta }{F}_{\mathcal{R}}=\sqrt{\frac{N}{V}}{\left[{\int }_{r_c}^{\mathrm{\infty }}{d}^{3}r{|\mathrm{\nabla }{\varphi }_{\mathcal{R}}(\mathbf{r})|}^2\right]}^{1/2},$$

where ${\varphi }_{\mathcal{R}}(\mathbf{r})$ is the short-range part of the chosen potential. In our example case of a Yukawa potential we have

$${\varphi }_{\mathcal{R}}(r)=\frac{Q^2}{2r}[e^{-\kappa r}\text{erfc}(\alpha r-\frac{\kappa }{2\alpha })+e^{\kappa r}\text{erfc}(\alpha r+\frac{\kappa }{2\alpha })],$$

where $\kappa ,\alpha$ are the dimensionless screening parameter and Ewald parameter respectively and, for the sake of clarity, we have a charge $Q=Ze/\sqrt{4\pi {ϵ}_0}$ with an ionization state of $Z=1$. Integrating this potential, and neglecting fast decaying terms, we find

$$\mathrm{\Delta }{F}_{\mathcal{R}}\simeq 2Q^2\sqrt{\frac{N}{V}}\frac{e^{-{\alpha }^2{r}_c^2}}{\sqrt{r_c}}{e}^{-{\kappa }^2/4{\alpha }^2}.$$

On the other hand $\mathrm{\Delta }{F}_{\mathcal{F}}$ is calculated from the following formulas

$$\mathrm{\Delta }{F}_{\mathcal{F}}=\sqrt{\frac{N}{V}}\frac{Q^2\chi }{\sqrt{V^{1/3}}}$$

$${\chi }^2{V}^{2/3}=\left(\sum _{\mathbf{k}\ne 0}{G}_{\mathbf{k}}^2|\mathbf{k}{|}^2\right)-\sum _{\mathbf{n}}\left[\frac{{(\sum _{\mathbf{m}}{\hat{U}}_{\mathbf{k}\mathbf{+}\mathbf{m}}^2{G}_{\mathbf{k}\mathbf{+}\mathbf{m}}{\mathbf{k}}_{\mathbf{n}}\cdot {\mathbf{k}}_{\mathbf{n}\mathbf{+}\mathbf{m}})}^2}{{\left(\sum _{\mathbf{m}}{\hat{U}}_{{\mathbf{k}}_{\mathbf{n}\mathbf{+}\mathbf{m}}}^2\right)}^2|{\mathbf{k}}_{\mathbf{n}}{|}^2}\right].$$

This is a lot to take in, so let’s unpack it. The first term is the RMS of the force field in Fourier space obtained from solving Poisson’s equation $-\mathrm{\nabla }\varphi (\mathbf{r})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$ in Fourier space. In a raw Ewald algorithm this term would be the PM part of the force. However, the P3M variant solves Poisson’s equation on a Mesh, hence, the second term which is non other than the RMS of the force obtained on the mesh. $G_{\mathbf{k}}$ is the optimal Green’s function which for the Yukawa potential is

$$G_{\mathbf{k}}=\frac{4\pi e^{-({\kappa }^2+{\left|\mathbf{k}\right|}^2)/(4{\alpha }^2)}}{{\kappa }^2+|\mathbf{k}{|}^2}$$

where

$$\mathbf{k}(n_x,n_y,n_z)={\mathbf{k}}_{\mathbf{n}}=(\frac{2\pi n_x}{L_x},\frac{2\pi n_y}{L_y},\frac{2\pi n_z}{L_z}).$$

${\hat{U}}_{\mathbf{k}}$ is the Fourier transform of the B-spline of order $p$

$${\hat{U}}_{{\mathbf{k}}_{\mathbf{n}}}={\left[\frac{\sin (\pi n_x/M_x)}{\pi n_x/M_x}\right]}^p{\left[\frac{\sin (\pi n_y/M_y)}{\pi n_y/M_y}\right]}^p{\left[\frac{\sin (\pi n_z/M_z)}{\pi n_z/M_z}\right]}^p,$$

where $M_{x,y,z}$ is the number of mesh points along each direction. Finally the $\mathbf{m}$ refers to the triplet of grid indices $(m_x,m_y,m_z)$ that contribute to aliasing. Note that in the above equations as $\kappa \to 0$ (Coulomb limit), we recover the corresponding error estimate for the Coulomb potential.

The reason for this discussion is that by inverting the above equations we can find optimal parameters $r_c,\alpha$ given some desired errors $\mathrm{\Delta }{F}_{\mathcal{R}\mathcal{,}\mathcal{F}}$. While the equation for $\mathrm{\Delta }{F}_{\mathcal{R}}$ can be easily inverted for $r_c$, such task seems impossible for $\mathrm{\Delta }{F}_{\mathcal{F}}$ without having to calculate a Green’s function for each chosen $\alpha$. As you can see in the second part of the output the time it takes to calculate $G_{\mathbf{k}}$ is in the order of seconds, thus, a loop over several $\alpha$ values would be very time consuming. Fortunately researchers have calculated an analytical approximation allowing for the exploration of the whole $r_c,\alpha$ parameter space [Dharuman et al., 2017]. The equations of this approximation are

$$\mathrm{\Delta }{F}_{\mathcal{F}}^{(\text{approx})}\simeq Q^2\sqrt{\frac{N}{V}}{A}_{\mathcal{F}}^{1/2},$$

$$A_{\mathcal{F}}\simeq \frac{3}{2{\pi }^2}\sum _{m=0}^{p-1}{C}_m^{(p)}{\left(\frac{h}{2}\right)}^{2(p+m)}\frac{2}{1+2(p+m)}\beta (p,m),$$

$$\beta (p,m)={\int }_0^{\mathrm{\infty }}dk{G}_k^2{k}^{2(p+m+2)},$$

where $h=L_x/M_x$ and the coefficients $C_m^{(p)}$ are listed in Table I of [Deserno and Holm, 1998].

Finally, by calculating

$$\mathrm{\Delta }{F}_{\text{tot}}^{(\text{apprx})}(r_c,\alpha )=\sqrt{\mathrm{\Delta }{F}_{\mathcal{R}}^2+(\mathrm{\Delta }{F}_{\mathcal{F}}^{(\text{approx})}{)}^2}$$

we are able to investigate which parameters $r_c,\alpha$ are optimal for our simulation.