Critical test of isotropic periodic sum techniques with group-based cut-off schemes

Truncation is still chosen for many long-range intermolecular interaction calculations to efficiently compute free-boundary systems, macromolecular systems and net-charge molecular systems, for example. Advanced truncation methods have been developed for long-range intermolecular interactions. Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. The former computes interactions of “atoms” inside the cut-off radius, whereas the latter computes interactions of “molecules” inside the cut-off radius. In this work, the effect of group-based cut-off is investigated for isotropic periodic sum (IPS) techniques, which are promising cut-off treatments to attain advanced accuracy for many types of molecular system. The effect of group-based cut-off is clearly different from that of atom-based cut-off, and severe artefacts are observed in some cases. However, no severe discrepancy from the Ewald sum is observed with the extended IPS techniques.

Molecular dynamics (MD) simulations are expected to give new insights on a molecular level into complex systems, but they incur massive computational costs. The most computationally expensive part of MD is the evaluation of the long-range interactions such as the electrostatic force. A distributed multipoles analysis gives the accurate description for spatial distribution of electric charge within a molecule and thus electrostatic structure of molecular systems 1-3 , but requires the enormous computational cost when considering electrostatic interaction calculations. Thus the electric charge of molecule is ordinary simplified using a few point charges that have constant values. However, although after the simplification, the computational cost of electrostatic interaction is still expensive to perform long-time MD simulations for large number of molecules. This cost can be reduced by employing truncation methods or lattice sum methods. The particle mesh Ewald (PME) sum is the de facto standard lattice sum method 4,5 , in which the long-range interaction is split into real and reciprocal space contributions. The reciprocal space can be calculated efficiently by introducing fast Fourier transforms (FFTs). The PME sum has an advanced computational efficiency, but its accuracy is guaranteed only in charge-neutral systems with periodic boundary conditions. This limitation makes it difficult to simulate some ionic systems without artificial counter ions for charge neutralization [6][7][8][9] . Thus, for net-charge systems and/or free boundary conditions, truncation has been mainly used as another approach [10][11][12][13][14][15][16][17] . Furthermore, a hierarchical tree algorithm that can be used for both free boundary and periodic boundary conditions is conceptually very similar to the truncation method, except for the hierarchical structure for the interaction calculation. Based on this similarity, several methods have been developed that combine the tree algorithm and truncation [18][19][20][21] . Switch/shift functions and Onsager's reaction field method are some examples of typical truncation methods, but these often cause severe artefacts in various systems [22][23][24][25][26][27][28][29][30] . For example, switching function may induce faster decay and oscillations in the velocity autocorrelation function of ions when the switching range is too short (Fig. 2 in ref. 22 ), it is well known that water transport properties are strongly influenced by shift functions (Fig. 5 in ref. 23 , Figs 1 and 5 in ref. 24 and Fig. 2 in ref. 25 ), and reaction field method with smoothing function causes unrealistic dipole-dipole correlations ( Fig. 5 and 9 in ref. 26 ). Therefore, the accuracy of truncation methods should be evaluated carefully. The isotropic periodic sum (IPS) method developed by Wu and Brooks 31 is an advanced truncation method that can effectively reduce the computational cost while maintaining adequate accuracy to estimate many types of molecular system. To date, IPS techniques have been applied to many systems: solids 32 , liquids 31,33-39 , liquid-vapour interfaces [38][39][40][41][42][43] , solid-liquid interfaces 44,45 , phospholipid monolayers 46 , proteins 31,47 , combined quantum-mechanics/molecular-mechanics methods 48 , solid-liquid crystalline phase-transitions 49,50 and constant-pH MD simulations [14][15][16] . Those studies suggest that IPS techniques can be applied to many macromolecular systems and that they give reasonable results. To attain better computational accuracy than that of the original IPS method, a linear-combination-based IPS (LIPS) method 38,39 has been developed. In LIPS, several versions of pair potentials, LIPS-5th 38 , LIPS-SW 39 and so on were designed. Our previous study showed that LIPS-SW is almost as accurate as PME with fine grid spacing (<0.1 nm) when used to estimate solid-liquid crystalline phase-transition temperatures 50 . Hence, the LIPS method is a possible way to extend the capability of the IPS technique.
Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. In atom-based cut-off, only the non-bonded interactions of atoms within a given cut-off length are considered. In short, the effects of "atoms" outside the cut-off length are simply truncated. It is easily seen that the scheme could cause appreciable artefacts when considering multipole moments of polar molecules. Moreover, several truncation methods cannot be used with atom-based cut-off because of its theoretical background. For example, Onsager's reaction field method requires that the outside of the cut-off radius be filled by a macroscopic bulk structure that consists of polar molecules. If an atom-based cut-off scheme was used for the reaction field, unexpected charged atoms would appear outside the cut-off radius. Therefore, the reaction field is normally used only with a group-based cut-off scheme. Group-based cut-off is also the standard cut-off scheme, implemented in several molecular simulation packages such as GROMACS [51][52][53][54][55] . In group-based cut-off, non-bonded interactions are truncated by the distance between charge groups that are defined considering charge-neutrality, functional group and residue, for example. Note that the computational accuracy of group-based cut-offs may depend on the definition of groups. The user has to define groups to attain accurate and efficient computations. Considering the charge-neutrality of groups is one of the possible ways. For example, one water molecule is proper as one group in most cases. For biological molecules, the interaction between charge-neutral groups are smoothly and continuously truncated for stabilizing the simulations [56][57][58][59] . From the perspective of molecular chemistry, this treatment is more acceptable than atom-based cut-off. Furthermore, for interaction calculations, group-based cut-off is computationally cheaper for searching interaction distances. In short, there are fewer pairs for groups than for atoms. Given its aforementioned advantages, there are many cases that group-based cut-off is more appropriate than atom-based cut-off for use with complex chemical/biological systems. However, a serious defect has been reported regarding the use of group-based cut-off for molecular systems, albeit very simple polar ones. A strange layered structure has been reported in bulk water systems [28][29][30] . The cause of this defect was investigated systematically, and a severe artefact around the cut-off distance was observed in the Kirkwood factor G k (r) 28,30 . This shows clearly that group-based cut-off causes serious artefacts in dipole-dipole correlations and stabilizes the anomalous layer structure in bulk water systems. This effect was also observed for switch/shift functions and the reaction field method with group-based cut-off 28,30 . Importantly, these defects were clearly reduced when using atom-based cut-off 37 . Thus we reason that group-based cut-off is more prone to erroneous behaviour than atom-based cut-off, and careful investigation is required.
As mentioned above, cut-off schemes are very important for all truncation methods. Because of how it is designed, IPS techniques are commonly used with atom-based cut-off. However, for IPS techniques to be applicable to wider varieties of simulation systems, conditions and packages, their capability with group-based cut-off should be evaluated. For example, a water molecule, which is abundant in many simulations, is treated as one charge group in most cases. Therefore, we carefully estimate herein the accuracy of IPS techniques with group-based cut-off for MD simulations of bulk water and water-vapour interfacial systems.

IPS methods.
The fundamental IPS concept starts with dividing the potential energy U i of a non-bonded interaction for particle i into two parts: interactions within a local region Ω i and long-range interactions outside Ω i : where r j is the position of particle j, r ij is the position vector from particle i to particle j and u(r ij ) is the potential that describes interaction between atoms based on their distance of separation. In typical cut-off methods, φ(r ij , Ω i ) is simply ignored or is applied to the long-range correction term based on continuum approximations. In IPS techniques, this term is expressed by a function of r c that reflects the effect from IPS image particles 31 . To date, several different types of IPS method have been developed, for example the IPS method for non-polar systems (IPSn) 31,47 , the IPS method for polar systems (IPSp) 33 , the LIPS method with a fifth-order cut-off boundary condition (LIPS-5th) 38 and the novel periodic reaction field method (LIPS-SW) 39 . The IPSn 31,47 and IPSp 33 are variations of the original IPS method, developed to calculate non-polar and polar molecular systems, respectively. However, the above two methods have difficulty to estimate homogeneous or heterogeneous polar molecular systems. The LIPS-5th 38 and LIPS-SW 39 are variations of the LIPS method, developed to enhance the capability of the original IPS method. These show advanced accuracy for estimating homogeneous and heterogeneous molecular systems regardless of polarity of molecules. For example, the LIPS-SW shows almost the same accuracy as the PME with fine grid spacing (<0 0.1 nm) when used to estimate solid-liquid crystalline phase transition temperatures 50 . Note that all the above characters of IPS techniques have been discovered by combinational use with the atom-based cut-off, not with the group-based cut-off. Herein, for efficient computation in GROMACS, all the IPS techniques are implemented on "tabulated interaction functions", a basic function of GROMACS to calculate many types of interactions without editing the source code. Therefore, the observed computational costs of the methods were almost equal to each other. Other implementations can make other results on computational efficiency. Note that detailed information about the tabulated interaction functions and numerical formulation of IPS methods are given in the supporting information.
Simulation system and conditions. MD simulations of bulk water and water-vapour interfacial systems were performed to investigate the effect of group-based cut-off scheme for the IPS techniques. The molecular-dynamics simulation package GROMACS 4.5.5 (double-precision version) 55 was used for all simulations in this study.
The bulk system consisted of 6192 water molecules, and the extended simple point charge (SPC/E) model was used for water molecules. Three-dimensional periodic boundary conditions were used to represent the bulk water structure. The bond lengths of the atoms in each water molecule were constrained using the SHAKE method with the relative tolerance 1e-6. The Verlet velocity algorithm was used with a time step of 2 fs, which is the most common time step for water molecule with the constraint method, and known that energy drift is prevented. The simulations were performed under a constant particle-number, volume and temperature (NVT) ensemble with the Nosé-Hoover thermostat technique. The time constant for coupling 0.04 ps is given and the thermostat is updated every time step. The density and temperature of the system were fixed to 0.997 g/cm 3 and 298.15 K, respectively. Non-bonded interactions were truncated with the group-based cut-off scheme. The non-bonded pair list is updated every time step to eliminate effect from the twin-range cut-off method. The cut-off radius of the Lennard-Jones (LJ) potential was set to 1.2664 nm, which is 4.0 in LJ length units for every condition. To treat the electrostatic interactions, the PME summation method, the IPSn method, the IPSp method, the LIPS-5th method and the LIPS-SW method were used. For the PME method, we chose very accurate conditions so that the interpolation was eighth order, namely a grid spacing shorter than 0.1 nm. For the IPS methods, the cut-off radius r c of the electrostatic potential was changed from 1.2 to 2.8 nm in 0.2 nm increments. Note that LJ and/or electrostatic interactions (PME, IPSn, IPSp, LIPS-5th and LIPS-SW) are computed with group-based cut-off scheme in this study. The results for atom-based cut-off scheme have been given in previous researches 38,39 . The system was equilibrated for 1 ns for each condition, and then sampling calculations were performed for 1 ns to calculate equilibrium averages. To observe the dynamic characteristics of the systems, the self-diffusion coefficient D was calculated by the Einstein relation or the Green-Kubo formula as follows: where t is the time, r i (t) is the position of particle i and <...> N denotes the particle average. To investigate the configuration of water, the radial distribution function g(r), the distance dependence of the Kirkwood factor G k (r) and the radial orientation h OO (r) of the radial distribution of dipole ordering were calculated to compare the accuracy of each method. The radial distribution function g(r) and the distance dependence of the Kirkwood factor G k (r) were given by where n i (r) is the number of molecules in the region between r and r + Δr from molecule i, and u i and u j are the normalized dipole moments of molecules i and j, respectively <...> e denotes an ensemble average at the equilibrium state. The radial orientation h OO (r) of the radial distribution of dipole ordering was given by where g OO (r) is the oxygen-oxygen radial distribution function and θ r cos ( ) e is the radial orientation of the dipole, given by For the water-vapour interfacial system, 10,976 SPC/E water molecules were set inside a cubic simulation box of 10.1 nm on each side. The x and y axes were tangential to the interface and the z axis was normal to the interface. The non-bonded LJ interactions were truncated at 2.5328 nm and the cut-off radius was chosen as 5 nm for IPS methods. All other conditions were the same as those described for the bulk water case. The interfacial system had to be observed for a long time, so a production run was performed for 10 ns for each condition and a further 10 ns simulation was run to calculate equilibrium averages. To assess the accuracy for the interfacial system, the density profile ρ z ( ) e and the electrostatic potential profile ψ(z) were calculated with respect to the surface normal (z axis). ψ(z) was given by double integration of the Poisson equation: where ε 0 is the vacuum permittivity and ρ c (z) is the charge density profile for the z direction. ψ(0) indicated the electrostatic potential for vacuum in liquid-vapour systems.

Results
Bulk water. In the following section, two cut-off schemes are described using "-atom" or "-group", for example, IPSn method with atom-based cut-off is labelled as IPSn-atom. Figure 1 shows the potential energies calculated using the PME, IPSn, IPSp, LIPS-5th and LIPS-SW methods. The simulations were performed with several different values of the cut-off radius r c to investigate the effect of r c on the bulk properties. IPSn-group overestimated the potential energy with shorter cut-off and converges to PME at cut-off radii longer than r c = 2.6 nm. By contrast, it has been reported that IPSn-atom underestimates the potential energy for bulk water systems 38,39 .
These facts indicate that the effect of group-based cut-off on IPSn differs completely from that of atom-based cut-off. Compared to other methods, IPSp-group converges more slowly and deviates more from PME at r c = 2.8 nm. This slow convergence has been reported in the case of IPSp-atom 38,39 . LIPS-5th-group and LIPS-SW-group converge to PME at cut-off radii longer than r c = 1.8 nm. These results show that LIPS-5th and LIPS-SW estimate the potential energy successfully despite using group-based cut-off. Figure 2 shows the self-diffusion coefficient D per molecule with different values of cut-off radius for each method. IPSn-group highly overestimates D regardless of the cut-off distance. By contrast, IPSn-atom tends to underestimate D for bulk water systems 38,39 . Corresponding to the results for potential energy described above, the effects of the two cut-off schemes on IPSn clearly differ. Therefore, it is shown that the choice of cut-off scheme strongly affects not only static properties but also dynamic properties. Other IPS techniques give similar values to Figure 1. Potential energies for bulk water system calculated using the PME, IPSn, IPSp, LIPS-5th and LIPS-SW methods with group-based cut-off. IPSn overestimates potential energy with shorter cut-off and converges to PME at r c = 2.6 nm; LIPS-5th and LIPS-SW converge to PME at r c = 1.8 nm. By contrast, IPSp does not converge within 1.2 nm ≤ r c ≤ 2.8 nm.   those of PME and successfully estimate D with adequate accuracy at 1.2 nm ≤ r c ≤ 2.8 nm. This indicates clearly that IPSp, LIPS-5th and LIPS-SW estimate the dynamic characteristics of bulk water systems successfully despite using group-based cut-off. Figure 3 shows the oxygen-oxygen radial distribution function g(r) calculated with the PME, IPSn, IPSp, LIPS-5th and LIPS-SW methods at r c = 2.0 nm. In Fig. 3 IPSn-group deviates more from PME than do the other methods. The deviation is larger than IPSn-atom 38,39 , and it hardly diminishes even when the cut-off radius is increased. By contrast, the IPSp-group, LIPS-5th-group and LIPS-SW-group forms of g(r) for agree well with the PME one. This indicates clearly that IPSp, LIPS-5th and LIPS-SW produce bulk water structures similar to that from PME despite using group-based cut-off. Figure 4 shows the distance dependence of the Kirkwood factor G k (r), which indicates the dipole-dipole correlation of bulk water systems. It is well known that G k (r) is strongly affected by the choice of cut-off treatment. Thus the distance dependence of the Kirkwood factor is useful for evaluating the choice of IPS technique and cut-off scheme on the dipole-dipole correlation. In Fig. 4, G k (r) calculated with IPSn-group, IPSp-group, LIPS-5th-group and LIPS-SW-group are given and compared with that of PME. The IPSn-group G k (r) deviates appreciably from the PME one and there is also a large fluctuation near r c . This highlights the problem of using IPSn-group. The observed deviation is much larger than IPSn-atom, as reported previously 38,39 . The tendency of G k (r) is similar to that of the reaction field method with group-based cut-off 28 , namely that G k (r) decreases sharply near r = r c . This indicates a serious defect in dipole-dipole correlation due to the poor treatment of the cut-off boundary condition for multipoles. By contrast, it is observed that G k (r) calculated by IPSp-group, LIPS-5th-group and LIPS-SW-group is almost the same as that calculated by PME despite using group-based cut-off. The dipole-dipole correlation can be estimated reasonably using IPS techniques with group-based cut-off except when using IPSn.
To observe dipole ordering of water molecules, the radial distribution h OO (r) of the dipole ordering for water molecules was calculated. Figure 5 shows h OO (r) calculated with the PME, IPSn, IPSp, LIPS-5th and LIPS-SW methods for r c = 2.0 nm. The IPSn-group h OO (r) deviates more from the PME one than others. The deviation is larger than that IPSn-atom 38,39 , it hardly diminishes even when the cut-off radius is increased. By contrast, it is observed that h OO (r) calculated by IPSp-group, LIPS-5th-group and LIPS-SW-group is almost the same as that calculated by PME. The dipole ordering of water molecules can be estimated reasonably using IPS techniques with group-based cut-off except when using IPSn.
As shown above, using group-based cut-off strongly affects the accuracy of the IPSn method. One should note that this effect can be avoided by the practical technique for making the tabulated potential, namely the introduction of a pseudo cut-off radius for group-based cut-off. When the pseudo cut-off radius ⁎ r c is set to be larger than the cut-off radius r c for IPSn, the potential at ⁎ < < r r r c c is zero regardless of using group-based cut-off. Thus this treatment is numerically equal to atom-based cut-off with cut-off radius r c (see supporting information).
Water-vapour interface. Figure 6 shows the electrostatic potential profile ψ(z) for the water-vapour interfacial system calculated using the PME, IPSn, IPSp, LIPS-5th and LIPS-SW methods. The IPSn-group, LIPS-5th-group and LIPS-SW-group forms of ψ(z) agree well with the PME one, while that of IPSp-group has a large discrepancy at the water slab. This discrepancy of IPSp-group is almost the same as IPSp-atom 38,39,41 . Interestingly, IPSn-group successfully estimates the physical properties of the water-vapour interfacial system but estimates the dipole-dipole correlation poorly for the bulk water system because of the differing conditions near the cut-off boundary. In our simulation, the cut-off boundaries of the water-vapour interfacial system were on the vapour phase whereas the boundaries of the bulk water system were on the water slab. The cut-off boundary has very Figure 3. Oxygen-oxygen radial distribution function g(r) for bulk water system calculated using PME, IPSn, IPSp, LIPS-5th and LIPS-SW with group-based cut-off at r c = 2.0 nm. The oxygen-oxygen g(r) of all methods except IPSn agree well with that of PME. Inset is an enlarged view of the main graph. . Distance-dependent Kirkwood factor G k (r) for bulk water system calculated using PME, IPSn, IPSp, LIPS-5th and LIPS-SW with group-based cut-off. For IPSn, G k (r) deviates appreciably from PME and fluctuation is observed near r c . The artificial configuration of G k (r) obtained with the IPSn method is not observed for IPSp, LIPS-5th and LIPS-SW. Figure 5. Radial distribution function h OO (r) of dipole ordering for bulk water system calculated using PME, IPSn, IPSp, LIPS-5th and LIPS-SW with group-based cut-off at r c = 2.0 nm. All forms of h OO (r) agree well with the PME one except when using IPSn. Inset is an enlarged view of the main graph. little effect because the vapour phase contains few water molecules. Consequently, the defect of IPSn-group is concealed.

Conclusions
The accuracy of IPS techniques with group-based cut-off was evaluated for bulk-water and watervapour-interfacial systems. The results of using IPSn for the bulk water system highlighted the fact that the effects of atom-and group-based cut-off schemes differ clearly. The potential energy and self-diffusion coefficient were overestimated using IPSn-group but underestimated using IPSn-atom. The liquid structure calculated by IPSn-group had defects that were more serious than those calculated by IPSn-atom. For G k (r), the deviation from the results of PME was much larger than that of IPSn-atom. This indicated a serious defect in dipole-dipole correlation due to the poor treatment of the cut-off boundary condition for multipoles. Importantly, our results showed that the aforementioned problems could be avoided by using other IPS techniques. The IPSp, LIPS-5th, and LIPS-SW methods have advanced cut-off boundary treatments that were quite effective despite using group-based cut-off. Therefore, the above methods should be chosen for homogeneous polar molecular systems, but one must pay attention when using IPS techniques for heterogeneous polar molecular systems. For the watervapour interfacial system, the calculated density and electrostatic potential profiles showed that IPSn-group is reasonably accurate without any improvement of the cut-off boundary conditions. This is because the effect of the cut-off boundary is concealed by the vapour phase, which contains few water molecules. Therefore, to estimate both homogeneous and heterogeneous polar molecular systems reasonably, the LIPS-5th and LIPS-SW perform better than other IPS techniques.
The difference between IPSn and other methods is how they treat the cut-off boundary. Because group-based cut-off requires smooth truncation at the cut-off boundary, the treatment of the higher-order derivatives of the potential energy is more critical than it is with atom-based cut-off. IPSn considers only one condition for the first derivative, therefore it may cause unphysical energy barriers at the cut-off boundary. IPSp avoids these problems by introducing a counter-charge assumption at the cut-off boundary. This makes the second and third-order differentials zero at the cut-off boundary condition, but the assumption causes poor estimation of the physical properties of interfacial water systems. In the extended IPS theory, LIPS-5th and LIPS-SW ensure that higher-order derivatives are zero on the cut-off boundary to remove unphysical energy barriers. Therefore they estimate both bulk and interfacial water systems successfully. This study shows that the aforementioned features of IPS techniques are not changed by group-based cut-off except for IPSn. One should note that the IPSm method 60 , an improved IPSn method that includes multipole interactions, may be able to overcome the above problem of IPSn.
The present work shows the robustness of the LIPS method for estimating multipolar molecular systems accurately. Therefore, LIPS methods have great potential to be combined successfully with the hierarchical tree algorithm. This may improve the previously reported IPS/tree method 20 . Combining LIPS-SW and the tree algorithm or fast multipole method may be particularly useful for computing various systems such as free boundary systems, periodic boundary systems and net-charge molecular systems. Figure 6. Electrostatic potential profile ψ(z) for water-vapour interfacial system calculated using PME, IPSn, IPSn, IPSp, LIPS-5th and LIPS-SW with group-based cut-off. The IPSn, LIPS-5th and LIPS-SW forms of ψ(z) agree well with the PME one, but IPSp underestimates ψ(z) appreciably at the water slab.