Is the Solution Activity Derivative Sufficient to ... - ACS Publications

Subscriber access provided by University of Newcastle, Australia

Article

Is the solution activity derivative sufficient to parameterize ionion interactions in aqueous solution?-ions for TIP5P water Vahid Satarifard, Sadra Kashefolgheta, Ana Vila Verde, and Andrea Grafmueller J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b01229 • Publication Date (Web): 10 Apr 2017 Downloaded from http://pubs.acs.org on April 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Is the solution activity derivative sufficient to parameterize ion-ion interactions?-ions for TIP5P water Vahid Satarifard,†,‡ Sadra Kashefolgheta,†,‡ Ana Vila Verde,∗,† and Andrea Grafmüller∗,† Theory and Bio-Systems, Max Planck Institute of Colloids and Interfaces, Potsdam, Germany E-mail: [email protected]; [email protected]

1

Abstract

Biomolecular processes involve hydrated ions, and thus molecular simulations of such processes require accurate force-field parameters for these ions. In the best force-fields, both ion-water and anion-cation interactions are explicitly parameterized. First, the ion Lennard-Jones parameters are optimized to reproduce, e.g., single ion solvation free energies; then ion-pair interactions are adjusted to match experimental activity or activity derivatives. Here we apply this approach to derive optimized parameters for concentrated NaCl, KCl, MgCl2 and CaCl2 salt solutions, to be used with the TIP5P water model. These parameters are of interest because of a number of desirable properties of the TIP5P water model, especially for the simulation of carbohydrates. The results show, that this state of the art approach is insufficient, because the activity derivative often reaches ∗ To

whom correspondence should be addressed and Bio-Systems, Max Planck Institute of Colloids and Interfaces, Potsdam, Germany ‡ Contributed equally to this work † Theory

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a plateau near the target experimental value, for a wide range of parameter values. The plateau emerges from the interconversion between different types of ion pairs, so parameters leading to equally good agreement with the target solution activity or activity derivative yield very different solution structures. To resolve this indetermination, a second target property, such as the experimentally determined ion-ion coordination number, is required to uniquely determine anion-cation interactions. Simulations show that combining activity derivatives and coordination number as experimental target properties to parameterize ion-ion interactions, is a powerful method for reliable ion-water force field parameterization, and gives insight into the concentration of contact or solvent shared ion pairs in a wide range of salt concentrations. For the alkali and halide ions Li+ , Rb+ , Cs+ , F− , Br− , I− , we present ion-water parameters appropriate at infinite dilution only.

2

Introduction

All biological processes and many environmental and geological ones take place in aqueous solutions which contain inorganic ions. The ions play important roles in many relevant processes. For example, structure and function of biomolecules such as proteins, nucleic acids and lipid bilayer membranes are regulated or affected by interactions with ions, 1–3 and chemical reactions in aerosol droplets critically depend on their salt concentration. 4 Molecular modeling can give unprecedented insight into these processes, because it can directly probe atomic or molecular degrees of freedom and allows for the calculation of ensemble- or time-averaged observables which can be directly compared with experiment. At present, most systems of interest can only be simulated by approximating interactions between atoms or ions via classical potential energy expressions, termed force fields. Only this way can the large length-scales and long time-scales that characterize these processes be accessed. Quantum mechanical calculations, which contain fewer approximations, are still limited to systems of only tens of atoms and picosecond time scales. An adequate force field must correctly reflect the balance between water-water, water-solute and solute-solute interactions. Moreover, for many applications it is indispensable that the force


Page 2 of 34

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


field captures the intrinsic directionality that characterizes short-range interactions in hydrogenbonded systems; in contrast, explicit consideration of polarizability is often unnecessary and fixedcharge models perform well under many conditions. 5–11 Short-range directionality naturally arises in all-atom models, where each atom is explicitly represented, and many all-atom, fixed-charge force fields have been developed for biomolecules, water and ions. 10 Despite several decades of development of all-atom force fields, optimizing the balance between all interactions remains an on-going effort which, because of its magnitude, proceeds stepwise and iteratively. Many force field were originally developed to reproduce properties that mostly reflect selfinteractions and interactions with solvent molecules, such as density and enthalpy or free energy of solvation, as well as structural properties of biomolecules. However, matching single molecule properties does not reflect the balance of solute and solvent interactions, which may require additional tuning. 12–14 Furthermore, biomolecular force field parameters are typically developed using a particular water model. Replacing that water model by a more recent one, which offers better descriptions of the structure and the dynamics of water, 10 is sometimes acceptable without further modifications, 15 using careful consideration. However, to take full advantage of the improved water model, further tuning of the protein and ion force fields 16–19 is indispensable. Of the various classical, all-atom, water models with functional forms that are in principle compatible with the most popular force fields for biomolecules, 10 the TIP5P model is of particular interest because it reproduces a number of properties of water (e.g., the liquid density, temperature of maximum density and diffusivity as a function of temperature and pressure) substantially better than the TIP3P model used to develop the widely used AMBER and CHARMM force fields for proteins. 8,9 Recent studies show that the TIP5P water model

20

is of particular interest for simula-

tions using the GLYCAM06h force field for sugars: this water model results in particularly good agreement with experimental data such as diffusion coefficients, free energies of hydration and osmotic pressure for solutions of disaccharides. 11,14 However, widespread usage of the TIP5P water model in simulations of biomolecular systems is currently hindered by the absence of a force field



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for inorganic ions that has been optimized for this water model: to the best of our knowledge, the closest parameter set that has been presented was developed for the TIP5P/Ew water model and is limited to halide ions and sodium. 21 Here we specifically address the need for parameters for the alkali, alkali-earth and halide ions that are compatible with the TIP5P water model. Our aim is to develop a force field for ions that adequately reflects the balance of ion-water and ion-ion interactions. To achieve this aim, we separately optimize the Lennard-Jones parameters for ion-water and anion-cation interactions to reproduce experimental observables that reflect those interactions. Possible targets for parameterization are for instance the free energy or the entropy of solvation, 18,19,22–25 which intrinsically depend on ion-water interactions, and the solution activities, activity derivatives or radial distribution functions at finite salt concentrations 26–32 which necessarily depend on ion-ion interactions. Here we chose to target the solvation free energy of ions at infinite dilution to develop parameters for Li+ , Na+ , K+ , Rb+ , Cs+ , F− , Cl− , Br− , I− , Mg2+ and Ca2+ ions. In addition, for chosen salts (NaCl, KCl, MgCl2 and CaCl2 ) ion-pair interactions are adjusted based on the molar activity derivative at 0.5 m concentration where m stands for molal concentration. Both properties can be accurately determined in both experiment and simulation. These target properties have been used by others to develop ion force field parameters for other water models, 19,23,28,31,33 but never for TIP5P. Our results show that adjusting ion-pair interactions based solely on the activity derivative – a state-of-the-art approach to optimize these interactions – is often insufficient, and that to obtain parameters that yields the correct solution structures requires that other experimental properties reflecting ion-ion interactions are also used for parameterization. Here, we chose the ion-ion coordination numbers, as an additional criterion where required, as this quantity can uniquely distinguish between types of ion pairs, and thus the correct solution structure. In addition, the approach followed here improves the previous applications in several aspects. Firstly, using solvation free energies for parameterization has in the past been problematic, 29 because of the large uncertainty associated with experimental values of single-ion solvation free energies, and the dramatic dependency of the calculated single ion solvation free energies on the


Page 4 of 34

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulation methodology, mostly associated with the treatment of electrostatic interactions in MD simulations. Therefore, often the difference of intrinsic free energies of hydration relative to a reference ion is used as the target property 23,28 to parameterize ion-water interactions. In those studies, the parameters for the reference ion – Cl− – were obtained by other means. This leaves the uncertainty related to the hydration free energy of the reference ion undetermined. In contrast, here we use the actual intrinsic free energies of hydration, and apply all the necessary correction terms to make the calculated free energies directly comparable to experiment. We minimize any systematic bias of the results that could derive from the uncertainty in the reference free energy value of the proton by using the recommended value following a thoroughly review of the literature. 29 Though other works have also used the actual hydration free energies as the target for parameterization, some correction terms – e.g. those related to P-summation 29 –, to the calculated intrinsic hydration free energies are often not applied. 19 Secondly, in the calculation of the activity derivatives used for the optimization of the ionion parameters we apply the recently developed version of the Kirkwood-Buff integrals that is appropriate for simulations in closed systems, 34 in contrast to many other works, which apply the original integrals developed for open systems. 27,28,31,33 Finally, we examine the issue of force field transferability in the context of the TIP5P water model. 20 This issue is important, because many more force fields are available for both charged and non-charged species in simpler water models than in TIP5P water, but TIP5P has substantial advantages over other models, as discussed above. Our results demonstrate that force fields for charged or partially charged species developed for other water models will in general not be transferable for TIP5P water. Transferability is only possible for species with low charge density. To evaluate the quality of the models, we investigate the ion-water radial distribution function (RDF) and the number of water molecules in the first hydration shell of the ions, the single ion and salt diffusivities and the residence times of water molecules in the first and second hydration shells of the ions, and compare the calculated values with experiment where possible.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

Methods

We perform four different types of simulations: (i) free energy perturbation simulations of single ions in water, to determine the free energies of solvation of the ions; (ii) molecular dynamics simulations of single ions in water, used to calculate single ion properties for each of the ions considered in this work; (iii) simulations of NaCl, KCl, MgCl2 and CaCl2 solutions at 0.5 m concentration, used to calculate the activity derivatives of these solutions; and simulations of NaCl at 6.18 m and KCl at 4.56 m, to calculate anion-cation radial distribution functions and coordination numbers. The TIP5P water model 20 is used throughout this work. Simulations are performed using the molecular dynamics package Gromacs version 5.0. 35 The initial configurations for the free energy perturbation simulations and the simulations of salt solutions (i and iii) are assembled with the editconf module in Gromacs and the Packmol package. 36 The initial configuration used for the molecular dynamics simulations of single ions (ii) is the final configuration of the fully coupled state of the free energy perturbation simulations. The simulation boxes are cubic and periodic boundary conditions are applied in all directions. Energy minimization is performed by the steepest descent algorithm. Non-bonded potentials are cut off at 1.4 nm; electrostatic interactions beyond that cutoff distance are calculated using the particle mesh Ewald (PME) method. 37 The leap-frog algorithm 38 is employed to integrate the equations of motion. Simulations are performed at 298.15 K; the temperature is kept constant using the Nose-Hoover scheme 39,40 with a time constant of 0.1 ps. Simulations in the NPT ensemble are done with a target pressure of 1 bar, using the Parrinello-Rahman 41 approach with a time constant of 0.5 ps. All systems are equilibrated for 100 ps in the NVT ensemble followed by 100 ps in the NPT ensemble. The production simulations for free energy calculations are carried out for 1 ns in the NPT ensemble. Production runs of the molecular dynamics simulations of single ions in water are done in the NVT ensemble and last 25 ns, with configurations saved every 1 ps. The NVT ensemble is chosen in this case to avoid unnecessary perturbations introduced by the barostat, which might affect the observables determined from these simulations; prior work has shown that using the Nose Hoover thermostat yields transport properties indistinguishable from 6 ACS Paragon Plus Environment

Page 6 of 34

Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


those calculated in the NVE ensemble. 42 The molar activity derivatives are calculated from the radial distribution functions using the Kirkwood-Buff (KB) theory modified for the canonical ensemble. 43,44 More details about the method can be found in the supporting information (SI). For the calculation of the activity derivatives, salt solutions are equilibrated in the NPT ensemble for 10 ns. For each salt, we select three distinct configurations from the equilibration run as the starting configurations for 3 independent NVT production runs. The starting configurations have a box size close to the average box size during equilibration. Each production run lasts 50 ns, totaling 150 ns simulation time for each salt, with configurations saved every 2 ps. To ensure converged results, values obtained from 50 ns trajectory blocks are compared (see table S2).

3.1 3.1.1

Calculation of solvation free energies Definitions

The different contributions to the solvation free energies are briefly summarize here. A full discussion can be found in the SI. The real free energy of solvation of single ions corresponds to the process of bringing the ion from its standard state in vacuum to its standard state in the liquid phase and thus includes two contributions: a surface term and an intrinsic term which arises from all ion-solvent interactions in the bulk of the solvent. To avoid uncertainties associated with deconvoluting the contributions to experimentally measured free energies, single-ion solvation free energies are often tabulated in relative, so-called conventional scales, instead of absolute (real or intrinsic) ones. Conventional scales list the free energy values relative to a chosen reference ion, typically the proton. An exhaustive review 29 indicates that the best estimate of the intrinsic solva+

+

H , of the proton is that by Tissandier et al., 45 ∆G H = −1104.5 kJ·mol−1 . tion free energy, ∆Gslv slv

This is the value used in this study, to convert conventional to real scales. The experimental target free energy values listed below are standard intrinsic solvation free energies, i.e. they reflect the transfer of an ion from its standard state in the ideal gas phase to its standard state in solution. The standard state in solution is an imaginary 1 M solution with the 7 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

properties of an infinitely dilute system; the standard state in the gas phase is that of an ideal gas at 1 bar and 298.15 K.

3.1.2

Simulation approach

We calculate the intrinsic Gibbs solvation free energy of each ion, using the Bennett Acceptance Ratio (BAR) 46 method implemented in the g_bar module of Gromacs. The ion, initially in the vacuum state, i.e. fully decoupled from the solvent, is inserted into the solvent box, by scaling the Lennard-Jones and electrostatic interactions of the ion with water, using the path coordinates λ LJ ∈ [0, 1] and λCoul ∈ [0, 1], respectively, and using LJ soft-core potentials. The free energy associated with cavity formation, ∆Gcav , was calculated by dividing each path coordinate into 21 steps

1

raw , was calculated by dividing and the charging free energy, ∆Gchg

each path coordinate into 21 equally spaced steps to ensure ensemble overlap. At each value of λ LJ and λCoul , the system is equilibrated for 200 ps and statistics are collected for 1 ns, with a frequency of 1 ps. To ensure that this is sufficient to produce reliable results, trajectories were divided into blocks. The simulated solvation free energies from the different time segments and using a different thermostat (listed in the SI, Table S1) are in good agreement, indicating that the sampling is sufficient for converged results and the results are robust with respect to the system setup. The solvation free energy calculated this way does not include a contribution of crossing the macroscopic gas-liquid interface and is thus intrinsic.

3.1.3

Correction Terms

To address the method dependence of calculated free energy values, Hünenberger et al., after a decade-long endeavor, have shown that the raw intrinsic free energy calculated from simulation can be corrected a posteriori to compensate for methodology-dependent artifacts. 29,47–50 They have presented a complete set of correction terms, derived through analytical models of continuum electrostatics, that must be added to the calculated observables to obtain single ion solvation prop1λ

LJ

= 0.00, 0.06, 0.12, 0.18, 0.24, 0.30, 0.36, 0.42, 0.46, 0.50, 0.52, 0.54, 0.56, 0.58, 0.60, 0.64, 0.68, 0.72, 0.76, 0.80, 1.00


Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


erties that are methodology-independent and can be compared with experiment. 48,49 To remain consistent with that work, we use their terminology for the different terms. 29,49 Here a type-C1 correction term was applied to correct for effects of the P-summation scheme, and a correction term of type-D to minimize the errors arising from the inaccuracy in the permittivity of the solvent model. A type-B correction term to eliminate the contribution of the interaction of the ion with its periodic images was not required, as the GROMACS package 35 already excludes the self-energy term in the decoupled state, and its magnitude in the coupled state is only of order 1-2 kJ·mol−1 , because the dielectric constant of the solvent is high. In addition, the experimental solvation free energy values include an entropic contribution that arises from changes in volume between the standard states in the gas and the liquid phase, which requires a standard state correction term ∆Gstd = 7.95 kJ·mol−1 . A full description of the above correction terms terms is given in the SI; a full account of all terms has been discussed in detail elsewhere. 29,49

3.2

Optimization of Lennard-Jones parameters

The non-bonded Lennard-Jones (LJ) potential form (6) !

(12)

V (rij ) =

(12)

is used, where Cij

(6)

and Cij

Cij

rij12

−

Cij

rij6

(1)

are LJ potential parameters describing the interactions between

species i and j, and the superscripts between parenthesis are identifiers and not mathematical exponents. The r12 dependence of short range interactions is arbitrary, but the r6 dependence of long range interactions is physically motivated. 51 We take advantage of this physical basis to obtain the homoionic C (6) coefficients for each ion species using the Slater-Kirkwood equation for homoionic dispersion interactions 50,52 (6)

Cii = CSK α˜ 3/2 n1/2


(2)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 34

where CSK = 23.97 kJ·mol−1 · nm3/2 , α˜ = α/(4πe0 ) is the static dipole polarizability volume of the ion, α is the static dipole polarizability, and n is the effective number of electrons for the isoelectronic noble gas. The static polarizability values for the cations are taken from the gas phase polarizabilities; 53–55 for the anions, the values for ion polarizabilities in solution are used. 56 The effective number of electrons are estimated based on the noble gas polarizabilities and dispersion coefficients. 57,58 The C (12) coefficients do not have any physical significance and can be used as tuning parameters to reproduce experimental observables. Here, we parameterize the C (12) coefficient associated (12)

with ion-water interactions (CiO ) to reproduce the experimental values of the intrinsic solvation free energy, and the C (12) coefficient associated with anion-cation interactions for some of the salts to reproduce the activity derivative of solutions at 0.5 m concentration. (n)

The Lennard-Jones coefficients, Cij

(where n = 6, 12) for interactions between different

species that are not explicitly parameterized are obtained based on the geometric average of the (n)

(n)

corresponding self-interaction parameters, 20 Cii and Cjj , i.e., (n) (n) (n) 1/2 Cij = Cii Cjj (12)

Equation 3 is also used to obtain the Cii

(3)

term for the interaction between two identical ions

(n)

(n)

from the C (12) coefficients for water-water (COO ) and ion-water (CiO ) interactions. To reproduce the experimental solvation free energies in simulation, only the repulsion co(12)

efficient CiO

(12)

of ion-water LJ interactions is tuned. The initial choice of CiO

parameters for

monovalent ions is taken from reference 49. For the divalent ions, Mg2+ and Ca2+ , the initial values are taken from reference 28. We perform multiple free energy simulations for each ion, for (12)

(12)

a range of CiO values. The optimal CiO for each ion is that which yields an intrinsic free energy (12)

of solvation identical to experiment. The initial choice of the anion-cation Cij

parameters that

are optimized against the activity derivative is the geometric average obtained with Equation 3, of the associated homoionic parameters for ions i and j.


Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12)

Other optimization strategies could have been used, namely scanning both Cij

(6)

and Cij pa-

rameters to simultaneously satisfy both the hydration free energy and the activity derivatives. 28 We opted not to do so because that approach requires the simultaneous optimization of parameters for the different ion pairs, and thus a much larger number of simulations, which would require reduced simulation times and lower statistical accuracy in the hydration free energies and the activity derivatives. Our chosen approach guarantees that the calculated observables have low associated statistical uncertainty, and ensures that the target properties are satisfied within the estimated limits of experimental uncertainty. Even though our approach appears less general than a 2-dimensional scan of the Lennard-Jones parameters, because we rely on a phenomenological theory 50,52 to es(6)

timate the Cij , the results below demonstrate that it performs well: we reproduce not only the target parameterization properties, but also other properties of dilute and concentrated electrolyte solutions.

3.3

Analysis of static and dynamic properties

We evaluate the quality of our optimized ion-water parameters by calculating a number of properties to be compared to experimental data: the ion-Ow RDF and corresponding coordination number (Cn ) of water oxygen atoms as a function of the distance from the ion, the self-diffusion coefficient 1 ) of water oxygens in the first hydration (D) of ions and selected salts and the residence time (τres

shell of the ions. The experimental range/values for the RDF (ref. 59), Cn (ref. 60), D (ref. 58) 1 (ref. 61) were obtained from literature and are tabulated below. and τres

The self-diffusion coefficient, Di , of ions, is determined from the well-known Einstein relation

limt→∞ |ri (t) − ri (0)|2 = 6Di t, 62 where ri denotes the position of an ion at the indicated time and the average is over all time origins. The MSD is averaged over 25 ns and the linear fit used to calculate the diffusion coefficient is done from 1 to 1000 ps of the MSD ensuring that the MSD is within the diffusive regime and has low statistical noise. The diffusion coefficient for a salt, Dsalt ,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

is calculated from the diffusion coefficients of its cations, D+ , and anions, D− , by the relation: 58

Dsalt = (z+ + |z− |) D+ D− / (z+ D+ + |z− | D− )

(4)

where z indicates the charge number (z = q I /|e|) of the cation or anion. To calculate the residence time of water molecules in the first hydration shell of the ions, we use the stable states picture of chemical reactions proposed by Laage and Hynes, 63 which has already been successfully used to estimate residence times of water around ions. 22 Calculation details can be found in the SI.

4

Results and discussion

4.1

Optimizing ion-water interaction parameters (12)

As mentioned above, to optimize the CiO parameters we perform free energy perturbation simu(12)

lations for each ion with multiple values of CiO . We find that the free energy of hydration depends (12)

linearly on CiO for all ions (shown in Figure S1. for Mg2+ ), within a narrow interval near the target experimental value. Thus, the data points are fitted using a linear function and the optimized (12)

values of the repulsive coefficients CiO are those that yield ∆Gsim = ∆Gexp . The optimized parameters for the interaction between ions and the water oxygens are shown in sim ) calculated from simulations using Table 1. The standard intrinsic solvation free energies (∆Gslv

these parameters are given in Table 2. The contribution of each correction term for the calculated exp

solvation free energy, as well as the reference experimental solvation free energies (∆Gslv ) are also shown in the same table. In all cases, our optimized parameters yield free energies of hydration that agree with the experimental target value, within the given uncertainty. To assess the extent to which ion-water parameters are transferable between water models, we compare our ion-water interaction parameters with parameter sets which were optimized for other water models using a methodology close to ours. This is the case for the ion-water interaction


Page 13 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Optimized ion-Ow LJ parameters. Ion Li+ Na+ K+ Rb+ Cs+ F− Cl− Br− I− Mg2+ Ca2+

(6)

CiO − 3 [10 kJ·mol−1 ·nm6 ] 0.1107235 0.4376606 1.7719211 2.7449745 4.3637442 2.2344670 5.5707096 7.4031937 10.373364 0.2426809 1.1917846

(12)

CiO − 7 [10 kJ·mol−1 ·nm12 ] 1.0962 5.4775 22.170 41.120 75.510 14.789 100.96 172.827 385.05 1.6165 10.411 (6)

parameters developed by Reif and Hünenberger 49 for monoatomic ions in SPC water, where CiO (12)

parameters were also obtained using Equations 2 and 3, and CiO was parameterized to reproduce (6)

free energies of solvation. Our CiO parameters are almost identical to those of Reif and Hünen(6)

berger, as expected, because the self-interaction parameter COO between two water oxygens differs only slightly between the two water models. However, the present parameter set and that by Reif (12)

and Hünenberger differ strongly in the CiO parameters for many of the ions, as shown in Table 3. (12)

(12)

For example, the parameter obtained for Li+ is substantially larger here: CLiO = 1.80CLiO,RH , where the subscript RH indicates the parameter from Reif and Hünenberger. Only for the K + , (12)

Rb+ and I − ions are the CiO parameters obtained by Reif and Hünenberger almost identical to ours. Interaction parameters for Mg2+ or Ca2+ and SPC/E water have been optimized by Netz et (6)

(12)

al., 28 also based on free energies of hydration. For those interactions, both CiO and CiO differ substantially from our parameters, as shown in Table 4. (12)

Thus, the parameters, and especially CiO , are not in general easily transferable between the SPC and TIP5P water models. To further understand the origin of this lack of transferability, we look at the differences between the two water models. In TIP5P water, the lone pairs of the oxygen atom are explicitly represented by two dummy atoms carrying negative charge placed at a distance 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

sim ) using Table 2: Standard intrinsic solvation free energies of ions calculated from simulation (∆Gslv exp optimized ion-water LJ parameters, and from experiment 45,64 (∆Gslv ). The various correction terms used to calculate the standard intrinsic solvation free energy from simulation are also shown (details are presented in SI).

Ion Li+ Na+ K+ Rb+ Cs+ F− Cl− Br− I− Mg2+ Ca2+

LS raw ∆Gchg ∆GC1 ∆GD h Li − 1 − 1 nm kJ·mol kJ·mol kJ·mol−1 4.988 -526.2 ±0.2 -19.721 0.162 4.988 -425.5 ±0.2 -19.718 0.142 4.988 -355.9 ±0.2 -19.714 0.127 4.989 -329.2 ±0.2 -19.711 0.118 4.989 -305.6 ±0.2 -19.709 0.114 4.988 -461.9 ±0.2 19.717 0.137 4.989 -341.6 ±0.2 19.707 0.111 4.989 -316.7 ±0.2 19.705 0.107 4.990 -282.7 ±0.1 19.698 0.098 4.987 -1907.2 ±0.5 -39.444 0.679 4.987 -1590.4 ±0.4 -39.436 0.569

∆Gcav ∆Gstd − 1 kJ·mol kJ·mol−1 9.45±0.07 7.95 13.14±0.10 7.95 13.19±0.07 7.95 13.62±0.08 7.95 12.45±0.14 7.95 4.77±0.07 7.95 10.80±0.14 7.95 12.08±0.14 7.95 17.69±0.12 7.95 9.06±0.12 7.95 10.61±0.08 7.95

exp

sim ∆Gslv ∆Gslv − 1 kJ·mol kJ·mol−1 -528.4±0.3 -529.4 ±8 -424.0±0.3 -423.8±8 -354.4±0.3 -352.0±8 -327.2±0.2 -329.3±8 -304.8±0.3 -306.1±8 -429.3±0.2 -428.8±8 -303.0±0.3 -304.2±8 -276.9±0.3 -277.4±8 -237.3±0.3 -240.0±8 -1928.9±0.6 -1931.4±8 -1610.7±0.5 -1608.3±8

(12)

Table 3: Factor, f , by which the CiO determined in this work differ from those determined by (12)

(12)

Reif and Hünenberger 49 for interactions with SPC water (CiO = f × CiO,RH ). Ion f

Li+ Na+ 1.80 1.25

K+ 1.0

Rb+ 1.0

Cs+ 0.90

F− 0.575

Cl− 0.80

Br− 0.865

I− 1.0

Table 4: Factor, f , by which the parameters determined in this work differ from those determined (n) (n) by Fyta and Netz for interactions with SPC/E water (CiO = f × CiO,FN , where FN indicates the parameters from reference 28 ). ion f C12 iO f C6

iO

Mg2+ 1.78 0.51

Ca2+ 1.50 0.81

of 0.7 Å from the oxygen. This leads to a greater separation between the sites carrying charges in the TIP5P model. Therefore, while both water models have approximately the same dipole moment, the partial charges in the TIP5P model are almost by a factor two smaller, than those in SPC. In addition, the Lennard-Jones excluded volume interaction is carried by the water oxygen, so that the charges carried by the lone pairs are placed further to the outside of the water molecule. Looking at the scaling factors in Table 3, we observe that for cations f >1, i.e. a stronger repulsive interaction is required for TIP5P than for SPC, whereas for anions f 1 factors introduce more repulsion to compensate for this effect. The result of the stronger repulsion in the optimized parameters is a small shift of the first peak in the cation-lone pair radial distribution functions (shown in Figure S3 for Li+ and Na+ ) to larger distances with the optimized parameters. For anions, on the other hand, the lower charge on the hydrogen atoms, +0.241e in TIP5P compared to +0.41e in SPC means a weaker electrostatic attraction in TIP5P, which has to be compensated by the weaker repulsion indicated by f

Is the Solution Activity Derivative Sufficient to ... - ACS Publications

Recommend Documents