Lattice-point Lagrangian interpolation with PME algorithm

Technical background

In a previousblog (loanword)In , we introduced the basic origin and representation of Lagrange interpolation. Here we are going to introduce an application scenario of Lagrange interpolation: lattice-point Lagrange interpolation. The advantage of this scenario is that if we want to sum or integrate over the whole space of real numbers, the computation grows with the shape of the variables. For example, in molecular dynamics simulationsCalculate the electrostatic potential energy, just calculating the charge distribution function is a\(O(N^2)\)The amount of computation in which\(N\)represents the number of point charges. And if we discretize the space, divide it into a series of lattice points, and interpolate the constant number of neighboring lattice points, then our summation computation can be reduced to\(O(N)\)。

lattice-point Lagrangian interpolation

Given a function\(y=f(x-x_r)\), which we can interpolate to the nearest 4 integer grid points:\(\lfloor x_r\rfloor-1.5,\lfloor x_r\rfloor-0.5,\lfloor x_r\rfloor+0.5,\lfloor x_r\rfloor+1.5\), according to the Lagrangian interpolation form there:

\[y_{interp}=c_1(x)f(\lfloor x_r\rfloor-1.5-x_r)+c_2(x)f(\lfloor x_r\rfloor-0.5-x_r)+c_3(x)f(\lfloor x_r\rfloor+0.5-x_r)+c_4(x)f(\lfloor x_r\rfloor+1.5-x_r) \]

provided that\(\lfloor x_r\rfloor\)The nearest center point is the origin, i.e.\(\lfloor x_r\rfloor=0\), then its coefficient has:

\[\begin{align*} c_1(x)&=\frac{(x-\lfloor x_r\rfloor+0.5)(x-\lfloor x_r\rfloor-0.5)(x-\lfloor x_r\rfloor-1.5)}{-6}=\frac{1}{48}(-8x^3+12x^2+2x-3)\\ c_2(x)&=\frac{(x-\lfloor x_r\rfloor+1.5)(x-\lfloor x_r\rfloor-0.5)(x-\lfloor x_r\rfloor-1.5)}{2}=\frac{1}{16}(8x^3-4x^2-18x+9)\\ c_3(x)&=\frac{(x-\lfloor x_r\rfloor+1.5)(x-\lfloor x_r\rfloor+0.5)(x-\lfloor x_r\rfloor-1.5)}{-2}=\frac{1}{16}(-8x^3-4x^2+18x+9)\\ c_4(x)&=\frac{(x-\lfloor x_r\rfloor+1.5)(x-\lfloor x_r\rfloor+0.5)(x-\lfloor x_r\rfloor-0.5)}{6}=\frac{1}{48}(8x^3+12x^2-2x-3) \end{align*} \]

Its image is roughly shown below (image from reference link 1):

For multidimensional lattice-point Lagrange interpolation, it is a cross-multiplication relation with the image:

Truncation of remote interaction terms

We apply this lattice point Lagrangian interpolation obtained above to the calculation of the electrostatic potential energy. In a previousblog (loanword)One of the electrostatic potential calculations presented is a charge distribution function that looks like this:

\[s(\mathbf{k})=|S(\mathbf{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}\mathbf{r}_l} \]

included among these\(S(\mathbf{k})=\sum_{i=0}^{N-1}q_ie^{j\mathbf{k}\mathbf{r}_i}=\sum_{i=0}^{N-1}q_ie^{j\mathbf{k}_xx_i}e^{j\mathbf{k}_yy_i}e^{j\mathbf{k}_zz_i}\). Replacing these latter exponential terms with lattice point Lagrangian interpolation yields:

\[S(\mathbf{k})=\sum_{i=0}^{N-1}q_i\sum_{x,y,z}\left[c_1(x)f(\lfloor x_i\rfloor-1.5-x_i)+c_2(x)f(\lfloor x_i\rfloor-0.5-x_i)+c_3(x)f(\lfloor x_i\rfloor+0.5-x_i)+c_4(x)f(\lfloor x_i\rfloor+1.5-x_i)\right]\left[c_1(y)f(\lfloor y_i\rfloor-1.5-y_i)+c_2(y)f(\lfloor y_i\rfloor-0.5-y_i)+c_3(y)f(\lfloor y_i\rfloor+0.5-y_i)+c_4(y)f(\lfloor y_i\rfloor+1.5-y_i)\right]\left[c_1(z)f(\lfloor z_i\rfloor-1.5-z_i)+c_2(z)f(\lfloor z_i\rfloor-0.5-z_i)+c_3(z)f(\lfloor z_i\rfloor+0.5-z_i)+c_4(z)f(\lfloor z_i\rfloor+1.5-z_i)\right] \]

With the functional form, we can abbreviate\(S(\mathbf{k})\)for a summation over three-dimensional spatial lattice points:

\[S(\mathbf{k})=\sum_{i=0}^{N-1}q_i\sum_{m_x=\lfloor x_{min}\rfloor-1.5}^{\lfloor x_{max}\rfloor+1.5}\sum_{m_y=\lfloor y_{min}\rfloor-1.5}^{\lfloor y_{max}\rfloor+1.5}\sum_{m_z=\lfloor z_{min}\rfloor-1.5}^{\lfloor z_{max}\rfloor+1.5}c_{m_x}(m_x)e^{jk_xm_{x}}c_{m_y}(m_y)e^{jk_ym_{y}}c_{m_z}(m_z)e^{jk_zm_{z}} \]

Take the coefficient terms separately again:

\[Q(m_x,m_y,m_z)=\sum_{i=0}^{N-1}q_ic_{m_x}(m_x)c_{m_y}(m_y)c_{m_z}(m_z) \]

Here.\(Q\)is actually a shape as\((N_x,N_y,N_z)\)The tensor of the\(m_x,m_y,m_z\)corresponds to the tensor index of a particular lattice point, and the tensor elements corresponding to each index are computed by the coefficient function, and with such a concept, rewrite the\(S(\mathbf{k})\)The function of the

\[S(\mathbf{k})=\sum_{m_x=\lfloor x_{min}\rfloor-1.5}^{\lfloor x_{max}\rfloor+1.5}\sum_{m_y=\lfloor y_{min}\rfloor-1.5}^{\lfloor y_{max}\rfloor+1.5}\sum_{m_z=\lfloor z_{min}\rfloor-1.5}^{\lfloor z_{max}\rfloor+1.5}Q(m_x,m_y,m_z)e^{j\mathbf{k}\cdot\mathbf{m}} \]

We will see that this interpolation comes out\(S(\mathbf{k})\)function actually computes the tensor\(Q\)exist\(\mathbf{k}\)at the Fourier transform, then it is possible to further abbreviate\(S(\mathbf{k})\)The form of the

\[S(\mathbf{k})=VF_{\mathbf{k}}^{*}(Q)(m_x,m_y,m_z) \]

included among these\(F^{*}\)denotes the inverse Fourier transform.\(V\)denotes the inverse Fourier transform normalization constant. Following the previous 4-lattice-point Lagrangian interpolation, at this point the obtained\(S(\mathbf{k})\)value is a tensor of shape (4, 4, 4), which means that each of the 64 lattice points for the inverted lattice vector\(\mathbf{k}\)contribution (the effect of interpolated individual point charges). Then similar can be obtained:

\[s(\mathbf{k})=VF_{\mathbf{k}}^{*}(Q)(m_x,m_y,m_z)F_{\mathbf{k}}(Q)(m_x,m_y,m_z)=V|F_{\mathbf{k}}(Q)(m_x,m_y,m_z)|^2 \]

Substituting into the Ewald form of the long-range interaction term (see thiswritings) is available in:

\[\begin{align*} E^L&=\frac{1}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}s(\mathbf{k})\\ &=\frac{V}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}|F_{\mathbf{k}}(Q)(m_x,m_y,m_z)|^2 \end{align*} \]

This is the trick of calculating the long-range interaction potential in the Particle-Mesh-Ewald method calculations. Since\(\mathbf{k}\)space cannot converge quickly, then reducing the computational complexity of the charge distribution term can likewise have the effect of substantial computational savings.

Truncation of the short-range interaction term

In the previous article on Ewald summation we described the different forms of convergence of the electrostatic potential energy calculation after splitting it into long-range, short-range and self-interaction terms. The simplification of the computational effort has been accomplished by the long range interaction term through the above sections, and there is also a short range interaction term\(E^{S}\), we know that the spacing of the short-range interaction terms with respect to the atomic real space is rapidly convergent, and we have already computed the nearest-neighbor table for a given cutoff truncation once in the calculation of the LJ potential energy. It is then easy to consider introducing the concept of a nearest-neighbor table and use this nearest-neighbor table directly to make a cutoff of the short-range interaction term for the electrostatic potential energy. Thus the short-range interaction term can be written as:

\[\begin{align*} 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)\\ &\approx \sum_{i,j\in \{Neigh\}}\frac{q_iq_j}{4\pi\epsilon_0|\mathbf{r}_j-\mathbf{r}_i|}Erfc\left(\frac{|\mathbf{r}_j-\mathbf{r}_i|}{\sqrt{2}\sigma}\right) \end{align*} \]

A presupposition here is that\(d_{cutoff}<<L_{pbc}\), so omitting the periodic box within the other boxes of the\(i\)Charge on the center box\(\mathbf{r}_i\)Role of the Branch.

Particle-Mesh-Ewald

Based on the approximate remote interaction terms and after the short-range interaction terms obtained in the above section, we can rewrite the total electrostatic potential energy in the PME (Particle-Mesh-Ewald) algorithm:

\[\begin{align*} E&=E^S+E^L-E^{self}\\ &=\sum_{i,j\in \{Neigh\}}\frac{q_iq_j}{4\pi\epsilon_0|\mathbf{r}_j-\mathbf{r}_i|}Erfc\left(\frac{|\mathbf{r}_j-\mathbf{r}_i|}{\sqrt{2}\sigma}\right)\\ &+\frac{V}{2k_xk_yk_z\epsilon_0}\sum_{|\mathbf{k}|>0}\frac{e^{-\frac{\sigma^2 k^2}{2}}}{k^2}|F_{\mathbf{k}}(Q)(m_x,m_y,m_z)|^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 presents the fundamentals of reducing the computational complexity of the electrostatic potential energy term in molecular force fields using the Particle Mesh Ewald algorithm based on lattice-point Lagrangian interpolation. The computation of electrostatic potential energy is split into a remote interaction term and a short-range interaction term in the Ewald summation framework, where the short-range interaction term converges rapidly about the point charge spacing in real space, while the remote interaction term converges slowly in inverse easy space. Therefore, in the calculation of remote interactions, the interpolation method can be used to reduce the computational complexity of individual inverted-easy lattice points, thus allowing the overall remote interaction term calculation to converge quickly as well.

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

Author ID: DechinPhy