Article Truncation Effects of Shift Function Methods in Bulk Water Systems

A reduction of the cost for long-range interaction calculation is essential for large-scale molecular systems that contain a lot of point charges. Cutoff methods are often used to reduce the cost of long-range interaction calculations. Molecular dynamics (MD) simulations can be accelerated by using cutoff methods; however, simple truncation or approximation of long-range interactions often offers serious defects for various systems. For example, thermodynamical properties of polar molecular systems are strongly affected by the treatment of the Coulombic interactions and may lead to unphysical results. To assess the truncation effect of some cutoff methods that are categorized as the shift function method, MD simulations for bulk water systems were performed. The results reflect two main factors, i.e., the treatment of cutoff boundary conditions and the presence/absence of the theoretical background for the long-range approximation.


Introduction
In the calculation of thermodynamic, structural and dynamical properties by molecular dynamics (MD) simulations, the effect of long-range interactions is an important issue.Long-range interactions on the periodic boundary conditions (PBCs) can be calculated using the Ewald sum or cutoff methods.The Ewald sum [1] is the key standard method used in calculations involving long-range interactions with the periodic boundary condition.In this method, the total energy is split into real and reciprocal space contributions.Calculation of the Ewald sum involves three problems, the first being that the reciprocal part is computationally expensive.Particle mesh Ewald (PME) [2,3] reduces computational cost for the reciprocal part by using fast Fourier transform (FFT); however, FFT has problems, becoming a cause of a strong bottle neck in massively parallel computers [4].The second is the inherent periodicity, which can develop artifacts [5][6][7][8][9][10][11][12][13][14].The third is that the thermodynamic limit is unclear [15].Notwithstanding the three problems, the Ewald sum and PME are methods of choice that most appropriately represent the long-range interactions.
Cutoff methods are often used to accelerate long-range interaction calculations.The interactions between molecular pairs only with a distance shorter than a given cutoff length are considered, and effects from more distant pairs are truncated or approximated.The plain cutoff, the cutoff with the switch/shift function and the reaction field (RF) method are some examples of typical cutoff methods for Coulombic systems.Simulations can be accelerated by using cutoff methods, but simple truncation or approximation of long-range interactions have serious defects for various systems.In Lennard-Jones (LJ) fluid systems, long-range interactions do not have a prominent effect on transport properties [16,17], but phase equilibria and interfacial properties change drastically [18][19][20].In water systems, the electrostatic interactions dominate the physical properties, and truncation or continuum approximation may lead to unphysical results.A lot of cutoff methods applied to the Coulombic interaction offer insufficient accuracy [21][22][23][24][25][26][27][28][29][30][31], and all of the results are highly sensitive to the cutoff distance.Cutoff methods are also applied to macromolecular systems [12][13][14][32][33][34][35][36][37][38][39][40][41][42][43], and many results indicate that the aforementioned approximations have difficulties when estimating these systems.On the other hand, advanced cutoff-like methods have been developed to avoid the aforementioned difficulties and to accelerate long-range interaction calculations.Wolf et al. [44] developed a method to calculate electrostatic interactions, which is simpler than the Ewald sum.They took into account charge neutrality in the cutoff sphere and discovered that the electrostatic potential of condensed phases seems to have short range behavior.The modified method developed by Fennel and Gezelter [45] could reproduce some thermodynamic properties of homogeneous systems obtained by the Ewald sum.However, the method can hardly estimate heterogeneous systems [46,47].Wu and Brooks developed the isotropic periodic sum (IPS) method [48][49][50].The IPS method is a method that can calculate contributions from the infinite periodic structure without reciprocal space calculations.Some reports on the accuracy of the IPS method of homogeneous [48,[50][51][52][53][54] and heterogeneous systems [47,50,[55][56][57] show that the method yields estimates in good agreement with the results of the Ewald sum.Improved methods were developed to speed up calculations for large-scale systems [58] and to improve the accuracy for homogeneous and heterogeneous systems [59,60].
Some cutoff methods and the aforementioned two cutoff-like methods can be regarded as the shift function method, which produces the pairwise-potential shifted from the original potential function (e.g., Coulombic interaction) by any theoretical or other requirements.To discover the relation between truncation effects and approximation treatments of long-range interactions, in this work, we focused on shift function methods.MD simulations of bulk water systems were carried out for evaluating the truncation effects of the potential energy, self-diffusion, radial distribution function and the dipole-dipole correlation.The results reflect two main factors, i.e., the treatment of cutoff boundary conditions and the presence/absence of the theoretical background for long-range approximation.

Experimental
MD simulations for bulk water systems were conducted to examine the truncation effects of shifted potentials, and physical properties were compared with those from the simulations of the Ewald sum.For shift function methods, CHARMm-shift [61], Ohmine-shift [62], the dumped-shifted-force potential of the Wolf method (Wolf-DSF) [45], the RFmethod with an infinite dielectric constant (RF-metal), the IPS method for non-polar systems (IPSn) [48], the IPS method for polar systems (IPSp) [50] and the linearcombination-based IPS (LIPS) method with a fifth-order cutoff boundary condition (LIPS-fifth) [59] were chosen.CHARMm-shift is used in CHARMm [61] for shifting Coulombic and LJ interactions.The cutoff boundary conditions are considered until the first-order differential of the interaction potential (first-order cutoff boundary condition).Ohmine-shift is originally one of the switching function methods.The method provides shifted pairwise-potential, if the switching point is set to zero.This potential has second-order cutoff boundary conditions.Wolf-DSF was developed by Fennel and Gezelter [45].The charge neutrality assumption inside the cutoff sphere is the basic concept of the original Wolf method and Wolf-DSF.In the bulk water systems, the α-parameter of Wolf-DSF is set to 0.2 nm −1 [45].It should be noted that better α-parameters for any other systems potentially exist [46,47].RF-metal is the RF method with an infinite dielectric constant.In the RF theory, the Coulombic interaction can be modified for homogeneous systems, by assuming a constant dielectric environment beyond the cutoff sphere.Originally, the dielectric constant of the RF method should be set to realistic value.However, some results indicate that the RF method with an infinite dielectric constant is best for estimating bulk water systems [30].Therefore, we set the dielectric constant of the RF method to the infinite value, like a bulk metal.IPSn and IPSp are two different versions of the IPS method.IPSn is applied to calculations for point charges, whereas IPSp calculates polar molecules.The IPS method assumes the isotropic periodic structure outside of the cutoff sphere.The contribution from this structure (periodic reaction field) determines the shape of the IPS potentials.LIPS-fifth is the potential produced by the improved IPS method, called the LIPS method.The LIPS method is based on the extended IPS theory that provides the design procedure of the periodic reaction fields.
In order to clarify the truncation effect of the Coulombic interaction for shift function methods, it was only applied for the Coulombic interaction, and a cutoff method was used for the LJ interaction.The cutoff radius of the LJ interaction was set as 1.2664 nm, which is 4.0 in LJ length units.For the shift function methods, the cutoff radius, r c , of the Coulombic interaction was changed from 1.2 nm to 2.0 nm by 0.2 nm increments.In this simulation, the extended simple point charge (SPC/E) model [63] was used for water molecules.The velocity Verlet algorithm [64] was used with three-dimensional periodic boundary conditions along with a time step of 2 fs.The atoms in a water molecule were constrained by the RATTLEalgorithm [65].The simulation was performed in a constant particle-number, volume and temperature ensemble with the Nosé-Hoover thermostat [66][67][68], where the number of water molecules was 6,192, the density was 0.997 cm 3 and the temperature was 298.15 K.After equilibrating the system, a total of 5 × 10 5 time steps (1 ns) were carried out for each cutoff radius of the shift function methods.The potential energy, U , the self-diffusion coefficient, D, the radial distribution function, g(r), the distance dependence of the Kirkwood factor, G K (r), and the radial distribution of the dipole ordering, s(r), were calculated.We calculated the self-diffusion coefficient for the transport coefficients.The self-diffusion coefficient can be determined either by the Einstein relation or the Green-Kubo formula, which are basically equivalent formulas.Here, we used the Einstein relation: where t is the time, r i (t) is the position of particle i and ⟨• • • ⟩ N denotes the particle average.The slope of the mean-squared displacement of the diffusing particle in the long-time limit is calculated for the diffusion coefficient.The radial distribution function, the distance dependence of the Kirkwood factor and the radial distribution of dipole ordering were calculated for the configuration of water.These properties are given as a function of the distance between two water molecules, denoted r.The conventional expressions give: and: where n i (r) is the number of molecules that exist in the region between r and r + ∆r from molecule i. u i and u j are the normalized dipole moments of molecules i and j, respectively, while ⟨• • • ⟩ e signifies an equilibrium ensemble average.
All of above properties calculated from the shift function methods were compared with that from the Ewald sum.For the Ewald sum, the cutoff radius for the real part was 2.8 nm.The α-parameter was determined by the following equation: where δ is a small number, which indicates the convergence of real space potentials in the Ewald sum.δ was 10 −6 , so αr c = (6 ln 10) 1/2 .

Potential Energy
The thermodynamic properties for the shift function methods and Ewald sum were calculated by potential energies.Figure 1 shows the potential energy per molecule with different cutoff radii, for the shift function methods and the Ewald sum.The results from CHARMm-shift, Ohmine-shift and Wolf-DSF are far from that of the Ewald sum.In contrast, RF-metal, IPSn, IPSp and LIPS-fifth are close to that of the Ewald sum.IPSp had intermediate values between the former and latter.The error of the potential energy for each method decreases by an increment of the cutoff radius, except for the case of the Wolf-DSF.The fastest decline was observed in the case of LIPS-fifth; the error was roughly in proportion to r −4  c .RF-metal and LIPS-fifth at r c = 2.0 nm achieved the smallest error and had the same value as that of the Ewald sum within 0.02%.To examine the cutoff radius tendency of the potential energy thoroughly, we plotted the error of the potential energy calculated with the shift function methods against that determined with the Ewald sum, as shown in Figure 2. The error of the potential energy for each method decreases by an increment of the cutoff radius, except for the case of the Wolf-DSF.The fastest decline was observed in the case of LIPS-fifth; the error was roughly in proportion to r −4  c .RF-metal and LIPS-fifth at r c = 2.0 nm achieved the smallest error and had the same value as that of the Ewald sum within 0.02%.It is clearly shown that CHARMm-shift, Ohmine-shift and Wolf-DSF poorly estimated the potential energy for bulk water systems.The reason for this is related to the presence/absence of the theoretical background for contributions outside the cutoff sphere.CHARMm-shift and Ohmine-shift do not have the theory that justifies their shifting procedure.Wolf-DSF explains its own truncation treatment by the charge neutrality assumption inside the cutoff sphere, but contributions from outside are not considered.In contrast, RF-metal, IPSn and LIPS-fifth are, respectively, based on a definite theory that considers the contributions from long-range interactions.Therefore, RF-metal, IPSn and LIPS-fifth have much better accuracy for estimating the potential energy.IPSp had intermediate values between the former and latter.This seems to be related to the counter-charge assumption of the IPSp.The counter-charge effect assumed at the cutoff boundary may partially interrupt long-range contributions.

Self-Diffusion Coefficient
We calculated the self-diffusion coefficient for the Figure 3 shows the self-diffusion coefficient for shift function methods and the Ewald sum.The results from CHARMm-shift could not have a similar value to that of the Ewald sum at 1.2 nm ≤ r c ≤ 2.0 nm.Other methods seem to estimate the selfdiffusion coefficient with an adequate accuracy.To examine the cutoff radius tendency of the self-diffusion coefficient thoroughly, we plotted the error of the self-diffusion coefficient calculated with the shift function methods against that determined with the Ewald sum, as shown in Figure 4.The convergence of the IPSp and LIPS-fifth is much faster than other methods.For IPSp, the self-diffusion coefficient is saturated at r c ≥ 1.6 nm, and the saturated value is almost the same as that of the Ewald sum (within 0.35%).For LIPS-fifth, the self-diffusion coefficient is saturated at r c ≥ 1.4 nm, and the saturated value is almost the same as that of the Ewald sum (within 0.36%).The difference of the cutoff radius tendency comes from the treatment of the cutoff boundary conditions and long-range interactions.The results show that the improvement of the cutoff boundary conditions or long-range interaction treatment strongly affects the accuracy of the self-diffusion coefficient.In CHARMm-shift, both treatments are insufficient.It has a first-order cutoff boundary condition and does not have any theoretical background for long-range interaction treatment.Ohmine-shift had improved accuracy in comparison with that of CHARMm-shift, even if it merely comes from an advantage on the cutoff boundary condition.RF-metal and IPSn consider the first-order cutoff boundary condition and the adequate treatment for long-range contributions.Therefore, these two methods have similar accuracy to that of Ohmine-shift.Wolf-DSF also had similar accuracy to the result of Ohmine-shift, despite the absence of the theoretical background for the long-range interaction treatment.Strictly, Wolf-DSF has a first-order cutoff boundary condition, but it can be regarded as an infinite-order cutoff boundary condition under the certain value of alpha.This is the reason for the results of Wolf-DSF.In IPSp, a faster convergence of errors were observed, because it has a third-order cutoff boundary condition and the long-range interaction treatment.LIPS-fifth achieved the fastest convergence.It has a fifth-order cutoff boundary condition and a reliable background for the long-range interaction treatment.Figure 4.The error of the self-diffusion coefficient calculated with the shift function methods against that determined with the Ewald sum.The convergence of the IPSp and LIPS-fifth is much faster than other methods.For IPSp, the self-diffusion coefficient is saturated at r c ≥ 1.6 nm and the saturated value is almost the same as that of the Ewald sum (within 0.35%).For LIPS-fifth, the self-diffusion coefficient is saturated at r c ≥ 1.4 nm, and the saturated value is almost the same as that of the Ewald sum (within 0.36%).

Radial Distribution Function
To examine the structure around a molecule for shift function methods, the radial distribution function, g(r), was calculated.Figure 5 shows the oxygen-oxygen, g(r), of the water molecule for shift function methods at r c = 2.0 nm and for the Ewald sum.In Figure 5, CHARMm-shift, Ohmine-shift, RF-metal and IPSn have notable deviations from the result of the Ewald sum.On the other hand, Wolf-DSF, IPSp and LIPS-fifth provided adequate accuracy.The oxygen-hydrogen and hydrogen-hydrogen, g(r), have very similar behavior in comparison to oxygen-oxygen in Figure 5.To examine the decrease of the deviation for r c thoroughly, we plotted the root mean square deviation (RMSD) of the oxygen-oxygen, g(r), for each shift function method against the Ewald sum at different cutoff radii in Figure 6(a).The RMSDs of the oxygen-hydrogen and hydrogen-hydrogen, g(r), are also plotted in Figure 6(b) and (c), respectively.The difference of the cutoff radius tendency is affected strongly by the treatment of the cutoff boundary conditions.In CHARMm-shift, RF-metal and IPSn that have a first-order cutoff boundary condition, the deviation decreases roughly in proportion to r −2.5 c , r −3 c and r −3 c , respectively.Ohmine-shift has a second-order cutoff boundary condition, and the RMSD of g(r) declines roughly in proportion to r −4  c .The RMSD of LIPS-fifth has a similar tendency with these shift function methods for cutoff radii, but a faster decline is observed.The RMSD of LIPS-fifth decreases roughly in proportion to r −6  c .Furthermore, LIPS-fifth gives accurate estimations of g(r); the RMSD converges at r c ≥ 2.0 nm.Converged values of RMSD for LIPS-fifth are most accurate.On the other hand, the RMSDs of Wolf-DSF and IPSp have an adequate accuracy in any cutoff radius.The charge neutrality and counter-charge assumptions of Wolf-DSF and IPSp, respectively, seem to work better for bulk water systems.Figure 6.The RMSDs of (a) the oxygen-oxygen, (b) the oxygen-hydrogen and (c) the hydrogen-hydrogen radial distribution function for the shift function method against the Ewald sum at different cutoff radii.In CHARMm-shift, RF-metal and IPSn that have a first-order cutoff boundary condition, the deviation decreases roughly in proportion to r −2.5 c , r −3 c and r −3 c , respectively.Ohmine-shift has a second-order cutoff boundary condition, and the RMSD of g(r) declines roughly in proportion to r −4  c .The RMSD of LIPS-fifth has a similar tendency with these shift function methods for cutoff radii, but a faster decline is observed.The RMSD of LIPS-fifth decreases roughly in proportion to r −6  c .Furthermore, LIPS-fifth gives accurate estimations of g(r); the RMSD converges at r c ≥ 2.0 nm.Converged values of RMSD for LIPS-fifth are most accurate.The RMSDs of Wolf-DSF and IPSp have an adequate accuracy in any cutoff radius.We focused on the distance dependence of the Kirkwood factor, G K (r), where one can see the dipole-dipole correlation of bulk water systems.G K (r) has a strong cutoff radius effect, and the influence of the interaction treatment is quantitatively-expressible by the shape of G K (r).An evident shortcoming of the cutoff-like method appears for the G K (r) value in bulk water systems.Thus, G K (r) of various cutoff radii were calculated using the shift function methods to evaluate the truncation effect of the dipole-dipole correlation.
Figure 7 shows the shape of G K (r) determined using the shift function methods and the Ewald sum.It is clearly seen that G K (r) calculated with CHARMm-shift, Ohmine-shift, RF-metal and IPSn fluctuate near r c as in g(r), and this fluctuation still remains in spite of the increment of the cutoff radius.The artificial configuration of Ohmine-shift was smaller than that of the other three methods.The defect of G K (r) for these above shift function methods was not seen in Wolf-DSF, IPSp and LIPS-fifth.Wolf-DSF, IPSp and LIPS-fifth can estimate G K (r) more adequately than other shift function methods.The result of s(r) shows the radial distribution of the dipole ordering for water molecules.Figure 8 presents s(r) calculated with the shift function methods at r c = 2.0 nm, along with that determined by the Ewald sum for comparison.s(r) calculated with the CHARMm-shift, Ohmine-shift, RF-metal and IPSn fluctuate near r c , like g(r), despite the long cutoff radius.Wolf-DSF, IPSp and LIPS-fifth did not calculate any singular configurations of s(r), like that for the other shift function methods.These three methods can estimate s(r) with adequate accuracy.
The aforementioned characteristics of the truncation effect in the dipole-dipole correlation for the shift function methods are affected strongly by the treatment of the cutoff boundary condition, like the case of g(r).

Conclusions
To assess the truncation effect of some cutoff methods that are categorized as the shift function method, MD simulations for bulk water systems were performed.The results reflect mainly two main factors, i.e., the treatment of cutoff boundary conditions and the presence/absence of the theoretical background for long-range approximation.
The difference of estimated value of the potential energy is related to the presence/absence of the theoretical background for contributions outside the cutoff sphere.CHARMm-shift, Ohmine-shift and Wolf-DSF poorly estimated the potential energy, because these methods do not have a reliable theory that justifies their shifting procedure.In contrast, RF-metal, IPSn and LIPS-fifth are, respectively, based on a definite theory, which considers contributions from long-range interactions.RF-metal, IPSn and LIPS-fifth have much better accuracy for estimating the potential energy.The fastest decline was observed in the case of LIPS-fifth; the error was roughly in proportion to r −4  c .RF-metal and LIPS-fifth at r c = 2.0 nm achieved the smallest error and had the same value as that of the Ewald sum within 0.02%.
For estimating the self-diffusion coefficient, the difference of the cutoff radius tendency comes from the treatment of cutoff boundary conditions and long-range interactions.In IPSp, a faster convergence of errors was observed, because it has a third-order cutoff boundary condition and the long-range interaction treatment.For IPSp, the self-diffusion coefficient is saturated at r c ≥ 1.6 nm, and the saturated value is almost the same as that of the Ewald sum (within 0.35%).LIPS-fifth achieved the fastest convergence in this work.It has a fifth-order cutoff boundary condition and reliable treatment for long-range interactions.For LIPS-fifth, the self-diffusion coefficient is saturated at r c ≥ 1.4 nm, and the saturated value is almost the same as that of the Ewald sum (within 0.36%).
The truncation effect in the radial distribution function mainly reflects the treatment of cutoff boundary conditions.In CHARMm-shift, RF-metal and IPSn, which have a first-order cutoff boundary condition, the deviation decreases roughly in proportion to r −2.5 c , r −3 c and r −3 c , respectively.Ohmine-shift has a second-order cutoff boundary condition, and the RMSD of g(r) declines roughly in proportion to r −4  c .The RMSD of LIPS-fifth has a similar tendency with these shift function methods for cutoff radii, but a faster decline is observed.The RMSD of LIPS-fifth decreases roughly in proportion to r −6  c .Furthermore, LIPS-fifth gives accurate estimations of g(r); the RMSD converges at r c = 2.0 nm.Converged values of RMSD for LIPS-fifth are most accurate.On the other hand, the RMSDs of Wolf-DSF and IPSp have an adequate accuracy in any cutoff radius.The charge neutrality and counter-charge assumptions of Wolf-DSF and IPSp, respectively, seem to work better for bulk water systems.
The cutoff radius effect in the dipole-dipole correlation is very similar to that of g(r).G K (r), and s(r) calculated with CHARMm-shift, Ohmine-shift, RF-metal and IPSn fluctuate near r c , as in g(r).The artificial configuration of Ohmine-shift was smaller than that of the other three methods.The defect of G K (r) and s(r) for these above shift function methods was not seen in Wolf-DSF, IPSp and LIPS-fifth.Wolf-DSF, IPSp and LIPS-fifth can estimate the dipole-dipole correlation more adequately than other shift function methods.
Overall, the shift function method that has a higher-order cutoff boundary condition and a reliable theoretical background for long-range interaction treatments achieves better accuracy.For estimating the potential energy and self-diffusion, LIPS-fifth is the most accurate shift function method.In the estimation of the radial distribution function, Wolf-DSF and IPSp have good accuracy with relatively short cutoff radii, and LIPS-fifth becomes most accurate at r c = 2.0 nm.

Figure 1 .
Figure 1.Potential energy for the shift function methods and the Ewald sum.The results from CHARMm-shift, Ohmine-shift and the Wolf method (Wolf-DSF) are far from that of the Ewald sum.In contrast, RF-metal, the isotropic periodic sum for non-polar systems (IPSn), for polar systems (IPSp) and the linear-combination-based IPS (LIPS)-fifth are close to that of the Ewald sum.

Figure 2 .
Figure2.The error of the potential energy calculated with the shift function methods against that determined with the Ewald sum.It is clearly shown that CHARMm-shift, Ohmine-shift and Wolf-DSF poorly estimated the potential energy for bulk water systems.In contrast, RF-metal, IPSn and LIPS-fifth have much better accuracy for estimating the potential energy.IPSp had intermediate values between the former and latter.The error of the potential energy for each method decreases by an increment of the cutoff radius, except for the case of the Wolf-DSF.The fastest decline was observed in the case of LIPS-fifth; the error was roughly in proportion to r−4  c .RF-metal and LIPS-fifth at r c = 2.0 nm achieved the smallest error and had the same value as that of the Ewald sum within 0.02%.

Figure 3 .
Figure 3.The self-diffusion coefficient for the shift function methods and the Ewald sum.The results from CHARMm-shift could not have a similar value to that of the Ewald sum at 1.2 nm ≤ r c ≤ 2.0 nm.Other methods seem to estimate the self-diffusion coefficient with an adequate accuracy.

Figure 5 .
Figure 5.The oxygen-oxygen radial distribution function of the water molecule for the shift function methods at r c = 2.0 nm and for the Ewald sum.CHARMm-shift, Ohmine-shift, RF-metal and IPSn have notable deviations from the result of the Ewald sum.On the other hand, Wolf-DSF, IPSp and LIPS-fifth provided adequate accuracy.The oxygen-hydrogen and hydrogen-hydrogen, g(r), have very similar behavior in comparison to oxygen-oxygen (figures not shown).

Figure 7 .
Figure 7. Distance dependence of the Kirkwood factor for the shift function methods and the Ewald sum.It is clearly seen that G K (r) calculated with CHARMm-shift, Ohmine-shift, RF-metal and IPSn fluctuate near r c as in g(r), and this fluctuation still remains in spite of the increment of the cutoff radius.The artificial configuration of Ohmine-shift was smaller than that of the other three methods.The defect of G K (r) for these above shift function methods was not seen in Wolf-DSF, IPSp and LIPS-fifth.Wolf-DSF, IPSp and LIPS-fifth can estimate G K (r) more adequately than other shift function methods.

Figure 8 .
Figure 8.Radial distributions of dipole ordering calculated with the shift function methods at r c = 2.0 nm and for the Ewald sum.s(r) calculated with the CHARMm-shift, Ohmine-shift, RF-metal and IPSn fluctuate near r c , like g(r), despite the long cutoff radius.Wolf-DSF, IPSp and LIPS-fifth did not calculate any singular configurations of s(r), like that for the other shift function methods.These three methods can estimate s(r) with adequate accuracy.