Theory

Potential Energy Slicing

This section contains the formulation used to implement the SlicedNonbondedForce class. For simplicity, the equations presented in here follow the conventional 1-based indexing, whereas the class API follows Python’s 0-based indexing.

Slicing the potential energy starts with partitioning all N particles of a system into n disjoint subsets or “colors”. This procedure results in \(n(n+1)/2\) slices, distinguished by order-invariant pairs of subset indices. Slice I, J comprises all pairs involving a particle in subset I and a particle in another (or the same) subset J.

Example: three particle subsets (●●●) form six energy slices

It is straightforward to slice pairwise Lennard-Jones interactions, as well as the direct-space, self-energy, and exclusion/exception parts of the Ewald summation. Therefore, the present section focuses on the reciprocal-space and background energy parts.

Sliced Background Energy

In the case of systems with non-zero net charge, it is necessary to add the energy of interaction with a neutralizing background charge. The background energy is given by [7]

\[E_{bg} = \frac{Q^2}{8\epsilon_0 V \alpha^2}.\]

where \(Q = \sum_{i=1}^N q_i\) is the total charge of the system. We can split the total charge into the sum of the charges in each subset,

\[Q = \sum_{I=1}^n Q_I,\]

where \(Q_I = \sum_{i \in I} q_i\) is the charge of subset I.

The goal is to express the background energy as

\[E_{bg} = \sum_{I=1}^n \sum_{J=I}^n \lambda^{elec}_{I,J} E^{bg}_{I,J},\]

where \(E^{bg}_{I,J}\) is the background contribution of Slice I, J. For this, we rewrite the square of the total charge as

\[Q^2 = \sum_{I=1}^n \sum_{J=1}^n Q_I Q_J.\]

Therefore, the background energy of a slice I, J is given by

\[E^{bg}_{I,J} = (2 - \delta_{I,J}) \frac{Q_I Q_J}{8\epsilon_0 V \alpha^2},\]

where \(\delta_{I,J}\) is the Kronecker delta. The prefactor \(2 - \delta_{I,J}\) allows us to run the summation over J only for \(J \geq I\).

Sliced Reciprocal-Space Energy

For a simulation box with edge matrix \(\mathbf L\) containing N particles under periodic boundary conditions, the standard reciprocal part of the electrostatic potential energy can be expressed as

\[E_{rec} = \frac{2\pi}{\epsilon_0 V} \sum_{\mathbf n \neq \mathbf 0} \frac{e^{-\frac{k^2}{4\alpha^2}}}{k^2} \sum_{i=1}^N \sum_{j=1}^N q_i q_j e^{\text{ⅈ} {\mathbf k}\cdot({\mathbf r}_i - {\mathbf r}_j)},\]

where \(\epsilon_0\) is the vacuum permittivity, \(V\) is the box volume, \(\alpha\) is the Ewald splitting parameter, \(\mathbf n \in \mathbb Z^3\) is an integer lattice vector, \(\mathbf k = 2\pi \mathbf L^{-1}{\mathbf n}\) is a reciprocal space wave vector, \(k = \|\mathbf k\|\) is the norm of \(\mathbf k\), and \(\text{ⅈ} = \sqrt{-1}\) is the imaginary unit. As the summations over indices i and j run for all particles, we can write

\[E_{rec} = \frac{2\pi}{\epsilon_0 V} \sum_{\mathbf n \neq \mathbf 0} \frac{e^{-\frac{k^2}{4\alpha^2}}}{k^2} \Bigg(\sum_{i=1}^N q_i e^{\text{ⅈ} {\mathbf k}\cdot{\mathbf r}_i} \Bigg) \Bigg(\sum_{j=1}^N q_j e^{-\text{ⅈ} {\mathbf k}\cdot{\mathbf r}_j}\Bigg).\]

The first term inside parentheses is the charge structure factor

\[S(\mathbf k) = \sum_i q_i e^{\text{ⅈ} {\mathbf k}\cdot{\mathbf r}_i}\]

and the second one is its complex conjugate \({\overline S}(\mathbf k)\), so that we can write

\[E_{rec} = \frac{2\pi}{\epsilon_0 V} \sum_{\mathbf n \neq \mathbf 0}\frac{e^{-\frac{k^2}{4\alpha^2}}}{k^2} |S(\mathbf k)|^2.\]

The goal is to express the reciprocal-space energy as

\[E_{rec} = \sum_{I=1}^n \sum_{J=I}^n \lambda^{elec}_{I,J} E^{rec}_{I,J},\]

where \(E^{rec}_{I,J}\) is the reciprocal-space contribution of Slice I, J.

The structure factor of a subset I can be computed as

\[S_I(\mathbf k) = \sum_{i \in I} q_i e^{\text{ⅈ} {\mathbf k}\cdot{\mathbf r}_i}.\]

Since we have disjoint subsets, this makes \(S(\mathbf k) = \sum_{I=1}^n S_I(\mathbf k)\). Substituting into \(S(\mathbf k) {\overline S}(\mathbf k)\) and expanding the product, all imaginary terms cancel out, as expected. The reciprocal-space energy then becomes

\[E_{rec} = \frac{2\pi}{\epsilon_0 V} \sum_{\mathbf n \neq \mathbf 0} \frac{e^{-\frac{k^2}{4\alpha^2}}}{k^2} \sum_{I=1}^n \sum_{J=1}^n S_I(\mathbf k) \cdot S_J(\mathbf k),\]

where \(x \cdot y = \text{Re}(x)\text{Re}(y) + \text{Im}(x)\text{Im}(y)\) for two complex numbers x and y. It is now clear that the energy of a slice I, J should be calculated by

\[E^{rec}_{I,J} = (2-\delta_{I,J}) \frac{2\pi}{\epsilon_0 V} \sum_{\mathbf n \neq \mathbf 0} \frac{e^{-\frac{k^2}{4\alpha^2}}}{k^2} S_I(\mathbf k) \cdot S_J(\mathbf k).\]

Sliced PME Implementation

Each subset-specific structure factor \(S(\mathbf k)\) can be calculated using the smooth Particle Mesh Ewald method [8]. This requires n FFT calculations for evaluating all slices. In the case of forces, n additional inverse FFT calculations are necessary. Therefore, in a single processor and for a fixed number of particles, the computational effort scales approximately linearly with the number of subsets (not the number of slices). In a GPU, as with the OpenCL and CUDA instances of openmm.Platform, we take advantage the capability of either package vkFFT or cuFFT to speed-up parallel FFT calculations via batched transforms.

Inner Force Group Energy

When an inner force group is enabled, the inner force of a pair is \(S(r)\) times the bare Lennard-Jones plus Coulomb force, where \(S\) is the quintic switch of the force-directed switching scheme of [6]. Writing the bare potential as \(V(r) = \sum_n c_n r^{-n}\) with \(n \in \{12, 6, 1\}\), the energy whose gradient is that force, shifted to vanish at the inner cutoff \(r_c\), is \(\hat V(r) = \sum_n c_n e_n(r)\) with

\[e_n(r) = n \int_r^{r_c} S(s)\, s^{-(n+1)}\, ds \quad (r_s < r < r_c), \qquad e_n(r) = r^{-n} - K_n \quad (r \le r_s), \qquad K_n = r_s^{-n} - e_n(r_s),\]

where \(r_s\) is the inner switching distance. Expanding \(S\) in powers of \(s\) and integrating term by term gives a closed form whose terms are much larger than their sum: the expansion coefficients grow like \((r_s/h)^k\) with \(h = r_c - r_s\), and the integration constants involve integrals from \(r_s\). In single precision this loses most of the significant digits even for a 0.5–0.8 nm switching range, and in double precision it fails for narrow ranges or small \(r_s\).

The switch is the smoothstep polynomial \(P(u) = 10u^3 - 15u^4 + 6u^5\) of \(u = (r_c - s)/h\). In the variable \(y = 1 - s/r_c\), that is \(S = P(\gamma y)\) with \(\gamma = r_c/h\), and \(s^{-(n+1)} ds = -r_c^{-n} (1-y)^{-(n+1)} dy\), so that

\[e_n(r) = n\, r_c^{-n} \left( 10\gamma^3 J_3 - 15\gamma^4 J_4 + 6\gamma^5 J_5 \right), \qquad J_k = \int_0^Y y^k (1-y)^{-(n+1)}\, dy, \qquad Y = 1 - \frac{r}{r_c}.\]

The only cancellation left is the alternation of the three smoothstep coefficients, which is bounded by a factor of about 30 because \(\gamma y \le 1\) over the integration range. For the Lennard-Jones powers, \(n \ge 6\), the substitution \(x = y/(1-y)\) turns each \(J_k\) into a polynomial with positive coefficients in \(x = (r_c - r)/r\),

\[J_k = \int_0^x t^k (1+t)^{n-k-1}\, dt = \sum_{j=0}^{n-k-1} \binom{n-k-1}{j} \frac{x^{k+j+1}}{k+j+1},\]

so that \(e_{12}\) and \(e_6\) are \(x^4\) times polynomials of degrees 8 and 2. For the Coulomb power, \(n = 1\), the integral has a logarithm, and the three \(J_k\) share it:

\[e_1(r) = \frac{\gamma^3}{r_c} \left[ \frac{Y^4 (10 - 15\gamma Y + 6\gamma^2 Y^2)}{1-Y} - 30 \left( (\gamma-1)^2 L_6(Y) + (1-2\gamma)\frac{Y^5}{5} + \frac{Y^4}{4} \right) \right], \qquad L_6(Y) = \sum_{j \ge 6} \frac{Y^j}{j} = -\ln(1-Y) - \sum_{j=1}^{5} \frac{Y^j}{j}.\]

The function \(L_6(Y)/Y^6 = \sum_{m \ge 0} Y^m/(m+6)\) is the Lerch transcendent \(\Phi(Y, 1, 6)\), equivalently \({}_2F_1(1, 6; 7; Y)/6\), and \(L_6(Y)\) itself is the incomplete beta function \(B_Y(6, 0)\). None of these is available in the GPU math libraries, and the logarithm, the only elementary function in the family, requires the subtraction above, which cancels for small \(Y\). Since the function does not depend on the switching range, the GPU kernels evaluate it from a fixed Chebyshev expansion on \(0 \le Y \le 0.8\), accurate to \(5 \times 10^{-10}\), and from the logarithm above that, where the subtraction loses at most one digit. The Reference platform sums the series for \(Y < 0.7\) and uses \(-\ln(1-Y)\) minus the leading terms above that. With these expressions, \(e_{12}\) and \(e_6\) are accurate to a few units in the last place and \(e_1\) to about \(3 \times 10^{-6}\) in single precision, for any switching range. In double precision, \(e_{12}\) and \(e_6\) are exact to rounding on all platforms and \(e_1\) on the Reference platform, while the Chebyshev expansion limits \(e_1\) on the GPU platforms to an absolute error of about \(3 \times 10^{-10}\) (about \(10^{-9}\) relative). These errors affect only the energies: the inner force is computed directly as \(S(r)\) times the bare force. The Reference platform and the constants of the GPU kernels are computed by the same code, in InnerSplit.h.

In a decapped slice (Capped Lennard-Jones potential), the Lennard-Jones potential below \(\sigma\) is not the bare one. Since \(\sigma\) cannot exceed \(r_s\), where \(S = 1\), the inner group takes the whole interaction there: the capped force with energy \(u_{cap}(r) - K_{LJ}\), where \(K_{LJ} = 4\epsilon(\sigma^{12} K_{12} - \sigma^6 K_6)\) is the constant the remainder reports for \(r \le r_s\), and all of the decapping term. The remainder is the same as for the bare potential. Coulomb decapping, which caps the \(1/r\) term of the pairs of a slice in the same way, is handled likewise: the inner group takes the capped Coulomb force with energy \(q_i q_j [c(r/\sigma)/\sigma - K_1]/(4 \pi \epsilon_0)\) and all of the Coulomb decapping term. Below \(r_s\) the remainder is \(q_i q_j [K_1 - \mathrm{erf}(\alpha r)/r]/(4 \pi \epsilon_0)\) with Ewald, PME and LJPME, evaluated from the series of \(\mathrm{erf}(\alpha r)/r\) at small \(\alpha r\) so that it stays finite and accurate as \(r \to 0\), and a constant plus the reaction-field term otherwise.