xi3: Tuning Program for the Ewald Summation

xi3 in RYUON-stokes package is a tuning program for ewald summation parameter $\xi$.

For simplicity, here we show the equation for F version.

$$6\pi\mu a \bm{U}^\alpha = \bm{F}^\alpha + \sum_{\gamma} {\sum_{\beta}}' \bm{M}(\bm{r}_{\beta\alpha} + \bm{r}_\gamma) \cdot\bm{F}^\beta$$ (2.20)

$$= \bm{F}^\alpha$$ (2.21)

$$+ \sum_{\gamma} {\sum_{\beta\neq\alpha}}' \bm{M}^{\rm real}(\bm{r}_{\beta\alpha} + \bm{r}_\gamma) \cdot\bm{F}^\beta$$ (2.22)

$$+ \frac{1}{V} \sum_{\lambda} \sum_{\beta = 1}^{N} \cos(\bm{k}_\la... ...}_{\beta\alpha})\ \tilde{\bm{M}}^{\rm recip}(\bm{k}_\lambda) \cdot\bm{F}^\beta$$ (2.23)

$$- \bm{M}^{\rm recip}(\bm{r} = \bm{0}) \cdot\bm{F}^\alpha ,$$ (2.24)

where $\gamma$ is a index of lattice vector $\bm{r}_\gamma = (n_x l_x, n_y l_y, n_z l_z)^\dagger$ for arbitrary integers $(n_x, n_y, n_z)$, $\bm{r}_{\beta\alpha} = \bm{x}^\alpha - \bm{x}^\beta$. The prime on the summation of particle $\beta$ in real-space lattice sum means that the self part ( $\beta = \alpha$) is excluded for the primary cell $\bm{r}_{\gamma_0} = \bm{0}$. The mobility matrix $\bm{M}(\bm{r})$ is so called Rotne-Prager-Yamakawa tensor given by

$$\bm{M}(\bm{r}) = \frac{3a}{4} \left( 1 + \frac{a^2\nabla^2}{6} \right)^2 \bm{J}(\bm{r})$$ (2.25)

$$= \frac{3}{4} \left[ \left( \frac{a}{r} + \frac{a^3}{3r^3} \right) \bm{I} + \left( \frac{a}{r} - \frac{a^3}{r^3} \right) \frac{\bm{rr}}{r^2} \right]$$ (2.26)

$$= \left( \frac{3a}{4} + \frac{a^3}{4} \nabla^2 \right) \left( \bm{I} \nabla^2 - \bm{\nabla\nabla} \right) r .$$ (2.27)

Because it is long-range interaction proportional to $1/r$, we need to take huge number of periodic images. The trick of Ewald summation technique is that using the identity

$${\rm erfc}(x) + {\rm erf}(x) \equiv 1,$$ (2.28)

we split $\bm{M}(\bm{r})$ into two parts, $\bm{M}^{\rm real}$ and $\bm{M}^{\rm recip}$ as

$$\bm{M}^{\rm real} = \left( \frac{3a}{4} + \frac{a^3}{4} \nabla^2 \right) \left( \bm{I} \nabla^2 - \bm{\nabla\nabla} \right) r\ {\rm erfc} (\xi r),$$ (2.29)

$$\bm{M}^{\rm recip} = \left( \frac{3a}{4} + \frac{a^3}{4} \nabla^2 \right) \left( \bm{I} \nabla^2 - \bm{\nabla\nabla} \right) r\ {\rm erf} (\xi r) .$$ (2.30)

Because the error function ${\rm erfc}(x)$ decays exponentially, we can truncate the lattice summation of $\gamma$ for $\bm{M}^{\rm real}$ at some point. Similarly, in reciprocal space, $\tilde{\bm{M}}^{\rm recip}(\bm{k})$ decays exponentially in $k$ so that we can truncate the lattice sum of $\lambda$ at some point.

Mathematically, the parameter $\xi$ introduced in Eqs. (2.29) and (2.30) is arbitrary. It specify how fast $\bm{M}^{\rm real}(\bm{r})$ decays and how slow $\tilde{\bm{M}}^{\rm recip}(\bm{k})$ decays. Practically, on the other hand, we can choose the optimal value of $\xi$ so that the calculation time is minimal.