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

\[\Delta F_{\textrm{tot}} = \sqrt{ \Delta F_{\mathcal R}^2 + \Delta F_{\mathcal F}^2 }\]

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

\[\Delta F_{\mathcal R} = \sqrt{\frac{N}{V} } \left [ \int_{r_c}^{\infty} d^3r \left | \nabla \phi_{\mathcal R}( \mathbf r) \right |^2 \right ]^{1/2},\]

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

\[\phi_{\mathcal R}(r) = \frac{Q^2}{2r} \left [ e^{- \kappa r} \text{erfc} \left( \alpha r - \frac{\kappa}{2\alpha} \right ) + e^{\kappa r} \text{erfc} \left( \alpha r + \frac{\kappa}{2\alpha} \right ) \right ],\]

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 \epsilon_0}\) with an ionization state of \(Z = 1\). Integrating this potential, and neglecting fast decaying terms, we find

\[\Delta F_{\mathcal R} \simeq 2 Q^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 \(\Delta F_{\mathcal F}\) is calculated from the following formulas

\[\Delta F_{\mathcal F} = \sqrt{\frac{N}{V}} \frac{Q^2 \chi}{\sqrt{V^{1/3}}}\]

\[\chi^2V^{2/3} = \left ( \sum_{\mathbf k \neq 0} G_{\mathbf k}^2 |\mathbf k |^2 \right ) - \sum_{\mathbf n} \left [ \frac{\left ( \sum_{\mathbf m} \hat{U}_{\mathbf{k + m}}^2 G_{\mathbf{k+m}} \mathbf{k_n} \cdot \mathbf{k_{n + m}} \right )^2 }{ \left( \sum_{\mathbf m} \hat{U}_{\mathbf{k_{n+m}}}^2 \right )^2 |\mathbf{k_{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 \(-\nabla \phi(\mathbf r) = \delta( \mathbf r - \mathbf r')\) 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_n} = \left ( \frac{2 \pi n_x}{L_x}, \frac{2 \pi n_y}{L_y}, \frac{2 \pi n_z}{L_z} \right ).\]

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

\[\hat U_{\mathbf{k_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 \rightarrow 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 \(\Delta F_{\mathcal {R,F}}\). While the equation for \(\Delta F_{\mathcal R}\) can be easily inverted for \(r_c\), such task seems impossible for \(\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

\[\Delta F_{\mathcal F}^{(\textrm{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^{\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

\[\Delta F_{\textrm{tot}}^{(\textrm{apprx})}( r_c, \alpha) = \sqrt{ \Delta F_{\mathcal R}^2 + ( \Delta F_{\mathcal F}^{(\textrm{approx})} ) ^2 }\]

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