Single-ion thermodynamics from first principles. Calculation of the

Subscriber access provided by Kaohsiung Medical University

Thermodynamics

Single-ion thermodynamics from first principles. Calculation of the Absolute Hydration Free Energy and Single-Electrode Potential of Aqueous Li using ab initio Quantum Mechanical/ Molecular Mechanical Molecular Dynamics Simulations. +

Niko Prasetyo, Philippe Henry Hünenberger, and Thomas S. Hofer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00729 • Publication Date (Web): 04 Oct 2018 Downloaded from http://pubs.acs.org on October 6, 2018

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 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 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.

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 55 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

Submitted to Journal of Chemical Theory and Computation - Revision 1

Single-Ion Thermodynamics from First Principles: Calculation of the Absolute Hydration Free Energy and Single-Electrode Potential of Aqueous Li+ using ab initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations. Niko Prasetyo1,2,3 , Philippe H. H¨ unenberger4 and Thomas Hofer* 1 1

Theoretical Chemistry Division, Institute of General, Inorganic and Theoretical Chemistry, University of Innsbruck, Innrain 80-82, A-6020 Innsbruck, Austria

2

Austria-Indonesia Centre (AIC) for Computational Chemistry, Universitas Gadjah Mada, Sekip Utara, Yogyakarta 55281, Indonesia

3

Department of Chemistry, Faculty of Mathematics and Natural Sciences, Universitas Gadjah Mada, Sekip Utara, Yogyakarta 55281, Indonesia 4

Laboratorium f¨ ur Physikalische Chemie, ETH Z¨ urich, ETH-Hönggerberg, HCI, CH-8093 Z¨ urich, Switzerland

September 21, 2018

*

[email protected]

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 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

Abstract A recently proposed thermodynamic integration (TI) approach formulated in the framework of quantum mechanical/molecular mechanical molecular dynamics (QM/MM MD) simulations is applied to study the structure, dynamics and intrinsic single-ion hydration free energy ∆s G◦M + ,wat of the Li+ ion at a correlated ab initio level of theory. ◦ Based on the results, the absolute intrinsic hydration free energy GH + ,wat of the pro-

ton, the surface electric potential jump χ◦wat upon entering bulk water, and the abso◦ lute redox potential VH + ,wat of the reference hydrogen electrode, are calculated to be

-1099.9±4.2 kJ·mol−1 , 0.13±0.08 V and 4.28±0.04 V, respectively, in excellent agreement with the values recommended by H¨ unenberger and Reif based on an extensive evaluation of the available experimental data (-1100±5 kJ·mol−1 , 0.13±0.10 V and 4.28±0.13 V). The simulation results for Li+ are also compared to those for Na+ and K+ from a previous study in terms of relative hydration free energies ∆∆s G◦M + ,wat and relative elec◦ trode potentials ∆VM The calculated values are found to agree extremely well + ,wat .

with the experimental differences in conventional hydration free energies ∆∆s G•M + ,wat ◦ and standard redox potentials ∆∆H VM + ,wat . The level of agreement between simulation

and experiment, which is quantitative within error bars, underlines the substantial accuracy improvement achieved by applying a highly demanding QM/MM approach at the resolution-of-identity second-order Møller-Plesset perturbation theory (RIMP2) level over calculations relying on purely molecular mechanical or density functional theory (DFT) descriptions. The detailed analysis of the structural and dynamical properties of the Li+ hydrate indicates that a correct description of the solvation structure is achieved as well at this level of theory. The consideration of the QM/MM potential energy components also shows, that the partitioning of the system into QM and MM zones does not induce any significant energetic artefact for the system considered.

1

Introduction Single-ion thermodynamics, i.e. the sub-field of this discipline concerned with the proper-

ties of atoms and molecules carrying a net charge viewed as separate entities, is still an area of active research and debate nowadays. 1–3 The reason for this is related to the electroneutrality constraint on macroscopic matter at equilibrium, 4 a consequence of the strength and


Page 2 of 55

Page 3 of 55 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


long-ranged nature of the Coulombic forces. 5 Despite immense efforts over more than one century to integrate single–ion processes into the thermodynamics framework, a number of challenging questions remain unresolved to date. 1 A prominent example is the decomposition of the hydration free energy ∆s G◦wat of a neutral salt into anionic and cationic contributions. 2 A closely related problem is that of determining ◦ the absolute potential Vwat of an aqueous Galvanic element, i.e the partitioning of the elec-

tromotive force of a cell into absolute anodic and cathodic single–electrode potentials. 6,7 A third one is the evaluation of the electric potential jump χ◦wat felt by a test charge upon transfer from vacuum/air into bulk water. 1 Although these single–ion quantities can be calculated from each other via simple relations involving well-established quantities and ultimately relate ◦ 1 (Fig. 1a), they cannot to the absolute intrinsic hydration fee energy GH + ,wat of the proton

be determined unambiguously from experiment due to the electroneutrality constraint. For this reason these key quantities have been termed elusive. 1 Three main strategies have been employed to overcome these obstacles in practical applications. The simplest approach is to ignore the problem by introducing conventional (relative) scales involving a standard reference. One example is the conventional scale of single-ion hydration ◦ free energies ∆s G•wat , defined by arbitrarily setting the hydration free energy GH + ,wat of the

proton to zero. 8–12 Alternative scales involve setting the formation parameters of the hydrated proton either to zero 12–15 or to its gas-phase value. 12,13 Another example is the electrochem◦ , defined by arbitrarily setting the single–electrode ical series of redox potentials ∆H Vwat ◦ 6,7 The resulting conventional potential VH + ,wat of the reference hydrogen electrode to zero.

data is entirely predictive in terms of electroneutral processes (e.g. hydration free energy of a salt, hydration free energy differences between two isovalent cations or two isovalent anions, electromotive force of a Galvanic cell), but does not enable to quantify directly the magnitude of the hydration forces for a single ion or the free–energy associated with a particular half-cell reaction. The second strategy to overcome the electroneutrality constraint is to introduce an extrathermodynamic assumption (ETA), i.e a hypothesis aiming at partitioning the anionic/anode and cationic/cathode contributions in a specific experimental situation using arguments based



on physical reasoning. 1 For example, the tetraphenylarsonium–tetraphenylborate (TATB) ETA assumes that the hydration free energies of these two ions are equal due to their similar sizes, geometric properties and relative hydrophobicities. 16 Another popular ETA is the cluster-pair assumption (CPA), which relies on extrapolating the properties of ion-water clusters upon increasing the number of solvent molecules, and assuming that the solvation of the cluster in the bulk becomes independent on the charge sign in the limit of very large clusters. 14,15,17 However, because the involved assumptions cannot be formally proved based on purely experimental data, the relative validities of specific ETAs is still extensively debated ◦ in the literature. 18–26 The hydration free energy GH + ,wat of the proton is commonly used as a

reference quantity to compare the results of using different ETAs. 14,15,17,25–30 Fig. 1b shows ◦ 98 estimates for GH + ,wat derived using various experiments and ETAs as a function of the

year of publication (reproduced from Tab. 5.19 of Ref. 1 ). These values cover a range of about 250 kJ·mol−1 and did not convergence to a consensus with time. An irksome consequence of this spread and lack of convergence is that a wide spectrum of theoretical “predictions” for ◦ GH + ,wat can be “validated” against some “experimental” value. For example in Fig. 1b (see

the bar at the right of the graph), every “prediction” located in the grey area is within 1% error of an “experimental” estimate. The third strategy to overcome the electroneutrality constraint is to rely on theoretical models, 31–36 where the consideration of an individual charged species in a bulk environment is perfectly feasible. 37–46 In this context the Born solvation model, 47 formulated in the framework of continuum electrostatics (CE), 48 represents the simplest and best established basis for intuitive reasoning and further refinements. This model describes the solvation of a nonpolarizable spherical ion in a homogeneous dielectric environment, and leads to hydration free-energy estimates that are intrinsic, i.e. omitting the contribution resulting from crossing the surface of the solvent. Furthermore, it disregards the granularity of the bulk liquid, i.e. the potential in the bulk of the non-polarized solvent (no ion or zero ionic charge) is zero. This setting implicitly represents an ETA, which was referred to previously 3 as the Born ETA, and effectively underlies most of the ETAs employed to date in the experimental literature. 1,3 However, due to an oversimplified representation of short-ranged ion-solvent interactions and


Page 4 of 55

Page 5 of 55 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


their dependence on ambiguous (experimentally elusive) ion-size parameters, hydration free energies calculated using the Born model are subject to very large inaccuracies. 1 Fortunately, modern computational approaches nowadays enable the application of models that go beyond the simple description of ion-solvent interactions afforded by CE. 1,2,49–51 The three main challenges in the application of these advanced models to the calculation of single-ion hydration free energies are: 3 (i) preserving the compatibility with the Born ETA; (ii) improving the accuracy in the description of short-range ion-solvent interactions; (iii) removing the arbitrariness in the definition of the ion-size parameters. One of the most general and reliable methods to calculate free-energy differences between two states of a molecular system is thermodynamic integration 52 (TI) implemented in the framework of molecular dynamics 31–34 (MD) simulations. In the TI approach the Hamiltonian of the hydrated ionic system is changed from that of a non-interacting mass point in the pure liquid (state A, λ=0) to that of the fully interacting hydrated ion (state B, λ=1). Based on a set of independent MD TI simulations at fixed λ values between 0 and 1, the hydration free energy is calculated by numerical integration over the ensemble average of the Hamiltonian λderivative. The advantage of this approach over simpler cluster calculations based on a limited set of low–energy configurations 23,26,53–57 is that it relies on correct Boltzmann sampling at a clearly defined state-point (e.g. 298.15 K and 1 bar) and properly accounts for bulk hydration by including a large number of explicit water molecules under periodic boundary conditions (Fig. 2a), possibly complemented by corrections for finite-size effects. 1,58–61 The accuracy of the forces acting on the atoms is a key factor for the accuracy of an MDbased TI calculation. In contrast to the common perception that a monoatomic ion in water is a “simple” chemical system, the ion-water interactions are already fairly complex. 38–40,62 Classical molecular mechanical (MM) models 35,36 clearly represent a large improvement over the CE description afforded by the Born model. However, they depend on arbitrary parameters for defining the ion size and typically fail to replicate the results of more demanding simulation techniques based on quantum mechanics (QM), 3,63,64 even when these ion-size parameters are optimized. 3,63–65 In the context of ionic solvation, the neglect of charge transfer 50,66–72 represents a particularly serious shortcoming of MM models. A number of studies



indeed suggest that effective charges on the order of +0.7 to +0.9e might be more appropriate to account for the QM electron density distribution near an aqueous monovalent cation, invoking either charge-transfer 3,50,66–73 or polarization 74–76 arguments. In the former case, however, it is not entirely obvious how to distribute the rest of the charge to the nearby water molecules in the context of a MM simulation protocol. Still, fixing the charge of the ion in the MM model to its nominal charge of +1e amounts to entirely neglecting these effects. Another key shortcoming of the MM-based description is its dependence on an elusive ion-size definition, 3,77 i.e. on a specific choice of ion-solvent Lennard-Jones interaction parameters. As a result of this dependence MM models are affected by the same inherent ambiguity as the Born model, because the calculated ∆s G◦wat values strongly depend on a specific choice for these Lennard-Jones parameters. 1 There exist two main approaches for incorporating a QM description of the ion-solvent interaction into an MD simulation while keeping the computational effort tractable. The first one is Car–Parrinello/Born-Oppenheimer molecular dynamics 78,79 (CPMD or BOMD). Here, the system size is reduced to a minimum (typically the ion plus about 30 – 100 water molecules), the QM level is limited to periodic density functional theory (DFT) with a generalised gradient approximation (GGA), 80,81 and the sampling is restricted to short timescales (typically about 10 to 20 ps sampling). Further limitations are an ambiguous definition of the reference potential 82–88 and approximations in the treatment of electron correlation. The latter approximations may largely affect the description of the properties for both bulk water 89–97 and hydrated ions. 98–101 Semi-empirical dispersion corrections 102,103 are also frequently applied, which re-introduce a dependence on arbitrary ion-size parameters. Despite these shortcomings, CPMD/BOMD approaches have been applied to study the hydration free energies of ions. 104 However, since no details about the hydration structure and dynamics have been provided, it remains difficult to assess the accuracy of these calculations based on experimentally accessible properties. The second approach relies on a hybrid QM/MM description. 105–109 Here, the chemically most relevant part of the system, namely the ion and its first hydration shell, is described at the QM level of theory, while MM interaction potentials are considered sufficiently accurate


Page 6 of 55

Page 7 of 55 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


to model the ion-solvent and solvent-solvent interactions beyond the first shell (Fig. 2a). The QM/MM approach is not limited to a particular subset of QM methods. In particular, it is possible to apply higher levels of theory such as resolution-of-identity second-order MøllerPlesset perturbation theory (RIMP2) 110–113 to take electron correlation effects into account. Comparing the hydration properties of Na+ and K+ obtained at the DFT and MP2 levels, the authors of Ref. 43 concluded that the application of the latter greatly improves the accuracy of results. 43 A similar conclusion was reached in Ref., 104 highlighting that the corrections based on RIMP2 are an important requirement to improve the hydration free energies calculated by DFT-BOMD simulations. The convergence of energy and free-energy profiles in QM/MM simulations of proton transfer in DNA was also investigated in Ref., 114 which concludes that the selection of an accurate QM/MM method such as MP2, of a chemically intuitive selection of the QM region, and of a proper simulation protocol are important requirements for accuracy next to the size of the QM region. By focusing the RIMP2 treatment on the critical interactions in a QM/MM approach, we recently showed in the context of hydrated Na+ and K+ 3 that long equilibration and simulation times can be achieved for comparatively large ion-solvent systems, which in turn leads to well-converged free-energy estimates with narrow confidence intervals. Another crucial element for the calculation of accurate single-ion hydration free energies at the QM/MM level is to include appropriate corrections terms for finite-size and cutoff effects, for standard-state conversion, as well as for enforcing compatibility with the Born ETA. These correction terms are described in details elsewhere. 1,3,58–61 They are easy and inexpensive to evaluate in a post-processing step, and do not lead to any overhead when executing the simulations. The application of this combined QM/MM MD/TI simulation strategy in Ref. 3 provided a remarkably accurate description of the hydrated Na+ and K+ ions at the structural, dynamical and energetic level. In particular, the calculated difference ∆∆s G◦M + ,wat in hydration free energies (73.9 kJ·mol−1 ; K+ minus Na+ ) agreed remarkably well with the recommended experimental value 1 (71.7 kJ·mol−1 ). Based on the calculated ∆s G◦M + ,wat for the two ions, ◦ ◦ ◦ estimates for the quantities GH + ,wat , VH + ,wat and χwat were also provided, and found to also



Page 8 of 55

agree very well with the recommended values of Ref. 1 for these elusive quantities. The goal of the present article is to report a third set of calculations performed using the same methodology, now applied to the aqueous Li+ ion. The hydration structure and dynamics of this comparatively smaller ion is compared to that of Na+ and K+ . The consideration of a third ion also represents an important validation step for the proposed QM/MM MD/TI methodology. It provides two new ∆∆s G◦M + ,wat values for comparison with experiment, as ◦ ◦ ◦ well as a new set of estimates for the elusive quantities GH + ,wat , VH + ,wat and χwat .

2

Computational Method The following sections summarise the key aspects of the employed methodology. A more

detailed description can be found in Ref. 3

2.1

Thermodynamic Integration

The basic thermodynamic cycle underlying the calculation of a single-ion hydration free energy via TI is shown in Fig. 2b (the cycle actually employed in the present work is slightly more complicated, see section 2.2). The isolated cation M+ (vac) is transformed into a noninteracting point mass X(vac), transferred into aqueous solution as X(aq), and transformed to the hydrated ion M+ (aq). For the standard hydration free energy ∆s G◦M + ,wat , this results in the equation

∆s G◦M + ,wat = ∆G1 + ∆G3 + ∆G◦std ,

(1)

where ∆G1 and ∆G3 correspond to the X(vac)← M+ (vac) and X(aq)→ M+ (aq) processes, ∆G◦std is a standard-state correction, and the process X(vac)→X(aq) has been omitted as it is associated with a free-energy change ∆G2 = 0. In principle, the sum ∆G1 + ∆G3 could be determined computationally using TI, by introducing a linear combination of the Hamiltonians HA and HB corresponding to the initial and final states using a coupling parameter λ changing from zero to one 52


Page 9 of 55 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


H(λ) = (1 − λ)HA + λHB = VSS + λVIS ,

(2)

where VSS and VIS represent the solvent-solvent and ion-solvent potential-energy terms, the latter being gradually “switched on” upon increasing λ. The potential-energy VII corresponding to the QM energy of the ion has been omitted, as it cancels out between ∆G1 and ∆G3 . In addition, since the masses of all particles remain unmodified, no contribution to the free energy arises from a change in the kinetic energy. By performing independent simulations for a set of suitably spaced λ-values between 0 and 1, the corresponding free energy is then obtained using the TI formula 52 Z1 ∆GAB = GB − GA =

dλ

∂H(λ) ∂λ

0

Z1

dλ hVIS iN P T λ ,

= NP T λ

(3)

0

where h. . .iN P T λ denotes ensemble averaging at a given λ-value under isothermal-isobaric (N P T ) conditions.

2.2

QM/MM TI

The direct determination of ∆G1 + ∆G3 in Fig. 2 and Eq. 1 via QM/MM MD/TI is not practical. First, it would require two QM calculations per MD step, one in the presence of the ion and one in its absence. Second, for simulations close to the initial state, the ionsolvent interaction energy and forces resulting from HB are strongly reduced, so that ion and solvent may sometimes overlap in space. As a result, the free energy derivative ∂H/∂λ in Eq. 3 presents large fluctuations and the corresponding ensemble average fails to converge even for very long simulations (end-point catastrophe). While this problem can be alleviated in MM TI simulations by the introduction of a λ-dependent soft-core potential, 115–119 there exists no equivalent formulation for QM methods, although grand-canonical approaches 120–124 represent a promising strategy in this context. In the present work, an alternative route for calculating ∆G1 + ∆G3 is thus followed, which involves splitting the path into three successive steps, namely cavitation, charging and quantization 3,125–128 (Fig. 3). The sum of the corresponding free-energy changes increased 9 ACS Paragon Plus Environment


Page 10 of 55

by the contributions ∆Gcor of the correction scheme 58–61 and a standard-state adjustment ∆G◦std (sections 2.4 and 2.5), yield the absolute intrinsic hydration free energy of the cation at QM/MM level

∆s G◦wat = ∆Gcav + ∆Gchg + ∆Gqua + ∆Gcor + ∆G◦std

(4)

The cavitation free energy ∆Gcav corresponds to the reversible work for inserting an uncharged Lennard-Jones particle MLJ into the solvent. It is calculated by TI at MM level using 33 equidistant λ-points. The charging free energy ∆Gchg corresponds to the reversible work for charging the Lennard-Jones particle leading to the fully established ion in solution. It is also calculated by TI at the MM level using 9 equidistant λ-points. For both ∆Gcav and ∆Gchg , the sampling time per λ-point is in the range of 2 to 3 ns after 0.5 ns equilibration. This sampling time is adjusted independently for each λ-point to achieve similar statistical errors in the respective Hamiltonian derivatives over the entire λ-range. For ∆Gchg , a linear Hamiltonian coupling scheme is applied. To avoid the occurrence of an end-point catastrophe at λ = 0, the calculation of ∆Gcav relies on a soft-core Hamiltonian coupling scheme 116 with the corresponding parameter αLJ set to 0.5. Note that uncharacteristically long simulations times are employed here to achieve small error bars. 3 The quantization free energy ∆Gqua corresponds to the reversible work for converting the MM representation of the hydrated ion to a QM/MM representation. The setup and execution of the corresponding TI calculation represents the main methodological and computational challenge of the present approach. At each simulation timestep, the energy of the simulation box is calculated twice, once at the MM level and a second time at a QM/MM level. 125,126 The corresponding energies (and forces) are then combined using a linear Hamiltonian coupling scheme

HM M →QM/M M = (1 − λ)HM M + λHQM/M M .

(5)

To make energies of the MM and QM/MM descriptions compatible, a common zero energy level E0QM has to be selected. This is done by anchoring the interaction energy to zero at both


Page 11 of 55 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


levels for the situation in which the interaction partners (ions, water molecules) are at infinite QM QM separation. In the QM case this means that the energies Eion and Ewat corresponding to the

isolated ion and the energy-minimised water molecule have to be subtracted from the total QM energy Etot of the QM region according to

QM QM QM QM V QM = Etot − E0QM = Etot − Eion + nwat · Ewat ,

(6)

where nwat is the number of water molecules included in the QM treatment. The resulting potential energy V QM accounts for the ion-water and water-water interactions, including contributions resulting from the shift in the electron density due to polarization and charge transfer as well as from changes in the molecular geometry. The definition of this joint zero energy level also cancels out the energy contribution of the isolated ion, i.e. ∆Gqua implicitly encompasses the contribution of ∆G1 in Fig. 2b and Eq. 1, the reason for which this term no QM/M M

longer appears in Eq. 4. Addition of the classical potential-energy contribution Vcla

,

composed of the QM/MM coupling energy V QM ↔M M and the interaction energy V M M of the QM/M M

water molecules within the MM region, yields the total QM/MM potential energy Vtot

QM/M M

Vtot

QM/M M

= V QM + Vcla

= V QM + V QM ↔M M + V M M .

(7)

As already noted in previous studies, 129–131 the convergence of MM-to-QM/MM perturbation calculations can be greatly enhanced by selecting MM parameters leading to a description of the system that is as close as possible to the QM one. Here, to improve the convergence of the ∆Gqua calculation and thus reduce its computational cost to a minimum (see Suppl. Mat. Section I and Fig. S1). The required number of λ-points can then be lowered to a minimum of 5 and the convergence of ∆Gqua is simultaneously largely improved. 3 It is important to stress, however, that the ion-solvent Lennard-Jones interaction parameters, which are physical parameters in a MM calculation, become merely numerical parameters in the present QM/MM calculations, i.e. they influence the convergence properties but no longer the result at full convergence. In this work the term ∆Gqua quantization term is calculated using 9 equidistant λ-points



with sampling times of 0.3 ns after equilibration periods of 0.02 ns.

2.3

QM/MM Simulation Setup

The hybrid QM/MM approach 105–109 partitions the system into two regions (Fig. 2a). The QM region consists of the ion and its first hydration shell, which can encompass a variable number of water molecules nwat along the simulation. The MM region includes a sufficient number of additional water molecules further away from the ion to adequately represent the bulk liquid. In the present simulations 2000 water molecules are included in total (both zones). The coupled simulation of the two zones follows the previously described 3,37–39,132–139 quantum mechanical charge field (QMCF) MD approach. The rigid three-sites extended simple point charge (SPC/E) water model 140 is employed for the representation of the MM water molecules. The QMCF MD approach relies on electrostatic embedding 141–143 where the partial charges of atoms in the MM zone induce an external Coulomb-potential in the QM Hamiltonian during the electronic structure calculation. This additional term has a similar form as the nucleielectron interactions, but with the nuclear charges Z replaced by the MM charges q in the corresponding operator. 141–143 Since the latter is a one-electron operator, this inclusion only incurs a moderate increase in computational cost. Electrostatic embedding has been shown recently 144 to compare favourably with other embedding approaches in terms of convergence with respect to the size of the QM region. The electrostatic influence of the QM subsystem onto the MM region must also be considered carefully, so as to avoid imbalances. To remain compatible with the effective parametrization of the SPC/E model 140 and the reaction-field (RF) scheme 145,146 applied in the MM treatment, comparable estimates for the effect of this influence is accounted for by means of QM-derived atomic partial charges. 3 Despite their known shortcoming of being sensitive to the employed basis set, Mulliken charges are used here in order to satisfy this compatibility requirement. 147 Note that the RF correction applied in the MM region is also applied in the QM zone on the basis of these charges.


Page 12 of 55

Page 13 of 55 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


2.4

Correction scheme

The charging free energy ∆Gchg is highly sensitive to a number of methodological choices entering into its calculation. 1,58–61 These are due to finite-size, approximate-electrostatics and potential-summation errors occurring in the MM simulations. The corresponding correction term ∆Gcor is obtained at MM/CE level, and consists of four terms A-D as formulated by Kastenholz, Reif and H¨ unenberger. 1,58–61,148,149 For the present simulations of a periodic system using RF electrostatics with an atom-based cutoff, 150 these contributions correct for: (A) The deviation of the solvent polarization around the ion relative to the polarization in an ideal Coulombic system and the omitted interactions of the ion with the polarized solvent beyond the cutoff. (B) The deviation of the solvent polarization around the ion in the finite periodic system relative to the polarization in a macroscopic system. (C) The deviation of the solvent-generated electric potential at the ion site as calculated from the simulated trajectory relative to the correct electric potential, a consequence of the application of an atom-based summation scheme (P-summation) instead of a molecule-based one (M-summation) for this potential. (D) An inaccurate dielectric permittivity of the employed solvent model. The correction scheme is formulated to produce charging free energies which are: (i) consistent across a wide range of simulation protocols; 59,61 (ii) intrinsic, 1,2,59 i.e. accounting for bulk solvation effects only and excluding the work of crossing the water surface; (iii) compatible with the Born ETA. The evaluation of these corrections terms in the present calculations is detailed in Suppl. Mat. Section S.II. The involved equations are analytical, relying on parameters that were fitted against the results of more complicated continuumelectrostatics calculations, 60 and their evaluation is simple and computationally inexpensive.



2.5

Page 14 of 55

Standard state conversion

Since in the bbmeT -convention 1 the hydration free energies are tabulated at a reference pressure P ◦ of 1 bar and a reference molality b◦ of 1 mol·kg−1 for the gas-phase and ideal solute reference states, respectively, a constant contribution

∆G◦std

b◦ ρwat = RT ln RT P◦

(8)

has to be added, where R is the ideal gas constant, T the reference temperature (298.15 K) and ρwat the density of water at P ◦ and T . The resulting value is 7.95 kJ·mol−1 for water at room temperature.

2.6

Confidence intervals

Error estimates σλ on the average free-energy derivative h∂H/∂λi of Eq. 6 were obtained by block averaging. 31,151 The total error estimate σ∆G on a calculated ∆G is then evaluated using numerical integration as

σ∆G

1/2  1 Z =  dλ σλ2  .

(9)

0

This error estimation is performed separately for ∆Gcav , ∆Gchg and ∆Gqua , the corresponding error bar on ∆s G◦wat being obtained via error propagation. 152,153 To increase the confidence interval to 99%, the result is multiplied by a t-factor of 3.0 appropriate for large data sets. 152,153

2.7

Simulation Parameters

The QM/MM MD/TI calculations are performed at the RIMP2 level 111–113 using Def2SVP 154 and DZP 155 basis sets for Li+ and water, respectively. The corresponding auxiliary basis sets are used for all QM atoms. 156 Orbitals below -3.0 a.u. are considered frozen, which in the present case only applies to the 1s orbitals of the oxygen atoms, in-line with the recommendation of Ref. 157 The radius of the QM region is set to 0.31 nm, so as to include


Page 15 of 55 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


the ion and its entire first solvation shell. To ensure a smooth transition of water molecules between the QM and MM regions, a smoothing zone is introduced from 0.29 to 0.31 nm where the atomic forces are progressively switched between their QM and MM values. 37–39,135,139 All QM calculation were carried out using the TURBOMOLE 158,159 package. The structural visualisations were performed using the VMD program 160 .

2.8

Calculation of elusive quantities

After determining the standard intrinsic hydration free energy ∆s G◦M + ,wat of an ion M + via ◦ QM/MM MD/TI simulations, the intrinsic solvation free energy GH + ,wat of the proton can be

calculated by comparison to the (non-elusive 1 ) conventional solvation free energy ∆s G•M + ,wat of the ion: 1

◦ ◦ • GH + ,wat = ∆s GM + ,wat − ∆s GM + ,wat .

(10)

◦ The absolute half-cell potential VH + ,wat of the reference hydrogen electrode and the air-liquid

interfacial potential χ◦wat can then be calculated as

◦ VH + ,wat =

and χ◦wat =

◦ ◦ GH + ,wat + ∆f GH + ,vac

(11)

F ◦ ˆ◦ GH + ,wat − GH + ,wat

F

,

(12)

◦ where ∆f G◦H + ,vac is the formation free energy of the gas-phase proton and GˆH + ,wat its stan-

dard real solvation free energy (both non-elusive 1 ), and F denotes the Faraday constant ◦ (96485.3399 C·mol−1 ). One may also calculate the absolute single-electrode potential VM + ,wat

for the ion as

◦ VM + ,wat

=

∆s G◦M + ,wat + ∆f G◦M + ,vac F

,

(13)

where ∆f G◦M + ,vac is the gas-phase formation free energy of the ion. The difference in ∆s G◦M + ,wat between two ions can be compared to the experimental differ15 ACS Paragon Plus Environment


ence ∆∆s G•M + ,wat between the corresponding conventional hydration free energies ∆s G•M + ,wat . ◦ Alternatively, the difference in VM + ,wat between two ions can be compared to the difference ◦ ◦ ∆∆H VM + ,wat between the corresponding standard redox potentials ∆H VM + ,wat . The values

for all non-elusive quantities required in Eqs. 10 to 13 or used in the indicated comparisons are taken from Ref. 1 and summarized in Tab. 1.

3

Result and Discussion

3.1

Structure and Dynamics

This section discusses the structural and dynamical properties of hydrated Li+ obtained from the MM and QM/MM MD simulations at full ionic charge. Although the QM/MM representation essentially provides a first-principle description of the system, the MM model relies on Lennard-Jones interaction parameters controlling the ion size and its dispersion interactions. These have been pre-optimized against the QM properties as described in Suppl. Mat. S.I. Validation of the QM/MM description in terms of structure and dynamics is essential because the hydration free energy is a single number which may coincidentally agree with experiment even if the underlying physical description is inaccurate. The comparison of key properties between MM and QM/MM models as well as with experimental data is presented in Tab. 2. The characteristic features of the ion-oxygen distance (peak positions and coordination numbers) are obtained from the associated radial distribution functions (RDFs) depicted in Fig. 4a. At both the MM and QM/MM levels, two well-defined hydration shells are identified. The first–shell peaks are located at nearly identical distances r1 of 0.197 and 0.196 nm, respectively, in excellent agreement with the experimental estimates of 0.194– 0.198 nm derived from scattering measurements 161–164 as well as with previous QM estimates of 0.194–0.196 nm. 165–167 In contrast, the second–shell peak is found at a larger distance of 0.424 nm in the MM simulation compared to 0.408 nm for its QM/MM counterpart. The latter is also closer to the former theoretical estimates. 166,167 Besides this, the main difference between the two representations is in the shape of the first peak, which is significantly broader


Page 16 of 55

Page 17 of 55 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


and flatter in the QM/MM case. The corresponding coordination number distributions (CNDs) are displayed in Figs. 4b and 4c. Only small differences are observed for the first shell, with a slight trend towards lower coordination numbers at the QM/MM level. This trend is more pronounced for the second shell CNDs, with a shift of the CND by about one unit. The resulting average first- and second-shell coordination numbers are 4.05 and 13.7 in the QM/MM simulation, compared to 4.11 and 14.5 in the MM case. A value close to four for the first-shell is expected for a tetrahedral arrangement of the water molecules, and supported by the experimental results. 161–164 Analysis of the water arrangement within the first shell via the oxygen-oxygen three-body distribution 168 is shown in Fig. 5a. The function is very similar at the MM and QM/MM levels, and reveals unimodal distribution of the oxygen-oxygen distances within the first shell, with a peak distance of about 0.32 nm. Compared to the first-shell peak distance of 0.196 nm of the ion-oxygen RDF, the ratio of 1.63 agrees very well with the value of 2×sin(54.2◦ ) = 1.62 expected for an ideal tetrahedral arrangement. The vibrational power spectrum of the ion-oxygen stretching vibration is shown in Fig. 5b. It points towards a significant difference in the effective ion-oxygen interaction strength at the MM and QM/MM levels, with peaks at 605 and 471 cm−1 , respectively, corresponding to effective force constants of 97.4 and 63.8 N·m−1 . The value obtained in the QM/MM case agrees well with the experimental estimates of 500–520 cm−1 obtained for lithium-oxygen complexes in silicon matrices. 169,170 The force–constant reduction observed when changing from the MM to a QM/MM description is in-line with the decrease in the peak height observed for the ion-oxygen RDF. These differences also correlate with a reduction of the mean residence time (MRT) τ of first shell water molecules, evaluated by monitoring exchange events with an allowed excursion time 171 t∗ of 0.5 ps. The MRT decreases from 27.2 ps in the MM simulation to 19.0 ps in the QM/MM case. The associated rate coefficient, 171 R measuring the average number of attempts to achieve one successful solvent exchange (i.e. lasting longer than t∗ ), also simultaneously decreases from 4.08 to 2.84.



Page 18 of 55

Finally, to highlight the importance of ion-water charge transfer, the time series and the distribution of the ionic Mulliken charge along the QM/MM simulation are shown in Figs. 5c and 5d, respectively. The most probable charge is +0.81 e, in good agreement with the value of +0.83 e obtained in an earlier force–matching study. 73 Since the QMCF scheme 37–39,135,139 is designed to incorporate this effect in the QM/MM coupling, the charge reduction causes a decrease in the second–shell coordination number compared to the MM simulation, even though this part of the system is included in the MM region. Most of the trends observed in Tab. 2, Fig. 4 and Fig. 5 upon changing the description from MM to QM/MM are qualitatively the same as those observed from Na+ and K+ in our previous study, 3 namely: (i) broader first shell RDF peak; (ii) slightly more distant second shell; (iii) reduced average coordination number; (iv) decreased first shell residence time. For Na+ and K+ , the change also led to noticeable differences in the oxygen-oxygen threebody distribution function, indicative of a more labile first-shell geometry. For Li+ , the simple tetrahedral arrangement is probably too stable for exhibiting such a dependence. In the opposite, the significant reduction of the ion-oxygen vibrational frequency (interaction strength) upon application of the QM/MM description appears to be specific to Li+ and is probably also related to the tetrahedral arrangement of the first shell water molecules. No corresponding change was observed 3 for Na+ and K+ , presumably because model-dependent interaction changes can be accommodated by adjustments in the first-shell geometry. For the three ions considered, all the observed trends can be rationalised by the neglect of charge-transfer effects in the MM model, i.e. the fact that the model relies on a nominal charge of +1 e that is significantly higher than the effective charge observed in the QM/MM simulations (about 0.6-0.8 e for the ions considered). A further critical aspect when performing QM/MM MD simulations is to ensure that the potential energy of the system is not affected by artefacts related to the presence of the QM/MM interface. The analysis of potential-energy terms carried out in our previous work 3 QM/M M

on Na+ and K+ showed that the total QM/MM potential energy Vtot

in Eq. 7 does not

show any spurious correlation with the number of water molecules nwat included in the QM region. This provides strong support for the adequacy of the applied QMCF MD simulation


Page 19 of 55 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


protocol. 37–39,135,139 A similar analysis carried out for the present simulation of aqueous Li+ is included in Suppl. Mat. section S.III, and leads to the same conclusion.

3.2

Absolute Intrinsic Hydration Free Energy

The curves of the average Hamiltonian λ-derivative as a function of the coupling parameter λ for the TI calculations of the cavitation, charging and quantization free energies are shown in Fig. 6. The resulting free-energy changes after integration are reported in Tab. 3 together with associated error bars corresponding to a confidence interval of 99%. The values of the terms ∆Gcor and ∆G◦std as well as the results for the absolute hydration free energy ∆s G◦M + ,wat at the MM and QM/MM levels are also reported. The corresponding values for aqueous Na+ and K+ are included for comparison. In order to provide a conservative error estimate, the error bars on the individual steps have simply been added. An alternative error estimate obtained by propagation of the variance, i.e. as the square-root of the sum of the variances, is provided in Suppl. Mat. Tab. S.3. The free-energy derivative for the cavitation step is shown in Fig. 6a. Due to the application of a soft-core coupling scheme, 116 the dependence on λ is non-monotonic, with a peak at λ = 0.34. The associated free-energy contribution ∆Gcav amounts to 6.0±0.6 kJ·mol−1 , highlighting the endergonic nature of the cavitation process. The corresponding charging process is shown in Fig. 6b, and presents a dependence of the Hamiltonian derivative on λ that is much closer to linearity, as would be expected for a solvent behaving as a perfect dielectric medium. The limited deviation relative to a straight line is highlighted in the figure. The resulting charging free-energy contribution ∆Gchg is −407.7±0.5 kJ·mol−1 , reflecting the substantial exergonic electrostatic contribution dominating the ion-solvent interaction. To correct the latter contribution for finite-size, approximate-electrostatics and potentialsummation errors in a way compatible with the Born ETA, a correction term ∆Gcor of -143.6±1.0 kJ·mol−1 has to be added. The detailed calculation of this term is described in Suppl. Mat. section S.II. Finally, the standard–state conversion term ∆G◦std of 7.95 kJ·mol−1 has to be added. 172 Summation of the contributions ∆Gcav , ∆Gchg , ∆Gcor and ∆G◦std yields 19 ACS Paragon Plus Environment


the absolute intrinsic hydration free energy at MM level, amounting to -537.3±2.1 kJ·mol−1 . The free-energy derivative curve for the quantization process is shown in Fig. 6c. It also shows an essentially linear dependence, resulting in a free-energy change ∆Gqua of 12.3±0.5 kJ·mol−1 . Thus, the MM simulation somewhat overestimate the strength of the ion-solvent interactions, in-line with the observations made in the previous section regarding the differences in structural and dynamical properties when comparing the MM and QM/MM descriptions. The estimate for the absolute intrinsic hydration free energy ∆s G◦M + ,wat for aqueous Li+ at the QM/MM level is -525.0±2.6 kJ·mol−1 , in remarkable agreement with the value recommended by H¨ unenberger and Reif, 1 namely -525.1±6.7 kJ·mol−1 . A similar level of agreement was previously 3 observed for Na+ and K+ , with calculated intrinsic hydration free energies of -417.3±3.1 and -343.4±2.3 kJ·mol−1 , respectively, also agreeing well with the recommended values 1 of -419.5±6.9 and -347.8±5.9 kJ·mol−1 . However, comparison with non-elusive quantities provides even stronger support for the calculated values. This is done here against the differences between the conventional hydration free energies ∆s G•M + ,wat of the three ions and between their relative redox potentials ◦ ∆H VM + ,wat . Both differences can be determined unambiguously and accurately via experi-

mental measurements. The comparison between the calculated values at the MM or QM/MM ◦ levels and the experimental ∆∆s G•wat and ∆∆H Vwat values for the changes Li+ → Na+ ,

Na+ → K+ and Li+ → K+ is performed in Tab. 4. In all cases the QM/MM estimates agree quantitatively with the experimental data, i.e. within the corresponding cumulated error bars, which is not always the case for the MM estimates. It should be stressed that the relatively small magnitude of ∆Gqua and thus, the acceptable accuracy of the MM estimates, largely results from the fitting of the Lennard-Jones interaction parameters of the MM model against structural properties from previous QM-based simulations. If desired, further fitting could easily achieve a situation in which ∆Gqua is zero. A small ∆Gqua does not imply that the MM model is accurate but only that its parameters have been well calibrated. As seen from the detailed analysis of structure and dynamics of the three hydrated ions considered, the MM model is not physically very accurate, because it is pair-wise additive and neglects charge-transfer effects.


Page 20 of 55

Page 21 of 55 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.3

Calculation of elusive quantities

◦ ◦ ◦ When deriving values for the elusive quantities GH + ,wat , VH + ,wat and χwat based on calcu-

lated intrinsic single ion hydration free energies ∆s G◦M + ,wat , each investigated system leads to an independent set of estimates. The values for the elusive quantities calculated based on the ∆s G◦M + ,wat from the QM/MM MD/TI simulations are reported in Tab. 5. The results ◦ are highly consistent for the three ions. The absolute intrinsic hydration free energy GH + ,wat

of the proton is in the range -1099.9 to -1095.6 kJ·mol−1 in excellent agreement with the recommended value of 1100±5 kJ·mol−1 from H¨ unenberger and Reif. 1 The same applies to ◦ the absolute VH + ,wat potential of the reference hydrogen electrode, with values in the range

4.28 to 4.33 V, also in excellent agreement with the estimate of 4.28±0.13 V from H¨ unenberger and Reif. 1 Similarly, the air-to-liquid surface electric potential χ◦wat is in the range 0.09 to 0.13 V, again in agreement with the recommended value of 0.13±0.10 V from Ref. 1 Since χ◦wat is obtained as the difference between two large quantities (Eq. 12), the associated error bars are comparably large, i.e. approximately 60 to 80%. However, considering the range of experimental estimates 1 for χ◦wat from -1.15 to +1.4 V, the calculated values can still be regarded as highly consistent. In all cases, the agreement between the theoretical and experimentally-inferred values is quantitative, i.e. within the corresponding cumulated error bars. The hydration structures of the three different ions show large variation. For example, the average first-shell coordination numbers for Li+ , Na+ and K+ are 4.1, 5.5 and 7.0, respectively, with average ion-oxygen distances of 0.196, 0.235 and 0.277 nm, and associated hydration free energies of -525.0, -417.3 and -343.4 kJ·mol−1 . In spite of these large differences, the values for the elusive quantities calculated for the three ions are very similar, which provides a strong validation argument for the overall consistency of the present QM/MM MD/TI methodology.

4

Conclusion The application of the QM/MM MD/TI methodology to aqueous Li+ yields very consis-

tent results in excellent agreement with those of our previous investigation 3 concerning Na+



Page 22 of 55

and K+ . The structural and dynamical properties of the hydrated Li+ ion observed in the simulation are also found to agree very well with available experimental and theoretical data, with an average coordination number of 4 (vs. 5.5 and 7.0 for Na+ and K+ , respectively) and a first-shell ion-oxygen distance of 0.196 nm (vs. 0.235 and 0.277 nm for Na+ and K+ , respectively). As observed previously 3 for Na+ and K+ the QM/MM description represents a significant improvement in terms of hydration structure and dynamics, even when the MM Lennard-Jones interaction parameters have been carefully pre-optimized. The main issue in the MM representation is the total neglect of charge transfer effects. By using an ionic charge equal to the its nominal charge of +1 e, the MM model leads to an overstructured RDF first-shell peak, higher average coordination numbers, a longer mean ligand residence time in the first shell, the occurrence of a more pronounced second-shell RDF peak at a shorter distance and a too rigid/static nature of the coordination complex. Interestingly, three recent approaches for improving the description of ion-water interactions at the MM level rely in effect on a weakening of the ion-solvent Coulombic interactions. These are (i) inductiontype 68–71,100,101,173 models, in which charge-transfer is included explicitly; (ii) the molecular dynamics in electric continuum (MDEC) approach, 74–76 using a scaled effective ion charge 174 (by a factor of about 0.7-0.9 in-line with the results of the QM/MM MD TI simulations) and (iii) multi-site ions, 175–181 redistributing the ionic charge over virtual off-center sites within the ion. The intrinsic single-ion hydration free energy (-525.0±2.6 kJ·mol−1 ) as well as the derived ◦ values for the elusive single-ion reference quantities GH + ,wat ,

◦ VH and χ◦wat + ,wat

(−1099.9±4.2 kJ·mol−1 , 4.28±0.04 V and 0.13±0.08 V, respectively) are again in remarkable agreement with the recommended values of H¨ unenberger and Reif 1 (-525.1 kJ·mol−1 , −1100.0±5.0 kJ·mol−1 , 4.28±0.12 V and 0.13±0.10V, respectively). The values of these three quantities have been debated for more than a century and it is likely that the debate will continue. This is because these numbers are assumption-bound, the present estimates referring to what is called here the Born ETA. In contrast, the differences in hydration free energies and redox potentials along the alkali series are not elusive, i.e they are unambiguous and accurately known experimentally. The fact that the QM/MM MD/TI calculations


Page 23 of 55 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


reproduce these differences very well (at least now for Li+ , Na+ and K+ ) certainly represents the strongest validation of the underlying physical model. Here also, model-based predictions have been available for more than a century. However, the present calculations provide the first estimates that are not bound to elusive ion-size parameters, but are entirely determined by the QM physics of the ion-solvent interactions. ◦ + to Na+ and then to A slight systematic drift in the calculated GH + ,wat values from Li

K+ can nevertheless be identified (Tab 5). Although the three estimates agree within their error bars, it is unlikely that the observed trend is entirely due to statistical errors. The main residual error source in the present calculations is probably the application of rigidbody dynamics. Because water molecules exchange between the first and second hydration shells along the simulations, the same constraints to the molecular geometry have to be applied to the solvent in the QM and MM zones (here, full rigidity in the SPC/E geometry). On the one hand, the choice of a fully rigid water model is justified by the fact that the intramolecular vibrations are not excited at room temperature 182,183 (quantum-mechanical ground state), in which case a constraint (no heat-capacity contribution) represents a better approximation than a harmonic oscillator (heat capacity contribution kB ). On the other hand, the choice of a specific geometry for the rigid model that is appropriate for the bulk liquid neglects the effect of possible geometric distortions in the neighbourhood of the ion. Such distortions are expected to play a more significant role when the ions is small 63,104,184,185 or/and highly charged 64,186 (stronger field), and when it is negatively charged 133,136–138,187 (hydrogen-bonding). Because the Li+ , Na+ and K+ ions considered here are monovalent and positively charged, first-shell geometric polarisation effects are anticipated to be small. The residual trend along the series nevertheless correlates with the ion size, i.e also with the expected magnitude of the geometric distortion likely, to be more pronounced for Li+ and least pronounced for K+ . We are currently investigating this issue by adjusting the QM/MM MD/TI methodology to the use of a flexible water model. Although flexible water models have already been applied in previous QM/MM studies of ionic solvation, 132–134,136–138,188–191 the high frequency of the vibrational modes requires a reduction of the MD timestep from 2 fs to about 0.2-0.5 fs, which



corresponds to an increase of the computational cost by a factor of 4-10. However, the consideration of molecular flexibility will probably represent an crucial unavoidable requirement for the extension of the present simulation methodology to multivalent or/and anionic solutes.

5

Acknowledgement NP gratefully acknowledges Financial support from a Ph.D grant of the Indonesia Endow-

ment Fund for Education of the Ministry of Finance Republik Indonesia (LPDP Kementrian Keuangan RI; grant S-2233/LPDP.3/2016). Research reported in this publication was jointly supported by the ASEAN-European Academic University Network (ASEA-UNINET), the Austrian Federal Ministry of Science, Research and Economy, and the Austrian Agency for International Cooperation in Education and Research (OeAD-GmbH). The computations presented have been performed using the HPC infrastructure of the University of Innsbruck.

6

Supporting Information The optimization of the Lennard-Jones potential parameters, the calculation of the cor-

rection terms ∆Gcor and its dependence with respect to the ionic radius RI , an analysis of the energy consistency of the QM/MM MD simulations and alternative error estimates are provided in the Supporting Infomation. This information is available free of charge via the Internet at http://pubs.acs.org

References [1] H¨ unenberger, P. H.; Reif, M. M. Single-ion solvation: Experimental and theoretical approaches to elusive thermodynamic quantities, 1st Ed.; Royal Society of Chemistry: London, UK, 2011. [2] Vlcek, L.; Chialvo, A. A. Single-ion hydration thermodynamics from clusters to bulk solutions: Recent insights from molecular modeling. Fluid Phase Equilib. 2016, 407, 58–75.


Page 24 of 55

Page 25 of 55 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] Hofer, T. S.; H¨ unenberger, P. H. Absolute proton hydration free energy, surface potential of water, and redox potential of the hydrogen electrode from first principles: QM/MM MD free-energy simulations of sodium and potassium hydration. J Chem. Phys. 2018, 148, 222814/1–222814/28. [4] Galembeck, F.; Burgo, T. A. L. Chemical Electrostatics: New Ideas on Electrostatic Charging: Mechanisms and Consequences, 1st Ed.; Springer International Publishing: Cham, Switzerland, 2017. [5] Stone, A. J. The Theory of Intermolecular Forces; Oxford University Press: Oxford, UK, 1995. [6] Rieger, P. H. Electrochemistry, 2nd Ed.; Springer Science + Business Media: Dordrecht, 1994. [7] Hamann, C. H.; Hamnett, A. Electrochemistry; Wiley-VCH: Weinheim, 2007. [8] Fajans, K. Löslichkeit und Ionisation vom Standpunkte der Atomstruktur. Naturwiss. 1921, 9, 729–738. [9] Baughan, E. The heat of hydration of the proton. J. Chem. Soc. 1940, 1403–1403. [10] Benjamin, L.; Gold, V. A table of thermodynamic functions of ionic hydration. Trans. Faraday Soc. 1954, 50, 797–799. [11] Halliwell, H. F.; Nyburg, S. C. Enthalpy of hydration of the proton. Trans. Farad. Soc. 1963, 59, 1126–1140. [12] Rosseinsky, D. R. Electrode potentials and hydration energies. Theories and correlations. Chem. Rev. 1965, 65, 467–490. [13] Noyes, R. M. Thermodynamics of ion hydration as a measure of effective dielectric properties of water. J. Am. Chem. Soc. 1962, 84, 513–522. [14] Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle Jr., T. R. The proton’s absolute aqueous enthalpy 25 ACS Paragon Plus Environment


and Gibbs free energy of solvation from cluster-ion solvation data. J. Phys. Chem. A 1998, 102, 7787–7794. [15] Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle Jr., T. R. Correction to “The proton’s absolute aqueous enthalpy and Gibbs free energy of solvation from cluster-ion solvation data.” [J. Phys. Chem. A 102, 7787-7794 (1998)]. J. Phys. Chem. A 1998, 102, 9308–9308. [16] Marcus, Y. The thermodynamics of solvation of ions. Part 4. Application of the tetraphenylarsonium tetraphenylborate (TATB) assumption to the hydration of ions and to properties of hydrated ions. J. Chem. Soc. Farad. Trans. I 1986, 83, 2985–2992. [17] Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous solvation free energies of ions and ion-water clusters based on an accurate value for the absolute aqueous solvation free energy of the proton. J. Phys. Chem. B 2006, 110, 16066–16081. [18] Coetzee, J. F.; Sharpe, W. R. Solute-solvent interactions. VI. Specific interactions of tetraphenylarsonium, tetraphenylphosphonium, and tetraphenylborate ions with water and other solvents. J. Phys. Chem. 1971, 75, 3141–3146. [19] Kim, J. I. Preferential solvation of single ions. A critical study of the Ph4 As-Ph4 B assumption for single ion thermodynamics in amphiprotic and dipolar-aprotic solvents. J. Phys. Chem. 1978, 82, 191–199. [20] Schurhammer, R.; Wipff, G. About the TATB hypothesis: solvation of the AsΦ+ 4 and BΦ− 4 ions and their tetrahedral and spherical analogues in aqueous/nonaqueous solvents and at water-chloroform interfaces. New. J. Chem. 1999, 23, 381–391. − [21] Schurhammer, R.; Wipff, G. Are the hydrophobic AsΦ+ 4 and BΦ4 ions equally solvated?

A theoretical investigation in aqueous and nonaqueous solutions using different charge distributions. J. Phys. Chem. A 2000, 104, 11159–11168. [22] Scheu, R.; Rankin, B. M.; Chen, Y.; Jena, K. C.; Ben-Amotz, D.; Roke, S. Charge


Page 26 of 55

Page 27 of 55 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


asymmetry at aqueous hydrophobic interfaces and hydration shells. Angew. Chem. Int. Ed. 2014, 53, 9560–9563. [23] Asthagiri, D.; Pratt, L. R.; Ashbaugh, H. S. Absolute hydration free energies of ions, ion-water clusters and quasichemical theory. J. Chem. Phys. 2003, 119, 2702–2708. [24] Vlcek, L.; Chialvo, A. A.; Simonson, J. M. Correspondence between cluster-ion and bulk solution thermodynamic properties: On the validity of the cluster-pair-based approximation. J. Phys. Chem. A 2013, 117, 11328–11338. [25] Pollard, T.; Beck, T. The thermodynamics of proton hydration and the electrochemical surface potential of water. J. Chem. Phys. 2014, 141, 18C512/1–18C512/10. [26] Pollard, T. P.; Beck, T. L. Quasichemical analysis of the cluster-pair approximation for the thermodynamics of proton hydration. J. Chem. Phys. 2014, 140, 224507/1– 224507/11. [27] Ishikawa, A.; Nakai, H. Quantum chemical approach for condensed-phase thermochemistry(III): Accurate evaluation of proton hydration energy and standard hydrogen electrode potential. Chem. Phys. Lett. 2016, 650, 159–164. [28] Rossini, E.; Knapp, E. Proton solvation in protic and aprotic solvents. J. Comput. Chem. 2016, 37, 1082–1091. [29] Zhang, H.; Jiang, Y.; Yan, H.; Yin, C.; Tan, T.; van der Spoel, D. Free-energy calculations of ionic hydration consistent with the experimental hydration free energy of the proton. J. Phys. Chem. Lett. 2017, 8, 2705–2712. [30] Bartmess, J. E. Thermodynamics of the electron and the proton. J. Phys. Chem. 1994, 98, 6420–6424. [31] Allen, M.; Tildesley, D. Computer simulation of liquids; Oxford University Press: New York, 1987. [32] Tuckerman, M. E. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: New York, USA, 2010. 27 ACS Paragon Plus Environment


[33] Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic press: San Diego – London, 2002. [34] Sadus, R. J. Molecular Simulation of Fluids:

Theory, Algorithms, and Object-

Orientation; Elsevier: Amsterdam, 1999. [35] Leach, A. R. Molecular Modelling; Prentice-Hall: Harlow, 2001; Vol. 2nd Edition. [36] Jensen, F. Introduction to Computational Chemistry; John Wiley & Sons Ltd.: Chichester, 2006; Vol. 2nd Edition. [37] Rode, B. M.; Hofer, T. S.; Randolf, B. R.; Schwenk, C. F.; Xenides, D.; Vchirawongkwin, V. Ab initio quantum mechanical charge field (QMCF) molecular dynamics: A QM/MM - MD procedure for accurate simulations of ions and complexes. Theor. Chem. Acc. 2006, 115, 77–85. [38] Hofer, T.; Pribil, A.; Randolf, B.; Rode, B. Ab initio quantum mechanical charge field molecular dynamics - A nonparametrized first-principle approach to liquids and solutions. Adv. Quantum Chem. 2010, 59, 213–246. [39] Hofer, T.; Rode, B.; Pribil, A.; Randolf, B. Simulations of liquids and solutions based on quantum mechanical forces. Adv. Inorg. Chem. 2010, 62, 143–175. [40] Weiss, A. K. H.; Hofer, T. S. Exploiting the capabilities of quantum chemical simulations to characterise the hydration of molecular compounds. RSC Adv. 2013, 3, 1606–1635. [41] Bhattacharjee, A.; Weiss, A. K. H.; Artero, V.; Field, M. J.; Hofer, T. S. Electronic structure and hydration of tetramine cobalt hydride complexes. J. Phys. Chem. B 2014, 118, 5551–5561. [42] Hofer, T.; Tirler, A. Structure and dynamics of the uranyl tricarbonate complex in aqueous solution: Insights from quantum mechanical charge field molecular dynamics. J. Phys. Chem. B 2014, 118, 12938–12951.


Page 28 of 55

Page 29 of 55 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


[43] Rowley, C.; Roux, B. The solvation structure of Na+ and K+ in liquid water determined from high level ab initio molecular dynamics simulations. J. Chem. Theory Comput. 2012, 8, 3526–3535. [44] Rempe, S.; Pratt, L. The hydration number of Na+ in liquid water. Fluid Phase Equilib. 2001, 183-184, 121–132. [45] Lamoureux, G.; Roux, B. Absolute hydration free energy scale for alkali and halide ions established from simulations with a polarizable force field. J. Phys. Chem. B 2006, 110, 3308–3322. [46] Whitfield, T. W.; Varma, S.; Harder, E.; Lamoureux, G.; Rempe, S.; Roux, B. Theoretical study of aqueous solvation of K+ comparing ab initio, polarizable, and fixed-charge models. J. Chem. Theory Comput. 2007, 3, 2068–2082. [47] Born, M. Volumen und Hydrationswärme der Ionen. Z. Phys. 1920, 1, 45–48. [48] Honig, B.; Nicholls, A. Classical electrostatics in biology and chemistry. Science 1995, 268, 1144–1149. [49] Nezbeda, I.; Mouˇcka, F.; Smith, W. R. Recent progress in molecular simulation of aqueous electrolytes: force fields, chemical potentials and solubility. Mol. Phys. 2016, 114, 1665–1690. [50] Pollard, T. P.; Beck, T. L. Toward a quantitative theory of Hofmeister phenomena: From quantum effects to thermodynamics. Curr. Opin. Colloid Interface Sci. 2016, 23, 110–118. [51] Li, P.; Merz Jr., K. M. Metal ion modeling using classical mechanics. Chem. Rev. 2017, 117, 1564–1686. [52] Kirkwood, J. G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300– 313. [53] Beck, T.; Paulaitis, M.; Pratt, L. The potential distribution theorem and models of molecular solutions; Cambridge University Press, Cambridge, UK: Cambridge, 2007. 29 ACS Paragon Plus Environment


[54] Rogers, D. M.; Beck, T. L. Modeling molecular and ionic absolute solvation free energies with quasichemical theory bounds. J. Chem. Phys. 2008, 129, 134505/1–134505/15. [55] Beck, T. L. Hydration free energies by energetic partitioning of the potential distribution theorem. J. Stat. Phys. 2011, 145, 335–354. [56] D.M. Rogers, L. P., D. Jiao; Rempe, S. Structural models and molecular thermodynamics of hydration of ions and small molecules. Annu. Rep. Comput. Chem. 2012, 8, 71–127. [57] M. Isegawa, F. N.; Pantazis, D. Ionization energies and aqueous redox potentials of organic molecules: Comparison of DFT, correlated ab initio theory and pair natural orbital approaches. J. Chem. Theory Comput. 2016, 12, 2272–2284. [58] Kastenholz, M. A.; H¨ unenberger, P. H. Computation of methodology-independent ionic solvation free energies from molecular simulations: I. The electrostatic potential in molecular liquids. J. Chem. Phys. 2006, 124, 124106/1–124106/27. [59] Kastenholz, M. A.; H¨ unenberger, P. H. Computation of methodology-independent ionic solvation free energies from molecular simulations: II. The hydration free energy of the sodium cation. J. Chem. Phys. 2006, 124, 224501/1–224501/20. [60] Reif, M. M.; H¨ unenberger, P. H. Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities and expansivities of solvated ions. J. Chem. Phys. 2011, 134, 144103/1–144103/30. [61] Reif, M. M.; H¨ unenberger, P. H. Computation of methodology-independent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones parameter sets for the alkali and halide ions in water. J. Chem. Phys. 2011, 134, 144104/1– 144104/25. [62] Hofer, T. S.; Tran, H. T.; Schwenk, C. F.; Rode, B. M. Characterization of dynamics


Page 30 of 55

Page 31 of 55 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


and reactivities of solvated ions by ab initio simulations. J. Comput. Chem. 2004, 25, 211–217. [63] D’Incal, A.; Hofer, T. S.; Randolf, B. R.; Rode, B. M. Be(II) in aqueous solution—an extended ab initio QM/MM MD study. Phys. Chem. Chem. Phys. 2006, 24, 2841–2847. [64] Hofer, T. S.; Randolf, B. R.; Rode, B. M. The influence of quantum forces on molecular dynamics simulation results for hydrated aluminium(III). Chem. Phys. Lett. 2006, 422, 492–495. [65] Vchirawongkwin, V.; Hofer, T. S.; Randolf, B. R.; Rode, B. M. Tl(I)-the strongest structure-breaking metal ion in water? A quantum mechanical/molecular mechanical simulation study. J. Comput. Chem. 2007, 28, 1006–1016. [66] M. Dal Peraro, P. C., S. Raugei; Klein, M. Solute-solvent charge transfer in aqueous solution. Chem. Phys. Chem. 2005, 6, 1715–1718. [67] Z. Zhao, D. R.; Beck, T. Polarization and charge transfer in the hydration of chloride ions. J. Chem. Phys. 2010, 132, 014502/1–014502/10. [68] Soniat, M.; Rick, S. The effects of charge transfer on the aqueous solvation of ions. J. Chem. Phys. 2012, 137, 044511/1–044511/9. [69] Yao, Y.; Kanai, Y.; Berkowitz, M. L. Role of charge transfer in water diffusivity in aqueous ionic solutions. J. Phys. Chem. Lett. 2014, 5, 2711–2716. [70] Y. Yao, M. B.; Kanai, Y. Communication: Modeling of concentration dependent water diffusivity in ionic solutions: Role of intermolecular charge transfer. J. Chem. Phys 2015, 143, 241101. [71] M. Soniat, L. H.; Rick, S. Charge transfer models of zinc and magnesium in water. J. Chem. Theory Comput. 2015, 11, 1658–1667. [72] M. Ahmed, J. M., A.K. Singh; Sarkar, S. Water in the hydration shell of halide ions has significantly reduced Fermi resonance and moderately enhanced Raman cross section in the OH stretch regions. J. Phys. Chem. B 2013, 117, 9728–9733. 31 ACS Paragon Plus Environment


[73] Li, J.; Wang, F. Accurate prediction of the hydration free energies of 20 salts through adaptive force matching and the proper comparison with experimental references. J . Phys. Chem. B 2017, 121, 6637–6645. [74] Leontyev, I.; Stuchebrukhov, A. Electronic continuum model for molecular dynamics simulations. J. Chem. Phys. 2009, 130, 085102/1–085102/8. [75] Leontyev, I.; Stuchebrukhov, A. Electronic polarizability and the effective pair potentials of water. J. Chem. Theory Comput. 2010, 6, 3153–3161. [76] Leontyev, I.; Stuchebrukhov, A. Accounting for electronic polarization in nonpolarizable force fields. Phys. Chem. Chem. Phys. 2011, 13, 2613–2626. [77] Patra, M.; Karttunen, M. Systematic comparison of force fields for microscopic simulations of NaCl in aqueous solutions: diffusion, free energy of hydration, and structural properties. J. Comput. Chem. 2004, 25, 678–689. [78] Car, R.; Parrinello, M. Unified approach for molecular-dynamics and density functional theory. Phys. Rev. Lett. 1985, 55, 2471–2474. [79] Hutter, J. Car–Parrinello molecular dynamics. WIREs Comput. Mol. Sci. 2012, 2, 604–. [80] Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory, 2nd Edition; Wiley–VCH: Weinheim, 2002. [81] Sholl, D. S.; Steckel, J. A. Density Functional Theory - a practical introduction; Wiley: Hoboken, 2009. [82] Kuo, I.-F.; Mundy, C. An ab initio molecular dynamics study of the aqueous liquidvapor interface. Science 2004, 303, 658–660. [83] McGrath, M. J.; Siepmann, J. I.; Kuo, I. W.; Mundy, C. J. Vapor-liquid equilibria of water from first principles: comparison of density functionals and basis sets. Mol. Phys. 2006, 104, 3619–3626.


Page 32 of 55

Page 33 of 55 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


[84] Leung, K. Surface potential at the air-water interface computed using density functional theory. J. Phys. Chem. Lett. 2010, 1, 496–499. [85] Kathmann, S. M.; Kuo, I. W.; Mundy, C.; Schenter, G. Understanding the surface potential of water. J. Phys. Chem. 2011, B 115, 4369–4377. [86] M. Sulpizi, M. S., M. Salanne; Gaigeot, M.-P. Vibrational sum frequency generation spectroscopy of the water liquid-vapor interface from density functional theory-based molecular dynamics simulations. J. Phys. Chem. Lett. 2013, 4, 83–87. [87] Remsing, R. C.; Baer, M.; Schenter, G. K.; Mundy, C. J.; Weeks, J. D. The role of broken symmetry in solvation of a spherical cavity in classical and quantum water models. J. Phys. Chem. Lett. 2014, 5, 2767–2774. [88] Duignan, T. T.; Baer, M.; Schenter, G. K.; Mundy, C. Electrostatic solvation free energies of charged hard spheres using molecular dynamics with density functional theory interactions. J. Chem. Phys. 2017, 147, 161716/1–161716/11. [89] Schwegler, E.; Grossman, J.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. II. J. Chem. Phys. 2004, 121, 5400–5409. [90] Sit, P.-L.; Marzari, N. Static and dynamical properties of heavy water at ambient conditions from first-principles molecular dynamics. J. Chem. Phys. 2005, 122, 204510/1– 204510/9. [91] Del Ben, M.; Schönherr, M.; Hutter, J.; VandeVondele, J. Bulk liquid water at ambient temperature and pressure from MP2 theory. J. Phys. Chem. Lett. 2013, 4, 3753–3759. [92] Del Ben, M.; Sch¨ utt, O.; Wentz, T.; Messmer, P.; Hutter, J.; VandeVondele, J. Enabling simulation at the fifth rung of DFT: Large scale RPA calculations with excellent time to solution. Comput. Phys. Commun. 2015, 187, 120–129. [93] Del Ben, M.; Hutter, J.; VandeVondele, J. Probing the structural and dynamical prop-



erties of liquid water with models including non-local electron correlation. J. Chem. Phys. 2015, 143, 054506/1–054506/12. [94] Willow, S. Y.; Salim, M. A.; Kim, K. S.; Hirata, S. Ab initio molecular dynamics of liquid water using embedded-fragment second-order many-body perturbation theory towards its accurate property prediction. Sci. Rep. 2015, 5, 14358/1–14358/14. [95] Willow, S. Y.; Zeng, X. C.; Xantheas, S. S.; Kim, K. S.; Hirata, S. Why is MP2-water “cooler” and “denser” than DFT-water. J. Phys. Chem. Lett. 2016, 7, 680–684. [96] M.J. Gillan, D. A.; Michaelides, A. Perspective: How good is DFT for water. J. Chem. Phys. 2016, 144, 130901/1–130901/33. [97] White, A.; Knight, C.; Hocky, G.; Voth, G. Communication: Improved ab initio molecular dynamics by minimally biasing with experimental data. J. Chem. Phys. 2017, 146, 041102/1–041102/5. [98] Soniat, M.; Rogers, D. M.; Rempe, S. B. Dispersion- and exchange-corrected density functional theory for sodium ion hydration. J. Chem. Theory Comput. 2015, 11, 2958– 2967. [99] Ikeda, T.; Boero, M. Role of van der Waals corrections in first principles simulations of alkali metal ions in aqueous solutions. J. Chem. Phys. 2015, 143, 194510/1–194510/7. [100] Arismendi-Arrieta, D. J.; Riera, M.; Bajaj, P.; Prosmiti, R.; Paesani, F. i-TTM model for ab initio-based ion-water interaction potentials. I. Halide-water potential energy functions. J. Phys. Chem. B 2016, 120, 1822–1832. [101] Riera, M.; Götz, A. W.; Paesani, F. The i-TTM model for ab initio-based ion-water interaction potentials. II. Alkali metal ion-water potential energy functions. Phys. Chem. Chem. Phys. 2016, 18, 30334–30343. [102] Grimme, S. Accurate description of van der waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463–1473.


Page 34 of 55

Page 35 of 55 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


[103] Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion contribution. J. Comput. Chem. 2006, 27, 1787–1799. [104] Duignan, T. T.; Baer, M. D.; Schenter, G. K.; Mundy, C. J. Real single ion solvation free energies with quantum mechanical simulations. Chem. Sci. 2017, 8, 6131–6140. [105] Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbenium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227. [106] Aaqvist, J.; Warshel, A. Simulation of enzyme reactions using valence bond force fields and other hybrid quantum/classical approaches. Chem. Rev. 1993, 93, 2523. [107] Field, M. J.; Bash, P. A.; Karplus, M. A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. J. Comput. Chem. 1990, 11, 700. [108] Lyne, P. D.; Hodoscek, M.; Karplus, M. A hybrid QM-MM potential employing HartreeFock or density functional methods in the quantum region. J. Phys. Chem. A 1990, 103, 3462. [109] Warshel, A. Molecular Dynamics Simulations of Biological Reactions. Acc. Chem. Res. 2002, 35, 385. [110] Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 618–622. [111] H¨ aser, M.; Weigand, F. RI-MP2: First derivatives and global consistency. Theor. Chem. Acc. 1997, 97, 331–340. [112] C. H¨ attig, A. H.; K¨ ohn, A. Distributed memory parallel implementation of energies and gradients for second-order Møller-Plesset perturbation theory with the resolutionof-the-identity approximation. Phys. Chem. Chem. Phys. 2006, 8, 1159–1169. [113] Weigend, F.; H¨ aser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143–152. 35 ACS Paragon Plus Environment


[114] Das, S.; Nam, K.; Major, D. T. Rapid convergence of energy and free energy profiles with quantum mechanical size in quantum mechanical−molecular mechanical simulations of proton transfer in DNA. J. Chem. Theory Comput. 2018, 14, 1695–1705. [115] Simonson, T. Free energy of particle insertion. An exact analysis of the origin singularity for simple liquids. Mol. Phys. 1993, 80, 441–447. [116] Beutler, T. C.; Mark, A. E.; van Schaik, R.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529–539. [117] Buelens, F.; Grubm¨ uller, H. Linear-scaling soft-core scheme for alchemical free energy calculations. J. Comput. Chem. 2012, 33, 25–33. [118] Gapsys, V.; Seeliger, D.; de Groot, B. L. New soft-core potential function for molecular dynamics based alchemical free energy calculations. J. Chem. Theory Comput. 2012, 8, 2373–2382. [119] Pham, T.; Shirts, M. Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformations in solution phase. J. Chem. Phys. 2012, 136, 124120/1–124120/14. [120] Tavernelli, I.; Vuilleumier, R.; Sprik, M. Ab initio molecular dynamics for molecules with variable numbers of electrons. Phys. Rev. Lett. 2002, 88, 213002/1–213002/4. [121] Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox potentials and acidity constants from density functional theory based molecular dynamics. Acc. Chem. Res. 2014, 47, 3522–3529. [122] Todorova, M.; Neugebauer, J. Extending the concept of defect chemistry from semiconductor physics to electrochemistry. Phys. Rev. Appl. 2014, 1, 014001/1–014001/15. [123] F. Ambrosio, G. M.; Pasquarello, A. Redox levels in aqueous solution: Effect of van der Waals interactions and hybrid functionals. J. Chem. Phys. 2015, 143, 244508/1– 244508/15. 36 ACS Paragon Plus Environment

Page 36 of 55

Page 37 of 55 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


[124] Frenzel, J.; Meyer, B.; Marx, D. Bicanonical ab initio molecular dynamics for open systems. J. Chem. Theory Comput. 2017, 13, 3455–3469. [125] Vaidehi, N.; Wesolowski, T. A.; Warshel, A. Quantum-mechanical calculations of solvation free energies. A combined ab initio pseudopotential free-energy perturbation approach. J. Chem. Phys. 1992, 97, 4264–4271. [126] Gao, J. Absolute free-energy of solvation from Monte Carlo simulations using combined quantum and molecular mechanical potentials. J. Phys. Chem. 1992, 96, 537–540. [127] Woods, C. J.; Manby, F. R.; Mulholland, A. J. An efficient method for the calculation of quantum mechanics/molecular mechanics free energies. J. Phys. Chem. 2008, 128, 014109/1–014109/8. [128] Shaw, K. E.; Woods, C. J.; Mulholland, A. J. Compatibility of quantum chemical methods and empirical (MM) water models in quantum mechanics/molecular mechanics liquid water simulations. J. Phys. Chem. Lett. 2010, 1, 219–223. [129] E.R. Pinnick, A. R., C.E. Calderon; Wang, F. Achieving fast convergence of ab initio free energy perturbation calculations with the adaptive force-matching method. Theor. Chem. Acc. 2012, 131, 1146/1–1146/11. [130] N.V. Plotnikov, S. K.; Warshel, A. Paradynamics: An effective and reliable model for ab initio QM/MM free-energy calculations and related tasks. J. Phys. Chem. B 2011, 115, 7950–7962. [131] E.C. Dybeck, N. S.; Shirts, M. Effects of a more accurate polarizable Hamiltonian on polymorph free energies computed efficiently by reweighting point-charge potentials. J. Chem. Theory Comput. 2016, 12, 3491–3505. [132] Hofer, T. S.; Randolf, B. R.; Shah, S. A. A.; Rode, B. M.; Persson, I. Structure and dynamics of the hydrated palladium(II) ion in aqueous solution. A QMCF MD simulation and EXAFS spectroscopic study. Chem. Phys. Lett. 2007, 445, 193.



[133] Vchirawongkwin, V.; Persson, I.; Rode, B. M. Structure and Dynamics of Sulfate Ion in Aqueous Solution – an ab initio QMCF MD Simulation and Large Angle X-ray Scattering Study. J. Phys. Chem. B 2007, 111, 4150. [134] Hofer, T. S.; Randolf, B. R.; Rode, B. M.; Persson, I. The hydrated platinum(II) ion in aqueous solution—a combined theoretical and EXAFS spectroscopic study. Dalton Trans. 2009, 1512. [135] Hofer, T. S.; Randolf, B. R.; Rode, B. M. In Solvation Effects on Molecules and Biomolecules; Canuto, S., Ed.; Challenges and Advances in Computational Chemistry and Physics; Springer: Heidelberg, 2008; Vol. 6; Chapter Molecular Dynamics Simulation Methods including Quantum Effects. [136] Eklund, L.; Hofer, T. S.; Pribil, A. B.; Rode, B. M.; Persson, I. On the structure and dynamics of the hydrated sulfite ion in aqueous solution – an ab initio QMCF MD simulation and large angle X-ray scattering study. Dalton Trans. 2012, 41, 5209. [137] Eklund, L.; Hofer, T. S.; Weiss, A. K. H.; Tirler, A. O.; Persson, I. Structure and water exchange of the hydrated thiosulfate ion in aqueous solution using QMCF MD simulation and large angle X-ray scattering. Dalton Trans. 2014, 43, 12711–12720. [138] Eklund, L.; Hofer, T. S.; Persson, I. Structure and water exchange dynamics of hydrated oxo halo ions in aqueous solution using QMCF MD simulation, large angle X-ray scattering and EXAFS. Dalton Trans. 2015, 44, 1816–1828. [139] Hofer, T. S. Perspectives for hybrid ab initio/molecular mechanical simulations of solutions: from complex chemistry to proton-transfer reactions and interfaces. Pure Appl. Chem. 2014, 86, 105–117. [140] Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269–6271. [141] Lin, H.; Truhlar, D. G. QM/MM: what have we learned, where are we, and where do we go from here? Theor. Chim. Acta 2007, 117, 185. 38 ACS Paragon Plus Environment

Page 38 of 55

Page 39 of 55 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


[142] Laio, A.; VandeVondele, J.; Rothlisberger, U. A Hamiltonian electrostatic coupling scheme for hybrid Car–Parrinello molecular dynamics simulations. J. Chem. Phys. 2002, 116, 6941. [143] Voloshina, E.; Gaston, N.; Paulus, B. Embedding procedure for ab initio correlation calculations in group II metals. J. Chem. Phys. 2007, 126, 134115. [144] Roßbach, S.; Ochsenfeld, C. Influence of coupling and embedding schemes on QM size convergence in QM/MM approaches for the example of a proton transfer in DNA. J. Chem. Theory Comput. 2017, 13, 1102−–1107. [145] Barker, J.; Watts, R. Monte Carlo studies of the dielectric properties of water-like models. Mol. Phys. 1973, 26, 789–792. [146] Tironi, I.; Sperb, R.; Smith, P.; van Gunsteren, W. A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 1995, 102, 5451–5459. [147] Martin, F.; Zipse, H. Charge distribution in the water molecule. A comparison of methods. J. Comput. Chem. 2005, 26, 97–105. [148] Dahlgren, B.; Reif, M. M.; H¨ unenberger, P. H.; Hansen, N. Calculation of derivative thermodynamic hydration and aqueous partial molar properties of ions based on atomistic simulations. J. Chem. Theory Comput. 2012, 8, 3542–3564. [149] Reif, M. M.; H¨ unenberger, P. H. Origin of asymmetric solvation effects for ions in water and organic solvents investigated using molecular dynamics simulations: The Swain acity-basity scale revisited. J. Phys. Chem. B 2016, 120, 8485–8517. [150] H¨ unenberger, P. H.; van Gunsteren, W. F. Alternative schemes for the inclusion of a reaction-field correction into molecular dynamics simulations: Influence on the simulated energetic, structural, and dielectric properties of liquid water. J. Chem. Phys. 1998, 108, 6117–6134. [151] Flyvbjerg, H.; Petersen, H. Error estimates on averages of correlated data. J. Chem. Phys. 1989, 91, 461–466. 39 ACS Paragon Plus Environment


[152] Taylor, J. An introduction to error analysis - The study of uncertainties in physical measurements Edition 2.; University Science Books: , Sausalito, CA, USA, 1997. [153] Lawrenz, M.; Baron, R.; McCammon, J. Independent-trajectories thermodynamicintegration free-energy changes for biomolecular systems: Determinants of H5N1 avian influenza virus neuraminidase inhibition by peramivir. J. Chem. Theory Comput. 2009, 5, 1106–1116. [154] Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. [155] Dunning, T. Gaussian basis functions for use in molecular calculations. Contraction of (12s9p) atomic basis sets for the second row atoms. Chem. Phys. Lett. 1970, 7, 423–427. [156] Weigend, F.; M. Haser, H. P.; Ahlrichs, R. RI-MP2: Optimized auxiliary basis sets and demonstration of efficiency. Chem. Phys. Lett. 1998, 294, 143–152. [157] Li, J.; Wang, F. The effect of core correlation on the MP2 hydration free energies of Li+ , Na+ , and K+ . J. Phys. Chem. B 2016, 120, 9088–9096. [158] Ahlrichs, R.; B¨ ar, M.; Häser, M.; Horn, H.; K¨ olmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165–169. [159] von Arnim, M.; Ahlrichs, R. Performance of parallel TURBOMOLE for density functional calculations. j. Comput. Chem. 1998, 19, 1746–1757. [160] Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14, 33–38. [161] M¨ ahler, J.; Persson, I. A study of the hydration of the alkali metal ions in aqueous solution. Inorg. Chem. 2011, 51, 425–438.


Page 40 of 55

Page 41 of 55 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


[162] Bouazizi, S.; Nasr, S. Local order in aqueous lithium chloride solutions as studied by X-ray scattering and molecular dynamics simulations. J. Mol. Struct. 2007, 837, 206– 213. [163] Mason, P.; Ansell, S.; Neilson, G.; Rempe, S. Neutron scattering studies of the hydration structure of Li+ . The Journal of Physical Chemistry B 2015, 119, 2003–2009. [164] Ansell, S.; Barnes, A. C.; Mason, P. E.; Neilson, G. W.; Ramos, S. X-ray and neutron scattering studies of the hydration structure of alkali ions in concentrated aqueous solutions. Biophys. Chem. 2006, 124, 171–179. [165] Lyubartsev, A. P.; Laasonen, K.; Laaksonen, A. Hydration of Li+ ion. An ab initio molecular dynamics simulation. J. Chem. Phys. 2001, 114, 3120–3126. [166] Loeffler, H. H.; Mohammed, A. M.; Inada, Y.; Funahashi, S. Lithium(I) ion hydration: a QM/MM-MD study. Chem. Phys. Lett. 2003, 379, 452–457. [167] Tongraar, A.; Liedl, K. R.; Rode, B. M. The hydration shell structure of Li+ investigated by Born–Oppenheimer ab initio QM/MM dynamics. Chem. Phys. Lett. 1998, 286, 56– 64. [168] Bhattacharjee, A.; Hofer, T. S.; Rode, B. M. Local density corrected three-body distribution functions for probing local structure reorganization in liquids. Phys. Chem. Chem. Phys. 2008, 10, 6653–6657. [169] Chrenko, R.; McDonald, R.; Pell, E. Vibrational spectra of lithium-oxygen and lithiumboron complexes in silicon. Phys. Rev. 1965, 138, A1775. [170] Andrews, L. Infrared spectrum, structure, vibrational potential function, and bonding in the lithium superoxide molecule LiO2 . J. Chem. Phys. 1969, 50, 4288–4299. [171] Hofer, T. S.; Tran, H. T.; Schwenk, C. F.; Rode, B. M. Characterization of dynamics and reactivities of solvated ions by ab initio simulations. J. Comput. Chem. 2004, 25, 211–217.



[172] Marcus, Y. Thermodynamics of solvation of ions. Part 5.—Gibbs free energy of hydration at 298.15 K. J. Chem. Soc., Faraday Trans. 1991, 87, 2995–2999. [173] Lee, A.; Rick, S. The effects of charge transfer on the properties of liquid water. J. Chem. Phys. 2011, 134, 184507/1–184507/9. [174] Kann, Z.; Skinner, J. A scaled-ionic-charge simulation model that reproduces enhanced and suppressed water diffusion in aqueous salt solutions. J. Chem. Phys. 2014, 141, 104507/1–104507/7. [175] qvist, J. A.; Warshel, A. Free energy relationships in metalloenzyme-catalyzed reactions. Calculations of the effects of metal ion substitutions in staphylococcal nuclease. J. Am. Chem. Soc. 1990, 112, 2860–2868. [176] Pang Y.-P., Xu K.; Yazal J. E.; Prendergast, F. Successful molecular dynamics simulation of the zinc-bound farnesyltransferase using the cationic dummy atom approach. Protein Sci. 2000, 9, 1857–1865. [177] Saxena, A.; Sept, D. Multisite ion models that improve the coordination and free energy calculations in molecular dynamics simulations. J. Chem. Theory Comput. 2013, 9, 3538–3542. [178] Duarte, F.; Bauer, P.; Barrozo, A.; Amrein, B. A.; Purg, M.; ˚ Aqvist, J.; Kamerlin, S. C. L. Force field independent metal parameters using a nonbonded dummy model. J. Phys. Chem. B 2014, 118, 4351–4362. [179] Q. Liao, S. K.; Strodel, B. Development and application of a nonbonded Cu2+ model that includes the Jahn-Teller effect. J. Phys. Chem. Lett. 2015, 6, 2657–2662. [180] F. Sessa, L. G., P. D’Angelo; Migliorati, V. Hidden hydration structure of halide ions: An insight into the importance of lone pairs. J. Phys. Chem. B 2015, 119, 15729–15737. [181] Y. Jiang, H. Z.; Tan, T. Rational design of methodology-independent metal parameters using a nonbonded dummy model. J. Chem. Theory Comput. 2016, 12, 3250–3260.


Page 42 of 55

Page 43 of 55 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


[182] Kräutler, V.; Van Gunsteren, W. F.; H¨ unenberger, P. H. A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. j. Comput. Chem. 2001, 22, 501–508. [183] Vesely, F. J. Of pendulums, polymers, and robots: Computational mechanics with constraints. Am. J. Phys. 2013, 81, 537–544. [184] Loeffler, H. H.; Rode, B. M. The hydration structure of the lithium ion. J. Chem. Phys. 2002, 117, 110. [185] Loeffler, H. H.; Inada, Y.; Funahashi, S. Water exchange dynamics of lithium(I) ion in aqueous solution. J. Phys. Chem. B 2006, 110, 5690–5696. [186] Hofer, T. S.; Weiss, A. K. H.; Randolf, B. R.; Rode, B. M. Hydration of highly charged ions. Chem. Phys. Lett. 2011, 512, 139. [187] Pribil, A. B.; Hofer, T. S.; Randolf, B. R.; Rode, B. M. Structure and dynamics of phosphate ion in aqueous solution: an ab initio QMCF MD study. J. Comput. Chem. 2008, 29, 2330–2334. [188] Tongraar, A.; Liedl, K. R.; Rode, B. Born-Oppenheimer ab initio QM/MM dynamics simulations of Na+ and K+ in water: From structure making to structure breaking effects. J. Phys. Chem. A 1998, 102, 10340–10347. [189] Azam, S.; Hofer, T.; Randolf, B.; Rode, B. Hydration of sodium(I) and potassium(I) revisited: A comparative QM/MM and QMCF MD simulation study of weakly hydrated ions. J. Phys. Chem. A 2009, 113, 1827–1834. [190] Bhattacharjee, A.; Pribil, A. B.; Randolf, B. R.; Rode, B. M.; Hofer, T. S. Hydration of Mg2+ and its influence on the water hydrogen bonding network via ab initio QMCF MD. Chem. Phys. Lett. 2012, 536, 39–44. [191] Hofer, T. S.; Randolf, B. R.; Rode, B. M. Influence of polarization and many body quantum effects on the solvation shell of Al(III) in dilute aqueous solution-extended ab initio QM/MM MD simulations. Phys. Chem. Chem. Phys. 2005, 7, 1382–1387. 43 ACS Paragon Plus Environment


◦ Figure 1: The absolute standard intrinsic hydration free energy GH + ,wat of the proton as a central quantity in single-ion thermodynamics. Panel (a) illustrates the connection between ◦ ◦ + GH + ,wat and the absolute intrinsic hydration free energy ∆s GM + ,wat of a cation M , the ab◦ + solute electrode potential VM + ,wat of the Galvanic element involving M and the air-liquid interfacial potential χ◦wat . The equations connecting these quantities are Eqs. 10-13. Panel ◦ (b) shows experimentally-inferred estimates for GH + ,wat as a function of the year of publication. The graph is based on the data of Tab. 5.19 of Ref. 1 The dashed lines indicate the recommended value of -1100 ± 5 kJ·mol−1 obtained from a thorough evaluation of this data ◦ −1 , any set. 1 Since the estimates for GH + ,wat cover a range of about -1225 to -975 kJ·mol ◦ prediction for GH + ,wat located in the grey area of the right bar is within a range of 1% with respect to one of the experimental estimates indicated by the horizontal black lines.


Page 44 of 55

Page 45 of 55 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


Figure 2: Simulated system and basic thermodynamic cycle relevant for the calculation of a single-ion hydration free energy. Panel (a) shows the simulation setup. The snapshot of the periodic simulation cube shows the Li+ ion and 2000 water molecules, highlighting the QM region encompassing the ion and its first hydration shell. Panel (b) shows the basic thermodynamic cycle employed for an ion M+ . The ion M+ (vac) is converted to a noninteracting point mass X(vac), transferred into aqueous solution as X(aq), and transformed to the fully hydrated ion M+ (aq). The calculated value is intrinsic, i.e. excluding the crossing of a vacuum-liquid interface. Note that this cycle is not employed directly in this work, but rather the cycle of Fig. 3.



Figure 3: Thermodynamic cycle applied in the present QM/MM MD/TI calculations. The calculation is split into three successive steps, namely cavitation, charging and quantization. Together with the correction term ∆Gcor and a standard–state adjustment ∆G◦std the sum of the corresponding free-energy contributions yields the absolute intrinsic single-ion hydration free energy ∆s G◦M + ,wat (Eq. 4). The zero energy level defined in Eq. 6 is applied throughout.


Page 46 of 55

Page 47 of 55 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


Figure 4: Lithium-oxygen radial distribution functions (RDFs) and corresponding coordination-number distributions (CNDs). The RDFs (panel a) and CNDs (panels b and c) are obtained from the QM/MM (black) and MM (red) MD simulations of the aqueous Li+ ion at full ionic charge.



Figure 5: Structural and dynamical characteristics of the aqueous Li+ ion. The results are shown for the QM/MM (black) and MM (top panels only, red) simulations of the ion at full ionic charge. The three-body correlation function g 3 (r) (panel a) represents the oxygenoxygen distance distribution within the first shell. The power spectrum of the lithium-oxygen stretching vibration (panel b) was obtained via Fourier transform of the respective corresponding autocorrelation function. The time evolution (panel c) and distribution (panel d) of the ionic partial charge qI was obtained from a Mulliken population analysis.


Page 48 of 55

Page 49 of 55 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


Figure 6: Hamiltonian-derivative curves corresponding to the TI calculation of the cavitation (panel a), charging (panel b) and quantization (panel c) free-energies for aqueous Li+ at 298.15 K and 1 bar. Shown are the average Hamiltonian derivative (black), the corresponding running integrals (red) and the associated error bars (blue). In panel (b), the dashed line connecting the initial and final states of the charging process highlights the slight deviation from a linear dependence.



∆f G◦M + ,vac ∆s G•M + ,wat ◦ ∆H VM + ,wat ◦ Gˆ + H ,wat

(kJ·mol−1 ) (kJ·mol−1 ) (V) (kJ·mol−1 )

H+ 1513.32 0 0 -1087.0±3.2

Li+ 644.98 574.9±1.7 -3.04±0.01

Na+ 570.71 680.5±1.9 -2.72±0.01

Page 50 of 55

K+ 477.63 752.2±0.9 -2.94±0.01

Ref 1 Tab. 5.9 Ref 1 Tab. 5.22 Ref 1 Tab. 5.25 Ref 1 Tab. 5.23

Table 1: Standard experimental values for the non-elusive quantities employed in Eqs. 10 to 13 and in corresponding comparisons to experiment, at 298.15 K and 1 bar. The quantities for an ion M+ (or the proton H+ ) are the gas-phase formation free energy ∆f G◦M + ,vac , the ◦ conventional hydration free energy ∆s G•M + ,wat and the redox potential ∆H VM + ,wat along with ◦ the standard real solvation free energy Gˆ + of the proton. H ,wat


0.41 0.40 0.40

0.194 0.195 0.193 0.196

4.1 4.2 4.1 4.05

Nc,1 4.11 4.05 4.0 4.0-5.0 4.0

13.0 16.0

Nc,2 14.50 13.70

τ (ps) 27.2 19.0

R 4.08 2.84

k (Nm−1 ) 97.4 63.8

78.4 71.9

ν (cm−1 ) 605 471

522.1 500

165

166

166

167

170

169

164

162

163

161

Refs this work this work

Table 2: Selected structural and dynamical properties of the hydrated Li+ ion from theory and experiment: The symbols r1 and r2 denote the lithium-oxygen distances corresponding to the first- and second-shell peaks in the radial distribution function, NC,1 and NC,2 being the corresponding average coordination numbers. The mean residence time τ for water molecules in the first solvation shell was obtained by monitoring exchange events 62 with an allowed excursion time t∗ of 0.5 ps. The quantity ν is the frequency of the ion−oxygen stretch vibration obtained via Fourier transform of the velocity autocorrelation function, k being the associated effective force constant.

r2 (nm) 0.424 0.408

r1 (nm) 0.197 0.196 0.196 0.196 0.198 0.194

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

method MM RIMP2 QM/MM LAXS NS XRS NS IR IR HF QM/MM HF QM/MM B3LYP QM/MM CPMD

Page 51 of 55 Journal of Chemical Theory and Computation

51

ACS Paragon Plus Environment


∆G in kJ·mol−1 ∆Gcav ∆Gchg ∆Gcor ∆G◦std ∆Gqua ∆s G◦M + ,wat MM ∆s G◦M + ,wat QM/MM

Li+ 6.0±0.6 -407.7±0.5 -143.6±1.0 7.95 12.3±0.5 -537.3±2.1 -525.0±2.6

Na+ 10.8±0.4 -300.4±0.9 -143.5±1.0 7.95 7.9±0.8 -425.2±2.3 -417.3±3.1

Page 52 of 55

K+ 15.1±0.3 -223.2±0.6 -143.3±1.0 7.95 0.1±0.5 -343.5±1.8 -343.4±2.3

Table 3: Contributions to the standard absolute intrinsic hydration free energy ∆s G◦wat of aqueous Li+ obtained from the individual steps in the thermodynamic cycle of Fig. 3 under standard conditions (298.15 K and 1 bar; 1 molal ideal-solute for the hydrated ion vs. 1 bar ideal-gas for the gas-phase ion) and calculated at the MM or QM/MM levels. The quantities ∆Gcav , ∆Gchg and ∆Gqua are obtained from TI, the corresponding Hamiltonian derivative curves being shown in Fig. 6. The contributions to ∆Gcor are reported in the Suppl. Mat. Tab. S.2. The standard-state correction ∆G◦std is according to Eq. 8. The individual error bars are reported based on a confidence interval of 99%, the error bar for ∆s G◦wat being estimated as the sum of the individual components. An alternative estimate obtained by propagating the corresponding variances (i.e. square-root of the sum of squared errors) is provided in the Suppl. Mat. Tab. S.3. The corresponding quantities for aqueous Na+ and K+ obtained in a previous study 3 have also been included for comparison.


Page 53 of 55 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


∆s G◦M + ,wat ∆s G◦M + ,wat ∆s G◦M + ,wat ◦ VM + ,wat ◦ VM + ,wat ◦ VM + ,wat ∆∆s G◦M + ,wat ∆∆s G◦M + ,wat ∆∆s G•M + ,wat ◦ ∆VM + ,wat ◦ ∆VM + ,wat ◦ ∆∆H VM + ,wat

MM QM/MM Rec MM QM/MM Rec MM QM/MM Exp MM QM/MM Exp

(kJ·mol−1 ) (kJ·mol−1 ) (kJ·mol−1 ) (V) (V) (V)

Li+ -537.3±2.1 -525.0±2.6 -525.1±6.7 1.12±0.02 1.24±0.03 1.24±0.07

Na+ -425.2±2.3 -417.3±3.1 -419.5±6.9 1.51±0.02 1.59±0.03 1.57±0.07

K+ -343.5±1.8 -343.4±2.3 -347.8±5.9 1.39±0.02 1.39±0.02 1.35±0.06

(kJ·mol−1 ) (kJ·mol−1 ) (kJ·mol−1 ) (V) (V) (V)

Li+ → Na+ 112.1±4.4 107.7±5.7 105.6±3.7 0.39±0.04 0.35±0.06 0.32±0.02

Na+ → K+ 81.7±4.1 73.9±5.4 71.7±2.8 0.12±0.04 0.20±0.05 0.22±0.02

Li+ → K+ 193.8±3.8 181.6±4.8 177.3±2.6 0.27±0.04 0.15±0.05 0.10±0.02

Table 4: Comparison between the calculated and experimental values for the standard hydration and redox properties of aqueous Li+ (this work) as well as Na+ and K+ (based on Ref. 3 ) at 298.15 K and 1 bar. The relevant quantities for an ion M+ are the absolute hydration free ◦ ◦ energy ∆s G◦M + ,wat , the single-electrode potential VM + ,wat , the redox potential ∆H VM + ,wat (relative to the standard hydrogen electrode) and the conventional hydration free energy ∆s G•M + ,wat . An additional ∆-prefix indicates a difference between two ions. The ∆s G•M + ,wat ◦ and ∆H VM + ,wat values for the three ions can be found in Tab. 1. Recommended (Rec) and experimental (Exp) values are from Ref. 1 Values calculated at the MM and QM/MM levels are reported from Tab. 3. All error estimates are obtained by summation of the error bars of the individual free energy contributions. An alternative estimate obtained by propagating the corresponding variances (i.e. square-root of the sum of squared errors) is provided in the Suppl. Mat. Tab. S.4.



◦ GH + ,wat ◦ VH + ,wat χ◦wat

(kJ·mol−1 ) (V) (V)

Rec. -1100±5 4.28±0.13 0.13±0.10

Li+ -1099.9±4.2 4.28±0.04 0.13±0.08

Na+ -1097.8±5.0 4.31±0.05 0.11±0.09

Page 54 of 55

K+ -1095.6±3.2 4.33±0.03 0.09±0.07

Table 5: Comparison between calculated and experimentally-inferred values for the three fundamental elusive quantities of single-ion thermodynamics at 298.15 K and 1 bar. The ◦ quantities are the absolute intrinsic hydration free energy GH + ,wat , the absolute potential ◦ VH + ,wat of the reference hydrogen electrode and the surface electric potential jump χ◦wat upon entering bulk water. The calculated values are derived from the QM/MM MD/TI simulation results for aqueous Li+ (this work) as well as Na+ and K+ (based on Ref. 3 ). All error estimates are obtained by summation of the error bars of the individual contributions. A corresponding estimate of the error bars obtained by propagating the respective variances (i.e. square-root of the sum of squared errors) is provided in the Suppl. Mat. Tab. S.5.


Page 55 of 55 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


ACS Paragon Plus Environment

Single-ion thermodynamics from first principles. Calculation of the

Recommend Documents