Performance Evaluation of the Zero-Multipole Summation Method in Modern Molecular Dynamics Software

TL;DR

This paper compares the zero-multiple summation method (ZMM) with the smooth particle mesh Ewald method (SPME) in molecular dynamics, showing ZMM's competitive or superior performance for electrostatic calculations.

Contribution

It provides a practical performance evaluation of ZMM, highlighting its efficiency and potential as a fast alternative for electrostatic computations in molecular dynamics.

Findings

ZMM performance is comparable or better than SPME.

ZMM with quadrupole or octupole cancellation is highly efficient.

ZMM is suitable for fast electrostatic potential calculations.

Abstract

We evaluate the practical performance of the zero-multiple summation method (ZMM), a method for approximately calculating electrostatic interactions in molecular dynamics simulations. The performance of the ZMM is compared with that of the smooth particle mesh Ewald method (SPME). Even though the ZMM uses a larger cutoff distance than the SPME does, the performance of the ZMM is found to be comparable to or better than that of the SPME. In particular, the ZMM with quadrupole or octupole cancellation and no damping factor is an excellent candidate for the fast calculation of electrostatic potentials.

Equations

$$V_{\mathrm{ZM},ij}=\frac{q_i{q}_j}{4\pi {ϵ}_0}{U}_{\mathrm{ZM}}\left(r_{ij}\right)$$

$$U_{\mathrm{ZM}}\left(r\right)=\{\begin{array}{ll}\frac{\mathrm{erfc}\left(\alpha r\right)}{r}+\sum _{k=0}^{\ell }{c}_k{r}^{2k}& \left(r\le r_c\right)\\ 0,& \left(r>r_c\right)\end{array}$$

$$\mathrm{erfc}\left(z\right)=1-\frac{2}{\sqrt{\pi }}{\int }_0^z\exp \left(-x^2\right)\mathrm{d}x,$$

$$V_{\mathrm{SPME}}=V_{\mathrm{real}}+V_{\mathrm{recp}},$$

$$V_{\mathrm{real}}=\sum _{i=1}^N\sum _{j=i+1}^N{V}_{\mathrm{real},ij}=\sum _{i=1}^N\sum _{j=i+1}^N\frac{q_i{q}_j}{4\pi {ϵ}_0}\frac{\mathrm{erfc}\left(\beta r_{ij}\right)}{r_{ij}}.$$

$$V_{\mathrm{recp}}=\frac{1}{2}\sum _{k_1=0}^{K_1-1}\sum _{k_2=0}^{K_2-1}\sum _{k_3=0}^{K_3-1}D\left(k_1,k_2,k_3\right){\left|\mathcal{F}\left(k_1,k_2,k_3\right)\right|}^2,$$

$$\mathcal{F}\left(k_1,k_2,k_3\right)=\sum _{m_1=0}^{K_1-1}\sum _{m_2=0}^{K_2-1}\sum _{m_3=0}^{K_3-1}Q\left(m_{1,}{m}_2,m_3\right)$$

$$\times \exp \left[-2\pi \sqrt{-1}\left(\frac{m_1{k}_1}{K_1}+\frac{m_2{k}_2}{K_2}+\frac{m_3{k}_3}{K_3}\right)\right].$$

$$\frac{\mathrm{erfc}\left(\alpha r\right)}{r}=\frac{1}{r}-\frac{\mathrm{erf}\left(\alpha r\right)}{r}=\frac{1}{r}-\alpha \cdot \frac{\mathrm{erf}\left(\alpha r\right)}{\alpha r},$$

$$s=\frac{(\text{Time to run a single MD step})}{(\text{\# of atoms in the system})/(\text{\# of CPU cores})}.$$

$$⟨\Delta V_{ij}⟩$$

$$\approx 4\pi {\left(r_{\mathrm{nl}}+{\sigma }_{ij}\right)}^2{\rho }_j{\int }_{-\infty }^{r_c}\mathrm{d}{r}_t{\int }_{r_{\mathrm{nl}}}^{\infty }\mathrm{d}{r}_0\left[\sum _{k=1}^{n_k}{V}_{ij}^{\left(k\right)}\left(r_c\right)\cdot \frac{1}{k!}{\left(r_t-r_c\right)}^k\right]$$

$$\times G\left(\frac{r_t-r_0}{{\sigma }_{ij}}\right),$$

Full text

††journal: Computer Physics Communications

1 Introduction

In molecular dynamics (MD) simulations, molecules are modeled in a computer using interaction potentials between atoms, and their dynamics are simulated by numerically solving the Newtonian equations of motion. In MD simulations, a significant amount of computation time is consumed during the evaluation of electrostatic interactions. Because accurate evaluation of electrostatic interactions is vital for an accurate MD simulation, computationally efficient yet reasonably accurate ways of modeling and approximating electrostatic interactions have been of much interest in the development of MD methods. Such methodologies include the particle–particle particle–mesh method (PPPM),[11] smooth particle-mesh Ewald method (SPME),[10, 12] Gaussian-split Ewald method,[45] and fast multipole method,[21, 14] to name a few. Among them, the PPPM and SPME are arguably the most commonly used techniques in current MD simulations. Current state-of-art MD simulation packages such as AMBER,[9] CHARMM,[8, 7, 29] GROMACS,[6, 36, 1] LAMMPS,[42] and NAMD[39] implement either (or both) of these algorithms. Roughly speaking, the PPPM and SPME accelerate the Ewald summation by splitting the electrostatic interaction into a short-ranged real part and a long-ranged reciprocal part, where the former is computed directly and the latter is computed approximately by convolution via a three-dimensional fast Fourier transform (FFT). Unfortunately, the FFT requires several global communication steps for the “butterfly” exchange of intermediate values. In modern computers, communication between processors is much costlier than computations, and thus methods that can avoid FFT offer huge advantages particularly when executed in parallel.

Historically, before the development of these efficient mesh-based methods, electrostatic interactions between atoms were typically calculated by only accumulating the interactions within some cutoff distance. These methods are called as cutoff-based methods.[18] Cutoff-based methods inherently do not require the FFT, and therefore, have a good scalability even in massively parallel environments. However, these cutoff-based methods have been mostly abandoned because of serious artifacts. For example, artifacts originated from layer structure in liquid water systems were reported when the plain-cutoff[52, 53, 47] or the switching force[53] is used. Similarly, properties of several water models deviate between the reaction-field method and the Ewald-based methods.[4, 24] With the introduction of efficient Ewald-based methods such as the PPPM and the SPME, cutoff-based methods have been rarely used in modern molecular simulations.

Recently, several cutoff-based methods[51, 50, 23, 46, 16, 17] that can avoid such artifacts have been proposed. The zero-multipole summation method (ZMM) is one of such examples.[16] The ZMM utilizes the neutralization of electrostatic multipoles, which is a reasonable assumption in solvents such as water-containing systems, to compute electrostatic interactions. The electrostatic energies in the ZMM have been shown to match those in the SPME to within a small deviation,[17] and both the static and dynamic properties of bulk water computed by the ZMM have been shown to agree with the SPME,[49] which makes the ZMM a promising candidate for computing approximate electrostatic potentials. The drawbacks of the ZMM are that it may require a slightly larger cutoff in real-space particle-particle interactions than does the SPME for sufficiently accurate computations. For example, Wang et al. [49] used a cutoff distance of 1.2 nm for validation of the ZMM, whereas in the SPME the real-space cutoff distance may vary and is typically set to 1.0 nm or shorter for good performance. The ZMM’s practical performance properties were thus not apparent. Although it was preliminarily reported that the zero-dipole summation method,[31] a special case of the ZMM, runs with performance comparable to that of the SPME, there have been no extensive software engineering efforts to optimize the performance of both methods and no comparison of the ZMM using higher order multipole parameters. Especially, in the case of the SPME, it has been reported that engineering efforts are crucial for high-performance implementations.[40, 39, 27, 43, 36] A rigorous comparison with the state-of-art implementation of the SPME is thus needed.

In the present article, we report the performance of the ZMM implemented in a molecular dynamics software package, GROMACS.[6, 25, 36, 1] GROMACS has been implemented with a mechanism to separate real and reciprocal ranks in the SPME calculation, 2D-decomposition of 3D FFT calculation for the reciprocal calculation, SIMD-optimized non-bonded interaction calculations, and fast short-ranged Ewald kernel to accelerate the SPME calculation.[27, 37, 1, 36] It thus contains one of the state-of-the-art implementations of the SPME. We compare the performance of the ZMM with that of the SPME to elucidate the differences at a practical level.

In sections 2.1 and 2.2, we briefly review the formulae and the reported properties of the ZMM and the SPME, respectively. We then explain the performance optimization of the ZMM and the SPME in sections 2.3 and 2.4, respectively. The system setup is shown in section 2.5, and the accuracy measurement is reported in section 2.6. The protocol of the performance measurement, which is the main issue of the study, is explained in section 2.7. Results and discussion for the performance measurement follows in section 3. Finally, we conclude the current work in section 4.

2 Methods

2.1 Zero-Multipole Summation Method

2.1.1 Formula

The ZMM assumes electric multipole neutralization in a system for fast computation of the electrostatics. More rigorously, for any particle $i$ ($1\leq i\leq N$ where $N$ is the number of particles), the ZMM assumes that there exists a particle set $\mathcal{M}_{i}^{\left(\ell\right)}$ such that the $\ell$th and lower order electrostatic multipoles111We refer to charges as the 0th order multipoles, dipoles as the 1st order multipoles, quadrupoles as the 2nd order multipole, and so on. are neutralized for $\mathcal{M}_{i}^{\left(\ell\right)}\cup\left\{i\right\}$, and the particles in $\mathcal{M}_{i}^{\left(\ell\right)}$ are within distance $r_{c}$ from particle $i$ (hereinafter this assumption is called the neutralization principle). The particles in $\mathcal{M}_{i}^{\left(\ell\right)}$ are further assumed to nearly fill the (cutoff) sphere around $i$ to a radius of $r_{c}$. Under these assumptions, the electrostatic interaction between particle $i$ and particle $j$ is suitably reformulated if their distance is less than $r_{c}$, while that is approximately replaced by a correction term if their distance is farther than $r_{c}$.

The potential function in the ZMM is then formed to approximate the electrostatic interaction. The obtained electrostatic interaction potential $V_{\mathrm{ZM},ij}$ between particle $i$ and $j$ is

[TABLE]

where $\epsilon_{0}$ is the dielectric constant of the vacuum, $q_{i}$ and $q_{j}$ are the charges of atoms $i$ and $j$ respectively, $r_{ij}$ is the distance between particles $i$ and $j$, $r_{c}$ is the cutoff distance, $\alpha$ is the real-space damping factor, $\mathrm{erfc}\left(\cdot\right)$ is the complementary error function defined as

[TABLE]

and $c_{k}$ are constants that depend on $\ell$, $\alpha$, and $r_{c}$. Because only the interactions within the cutoff radius $r_{c}$ are considered, the time complexity to run the ZMM per atom is expected to be constant, thus the total time to compute the interactions scales with $O\left(N\right)$.

2.1.2 Accuracy

This study focuses on the computational performance of the ZMM, since its accuracy has already been clarified in previous studies. We briefly review those results in this section. The accuracy of the ZMM is clear when applied to systems that are governed purely by electrostatic interactions. The method was used to calculate the Madelung energy of ionic crystal systems[16]. In these pure electrostatic systems, the neutralization principle works well, and the accuracy is improved by increasing the multipole order $\ell$. The ZMM gives a high accuracy, and the absolute error rate from the exact value is, for example, of the order $10^{-6}$ at around $r_{c}=4a$, where $a$ is the nearest-neighbor ionic distance in NaCl crystal. This is significantly more accurate than conventional real-space methods, including the cubic minimum image convention, and methods by Harrison [22] and Gaio and Silvestrelli [19].

The best accuracy of the ZMM in liquids is attained at $\ell=2$ or 3 with $\alpha\leq 1$ nm*-1*. Ionic liquids and pure water systems are governed by both electrostatic interactions and the thermodynamic environment.[17] In these systems, neutralized states are dominant but not always realized at a higher order multipole level due to thermal fluctuations, and so a moderate value of multipole order is efficient. The accuracies in the TIP3P system have been investigated in detail.[49] Thermodynamic quantities, including the chemical potential, heat capacity, and isothermal compressibility, and dynamical quantities, including the diffusion constant and viscosity, obtained by the ZMM with $r_{c}=1.2$ nm agreed with the results of a highly optimized SPME computation within statistical uncertainty at a confidence level of 95%. The dielectric properties were carefully investigated using long-time simulations[20] taking into account the correction formula[34] and system size effects.[47] The dielectric constant agreed with the SPME result, and the Kirkwood-G factor, which is very sensitive to the electrostatic treatment, also showed accurate results beyond two or three times the cutoff length and they are qualitatively excellent compared with the reaction field method with the group-based cutoff. The density deviation between the ZMM with $r_{c}=1.2$ nm and the optimized SPME was of the order of 0.1%. The radial distribution error in the ZMM was less than 0.2% with $\ell=2$, and was considerably smaller with $\ell=3$. The ZMM with $\ell=1$ has been applied to biomolecule simulations, where the systems are filled with water molecules in realistic explicit solvent studies. In explicit solvated systems of DNA[5] and a membrane protein,[31] the accuracies of energetic and dynamical properties for each scale level, such as molecule, residue, or nucleotide, have been confirmed by comparison with the SPME method.

2.2 Smooth Particle Mesh Ewald Method

2.2.1 Formula

The SPME method approximates the total electrostatic interaction energies in a periodic boundary condition as a sum of two energies:

[TABLE]

where $V_{\mathrm{real}}$ is the real-space electrostatic energy, and $V_{\mathrm{recp}}$ is the reciprocal-space electrostatic energy. The real-space electrostatic energy $V_{\mathrm{real}}$ is represented as a sum of short-ranged interaction potentials as

[TABLE]

Because $V_{\mathrm{real},ij}$ decreases rapidly as $r_{ij}$ becomes large for typically used value of $\beta$, the real-space electrostatic energy can be computed by simply summing interaction potentials from the atoms within a cutoff distance. As in the case of the ZMM, we denote this cutoff length as $r_{c}$.

The reciprocal-space electrostatic energy $V_{\mathrm{recp}}$ is calculated by grid-wise summation with a grid size $K_{1}\times K_{2}\times K_{3}$,

[TABLE]

where $D$ is a modified Green’s function in reciprocal-space and $\mathcal{F}$ is a Fourier-transformed grid charge distribution $Q$:

[TABLE]

The formulations of $D$ and $Q$ are described elsewhere.[12]

2.2.2 Computational Cost

The Fourier transform (7) is typically computed by a parallel three-dimensional FFT. Although the time complexity of the FFT calculation per process scales with $O\left(K_{1}K_{2}K_{3}P^{-1}\log K_{1}K_{2}K_{3}\right)$, where $P$ is the number of processors, and the total communication amount scales with $O\left(K_{1}K_{2}K_{3}P^{-1}\right)$, the number of communication messages per process grows as $O\left(P\right)$.[38] Thus, even with a fixed grid size, the total amount of time needed to initiate communication increases as the number of processors increases until it eventually becomes the dominant part of the calculation for large $P$. This makes the three-dimensional FFT costly in terms of the wall-clock time to complete the calculation in a parallel environment.

To reduce the increased communication cost, many SPME implementations separate processes into two groups, namely, real-space processes and reciprocal-space processes.[40, 27] The former evaluates $V_{\mathrm{real}}$ and the corresponding forces, and sends atomic coordinates to the reciprocal-space processes. The latter computes the reciprocal space energy and forces, and sends them back to the real-space nodes. In this way, the total number of processes involved in the FFT can be decreased. For large $P$, communication between the real-space and reciprocal-space processes, as well as intra-communication within the reciprocal process group, may occupy a significant amount of time. We will illustrate this problem in section 3.

2.3 Implementation and Optimization of the ZMM

We implemented the ZMM in GROMACS 5.0.7. Only the cases with $\ell=2,3$ were implemented for performance comparison. Since the ZMM with $\alpha=0$ can be computed efficiently because $\mathrm{erfc}\left(0\cdot r\right)/r=1/r$, we implemented the specialized version, too, as it can significantly decrease the number of floating point operations during the force and potential evaluation. To compute $\mathrm{erfc}\left(\alpha r\right)/r$ quickly, we reused a numerical trick from the SPME real-space calculation in GROMACS:

[TABLE]

where $\mathrm{erf}\left(\cdot\right)=1-\mathrm{erfc}\left(\cdot\right)$ is the error function. The term $\mathrm{erf}\left(x\right)/x$ is then computed using a rational polynomial approximation. The corresponding force is computed similarly from another rational polynomial approximation. Using these techniques, a computational kernel for the ZMM electrostatic force with $\ell=2$, $\alpha=0$ was implemented with 33 floating point operations (flops) per interaction, and that for the ZMM with $\ell=2$, $\alpha\neq 0$ was implemented with 65 flops. They are only 2 and 4 extra flops compared to the abrupt cutoff (31 flops)222The electrostatic force kernel for the abrupt cutoff in GROMACS shares the code with the one for the reaction field method. Due to this, the flops required for the abrupt cutoff kernel is increased from 29 flops to 31 flops. and the SPME real-space calculation (61 flops), respectively. In GROMACS, highly optimized computational kernels are implemented using the SIMD (single-instruction, multiple-data) instructions of modern processors. Thus, we implemented a SIMD version of the ZMM to maximize the computational performance. The implementation is presented in the supporting information and also at the website of one of the authors.[44]

As stated in the introduction, the bulk of MD computation time is consumed by the electrostatic interaction. In the innermost loop of the MD calculation in GROMACS, interactions between the clusters of atoms are considered.[37] To compute the interaction efficiently, the program must keep a list of neighboring atom clusters – which is often called as the (Verlet) neighbor list.[48, 3, 37] The neighbor lists are typically constructed to store a list of atoms within the distance that is the sum of the cutoff distance and some buffer length; that is, $r_{\mathrm{nl}}=r_{c}+r_{b}$, where $r_{\mathrm{nl}}$ is the total length of the neighbor list and $r_{b}$ is the buffer length. The neighbor list is updated periodically to save computational time with a certain interval. The choice of $r_{b}$ is crucial to balance the accuracy and speed of the computation. For accurate computation, $r_{b}$ must be sufficiently large so that all atoms within the cutoff distance are included in the neighbor list until the next list update. Meanwhile, large $r_{b}$ leads to a performance degradation due to the increased calculation cost. In GROMACS 5, the buffer length $r_{b}$ is determined so that the expected energy drift in the simulation does not exceed the limit.[2, Section 3.4.2] This scheme is often critical to attain a good performance, and we also utilized the scheme for the ZMM. The details are described in A.

We also optimized the neighbor list update interval. The performance of the ZMM is measured with a neighbor list update interval of 10, 20, or 40 steps, and the fastest parameter was selected. Although GROMACS supports hybrid parallelism with both OpenMP and MPI, preliminary tests showed that flat MPI execution is better than execution with OpenMP under most conditions. We therefore did not use OpenMP parallelization for the ZMM.

The following ZMM parameters were chosen for the comparison. (1) $\ell=2,$ $r_{c}=1.2$ nm and $\alpha=0.0$ nm*-1*; (2) $\ell=3$, $r_{c}=1.2$ nm and $\alpha=0.0$ nm*-1*; (3) $\ell=2$, $r_{c}=1.2$ nm and $\alpha=1.0$ nm*-1*; (4) $\ell=2$, $r_{c}=1.4$ nm and $\alpha=0.0$ nm*-1*. In the present research, we named these settings (1) ZQuad, (2) ZOct, (3) ZQuad-$\alpha$1, and (4) ZQuad-1.4. ZQuad and ZOct represent the ZMM runs with the most-typical parameters. ZQuad-$\alpha$1 and ZQuad-1.4 were tested to show the performance difference with a non-zero dumping factor and with an increased cutoff length. We also benchmarked conditions with $\ell=2,$ $r_{c}=1.2$ nm and different $\alpha$ values, $\alpha=0.5$ nm*-1* and $\alpha=1.5$ nm*-1* in a preliminary analysis for the case of 10,000 TIP3P molecules (see section 2.5). The results showed that the performance differences in comparison with condition (3) were small (within $\pm$3%), and thus these conditions were not examined in larger systems.

2.4 Optimization of the SPME

The SPME performance was optimized by tweaking the SPME run parameters[32] for careful comparison of the ZMM and the SPME. Fundamental SPME parameters were taken from the default values in GROMACS. The minimum grid distance is 0.12 nm and the SPME spline order was set to 4. The following parameters were considered in the optimization: (1) the SPME was performed with separate reciprocal processes or without separate reciprocal processes (i.e., the real and reciprocal calculations were performed on all processes), (2) various numbers of real-reciprocal process number ratios (the number of reciprocal processes were searched within 15%–50% of the real-space cores), (3) runs were performed with and without OpenMP parallelization with two threads, (4) the compensation between the real-space cutoff distance and the reciprocal grid size (up to 20% increase in the cutoff distance and the grid size) was considered, and (5) the neighbor list was updated every 10, 20, or 40 steps. Some of the parameters were optimized using a tool gmx tune_pme in GROMACS. The fastest combination of (1)-(4) was first searched, and then the best neighbor list update interval (5) was searched to find the optimal condition for the SPME. We set the minimum cutoff distance to 1.0 nm, because the Lennard-Jones interaction cutoff distance must be equal to the electrostatic cutoff for the optimal performance in GROMACS 5, and lowering the cutoff below 1.0 nm showed non-negligible stability differences for the case of proteins.[41] We call the SPME calculation with the optimal parameter set SPME-1.0. We also tried the Lennard-Jones Particle Mesh Ewald (LJ-PME)[12] to have a shorter real-space cutoff, but preliminary analysis showed that the combination of the SPME and the LJ-PME with 0.9 nm cutoff was slower than SPME-1.0 (data not shown), and thus we did not include the LJ-PME in the comparison. In addition to the SPME-1.0, we also performed an SPME simulation with 1.2 nm minimum cutoff for comparison with the ZMM, which we call SPME-1.2.

2.5 Materials

We investigated water systems in detail for fundamental considerations and also protein systems in addition.

Water systems

Since water is necessary for biomolecular simulations and since the accuracy of the ZMM in TIP3P water systems was revealed in Ref. [49], detailed investigation of the computational performance of the ZMM in the TIP3P system is a current principal interest. We prepared four systems with TIP3P[30] water boxes containing 10,000, 20,000, 40,000, and 80,000 molecules, respectively, meaning 30,000, 60,000, 120,000, and 240,000 atom systems, respectively.

Protein systems

As an application to more realistic systems, we benchmarked the ZMM calculation of two systems, Trp-cage and $\gamma\delta$TCR proteins. The Trp-cage is a mini-protein consisting of 10-residue amino acids, consisting of 15,596 atoms including water and ions. The $\gamma\delta$TCR is a 206-residue protein dimer (412 residues in total), solvated in water and ions consisting 93,794 atoms. The trp-cage was modeled as a triclinic periodic boundary system while the $\gamma\delta$TCR as a rhombic dodecahedron periodic boundary system.

2.6 Accuracy measurement

Prior to the performance measurement, we have carefully investigated the accuracy of the ZMM for the mainly concerned system, TIP3P water. Although the accuracy of the TIP3P system was revealed in Ref. [49], we verified the accuracies for the specific systems, including the ones with over $20,000$ molecules, prepared in the current performance study. First we measured the relative differences in both electrostatic energies and densities between the SPME and ZMM to confirm the correctness of the implementation. The mean electrostatic energy $\bar{E}_{\text{el,}\mathsf{M}}$ was measured as the mean of $1000$ total electrostatic energies of method M, $E_{\text{el,}\mathsf{M}}$, which were obtained from a $1$ ns NPT MD run at $310$ K and $1$ atm using method M after an equilibration. The relative difference in electrostatic energies is defined as $\left|\left(\bar{E}_{\text{el,ZMM}}-\bar{E}_{\text{el,SPME-1.2}}\right)/\bar{E}_{\text{el,SPME-1.2}}\right|$. The relative difference in the densities was obtained similarly. We confirmed that the relative differences of electrostatic energies and densities between the SPME-1.2 and the ZMM are less than $4.1\times 10^{-4}$ and $1.1\times 10^{-3}$, respectively, in all three ZMM parameters with $r_{c}=1.2$ nm (ZQuad, ZOct, ZQuad-$\alpha 1$) and all system sizes; these results are consistent to the reported deviations in previous literatures[17, 49] and in section 2.1.2. The standard deviations of the electrostatic energies with the SPME and the ZMM are presented in Table S1. In all the three ZMM parameters and the four system sizes, differences in $\bar{E}_{\text{el,}\mathsf{M}}$ between the SPME and the ZMM were within the margin of the standard deviations.

We further studied the accuracy about the electrostatic energy via the differences in the energy between two molecule conformations. This test can effectively capture the deviations among different methods, because it employs the pair of conformations so that the data are larger than that in the ordinary energy comparison employing every single conformation.[13] After an equilibration, NPT simulations at $310$ K and $1$ atm for the TIP3P water systems with different sizes were performed for $1$ ns with SPME-1.2 (this protocol is the same as that in the performance measurement described in section 2.7, except the simulation length), and conformation snapshots in the trajectories were recorded every $10$ ps. The electrostatic energies of these snapshots were then evaluated by SPME-1.2, ZQuad, ZQuad-$\alpha 1$, and ZOct. For all possible combinations of two snapshots, the energy differences were computed by SPME-1.2 and the ZMMs independently, which made $9900\times 3$ total data.

Figure S1 shows the correlation plots between SPME-1.2 and the individual ZMM for the $10,000$ TIP3P molecule system. The linear regression between the two methods shows an excellent agreement, indicated by $R^{2}\geq 0.9997$ for all ZMMs. This result was similar for the larger systems having $20000$, $40000$, $80000$ TIP3P molecules, as shown by the regression coefficients and correlation coefficients that are also close to unity (Table S2). These results indicate that the all ZMMs with the different parameters were equally accurate for the water systems used in the performance measurement, which is the main issue of this paper, as detailed below. For further investigation of accuracy for other physical quantities and for other materials, readers can consult the literature stated in section 2.1.2.

2.7 Performance Measurement

The performance of the ZMM was measured and compared with that of the SPME in the water systems and protein systems. The benchmark was taken on the Cray XC30 supercomputer at Kyoto University. Each node was composed of two Xeon E5-2695v3 CPUs having 14$\times$2 physical processor cores, and nodes were connected by the Aries interconnect. Benchmarks were taken with 1, 2, 4, 8, 12, and 16 nodes (28, 56, 112, 224, 336, and 448 cores respectively). Throughout the simulations, hyperthreading was not used. The program was compiled with GCC 4.9.2 and FFTW-3.3.3[15] was used as the FFT library. The simulation conditions were set up as follows: the time step for the simulation was set to 2.5 fs and hydrogens connected to heavy atoms were constrained to fixed lengths by LINCS.[26, 25] The SETTLE algorithm was used to simulate rigid TIP3P water molecules efficiently.[33] The systems were equilibrated under SPME-1.2 conditions for 2 ns NPT runs. Simulations were performed with the NPT ensemble using the Parrinello-Rahman barostat and Langevin dynamics. The system was kept at 310 K and 1 atm. The simulation was run for 200 ps to measure the performance.

3 Results and Discussion for the Performance

Figure 1 shows a comparison of the performance of the ZMM and SPME. ZQuad outperformed SPME-1.2 under all conditions. This is an expected result as the ZMM with $\alpha=0$ nm*-1* has a simpler form of the function, and thus requires a lower number of floating point operations. Furthermore, reciprocal convolution is not necessary in the ZMM, which also contributes favorably to the performance of the ZMM. ZQuad also showed comparable or better performance compared to SPME-1.0. This is because the performance loss from the longer cutoff length in ZQuad was compensated by the lack of reciprocal convolution part. The performance difference is more prominent for a larger number of cores. This is possibly due to the increased communication time in the SPME. Both in the 3D-FFT calculation and in the data exchange between the real and reciprocal processes, the relative cost for communication increases as the number of cores grows. This is reflected by the fact that in the SPME-1.0 simulation run with 30,000 atoms and 448 cores, 5.3% of the wall-clock time was consumed waiting for the real-space processes to receive the force information from the reciprocal processes, whereas only 0.7% of the wall-clock time was consumed by the SPME-1.0 run on 112 cores with the same system.

The performance advantage of the ZMM was even more prominent for $\ell=3$, as indicated by the ZOct in Fig. 1. ZOct exhibits a smoother potential function near $r=r_{c}$ than ZQuad, and this leads to a smaller energy drift by the neighbor list (see A). Consequently, a shorter neighbor list cutoff buffer length is required, and this is the main source of the performance increase. For example, in the 30,000 atoms case with a neighbor list update interval of 20 steps, the neighbor list distance $r_{\mathrm{nl}}$ of ZOct was 1.294 nm, while that of ZQuad was 1.350 nm. As a result, ZOct exhibited an increased speed over ZQuad of from 4% to 11% due to the lower number of interactions to compute.

ZQuad-$\alpha 1$ was noticeably slower than ZQuad. This is because the computation of $\mathrm{erfc}\left(\alpha r\right)/r$ is costly. The performance of ZQuad-$\alpha$1 was comparable with or slower than that of SPME-1.0 except for cases with 336 and 448 cores in the 10,000 TIP3P system. ZQuad-1.4 also exhibited similar performance to ZQuad-$\alpha 1$. From these results, it is not worth choosing these electrostatic potentials in terms of computation speed compared to SPME-1.0, at least for $\leq 224$ cores under the present benchmark conditions. When the number of cores is larger, these ZMM settings may become beneficial due to the lack of a reciprocal part, as is the case in the 30,000 atoms case. Overall, ZQuad-$\alpha$1 and ZQuad-1.4 exhibited a performance decrease of 12% to 24% and 17% to 25%, respectively, compared with ZQuad.

Judging from these results and the accuracy of the ZMM presented in the “Methods”, the ZMM with $\ell=2,3$, $r_{c}=1.2$ nm, and $\alpha=0$ nm*-1* is promising for fast simulations, particularly in medium- to high-level parallel environments.

We further analyzed the time required to process an atom to investigate the difference between the ZMM and the SPME. Figure 2 shows the time to process an atom in systems of different sizes. The average time $s$ to process a single atom in a single CPU core is given by the following equation.

[TABLE]

For all methods, the time per step per atom using 28 cores shown in Fig. 2 (a) remains almost constant irrespective of the number of atoms. This is consistent with the time complexity of the ZMM, $O\left(N\right)$, as explained in section 2.1. For the SPME, while the time to process each atom is expected to scale as $O\left(N\log N\right)$,[12] this was not observed in the present results. With the present range of $N\leq 240,000$ and 28 CPU cores, the reciprocal calculation part was only ${}<20$% of the total calculation time, and was thus dominated by other calculations in the MD simulation. In the case of using 448 cores in Fig. 2 (b), the time to process each atom decreased as the number of atoms increased. The increase in the number of atoms under a fixed number of cores indicates a decrease in inter-node communications relative to intra-node communications, if the latter cost is not significant. This is the case for the 448 cores. However, this is not the case for the case of 28 cores shown in Fig. 2 (a), where all communication is intra-node such that the time to process each atom is almost constant. This is also related to the fact that the average time when using 448 cores is longer than that when using 28 cores. In fact, the inter-node communication between processes for 448 cores is much slower than intra-node communications for 28 cores.

The increase in the average time for many cores is relatively significant for the SPME. In particular, the SPME-1.2 for 448 cores takes considerably longer time per step than the ZMM does. This also confirms the increased communication time for the reciprocal-space part in parallel environments. This point can be clarified by a direct comparison between SPME-1.2 and ZQuad-$\alpha$, since the real-space parts of SPME-1.2 and ZQuad-$\alpha$1 have comparable computational costs. For example, in the case of 30,000 atoms and 28 cores, SPME-1.2 required 2.44 $\mu$s to process an atom, only 8% slower than ZQuad-$\alpha$1 which took 2.26 $\mu$s on average. In contrast, in the case of 448 cores, the difference between the two increased to 34% (6.46 $\mu$s and 4.84 $\mu$s, respectively).

Figure 2 also indicates, as was observed in Fig. 1, that ZOct exhibited the highest performance among the methods considered, and that ZQuad-1.4 exhibited performance similar to that of ZQuad-$\alpha$1.

Figure 3 shows the performance of the ZMM and the SPME for the Trp-cage and $\gamma\delta$TCR protein systems. The performance characteristics of the ZMM compared to the SPME in the water systems are retained even with realistic setups, including inhomogeneous systems (both Figs. 3 (a) and 3 (b)) and the rhombic dodecahedral periodic boundary conditions (Fig. 3 (b)). In particular, the overall behavior and relative differences among the ZMM runs with different parameters and the SPME are similar for the 30,000-atom TIP3P system (Fig. 1 (a)) and the 15,000-atom Trp-cage system (Fig. 3 (a)), and are also similar for the 240,000-atom TIP3P system (Fig. 1 (b)) and the 100,000-atom $\gamma\delta$TCR system (Fig. 3 (b)). These results imply that the performance of the ZMM is mainly governed by the system size and is unrelated to the system details that are typical to bimolecular simulations. We thus conclude that the ZMM is also beneficial in terms of performance for real-world applications, including protein MD simulations.

As in the case of the SPME, there are several important parameters in the ZMM that are tunable for better performance, while they are out of the scope of the current study. First, testing the accuracy with respect to the neighbor list cutoff distance may be beneficial for pursuing further performance. As we have shown in the comparison between the ZQuad and ZOct methods, the performance of the ZMM was strongly affected by the neighbor list distance. The total drift energy is currently estimated from the sum of absolute values of Eq. (11) to give a conservative evaluation. In contrast to the Lennard-Jones interactions, electrostatic interactions may cancel each other out, and the drift energy may thus be overestimated. The development of a better drift energy estimation method may contribute to better performance of both the ZMM and the SPME.

Second, the ZMM with decreased cutoff lengths below $1.2$ nm, which were not examined in the current study, becomes an interesting topic for exploration. In this regard, the trade-off between the accuracy and the speed should be evaluated. The importance of such evaluations has been revealed in recent studies. In fact, Nozawa et al.[35] discussed deviations between the Ewald method and a family of periodic reaction field methods near the phase transition point in a low-charge-density polymer system. Also for the Ewald method, Hub et al.[28] investigated artifacts that appear in a potential of mean force calculation for a non-zero net charge membrane/water system under an application of the Ewald method. These studies suggest that the accuracy investigations also for the ZMM under various conditions including that away from equilibrium should be important.

4 Conclusion

We implemented the ZMM, a method for calculating electrostatic interactions, in the GROMACS molecular dynamics software package and compared its performance with the SPME. We showed that even though a larger cutoff length is required in the ZMM than in the SPME, the ZMM provides performance comparable to or better than that of the SPME. In the case of medium to high parallelism in particular, the ZMM is a strong competitive alternate to the SPME in terms of performance.

5 Acknowledgement

The authors would like to thank Dr. Kota Kasahara for providing the Trp-cage structure. Some of the computations were performed at the Supercomputer Center in Kyoto University, Kyoto, Japan, through the HPCI research project (project ID: hp150065, hp160165), and Research Center for Computational Science, Okazaki, Japan. This research was supported by a Grant-in-Aid for Scientific Research on Innovative Areas (No. JP221S0002) from the Japanese Ministry of Education, Culture, Sports, Science and Technology, and by the “Development of core technologies for innovative drug development based upon IT” from Japan Agency for Medical Research and Development (AMED).

Appendix A Verlet Buffer Distances

Given a potential function $V_{ij}\left(r\right)$, the expected drift $\left\langle\Delta V_{ij}\right\rangle$, caused by atom $j$ violating the neighbor list of atom $i$, is approximated as follows.[2]

[TABLE]

where $\rho_{j}$ is the number density of the atom $j$, $\sigma_{ij}$ is the standard deviation of the distance between freely moving (i.e., gas phase) particles $i$ and $j$ between updates, $V_{ij}^{\left(k\right)}\left(r_{c}\right)$ is the $k$th derivative of the potential function $V_{ij}\left(r\right)$ at $r=r_{c}$, $n_{k}$ is the order of the Taylor expansion, and $G\left(\cdot\right)$ is the normal distribution with zero mean and unit variance. From Eq. (10) to Eq. (11), $V^{(0)}_{ij}(r_{c})=0$ is assumed and the lower bound of $r_{t}$ was changed to negative infinity. We used Eq. (11) to determine an appropriate $r_{\mathrm{nl}}$ distance for both the ZMM and the SPME so that the expected sum of drifts from both the electrostatic and the Lennard-Jones potentials do not exceed the maximum allowed $\left|\left\langle\Delta V_{ij}\right\rangle\right|$. The value of maximum $\left|\left\langle\Delta V_{ij}\right\rangle\right|$ was set to 0.005 kJ / mol / ps / atom in this research (the default value in GROMACS 5). In GROMACS 5, the order of Taylor expansion $n_{k}$ is only up to $n_{k}=2$ for electrostatic interactions. Because the ZMM electrostatic potential satisfies $V_{ij}^{\left(k\right)}\left(r_{c}\right)=0$ for $k\leq\ell$, using $n_{k}=2$ underestimates the absolute value of drift to be zero with any buffer size for $\ell\geq 2$. We thus increased the order of the Taylor expansion for the ZMM to prevent this underestimation. The third order expansion $n_{k}=3$ was used for $\ell=2$, and $n_{k}=4$ was used for $\ell=3$.