Optimized Magnesium Force Field Parameters for Biomolecular Simulations with Accurate Solvation, Ion-Binding, and Water-Exchange Properties in SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D

Magnesium is essential in many vital processes. To correctly describe Mg2+ in physiological processes by molecular dynamics simulations, accurate force fields are fundamental. Despite the importance, force fields based on the commonly used 12-6 Lennard-Jones potential showed significant shortcomings. Recently progress was made by an optimization procedure that implicitly accounts for polarizability. The resulting microMg and nanoMg force fields (J. Chem. Theory Comput.2021, 17, 2530–2540) accurately reproduce a broad range of experimental solution properties and the binding affinity to nucleic acids in TIP3P water. Since countless simulation studies rely on available water models and ion force fields, we here extend the optimization and provide Mg2+ parameters in combination with the SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water models. For each water model, the Mg2+ force fields reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl2 solutions, and the binding affinity and distance to the phosphate oxygens on nucleic acids. We present two parameter sets: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. Depending on the water model, nanoMg yields accelerated water exchange in the range of 106 to 108 exchanges per second. The nanoMg parameters can be used to enhance the sampling of binding events, to obtain converged distributions of Mg2+, or to predict ion binding sites in biomolecular simulations. The parameter files are freely available at https://github.com/bio-phys/optimizedMgFFs.


INTRODUCTION
Magnesium is the second most abundant intracellular cation. 1 It is involved in more than 600 enzymatic reactions 2 and plays a crucial role in vital processes such as ATP hydrolysis, 3 cellular signaling, 2 or the catalytic activity of ribozymes. 4 The specific requirement for Mg 2+ is prevalent in nucleic acid systems where Mg 2+ is crucial for the stability, folding, and biological function. 5−9 In addition to its physiological relevance, Mg 2+ is important in DNA nanotechnology 10−12 and is a promising candidate for divalent batteries. 13,14 Due to the distinct role of Mg 2+ , the modeling of Mg 2+ by all-atom molecular dynamics (MD) simulations has received significant attention. 15−26 Despite tremendous efforts, Mg 2+ force fields based on the commonly used 12-6 Lennard-Jones interaction potential failed to provide a quantitative description (see Table S1) for three reasons. (i) No parameter combination existed that simultaneously reproduced the solvation free energy and the size of the first hydration shell. 20,21,23,27 (ii) The parameters yielded too low water exchange rates leading to unrealistically slow exchange kinetics such that transitions from outer to inner sphere binding and back could never be observed on the typical time scale of MD simulations. 28 (iii) The binding affinity of Mg 2+ to ion binding sites on biomolecules was overrated significantly. 19,25,29,30 Regarding these shortcomings, the immediate question arises why classical, nonpolarizable Mg 2+ force fields fail to provide an accurate description. Clearly, Mg 2+ ions strongly polarize their environment, and the lack of charge transfer and polarization effects in classical simulations likely leads to the observed deviations between experiments and simulations. Possibilities to provide improvement include the use of computationally more demanding polarizable force fields, 31 to include additional parameters in the interaction potential, 22,30,32 or to scale the charge of the Mg 2+ ion. 18,33 Another possibility is to include the effect of polarizability implicitly without introducing additional terms in the interaction potential. As shown in our previous work, 26 polarizability can be taken into account by two measures. First, optimizing the parameters of the 12-6 Lennard-Jones parameters in an enlarged parameter space captures the more attractive and long-ranged Mg 2+ −water interactions, thereby facilitating the simultaneous matching of the experimental solvation free energy, size of the first hydration shell, and water exchange rate. Second, by introducing scaling factors 34,35 in the standard combination rules for the Mg 2+ −Cl − and Mg 2+ −RNA interactions, deviations due to polarization are taken into account, and the experimental activity derivative and the binding affinity to the phosphate oxygen on nucleic acids is reproduced. The resulting force field parameters provide an efficient and highly accurate model for Mg 2+ in biomolecular simulations in combination with the TIP3P water. The TIP3P water model is frequently used in biomolecular simulations since protein and nucleic acid force fields were frequently optimized in the presence of TIP3P. However, to date, a large variety of different water models exist, some of which reproduce the physical properties of water better compared to TIP3P. 36 In principle, it is possible to combine more advanced water models and ion parameters with the force fields for biomolecules. In some cases, such a combination can leverage the strengths of the respective parameters and improve the agreement between simulated and experimental structures. 25,37−42 However, the combination of force field parameters does not guarantee that the physical properties, which were targeted in the first place, are reproduced. For metal cations, previous studies indicate limited transferability and different water models can have significant effects on the solvation free energy, the exchange kinetics, and even the reaction mechanism. 43−45 It is therefore crucial to assess the transferability of Mg 2+ parameters to different water models and optimize the force field parameters if necessary in order to ensure physically meaningful results.
The aim of our current work is to provide optimized Mg 2+ force field parameters in combination with five different water models. We focus on some of the most popular rigid nonpolarizable water models, namely, SPC/E, 46 TIP3P-fb, 47 TIP4P/2005, 48 TIP4P-Ew, 49 and TIP4P-D. 50 Our results illustrate that the transferability of the Mg 2+ parameters 26 developed in combination with TIP3P is limited. In particular, for the 4-site water models, significant deviations from the experimental properties are observed, in agreement with similar studies. 44 In order to provide accurate parameters for the different water models, we systematically derive Mg 2+ parameters that reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl 2 solutions, and the binding affinity and distance to nonbridging phosphate oxygens on nucleic acids, using our previously developed optimization procedure. 26 In order to provide consistent and robust results, Cl − was chosen as the reference ion, and its parameters were optimized in a preceding step. For each water model, we present two parameter sets: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. NanoMg yields accelerated water exchange in the range of 10 6 to 10 8 exchanges per second depending on the water model used. The nanoMg parameters are suited to accelerate the binding kinetics in biomolecular simulations and improve the sampling of ionic distributions. Our results reveal that the largest speed-up is obtained in combination with TIP3P or SPC/E.  σ ii , ε ii , σ io , and ε io are the ion−ion and ion−water LJ parameters. λ σ X and λ ε X are the scaling factors for the Lorentz−Berthelot combination rules (eq 2) for the interaction with Cl − or the RNA atoms, shown exemplary for the interaction between Mg 2+ and OP. Note that the scaling factors are only valid in combination with the Cl − parameters from Smith−Dang 54 for SPC/E water and those listed in Table 2 for the other water models, and the parmBSC0χ OL3 RNA parameters. 84−86 Values for parameters in TIP3P are taken from ref 26. pairwise interaction potential. We use the most common form of the Lennard-Jones (LJ) potential with a repulsive r −12 and an attractive r −6 term where q i and q j are the charges of atoms i and j, r ij is the distance between them, and ϵ 0 is the dielectric constant of the vacuum. The parameters σ ij and ε ij describe the LJ diameter and interaction strength, respectively. Following the procedure to optimize the force field parameters for Mg 2+ in TIP3P water from our previous work, 26 we systematically optimize the two parameters of the LJ potential and the combination rules for 5 different nonpolarizable rigid water models. From the different three site models, we chose the most commonly used SPC/E 46 and the more recent TIP3P-fb 47 49 and TIP4P-D. 50 TIP4P/2005 has gained popularity and is often quoted as the best nonpolarizable general-purpose model. 36 In addition to TIP3P, TIP4P-Ew is frequently used in the AMBER and CHARMM force fields. 30,51,52 TIP4P-D is one of the newer offsprings in the TIP4P family and was designed to improve water dispersion interactions. Partial charges, Lennard-Jones parameters, bond lengths, and bond angles of the different waters are shown in Section S1.2.
Since the standard combination rules do not take polarizability and charge transfer effects into account, we use scaling parameters in the Lorentz−Berthelot combination rules to describe the Mg 2+ −Cl − and the Mg 2+ −RNA interactions. 26 Following the work by Fyta and Netz 35 and previous similar studies, 29 where X denotes Cl − or the atoms of the RNA. Note that the solvation free energy, the structural properties of the first hydration shell, and the rate of water exchange remain unchanged upon changing the Mg 2+ −Cl − and the Mg 2+ − RNA combination rule. In this work four scaling parameters are introduced, λ σ RNA , λ ϵ RNA , λ σ Cl , and λ ϵ Cl , to optimize ion−ion interactions and Mg 2+ −DMP interactions. Note that, in the RNA system, the same scaling factors λ σ RNA and λ ϵ RNA are applied to all atoms of the RNA.
To compare the optimized parameters (Table 1) to the available force fields from the literature, we performed simulations using the 12-6 based parameters by Åquist, 15 Mamatkulov−Netz, 20 and Li−Merz 21 (HFE set) and the 12-6-4 based parameters by Li−-Merz 22 for SPC/E water. In addition, the parameters by Li−Merz 21,22,32 (12-6 and 12-6-4) for TIP4P-Ew and TIP3P-fb were used (data not shown). For SPC/E water, the Cl − parameters were taken from Smith− Dang. 54 For all other water models, the Cl − parameters were adjusted as described below.
Simulations with force fields of the 12-6 type were performed with GROMACS 55 (versions 2018, 2020). Simulations with force fields of the 12-6-4 type were done with AMBER 52 (version 2018) and PLUMED 56 (version 2.6.4) since GROMACS does not support 12-6-4 interaction potentials. Section S1.3 (Tables S3 and S4) lists all simulation setups. Similar to previous work, 19,26,30,57,58 we used dimethylphosphate (DMP) to optimize the ion−RNA interactions. The force field parameters for the DMP molecule are based on a parametrization with GAFF 59 and can be found in ref 26. The analysis was performed with built-in GROMACS 55 code and by using the MDAnalysis package 60,61 for python.
2.2. Optimization Procedure. The optimization procedure follows the same three step strategy described in ref 26, starting with an optimization of the ion−water interaction by performing a grid search in σ io −ε io space. Here, in a first step, all σ io −ε io parameter combinations matching the experimental solvation free energy ΔG solv , the Mg 2+ -oxygen distance in the first hydration shell R 1 , and the coordination number n 1 are selected.
Second, by calculating the rate of water exchange for the above-mentioned parameter combinations, we optimize the water exchange dynamics. For each water model two parameter sets were chosen: The microMg parameter sets yield water exchange on the microsecond time scale and reproduce the experimental rate within errors. The second parameter sets yield faster water exchange time scales providing exchange dynamics that range between 10 6 to 10 8 exchanges per second while simultaneously reproducing all thermodynamic and structural properties. For consistency with our previous work, 26 we refer to the second set as nanoMg, even though the nanosecond time scale could not be reached for all water models.
In the final step, we optimize the ion−ion and ion−RNA interactions by calculating activity coefficient derivatives and ion binding affinities. We performed a grid search in λ σ X and λ ε X parameter space (eq 2). The activity coefficient derivatives a cc are collected and the scaling factors λ σ Cl and λ ε Cl are selected by reproducing the experimental activity derivative over a broad range of MgCl 2 concentrations, using the Kirkwood−Buff theory. 62 To calculate the binding affinity of Mg 2+ toward one of the nonbridging phosphate oxygens of DMP, we used alchemical transformation calculations. Reproducing the experimental binding affinity ΔG b 0 and binding distance R b toward the phosphate oxygen, we select the scaling factors λ σ RNA , λ ε RNA .

Solvation Free Energy.
The solvation free energy is considered the most important thermodynamic property in the development of accurate force field parameters. 27 Yet, the proton solvation free energy can vary depending on the experimental sources. In most cases, the solvation free energy of the proton by Tissandier et al., 63 −1104 kJ/mol, or Marcus, 64 −1056 kJ/mol, are used. Since, depending on the exact choice, the estimates can vary by more than 50 kJ/mol, we take the more robust solvation free energy of neutral ion pairs into account. In the following, we use chloride as the reference ion as in previous works. 20,23,27 In order to obtain consistent results, appropriate corrections for finite size effects, compression, and interfacial crossing have to be applied. Finite size effects are taken into account via 65 where z is the valency, N A Avogardo's number, e the unit charge, ϵ 0 the vacuum permittivity, and R 1 the first peak of the ion−water radial distribution function. Hence, the effective ion radius, ζ ew = −2.837 279/L, is the Wigner potential, with L being the edge length of the simulation box in nanometers. 65 47,49,50 The correction term related to the compression of the gas is given by 27 and is independent from the choice of the water model, with p 0 = 1 atm being the pressure of ideal gas and p 1 = 24.6 atm corresponding to an ideal solution under pressure at a density of 1 mol/L. In the experiments, 63 the ions have to pass the air−water interface in order to enter into the aqueous phase. The corresponding free energy correction term is 27 Here, we chose the surface potential as ϕ surf = −0.527 V obtained for polarizable TIP4P-FQ water. 67,68 The reason is that this choice closely matches with the experimentally based estimation of −0.50 V. 69 Even more importantly, the interfacial crossing term almost exactly cancels the shift between the solvation free energies by Marcus 64 and by Tissandier 63 for different anions and cations. Consequently, we obtain consistent results if we compare the simulation results with the interfacial crossing term to the free energies by Tissandier (since those values include a full contribution from the bulk water surface potential 69 ) or if we compare the simulation results without interfacial crossing term to the bulk free energies by Marcus (which omit the surface potential 69 ).
Using the former approach, the solvation free energy is given by The solvation free energy of the neutral MgCl 2 ion pairs is given by 2.4. Parametrization of Cl − and Mg 2+ . Initially, we optimize the parameters of the reference chloride ion. For SPC/E water, we use the well established Smith−Dang parameters. 54 In particular, the calculated solvation free energy is ΔG solv Cl− = −306 kJ/mol 27 and closely matches the experimental results by Tissandier 63 (−304.2 kJ/mol). For the optimization of the parameters for Cl − in combination with TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water, we used two different approaches: In the first approach, we modified ε io starting from the Smith−Dang parameters 54 (σ io = 0.378 nm, ε io = 0.52 kJ/mol). We selected the values for ε io that reproduced the experimental results while keeping σ io fixed. In the second approach, we modified σ io starting from the Joung−Cheatham chloride parameters 70 for TIP4P-Ew water (σ io = 0.404 105 546 nm, ε io = 0.182 336 936 kJ/mol). We selected the values for σ io that reproduced the experimental results while keeping ε io fixed. As shown in Figure 1, the second approach allowed us to simultaneously reproduce ΔG solv Cl− and R 1 , the radius of the first hydration shell, for all water models. The parameters from the second approach were therefore used in the optimization of Mg 2+ . The resulting LJ parameters are listed in Table 2.
Subsequently, the solvation free energy for Mg 2+ was calculated on a grid with σ io = 0.195−0.22 nm and ε io = 0.02−28 kJ/mol. In addition, R 1 , n 1 , and the self-diffusion coefficient D 0 were calculated for all waters as described in ref 26. 2.5. Free Energy Profiles, Water Exchange Rate, Activity Derivative, and Binding Affinity. We used umbrella sampling 71,72 to calculate one-and two-dimensional free energy profiles, straightforward simulations to calculate the water exchange rate (see Section S1.5 for more details), Kirkwood−Buff theory 62 to calculate the activity derivative with respect to the natural logarithm of the number density of the Mg 2+ ions (see Section S1.6 for more details), and alchemical transformation to calculate the binding affinities. Additional details regarding the calculation of free energy profiles (Section S1.7) and binding affinities (Section S1.8) can be found in the Supporting Information of this work.

RESULTS AND DISCUSSION
In the following, we present the results from our parameter optimization of Mg 2+ in combination with the SPC/E, TIP3Pfb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water models. Our optimization procedure is done in three sequential steps 26 and allows us to reproduce a broad range of thermodynamic properties including the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl 2 solutions, and the binding affinity and distance to the phosphate oxygens of RNA (Tables 3 and 5). For each water model, we present two sets of optimal parameters: MicroMg yields water exchange on the microsecond time scale and matches the experimental exchange rate. NanoMg yields accelerated water exchange in the range of 10 6 to 10 8 exchanges per second depending on the water model used (Table 4).
3.1. Transferability of Ion Parameters between Different Water Models. We tested the transferability of the recently developed Mg 2+ parameters 26 microMg(TIP3P) and nanoMg(TIP3P), that were optimized in combination with the TIP3P water model. Figure 2 and Table S3 show that the transferability of the parameters to other water models is limited and deviations from the experimental solvation free energy are observed. While for the 3-site water models (TIP3P-fb and SPC/E) the deviations are small (3−4 kJ/mol), significant deviations up to 80 kJ/mol are observed for the 4site models (TIP4P/2005, TIP4P-Ew, and TIP4P-D). The influence of the different water models on the solvation free energy ΔG solv , the radius R 1 of the first hydration shell, and the coordination number n 1 is also reflected in the free energy isolines which display differences between 3-and 4-site water models ( Figure S3). Since the transferability of the Mg 2+ parameters is limited, in agreement with a similar previous study on metal cations, 44 we have systematically optimized the parameters in combination with 5 different water models.
3.2. Optimization of Solvation Free Energy, Radius of Hydration Shell, and Coordination Number. For all 5 water models, we selected parameter combinations for σ io and ε io that accurately reproduce ΔG solv , R 1 , and n 1 simultaneously ( Figure 3A and Table 3). Hereby, the key to a successful parametrization was to take polarizability into account implicitly by considering a larger range of possible LJ parameters. 26 As shown exemplary for SPC/E water in Figure  3B (for the other waters see Figure S4), the interaction potentials with our optimized parameter sets microMg and  49 and TIP4P-D 50 ), resulting in the converted parameter sets microMg(TIP3P|W) and nanoMg(TIP3P|W). Additional single-ion properties are shown in Table S7.  Figure S4.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article nanoMg are more long-ranged compared to the standard parameter sets from the literature 15,20,21 and similar to 12-6-4 potentials, which explicitly consider ion−dipole interactions via additional terms in the interaction potential. 22 Thus, our optimized parameters take polarization effects into account implicitly, and the interaction between Mg 2+ and water is more attractive and long-ranged compared to the standard force fields from the literature. This in turn allows us to reproduce the experimental results for ΔG solv , R 1 , and n 1 simultaneously and to achieve agreement comparable to 12-6-4 interaction potentials 22 ( Figure 3A) without introducing additional parameters in the functional form of the interaction potential (eq 1).

Optimization of the Water Exchange Rate.
In the next step of our optimization, we selected two sets of parameters from the data that reproduced ΔG solv , R 1 , and n 1 in the previous step based on the calculated water exchange rate. In order to provide an accurate estimate of the rate constant, the rate was calculated from 1 μs long straightforward simulations since transition state theory, albeit computationally less demanding, can overestimate the true rate by more than 2 orders of magnitude. 28 The rates for the two sets of parameters for all water models are shown in Figure 4A and Table 4. The first parameter sets, microMg, yield water exchange on the microsecond time scale in agreement with experimental results. 73,74 The second parameter sets, nanoMg, yield accelerated water exchange with 10 6 to 10 8 exchanges per second dependent on the water model used (Table 4). In all cases, the parameters for nanoMg were chosen such that they yield the highest rate while still reproducing all other experimental properties. Yet, the different properties of the water models, 36 in particular the dielectric constant of the models, lead to distinct differences in the parameter range that reproduces ΔG solv , R 1 , and n 1 ( Figure S3). Consequently, the maximum exchange rate that could be achieved differs for the different water models. NanoMg(SPC/E) yields the highest exchange rate with 10 8 exchanges per second, similar to the acceleration achieved for TIP3P. 26 NanoMg(TIP4P-Ew) yields the next highest rate with 10 7 s −1 followed by nanoMg(TIP3Pfb), nanoMg(TIP4P/2005), and nanoMg(TIP4P-D) with 10 6 exchanges per second (Table 4 and Figure 4A).
The improvement of the exchange rate compared to the standard sets from the literature is shown exemplary for SPC/E in Figure 4B (for the other waters, see Figure S5). The barrier height between the two stable states (i) and (ii), which corresponds to replacing one water molecule in the first hydration shell with a water molecule from the second shell (green and blue water molecules shown in Figure 4C), differs by more than 10 k B T for different force fields ( Figure 4B). The different barrier heights provide a qualitative explanation of the different water exchange rates obtained with different force fields from the literature. Force fields with slow exchange kinetics (Åquist, 15 Mamatkulov−Netz, 20 and Li−Merz 21 ) have high free energy barriers while force fields with a more attractive and long-ranged interaction potential (microMg and Li−Merz [12-6-4] 22 ) have lower energetic barriers and yield exchange rates close to the experimental result. Finally, nanoMg has the lowest barrier, resulting in the fastest exchange kinetics that is particularly useful for enhanced sampling of Mg 2+ binding, as will be discussed further below.
The microscopic mechanism of water exchange is, however, more complex than the one-dimensional free energy profiles suggest. The exchange involves the concerted motion of two water molecules 28,45 and is captured more realistically by twodimensional free energy profiles (Figures S6 and S7).
In addition, the self-diffusion coefficients D 0 were calculated for each water model. Unlike the results in TIP3P water, where the diffusion coefficient matched the experimental value without optimization, D 0 is underestimated for all water models investigated in the present study (Table 3). Likely, the better agreement in TIP3P is due to the higher self-diffusion constant of TIP3P. 36 3.4. Optimization of the Activity Derivative. Subsequently, we balance ion−water and ion−ion interactions by matching with experimental activity derivatives. 34,35,53,75,76 Similarly to previous work, 20,23,25,26,35 the standard combination rules (eq 2) (with λ σ,ε Cl− = 1.0) overestimate the anion− cation interaction and consequently underestimate the activity derivative over the full parameter range. Since the standard combination rules are valid only in idealized cases, polarization and charge transfer can lead to deviations. Introducing scaling factors in the combination rule allows us to take these effects into account and to provide closer agreement with Solvation free energy of neutral MgCl 2 ion pairs ΔG solv , Mg 2+ −oxygen distance in the first hydration shell R 1 , coordination number of the first hydration shell n 1 , self-diffusion coefficient D 0 , and a cc the activity derivative of a MgCl 2 solution at 0.25 M concentration. Values for parameters in TIP3P are taken from ref 26. experimental results without changing ΔG solv , R 1 , n 1 , D 0 , or k. In all cases, increasing the effective size, σ MgCl , and decreasing the cation−anion interaction energy, ε MgCl , reproduces the experimental activity coefficient derivative of MgCl 2 solutions over a broad concentration range ( Figure 5). In particular, we found one set of universal scaling factors (λ σ Cl = 1.59 and λ ε Cl = 0.1) that leads to reasonable agreement with the experimental values for all five water models (Table 3, Figure 5). Note, however, that, at high concentrations, individually adjusted combination rules for water models such as TIP4P-D or SPC/ E might lead to better agreement.
3.5. Optimization of the Binding Affinity to Nucleic Acids. So far, the optimization was based only on bulk ion properties. However, the transferability of such bulk-optimized ion parameters to describe the interactions of Mg 2+ with biomolecules turned out to be limited. 25,26,30 Therefore, in a last step, we optimized the interaction between Mg 2+ and nucleic acids by matching with the experimental binding affinity and distance to one of the nonbridging phosphate Kinetic rate coefficients k (eq S3) of microMg (blue) and nanoMg (green) of different water models (SPC/E, 46 TIP3P-fb, 47 TIP4P/ 2005, 48 TIP4P-Ew, 49 and TIP4P-D 50 ) and for TIP3P from the literature. 26 The gray vertical line indicates the experimental values 73,74 (Table 4). (B) One-dimensional free energy profiles as a function of the distance between Mg 2+ and the leaving water molecule R Mg−Ox for different force fields in combination with SPC/E water. (C) The snapshots show representative conformations in the two stable states: (i) Before exchange: Leaving water (shown in green) is part of the first and entering water (shown in blue) is part of the second hydration shell. (ii) After exchange: Leaving water is in the second hydration shell and the entering water molecule filled the void in the first hydration shell. The snapshots were taken using microMg(SPC/E).  73,74 is obtained from eq S3. The rate constant k is calculated from the number of transitions N for microMg and nanoMg and eq S3. The errors for N and k are obtained from block averaging. ΔF* is the free energy difference between the top and the first minimum ( Figure 4B). Values for parameters in TIP3P are taken from ref 26.  (Table 1) and from experiments. 88 The inset shows a representative snapshot of Mg 2+ in interaction with Cl − , including their first hydration shell. This snapshot was taken using microMg(SPC/E) and Cl − parameters from Smith−Dang. 54 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article oxygens. 77,78 Similar to previous work, 19,26,30,57,58 a dimethylphosphate (DMP) is used to mimic the phosphate backbone, the most important inner-sphere Mg 2+ binding sites on larger RNA molecules. 5,77,79−81 Here, the standard combination rules (eq 2 with λ σ,ε RNA = 1.0) significantly overestimate the binding affinity toward the phosphate oxygen reflecting that the Mg 2+ − RNA interactions are too attractive (data not shown). As before, this problem can be solved by increasing the effective diameter via λ σ RNA and simultaneously reducing the interaction energy via λ ε RNA ( Table 1). Note that, in this case, the scaling parameters for the different water models are similar but not identical since, unlike for chloride, the force field parameters of the RNA were not adjusted to the different water models. Figure 6 shows the binding affinity ΔG b 0 and binding distance R b for different water models and force field parameters from the literature. 19,20,22,30 Clearly, the bulkoptimized parameters (Mamatkulov−Netz(SPC/E), Allneŕ− Villa(TIP3P)) significantly overestimate ΔG b 0 and underestimate R b . The Panteva−York[m12-6-4] parameters 30 provide improvement. However, in the optimization the 4fold access of the phosphate oxygen binding site on the backbone compared to the nucleobase binding sites 77 was taken into account, which is only applicable in the context of the modular model used in the experimental study. The affinity is thus overrated. Also, the binding distance R b was not explicitly considered. With the optimized ion-RNA scaling factors for microMg and nanoMg, the results for ΔG b 0 and R b agree within errors with the experimental results for all 5 water models ( Figure 6, Table 5, and Tables S9 and S10).
The Mg 2+ force field parameters and the water model have a significant influence on the ligand exchange kinetics. 45 The effect is illustrated by the one-dimensional free-energy profiles as a function of the distance between Mg 2+ and the phosphate oxygen for different force fields and water models ( Figure  6B,C). The free energy profiles show two stable states (i) and (ii), corresponding to the inner-sphere conformation (direct contact between Mg 2+ and phosphate oxygen) and the watermediated outer-sphere conformation ( Figure 6D). The bulkoptimized parameters (Mamatkulov−Netz(SPC/E), Allneŕ− Villa(TIP3P)) significantly overestimate ΔG b 0 and underestimate R b as reflected by the shift and depth of the first minimum ( Figure 6B). For microMg and nanoMg, the innersphere and outer-sphere minima are identical, reflecting that the parameters have identical thermodynamic and structural properties. In addition, clear differences are observed for the barrier height, which corresponds in large part to the free energy cost of exchanging one water molecule from the first hydration shell with a water from the second shell. Force fields that underestimate the rate of water exchange (Mamatkulov− Netz(SPC/E), Allneŕ−Villa(TIP3P)) have high energetic barriers, Mg 2+ association/dissociation too slow, binding affinity too high, and binding distance too small ( Figure 6B). MicroMg and nanoMg provide significant improvement (Table  5). At the same time, the free energy barrier for nanoMg is up to 4 k B T lower compared to microMg depending on the water models used ( Figure 6B, Table 5).
Reproducing ΔG b 0 and R b is crucial to correctly describe the structure and thermodynamics of specifically bound Mg 2+ ions. Since microMg and nanoMg reproduce the properties of specifically bound ion as well as important bulk properties, both parameter sets are equally suited to calculate the distribution of Mg 2+ around biomolecules. In simulations that target the binding kinetics of the metal cations, the  49 or TIP4P-D. 50 The inset is a zoomed in view of the experimental area. The experimental values (gray) are taken from refs 77 and 78. (B, C) One dimensional free energy profiles along the distance between one of the nonbridging microMg parameter set should be used. To enhance the sampling of binding events, the nanoMg set can be used since the exchange kinetics is accelerated, thereby improving the sampling. However, the acceleration of the binding kinetics strongly depends on the water model used. For the largest speed-up, we recommend using nanoMg in combination with the TIP3P or SPC/E water model ( Figure 4A).

CONCLUSIONS
The distinct role of magnesium in biological systems has driven the development of force fields for molecular dynamics simulations. Despite considerable efforts, Mg 2+ force fields based on the 12−6 Lennard-Jones (LJ) potentials showed significant shortcomings in reproducing a broad range of structural, thermodynamic, and kinetic properties. Since Mg 2+ strongly polarizes its environment, the lack of polarizability and charge transfer effects in classical simulations likely causes such deviations. Recently, progress was made by an optimization procedure that implicitly accounts for polarizability. Considering an enlarged range of possible LJ parameters for Mg 2+ − water interactions and optimizing the Lorentz−Berthelot combination rule for Mg 2+ −Cl − and Mg 2+ −RNA interactions yielded simple and accurate force fields in TIP3P water. 26 However, a large variety of water models exist, and models with improved properties are increasingly used in biomolecular simulations. Since simulations rely on available water models and ion force fields that provide physically meaningful results when combined, 82 the transferability of the Mg 2+ parameters to other water models needs to be assessed and the parameters need to be optimized if necessary. Our results reveal that the Mg 2+ parameters developed in combination with TIP3P show limited transferability to SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D. While for the 3-site water models (TIP3P-fb and SPC/E) the deviations from the experimental solvation free energy are small, significant deviations are observed for the 4-site models (TIP4P/2005, TIP4P-Ew, and TIP4P-D).
To provide improvement, we have systematically developed improved Mg 2+ parameters. The optimized parameters (Table  1) reproduce the solvation free energy, the distance to oxygens in the first hydration shell, the hydration number, the activity coefficient derivative in MgCl 2 solutions, and the binding affinity and distance to the phosphate oxygens on nucleic acids (Tables 3, 5). In order to provide consistent and robust results for the solvation free energy and the activity derivative, Cl − was chosen as a reference ion and its parameters were optimized in a preceding step for each water model (Table 2). Here, an alternative and promising approach for future work is the simultaneous optimization of anion and cations parameters, which eliminates the necessity to select a reference ion. 83 In addition, including experimental binding affinities toward specific ion binding sites on biomolecules has turned out to be essential to capture the structure of specifically bound ions 26 as well as to reproduce the structural properties of large nucleic acids. 42 Two parameter sets are presented for each water model: The first sets, microMg, yield water exchange on the microsecond time scale in agreement with experimental results. 73,74 For the second sets, nanoMg, the parameters were chosen to yield the highest exchange rate possible while still reproducing all other thermodynamic and structural properties. Since the water models have different properties, 36 including different dielectric constants and diffusion coefficients, the maximum achievable exchange rate differs and ranges from 10 6 to 10 8 exchanges per second (Table 4).
In summary, the Mg 2+ parameters presented here provide simple, computationally efficient, and highly accurate models for the simulation of Mg 2+ ions in aqueous solutions of SPC/E, TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D water. For simulations targeting the kinetics of ion pairing or ion binding, the microMg is recommended. For simulations targeting the distribution of Mg 2+ around nucleic acids, proteins, or lipids the nanoMg is recommended to enhance the sampling of binding events. Here, the parameters for TIP3P 26 or SPC/E yield the largest enhancement.
Additional methods, including detailed lists on the different setups used in this study, and additional results: single-ion properties and their isosurfaces calculated to study the transferability of the original parameters microMg(TIP3P) and nanoMg(TIP3P) to other water models; additional Lennard-Jones interaction potentials and one-dimensional free energy profiles of the optimized parameters for TIP3P-fb, TIP4P/2005, TIP4P-Ew, and TIP4P-D that are not shown in the   Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article main text; two-dimensional free energy profiles of all 10 optimized Mg 2+ parameter sets; and binding affinity values obtained from alchemical transformation for all 10 optimized parameter sets (PDF)