{'fraction_non_alphanumeric': 0.03678764172896965, 'fraction_numerical': 0.014160518888943903, 'mean_word_length': 4.645053102291783, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Influence of Electrostatic Truncation on Simulations of Polarizable Systems

Different schemes for the treatment of long-ranged electrostatic interactions will be examined for water simulations using the polarizable fluctuating charge potential. Several different methods are compared, including Ewald sums, potential shifting, spherical truncation and reaction field corrections. For liquid water, properties such as the energy, pressure, dynamics and structure are more sensitive to the treatment of the long-ranged interactions with polarizable than with non-polarizable potentials.


INTRODUCTION
Simulating systems with long-ranged electrostatic interactions using periodic boundary conditions requires a treatment of the interactions beyond the central simulation box, either using Ewald sums, truncation or modification of the potential. The importance of the proper treatment of long-ranged forces has been demonstrated for many systems, including pure water [1][2][3][4][5], hydrophobic aggregation [6], ionic solutions [3], and proteins and peptides [7][8][9][10][11]. Nevertheless, many, if not most, simulations of aqueous systems are done with periodic boundary conditions but without using Ewald. There are two main motivations for avoiding Ewald. First, many seek to avoid the additional computational time that evaluating the Ewald sums requires, although algorithms such as the cell multipole [12] and particle mesh Ewald [13,14] have made Ewald efficient for large systems. Second, many seek to use periodic boundary conditions to avoid edge effects but eliminate Ewald in the hope that having no direct interactions between periodic images better represents infinite dilution for aqueous solutes. As stated by Allen and Tildesley, Ewald methods "will tend to overemphasize the periodic nature of the model fluid [15]." An alternative to Ewald and simple truncation are reaction field methods in which the volume outside the cut-off distance is treated as a dielectric continuum [15,16]. The influence of long-ranged interactions may be even more significant for polarizable systems since the Coulomb or dipole-dipole interactions will couple to the polarizable charges or dipoles. The recent widespread use of polarizable water models, sometimes without using Ewald [17][18][19][20][21][22][23][24][25][26][27], suggests that an examination of the effects of truncation on polarizable systems is necessary. The polarizable fluctuating charge model for water will be used [28] to examine how different truncation schemes and Ewald influence the structure and dynamics of pure water.
Six different simulation methods will be used.
B. Scaling by the complementary error function, S B (r)=erfc(λr). This is simply the real space part of Ewald and λ is set equal to 5/L as is fairly standard [29].
C. Scaling by erfc(λr), with the addition of the Ewald self-term and mean-field approximation for Fourier space term.

D. Scaling by
Loncharich and Brooks used n = 2 in Equation 1 [7]. However, References [3] and [4] found that n = 1 works better than n = 2 for pure water. The present simulations will use n = 1 and r cut =L/2, where L is the box length.
E. Nearest image truncation. Coulombic interactions between two-molecules are included only if the distance between center-of-masses is less than a cut-off distance, taken to be L/2.
F. Reaction field correction to truncation. The Coulomb interaction becomes where ǫ RF is the dielectric constant of the reaction field, r cut is the cut-off radius and r RF is the radius of the reaction field. Hünenberger and van Gunsteren have found that r RF =r cut is optimal [5] and so the present simulations will use r RF =r cut =L/2 and ǫ RF is set equal to 79. With r RF =r cut and ǫ RF much larger than 1, the reaction field interaction becomes a scaling function acting on the Coulomb interaction equal to 1-(3/2)(r/r cut ) +(1/2)(r/r cut ) 3 .
The different methods are plotted on Figure 1, which compares S(r)/r with the bare Coulomb interaction. Notice that the functions S B (r), S D (r) and the reaction field interaction are similar and modify the potential over its whole range. Notice also that at the cut-off, at 10, the Coulomb interaction is 0.1, only 20% lower than its value at hydrogen bonding distances of 2 and is still far from zero. Cut-off lengths are usually in the range of 8 to 10Å.
The periodic boundary conditions can either be applied to atoms individually or by groups. If applied by groups (for example, by molecules or residues), then the entire group is moved together so that the nearest image between group centers (for example, the center of mass) is used. If applied by atoms, then each atom is moved individually. The second method can create large dipoles by splitting molecules across a box length. On the other hand, applying the periodic boundary conditions to groups can suddenly shift entire molecules across the box, resulting in discontinuities. If the length scale of a molecule or residue (the distance from the group center to the outermost atom) is R then the distance between two atoms can suddenly change from |r+R| to |r-R|. These problems will not be present if the potential equals zero at L/2-R. For the methods in which the potential slowly approaches zero at L/2 (methods A-D and F), atom based periodic boundary conditions are used, since this was found to conserve energy better. (Although the differences between atom and group based periodic boundary conditions are not large for the Ewald simulation.) For nearest image truncation with a spherical cut-off, the treatment of nearest image makes a large difference and using atombased periodic boundary conditions leads to great instabilities in the system, so a molecular based method was used.
The effects of the various treatments of the long-range forces on the energy, pressure, dipole moment and pair correlation functions will be examined. In addition, the dependence of dynamical properties in terms of the translational and rotational diffusion constants will be examined. Convergence of these properties will be examined for two different system sizes: 256 water molecules (corresponding to L=19.7Å at a density of 1 g/cm 3 ) and 512 molecules (L=24.8Å).

POTENTIAL MODEL
The fluctuating charge (FQ) model is a polarizable potential model in which the partial charges on atomic sites are treated as variables which respond to changes in their environments [28]. The model gives accurate predictions for liquid state properties. At T=298 and P=1 atm, the dielectric constant is 79. The TIP4P-FQ model uses the geometry of the TIP4P water model and includes Lennard-Jones interactions between oxygen sites and three charge sites: two on the hydrogen atoms and one on the M-site 0.15Å from the oxygen atom [30]. The FQ model has additional interactions between charge sites on the same molecule. The energy is where E LJ is the sum of all Lennard-Jones interactions between oxygen sites, E Coulomb is the sum of all Coulomb interactions between different molecules where Q iα is the charge of atom α on molecule i and E pol is the difference in the molecular energy between the liquid and gas-phase, where E gp is the gas-phase energy,χ 0 α is the Mulliken electronegativity of the isolated atom and J αβ (r iα,iβ ) is the intramolecular interaction. The gas-phase energy is the minimized energy of the isolated molecule. With Ewald, E Coulomb becomes and there are two additional energy terms: the Fourier space term and the self-term, which corrects for including i = j terms in the same box in E F S , where erf(x) is the error function [15]. The screening parameter λ was set equal to 5/L, 256 lattice vectors were used in the Fourier space sum and conducting boundary conditions were used. There is a self energy for the reaction field method also, which is given by [5] Model C consists of the Coulomb energy from Equation 6 and the self term (Equation 8) plus an approximate Fourier space term with which was calculated from the Ewald simulation. These terms are constants for rigid molecules but by coupling to the charges they can influence the dynamics. For 256 water molecules and L=19.7Å, φ F S H =0.0125 kcal/mol/e and φ F S O =-0.0280 kcal/mol/e. The set of charges which minimize Equation 3 are the ground state charges, subject to a charge neutrality constraint on each molecule. Rather that solving for the charges exactly at each time step, the method treats them as dynamical variables, which are propagated in an extended Lagrangian formalism at a low temperature so as to remain near the potential energy minimum [28]. The extended Lagangian method introduces a new complication when using cut-offs. If the cutoffs introduce discontinuities in the potentials (as method E does) then propagating the charge degree-of-freedom becomes more difficult, because the forces on the charges are discontinuous. For the case of method E, in order to prevent the charge degrees-of-freedom from getting too hot and drifting away from the potential energy minimum, the exact set of charges was found every 1000 times steps. In all other cases, the charge temperature remains under 10 Kelvin for a 100 ps simulation. In addition, we found that the discontinuities in spherical truncation caused the system to gradually heat up (at about 3 K/ps) so these simulations were all done with a Nosé-Hoover temperature bath with a mass for the Nosé variable equal to 2.0 kcal/mol psec 2 [31,32]. All other simulations were done at constant E,V,N, except for equilibration (at T,V,N) and as noted. The simulations were done with a 1 fs time step and used SHAKE to enforce bond constraints [33]. The Lennard-Jones interactions were calculated only between the nearest periodic images. This too introduces some discontinuities into the energy and forces and sometimes switching functions are used for the Lennard-Jones interactions. However, for a box length of 20Å, the TIP4P-FQ Lennard-Jones interaction at half the box length is only 1x10 −3 kcal/mol. The data presented in the next section is from four separate 100 ps runs.

RESULTS
Energy. The energies for the different simulations are listed in Table 1. The parentheses give 95% confidence intervals. The Ewald results (A) do not show a system size dependence for the total energy. The individual contributions (E Coulomb , E F S and E self ) should be system size dependent, since they use a value of λ dependent on the box size, but the sum of these three terms should be size independent. The largest component of the Ewald terms is E Coulomb . The other terms E F S and E self make much smaller contributions. However, removing these terms, which would correspond to using a complementary error function shifting function (method B), gives a much different energy, which is higher by 1 kcal/mol for the 256 molecule system. For the larger system with a longer r cut the energy is closer to the Ewald result. If the self-term and a mean-field estimate of E F S (see Equations 10 and 11) are added back in (model C), the results are improved considerably. These results are essentially indistinguishable from the Ewald results.
The use of the shifting function (model D) gives energies similar to shifting by the complementary error function and the energies show a strong dependence on cut-off length. Other studies with non-polarizable water potentials have also found that using the shifted potential, S D (r), leads to an increase in the energy. For the SPC potential with r cut =9.3Å, the energy is 0.4 kcal/mol higher [4] and for the TIPS potential with r cut =8.0Å, the energy is 0.6 kcal/mol higher [3]. For the polarizable model with a similar cut-off distance (for 256 molecules, r cut =9.85Å) the difference in the energy is greater. Nearest image truncation (model E), on the other hand, overestimates the energy by almost a half a kcal/mol and does not improve with system size. For a non-polarizable models, the energy with spherical truncation also does not show much of a dependence on cut-off length [2,34] and overestimates the energy, but only by 0.1 kcal/mol [5]. Once again the differences between Ewald and other treatments is greater for the polarizable model. Another study using spherical truncation with r cut =10.5Å and the SPC-FQ model also finds a lower energy (-11.5 kcal/mol) [27] than the Ewald result (-9.9 kcal/mol) [28]. For the RPOL model of water, which treats polarizability using point inducible dipoles, spherical truncation (with r cut =9Å) only slightly overestimates the energy [35,36]. The cut-offs apparently are more severe for the charge-charge (1/r) interactions than for the dipole-dipole (1/r 3 ) interactions. In fact, the cut-offs introduce less errors for the RPOL model than for non-polarizable models.
The reaction field method (model F) gives very good agreement with the Ewald results for the energy. This agreement, and the improvement over the shifting function (method D), is remarkable considering the similarity of the treatment of the Coulomb interaction (see Figure 1). The reaction field interaction (Equation 2)  is greater than S D (r)/r so the electrostatic interactions are stronger, which leads to the lower energy. In addition, there is the self-term (Equation 9) which acts to slightly increase the magnitude of the charges and also lowers the energy. For non-polarizable water potentials, the energy using the reaction field method also agrees well with the energy using Ewald [5].
Pressure. The pressure is a balance of repulsive and attractive forces and so is sensitive to the treatment of electrostatic interactions ( Table 2). Listed are the total pressure, plus the contributions to the virial from the different interactions. The components do not add up to the total because there is the additional ideal gas part. Only the Ewald method gives the correct pressure. The reaction field method, which gave an accurate energy, does not do well for the pressure. The differences in the pressure between Ewald and other methods are more substantial for the polarizable model than for non-polarizable models [4,5]. Like the energy, the pressure with Ewald is size independent. A previous simulation using Ewald with a non-polarizable potential found that the energy and pressure shows no size dependence for systems of 64 or more molecules [37].
Dipole moment. The dipole moment, which in all cases is enhanced relative to the gas-phase value of 1.85 Debye, correlates very well with the energy. Treatments which give an accurate energies (models C and F) also have the same dipole moment as the Ewald method. Those with a larger dipole moment (model E) have a lower energy and those with a smaller dipole moment (models B and D) have a higher energy. It is the sensitivity of the charges that makes the proper treatment of long-ranged electrostatics more important for polarizable models. The method C, which did not work well, is improved considerably just by adding the constant terms which couple to the charges to give the right dipole moment.
Dynamical properties. Also listed in Table 3 are the translational diffusion con-stant and the NMR relaxation time, τ N M R , which gives the time scale for rotations around the axis connecting the hydrogen atoms [28,40]. Methods which have an accurate energy and dipole moment (C and F) have good transport properties. Methods with a higher energy and a lower dipole moment (B and D) have transport properties which are too fast. The diffusion constant is larger and τ N M R is smaller for these methods. For the model (E) with a lower energy and a larger dipole moment, the transport properties are too slow. For method E, constant temperature dynamics is necessary to avoid heating. Transport properties are sensitive to how the velocity rescaling is done. It is preferable to use constant E,V,N dynamics but the Nosé-Hoover method for constant temperature dynamics can reproduce diffusion properties well [41]. Constant temperature dynamics with the Nosé-Hoover method were run using Ewald as a check and the resulting diffusion constant and τ N M R were identical to the constant E,V,N results. In other studies with non-polarizable water potentials, the diffusion constant was found to be about the same as the Ewald result using spherical truncation [5], a reaction field [5] and the shifting potential [4].
Structure. The radial distribution functions are sensitive to the treatment of the long-ranged interactions. Figure 2 shows the oxygen-oxygen radial distribution function, g OO (r), for Ewald with 512 molecules and the shifted potential (method D) with 256 and 512 molecules. The Ewald g OO (r)'s do not show a size dependence (data not shown) but the shifted potential results do show a size dependence. The first peaks do not have as much structure as Ewald, which is consistent with the smaller charges that result using the shifted potential. Also, there is structure around the cut-off distances (9.85Å and 12.4Å for the different sized systems) due a. Reference [38] b. Reference [39] to truncation effects. For non-polarizable models, the agreement between method D and Ewald is much better, although there still remains structure around the cut-off distance [3,4]. The results using method B are similar to the method D results. They show a similar size dependence, but do not have peaks at the cut-off distances. The interactions are sufficiently modified near the cut-off distance so that there are no truncation effects (see Figure 1). The results using the reaction field method agree well with the Ewald results, except for the peaks near the cut-off distance (Figure 3). The reaction field method also gives peaks at the cut-off distance for non-polarizable potentials [5]. With spherical truncation, the radial distribution functions are similar to the Ewald results and are smooth at the cut-off distances (Figure 4). There is slightly more structure in the first peaks consistent with the larger charges. There is no system size dependence in the two spherical truncation results and the curves are indistinguishable from each other. The g OO (r) using method C is identical to the Ewald result and is not shown.
From the radial distribution functions, the pressure results can be understood. The simulation methods which produce the largest pressure difference (models D and E) are also those which have peaks in the g OO (r) at the cut-off region. These peaks will contribute to the virial and change the pressure.

CONCLUSIONS
For systems with polarizable charges, the long-ranged interactions are coupled to the charge distributions. Therefore, modifying the treatment of these interactions will modify the charges and the Coulomb interaction, q i q j /r, is changed not only by changing 1/r to S(r)/r, but also by changing q i and q j . The results presented here indicate that the properties of liquid water are more sensitive to the treatment of the Coulomb interactions for polarizable systems than for non-polarizable systems. The reaction field correction to truncation (method F) is found to work well for the energy and transport properties, although it does give a large negative pressure. The difference between the reaction field method and shifting the potential by S D (r) (Equation 1), which did not work nearly as well, are small (see Figure 1). Subtle differences in the treatment of the electrostatics can cause large differences in the results. Another truncation method which does not work well, S(r)=erfc(λr), can be made to give results almost identical to the Ewald results just by adding two constant terms which couple to the charges (compare methods B and C). These terms represent the self-term of Ewald (Equation 8) and a mean-field approximation to the Fourier-space term of Ewald (Equations 10 and 11). The fact that the meanfield approximation works so well suggests that fluctuations in the Fourier-space and the forces from this term are not important. Once these terms are added back into model B, the charges become equal to the charges with the Ewald method and all the other properties including the energy, pressure, dynamics and structure are accurately reproduced.