The Effect of Using a Twin-Range Cutoff Scheme for Nonbonded Interactions: Implications for Force-Field Parametrization?

Recently, concerns have been voiced regarding the validity of the GROMOS force fields, being parametrized using a twin-range cutoff scheme, in which longer ranged nonbonded forces and energies are updated less frequently than shorter ranged ones. Here we demonstrate that the influence of such a scheme on the thermodynamic, structural, and dynamic properties used in the parametrization of the GROMOS force fields is minor. We find root-mean-square differences of maximally 0.5 kJ/mol for the solvation free energy and heat of vaporization and of maximally 0.4% for the density. Slightly larger differences are observed when switching from a group-based to an atom-based cutoff scheme. In cases where the twin-range cutoff scheme does result in minor differences compared to a single-range cutoff these are well within the deviation from the experimentally measured values.

R ecently, concerns have been voiced regarding the validity of the GROMOS force fields for biomolecular simulation. 1 In particular, when using the GROMOS force field in the popular molecular dynamics program GROMACS (version 2019.3 and newer), the following warning is issued: "The GROMOS force fields have been parametrized with a physically incorrect multiple-time-stepping scheme for a twin-range cut-off. When used with a single-range cut-off (or a correct Trotter multiple-time-stepping scheme), physical properties, such as the density, might differ from the intended values. Check if molecules in your system are affected by such issues before proceeding." Here, we address the following questions: Is this warning warranted and does the use of this time-saving technique affect the parametrization of the GROMOS force field in a relevant way?
The use of force fields to describe the interaction energy between atoms and molecules in molecular dynamics simulation is a commonly used approach, and the validation of such simulations is a crucial matter. 2 Apart from the accuracy of the force field, a wide range of simulation parameter settings will influence the outcome of a simulation. Ideally, the outcome of a simulation using a given force field should be independent of the parameter settings that were used to parametrize the force field. Unfortunately, this situation is extremely difficult, if not impossible, to reach. It is certainly not currently possible as atomistic simulations do not reach macroscopic space and time scales. This means that, by necessity, a range of approximations, including treating electrostatic interactions as either lattice-sum or cutoff with reaction field, are required. For this reason, it is generally prudent to use simulation parameter settings that are similar to those used to derive the force-field parameters.
The GROMOS force field parameters have been derived using a twin-range (TR) cutoff scheme for the nonbonded interactions as a simulation-time-saving technique. Typically, a pair-list is generated after a fixed time interval (e.g., 10 fs). The forces and energies up to a short-range cutoff (typically 0.8 nm) are computed at every time step according to this pair-list. At pair-list updates, also after the chosen time interval, forces and energies up to a long-range cutoff (typically 1.4 nm) are computed and kept constant between updates. 3 The discontinuity this introduces to the forces will lead to additional noise in the properties of the system, which will be small if the longer ranged nonbonded forces can be expected to change little between the updates. This approximation was introduced at a time when computational power was limited with tests at the time, suggesting the increase in computational efficiency outweighed any possible loss in accuracy. We note that the choices of a specific time step, van der Waals cutoff, shifting function, the use of constraints, or the choice of a particular level of numerical precision or solvation model in a simulation are all based on a similar balance between computational efficiency and accuracy. 4 However, since the current GROMOS force field was parametrized using this approximation, it is reasonable to ask if this approximation affected the parametrization to such an extent that if pair-list and nonbonded forces were computed at every time step (single range, SR), the results of the simulation would no longer be valid.
Before presenting our findings we note that there have been a number of recent reports suggesting that the results obtained in simulations performed using the twin-range scheme and simulations performed using a single-range cutoff were essentially identical. No significant differences were observed in the area per lipid of a membrane, when the only difference in the simulation parameter settings was the use of a SR vs a TR cutoff scheme. 5 Similarly, the radius of gyration of a polyamidoamine dendrimer and the structural properties of a protein or a membrane in constant pH simulations were indistinguishable when using SR vs TR cutoffs. 6 However, differences were found when using an atomistic (AT) cutoff rather than a charge-group (CG)-based cutoff for nonbonded forces. 6 Also, the potentials of mean force between small solutes and carbon nanotube model systems were not significantly affected by the choice of SR vs TR cutoffs. 7 None of these studies directly addressed the question of Averages over three independent sets of simulations are reported with errors reported as standard deviations. b Reference 10. c Reference 9.   Table 1, with significant differences (p < 0.01) indicated by asterisks. whether the types of properties used in the parametrization of the GROMOS force field were affected by the use of a twinrange cutoff scheme. In a related and very extensive study, Gonsalves et al. examined the effect of various simulation parameter settings on the thermodynamic and transport properties of pure liquids. 8 Significant differences were observed, depending on the exact choice of parameters governing the computation of nonbonded interaction energies. However, a direct comparison of SR and TR using otherwise identical parameter settings was not considered.
The nonbonded parameters of the GROMOS force fields were parametrized against thermodynamic properties of small molecules. Specifically, the density and heat of vaporization of pure liquids and the solvation free energy of amino acid sidechain analogues in water and in cyclohexane. 9−11 Here, the hydration free energies for the neutral amino acid side-chain analogues used in the parametrization of the GROMOS force field are recalculated using SR and TR cutoff schemes. In addition, the effect of the cutoff scheme on key properties of a set of five liquids of different polarity and molecular complexity is also examined. All calculations were performed using both an AT-and a CG-based cutoff, to further clarify whether any of the key assumptions made during the parametrization of the force field affects the outcome of the simulations. The GROMOS force fields were initially parametrized using a charge-group cutoff scheme. This was used as it reduces the artifacts in the forces due to atoms at the cutoff distance. The main question we aim to address is whether the GROMOS force field parameters would have been affected, if a SR cutoff was used instead of the TR cutoff.
To address this question, four sets of simulations were performed, two with a SR cutoff scheme and two with a TR cutoff scheme. In both cases either a CG-based cutoff and pairlist or an AT cutoff and pair-list were used. All simulations were performed using GROMOS11. 12 Note that other simulation packages may not support all four combinations of cutoff schemes.
Analogous to our previous work, the amino acid side-chain analogues were generated by breaking the C α −C β bond in the naturally occurring neutral amino acids and adjusting the parameters of the united atom C β carbon. 10 The molecules were placed in a rectangular periodic box and solvated in SPC water. 13 The systems were then simulated at a constant reference temperature of 298.15 K and a reference pressure of 1 atm, using the weak-coupling scheme. 14 The coupling times were τ T = 0.1 ps and τ P = 0.5 ps. The isothermal compressibility was set to 4.575 × 10 −4 (kJ·mol −1 ·nm −3 ) −1 . Separate temperature baths were used for the solute and solvent. Bond lengths were constrained using the SHAKE algorithm, with a relative geometric accuracy of 10 −4 . 15 A time step of 2 fs was used in the leapfrog algorithm 16 to integrate the equations of motion. For the TR cutoff scheme, a shortrange cutoff of 0.8 nm and a long-range cutoff of 1.4 nm were used with updates of longer ranged forces and energies and the pair-list every 10 fs (five steps). The SR cutoff scheme involved a single cutoff of 1.4 nm, with an update of pair-list and forces and energies at every step. A reaction-field contribution to the energies and forces was added to the electrostatic interactions to account for a homogeneous medium with a relative dielectric constant of 61 17 beyond the 1.4 nm cutoff. 18 Solvation free energies were calculated by alchemically removing the nonbonded interactions of the amino acid sidechain analogues using a λ-coupling parameter approach and thermodynamic integration. 19 Simulations were performed at 21 equally spaced λ-values. At each λ-value, the systems were

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Letter equilibrated for 100 ps and data were collected for 1 ns. Every set of simulations was performed three times, using different random number seeds for the generation of initial velocities. Error estimates are reported as standard deviations over the three independent estimates. Five liquids of different polarity were also examined: hexane, butylamine (BAN), ethanol, dimethyl sulfoxide (DMSO), and water. Each system consisted of 1000 molecules. The systems were simulated at constant pressure, using the same simulation parameter settings as above. The only differences were that the isothermal compressibility for the pressure scaling algorithm and the relative dielectric constant for the reaction-field contributions to long-range interactions were set to values derived from experiment. 6 The translational and internal/ rotational degrees of freedom were coupled to separate temperature baths. Three independent simulations of 10 ns each were performed, using different initial velocities. The last 5 ns of every simulation was used to analyze the data. Error estimates are reported as the standard error of the mean over the three simulations.
Simulations in the gas phase were performed by simulating a single amino acid analogue or by placing the molecules at a distance of at least 50 nm from each other (no intermolecular interactions) for multimolecular systems. In these simulations a Langevin thermostat 20 involving stochastic forces was added to help the molecule to escape conformational local minima (friction coefficient γ = 91 ps −1 ). Diffusion constants (D) were estimated from the Einstein relation; rotational relaxation constants (τ 2 ) were estimated from the autocorrelation function of the second order Legendre polynomial of a molecular axis. To estimate the relative dielectric permittivity (ε r ), five additional simulations were performed in the  Table 1 shows the average solvation free energies using the different cutoff schemes. For comparison the literature values from the original parametrizations 9,10 are given, as well as the experimental data used in the parametrization. 22 Deviations between the literature values and the current values can be traced to longer simulation times and the use of more λ values in the thermodynamic integration in the current work. Absolute differences between the nonbonded cutoff schemes for TR vs SR and CG vs AT comparisons are given in Figure 1, along with an indication of significant differences (p < 0.01).
The differences between TR and SR are very small. The largest deviation is observed for the histidine side chain at 1.1 kJ·mol −1 when using an group-based cutoff. The top half of Table 2 shows the root-mean-square differences for comparisons between each of the columns in Table 1. A comparison between TR and SR amounts to 0.53 kJ·mol −1 (CG) or 0.52 kJ·mol −1 (AT). Using a two-tailed t test, the observed differences are significant only for Arg when using a groupbased cutoff and for Arg, His, and Trp when using an atombased cutoff. The differences between CG and AT are slightly larger, with a maximum difference of 2.4 kJ·mol −1 for Arg when using a twin-range cutoff scheme. Root-mean-square differences amount to 0.85 kJ·mol −1 (SR) and 0.92 kJ·mol −1 (TR). The t test reveals significant differences for seven (SR) and nine (TR) compounds. The lower half of Table 2 shows the average differences between the columns of Table 1. Comparison of these values to the RMSD values in the upper half reveals that the small differences are largely systematic. Overall, the largest difference in the simulations is observed for CG/TR vs AT/SR, with an RMSD of 1.32 kJ· mol −1 , where the two effects (CG to AT and TR to SR) seem to add up. For the comparison CG/SR vs AT/TR, the RMSD amounts to 0.52 kJ·mol −1 and systematic effects partially cancel. Table 3 shows selected properties for the pure liquids. Comparing the SR and TR schemes, the largest deviation in terms of the heat of vaporization amounts to 0.6 kJ·mol −1 for butylamine using AT. The root-mean-square deviation over these five liquids amounts to 0.24 kJ·mol −1 (CG) and 0.42 kJ· mol −1 (AT). The densities of all liquids are rather similar. The largest deviation in a SR vs TR comparison is seen for butylamine and amounts to 0.9% (AT). Root-mean-square deviations amount to 2.0 kg·m −3 (CG) and 3.0 kg·m −3 (AT) or 0.3% (CG) and 0.4% (AT). Differences regarding the comparison CG vs AT are slightly larger with root-meansquare differences of 3.8 kg·m −3 (TR) and 3.5 kg·m −3 (SR). Figure 2 shows the radial distribution functions and dipole correlation functions for the simulations of the polar liquids. As observed before, differences are found when comparing CG and AT simulations. 6,23 When comparing simulations performed with SR to those with TR, the curves for both of these structural properties are indistinguishable. A detailed comparison for group-based and atom-based cutoff schemes shows that artifacts at the cutoff appear for atom-based cutoff schemes. 6,24 Dynamic properties were computed for the pure liquids too. For the diffusion coefficient, τ 2 rotational relaxation time and the relative dielectric permittivity, the deviations between SR and TR are negligible. For the dielectric permittivity, the difference between AT and CG that was previously reported 6 is reproduced.
In conclusion, we have repeated the key simulations used in the parametrization of the GROMOS force fields for 18 amino acid side-chain analogues and five pure liquids of different polarities. We have shown that the effect of using a twin-range cutoff scheme as compared to a single-range cutoff on the thermodynamic, structural and dynamic properties is minor. All differences that were observed are well within the deviation from the experimentally measured values and are likely to be irrelevant in the vast majority of applications of biomolecular simulations. 24 On the basis of these results we can conclude that the parametrization of the GROMOS force fields has not been influenced by the use of the twin-range cutoff scheme to a large extent. While ideally all force fields should be parametrized without the use of any time-saving approximations, this has not been true in the past. At least in the case of the GROMOS family of force fields we can safely say that the use of a single-as opposed to a twin-range cutoff has had no effect on the parametrization in terms of the agreement with experiment. Thus, in codes where a GROMOS-type twinrange is not implemented users should use a group-based, single-range cutoff to obtain equivalent results.