Application of Ewald summation to molecular electrostatic potential energy calculations

Technical background

Calculating the electrostatic potential for periodic boundary conditions is often seen as one of the bottlenecks in molecular dynamics simulations. Formally it is easier, for example the electrostatic potential energy is if periodic boundary conditions are not considered:

\[E=\frac{1}{4\pi\epsilon_0}\sum_{i=0}^{N-2}\sum_{j=i+1}^{N-1}\frac{q_iq_j}{r_{ij}} \]

If periodic boundary conditions are considered, then the electrostatic potential energy becomes:

\[E=\frac{1}{4\pi\epsilon_0}\sum_{\mathbf{n}}\sum_{i=0}^{N-2}\sum_{j=i+1}^{N-1}\frac{q_iq_j}{\left|\mathbf{r}_{ij}+\mathbf{n}\mathbf{L}\right|}+\frac{1}{4\pi\epsilon_0}\sum_{|\mathbf{n}|>0}\frac{q_i^2}{\left|\mathbf{n}\mathbf{L}\right|},\mathbf{n}=(n_1,n_2,...,n_d)\in\mathbb{Z} \]

Not to mention the superposition of an infinite number of boxes, even a single\(O(N^2)\)of computational complexity, which is the reason for the difficulty in solving. In molecular dynamics, three ideas are used to simplify this computational effort: the Fourier Transform, Ewald Summation, and the Particle Mesh approach, and this paper deals with the part of the Fourier Transform that is computed with Ewald Summation.

periodic potential

First we understand the meaning of the electrostatic potential energy term in terms of physical concepts: a single charge\(q_i\)would create an electric potential in space\(V_i\)When another charge\(q_j\)be situated at\(V_i\)When in the corresponding electric field, one is subjected to electrostatic interactions with an energy of\(E_{ij}\). Since the electrostatic potential energy at infinity is 0, this also assumes that if we need to convert the\(q_j\)Pushing beyond infinity from the original position would require the\(E_{ij}\)of energy. This idea tells us that we can first calculate the potential\(V_i\)and then calculate the potential energy\(E_{ij}\), which is the same thing:

\[V_i(\mathbf{r})=\frac{1}{4\pi\epsilon_0}\frac{q_i}{\left|\mathbf{r}-\mathbf{r}_i\right|} \]

If one considers the upper periodic boundary condition that is:

\[V_i(\mathbf{r})=\frac{1}{4\pi\epsilon_0}\sum_{\mathbf{n}}\frac{q_i}{\left|\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L}\right|} \]

Because there is no way to calculate this infinitely many summation directly, it is only clear that the potential is periodic:\(V_i(\mathbf{r})=V_i(\mathbf{r}+\mathbf{n}\mathbf{L})\)。

Poisson's equation for an electrostatic field

Charge in space\(i\)exist\(\mathbf{r}\)The potential at is given by the Poisson equation:

\[\Delta V_i(\mathbf{r})=-\frac{\rho_i(\mathbf{r})}{\epsilon_0} \]

included among these\(\Delta\)is the Laplace operator.\(\rho_i\)denotes the charge density. If the case of a point charge is considered, then there is:\(\rho_i(\mathbf{r})=q_i\delta(\mathbf{r}-\mathbf{r}_i)\). In turn, write the Poisson equation in Euclidean space:

\[\nabla^2V_i(\mathbf{r})=-\frac{q_i\delta(\mathbf{r}-\mathbf{r}_i)}{\epsilon_0} \]

included among these\(\nabla^2V_i(\mathbf{r})=\frac{\partial^2V_i(\mathbf{r})}{\partial x^2}+\frac{\partial^2V_i(\mathbf{r})}{\partial y^2}+\frac{\partial^2V_i(\mathbf{r})}{\partial z^2}\). Consider again the periodic boundary conditions, with a point charge in each box\(q_i\), so the equation should be written as:

\[\nabla^2V_i(\mathbf{r})=-\sum_{\mathbf{n}}\frac{q_i\delta(\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L})}{\epsilon_0} \]

Poisson equation in Fourier space

Simultaneous Fourier transforms are applied to both sides of the Poisson equation in the above single-point form (for an understanding of Fourier transforms, seePreface Article 1respond in singingPreface Article 2, there is a more detailed introduction to the principle and related code implementation) there:

\[\int\left(\frac{\partial^2V_i(\mathbf{r})}{\partial x^2}+\frac{\partial^2V_i(\mathbf{r})}{\partial y^2}+\frac{\partial^2V_i(\mathbf{r})}{\partial z^2}\right)e^{-2\pi j\mathbf{k}\mathbf{r}}d\mathbf{r}=-\frac{q_i}{\epsilon_0}\int\delta(\mathbf{r}-\mathbf{r}_i)e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r} \]

Start by calculating one of the items on the left using the division integral:

\[\begin{align*} \int\frac{\partial^2V_i(\mathbf{r})}{\partial x^2}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}&=\left.\frac{\partial V_i(\mathbf{r})}{\partial x}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}\right|_{-\infty}^{\infty}+2\pi jk_x\int\frac{\partial V_i(\mathbf{r})}{\partial x}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r} \end{align*} \]

It is important to note that here\(\frac{\partial V_i(\mathbf{r})}{\partial x}=F_i(\mathbf{r})\)the physical meaning of the action, then we can take the first type of boundary conditions:

\[\lim_{x\rightarrow\infty}V_i(x)=0 \]

Such that by definition of differentiation there is:

\[\lim_{x\rightarrow\infty}\frac{\partial V_i(x)}{\partial x}=\lim_{x\rightarrow\infty}\lim_{\epsilon\rightarrow0}\frac{V_i(x+\epsilon)-V_i(x)}{\epsilon}=0 \]

This limit above substitutes for the\(V_i(\mathbf{r})\)The significance of the single-point form, which is located in the\(\mathbf{r}_i\)point charge at\(q_i\), the electric field generated at infinity is 0 for the charge at infinity\(q_{\infty}\)The force is also 0. The periodic boundary conditions are not considered here. Then there are:

\[\begin{align*} \int\frac{\partial^2V_i(\mathbf{r})}{\partial x^2}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}&=2\pi jk_x\int\frac{\partial V_i(\mathbf{r})}{\partial x}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}\\ &=2\pi jk_x\left[\left.V_i(\mathbf{r})e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}\right|_{x=-\infty}^{x=\infty}+2\pi jk_x\int V_i(\mathbf{r})e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}\right]\\ &=-4\pi^2k_x^2V_i(\mathbf{k}) \end{align*} \]

Similarly, the left-hand side of the Poisson equation can be obtained in the form:

\[\int\nabla^2V_i(\mathbf{r})e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}=-4\pi^2k^2V_i(\mathbf{k}) \]

And the right-hand side form requires the sampling properties of the Dirac function:

\[\int_{-\infty}^{\infty}\delta(t)f(t)dt=\int_{-\infty}^{\infty}\delta(t)f(0)dt=f(0)\int_{-\infty}^{\infty}\delta(t)dt=f(0) \]

To wit:

\[\int\delta(\mathbf{r}-\mathbf{r}_i)e^{-2\pi j\mathbf{k}\cdot\mathbf{r}}d\mathbf{r}=e^{-2\pi j\mathbf{k}\cdot\mathbf{r}_i} \]

The single-point Poisson equation in Fourier space is obtained:

\[4\pi^2k^2V_i(\mathbf{k})=\frac{q_i}{\epsilon_0}e^{-2\pi j\mathbf{k}\cdot\mathbf{r}_i} \]

reciprocal space

When it comes to Fourier space, it is natural to think of using the inverse eigenspace transformation of solid state physics, that is, treating a periodic box as a protocell. According to the inverse easy space (also called\(\mathbf{k}\)The (spatial) lattice vector (inverted lattice vector) definition has:

\[\mathbf{k}=m_1\mathbf{k}_1+m_2\mathbf{k}_2+m_3\mathbf{k}_3\\ \]

Among them:

\[k_i\cdot L_j=\delta_{i,j}\\ \delta_{i,j}=\left\{ \begin{matrix} 1, i=j\\ 0, i\neq j \end{matrix} \right. \]

Follow our usual rectangular periodic boundary conditions:

\[\begin{matrix} \mathbf{L}=(L_x,L_y,L_z)=\mathbf{L}_1+\mathbf{L}_2+\mathbf{L}_3\\ \mathbf{L}_1 = (L_x,0,0)\\ \mathbf{L}_2 = (0,L_y,0)\\ \mathbf{L}_3 = (0,0,L_z) \end{matrix} \]

It can be calculated:

\[\mathbf{k}_1=\frac{2\pi(\mathbf{L}_2\times\mathbf{L}_3)}{\mathbf{L}_1\cdot(\mathbf{L}_2\times\mathbf{L}_3)}= \frac{2\pi}{\Omega}(\mathbf{L}_2\times\mathbf{L}_3)=\frac{\pi L_yL_z}{\Omega L_x}\mathbf{L}_1=\frac{2\pi}{L_x^2}\mathbf{L}_1 \]

included among these\(\Omega\)denotes the volume of the periodic box, similarly:

\[\mathbf{k}_2=\frac{2\pi(\mathbf{L}_3\times\mathbf{L}_1)}{\mathbf{L}_2\cdot(\mathbf{L}_3\times\mathbf{L}_1)}= \frac{2\pi}{L_y^2}\mathbf{L}_2 \]

as well as

\[\mathbf{k}_3=\frac{2\pi(\mathbf{L}_1\times\mathbf{L}_2)}{\mathbf{L}_3\cdot(\mathbf{L}_1\times\mathbf{L}_2)}= \frac{2\pi}{L_z^2}\mathbf{L}_3 \]

After the inverse easy space transformation, the protocell volume increases from\(\Omega=L_xL_yL_z\)Becoming:\(\Omega^*=\frac{(2\pi)^3}{L_xL_yL_z}\). Because in the previous step of the Poisson equation in Fourier space we note that\(\mathbf{k}\)The front always carries a\(2\pi\), it may be useful here to use the definition of the inverse easy lattice vector for the\(\mathbf{k}\)to do a simplification of the form of the

\[\mathbf{k}=\left(\frac{2\pi}{L_x}, \frac{2\pi}{L_y}, \frac{2\pi}{L_z}\right) \]

In this way the Poisson equation in Fourier space can be abbreviated as:

\[V_i(\mathbf{k})=\frac{q_i}{\epsilon_0k^2}e^{-j\mathbf{k}\cdot\mathbf{r}_i} \]

included among these\(k^2=|\mathbf{k}|^2=\frac{4\pi^2m_1^2}{L_x^2}+\frac{4\pi^2m_2^2}{L_y^2}+\frac{4\pi^2m_3^2}{L_z^2}\)The transformation from real space to inverse easy space can be realized.

Decay function construction

For the form of the single-point potential after the Fourier transform described above, even if we have the entire\(\mathbf{k}\)space for integration, which is also a divergence result. So a very special idea is used here, proposed by Edwald, to divide the electrostatic energy term into a long-range interaction term and a short-range interaction term, which converge in inverse and real space, respectively, so that the electrostatic energy can be calculated exactly. In practice, there are different derivations, so we quote a more "mathematical" one (see link 1).

First we construct a decay function\(e^{-k^2t}\), a property of this decay function is:

\[\int_0^{\infty}e^{-k^2t}dt=\left.\left(-\frac{1}{k^2}e^{-k^2t}\right)\right|_0^{\infty}\\=\frac{1}{k^2} \]

This allows us to replace the Fourier-Poisson equation in the Fourier-Poisson equation with this form of integration\(\frac{1}{k^2}\)Item:

\[V_i(\mathbf{k})=\frac{q_i}{\epsilon_0}e^{-j\mathbf{k}\cdot\mathbf{r}_i}\int_0^{\infty}e^{-k^2t}dt \]

Since an integral form from 0 to infinity is used here, then we can implement a truncation that divides it into the sum of two integrals, if we are in the\(\eta\)Doing a truncation at the place, then there is:

\[V_i(\mathbf{k})=\frac{q_i}{\epsilon_0}e^{-j\mathbf{k}\cdot\mathbf{r}_i}\left(\int_0^{\eta}e^{-k^2t}dt+\int_{\eta}^{\infty}e^{-k^2t}dt\right) \]

The Short Term (ST) interaction is taken here:

\[V_i^S(\mathbf{k})=\frac{q_i}{\epsilon_0}e^{-j\mathbf{k}\cdot\mathbf{r}_i}\int_0^{\eta}e^{-k^2t}dt \]

as well as Long Term (LTP) interactions for:

\[V_i^L(\mathbf{k})=\frac{q_i}{\epsilon_0}e^{-j\mathbf{k}\cdot\mathbf{r}_i}\int_{\eta}^{\infty}e^{-k^2t}dt=\frac{q_i}{\epsilon_0k^2}e^{-j\mathbf{k}\cdot\mathbf{r}_i}e^{-\eta k^2} \]

Calculation of short-range action terms

As expected, after dividing the short-range action term and the long-range action term, a real-space convergent short-range interaction should be obtained, and we do an inverse Fourier transform of the short-range action:

\[\begin{align*} V_i^S(\mathbf{r})&=\frac{1}{k_xk_yk_z}\sum_{\mathbf{k}}V_i^S(\mathbf{k})e^{i\mathbf{k}\cdot\mathbf{r}}\\ &=\frac{q_i}{k_xk_yk_z\epsilon_0}\sum_{\mathbf{k}}e^{j\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}_i)}\int_0^{\eta}e^{-k^2t}dt\\ &=\frac{q_i}{k_xk_yk_z\epsilon_0}\int_0^{\eta}\sum_{\mathbf{k}}e^{j\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}_i)}e^{-k^2t}dt \end{align*} \]

It is clear that the summation term within the integral is a discrete Fourier transform of an exponential square function. And we can know that the normal distribution function\(f(\xi)=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{\xi^2}{2\sigma^2}}\)The Fourier and inverse Fourier transforms do not change the form of the distribution:

\[\begin{align*} F(k)&=\int f(\xi)e^{-jk\xi}d\xi\\ &=\frac{1}{\sqrt{2\pi}\sigma}\int e^{\frac{-\xi^2-2jk\xi\sigma^2}{2\sigma^2}}d\xi\\ &=\frac{1}{\sqrt{2\pi}\sigma}\int e^{\frac{-\xi^2-2jk\xi\sigma^2-(jk\sigma^2)^2+(jk\sigma^2)^2}{2\sigma^2}}d\xi\\ &=\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{(jk\sigma^2)^2}{2\sigma^2}}\int e^{-\frac{(\xi+jk\sigma^2)^2}{2\sigma^2}}d\xi \end{align*} \]

Since the integral term is just an integral in real space, it is still essentially an integral of a normal distribution function, and we know that the result of its integration is a constant, so there:

\[F(k)=e^{-\frac{k^2\sigma^2}{2}} \]

is also a normal distribution, except that its mean square deviation varies from\(\sigma\)It's turned into\(\frac{1}{\sigma}\), that is, its integral results in:

\[\int F(k)dk=\frac{\sqrt{2\pi}}{\sigma} \]

In the same way we can calculate that the inverse Fourier transform of a normal distribution function is still a normal distribution function with the mean square deviation of\(\frac{1}{\sigma}\)。
So back to the short-range static potential, let's start with a variable substitution\(t=\frac{\sigma^2}{2},\sigma\geq0\)then there are\(dt=\sigma d\sigma,\sigma_{t=\eta}=\sqrt{2\eta},\sigma_{t=0}=0\). At this point the short-range interaction potential is:

\[\begin{align*} V_i^S(\mathbf{r})&=\frac{q_i}{k_xk_yk_z\epsilon_0}\int_{0}^{\sqrt{2\eta}}\left(\sum_{\mathbf{k}}e^{j\mathbf{k}\cdot(\mathbf{r}-\mathbf{r}_i)}e^{-\frac{k^2\sigma^2}{2}}\right)\sigma d\sigma\\ &=\frac{q_i}{\epsilon_0}\int_{0}^{\sqrt{2\eta}}\frac{e^{-\frac{(\mathbf{r}-\mathbf{r}_i)^2}{2\sigma^2}}}{N_{coe}}\sigma d\sigma \end{align*} \]

here are\(N_{coe}\)is a constant used to normalize the normal distribution:

\[\begin{align*} N_{coe}&=\int_{\mathbf{r}}e^{-\frac{\left|\mathbf{r}-\mathbf{r}_i\right|}{2\sigma^2}}d\mathbf{r}\\ &=\int_{-\infty}^{\infty}e^{-\frac{(z-z_i)^2}{2\sigma^2}}\left[\int_{-\infty}^{\infty}e^{-\frac{(y-y_i)^2}{2\sigma^2}}\left(\int_{-\infty}^{\infty}e^{-\frac{(x-x_i)^2}{2\sigma^2}}dx\right)dy\right]dz\\ &=(2\pi\sigma^2)^{\frac{3}{2}} \end{align*} \]

So there:

\[\begin{align*} V_i^S(\mathbf{r})&=\frac{q_i}{\epsilon_0}\int_{0}^{\sqrt{2\eta}}\frac{e^{-\frac{(\mathbf{r}-\mathbf{r}_i)^2}{2\sigma^2}}}{(2\pi\sigma^2)^{\frac{3}{2}}}\sigma d\sigma\\ &=\frac{q_i}{(2\pi)^{\frac{3}{2}}\epsilon_0}\int_{0}^{\sqrt{2\eta}}\frac{e^{-\frac{(\mathbf{r}-\mathbf{r}_i)^2}{2\sigma^2}}}{\sigma^2} d\sigma \end{align*} \]

Replace it with a variable here:\(y=-\frac{|\mathbf{r}-\mathbf{r}_i|}{\sqrt{2}\sigma}\)Then there is:\(\sigma=-\frac{|\mathbf{r}-\mathbf{r}_i|}{\sqrt{2}y}\), whose differential transformations take the form:\(d\sigma=\frac{|\mathbf{r}-\mathbf{r}_i|}{\sqrt{2}y^2}dy\), whose boundaries are:\(y_{\sigma=\sqrt{2\eta}}=-\frac{|\mathbf{r}-\mathbf{r}_i|}{2\sqrt{\eta}},y_{\sigma=0}=-\infty\), substituting to get:

\[\begin{align*} V_i^S(\mathbf{r})&=\frac{q_i}{(2\pi)^{\frac{3}{2}}\epsilon_0}\int_{-\infty}^{-\frac{|\mathbf{r}-\mathbf{r}_i|}{2\sqrt{\eta}}}\frac{\sqrt{2}e^{-y^2}}{|\mathbf{r}-\mathbf{r}_i|} dy\\ &=\frac{q_i}{2\pi^{\frac{3}{2}}\epsilon_0|\mathbf{r}-\mathbf{r}_i|}\int_{-\infty}^{-\frac{|\mathbf{r}-\mathbf{r}_i|}{2\sqrt{\eta}}}e^{-y^2}dy\\ &=\frac{q_i}{2\pi^{\frac{3}{2}}\epsilon_0|\mathbf{r}-\mathbf{r}_i|}\int_{\frac{|\mathbf{r}-\mathbf{r}_i|}{2\sqrt{\eta}}}^{+\infty}e^{-y^2}dy \end{align*} \]

At this point an error function is used\(Erfc(y)=\frac{2}{\sqrt{\pi}}\int_{y}^{\infty}e^{-x^2}dx\)Substitute to make the substitution:

\[V_i^S(\mathbf{r})=\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i|}Erfc\left(\frac{|\mathbf{r}-\mathbf{r}_i|}{2\sqrt{\eta}}\right) \]

on account of\(\eta\)is a constant covariate that we introduce manually, and if we consider the\(\eta=\frac{\sigma^2}{2}\), then the form becomes:

\[V_i^S(\mathbf{r})=\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i|}Erfc\left(\frac{|\mathbf{r}-\mathbf{r}_i|}{\sqrt{2}\sigma}\right) \]

This form of short-range interaction potential represents a single charged particle inside a single box\(q_i\)exist\(\mathbf{r}\)The potential at, if periodic boundary conditions are considered, the form needs to be changed to:

\[V_i^S(\mathbf{r})=\sum_{\mathbf{n}}\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}Erfc\left(\frac{|\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}{\sqrt{2}\sigma}\right) \]

This form of interaction potential compared to the original form\(V_i(\mathbf{r})=\frac{1}{4\pi\epsilon_0}\sum_{\mathbf{n}}\frac{q_i}{\left|\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L}\right|}\)for example, an error function is used to make a truncation of the real space:

\[V_i^S(\mathbf{r})\approx\sum_{\mathbf{n'}}\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i+\mathbf{n'}\mathbf{L}|}Erfc\left(\frac{|\mathbf{r}-\mathbf{r}_i+\mathbf{n'}\mathbf{L}|}{\sqrt{2}\sigma}\right),1-Erfc\left(\frac{|\mathbf{r}-\mathbf{r}_i+\mathbf{n'}\mathbf{L}|}{\sqrt{2}\sigma}\right)<\epsilon \]

Thus, only a finite number of computational\(\mathbf{n'}\)of the periodic box is sufficient. Since the truncation here is at a distance of\(d_{\mathbf{n}}=|\mathbf{r}-\mathbf{r}_i+\mathbf{n}\mathbf{L}|\), one can prove the convergence of the short-range interaction potentials in real space by the D'Alembert discriminant (by generalizing the method of taking first a\(\delta>0\)）：\(\lim_{l_{\mathbf{n}}\rightarrow\infty}\frac{Erfc(\frac{l_{\mathbf{n}}+\delta}{\sqrt{2}\sigma})l_{\mathbf{n}}}{(l_{\mathbf{n}}+\delta)Erfc(\frac{l_{\mathbf{n}}}{\sqrt{2}\sigma})}=\lim_{l_{\mathbf{n}}\rightarrow\infty}\frac{Erfc\left(\frac{l_{\mathbf{n}+\delta}}{\sqrt{2}\sigma}\right)}{Erfc\left(\frac{l_{\mathbf{n}}}{\sqrt{2}\sigma}\right)}=e^{-\frac{\delta^2}{2\sigma^2}}\lim_{l_{\mathbf{n}}\rightarrow\infty}e^{-\frac{l_{\mathbf{n}}\delta}{\sigma^2}}=e^{-\frac{\delta^2}{2\sigma^2}}<1\). That is, the short-range interaction of potentials converges in real space. After obtaining the form of the short-range interaction potential, the potential energy of the short-range interaction can be further calculated:

\[E^S=\sum_{\mathbf{n}}\sum_{i=0}^{N-2}\sum_{j=i+1}^{N-1}\frac{q_iq_j}{4\pi\epsilon_0|\mathbf{r}_j-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}Erfc\left(\frac{|\mathbf{r}_j-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}{\sqrt{2}\sigma}\right)+\sum_{|\mathbf{n}|>0}\frac{q_i^2}{4\pi\epsilon_0|\mathbf{n}\mathbf{L}|}Erfc\left(\frac{|\mathbf{n}\mathbf{L}|}{\sqrt{2}\sigma}\right) \]

Calculation of long-range action terms

The long-range interaction potential form is obtained earlier as:

\[V_i^L(\mathbf{k})=\frac{q_i}{\epsilon_0 |\mathbf{k}|^2}e^{-j\mathbf{k}\mathbf{r}_i}e^{-\eta |\mathbf{k}|^2} \]

Also use the values in the short-range action term\(\eta=\frac{\sigma^2}{2}\).. Do an inverse Fourier transform back to real space:

\[V_i^L(\mathbf{r})=\frac{q_i}{k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-j\mathbf{k}\mathbf{r}_i}e^{-\frac{\sigma^2 |\mathbf{k}|^2}{2}}}{|\mathbf{k}|^2}e^{j\mathbf{k}\cdot\mathbf{r}} \]

Similarly one can prove the convergence of the formula according to the D'Alembert discriminant method. And with reference to the previous inverse easy space\(\mathbf{k}\)The definition has:

\[e^{j\mathbf{k}\cdot\mathbf{nL}}=e^{2j|\mathbf{n}|\pi}=1 \]

That is, the long-range interaction term can be omitted from the summation term of the periodic box, so the long-range interaction potential energy takes the form:

\[E^L=\frac{1}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}\sum_{i=0}^{N-1}q_ie^{-j\mathbf{k}\mathbf{r}_i}\sum_{l=0}^{N-1}q_le^{j\mathbf{k}\cdot\mathbf{r}_l} \]

The content of these last two summations takes the form of the inner product of two plane wave functions, the physical meaning of which is that by distributing a fixed charge in real space according to the probability amplitude to different lattice points in the inverse easy space, a charge distribution function in the inverse easy space can be defined:

\[S(\mathbf{k})=\sum_{i=0}^{N-1}q_ie^{j\mathbf{k}\cdot\mathbf{r}_i} \]

Then the writing of the long-range interaction potential can be further simplified:

\[E^L=\frac{1}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}\left|S(\mathbf{k})\right|^2 \]

One point to note is that although the long-range interaction potential is convergent here, it contains the point charge\(i\)treat (sb a certain way)\(\mathbf{r}_i\)The potential energy generated at (needs to be removed). And the long-range interaction form of the error function form that was previously available when calculating the short-range interaction:

\[\begin{align*} V_i^L(\mathbf{r})&=\sum_{\mathbf{n}}\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i-\mathbf{nL}|}-V_i^S(\mathbf{r})\\ &=\sum_{\mathbf{n}}\frac{q_i}{4\pi\epsilon_0|\mathbf{r}-\mathbf{r}_i-\mathbf{nL}|}Erf\left(\frac{|\mathbf{r}-\mathbf{r}_i-\mathbf{nL}|}{\sqrt{2}\sigma}\right) \end{align*} \]

Although not astringent, it would have been better to take a\(\mathbf{r}=\mathbf{r}_i,\mathbf{nL}=0\), which is also the charge self-interaction mentioned earlier, the long-range interaction form becomes:

\[V_i^L(\mathbf{r}_i)=\frac{q_i}{4\pi\epsilon_0}\sqrt{\frac{2}{\pi}}\frac{1}{\sigma} \]

The self-interaction term to be subtracted from the long-range interaction term in the inverse easy space can then be obtained as:

\[E^{self}=\frac{1}{2}\sum_{i=0}^{N-1}q_iV_i^L(\mathbf{r}_i)=\frac{1}{4\pi\epsilon_0}\frac{1}{\sqrt{2\pi}\sigma}\sum_{i=0}^{N-1}q_i^2 \]

total potential energy

After the previous calculations, we have obtained the short-range interaction term for real-space convergence, the long-range interaction term for inverse-easy-space convergence, and a self-interaction term that needs to be subtracted from the long-range interaction term, respectively, then the potential energy can be summarized as:

\[\begin{align*} E&=E^S+E^L-E^{self}\\ &=\sum_{\mathbf{n}}\sum_{i=0}^{N-2}\sum_{j=i+1}^{N-1}\frac{q_iq_j}{4\pi\epsilon_0|\mathbf{r}_j-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}Erfc\left(\frac{|\mathbf{r}_j-\mathbf{r}_i+\mathbf{n}\mathbf{L}|}{\sqrt{2}\sigma}\right)+\sum_{|\mathbf{n}|>0}\frac{q_i^2}{4\pi\epsilon_0|\mathbf{n}\mathbf{L}|}Erfc\left(\frac{|\mathbf{n}\mathbf{L}|}{\sqrt{2}\sigma}\right)\\ &+\frac{1}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}\left|S(\mathbf{k})\right|^2\\ &-\frac{1}{4\pi\epsilon_0}\frac{1}{\sqrt{2\pi}\sigma}\sum_{i=0}^{N-1}q_i^2 \end{align*} \]

Summary outline

This paper describes the Ewald summation method for calculating the electrostatic potential energy under periodic boundary conditions. The periodic electrostatic potential function is not a spatially convergent function, and by Ewald summation the electrostatic potential can be cut into short-range interactions and long-range interactions, and the two terms converge in the real space and the inverse easy space (or Fourier space, k-space, etc.), respectively. Further truncation can then be performed to obtain more accurate potential energy calculations with less cost.

This article was first linked to:/dechinphy/p/

Author ID: DechinPhy