Reaction-field electrostatics in molecular dynamics simulations: development of a conservative scheme compatible with an atomic cutoﬀ †

In molecular dynamics (MD) simulations of condensed-phase systems, straight-cutoﬀ truncation of the non-bonded interactions is well known to cause cutoﬀ noise and serious artifacts in many simulated properties. These eﬀects can be drastically reduced by applying the truncation based on distances between neutral charge groups (CG) rather than between individual atoms (AT). In addition, the mean eﬀect of the omitted electrostatic interactions beyond the cutoﬀ distance can be reintroduced using the reaction-field (RF) method, where the medium outside the cutoﬀ sphere is approximated as a dielectric continuum of permittivity equal to that of the solvent. The RF scheme is generally applied with CG truncation. This is justified for low solvent permittivities, where the RF correction is small and an AT truncation would lead to severe issues, just as in the straight-cutoﬀ case. However, it is less appropriate for solvents with high permittivities, where the RF correction acts as a physically motivated shifting function, and a CG truncation may in turn lead to artifacts and poorer energy conservation. In this study, we assess the impact of truncation artifacts considering the 57 organic liquids which were used in the calibration of the GROMOS-compatible 2016H66 force field. Combinations of shifting or switching schemes with RF-based electrostatic interactions as well as van der Waals (Lennard-Jones) interactions are then introduced to resolve the issues with AT truncation. These shifting and switching schemes have the following properties: (i) they bring the force but not the potential energy to zero at the cutoﬀ; (ii) as a result, they lead to a modification of the interaction that is comparatively small; (iii) they permit to conduct rigorously conservative simulations; (iv) the energies can easily be corrected back to the unmodified form, either on the fly or in a post-processing step. The mathematical formalism of these schemes is presented in detail, and their validation is performed using the 57 organic liquids.


Introduction
Classical atomistic simulation and, in particular, molecular dynamics (MD) is nowadays a key component in the investigation of (bio-)chemical processes in the condensed phase. [1][2][3][4][5][6] These simulations rely on a potential-energy function or force field, [7][8][9][10][11] constructed as a sum of covalent terms, accounting for interactions between close neighbors within the same molecule, and non-bonded terms, describing the through-space interactions between more remote neighbors as well as between different molecules. The non-bonded interactions include an electrostatic and a van der Waals component, the latter commonly represented using the Lennard-Jones (LJ) function. 12 For condensed-phase (bio-)molecular systems, the accuracy of a force field depends most critically on the representation of the dihedral-angle torsions [13][14][15][16] and of the non-bonded electrostatic and LJ interactions 11,[17][18][19][20][21] (see Tables 1 and 2 in ref. 20 for an overview of studies investigating the impact of different choices). The calculation of the non-bonded interactions also represents the computational bottleneck in MD simulations, and a trade-off must be established between accuracy and efficiency depending on the system and properties of interest.
Restricting the discussion to pairwise-additive force fields and to simulations carried out under periodic boundary conditions, the non-bonded interactions are typically split into short-and long-range components by application of a cutoff distance R c . A specific calculation scheme is then characterized by the following components: 19,20,[22][23][24] (A) the long-range approximation used for the electrostatic and LJ interactions beyond the cutoff; (B) the application modus of the cutoff splitting, i.e. the entities considered for the truncation; (C) the pairlisting procedure for the short-range interaction; (D) the short-range smoothing possibly applied to the interaction; (E) the force-field parameters used for both electrostatic and LJ interactions. These components are discussed in turn below.
For both electrostatic and LJ interactions, there are three main options for the long-range approximation (A), namely: (i) complete omission, leading to a straight-cutoff scheme; (ii) approximation as a continuous medium, leading to a mean-field scheme; (iii) approximation as a periodic medium, leading to a lattice-sum scheme.
For the cutoff distances typically used in simulations (ranging between 0.9 and 1.4 nm), and considering the large magnitude and long-range nature of the electrostatic interactions, abrupt truncation of these interactions in a straight-cutoff approach (Fig. 1a) leads to large errors in many simulated properties. 20,[25][26][27][28] For the LJ interactions, the straight-cutoff scheme is less problematic. 29 It mainly causes a cutoff-dependent underestimation (relative to the long-cutoff limit) of the densities, vaporization enthalpies and surface-tension coefficients, 21,30,31 due to the neglect of the predominantly dispersive (i.e. attractive) interactions beyond the cutoff. The straight-cutoff approach nevertheless remains the most common scheme for LJ interactions nowadays, as it is consistent with the calibration conditions of many popular force fields. This situation might, however, change in the near future.
In the lattice-sum approximation, an assumption of exact periodicity is made for the system simulated under periodic boundary conditions, which is considered to be an infinite pseudo-crystal. By convoluting point charges with a chargeshaping function, the calculation is split into a real-and a reciprocal-space component. 69 The range of the former component is made shorter than the cutoff, so that it can be calculated exactly 70 (or nearly exactly in the case of a Gaussian chargeshaping function) by pairwise summation. The latter component is long-ranged but periodic and smooth in real-space, so that it converges rapidly in reciprocal-space. For the electrostatic interactions, it can be evaluated either by direct summation over reciprocal-lattice vectors, as in the Ewald 71,72 method, or using grid-based fast Fourier transforms, as in the particleparticle-particle-mesh (P3M) 73,74 and particle-mesh Ewald (PME) [75][76][77][78] methods. Both approaches exist as well for the LJ interactions, namely LJ Ewald,30,72,[79][80][81][82] LJ P3M, 83,84 and LJ PME. 76,85,86 Note that these lattice-sum LJ schemes may be computationally expensive unless a strict geometric-mean combination rule 87 applies to the LJ parameters, 30,80 or such a rule is used as an approximation during the calculation. 86 The application of a lattice-sum approach is most appropriate for truly periodic systems (crystals, neglecting static and dynamic disorder), but only approximate for inherently non-periodic systems such as liquids and solutions. In practice, for the comparatively short-ranged LJ interactions, lattice-sum and mean-field approximations typically lead to comparably accurate results for bulk systems, 21,78 but generally not for interfacial systems, where lattice-sum is to be preferred. 30,56,83,85,86,88 Besides reducing the effect of artificial periodicity, there are three additional advantages to the use of mean-field over latticesum schemes for both electrostatic and LJ interactions: (i) for large systems, they enable an O(N) scaling of the computational Fig. 1 Schematic illustration of different treatments of the pairwise electrostatic interactions (similar considerations apply to LJ interactions as well). The cutoff distance is R c . The unmodified (physical) interaction function u(r) is shown for reference in all panels (blue dotted). (a) u(r) truncated with a straight-cutoff scheme (blue solid). (b) u(r) modified using a SH (orange) or a SW (pink) scheme. (c) u(r) corresponding to the RF scheme in its usual implementation for e RF = 5. (d) u(r) corresponding to the RF scheme for e RF = 5 but without off-setting the energy to zero at R c (dotted green line shifted up by a constant to the dashed green line), together with the proposed RF(+LJ)/AT/SH (orange) or RF(+LJ)/AT/SW (pink) modifications.
cost, compared to an O(N log N) scaling for FFT-based latticesum schemes; 74,89 (ii) they enable a pairwise (or groupwise) partitioning of the interaction energy without additional computational cost, whereas such a direct partitioning is not possible in lattice-sum schemes (i.e. it would require additional FFTs 69 ); (iii) for processes involving net-charge changes, the application of lattice-sum schemes is affected by intricate methodological issues. 90 A second characteristic of the non-bonded interaction scheme is the application modus (B) of the cutoff splitting, which may be based on distances between: (i) individual atoms, resulting in an atomic (AT) cutoff scheme; (ii) atom groups, resulting in a charge-group (CG) cutoff scheme. Here, chargegroups are small groups of close covalent neighbors within a molecule, for which the partial charges should add up to an integer value. Historically, charge groups truncation was found necessary to permit the simulation of polar liquids using a straight-cutoff scheme for the electrostatic interactions. [91][92][93][94]  truncation. In addition, the use of a CG-based short-range pairlist 95 represented a simple way to reduce memory requirement and computational cost of the pairlist generation. Note, however, that a well-implemented AT scheme relying on a group-based extended pairlist can reach a very comparable level of computational effort and memory usage. In practice, there are several reasons why modern force fields should probably favor AT over CG truncation: (i) if the potentialenergy function of a mean-field scheme does not exactly vanish at the cutoff, CG truncation leads to loss of energy conservation and heating, 91,92 as well as other artifacts; 28,52,[96][97][98] (ii) in the context of lattice-sum and of some mean-field schemes, the short-range potential-energy function actually vanishes at the cutoff, and an AT truncation then becomes the natural choice; (iii) the definition of neutral charge-groups within molecules generally requires charge adjustments that are rather arbitrary as well as difficult to automate in the context of force-field design. 99 The third characteristic of the non-bonded interaction scheme concerns the pairlisting procedure (C), which may rely on: (i) a pairwise summation that is exact at each simulation timestep; (ii) a summation dictated by an approximate pairlist. Here, common time-saving techniques include Verlet 95 and twin-range 100,101 schemes, where the pairlist is only recalculated at periodic intervals. These variants are in principle designed to increase efficiency while negligibly affecting the simulated properties. However, a number of cases have been reported where the implementation of the method and/or the associated parameters may cause serious artifacts in simulations. 21,102,103 The fourth important aspect concerns the possible shortrange smoothing of the interaction 22,23,96 (D), where the options are: (i) no smoothing applied; (ii) application of a shifting (SH) function; (iii) application of a switching (SW) function. The SH approach 22,23,[104][105][106][107][108][109] refers to the modification of the short-range interaction over the entire range from 0 to R c , so as to enforce specific continuity conditions at R c (e.g. vanishing potentialenergy function along with a number of its derivatives). The SW approach 104,107,108,110,111 refers to a similar modification applied only over a limited range [R c À d, R c ] with 0 o d o R c , and enforcing continuity conditions at R c À d as well. Switching has the advantage over shifting that the very-short-range interactions (below R c À d) remain unmodified, but the drawback that it may induce large forces in the range [R c À d, R c ] if d is chosen small. Historically, SH and SW schemes (with either CG or AT truncation) were applied for the electrostatic interactions to reduce noise and artifacts introduced by straight-cutoff truncation (Fig. 1b). However, the use of such ad hoc schemes without physical motivation may cause major disruptions in the simulated structural, thermodynamic and dynamical properties. 96,108,[111][112][113][114][115] The SW and SH approaches have also been applied to electrostatic interactions in complement to mean-field schemes, for example in the switched RF approach, 116 or in variants of the Wolf-damping method. 114,[117][118][119][120][121][122] In view of the small LJ forces at typical cutoff distances, the application of SH or SW schemes to the LJ interactions has been less frequently considered, 109,111,123 although it may still be useful if one wishes to enforce a strict energy conservation in simulations. Note that SH and SW schemes are irrelevant in the lattice-sum case, where the cutoff is only a numerical parameter and the level of noise is (nearly) exclusively controlled by the accuracy of the reciprocalspace evaluation.
Finally, the choice of appropriate force-field parameters (E) and the options selected for points A-D above (also including the cutoff distance R c itself, and the possible use of specific LJ combination rules) are inter-dependent. This is because force-field parameters are effective quantities, which partly compensate for approximations made in the functional-form representation. As a result, a force field is only expected to provide accurate results under simulation conditions that are close to the ones selected during its calibration. For example, the use of a mean-field or lattice-sum long-range component is not recommended for a force field that was calibrated with straight-cutoff, both in the context of LJ interactions 21,31,55,85,[124][125][126] and electrostatic interactions. 21,127,128 Although the cutoffdependence of the simulated properties may be reduced, the deviations from experiment are likely to increase, unless an appropriate reparametrization is performed.
The present work focuses on the implementation of one of the best established mean-field approaches for electrostatic interactions, namely the RF method. 28,32,33 In the context of molecular liquids, this scheme is underpinned by a particularly sound theoretical basis, as it rests on the Onsager model for a dipole at the center of a cavity immersed in a dielectric medium. 129 It is also the standard electrostatic scheme compatible with the GROMOS force field. [130][131][132] The RF scheme depends on the permittivity e RF of the dielectric medium outside the cutoff sphere, which is generally taken to be that of the solvent (or pure liquid) considered. In its usual implementation (see eqn (1) in the Theory section), the potential energy vanishes at the cutoff, but this is generally not the case for the force (Fig. 1c, green curve with non-zero derivative at R c ). In the limit e RF -1 (vacuum), the RF scheme becomes like a straight-cutoff scheme (Fig. 1a), except for the introduction of a constant offset in ÀR c À1 that brings the function to zero at R c .
In the limit e RF -N (conducting), the RF scheme becomes similar to a physically motivated SH scheme (Fig. 1b, orange curve), in which the potential energy and its first derivative vanish at R c . The RF scheme is generally applied with CG truncation, which is justified for low solvent permittivities. In this case, the RF correction is small and an AT truncation would lead to severe issues, just as in the straight-cutoff case (e.g. noise in O[R c À1 ] rather than O[R c À3 ] for pairs of neutral charge groups). However, a number of studies 28,98,101,133,134 have suggested that the use of a RF/AT scheme gives similar or improved results compared to the RF/CG scheme for aqueous systems.
The goal of this study is to assess the influence of the cutoff treatment on selected properties of pure liquids, as well as to eliminate the need for charge groups by introducing SH and SW schemes in combination with an AT truncation of the RF-based electrostatic interactions (and of the LJ interactions), which are appropriate for low as well as high values of e RF . This work follows analogous developments in the context of Wolf damping 119,120 and isotropic periodic sum 51 methods. The proposed RF(+LJ)/AT/SH and RF(+LJ)/AT/SW schemes have the following properties: (i) they bring the force but not the potential energy to zero at the cutoff; (ii) as a result, they lead to a modification of the interaction that is comparatively small; (iii) they permit rigorously conservative simulations; (iv) the energies can easily be corrected back to the unmodified form, either on the fly or in a post-processing step. The main idea is to shift/switch the force to zero at R c , and let the potential energy level-off to a constant (Fig. 1d). This results in an energy contribution involving all atom pairs at a distance larger than R c . However, this term is easy to calculate and induces neither a force nor a virial contribution. The benefit of this approach is that the continuity of the potential energy as well as that of its first (and second) derivative at R c can be achieved (conservative scheme), but with a limited perturbation of the original RF forces compared to a corresponding shifting/switching of the potential energy (Fig. 1b).
In this article, the mathematical formalism of these RF(+LJ)/ AT/SH and RF(+LJ)/AT/SW schemes is presented in detail. This is followed by their validation in terms of the densities and vaporization enthalpies of the 57 organic liquids considered in the calibration of the GROMOS-compatible 2016H66 force field. 135 The effects of different schemes on the radial distribution functions (RDFs) and dipole-dipole orientation correlation functions (DOCs) of some model liquids are also analyzed.

RF-based electrostatic interactions with AT or CG cutoff schemes
The Coulomb plus RF potential energy v RF,ij between two atoms i and j with charges q i and q j at a distance r ij = |r j À r i | is commonly implemented in the form 28,32,33,[130][131][132] (1) where e 0 is the permittivity of free space. This expression assumes that the atoms are at the centers of their respective spherical cutoff spheres of radii R c , filled by vacuum and surrounded by a homogeneous dielectric medium of relative permittivity e RF (Onsager model 129 ) that is exempt of free charges (no Debye screening 33 ). For simulations in solution (or pure liquids), the value of e RF is typically selected to be that of the solvent (or liquid). The potential energy v RF,ij and the associated interatomic force f RF,ij are shown in Fig. 2 for different values of e RF , assuming atomic charges of the same sign. Owing to the third term in eqn (1), v RF,ij is exactly zero at the cutoff distance R c , irrespective of the value of e RF (Fig. 2a). One can distinguish two limiting cases.
In the limit e RF -N (e.g. orange curve for e RF = 80, representative for water), the interatomic force f RF,ij nearly vanishes at the cutoff distance R c . Thus, in this highpermittivity regime, the RF interaction plays the role of a nearly perfect shifting function, a favorable situation to avoid cutoff noise and pair-distribution artifacts. However, in this regime, the potential energy rises again shortly beyond R c for charges of the same sign, i.e. the interatomic force becomes unphysically attractive (or repulsive for charges of opposite signs). Such a situation may occur when using CG truncation, 28 because interacting CG pairs that have their centers close to R c may involve some interatomic distances that actually exceed R c . In contrast, this situation is excluded with AT truncation.
In the limit e RF -1 (e.g. green curve for e RF = 2, representative for a non-polar organic solvent), the RF interaction becomes progressively closer to a Coulomb interaction with straight-cutoff truncation, i.e. the RF screening no longer acts like a shifting function. In this low-permittivity regime, an AT truncation scheme is expected to induce serious artifacts in the simulated properties. 28 These can be largely remedied by a CG-based truncation. The force f RF,ij as a function of r ij , which is the negative derivative of the potential energy and corresponds to the force exerted by atom i on atom j projected along the interatomic vector r ij = r j À r i .
The consequences of these features on simulated pairdistribution functions are illustrated in Fig. 3 for four model liquids, in the form of RDFs calculated using either a RF/CG or a RF/AT scheme for the electrostatic interactions. The results obtained using a lattice-sum scheme are shown for comparison. In addition, the corresponding DOCs and distance-dependent Kirkwood factors G(r) are depicted in Fig. 4 for chloroform and water.
For the three low-permittivity systems ( Fig. 3a-c), the use of AT truncation leads to an artificial peak in the RDF at R c , particularly pronounced for the CCl 4 model with high partial charges. The reason is that molecule pairs at a distance close to R c can take advantage of a strong potential-energy decrease when they adopt relative positions and orientations that maximize the number of oppositely-charged atom pairs just inside the cutoff and of like-charged atom pairs just outside the cutoff. In contrast, the RF/CG curves are smooth at R c and very close to the reference lattice-sum curves.
In the high-permittivity case of water (Fig. 3d), the AT scheme results in a small wiggle in the RDF around R c , which is essentially absent for the RF/CG scheme. However, an opposite difference is observed for the DOC of water (Fig. 4b), where the RF/CG scheme shows negative correlations below R c and positive ones above, which is also reflected by a dip in the corresponding G(r). In contrast, the RF/AT curve is very close to the reference lattice-sum one. This effect in the RF/CG case is probably related to the slight like-charge attraction between the hydrogen atoms beyond R c that can occur for molecule pairs in specific orientations with their charge-group centers (oxygen atoms) just below R c . Note that these observations on the differences between RF/CG and RF/AT for water were already reported previously in ref. 28 (see Fig. 2-4 therein; note that the labels for the CG/1 and AT/2 curves are mistakenly inverted in Fig. 2a; see also ref. 101). For CHCl 3 , the effect on the DOC is very small (Fig. 4a), likely due to the much smaller partial charges and dipole moment. The Kirkwood factors with latticesum and RF/CG are in good agreement up to R c and diverge afterwards. The use of an AT cutoff causes a slightly stronger alignment of the dipoles close to the cutoff distance compared to the lattice-sum reference.

Objective of the modification schemes
The observations of the previous section suggest that when applying RF-based electrostatic interactions in their usual implementation, 28,32,33,130-132 a CG truncation is to be preferred in the low-permittivity limit, whereas an AT truncation is more adequate for solvents with high permittivity. This is inelegant as well as problematic for intermediate or mixed situations. To address this issue, we introduce here a scheme, which combines RF-based electrostatic interactions with AT truncation, and is designed to fulfill the following requirements ( Fig. 1d): (i) the pairwise potential-energy function is unaltered at zero interatomic distance (in a limiting sense); (ii) it is continuous at the cutoff; (iii) it is constant (generally non-zero) at and beyond the cutoff; (iv) its first derivative vanishes at the cutoff (i.e. zero force at R c ); (v) its second derivative also vanishes at the cutoff (i.e. continuous force at R c ). A similar change is also explored for the LJ interactions, so as to achieve non-bonded interaction schemes ensuring a rigorous conservation of the total energy.
The proposed AT schemes will be further referred to as modified schemes (mod), to distinguish them from the standard (std) implementation. The modified schemes involve  either an additive polynomial shifting function or an additive spline switching function that levels off the force but no the potential energy to zero at R c . In order to be comparable with the standard energies, a correction can be added to the modified energies, either on the fly or as a post-processing step. These recalculated energies will be referred to as corrected (cor) energies.

Standard and modified schemes
For the ease of notation, the subscript o is introduced here to specify the interaction type, namely ''RF'' for the Coulomb plus reaction-field term, ''LJ6'' for the dispersive C 6 term of the Lennard-Jones interactions, and ''LJ12'' for the corresponding repulsive C 12 term. Defining pairwise topology-dependent prefactors as the standard potential-energy function can be written as where the three influence functions are defined as the RF offset term as and the RF self-term as In eqn (3) and (5), the notation j A PL(i) refers to the pairlist of an atom i, i.e. all the atoms j 4 i that are within R c of i according to the selected CG or AT cutoff criterion. The term u RF (r) in eqn (4) accounts for the first two terms of the RF potential-energy function v RF,ij (r ij ) of eqn (1). The contribution of the third term is gathered into the energy offset U off RF of eqn (5). The motivation for the self-term U slf RF of eqn (6) is discussed in ref. 138 (see Appendix B therein; see also ref. 131). Intuitively, this term may be interpreted as the reversible work required to individually charge the atoms when they are at infinite separation (excluding the Coulombic self-energy but including the RF contribution). Both U off RF and U slf RF are configuration-independent. These two terms only affect the energy of the system and induce neither atomic forces nor a virial contribution. However, it may be essential to include them in free-energy calculations involving alterations of the atomic partial charges and comparisons between media with different permittivities. For simplicity, the handling of excluded atom pairs (usually first and second covalent neighbors) has not been explicitly indicated in the above equations. As discussed elsewhere, 131 the entire LJ interaction but only the Coulomb component in r ij À1 of the RF interaction should be omitted from the terms of eqn (3) involving excluded atom pairs. Turning to the modified schemes, it is convenient to introduce a superscript O standing for the modification type, namely ''SH'' for shifting or ''SW'' for switching. For these modified schemes, eqn (3) is altered to Besides omitting U off RF and U slf RF , the modification involves the addition of a short-range (shr; inside the cutoff) and a longrange (lnr; outside the cutoff) component. These are given by and Here, the modification is specified by the influence function u O o (r) to be added to the original influence function u o (r) of eqn (4), and the notation is introduced for the modified function inside the cutoff. The short-range term will be designed to enforce the desired continuity conditions. The long-range term is required to ensure energy conservation. However, it induces neither a force nor a virial. Because it does not account for the long-range interaction energy in a physically meaningful way, this component will be removed when calculating corrected energies in a postprocessing step.
Since the summation in eqn (9) involves a mere constant, this equation can be converted to The first term involves all particle pairs but is constant throughout a simulation. The second term must be recalculated at each simulation timestep, but only involves the inexpensive summation of a distance-independent quantity over particle pairs within the cutoff. Following the objectives of Section 2.2, the conditions to along with In the switching case, conditions similar to eqn (12) are also imposed to enforce corresponding continuity constraints at the onset of the switching range.
The forms of selected functions u O o (r) satisfying these conditions are detailed in the two following sections for the shifting (O = SH, eqn (14)) and switching (O = SW, eqn (15)) schemes. Note that the indicated equations only specify these functions in the relevant range up to R c , as the contribution of ũ O o (r) beyond R c is already encompassed in the long-range term U lnr,O o . It should also be stressed that the retained functional forms are dictated by considerations related to the above continuity conditions and to computational efficiency, but they are by no means unique.

Shifting function
The SH scheme relies on a polynomial of the form which already satisfies eqn (13)  The exponents m o and n o are related to the distance scale of the modification. When the exponents are high, most of the damping of the force towards zero is operated over a small distance range below the cutoff, i.e. the decrease of the force becomes an increasingly close approximation to a step function. When the exponents are low, the force is shifted more progressively towards zero over the entire distance range from 0 to R c . The choice of even integer exponents is also advantageous in terms of computational efficiency. In this work, only pairs of even integers with n o = m o + 2 are considered. In addition, the condition m o Z 4 is imposed, because the choice m RF = 2 for the RF interaction would involve altering the quadratic term of the potential energy, i.e. changing the effective permittivity e RF . This option is nevertheless briefly discussed in the ESI. † As an illustration, the curves for the modified potential energy c RF,ij ũ SH RF and the associated interatomic force c RF,ij f SH RF are displayed in Fig. 5a and c for two variants with m RF -n RF combinations of 4-6 (purple) and 12-14 (orange), respectively. A more detailed analysis of the amount of change in the influence function introduced by different exponent combinations in the context of RF-based electrostatic interactions is provided in Appendix B. This analysis indicates that m o , the smaller of the two exponents, is the one that predominantly determines the shape of the modified influence function.

Switching function
The SW scheme relies on a fourth-order polynomial operating over the distance range from R c -d o to R c , namely where Y is the Heaviside step function (zero if its argument is negative, one otherwise). This expression can be identified with a cubic spline 139 for the forces, and already satisfies eqn (13)   Here also, these values can be precalculated.
As an illustration, the curves for the modified potential energy c RF,ij ũ SW RF and the associated interatomic force c RF,ij f SW RF are displayed in Fig. 5b and d for two variants with d scl RF = 0.78 (purple) and 0.36 (orange), respectively. Comparing to the SH scheme in Fig. 5a and c, one sees that for a sensible range of exponents m o and n o in SH or distance parameters d scl o in SW, both modification schemes can yield analogous influence functions.

Corrected energies
The potential energy U mod,O in eqn (7), with O = SH or SW, is the one relevant for the dynamics of the system when using the corresponding modification scheme. Since the schemes are applied with AT truncation and the potential energy is continuous (up to its second derivative) at the cutoff, strict energy conservation is ensured. However, the potential energy U std of the standard scheme in eqn (3)  For this reason, it is useful to calculate not only the modified but also the corrected (to standard) non-bonded potential energy while performing the dynamics or afterwards. This can be done straightforwardly using eqn (7) rewritten for U cor as This calculation induces little extra overhead during the simulation, considering that all the terms can be precalculated or are already calculated at each timestep. The only exception is the second term in eqn (11), but this is inexpensive. One may also consider the application of a correction similar to eqn (17) to the virial during the simulation. This option is explored in the ESI † and the corresponding results are only briefly discussed in the main article. Note that neither the standard nor the corrected energies encompass a long-range component for the LJ interactions, e.g. in the sense of a tail-correction, 1 in agreement with the calibration conditions for the GROMOS force fields. 21,135 Such a tail-correction could (and probably should) be added in the future, but this would require a consistent force-field reparametrization.

Simulation Details
The proposed modification schemes were implemented into the GROMOS program. 132 They were tested using a set of 57 small molecules, previously employed in the calibration of the GROMOS-compatible 2016H66 force field. 135 The same forcefield parameters were used in the present simulations. The list of compounds, along with key experimental properties, were taken from ref. 135 and are provided in Table S1 of the ESI. † A subset of 11 molecules (marked in the table) was selected as the calibration set, which was used to investigate the effect of parameter variations in the modification schemes. The full set of 57 molecules, referred to as the validation set, was then used to further test the schemes with appropriate parameter choices.
The calculations for each compound included 1 ns simulations in the liquid phase and 10 ns simulations in the gas phase. In both cases, the timestep was set to 2 fs, all bond lengths were constrained using SHAKE 140 with a relative geometric tolerance of 10 À4 , and the coordinates and energies were written to file every 1 ps.
The liquid-phase simulations were carried out using MD for a cubic box of 512 molecules under periodic boundary conditions in the isothermal-isobaric ensemble. The average temperature and pressure were maintained close to 298.15 K (unless indicated otherwise in Table S1 of the ESI †) and 1 bar using a Berendsen thermo-and barostat, 141 with coupling times of 0.1 and 0.5 ps, respectively, and a compressibility of 4.575 Â 10 À4 kJ À1 mol nm 3 for the barostat. The center-of-mass translational motion was removed every 2 ps. The simulations were started from pre-equilibrated configurations, and an additional 0.1 ns equilibration was performed for every new parameter setting. The gas-phase simulations were carried out using stochastic dynamics (SD) for a single isolated molecule with a friction coefficient of 91 ps À1 . The reference temperature was identical to that in the corresponding liquid simulation.
All simulations presented in the main article relied on a cutoff distance R c = 1.4 nm, a pairlist updated every timestep, the inclusion of long-range electrostatic interactions, and the omission of the long-range LJ component. These settings are compatible with the parametrization conditions of the 2016H66 force field. 21,135 The results of analogous simulations with R c = 1.2 nm instead are provided in the ESI. † Five schemes were considered for the long-range electrostatic interactions in the liquid-phase simulations: (i) a latticesum scheme (LS); (ii) the standard RF scheme with CG cutoff (RF/CG); (iii) the standard RF scheme with AT cutoff (RF/AT); (iv) SH-modified RF schemes with AT cutoff (RF/AT/SH); (v) SW-modified RF schemes with AT cutoff (RF/AT/SW). The lattice-sum scheme relied on the P 3 M method, 73,74 using a Gaussian charge-shaping function of width parameter 0.4 nm À1 and a grid of 64 Â 64 Â 64 points. For the RF schemes, the permittivity e RF was set to the experimental value for the compound (see Table S1 of the ESI †). For the LJ interactions in the first series of simulations, straight-cutoff truncation was used with a CG truncation when using RF/CG, or an AT truncation in all other cases.
After determining appropriate parameter choices for the SH and SW schemes applied to the electrostatic interactions, a second series of calculations considered the possibility of also modifying the LJ interactions with a SH or SW scheme. The additional schemes are labeled RF+LJ/AT/SH and RF+LJ/AT/SW. They rely on the same parameter choices for the RF modifications, selected as appropriate in the preceding phase.
The gas-phase simulations always relied on the same electrostatic and LJ scheme as used in the corresponding liquidphase simulation (including the possible application of the SH or SW modifications), with only two differences. The RF permittivity e RF was set to 1, and the lattice-sum case was mapped to exactly Coulombic interactions using straight-cutoff with AT truncation. Note that for the small molecules considered here, intramolecular interatomic distances larger than R c never occur in the simulations.
In addition, simulations for four model liquids were performed to assess the magnitude of possible cutoff artifacts (shown in Fig. 3, 4, 5, 8 and 9). The corresponding topologies are based on the standard GROMOS models for chloroform 137 and carbon tetrachloride, 136 as well as the simple-point-charge (SPC) water model. 93 For CCl 4 , the original model does not involve partial charges. Thus, charges of +0.4 e or +0.86e were introduced on the carbon atom (and neutralized by negative charges on the chlorine atoms) to test the effect of a molecular quadrupole. These simulations were carried out in the canonical ensemble at 298.15 K and densities set close to the corresponding experimental values. For water, they involved 2200 molecules instead of 512. For chloroform and water, the simulations lasted 10 ns to obtain converged DOCs and Kirkwood factors. For the other two solvents 1 ns simulations were performed. Besides this, the same parameters were used as for the other liquid simulations.

Analyses
The analysis of the trajectories was performed using the GROMOS++ analysis programs, 142 and involved the calculation of the liquid densities r liq , vaporization enthalpies DH vap , static relative dielectric permittivities e, and self-diffusion coefficients D. For the four model liquids, RDFs g(r), DOCs c(r), and Kirkwood factor G(r) were calculated.
The liquid density r liq as well as the total potential energies U liq and U gas were extracted from the trajectories using the GROMOS++ ene_ana program. The potential energies were corrected according to eqn (17), and the result used to calculate the vaporization enthalpy as where N = 512 is the number of molecules in the liquid simulations, R the universal gas constant, and T the absolute temperature.
The static relative dielectric permittivities e were determined using the applied electric-field method. 143,144 To this purpose, an external electric field E ext z ranging between 0.01 and 1 e nm À2 was applied along the z-direction of the computational box. For the lattice-sum calculations, the local field determining the forces on the individual charges is equal to E ext z . However, for the RF/CG calculations, this local field must be reduced by a factor 3e RF /(1 + 2e RF ) (see eqn (14) in ref. 143). This arises as a consequence of the non-Coulombic nature of the integrated dipole-dipole interaction tensor. For RF/AT, it turns out that the integration of this tensor leads to a result that is very close to the lattice-sum rather than the RF/CG one (see Appendix B in ref. 28). Consequently, for RF/AT, the local field determining the forces on the individual charges is set equal to E ext z , as in the lattice-sum case. The resulting average box dipole moment hM z i in the z-direction was evaluated from 0.5 ns simulations at each selected field strength. The permittivity e was then calculated using the GROMOS++ program eps_field according to where hVi is the average box volume. Following ref. 143, in order to avoid saturation effects, the range considered for E ext z was subsequently restricted by the condition (12k B pe 0 ) À1 mE ext z o 50, where m stands for the molecular dipole moment and k B for the Boltzmann constant. The reported values of e were obtained by averaging over the results of eqn (19) at the different field strengths satisfying this condition. For flexible molecules, m is configuration-dependent and was approximated by averaging over all molecules in the final configuration of the simulation. The permittivities of the liquids contained in the calibration set, which were obtained using the external electrc-field method, were also compared to the values calculated using corresponding fluctuation formulae in the ESI. † The two methods have been found to be in good agreement.
The self-diffusion coefficients D were calculated based on the liquid simulations using the GROMOS++ program diffus, which relies on the Einstein relation 145,146 In practice, D is evaluated from the slope of the atomic meansquare displacements with respect to time t in the longtime limit. The RDFs g(r) were calculated with the GROMOS++ program rdf according to where hN I i r,Dr is the average number of atoms of type I at a distance between r and r + Dr from the reference atom, and r I is the bulk number density of atoms I. Finally, the DOCs c(r) were calculated according to where f ij is the angle between the dipole moment vectors l of molecules i and j at a distance r. Distance-dependent Kirkwood factors G(r) were obtained using, ð r 0 r 02 gðr 0 Þcðr 0 Þdr (23) using the RDF g(r), the DOC c(r), and the bulk number density Z.

Modification schemes for RF-based electrostatic interactions
The deviations in the simulated r liq and DH vap values over the calibration set relative to the lattice-sum reference are illustrated in Fig. 6 for the RF/CG, RF/AT, RF/AT/SH and RF/AT/SW schemes with a cutoff distance R c = 1.4 nm. A number of different parameter combinations are compared for the SH and SW schemes applied to the RF-based electrostatic interactions. The corresponding results for R c = 1.2 nm are shown in Fig. S1 of the ESI. † For the SH scheme, the lowest exponents considered are m RF = 4 and n RF = 6 (4-6 SH scheme), because the choice m RF = 2 would correspond to changing the effective permittivity e RF . The choice of the highest exponent combination (12-14 SH scheme) is based on the onset of significant deviations from the lattice-sum curves in the RDFs of the model liquids ( Fig. 5e and f). For the SW scheme, the entire accessible range (0,1] in terms of the scaled switching distance d scl RF is considered (i.e. 0.05 to 1.0). As can be seen in Fig. 6, the 4-6 SH and the 0.75 SW schemes appear to reproduce the lattice-sum results most accurately with R c = 1.4 nm. The same holds when considering a 1.2 nm cutoff instead (see Fig. S1 of the ESI †). These specific parameter combinations are thus retained as the most appropriate for the implementation of smooth RF/AT/SH and RF/AT/SW schemes. For these specific schemes, most of the deviations in r liq and DH vap relative to lattice-sum are within one standard deviation of the purely statistical error affecting the calculated properties. Since the smoothing corrects for discontinuities at the cutoff but not for errors related to the instantaneous inhomogeneity/anisotropy of the medium outside the cutoff sphere, the above observation suggests that the errors related to the cutoff discontinuities are more pronounced than those related to the continuum approximation, at least for simple bulk systems and for the longest cutoff distance of 1.4 nm.
The validation of the selected 4-6 SH and 0.75 SW schemes in the context of RF/AT electrostatics on the full set of 57 organic compounds is shown in Fig. 7 for R c = 1.4 nm (and in Fig. S2 of the ESI † for R c = 1.2 nm). In addition to r liq and DH vap , the static relative dielectric permittivities e and the selfdiffusion coefficients D are also compared to the lattice-sum results. For both cutoff values, a direct comparison to the experimental data is also shown in Fig. S3 and S4 of the ESI. † However, the comparison against experiment is viewed here as less relevant than the comparison against lattice-sum, because it mixes errors from the force-field parametrization with those resulting exclusively from the non-bonded scheme. Fig. 6 Distributions of the absolute errors in the liquid density r liq (top) and vaporization enthalpy DH vap (bottom) over the calibration set for various RF-based electrostatic interaction schemes: RF/CG, RF/AT, RF/AT/SH and RF/AT/SW. All schemes rely on straight-cutoff truncation of the LJ interactions (CG when using RF/CG, AT otherwise). The errors are determined relative to the lattice-sum (LS) scheme along with straight-cutoff AT truncation of the LJ interactions as a reference. Different parameters are considered for the SH scheme (exponents m RF and n RF ) and the SW scheme (scaled switching distance d scl RF ). The cutoff distance is R c = 1.4 nm. The blue area sketches the error distribution, the central blue bar indicates the mean absolute error (MAE), the red circle shows the root-mean-square error (RMSE), and the gray area represents one standard deviation of the statistical error affecting the MAE. The latter value is obtained by evaluating the statistical error in each simulation by block averaging, and propagation onto an error on the MAE considering the set of molecules.
As was observed for the calibration set (Fig. 6), the 4-6 SH and 0.75 SW schemes substantially reduce the deviations in r liq and DH vap relative to lattice-sum for the validation set (Fig. 7a-c). Note, however, that the agreement with lattice-sum was already good for the original RF/CG and RF/AT schemes.
The static relative dielectric permittivities e (Fig. 7d) calculated using the RF/AT/SH and RF/AT/SW schemes remain close to those obtained with the RF/AT scheme. These values also agree very well with the reference lattice-sum values. The correlation of the RF/CG permittivities with lattice-sum is significantly worse compared to the schemes using an AT cutoff. For a few compounds, the present permittivity results are compared in the ESI † to values calculated using corresponding fluctuation formulae, and found to also agree very well.
The self-diffusion coefficients D (Fig. 7e) obtained using RF/AT tend to be somewhat lower compared to the results obtained using lattice-sum and the three other RF schemes. This reduced diffusion likely results from the preferential positioning of molecule pairs close to the cutoff distance R c , which appears in the form of an artificial peak in the RDF (Fig. 5e and f). The modification schemes eliminate these RDF artifacts and, as a result, the diffusive behavior in RF/AT/SH and RF/AT/SW matches very well the lattice-sum and RF/CG results.
As mentioned above, the choice m RF = 2 and n RF = 4 was not considered for the SH scheme. However, if one accepts to use a SH scheme corresponding to an altered effective permittivity e RF , the simplest alternative is to use the unmodified RF interaction with e RF -N. In this case, the corrected energy U cor is not evaluated using eqn (17), but by re-calculating the standard potential energy according to the actual permittivity e RF . For completeness, this option was also explored and the corresponding results are shown in Fig. S5 of the ESI † for R c = 1.4 nm (and Fig. S6 of the ESI † for R c = 1.2 nm). Even this simple modification leads to a significant improvement of the agreement with the lattice-sum results, although the 4-6 SH scheme still performs slightly better. Another option to consider is to also apply a correction similar to that of eqn (17) to the virial during the simulation. The results without and with such a correction in the context of the RF/AT/SH scheme are compared in Fig. S7 of the ESI † for R c = 1.4 nm (and Fig. S8 of the ESI † for R c = 1.2 nm). The virial correction slightly worsens the results and was thus not adopted. In particular, the values for benzene (like CCl 4 , a quadrupolar molecule with no dipole moment) differ significantly, due to changes in the RDF at short-range, as shown in Fig. S9 of the ESI. † Finally, in order to verify that the final RF/AT/SH and RF/AT/ SW schemes are indeed exempt of the cutoff artifacts observed in the original RF/CG and RF/AT schemes, the RDFs, DOCs and Kirkwood factors G(r) analogous to those of Fig. 3 and 4 are shown in Fig. 8 and 9, along with the curves from lattice-sum as a reference. Except for the slight oscillation in the RDF for water observed previously for RF/AT, which is still present, the RF/AT artifacts in the RDFs as well as the RF/CG artifact in the DOC of water have both been entirely removed. It can be seen in Fig. 9a that in the low-permittivity regime, the distance-dependent Kirkwood factor still deviates from the lattice-sum reference, resembling the RF/AT case (Fig. 4a). A similar alignment of RF/AT with RF/AT/SH and RF/AT/SW is also observed for the dielectric permittivity ( Fig. 7 and Fig. S2 of the ESI †), which is a dipolar property as well. However, it can be argued that due to the poor correlation of static dielectric permittivities with experimental values irrespective of the cutoff treatment ( Fig. S3 and S4 of the ESI †), this finding is of little practical relevance.

Modification schemes for the LJ interactions
The deviations in the simulated r liq and DH vap values over the calibration set relative to the lattice-sum reference are illustrated in Fig. 10  As can be seen in Fig. 10, the deviations of both r liq and DH vap appear to follow systematic trends with R c = 1.4 nm (and R c = 1.2 nm, see Fig. S10 of the ESI †). Relative to the lattice-sum reference, the deviations in both quantities become more limited when the modification scheme acts mainly at large distances, i.e. just below the cutoff (high exponents for SH, short switching distance for SW). The reason for these pronounced trends resides in the fact that the LJ interactions close to the cutoff distance are dominated by the dispersive (attractive) component. As a result, the modification scheme slightly damps the cohesive forces in the system, more so when the modification extends into the shorter distance range. The effect is more significant for r liq , because the modification influences the virial, provoking a slight expansion of the system and a decrease in the density. It is less significant for DH vap , because the potential energy is corrected for the effect of the modification, i.e. the change only reflects the effect of the slightly altered density, not of the altered energetics.
For SW, the choice of a 0.05 SW scheme, which only becomes active at a distance of 0.025 nm from R c , leads to reasonably small deviations of r liq and DH vap , and was retained for the following. The SH scheme is less suited to achieve similar distance scales for the modification, because the  magnitude of the required exponents would give rise to numerical instabilities. For example, the 0.05 SW scheme could be approximated by a SH scheme where the lower exponent is about 200, which would lead to coefficients a LJ6,m and a LJ6,n on the order of 10 À30 . For the subsequent validation, the 36-38 SH scheme was selected as a good compromise.
The validation of the selected 36-38 SH and 0.05 SW schemes for LJ interactions (along with 4-6 SW or 0.75 SH schemes for RF-based electrostatic interactions) on the full set of 57 organic compounds is shown in Fig. 11 for R c = 1.4 nm (and for R c = 1.2 nm in Fig. S11 of the ESI †).
As was observed for the calibration set (Fig. 10), a slight systematic underestimation of r liq and DH vap is observed. As expected, this trend is more pronounced for the SH compared to the SW scheme (Fig. 11a-c). The self-diffusion coefficients D (Fig. 11d) seem to be sensitive to the perturbation introduced in the potential-energy function, but the similar least-square fits to the data indicate that this might be a statistical rather than a systematic effect.
Here again, it is of interest to consider the application of a correction similar to that of eqn (17) to the virial during the simulation. Although this option was not retained in the case of the RF-based electrostatic interactions, in the context of the LJ interactions, it may provide a possibility to counteract the systematic decrease in density resulting from the modification of the potential-energy function. The results without and with such a virial correction, applied only for the LJ interactions in the context of the RF+LJ/AT/SH and RF+LJ/AT/SW schemes are compared in Fig. S12 of the ESI † for R c = 1.4 nm (and for R c = 1.2 nm in Fig. S13 of the ESI †). This virial correction noticeably improves the results. For this reason, one might consider adopting it in future work.

Conclusion
The RF approach for calculating long-range electrostatic interactions in MD simulations represents a common alternative to the use of lattice-sum methods. Compared to the latter methods, it limits the effect of artificial periodicity in nonperiodic systems (liquids and solutions), it involves a more favorable scaling of the computational cost with system size (linear), it enables a pairwise (or groupwise) partitioning of the interaction energy without additional overhead, and it alleviates intricate methodological issues for processes involving net-charge changes. The RF scheme is commonly applied with a CG truncation, to avoid the occurrence of severe artifacts in the context of low-permittivity solvents. However, AT truncation may be preferable when considering high-permittivity solvents, where it reduces the cutoff noise and leads to an improved representation of dipole-dipole correlations. Furthermore, CG truncation relies on the definition of neutral charge groups within molecules, generally requiring charge adjustments that are rather arbitrary as well as difficult to automate in the context of force-field design.
To address these issues, simple modification schemes are proposed here, which are compatible with an AT truncation for both RF-based electrostatic as well as LJ interactions. The objective of these modifications is to ensure that the pairwise force and its derivative vanish at the cutoff distance. However, instead of leveling the potential energy to zero, a constant contribution is added for all atom pairs at a distance larger than the cutoff. The calculation of the latter component is inexpensive and the modified potential energy permits to conduct rigorously conservative simulations. The energies can easily be corrected back to the unmodified form, either on the fly or in a post-processing step. Two different types of modification were considered.
The SH schemes rely a two-terms polynomial added to the pairwise interaction energy. The effective range of the shifting is determined by the selected even integer exponents m o and n o (with n o = m o + 2, m o Z 4). Due to the restriction to small even integer exponents, the calculation of the modified forces induces only limited overhead in the pairwise interaction calculation, provided that the exponents are pre-selected and hard-coded in the program.
The SW schemes rely on a cubic spline function added to the pairwise force. The perturbation range is determined by a scaled switching distance d scl o , with d scl o in the range (0,1]. The computational overhead of this scheme may be slightly higher, considering that its straightforward implementation adds two conditional (if) statements in the pairwise interaction calculation (or only one if the same switching range is used for the electrostatic and LJ interactions). Alternatively, the use of three separate range-based pairlists would permit to avoid these conditional statements, at the cost of more implementation efforts. In both cases, the scaled switching distance could be given as an input parameter to the simulation, and does not need to be hard-coded.
In practice, using the GROMOS code and RF/AT with 1.4 nm cutoff as a reference, the introduction of a shifting function in RF/AT/SH is found to increase the computational cost of the non-bonded interaction calculation by less than 1%, while the introduction of a switching function in RF/AT/SW increases it by about 8%. In contrast, the reciprocal-space evaluation involved in a lattice-sum P3M calculation (FFT operations) with 64 grid points in each dimension is much more expensive, with an increase in the calculation cost by a factor 2.8.
Simulations involving a set of 11 organic liquids were used to calibrate appropriate distance parameters for both schemes, separately for the RF-based electrostatic and the LJ interactions. Further simulations considering 57 organic liquids were then used to validate the retained SH and SW schemes in terms of liquid density, vaporization enthalpy, relative dielectric permittivity, and diffusion constant.
For the electrostatic interactions, an appropriate modification improves or, at least, does not mitigate the agreement with lattice-sum for all the monitored properties. For the permittivity calculation, it should be noted that the same equations are to be applied for the (modified) RF/AT schemes as for lattice-sum, whereas the equations differ in the RF/CG case. Both the SH and SW schemes can be tuned to give very similar modifications, and are in principle equally adequate for RF-based electrostatic interactions.
For the LJ interactions, even the least disruptive modifications increase the deviations relative to the lattice-sum reference. This is due to the slight damping of the dispersive (cohesive) interactions, and a resulting decrease in density. The effect is more pronounced for the density itself compared to the vaporization enthalpy. Due to the use of corrected energies, the change in vaporization enthalpy only reflects the effect of the altered density, not of the altered energetics. Correcting the virial similarly during the simulation noticeably improves the results for both quantities. One might thus consider adopting this strategy in future work. In the context of LJ interactions, the SW scheme leads to smaller deviations relative to lattice-sum compared to the SH scheme, because it enables shorter modification ranges without giving rise to numerical problems.
The influence of the modification schemes introduced in this work on the simulated properties are limited, suggesting that these schemes can be used together with the GROMOS force fields without requiring a reparametrization and without incurring a pronounced accuracy loss. Further tests on more complex (multi-component, biomolecular) systems will be performed, also considering the occurrence of charged species (ions, ion pairs, charged residues). It is important to stress that the modifications introduced here for the LJ interactions do not account for long-range LJ interactions. In the future, a proper LJ tail correction could (and probably should) be added, which will require a consistent force-field reparametrization.

Conflicts of interest
There are no conflicts to declare.

C. Coefficients of the switching function
Given a choice for the range parameter d o , analytical expressions are provided here for the five coefficients b o,0 to b o,4 in eqn (15). These coefficients also depend on the cutoff radius R c and the RF permittivity e RF . The expressions are determined by the conditions of eqn (12), along with corresponding requirements of continuity up to the second derivative at R c -d o , and reported here for o A {RF,LJ6,LJ12}. The coefficients are given by Note that the numerators are always negative, so that the switching distances d mid   (30) is shown as a function of the exponents m RF and n RF in eqn (14). The graphs rely on a cutoff distance R c = 1.4 nm and consider four different permittivities e RF . The lower triangle as well as the diagonal are irrelevant given the selected condition 0 o m o o n o specified in eqn (14).