Polarizable Force Field for Molecular Ions Based on the Classical Drude Oscillator
Development of accurate force field parameters for molecular ions in the context of a polarizable energy function based on the classical Drude oscillator is a crucial step towards an accurate polarizable model for modeling and simulations of biological macromolecules. Towards this goal we have undertaken a hierarchical approach in which force field parameter optimization is initially performed for small molecules for which experimental data exists that serve as building blocks of macromolecular systems. Small molecules representative of the ionic moieties of biological macromolecules include the cationic ammonium and methyl substituted ammonium derivatives, imidazolium, guanidinium and methylguanidinium and the anionic acetate, phenolate and alkanethiolates. In the present work, parameters for molecular ions in the context of the Drude polarizable force field are optimized and compared to results from the nonpolarizable additive CHARMM General Force Field (CGenFF). Electrostatic and Lennard-Jones parameters for the model compounds are developed in the context of the polarizable SWM4-NDP water model, with emphasis on assuring that the hydration free energies are consistent with previous reported parameters for atomic ions. The final parameters are shown to be in good agreement with the selected QM and experimental target data. Analysis of the structure of water around the ions reveal a substantial differences between the Drude and additive force fields indicating an important role of polarization in dictating the molecular details of aqueous solvation. The presented parameters represent the foundation for the charged functionalities in future generations of the Drude polarizable force field for biological macromolecules as well as for drug-like molecules.
1. Introduction<!>2. Computational Methods<!>Method 1<!>Method 2<!>3.1 Parametrization strategy<!>3.2 Electrostatic and LJ Parameters<!>3.3 Internal geometries and vibrational spectra<!>3.4 Free Energies of Hydration<!>3.5 Ion-Water coordination structure<!>4. Summary
<p>Classical molecular dynamics (MD) simulations based on empirical force fields (FFs) are widely used in a broad range of scientific fields such as chemistry, biochemistry, biophysics, engineering, and materials science. In particular, MD simulations may be used to investigate the structure and dynamics of complex biological processes at the atomic level under a wide variety of conditions. Examples include studies of ligand binding, enzymatic-reaction mechanisms, protein or RNA folding and protein-protein interactions.1</p><p>Central to the quality of the results from MD-based technologies is the accuracy of the underlying force field used to calculate the energies and forces acting on the system. As a result of many years of careful refinement, current additive biomolecular FFs have achieved a level of quality allowing them to be used to study dynamics, complex interactions and in pharmacological applications, among others.2 However, limitations are present in the additive models. A simple example is the inability of such models to correctly treat the distribution of atomic ions at air-water interfaces,3–5 an issue that was overcome by the explicit inclusion of polarizability in the force field. Thus, it is clear that further improvements in the accuracy of macromolecular FFs will likely require the inclusion of electrostatic polarization, as induced fields from the solvent environment, including ions as well as the molecules of interest themselves will affect the electrostatic interactions.6–10</p><p>Development of the Drude polarizable FF in the context of CHARMM,11, 12 subsequently extended to NAMD,13, 14 GROMACS,15, 16 and OpenMM,17 started in 2001 with optimization of the first water model (SWM4-DP).18 A later SWM4 water model was developed which included a negatively charged Drude particle, becoming the standard polarizable water model in the Drude polarizable FF (SWM4-NDP)19 and this has been followed by a six-point water model, SWM6-NDP,20 which showed improvements in a number of properties, including the treatment of water clusters, though at an additional computational cost over the SWM4-NDP model. Development of the Drude polarizable FF continued with parametrization of small molecules covering the functional groups encountered in biomolecules, including alkanes,21 alcohols,22 aromatics,23 among others.24–29 Development of parameters for N-methylacetamide (NMA) motivated extensions of the Drude energy function to include anisotropic polarizabilities and atom-based Thole scale factors,30 with the NMA model subsequently updated.31 More recently, Drude models for proteins,32 selected lipids33, 34 and carbohydrates28, 35–37 and DNA38, 39 have been presented. These are referred to as the Drude-2013 FF, with a second generation DNA force field,40, 41 along with the required parameters for RNA (Lemkul and MacKerell, Work in progress), referred to as Drude-2017 recently presented.</p><p>The molecular ions included in Drude-2013 FF were optimized to be consistent with published Drude parameters for atomic ions,42 though the details of the molecular ion parametrization were not presented. Moving ahead as part of the second generation Drude force field we have undertaken additional optimization of the Drude molecular ions. In the present study we provide a description of the parametrization of selected small molecule ions. Results on the performance of the FF in reproducing intramolecular geometries and vibrational spectra, interactions with water, dipole moments and molecular polarizability tensors based on quantum mechanical (QM) data, and experimental hydration free energies of the molecular ions are presented and discussed. The molecular ions considered in this study are listed in Table 1 along with their abbreviations. These were selected to cover the common ionic moieties in biological macromolecules at physiological pH, as well as the charged functionalities with pKa values close to physiological pH. Parameters for phosphate containing species will be presented as part of our ongoing efforts in the development of next generation nucleic acid and lipid force fields. Details of the Drude energy function and general parameter optimization strategies have been reported elsewhere43, 44 and are not presented in this report.</p><!><p>QM calculations were performed with the PSI445 and Gaussian 03 programs.46 Gas phase geometry optimizations were done at the MP2/aug-cc-pVDZ level. Vibrational frequencies were calculated using the optimized geometry with the same model chemistry. A scaling factor 0.959, obtained from the NIST Computational Chemistry Comparison and Benchmark Database,47 was used to account for limitations in the level of theory.48 A symbolic potential energy distribution (PED) analysis was performed as proposed by Pulay et al.49 using the MOLVIB module in CHARMM12 to obtain the vibrational frequencies and assignments. QM electrostatic potential (ESP) calculations were performed on the gas-phase MP2/aug-cc-pVDZ optimized geometries using the B3LYP hybrid functionals50–52 and the correlation-consistent double-ζ Dunning aug-cc-pVDZ basis set.53 Interaction energies of the model compounds with water were performed using the QM gas phase model compound and SWM4-NDP or TIP3P internal geometries at the MP2/cc-pVQZ model chemistry including the basis set superposition error correction of Boys and Bernardi.54 The calculations involved scans of the water to model compound distances in selected orientations (Figure 1 and Figure S1–S6) with the monomer geometries constrained to those obtained in the gas phase.</p><p>Empirical FF calculations were performed with the program CHARMM.12 All calculations were performed with the SWM4-NDP19 water model as this is the default water model in the Drude force field due, in part, to its computational efficiently over the SWM6 model.20 Additive force field calculations used the CHARMM General Force Field (CGenFF) parameters55 along with the CHARMM TIP3P water model.56 Gas phase optimizations included all possible nonbond interactions and were performed to a convergence threshold of 10−5 kcal·mol−1·Å−1 for the forces. Interactions with water used gas phase geometries in interaction orientations identical to those used in the QM calculations with the interaction energies obtained by taking the energy difference between the water-model complex and the water and model compound monomers, following relaxation of the Drude particles via a minimization by steepest-descent (SD) algorithm for 200 steps followed by an adopted-basis Newton-Raphson (ABNR) algorithm for 500 steps to a final gradient of 10−5 kcal·mol−1·Å−1 with the atomic positions restrained with a force constant of 107 kcal/mol/Å2.</p><p>The hydration free energy of molecular ions was calculated through alchemical free energy perturbation (FEP) simulations.57 Final hydration free energies were based on three individual calculations, as described below. In addition, to study the structure and dynamics of the hydrated ions to recover the full layered structure of the water, the cubic box of 250 water molecules and a single ion was submitted to a 10 ns simulation for each molecular ion for analysis of the radial distribution functions and 3D spatial probability distributions of water.</p><p>Simulations of the ions in aqueous solution were performed with cubic periodic boxes containing 250 water molecules and a single molecular ion restrained to the center of the simulation box through a force constant of 0.5 kcal/(mol·Å2) on the center of mass of all solute real atoms. Simulations were performed at 298.15 K and 1 atm pressure, corresponding to the experimental conditions. MD simulations of the Drude models were performed using the extended Lagrangian formalism implemented in CHARMM,18, 58 with a mass of 0.4 AMU subtracted from the polarizable atoms (ie. nonhydrogens) and assigned to the Drude particles. The time step used in all simulations was 1 fs for the Drude model and 2 fs for the additive model. Long-range electrostatic interactions were computed using particle mesh Ewald summation59 with a real space cutoff of 12 Å, Ewald splitting parameter of 0.34 Å−1, a grid spacing of ~1.0 Å, and a sixth-order interpolation of the charge to the grid was used. Lennard-Jones (LJ) 6–12 term parameters were truncated at 12 Å with switch smoothing60 used for the Drude and force switch smoothing61 for the additive applied between 10 and 12 Å. Simulations in the isothermal–isobaric ensemble used a Nosé–Hoover thermostat62, 63 and the modified Andersen–Hoover barostat of Martyna et al.64 The internal geometry of the SWM4-NDP water molecule and covalent bonds involving hydrogens were kept rigidly fixed using the SHAKE/Roll and RATTLE/Roll algorithm.65, 66 Gas-phase simulations were carried out using Langevin Dynamics with a friction constant of 5 ps−1.</p><p>The computation of hydration free energies of charged species is more complex than with neutral compounds because it must also account for the contribution arising from the Galvani potential Φ originating from the vacuum/water interface.67–69 Accordingly, we define the 'real' hydration free energy 
ΔGhyd0,real as the total reversible work to physically transfer a charged species from the gas phase to the bulk solvent phase (crossing the physical air/water interface).69–71</p><p>In equation 1, z is the charge of the species, F is the Faraday constant (23.06 kcal/mol/V), Φ is the phase potential of water relative to vacuum or the electrostatic Galvani potential at the vacuum/water interface, and 
ΔGhyd0,intr is referred to as the 'intrinsic' free energy. Although the real free energy 
ΔGhyd0,real is physically invariant, both Φ and 
ΔGhyd0,intr depend on the choice of mathematical convention. While Φ cannot be measured experimentally by a physical process, "internal" or "external" Galvani potentials may be defined mathematically unambiguously by specifying the convention used in the calculation (P- or M-sum).68 As discussed by Kastenholz and Hünenberger,72, 73 the convention for summing up the contributions of solvent charges to the electrostatic potential in the aqueous phase can be made based on point charges within entire solvent molecules (M-sum) or on the basis of individual point charges (P-sum). In the present work, the calculated free energies rely on periodic boundary conditions (PBC) and the particle mesh Ewald (PME) method carried out within the P-sum convention, and the implicit reference phase potential Φ of the liquid is the "internal" Galvani potential. The results from FEP simulations with PBC/PME are the intrinsic free energies, 
ΔGhyd0,intr, which then need to be corrected by zFΦ (with Φ calculated within the same P-sum convention) to yield the physically meaningful 
ΔGhyd0,real. For the polarizable SWM4-NDP water model, the P-sum internal Galvani potential Φ is equal to −545 mV19 (negative in the liquid phase relative to vacuum), giving rise to an energy shift of zFΦ being ±12.6 kcal·mol−1 for the monovalent cations/anions.19 For the additive TIP3 water model, Φ is equal to −500 mV such that zFΦ is ±11.5 kcal/mol.69</p><p>Following the previously reported Deng and Roux protocol,74 the total solvation free energy (ΔGtot) is calculated as a sum of the differences between gas and aqueous-phase system vdW repulsive (ΔGrep), vdW dispersive (ΔGdisp) and electrostatic (ΔGelec) terms as shown in equation 2.</p><p>In the free energy perturbation protocol, the gas and aqueous-phase simulations used the same coupling parameter schemes. The repulsive contribution (ΔGrep) was computed using a soft-core scheme74 with the staging parameter, s, set to 0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0. The electrostatic (ΔGelec) and dispersive (ΔGdisp) components were calculated using the standard linear coupling scheme with coupling parameters λ and ξ, varying from 0 to 1 with the increments of 0.1. For each window, a 50 ps equilibration and a 400 ps production simulation were performed. The weighted histogram analysis method (WHAM)75 was used for computing nonpolar contributions to the free energy difference while the electrostatic contributions were computed using thermodynamic integration (TI). A long-range correction (LRC)76 was then included in the calculation of 
ΔGhyd0,intr, where the LRC is the difference when the LJ solvent-solute interaction energies use cutoff schemes of 12 and 50 Å. Thus, the 
ΔGhyd0,real are computed from the sum of 
ΔGhyd0,intr, the entropic contribution of +1.9 kcal·mol−1,42, 69 and zFΦ. The real free energy 
ΔGhyd0,real is commonly referred to as "absolute" solvation free energy to draw a clear distinction with relative free energy between two different ions, although in physical terms, it corresponds to the reversible thermodynamic work to move the ion from vacuum to the bulk water phase. While such an absolute solvation free energy can be calculated unambiguously, in practice, only relative free energies of ions of the same charge or of neutral salts can actually be measured experimentally. In spite of this, experimental values of charged species are often reported as if they are absolute values, although they always depend on an arbitrary offset that may differ between the various reports. The offset is sometimes reported in terms of the absolute hydration free energy of a proton, but this information is not always provided explicitly. The confusion has been further compounded at times by the concepts of "intrinsic" and "real" solvation free energy (see discussion above). Accordingly, one cannot rely naïvely on absolute experimental values because there is no unambiguous reference scale for charged species. All experiments have to utilize some reference scale. Thus, special caution is needed when combining experimental values taken from different studies; a simple subtraction of two absolute values from two different studies does not yield a valid relative free energy since both studies may not have assumed the same implicit offset. This issue is important in force field development because the absolute hydration free energies of the different ions must be internally consistent, for example, to accurately account for the relative binding and ion-pairing affinities. To resolve this issue, let us assume that one knows from experiment the absolute experimental solvation free energy of a negatively charged compound with a charge of −1. There are essentially two methods to combine the experimental values from different studies to obtain the target data for the present optimization. These take advantage of hydration free energy of Cl− and Na+ and reference H+ used in the previous development of the Drude monoatomic ions42 (Table 2), which will be used to normalize the hydration free energies for the molecular ions, Y+/−, reported in the literature.77–81 Method 1 is used when the reference is a co-ion of the same charge (eg. Cl− as the reference in the case of CH3COO−, or Na+ in the case of NH4+). Method 2 is used when the reference is based on a counterion for forming a neutral salt (eg. Na+ or H+ in the case of CH3COO−, or Cl− in the case of NH4+). We note that Hofer and Hünenberger have recently developed a QM/MM based approach to calculate the accurate proton solvation free energy,82 though this result is not relevant to the present study.</p><!><p>Let us assume the molecule is CH3COO−. If the literature reports ΔGhyd(Cl−ref) and ΔGhyd(CH3COO−ref), then the difference of ΔGhyd(Cl−ref) and ΔGhyd(CH3COO−ref) will be used to determine the absolute scale. As ΔGhyd(Cl−Drude) is −78.4 kcal/mol,42 the target ΔGhyd,target value for CH3COO− will be −78.4 kcal/mol + ΔGhyd(CH3COO−ref) − ΔGhyd(Cl−ref). The same method also applies to cationic molecules. This calculation is described in equation 3.</p><!><p>Again let us assume the molecular ion is acetate and the literature reports ΔGhyd(Na+ref) and ΔGhyd(CH3COO−ref). Then, the sum of ΔGhyd(Na+ref) and ΔGhyd(CH3COO−ref) will be used to determine the absolute scale. As the ΔGhyd(Na+Drude) is −96.3 kcal/mol,42 the target ΔGhyd,target value for CH3COO− will be ΔGhyd(Na+ref) + ΔGhyd(CH3COO−ref) − (−96.3 kcal/mol). The same method also applies to cationic molecules. This calculation is described in equation 4.</p><p>In equations 3 and 4, ΔG(Y−/+ref) is the hydration free energies for the molecular ions, Y+/−, reported in the literature, ΔG(X−/+ref) is the hydration free energies for ions, X+/−, reported in the literature, ΔG(Y−/+
target) is the target hydration free energies for the molecular ions obtained from the Drude force field and ΔG(X−/+
Drude) is the hydration free energies for the ions obtained from Drude or reference H+ used in the development of the Drude monoatomic ions (Table 2).</p><p>Finally, while the shifting of the relative scales is simple in the case of monovalent ions, it requires more caution in the case of divalent ions. In this case, forming the neutral salt is the safest way to avoid mistakes. For example, assume that your compound (X2−) is negative with a charged of −2 with a reported hydration free energy of 
ΔGX2-exp in the literature. If the report provides the hydration free energy of another cation like Na+, 
ΔGNa+exp, then one can use this to determine the absolute scale of 
ΔGX2-FF for X2− by forming a neutral salt with two Na+,</p><p>This reasoning can be generalized to charged species of arbitrary valence. Based on the above correction schemes the target hydration free energies along with the various correction values are shown in Table 3. Analogous data for the additive CGenFF force field is presented in Table S1 of the supporting information.</p><!><p>Development of parameters for the Drude FF in CHARMM follows the previously described iterative methodology.20,23 For example, the initial electrostatic, internal and LJ parameters are used to calculate the free energies of hydration. Based on those values as well as considering the remaining target data, changes in, for example, the electrostatic parameters are performed followed by analysis of the remaining target data, re-evalution of the free energy of hydration and so on. Accordingly, when any of the parameters were adjusted all remaining target data was recalculated and additional adjustments of the parameters made was required. The following results are based on the final set of parameters, which will be accessible from the MacKerell laboratory web site (http://mackerell.umaryland.edu/CHARMM_ff_params.html), with future releases of the program CHARMM and from the Drude Prepper module in the CHARMM-GUI.83</p><p>Electrostatic parameters in the polarizable Drude model include partial atomic charges, atomic polarizabilities, and atom-based Thole scale factors.84, 85 In addition, selected hydrogen bond acceptor atoms include anisotropic corrections to their polarizabilities as well as virtual particles representative of lone pairs (LP).85 Initial determination of the electrostatic parameters is based on a series of ESPs surrounding the model compound determined using QM density functional theory on a set of specified grid points, differing by the presence and location of a perturbing charge in the environment surrounding the molecule.10, 84 Electrostatic parameters are then fitted to minimize the difference between the QM and Drude ESP maps. Initial guesses for the electrostatic parameters were taken from the neutral analogues of each molecular ion. Initial polarizabilities were the atomic hybrid polarizability values published by Miller.86 Subsequently, iterative optimization of the electrostatic parameters, in conjunction of the optimization of the remaining parameters (see below), involved manual adjustments targeting dipole moments, molecular polarizabilities, interactions with water and hydration free energies.</p><p>Optimization of the remaining components of the FF, namely the internal and LJ terms, were adjusted as described previously for the additive CHARMM FF.87, 88 Internal parameters were optimized by targeting experimental and QM data that included geometries, vibrational spectra and conformational energies. Optimization of LJ parameters, which model vdW forces in the Drude FF, rely on reproduction of experimental hydration free energies (Table 3). In addition, QM minimum interaction energies and distances of the model compounds with water are used as LJ target data (Table S2).</p><!><p>Partial atomic charges, atomic polarizabilities, and Thole scale factors were initially optimized targeting QM ESPs. Target molecular polarizabilities were not scaled in the cations and scaled by 0.6 to 0.7 for the anions, a value similar to that used for the SWM4-NDP water model19 and for the Drude alcohols.22 In ACET, MES and ETS the positions of LPs were manually adjusted to improve agreement between the Drude FF and QM calculated interactions with water. With PHET the quality of the model in the absence of the LPs was deemed adequate based on reproduction of QM interactions with water and the free energy of hydration such that LPs are not included in that model.</p><p>Initial values for the LJ parameters were based on neutral analogs of the molecular ions with adjustments first based on the reproduction of QM minimum interaction energies and distances between the ions and water. The interaction orientations between water and the protonated amines are shown in Figure 1 along with the QM, Drude and additive CGenFF interaction energy surfaces. Similar images and data for all the molecules considered in this study are shown in Figures S1 to S6 of the supporting information. It should be noted that the interactions were treated individually (e.g. as monohydrates). A summary of the model compound-water interaction energy results is shown in Table 4, with the data for the individual interactions presented in Table S2 of the supporting information. Overall, good agreement between QM and MM reported minimum interaction distance and energies are observed. The minimum interaction distances are systematically longer in the Drude model and the interaction energies are slightly less favorable. The average distance and energy differences are better with the additive CGenFF force field though the variation is larger, which is also seen in the larger AUE and RMSD values with the additive model. The overall level of agreement with the Drude model, whose parameters were optimized as part of the present study, was partially compromised due to the LJ parameters also being optimized to reproduce the hydration free energies. Results of the hydration free energy analysis are presented below.</p><p>Dipole moments of the model compounds were also considered during optimization of the electrostatic parameters. Presented in Table 5 are the QM dipole moments along with Drude and additive values. In general, the Drude values are in good agreement with the QM data. The differences with the additive model are systematically larger as expected given that the QM data was not considered when optimizing that model. We note that the impact of the dipole moment on interactions with the environment is less important with charged species vs. polar, neutral compounds as the monopole on the ions dominates such interactions.</p><p>An important quality of a polarizable FF is to accurately represent the molecular polarizability tensors, which determine the magnitude of the induced dipoles in the presence of an electric field. Thus, the polarizability tensors of all model compounds were computed and compared to QM values. Only the diagonal elements of the polarizability tensors are non-zero due to the alignment of each molecular ion along the axis of inertia. The components of the polarizability tensor along each coordinate axis were computed from finite-difference of the respective dipole moment component, x, y or z, in the presence and absence of an applied electrical field of 1010 V/m, respectively, parallel to the x, y, or z directions. Comparison of the QM reference values and equivalent quantities calculated with the polarizable Drude model is presented in Table 6. Overall, the discrepancies between the polarizable Drude model and QM values are small, evidencing the ability of the polarizable model to reproduce subtle aspects of the electrostatic behavior of the molecules under study. With polar neutral species in the Drude force field it has been shown that the polarizabilities need to be scaled by values down to 0.6 in order to reproduce the dielectric constants of pure solvents.27 In the present study the interactions with water, reproduction of dipole moments and reproduction of experimental hydration free energies were used to optimize electrostatic parameters. This approach revealed that there was no need to scale the polarizabilities of ammonium and the methylated derivatives of ammonium as good reproduction of the hydration free energies was achieved (see below). However, polarizability scaling was required with guanidinium, methylguanidinium and the anions. The scale factor of 0.85/0.86 used with GUAN, MGUAN and ACET falls in the range of scaling values obtained with the polar neutral species. With PHET and the sulfur containing compounds scale factors of 0.6 were required. This value is consistent with that determined for neutral sulfur containing compounds,27 though the value of 0.64 for PHET is lower the values in the range of 0.7 used for neutral alcohols.22</p><!><p>Parameters for the bonds, valence angles and dihedral angles were optimized to reproduce QM optimized geometries and vibrational spectra of the studied molecular ions. Internal parameters for the methyl groups in ACET and the methylated derivatives of ammonium and guanidinium were transferred from the alkanes.21 Target data for the geometries included the optimized structure of the compounds at the MP2/aug-cc-pVDZ level. RMS differences of the bonds and angles with respect to the target data are shown in Table 7. The absolute values of all bond distances and angles and comparison with the QM data are shown in Table S3 of the Supporting Information. The agreement between calculated and target values is in general good, falling within the targeted differences of 0.02 Å for bonds and ~2° for angles, a level of agreement that is commonly used for the CHARMM FFs.2, 87, 89</p><p>One of the distinctive features of the CHARMM additive and Drude force fields is the use of reproduction of QM vibrational frequencies within the harmonic approximation to optimize the empirical force constants. Emphasis was placed on the magnitude of the frequencies along with the reproduction of the assignments. All vibrational frequencies and assignments for the model compounds are shown in Table S4 of the Supporting Information. Overall, the agreement for the lowest frequency vibrational modes between the QM and MM models is quite good, indicating that the Drude force field will satisfactorily treat structural variations of the ions that occur during MD simulations.</p><!><p>Central to the quality of a force field for ions is its accuracy with respect to estimation of hydration free energies, with those estimates based on common offset of the free energies as discussed above. Proper offset of the hydration free energies is essential to assure that the relative hydration free energies of all ions in the force field are correct. Statistical analysis of the results for the presented Drude model and CGenFF are shown in Table 8 with the results for the individual ions shown in Table 9. With the exception of NH4, the calculated hydration free energies for the Drude model are within 2 kcal/mol of the experimental target values. Substantially larger differences occur with CGenFF leading to the larger AVG, AUE and RMSD values with the additive model. This larger discrepancy is expected due to the additive force field parameters primarily being optimized based on QM interactions with water, where the QM calculations were performed at the HF/6-31G* model chemistry, consistent with that used for CGenFF and the entire C36 additive force field.2, 55 However, the ordering of the hydration free energies with the additive force field is quite acceptable, as is that for the Drude model as seen in Figure 2 and in the R2 values in Table 8. This indicates the potential utility of both force fields in predicting the relative hydration energies of the molecular ions in different chemical contexts.</p><p>With the Drude force field the discrepancy with NH4 (ammonium) is based on the optimization focusing on the methylated ammonium analogs (NC1 through NC4), such that the same LJ parameters are used for NH4 and NC1 and for NC2 through NC4. Specific optimization of LJ parameters for NC1 (ε = −0.025, Rmin/2 = 1.55) yielded a calculated ΔGhydreal of −72.6 kcal/mol, in excellent agreement with target value in Table 9. An improved ΔGhydreal for NH4 of −81.3 kcal/mol is obtained if LJ parameters (ε = −0.025, Rmin/2 = 1.76) specifically optimized for NH4 are used. Concerning NC2 through NC4, the ΔGhydreal values of NC2 and NC3 are not favorable enough while NC4 is too favorable, a compromise that was deemed acceptable in the context of a single N atom type for the three compounds. With GUAN and the methyl analog, MGUAN, emphasis was placed on the methylated species as this is a better model for the side chain of Arginine. However, the discrepancy with GUAN is less than 2.0 kcal/mol, which is considered acceptable in the context of the same atom types. Overall, the presented Drude parameters for the molecular ions yield calculated hydration free energies that are in excellent agreement with experiment for the biomolecular-related ions included in this study.</p><!><p>Given the difference in the underlying potential energy function in the Drude polarizable and CGenFF additive force fields, radial distribution functions (RDFs) for the ions with water were calculated to determine if difference in the water structure is present between the two models. RDFs were extracted from 10 ns simulations of the single ions immersed in a box of water, with the results presented in Figures S7–S8 for molecular cations and anions, respectively. The resulting RDFs, shown in Figures S7–S8 show the overall structures of water around the ions to be similar for the two models. However, differences are present and, interestingly, consistent trends in those difference are largely not present. For example with NC1, the first peak is higher with CGenFF while the opposite is true with NC2, NC3 and NC4. With GUAN and MGUAN, the RDFs are smooth at 3~4 Å with the Drude model (Figure S7). With the anions, subtle changes occur in the first peak and a trend occurs where the second peak is systematically larger with the Drude force field (Figure S8).</p><p>Given the subtle differences in the RDFs, more detailed analysis was undertaken involving the calculation of the 3D spatial probability distributions of water around the ions. Results in Figures 3–5 and S9–S13 do show obvious differences in most of the systems. For example, with NH4 and the methylated ammonium analogs (NC1 through NC4) the water distributions are less localized in the Drude model (Figure 3 and S9–S11). With GUAN and MGUAN, small water distributions occur above and below the planes of the molecules in the Drude model but are absent in additive model. With ACET, the water around the oxygens in the Drude model are more distributed and correspond to the lone pair directions (Figure 4a–b), which is likely associated with the use of lone pairs on the oxygens for better describing the directional hydrogen bonds. With PHET, the water distributions are more obvious on the external side of oxygen along the C-O bond in the Drude model (Figure 4c–d), showing that they are more favored around the negatively charged oxygen than in the additive model. Thus, while the two models yield generally similar solvation structures based on RDFs, differences are evident from the 3D spatial water distributions indicating the impact of the Drude model on solvation structure. Future efforts that compare calculated solvation structures with scattering studies of salt solutions are required to better evaluate the quality of the force fields with respect to their modeling of water structure.</p><!><p>Presented are classical Drude polarizable FF parameters of individual molecular ions that represent building blocks of larger molecules of biological interest such as peptides, lipids and carbohydrates. Contrary to neutral compounds, where extensive condensed phase properties are available, and against which newly developed parameters can be tested, for molecular ions such data is scarce. Moreover, when it is available, as with experimental Gibbs free energies of hydration, there are significant assumptions in deriving such data.94 The fact that ions cannot exist alone or undergo phase transfer without a counter-ion is the primary complication in obtaining this data. To overcome this issue, we offset the hydration free energies in a consistent fashion thereby assuring that the relative hydration free energies are consistent throughout the force field. In the case of the Drude force field this was done to assure consistency with previously reported atomic ions.42 Similar corrections were applied to the additive CGenFF force field, with results from the model presented for comparison with the developed Drude molecular ions.</p><p>With respect to model development, special care was taken in optimizing the electrostatic and LJ parameters to assure a proper balance of local and condensed phase properties. This included assuring that the balance of interactions with water in different model compound-water orientations were satisfactorily reproduced as were dipole moments and molecular polarizabilities. Care was also taken to assure that the intramolecular geometries and the vibrational frequency distributions were in good agreement with QM target data. Finally, the Drude model is shown to reproduce hydration free energies at a high level of accuracy, with improvements over the additive CGenFF evident, though it is anticipated that improvements in the additive model could be made with optimization of the model undertaken. These results indicate that the Drude polarizable model of molecular ions is of satisfactory quality for use in molecular modeling and simulation studies.</p>