Origin of Asymmetric Solvation Effects for Ions in ... - ACS Publications

May 12, 2016 - the asterisk is used to denote a static quantity or a quantity including the static ..... center-of-mass translation of the computation...
0 downloads 0 Views 1MB Size
Subscriber access provided by University of Sussex Library

Article

Origin of Asymmetric Solvation Effects for Ions in Water and Organic Solvents Investigated Using Molecular Dynamics Simulations: The Swain Acity-Basity Scale Revisited Maria M. Reif , and Philippe Henry Hünenberger J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b02156 • Publication Date (Web): 12 May 2016 Downloaded from http://pubs.acs.org on May 17, 2016

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

The Journal of Physical Chemistry B 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 74

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 Journal of Physical Chemistry

Submitted to J. Phys. Chem. B

Origin of Asymmetric Solvation Effects for Ions in Water and Organic Solvents Investigated Using Molecular Dynamics Simulations: The Swain Acity-Basity Scale Revisited Maria M. Reif‡ and Philippe H. H¨ unenberger†,∗ May 11, 2016



Laboratory of Physical Chemistry, ETH Z¨ urich, CH-8093 Z¨ urich, Switzerland



Physics Department (T38), Technische Universit¨at M¨ unchen, D-85748 Garching, Germany

Corresponding author:∗ Prof. Philippe H. H¨ unenberger, ETH Z¨ urich, Laboratorium f¨ ur Physikalische Chemie, ETH– H¨ onggerberg, HCI, CH-8093 Z¨ urich, Switzerland Phone: +41 44 632 5503, Fax: + 41 44 632 1039, e–mail: [email protected] Running title: Asymmetric solvation effects

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 The asymmetric solvation of ions can be defined as the tendency of a solvent to preferentially solvate anions over cations or cations over anions, at identical ionic charge magnitudes and effective sizes. Taking water as a reference, these effects are quantified experimentally for many solvents by the relative acity (A) and basity (B) parameters of the Swain scale. The goal of the present study is to investigate the asymmetric solvation of ions using molecular dynamics simulations, and to connect the results to this empirical scale. To this purpose, the charging free energies of alkali and halide ions, and of their hypothetical oppositely-charged counterparts, are calculated in a variety of solvents. In a first set of calculations, artificial solvent models are considered that present either a charge or a shape asymmetry at the molecular level. The solvation asymmetry, probed by the difference in charging free energy between the two oppositely-charged ions, is found to encompass a term quadratic in the ion charge, related to the different solvation structures around anion and cation, and a term linear in the ion charge, related to the solvation structure around the uncharged ion-sized cavity. For these simple solvent models, the two terms are systematically counteracting each other, and it is argued that only the quadratic term should be retained when comparing the results of simulations involving physical solvents to experimental data. In a second set of calculations, 16 physical solvents are considered. The theoretical estimates for the acity A are found to correlate very well with the Swain parameters, whereas the correlation for B is very poor. Based on this observation, the Swain scale is reformulated into a new scale involving an asymmetry parameter Σ, positive for acitic solvents and negative for basitic ones, and a polarity parameter Π. This revised scale has the same predictive power as the original scale, but it characterizes asymmetry in an absolute sense, the atomistic simulations playing the role of an extra-thermodynamic assumption, and is optimally compatible with the simulation results. Considering the 55 solvents in the Swain set, it is observed that a moderate basity (Σ between -0.9 and -0.3, related to electronic polarization) represents the baseline for most solvents, while a highly variable acity (Σ between 0.0 and 3.0, related to hydrogen-bond donor capacity modulated by inductive effects) represents a landmark of protic solvents.

2 ACS Paragon Plus Environment

Page 2 of 74

Page 3 of 74

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 Journal of Physical Chemistry

I

Introduction

The asymmetric solvation of ions by a given solvent can be defined as the tendency of this solvent to preferentially solvate anions over cations or cations over anions, at identical ionic charge magnitudes and effective sizes.1 By analogy with protolytic equilibria, the preferential affinity of a solvent for anions or for cations are referred to as its acity A˜ and its basity ˜ respectively. 2−5 As a logical extension of this nomenclature, one may say that a solvent B, ˜ and basitic when it preferentially is acitic when it preferentially solvates anions (A˜ > B) ˜ > A). ˜ solvates cations (B The above definition involves at least two implicit assumptions. If acity and basity are to be regarded as intrinsic properties of a given solvent, it requires that: (i) the concept of identical ionic sizes can be formalized, e.g. by means of ionic radii; (ii) the acity and basity of the solvent derived from its affinity for two isovalent, equisized and oppositely-charged ions are independent of the selected ion pair, given an appropriate normalization. The validity of these assumptions are, however, very difficult to assess a priori. The concept of ionic radius is already very ambiguous for monovalent and monoatomic ions 1 (see Table 5.4 and Figure 5.1 therein for a collection of 63 sets of radii proposed to date for the alkali and halide ions), and even more for polyvalent or polyatomic ones. 6−10 The same holds for the concept of affinity, as it relies on the highly non-trivial definition of single-ion solvation free energies. 1, 11

Furthermore, the definition also assumes that the ionic charge (including its sign) and the

ionic size are the only factors governing the magnitude of the solvation forces. This neglects, for example, the different polarizabilities 12−14 (and thus, ion-solvent dispersion interactions) as well as the different hydrogen-bond donor/acceptor capacities 15, 16 (and thus, ion-solvent hydrogen-bonding interactions) between different ions at identical charge and size. In spite of this fundamental ambiguity in their formal definition, acity and basity remain extremely useful qualitative concepts from a practical perspective. For example, by a judicious choice of solvent, experienced organic and analytical chemists often take advantage of asymmetric solvation effects to: 4 (i) alter the position of chemical equilibria, e.g. conformational, complexation, tautomeric or protolytic equilibria; (ii) alter reaction mechanisms and rates and, thereby, product or stereoisomer ratios in reactions under kinetic control; (iii) alter the spectroscopic properties of solute molecules, e.g. light absorption/emission spectra, nuclear magnetic resonance (NMR) or electron spin resonance (ESR) properties. In view of the great practical relevance of solvent effects, a number of empirical scales

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 4 of 74

have been proposed to characterize the solvation properties of liquids (see e.g. Refs. 3−5, 17−19 for overviews). These include in particular the Dimroth ET (30) scale 20, 21 (see also Ref. 12 ), the Mayer-Gutmann AN /DN scale 22−27 (see also Ref. 28 ), the Kamlet-Taft π ∗ /α/β scale 28− 32

(see also Refs. 12, 33 ), the Swain A/B scale, 2, 34 the Catal´an SA/SB/SdP /SP scale, 18 as

well as others.16,35−37 These scales encompass information on asymmetric solvation, along with other key solvent properties 12, 18, 19 such as polarity (including permittivity effects), dipolarity, induction, charge transfer, polarizability (including dispersive effects), hydrogen-bonding capacity (aprotic or protic nature, hydrogen-bond donor/acceptor potential). Furthermore, the information on asymmetry is only included in an indirect way, i.e. not directly based on thermodynamic data pertaining to ionic solvation (see, however, Refs. 32, 38, 39 ). Instead, the scales are constructed by investigating and fitting solvent effects on spectroscopic properties (solvatochromy 3,4 ) or parameters of chemical processes sensitive to the asymmetric solvation of the reactants, products, intermediates or transition states. The recent years have witnessed a renewed attention for these scales of solvation descriptors, in the context of both computer-based data analysis 35−37, 40 and computer-aided process optimization. 41, 42 In the present work, we focus on the Swain acity-basity scale 2, 34 as the one that appears to connect most directly to asymmetric solvation properties. This scale was developed by fitting 1080 experimental values pij concerning 77 types of spectroscopic, thermodynamic or kinetic properties measured in 55 solvents (plus 6 solvent mixtures) to an equation of the form pij = ai Aj + bi Bj + ci

,

(1)

where the i-index refers to a given type of property and the j-index to a given solvent (the latter subscript will be omitted in the following for simplicity). The data considered included equilibrium constants (logarithms) for solvent-solvent partitioning, rate constants (logarithms) and product ratios for solvolysis and other reactions, as well as various spectroscopic properties (light absorption/emission/fluorescence peaks, ESR splitting constants, NMR chemical shifts). The quality of the fit was found to be remarkable, with a linear correlation coefficient of 0.991, which suggested that all solvent effects could be accurately captured by only two solvent-dependent parameters. However, the corresponding parameter vectors are not uniquely defined by Eq. 1, which must be complemented by four critical (scale-setting) conditions. Based on the intuition that the two solvent parameters A and B should measure the acity and basity, respectively, of the solvent, Swain et al. introduced the following critical conditions: (i) A = B = 0 for hexane (assumed non-solvating); (ii) 4 ACS Paragon Plus Environment

Page 5 of 74

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 Journal of Physical Chemistry

A = B = 1 for water (taken as a reference); (iii) A = 0 for hexamethylphosphoric triamide (assumed non-solvating for anions); (iv) B = 0 for trifluoroacetic acid (assumed non-solvating for cations). The resulting A and B values for the 55 pure solvents, taken directly from Ref. 2 (see Table III therein), are reported in Table 1. A number of important comments must be made concerning the Swain acity-basity scale. First, although the fitting involved an extensive experimental data set, it may still be incomplete in terms of the solvation effects probed by the selected spectroscopic properties and chemical processes. 15 Second, given this data set, the fitting accuracy of Eq. 1 is unquestionable, but the assignment of specific A and B values is not unique, as it results from an interpretation in terms of acity and basity via the choice of the critical conditions. 34 Third, the Swain parameters also implicitly encompass information on the solvent polarity (both A and B are expected to increase with polarity) and hydrogen-bonding capacity (A is expected to be higher for hydrogen-bond donor solvents and B for hydrogen-bond acceptor solvents), as well as dipolarity and polarizability. Fourth, the scale provides no information on the ˜ but only on its absolute asymmetry of a solvent (as noted above by the symbols A˜ and B), asymmetry relative to water (as noted here by the symbols A and B). However, although water is characterized by A = B = 1 in the Swain scale, it is in reality a highly acitic solvent ˜ Principally due to hydrogen-bonding effects, it presents a higher affinity for with A˜ > B. anions compared to cations at identical charges and effective sizes, as evidenced by both experiment 43−47 and theory. 48−65 Thus, if a solvent with A > B is more acitic than water, a solvent with A < B may be either less acitic than water, or it may be basitic. This issue represents a major drawback of the Swain scale, which will be tentatively remedied in the present article. Asymmetric solvation effects have also attracted considerable attention in the theoretical community. Although the present article focuses on the role of solvent asymmetry in single-ion solvation (see below), the problem has also been investigated by theory in the context of the solvation of neutral polar 62,66,67 or polarizable 13,14 solutes, of the dynamic properties of ions, 47 , 50, 68, 69

of the interaction between ion pairs, 63, 70, 71 of the properties of ions at interfaces, 63

, 69, 70, 72−76

and of the electrostatic properties of biomolecules. 77−81 In the context of ionic

solvation, theoretical investigations present two major advantages relative to experiment: (i) it is possible to invert the charge of an ion without altering its other interaction properties; (ii) it is possible to consider exactly symmetric dipolar solvent molecules. In nature, because the net charge of an ion and the multipoles of a solvent molecules are intimately related to 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 6 of 74

the atomic numbers and electronic structures of the constituting atoms, ions with opposite charges and otherwise rigorously identical properties as well as perfectly symmetric dipolar solvents do not exist. These specific features of the theoretical approach have been taken advantage of in numerous studies. 49, 50, 61, 72, 82−84 The root of all modern theoretical approaches to single-ion solvation remains the Born model 85 (BO). Proposed in 1920, this model is based on continuum electrostatics and represents a solvated monoatomic ion as a rigid non-polarizable sphere immersed in a homogeneous medium with local and linear dielectric response. The resulting equation for the charging free energy can be derived analytically and reads ∆GBO chg = −

NA 2 −1 q (1 − ǫ−1 S ) RI , 8πǫo I

(2)

where qI is the ionic charge, RI the effective ionic radius, ǫS the static relative dielectric permittivity of the solvent, ǫo the permittivity of vacuum and NA Avogadro’s constant. It is immediately seen that the Born equation neglects asymmetric solvation effects (along with many others 1 ), because Eq. 2 is insensitive to the sign of qI . Accordingly, numerous early attempts to partition experimental hydration free energies measured for neutral sets of ions into absolute intrinsic single-ion quantities 1 (see Table 5.19 and Figure 5.7 therein for a collection of 98 estimates proposed to date for the absolute intrinsic hydration free energy of the proton, required to convert the conventional values into intrinsic ones) also neglected the effect of the water asymmetry. For example, a symmetric partitioning was performed based on a specific pair of monovalent and presumably equisized ions, e.g. Cs+ /I− in Ref.86 , K+ /F− in Refs.87 , K+ /Br− in Ref.88 (with a 0.63:0.37 partitioning ratio), Rb+ /Cl− in Ref. 89 , K+ /Cl− in Ref. 90 , or [Ph4 As]+ /[Ph4 B]− in Refs. 91−93 The latter choice is referred to as the tetraphenylarsonium-tetraphenylborate (TATB) assumption, and the discussion of its validity has attracted considerable attention in the literature. 43, 53, 55, 76, 94−96 In 1939, Latimer et al. 48 introduced a simple modification the Born equation to account for asymmetric hydration effects in an ad hoc fashion, and used it to derive absolute intrinsic single-ion hydration parameters taking these effects into account. The idea was to add an increment ∆RI to the (crystallographic) ionic radius RI in Eq. 2, chosen to be smaller for the anions than for the cations. Since then, a wealth of modified Born equations have been proposed that involve diverse theoretical justifications (see further below), but often result in expressions involving such asymmetric radius increments 97−104 (see also Refs. 105−114 for other types of modified Born equations). In a similar way, numerical continuum-electrostatics 6 ACS Paragon Plus Environment

Page 7 of 74

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 Journal of Physical Chemistry

calculations 115−120 involving solutes of arbitrary geometries and charge distributions, as well as corresponding generalized-Born model approximations, 121−123 have relied on empiricallyoptimized radii taking asymmetric hydration effects into account either implicitly 124 (i.e. during the radius optimization) or explicitly 10, 66, 67, 125−128 (e.g. by modifying the boundary conditions or making the radii depend on the atom or surface charges). If the empirical radius-increment approach has the merit of acknowledging the effect of water asymmetry, it does not provide insight into the physics of this asymmetry, introduces two relatively obscure solvent-dependent quantities (the selected increments), and has seldom been generalized to solvents other than water (see, however, Refs. 99, 129 ). From a theoretical perspective, the basic problem amounts to establishing a connection between the molecular structure of the solvent and its asymmetric solvation properties. A simple view of this connection, illustrated in Figure 1, is that solvation asymmetry is related to an asymmetry of the solvent molecule itself. 113 It can result either from a charge asymmetry (e.g. positively charged region close to the molecular surface, capable of approaching an anion at shorter distances), or from a shape asymmetry (e.g. protruding positively charged regions, capable of packing at a higher density in the solvation shell of an anion). Intuitively, one would expect more exposed positive charges to promote acity and more exposed negative charges to promote basity. To account for these specific features of the solvent molecule, the level of theory has to be increased beyond that of standard continuum electrostatics. This can be done using (quasi-)analytical liquid theories, or using molecular-mechanics (MM), quantummechanics (QM) or hybrid (QM/MM) approaches, as implemented in Monte Carlo (MC) or molecular dynamics (MD) simulations. In the context of (quasi-)analytical liquid theories, 130 asymmetric solvation effects have been investigated in particular using modified Born equations (MBE; see above), non-local electrostatics131−136 (NLE) in Refs.109,114,137−140 , the mean-spherical approximation141−144 (MSA) in Refs. 98−100, 113, 129, 145 , or integral-equation theories 146−151 (IET) in Refs. 51, 73, 152−158 The NLE approach, which features a non-local susceptibility (permittivity) kernel, and the MSA approach, which introduces correlations between hard-sphere solvent molecules, still imply a solvent response that is strictly linear in the ionic charge (and thus exempt of asymmetry; see, however, Ref. 104 for a non-linear model). Yet, they provide a theoretical justification for the application of radii increments, 130 and can be modified to account for the use of distinct increments for anions and cations. 129, 145 In contrast, the IET approach, which fully (albeit approximately) accounts for ion-solvent and solvent-solvent site-site correlations, is 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 8 of 74

non-linear as well as non-local. It is arguably the only one of the four approaches that is able to connect directly an information on the asymmetry of the solvent molecule itself (e.g. based on an interaction potential taken from a classical force field) to its asymmetric solvation properties. 113,154,157 The three other approaches must import this information indirectly, either from distribution functions or solvation free energies monitored in independent atomistic MC or MD simulations 66, 67, 101−103, 110, 111, 139, 140 or from experimentally-inferred intrinsic single-ion solvation free energies. 99, 100, 104, 109, 129, 135, 145 The second procedure is very common but also highly problematic since intrinsic single-ion parameters can only be derived based on experimental data by application of a model, 1, 11 so that one essentially ends up fitting a new model against an older one. Over the last three decades, much insight into asymmetric solvation effects has been gained from classical force-field based MC and MD simulations. One of the main findings in the context of the solvation of spherical ions is that the linear response implied in the Born model should be revised into a piecewise-affine (PA) response. 49, 54, 60, 64, 82, 83, 139 More precisely, the electric potential φ at the ion site as a function of the fractional ionic charge qI′ upon charging appears to obey a relation of the form   φ′ q ′ + I φP A (qI′ ) = φ∗ −  φ′ q ′ − I

if qI′ ≥ 0 if qI′ < 0

,

(3)

where φ∗ is a static potential, already present for the uncharged cavity and constant during the charging process (see Ref. 108 for an early MBE including such a term), while φ′+ and φ′− are (positive) steric potential derivatives, characterizing the distinct linear solvent responses to a positive or a negative charge, respectively. As a result, the charging free energy of an ion up to its full charge qI is given by ∆G∗chg = ∆G∗lin + ∆Gchg

 NA 2  φ′+ if qI ≥ 0 q = N A φ ∗ qI − 2 I  φ′ if q < 0 I −

.

(4)

It follows that the solvation asymmetry, measured by the difference in solvation free energy for two oppositely-charged but otherwise identical ions (cation minus anion and using qI > 0 to characterize the charge of the pair), is given by ∆∆G∗chg = ∆∆G∗lin + ∆∆Gchg = 2NA φ∗ qI +

NA ′ (φ− − φ′+ ) qI2 . 2

(5)

Note that in Eqs. 4 and 5, and throughout the rest of this article, the asterisk is used to denote a static quantity or a quantity including the static contribution. 8 ACS Paragon Plus Environment

Page 9 of 74

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 Journal of Physical Chemistry

Interestingly, a number of MC and MD studies60,64,77 have revealed that, at least to a good approximation, the static potential is essentially: (i) insensitive to the size and shape of the ionic cavity; (ii) homogeneous over the volume of this cavity. Thus, as φ∗ is also by definition independent of the solute charge distribution because it concerns the uncharged ionic cavity, it can essentially be regarded as a property of the solvent that is entirely independent of the nature of the solute. The consistent and accurate calculation of the parameters φ∗ , φ′+ and φ′− in Eqs. 4 and 5 based on MC or MD simulations is a non-trivial task, as the raw results are typically very sensitive to: (i) the approximate treatment of electrostatic interactions during the simulation, e.g. truncated-Coulomb, lattice-sum or reaction-field electrostatics; (ii) the boundary conditions employed for the simulation, e.g. finite droplet under fixed boundary conditions (FBC) or computational box under periodic boundary conditions (PBC); (iii) the method employed for calculating the potential at the ion site based on the solvent configurations, i.e. molecule-based sum (M-summation) or charged-based sum (P-summation). However, as shown by Kastenholz & H¨ unenberger 159, 160 (see also Refs. 1, 84, 161, 162 ), the raw simulation results can be corrected ex post so that methodological independence is achieved and consistent intrinsic single-ion solvation free energies are obtained irrespective of the simulation protocol 160, 161 (see also Refs. 163−165 for discussions of these issues and their correction in binding free-energy calculations involving charged species). The presence of the static potential φ∗ in Eq. 5 raises two important issues. The first one concerns its definition. Following others, 166−170 we have advocated the use of a M-summation convention for the calculation of single-ion solvation free energies, a choice which anchors the zero point for φ∗ in the the subvolume of the bulk liquid located outside the solvent molecules in each configuration. Although we have provided strong arguments supporting this choice in an article 159 of 2006, many recent (and otherwise excellent) studies 64, 77, 154, 171−173 still rely on an alternative P-summation convention, a choice which anchors the zero point for φ∗ in the entire volume of the bulk liquid, including the inside of the solvent molecules. For this reason, we restate our main arguments 159 in favor of M-summation in Appendix A of this article, in a qualitative and intuitive rather than formal manner. The second issue concerns the actual relevance of φ∗ in terms of experimentally observable spectroscopic properties or chemical processes, most of which being bound by the constraint of electroneutrality. 1, 11 Spectroscopic properties usually involve nuclear, rotational, vibrational or electronic excitation processes that do not change the net charge of a solute. And in 9 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

chemical processes, whenever a cation (or a set thereof) is created, an anion (or a set thereof) that has an exactly opposite charge is created simultaneously. If φ∗ accounts for a potential that is, at least as a first approximation, homogeneous over the solute cavity as well as independent of the solute shape, size and charge distribution (see above), it follows that the term ∆G∗lin in Eq. 4 cannot contribute to any experimental observable because it is linear in the solute charge. Exceptions may result from a residual dependence of φ∗ on the size and shape of the solute, from a possible heterogeneity of φ∗ over the solute cavity, or in situations where a spectroscopic measurement involves the absorption or emission of charged subatomic particles (e.g. measurement of electron affinities, 174 ionization potentials 175 or metal work functions 176 ), or where single ions are transferred between different media (e.g. Voltaic cell experiments 1, 130, 177, 178 ). But as a general principle, the above considerations suggest that ∆G∗lin is experimentally irrelevant. In the context of solvation asymmetry, it follows that an experimental scale like the Swain scale most likely characterizes the asymmetry in terms of ∆∆Gchg and not of ∆∆G∗chg . We will nevertheless further consider the two types of ˜ ∗ ) will refer to quantities including quantities. Symbols with an asterisk (e.g. ∆∆G∗chg , A˜∗ , B the linear term, which will be called intrinsic quantities. Symbols exempt of the asterisk (e.g. ˜ B) ˜ will refer to quantities excluding this term, which will be called effective ∆∆Gchg , A, quantities. The symbols A and B will be kept to refer to the Swain relative acity and basity values, which are presumably effective quantities. In the present study, charging free energies calculated from MD simulations are used to investigate solvation asymmetry. This is done first in the context of two simplified solvent models that are characterized by tunable charge or shape asymmetries, as illustrated in Figure 1: (i) a spherical anisotropically-charged model (SAM), which presents a charge asymmetry but is exempt of shape asymmetry (see Refs. 113, 157, 179 for prior usage of SAM-like models to represent the water molecule); (ii) a non-spherical isotropically-charged model (NIM), which presents a shape asymmetry but is essentially exempt of charge asymmetry (see Refs. 150, 155, 156

for prior usage of NIM-like models in the context of IET). Relative charging free energies

of the ions Na+ , K+ , F− , Rb+ and I− in these solvents are compared to those of hypothetical ions Na− , K− , F+ , Rb− and I+ with opposite charges and identical solute-solvent LennardJones interaction parameters. In a second step, the comparison is performed for 16 physical solvent models considering the Na+ and (hypothetical) Na− ions. The results of the latter simulations are compared with the A and B parameters of the Swain acity-basity scale. This comparison leads us to suggest a reformulation of the scale using alternative Σ (or Σ∗ ) and 10 ACS Paragon Plus Environment

Page 10 of 74

Page 11 of 74

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 Journal of Physical Chemistry

Π parameters, leading to the definition of an asymmetry-polarity scale instead.

II

Computational details

II.1

Molecular dynamics simulations

The MD simulations were performed using the GROMOS11 program, 180 as well as with a modified version of the GROMOS96 program181,182 implementing the particle-particle particlemesh 183, 184 (P3 M) lattice-sum method for the calculation of electrostatic interactions. They were all carried out under PBC based on cubic computational boxes. The simulations aiming at the calculation of an ionic charging free energy (Section II.2) involved computational boxes containing one ion (charge state qI′ ) and a number of NS solvent molecules (Table 2). The ions considered include the alkali ions Na+ , K+ and Rb+ (full ionic charge qI = +1 e) as well as the halide ions F− and I− (qI = −1 e), with ion-water Lennard-Jones interaction185 iw ) previously reoptimized 84 for the simple point charge (SPC) waparameters (C6iw and C12

ter model 186 against conventional single-ion hydration free energies by Tissandier et al. 187 + −1 for the standard intrinsic hydration free along with a value ∆G

hyd [H ] = −1100 kJ·mol

energy of the proton 1 (set L in Ref. 84 ). In addition, corresponding hypothetical ions with reversed charges, i.e. Na− , K− and Rb− (qI = −1 e) as well as F+ and I+ (qI = +1 e), were also considered, with the same ion-water Lennard-Jones interaction parameters as the corresponding physical ions but with opposite charges. The latter ions serve as a test for the effect of anisotropic solvation at constant effective ionic radius, and are not meant to provide realistic models for the corresponding physical (but rather uncommon 188−190 ) ions. The derivation of Lennard-Jones interaction parameters for ion-solvent interactions in the various solvents considered strictly followed a geometric-mean combination rule, 181, 191, 192 i.e. is = (C ii C ss )1/2 , where i denotes the ion and s a solvent interaction C6is = (C6ii C6ss )1/2 and C12 12 12 ii coefficients were initially derived from the C iw and C iw site. In particular, the C6ii and C12 6 12 ss appropriate for the SPC water values of Ref. 84 (see above) given the coefficients C6ss and C12 ww ) and subsequently recombined with the C ss and C ss parameters model 186 (C6ww and C12 6 12

associated with the different sites of the solvent models. Note that GROMOS distinguishes ss values for different types of interactions 181 (non-polar, polar-uncharged three alternative C12

and oppositely-charged) and, if available for a given atom of the solvent model, the second value was selected as appropriate for interactions involving an ion and a neutral solvent ii (different ions) as well as C ww and C ww molecule. The values of the parameters C6ii and C12 6 12

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 12 of 74

ss for (SPC water model 186 ) are reported in Table 3. The values of the parameters C6ss and C12

the different sites of the physical solvent models can be found in the corresponding literature references (see below). The solvent models considered are listed in Table 2, namely: (i) artificial solvent models referred to as spherical anisotropically-charged solvent models (SAMs); (ii) artificial solvent models referred to as non-spherical isotropically-charged solvent models (NIMs); (iii) water (H2O) described according to the SPC model; 186 (iv) the dimethylsulfoxide (DMS) model of Ref. 193 ; (v) the dimethylacetamide (DAMD) model of Ref. 194 ; (vi) the methanol (MET) model of Ref. 195 ; (vii) the ethanol (ETL) model of Ref. 194 ; (viii) the acetone (PPN) model of Ref. 194 ; (ix) the isopropanol (2PL) model of Ref. 194 ; (x) the pyridine (PYRI) model of Ref. 194 ; (xi) the acetic acid (ACA) model of Ref. 194 ; (xii) the ethylacetate (EAE) model of Ref. 194 ; (xiii) the butanamine (BAN) model of Ref. 194 ; (xiv) the chloroform (CHL) model of Ref. 196 ; (xv) the diethylether (DEE) model of Ref. 194 ; (xvi) the toluene (TOL) model of Ref. 194 ; (xvii) the benzene (BZN) model of Ref. 194 ; (xviii) the carbontetrachloride (CCL) model of Ref. 197 Note that the CCL model has no partial charges so that no simulations were required in practice for this solvent. A SAM (Figure 1) consists of a van der Waals sphere with the same Lennard-Jones interacww ), and encompassing a first (centered) tion parameters as a SPC water oxygen186 (C6w and C12

particle (charge −qs , mass 15.999 g·mol−1 ) and a second (off-center) particle (charge qs , mass 1.008 g·mol−1 ), the latter offset by a constrained distance b and not involved in LennardJones interactions. The two particles are excluded from any mutual non-bonded interaction. Twenty different SAMs were considered, with values of qs ranging from 0.1 to 0.5 e in steps of 0.1 e and values of b ranging 0.06 to 0.12 nm in steps of 0.02 nm. Note that a SAM with b = 0 or qs = 0 would induce no electrostatic solvation (uncharged Lennard-Jones sphere). A NIM (Figure 1) consists of two van der Waals spheres, each associated with a centered point charge. The first (unscaled) site (charge −qs , mass 15.999 g·mol−1 ) has the same ww ). The second Lennard-Jones interaction parameters as a SPC water oxygen186 (C6ww and C12

(scaled) site (charge qs , mass 15.999 g·mol−1 ) has the homoatomic Lennard-Jones interaction parameters C6 = α6 C6ww

and

ww C12 = α12 C12

,

(6)

where α is a radius-scaling factor. The geometric-mean combination rule 181, 191, 192 was also applied to deduce the Lennard-Jones interaction parameters for all interactions involving this second site. The two sites of the NIM are separated by a constrained distance b = 0.2 nm. 12 ACS Paragon Plus Environment

Page 13 of 74

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 Journal of Physical Chemistry

The two particles are excluded from any mutual non-bonded interaction. Twenty different NIMs were considered, with values of qs ranging from 0.1 to 0.5 e in steps of 0.1 e and values of α ranging 0.4 to 0.7 in steps of 0.1. Note that a NIM with qs = 0 would induce no electrostatic solvation (uncharged Lennard-Jones diatomic molecule) and a NIM with α = 1 would present no asymmetric solvation properties (purely dipolar molecule). Considering Figure 1, it is easily seen that a SAM presents a charge but no shape asymmetry, while a NIM presents a shape but essentially no charge asymmetry. The word “essentially” is added in the latter statement considering that the concept of charge symmetry (relative to the molecular envelope) is somewhat ambiguous in the presence of a shape asymmetry. Since the asymmetric solvation properties of a SAM or NIM with qs < 0 can easily be deduced from those of the corresponding SAM or NIM with an opposite qs , the simulations only involved SAMs and NIMs with qs > 0. However, the analysis considers both situations. For all MD simulations, the equations of motion were integrated using the leap-frog scheme 198 with a timestep of 2 fs. The fixed geometry of the rigid solvent models (SAMs, NIMs, H2O, DMS, MET, CHL and CCL) and the bond-length distances within otherwise flexible solvent models (all other solvents) were enforced by application of the SHAKE algorithm 199 with a relative geometric tolerance of 10−4 . The center of mass translation of the computational box was removed every 2 ps to avoid the accumulation of translational velocity caused by algorithmic noise. All simulations were performed at constant volume and temperature (canonical NVT ensemble). The temperature was maintained at 298.15 K by weak coupling to a heat bath 200 using a coupling time of 0.1 ps. The edge length L of the cubic computational box was selected in the following way. For the simulations involving physical solvents, it corresponded to that of a system containing the solvent and Na− , equilibrated at a constant pressure of 1 bar during 0.1 ns. Note that performing this equilibration with Na+ instead would result in a negligible change of L (data not shown). For the simulations involving the SAMs, a density ρ′S = 922 kg·m−3 was selected, corresponding to the equilibrium density of SPC water at 298.15 K and 1 bar corrected for the slightly lower molar mass of the SAM compared to water. For the simulations involving the NIMs, a density ρ′S = 1735 kg·m−3 was selected, likewise corresponding to the equilibrium density of SPC water corrected for the almost twofold molar mass of the NIM. In all cases, the box-edge length L and the number NS of solvent molecules in the box are connected to

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 14 of 74

the equilibrium density ρ′S of the solvent model via ρ′S = NA−1 MS L−3 (NS + 1)

,

(7)

where MS is the molar mass of the solvent and the ionic volume is approximated by that of one solvent molecule. 52 For the physical solvents, this equation was used to calculate ρ′S . For the SAMs and NIMs, it was used to calculate L. The values of NS , L and ρ′S are reported in Table 2. Three different electrostatic interaction schemes were employed during the simulations: (i) lattice summation with tinfoil boundary conditions 183, 184 (LS); (ii) cutoff truncation of Coulombic interactions including a Barker-Watts reaction-field correction 201, 202 using a molecule-based cutoff 203 (BM); (iii) straight truncation of Coulombic interactions using a molecule-based cutoff 203 (CM). The LS scheme (P3 M) was applied with 183, 184, 204−210 a spherical hat charge-shaping function of width 1.0 nm, a triangular-shaped-cloud assignment function, a finite-difference scheme of order two, and a grid-spacing of about 0.1 nm (simulations involving the SAMs and NIMs or the SPC water model) or about 0.06-0.09 nm (all other simulations). This choice of parameters led to a ratio of the root-mean-square error in the atomic forces to the average norm of these forces 184, 207, 208 of about 0.2-0.5% for the different systems. The BM and CM schemes were applied using a twin-range setup,211 with short- and long-range cutoff distances set to 0.8 and 1.4 nm, respectively, and an update frequency of 5 timesteps for the short-range pairlist and intermediate-range interactions. The BM scheme was applied with values ǫRF for the relative permittivity of the dielectric continuum surrounding the cutoff sphere set to the permittivity ǫ′S of the solvent model (Tables 2, 4 and 5; see also Suppl. Mat. Section S.I). All production simulations (different ions, ionic charge states, solvents and electrostatic schemes; Section II.2) were preceded by an equilibration period of 100 ps (SAMs, NIMs) or 500 ps (all other solvents). The configurations sampled along these simulations were written to file every 0.2 ps for subsequent structural analysis, while the corresponding energetic data (electric potential at the ion site) was written every timestep (SAMs, NIMs) or also every 0.2 ps (all other solvents).

II.2

Free-energy calculations

∗,raw Raw charging free energies ∆Gchg exempt of self-term contribution for all electrostatic

= 0 in the absence of schemes considered, 84, 160 i.e. calculated in such a way that ∆G∗,raw chg 14 ACS Paragon Plus Environment

Page 15 of 74

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 Journal of Physical Chemistry

solvent, were obtained as follows. ∗,raw For the SAMs and NIMs, ∆Gchg was calculated by reversibly charging a neutral ion-

sized cavity according to the method described by Hummer 52 and applied in the previous work by Kastenholz, Reif & H¨ unenberger. 84, 160 This approach, related to thermodynamic integration,212 requires the sampling of the average solvent-generated electric potential at the ion site and of its fluctuations for different ionic charge states qI′ from zero to the full ionic charge qI . Three charge states (qI′ /qI = 0.0, 0.5 and 1.0) were used for all ions. Simulations at each charge state were carried out for 1 ns after 0.1 ns equilibration. ∗,raw For the physical solvents, ∆Gchg was calculated by reversibly charging a neutral ion-sized

cavity using thermodynamic integration,212 which requires the sampling of the average solventgenerated electric potential at the ion site for different ionic charge states qI′ from zero to the full ionic charge qI . The calculations relied on 21 equispaced charge states (qI′ /qI = 0.0, 0.05, ..., 1.0) and the integration was performed using the trapezoidal rule. Here, the simulations at each charge state were carried out for 0.15 ns (rigid solvent models) or 0.5 ns (flexible solvent models) after 0.1 ns equilibration. The quantity ∆G∗,raw was calculated for all cations and anions (Table 3) in the SAMs chg and NIMs (Tables 4 and 5) using the LS scheme only. It was calculated for Na+ and Na− in the physical solvents (Table 2) using the BM and CM schemes (all solvents) as well as the LS scheme (rigid solvent models only).

II.3

Principle of the free-energy correction scheme

As shown by Kastenholz & H¨ unenberger 159, 160 (see also Refs. 1, 84, 161 ), the raw charging free energies ∆G∗,raw (Section II.2), which are typically very sensitive to the employed simulation chg methodology, can be corrected ex post so that methodological independence is achieved. The corresponding additive corrections may be evaluated using continuum electrostatics and analytical models, or empirically-fitted versions thereof, 161 and must account for (see Refs. 1, 84, 160, 161

and references therein):

(A) The deviation of the solvent polarization around the ion relative to the polarization in an ideal Coulombic system and the incomplete interaction of the ion with the polarized solvent, a consequence of possible approximations made in the representation of electrostatic interactions during the simulation (e.g. non-Coulombic interactions involving cutoff truncation, possibly with a shifting or switching function or a reaction-field cor-

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

rection). This type-A correction is evaluated for the given approximate electrostatic scheme used in the simulation, but in the idealized context of a macroscopic system. For cutoff-based schemes, it can be subdivided into a type-A1 correction, acting beyond the cutoff sphere of the ion, and a type-A2 correction, acting within this sphere. There is no such correction for LS schemes, which are Coulombic in the limit of infinite box sizes. (B) The deviation of the solvent polarization around the ion relative to the polarization in an ideal macroscopic system, a consequence of the use of a finite system during the simulation (e.g. droplet under FBC or computational box under PBC). This typeB correction is evaluated for the given approximate electrostatic scheme used in the simulation. (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 (see Appendix A), a consequence of the possible application of an inappropriate summation scheme (along with a possibly non-Coulombic potential function) for the contributions of individual solvent atomic charges to this potential (i.e. P-summation, considering individual solvent charges instead of M-summation, considering charges within individual solvent molecules), as well as of the possible presence of a constant offset in this potential (e.g. due to an interfacial potential jump at the surface of a FBC system or to the constraint of vanishing average potential in LS schemes). This type-C correction is evaluated for the given approximate electrostatic scheme and choice of boundary conditions used in the simulation, and can be subdivided into a type-C1 correction for improper potential summation and a type-C2 correction for the potential offset. (D) An inaccurate dielectric permittivity of the employed solvent model. When the corrected charging free energies ∆G∗chg are increased by an estimate ∆Gcav for the free energy of cavity formation and by a term ∆G

std to match the appropriate definition of standard states (possibly differing for gas phase and solution), the scheme leads to standard absolute single-ion solvation free energies that are consistent across a wide range of simulation protocols. 84, 160 These values are intrinsic 1, 11, 160 (rather than real), i.e. they account for bulk solvation effects only and exclude the work of crossing an air-liquid interface. In the present work, the terms ∆Gcav and ∆G

std were not calculated, because they do not contribute to 16 ACS Paragon Plus Environment

Page 16 of 74

Page 17 of 74

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 Journal of Physical Chemistry

asymmetric solvation effects, i.e. they are equal for anions and cations of identical sizes in a given solvent.

II.4

Application of the free-energy correction scheme

Following the principles of the correction scheme (Section II.3), the raw charging free energies ∆G∗,raw (Section II.2) were used to calculate corresponding corrected (methodologychg independent) values ∆G∗chg for the different ions and solvents, as 1, 84, 160, 161 ∆G∗chg = ∆G∗,raw chg + ∆Gcor

,

(8)

where ∆Gcor is the free-energy correction for the different (methodology-dependent) errors committed during the calculation (points A-D in Section II.3). In the specific case of the LS, BM and CM schemes, this correction term can be written 161 (see Eq. 39 therein) ∆Gcor

= +





∆GA1 |qI | , ǫ′S , RC

BM,CM

∆GA2 |qI | , RI , ǫ′S , RC

BM,CM

 + ∆GB |qI | , RI , L, ǫ′S [, RC ]BM LS,BM  + ∆GC1 qI , RI , L, NS [, ǫ′S , RC ]BM , γ ′S  LS + ∆GC2 qI , RI , L, χ′S , χ ˜′S  . + ∆GD |qI | , RI , ǫ′S , ǫS

(9)

The quantities involved in this expression are the full ionic charge qI (±1.0 e), the effective ionic radius RI (see below), the (long-range) cutoff distance RC (1.4 nm), the number NS of solvent molecules in the (cubic) computational box (Table 2), the corresponding box-edge length L (Table 2), the relative dielectric permittivity ǫ′S of the employed solvent model (Tables 2, 4 and 5) and the corresponding experimental value ǫS , the effective quadrupolemoment trace γ ′S of the solvent model (see below), and its air-liquid interfacial potential χ′S (in the air-to-liquid direction) along with the dependence χ ˜′S of this quantity on the inverse curvature radius for a concave interface (see below). For the BM scheme, it is assumed that the reaction-field permittivity ǫRF is equal to ǫ′S , as is the case in the present simulations (Section II.1). The terms ∆GA1 and ∆GA2 are only required (non-vanishing) when using cutoff-based schemes (CM or BM), to correct for the omission of interactions between the ion and the polarized solvent beyond the cutoff distance and to correct for the error in the solvent polarization within the cutoff sphere of the ion, respectively. The term ∆GA1 is analytical (Born 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 18 of 74

equation), = −(8πǫo )−1 NA qI2 (1 − ǫ′S ∆GBM,CM A1

−1

−1 )RC

.

(10)

The term ∆GA2 can in principle be evaluated numerically based on the method of Refs. 213 Here, an empirical (parameterized) equation from Ref. 161 (see Eqs. 42 and 43 therein) was used instead for the approximate evaluation of this term. The term ∆GB is required to correct for the finite size of the computational box. For the LS scheme, this term is given by −1

∆GLS = (8πǫo )−1 NA qI2 (1 − ǫ′S )L−1 B "     # 4π RI 2 16π 2 RI 5 − × αLS + 3 L 45 L

,

(11)

where αLS ≈ −2.837297 is the LS self-term constant.52,184,204,214 For the BM and CM schemes, this term can in principle be evaluated numerically based on the methods of Refs. 215,216 Here, an empirical (parameterized) equation from Ref. 161 (see Eqs. 20 and 50-53 therein) was used instead for the approximate evaluation of ∆GB . The term ∆GC1 is only required (non-vanishing) for the LS and BM schemes. It corrects the P-summation implied by these schemes to a proper M-summation, and is given by   4πRI3 LS ξS′ (12) ∆GC1 = −NA qI 1 − 3L3 or ∆GBM C1

2(ǫ′S − 1) = −NA qI 2ǫ′S + 1



R3 1 − 3I RC



ξS′

,

(13)

where ξS′ is the exclusion potential of the solvent model, which is a function of NS , L and γ ′S (see Eq. 15 below). In Eqs. 12 and 13, the factor between the parentheses is typically very close to unity, because RI ≪ L or RI ≪ RC in all practically relevant situations. For polar solvents with ǫ′S ≫ 1 the fractional prefactor in Eq. 13 is also very close to one. In these limiting situations, ∆GC1 = −NA qI ξS′ accounts for a shift in the zero of the potential (within the bulk solvent) by an amount that is equal to the exclusion potential ξS′ of the solvent model (see further below and Appendix A). Finally, the term ∆GC2 is only required (non-vanishing) for the LS scheme. It accounts for the fact that the boundary condition on the LS electric potential is that of a vanishing value over the box, rather than over the solvent-occupied region of this box. Even at zero ionic charge, the presence of a potential jump at the surface of the ionic cavity induces a small difference between the two types of averages, which must be corrected. The term ∆GC2 depends 18 ACS Paragon Plus Environment

Page 19 of 74

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 Journal of Physical Chemistry

on the solvent model, as well as on pressure and temperature, through the magnitude and curvature dependence of the surface potential of the liquid at a concave air-liquid interface. It is relatively difficult to calculate explicitly. 159 However, because this term is proportional to the ratio of the ion volume to the box volume according to Ref. 161 (see Eqs. 31 and 33 therein; note that the latter equation incorrectly contains a term in ∆GB ), its magnitude is very small for the systems considered here. For this reason, this term was neglected in the calculation of ∆Gcor . Finally, the correction of type D, related to the possibly inaccurate permittivity ǫ′S of the solvent model relative to the experimental permittivity ǫS of the liquid 161 (see Eq. 35 therein), is omitted in the present work because it is also typically very small (as long as ǫ′S −1 is reasonably close to ǫ−1 S ). Two parameters involved in Eq. 9 deserve special attention, namely the effective radius RI of an ion and the effective quadrupole moment trace γ ′S of a solvent model. The effective radius of an ion is not an unambiguously defined quantity. 1 However, previous investigations 84 have shown that the correction term ∆Gcor is only moderately sensitive to a specific choice for this parameter, within reasonable bounds. For all calculations, RI was defined here as in previous work 84, 161 by the arithmetic mean of the location Rrdf of the first peak in the ion-solvent (reference site) radial distribution function and the Goldschmidt in-crystal radius217 Rgld , based on the suggestion that this definition of RI represents a close-to-optimal Born radius in solution. 102 The choices of reference sites and corresponding values of Rrdf are reported in Tables 2, 4 and 5 for the different ions and solvents considered. The values of Rgld for the ions can be found in Table 3. The effective quadrupole moment trace γ ′S of the solvent model is required for the calculation of the correction term ∆GC1 . This term is proportional to the exclusion potential ξS′ of the solvent model. As detailed elsewhere160,161 (see also Appendix A), the appearance of these quantities in the correction term ∆GC1 arises from the fact that in the LS scheme, the electric potential at the ion site is measured relative to the average electric potential in the periodic computational box. Because this averaging region also includes the interior of the solvent molecules, the zero of the potential encompasses a spurious and arbitrary component related to the average potential generated by the solvent charges within the molecules themselves. This spurious offset is removed by application of the correction term ∆GC1 . A similar artifact arises in the BM scheme, leading to the requirement for an analogous correction. 159 For a fully rigid solvent model with a single Lennard-Jones interaction site, it is easily shown 218 19 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 20 of 74

that the exclusion potential is given as in Ref. 161 (see Eq. 21 therein) ξS′ = NA (6ǫo )−1 MS−1 ρ′S γ ′S

,

(14)

where γ ′S = γS′ is the quadrupole moment trace of the solvent model relative to its LennardJones interaction site, which is readily calculated from the knowledge of the geometry and charge distribution of the solvent model, MS is the molar mass of the solvent, and ρ′S the density of the solvent model in the simulated system. In the present simulations, the density is directly related to MS , NS and L via Eq. 7, so that one may write ξS′ = (6ǫo )−1 L−3 (NS + 1)γ ′S

.

(15)

For simplicity, Eqs. 14 and 15 are also retained in the situation of a rigid or flexible solvent model with multiple Lennard-Jones interaction sites, by introducing an effective quadrupole moment trace γ ′S . However, in this case, γ ′S is no longer trivially related to molecular properties of the solvent model, i.e. it becomes a collective property. This quantity must be evaluated numerically via MD simulation, e.g. by performing a sampling in the so-called orientational disorder limit 159 (ODL), by introducing biasing potentials (electric-potential restraint), 219 or by comparing CM and BM (or LS) raw electric potentials in neutral cavities (or corresponding raw charging free energies). The details of these calculations, which are not trivial, are provided in Suppl. Mat. Section S.II. The γ ′S estimates obtained using different approaches are compared in Suppl. Mat. Tables S.1 (NIMs) and S.2 (physical solvents). The estimates considered most reliable are reported in Tables 2, 4 and 5 for the different solvent models.

II.5

Analysis of the solvation asymmetry

The analysis of the solvation asymmetry was performed based on simulations using the LS scheme for the NIMs, the SAMs and the rigid physical solvent models (H2O, DMS, MET and CHL), or using the CM scheme for the flexible solvent models. Simulations using the BM or CM schemes for the rigid physical solvent models and using the BM scheme for the flexible solvent models were performed as well, to verify the methodological independence of the corrected charging free energies. The corrected charging free energies ∆G∗chg of the physical cations or anions (Na+ , K+ , F− , Rb+ and I− ) and of their hypothetical oppositely-charged counterparts (Na− , K− , F+ , Rb− and I+ ) were calculated in the SAMs and NIMs (all ions) as well as in the physical solvents (Na+ and Na− only). A list of all calculations performed 20 ACS Paragon Plus Environment

Page 21 of 74

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 Journal of Physical Chemistry

can be inferred from Tables 4 (SAMs, all with LS), 5 (NIMs, all with LS) and 6 (physical solvents). As mentioned earlier, the results for the SAMs and NIMs with qs < 0 can be deduced from the simulations involving the corresponding model solvents with qs > 0, and the results for CCL did not require any simulation as this model has zero partial charges. In quantitative terms, the absolute acity of a given solvent could be defined as the standard solvation free energy of a reference anion X − in the solvent, and the absolute basity as the corresponding value for the equisized cation X + . In this case, the asymmetry of the solvent in the sense of Eq. 5 corresponds to the free-energy difference ∆∆G∗chg = ∆G∗chg [X + ] − ∆G∗chg [X − ]

,

(16)

where ∆G∗chg [X + ] = ∆G∗lin [X + ] + ∆Gchg [X + ] = F Φ∗ + ∆Gchg [X + ]

(17)

and ∆G∗chg [X − ] = ∆G∗lin [X − ] + ∆Gchg [X − ] = −F Φ∗ + ∆Gchg [X − ]

,

(18)

F being the Faraday charge. Here, ∆G∗lin represents the static contribution arising from the interfacial potential already exerted by the solvent when the ion is uncharged, i.e. a contribution linear in the ionic charge, while ∆Gchg represents the steric contribution arising from the change in solvent orientation and packing around the ion occurring during its charging, i.e. a contribution quadratic in the ionic charge. Note that when the charging free energies of the anion and the cation are averaged, the linear contribution cancels out, i.e. 1 1 ∗ ∆Gchg = {∆Gchg [X + ] + ∆Gchg [X − ]} = {∆G∗chg [X + ] + ∆G∗chg [X − ]} = ∆Gchg 2 2

(19)

For simplicity, the notation ∆Gchg will be used for this average in the following. The quantity ∆∆G∗chg may also be partitioned into two contributions, as ∆∆G∗chg = ∆∆G∗lin + ∆∆Gchg

,

(20)

where ∆∆G∗lin = ∆G∗lin [X + ] − ∆G∗lin [X − ] = 2F Φ∗

(21)

and ∆∆Gchg = ∆Gchg [X + ] − ∆Gchg [X − ] = ∆∆G∗chg − ∆∆G∗lin

.

(22)

The term ∆∆G∗lin accounts for the fact that if a solvent generates a positive potential Φ∗ at the center of an uncharged ion-sized cavity, the contribution of this potential to the asymmetry ∆∆G∗chg will be positive, i.e. it will favor anion over cation solvation, the opposite 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 22 of 74

being true for a negative potential. However, this linear term does not arise from a response of the solvent to the ion charge, but solely from the interfacial properties of the pure solvent at contact with a cavity. In contrast, ∆∆Gchg accounts for the differential response of the solvent to the charging of a positive vs. a negative ion, after removal of the uncharged-cavity contribution. Upon reorientational response to the charging, a solvent that is capable of increasing more significantly the density of positive charges close to the surface of an anion than the density of negative charges close to the surface of a cation, will be characterized by a positive ∆∆Gchg contribution to ∆∆G∗chg , i.e. this term will favor anion over cation solvation, the opposite being true in the inverse situation. However, because they have very distinct physical origins, the contributions ∆∆G∗lin and ∆∆Gchg need not automatically be correlated with each other. The electric potential Φ∗ is required for the partitioning of the charging free energies into linear and quadratic contributions. In practice, it can be calculated based on the corresponding raw electric potential Φ∗,raw monitored during the simulation at zero ionic charge by applying a correction Φ∗cor that only involves the corresponding free-energy terms linear in qI , as described in Ref. 161 (see Eq. 32 therein; note that this equation incorrectly contains a term in ∆GB ) Φ∗ = Φ∗,raw + Φ∗cor = Φ∗,raw +

1 N A qI

∆GC1

.

(23)

For the CM scheme, ∆GC1 = 0 so that Φ∗,CM =0 . cor

(24)

For the LS scheme, using Eq. 12, one has Φ∗,LS cor

  4πRI3 ξS′ =− 1− 3L3

.

(25)

For the BM scheme, using Eq. 13, one has Φ∗,BM cor

2(ǫ′ − 1) = − ′S 2ǫS + 1

  RI3 1 − 3 ξS′ RC

.

(26)

The quantities ∆∆G∗chg or ∆∆Gchg are expected to provide absolute measures of the asymmetry, positive for an acitic solvent and negative for a basitic solvent. However, they depend in magnitude on the selected ion pair. To remove this dependence, one may try to connect the calculated charging free energies to the relative acities A and basities B introduced in the Swain scale. 2, 34 In this case, the relative acity A of a solvent S could be defined as the charging free energy of the reference anion X − in the solvent, divided by the 22 ACS Paragon Plus Environment

Page 23 of 74

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 Journal of Physical Chemistry

corresponding value in a reference solvent R, and the basity B as the corresponding quotient for the equisized cation X + , namely A∗ =

∆G∗chg,S [X − ] ∆G∗chg,R [X − ]

and

B∗ =

and

B=

∆G∗chg,S [X + ]

,

∆G∗chg,R [X + ]

(27)

or A=

∆Gchg,S [X − ] ∆Gchg,R [X − ]

∆Gchg,S [X + ] ∆Gchg,R [X + ]

.

(28)

Taking the reference solvent R to be water, the above definitions would be compatible with the Swain scale. More precisely, taking the definition of Eq. 28 as an example, it would correspond to applying Eq. 1 considering the process of charging the ion pair X − + X + (at infinite separation) in solvent S ∆Gchg,S [X − + X + ] = ∆Gchg,R [X − ]A + ∆Gchg,R [X + ]B + 0 .

(29)

Note that the above equations have been formulated in terms of charging free energies rather than standard solvation free energies, because the cavity-formation and standard-state correction terms do not contribute to the asymmetry. As mentioned in Section I, the linear component of the charging free energy should probably be disregarded as well and the fact that Eq. 29 considers the charging process of a neutral ion pair is in line with this suggestion. Nevertheless, the alternative definitions of intrinsic parameters A∗ and B ∗ (Eq. 27) and effective parameters A and B (Eq. 28) will both be further considered.

III III.1

Results and Discussion Model solvents

The properties of the SAMs and NIMs required for the evaluation of the free-energy correction terms (Section II.4) are reported in Tables 4 and 5, respectively. For the SAMs (Table 4), the polarity of the solvent, as estimated by its permittivity ǫ′S , increases systematically with both qs and b, approximately as the third to fourth power of the dipole moment qs b. When qs = 0.5 e and b = 0.12 nm, the solvent becomes so polar that it is in a glassy state, and its permittivity can no longer be evaluated accurately. Except for this model, the SAM series spans the range 1.3-109.5 in terms of ǫ′S . For these solvent models, which are rigid and characterized by a single Lennard-Jones interaction site, the effective quadrupole-moment trace γ ′S is given analytically by the quadrupole-moment trace γS′ = qs b2 . 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

For the NIMs (Table 5), the value of ǫ′S increases approximately quadratically with qs and is only weakly sensitive to α. Here also, the highly-polar solvents with qs = 0.5 e and α = 0.6 or 0.7 are in a glassy state, so that their permittivities cannot be evaluated accurately. Except for these two models, the NIM series spans the range 4.6-99.6 in terms of ǫ′S . For these solvent models, which are rigid but characterized by two Lennard-Jones interaction sites, the effective quadrupole-moment trace γ ′S must be evaluated numerically as detailed in Suppl. Mat. Section S.II. Three numerical approaches 159, 219 were applied to determine γ ′S , which give comparable results (see Suppl. Mat. Table S.1). The γ ′S from the method considered most reliable are reported in Table 5, and were used to calculate the correction terms. As expected, the γ ′S values for a given α are essentially proportional to qs . For a given qs , they also appear to depend only weakly on α in the range considered here. As a result, γ ′S is 2

approximately equal to qs b with b = 0.197 nm. This effective value b is essentially identical to the actual distance b = 0.2 nm between the two particles. The results of the simulations involving the ion pairs Na+ /Na− , K+ /K− , F+ /F− , Rb+ /Rb− and I+ /I− in the SAMs with different values of the parameters qs and b, and in the NIMs with different values of the parameters qs and α, are reported graphically in Figures 2 and 3, respectively. The corresponding numerical values and calculation details can be found in Suppl. Mat. Tables S.3 and S.4, respectively. In the figures, the ion series is presented in order of increasing ion size (see the effective radii σ ˜ iw in Table 3) and the graphs also include SAMs and NIMs with qs < 0. It is interesting to correlate this information with structural features concerning the average distance and radial orientation of the SAMs and NIMs around the anions, the uncharged cavities and the cations. The latter information is provided graphically in Figure 4, which displays the locations of the first peaks corresponding to the two solvent sites in the corresponding radial distribution functions, considering the SAMs and NIMs with qs > 0 around the smallest (Na− , Na0 , Na+ ) and the largest (I− , I0 , I+ ) cavities only. Considering the SAMs (Figure 2), the average charging free energy ∆Gchg over the cation and the anion is always negative. It is zero for qs = 0 or b = 0 (uncharged Lennard-Jones sphere), and increases in magnitude with increasing |qs |, with increasing b, and with decreasing ion size. The enhancement of solvation upon increasing |qI | or b is related to the increasing polarity of the SAM, in good agreement with the corresponding increase in the molecular dipole moment qs b and dielectric permittivity ǫ′S (Table 4). The effect of ion size is immediately rationalized based on the Born model (Eq. 2). The slight irregularities in the 24 ACS Paragon Plus Environment

Page 24 of 74

Page 25 of 74

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 Journal of Physical Chemistry

curves for the SAM with qs = ±0.5 e and b = 0.12 nm are probably related to the formation of a glassy state (see above). The intrinsic solvation asymmetry ∆∆G∗chg as defined by Eq. 16 is an odd and monotonically increasing function of qs (except for the I+ /I− pair), i.e., as expected, the SAMs for which the off-center charge (closest to the molecular surface) is positive tend to preferentially solvate anions over cations, the opposite being true for the SAMs with a negative off-center charge. The asymmetry also increases in magnitude with increasing b and with decreasing ion size. However, it is also seen that the dependence of ∆∆G∗chg on qs evolves from quasilinear (Na+ /Na− ) to highly damped (Rb+ /Rb− ) and even non-monotonic (I+ /I− , with an extremum at qs = ±0.2 e) upon increasing the ion size. The latter behavior is easily understood when considering the linear and quadratic (effective) contributions ∆∆G∗lin and ∆∆Gchg , respectively, to ∆∆G∗chg as defined by Eq. 20. The quadratic contribution ∆∆Gchg is typically the largest and follows similar trends as ∆∆G∗chg (see above), but without damping effect upon increasing the ion size. On the other hand, the linear contribution ∆∆G∗lin is an odd and monotononically decreasing function of qs , i.e. anticorrelated with ∆∆Gchg . Furthermore, although the magnitude of this term increases with b, the curves are nearly identical for all ion sizes. For the smallest ions, ∆∆Gchg clearly dominates over ∆∆G∗lin in ∆∆G∗chg . However, for the largest ions, the magnitudes of the two terms become comparable, especially at high |qs |. For example, for the I+ /I− pair, the two terms partly cancel out, leading to a very limited solvation asymmetry for all the SAMs, and to the observed non-monotonic dependence on |qs |. Nevertheless, for all ion sizes and SAM parameter combinations considered and in terms of the intrinsic asymmetry ∆∆G∗chg , the SAMs with qs > 0 (exposed positive charge) preferentially solvate anions and the SAMs with qs < 0 (exposed negative charge) preferentially solvate cations, in agreement with intuitive expectations based on Figure 1. The above trends are intimately connected with the average distance and radial orientation of the SAMs around an anion, an uncharged cavity and a cation, as illustrated in Figure 4 for the sodium-sized cavity in the case qs > 0. Close to an anion, the positivelycharged off-center site of the SAM points towards the ion and the distance from the ion to the central site is the smallest, due to the strong resulting electrostrictive force. Close to a cation, the off-center site points away from the ion and the distance to the central site is longer, consequence of a reduced electrostrictive force. Close to an uncharged cavity, the off-center site also points on average away from the ion (although less pronouncedly) and the 25 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

distance to the central site is even larger in the absence of electrostrictive force. The latter situation is illustrated in more details in Figure 5a, which shows the orientation correlation function cI1 (r) and the radial distribution functions gI1 (r) and gI2 (r) for the cavity Na0 and the SAM with b = 0.12 nm and qs = 0.4 e, where I is the ion, 1 the central site and 2 the off-center site (see also Suppl. Mat. Figure S.4 for the anion and cation cases as well as the iodine-sized cavity). In agreement with the observations based on Figure 4 and the negative potential of -0.53 V in the cavity, the cI1 curve is positive at the ion surface, indicative of a preferential orientation involving the off-center site pointing on average towards the bulk liquid. This specific residual orientation of the SAMs around a cavity is responsible for the negative ∆∆G∗lin value when qs > 0. This discussion is easily transposed to the case qs < 0 where the off-center site still points preferentially towards the bulk liquid, now resulting in a positive ∆∆G∗lin . Considering the NIMs (Figure 3), the average charging free energy ∆Gchg is again always negative. It is zero for qs = 0 (uncharged Lennard-Jones diatomic molecule) and increases in magnitude with increasing |qs |, with decreasing α, and with decreasing ion size. Note that NIMs with α = 1 (purely dipolar diatomic molecules; not considered here) would be solvating entities presenting ∆Gchg values similar to those of the corresponding NIMs with α = 0.7, but they would present no asymmetry whatsoever (∆∆G∗lin = ∆∆Gchg = 0). As for the SAMs, the above observations are easy to rationalize based on the trends observed in the dielectric permittivity ǫ′S (Table 5) and on the Born model (Eq. 2). The slight irregularities in the curves for the NIMs with qs = ±0.5 e and α = 0.6 or 0.7 are probably related to the formation of a glassy state (see above). The intrinsic solvation asymmetry ∆∆G∗chg is again an odd function of qs . However, here, it nearly systematically presents extrema at qs = ±0.1 e (except for α = 0.7). For low α and |qs | values, and more pronouncedly for the smallest ions, the NIMs for which the protruding charge is positive expectedly tend to preferentially solvate anions over cations. In contrast, for high α or |qs | values, and more pronouncedly for the largest ions, the opposite and rather counter-intuitive preference is observed. This irregularity arises from the cumulative effects of the linear and quadratic (effective) contributions ∆∆G∗lin and ∆∆Gchg , respectively, to ∆∆G∗chg . As was the case for the SAMs, ∆∆Gchg is a monotonically increasing function of qs . The dependence of this function on qs is more pronounced for the lowest α values and for the smallest ions. Also similarly to the SAMs, ∆∆G∗lin is a monotonically decreasing function of qs that is nearly independent of 26 ACS Paragon Plus Environment

Page 26 of 74

Page 27 of 74

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 Journal of Physical Chemistry

the ion size. Here, the dependence of this function on qs is also essentially insensitive to the second parameter α of the solvent model. However, in sharp contrast with the SAMs, it is the contribution ∆∆G∗lin that is typically the largest in magnitude (except for the smallest |qs | charges, α values and ion sizes), resulting in a generally opposite asymmetric solvation behavior. In terms of intrinsic asymmetry, most of the NIMs with qs > 0 (exposed positive charge) preferentially solvate cations and most of the NIMs with qs < 0 (exposed negative charge) preferentially solvate anions, which is clearly at odds with intuitive expectations based on Figure 1. The above trends are intimately connected with the average distance and radial orientation of the NIMs around an anion, an uncharged cavity and a cation, as illustrated in Figure 4 for the case qs > 0. Close to an anion, the positively-charged scaled site points towards the ion and the distance from the ion to the unscaled site is relatively large, due to the space occupied by the tightly packed unscaled sites in contact with the ion. Close to a cation, the scaled site points away from the ion and the distance to the unscaled site is relatively short, because these sites are now directly in contact with the ion. The situation is somewhat more complicated close to an uncharged cavity. Here, Figure 4 would suggest that the scaled site can be either closer (most often) or further away (especially for high ǫ′S ) than the unscaled site, while the distance from the ion to the unscaled site is intermediate between the anion and the cation situations. However, as illustrated in more details in Figure 5b for the cavity Na0 and the NIM with α = 0.4 and q = 0.2 e (see also Suppl. Mat. Figure S.5 for the anion and cation cases as well as the iodine-sized cavity), this observation is misleading. Although the first peak in gI2 indeed occurs at a closer distance (0.20 nm) compared to that in gI1 (0.29 nm), the cI1 curve is positive at the ion surface, indicative of a preferential orientation involving the scaled site pointing on average towards the bulk liquid, in agreement with the negative potential of -0.56 V in the cavity. The reason why the analysis of gI1 and gI2 instead of cI1 may be misleading is given in Appendix A in the context of the SPC water model 186 (see Figure 5c and Suppl. Mat. Figure S.6). The specific residual orientation of the NIMs around a cavity is responsible for the positive ∆∆G∗lin value when qs < 0 (and negative value when qs > 0). Although the SAMs and NIMs represent highly simplified models for physical solvent molecules, the above analysis already permits to draw a few general conclusions concerning asymmetric solvation effects. The intuitive picture of Figure 1 appears to be correct in terms of the effective asymmetry ∆∆Gchg . In terms of this quantity, a more exposed positive charge 27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

(either off-center or protruding) promotes acity, due to an increased ability to concentrate this positive charge close to an anion. Conversely, a more exposed negative charge promotes basity, due to an increased ability to concentrate this negative charge close to a cation. Furthermore, the close similarity between the shapes of the ∆∆Gchg curves for different ion pairs in Figures 2 and 3 suggests that this quantity can be normalized by an ion-dependent factor so as to become a unique property of the solvent. For the SAMs, the intrinsic asymmetries ∆∆G∗chg always present the same signs as the corresponding effective asymmetries ∆∆Gchg , and are thus also in line with intuitive expectations based on Figure 1. This is no longer systematically the case for the NIMs where, in most situations (except for the lowest α and |qs | values for the smallest ions), an opposite preference is observed. Furthermore, for both SAMs and NIMs, the dissimilarity between the shapes of the ∆∆G∗chg curves for different ion pairs in Figures 2 and 3 suggests that this quantity cannot be normalized by an ion-dependent factor so as to become a unique property of the solvent. The reason for this difference is the presence of the linear term ∆∆G∗lin entering into ∆∆G∗chg , which is generally smaller (SAMs) or generally larger (NIMs) in magnitude compared to ∆∆Gchg , appears to be nearly independent of the ion size, and presents opposite trends relative to ∆∆Gchg . The observation that the trends in ∆∆G∗lin and ∆∆Gchg are opposite for the simple solvent models considered here is easy to understand. A solvent with more exposed positive charge is expected to preferentially orient this exposed charge towards an anion and away from a cation, resulting in an acitic behavior in terms of ∆∆Gchg . However, in contact with an uncharged cavity, this solvent is expected to preferentially orient the exposed charge towards the solvent (more polar region) than towards the cavity (non-polar region), resulting in a basitic component in terms of ∆∆G∗lin . Taken together, the combined effect of the linear and quadratic terms can be viewed as resulting from a competition between the ion and the solvent for the most polar ends of the solvent molecules at their interface. At this point, the question arises whether the intrinsic quantity ∆∆G∗chg or the effective quantity ∆∆Gchg is the appropriate one to correlate with experimental estimates for the solvation asymmetry as encompassed e.g. in the Swain scale. 2,34 In Section I, we argumented in favor of the effective quantity based on the electroneutrality conservation principle. The analysis of the SAMs and NIMs provides a second argument. Figures 2 and 3 show that although both ∆∆Gchg (after an appropriate normalization to account for the different ionic sizes) and ∆∆G∗lin (nearly independent of the ionic size) can be made properties of the sol28 ACS Paragon Plus Environment

Page 28 of 74

Page 29 of 74

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 Journal of Physical Chemistry

vent only, their sum ∆∆G∗chg cannot. In other words, if the experimental asymmetry was measured by ∆∆G∗chg , it would be impossible to construct a scale such as the Swain scale where the parameters A and B characterize the solvent only, independently of the solute. To illustrate that the latter issue is also relevant in the context of a physical solvent, the quantities ∆∆Gchg and ∆∆G∗chg were also calculated for the five ion pairs considering the SPC water model 186 . The results can be found in Suppl. Mat. Table S.6. They show that 3 ∆∆G both ∆∆G∗lin and Reff chg are essentially independent of the ion pair, where Reff is the 3 used for normalization was chosen empirically; for mean of Rgld and Rrdf (the factor Reff

example the Latimer radius-increment model 48 suggests that the leading term of ∆∆Gchg −2 ). In contrast, ∆∆G∗chg evolves from positive for Na+ /Na− to negative for I+ /I− . is in Reff

Thus, if acity and basity were defined based on ∆∆G∗chg , one would have to conclude that water is acitic for the Na+ /Na− pair, but basitic for the I+ /I− pair.

III.2

Simulation results for the physical solvents

The properties of the physical solvent models required for the evaluation of the free-energy correction terms (Section II.4) are reported in Table 2. The exclusion potentials ξS′ of the models or, equivalently, their effective quadrupole-moment traces γ ′S (related to the latter via Eq. 14) can be determined analytically for rigid solvent models with a single Lennard-Jones interaction site. This applies here to the H2O and CCL models, where γ ′S is simply equal to the quadrupole-moment trace γS′ of the model relative to this site. For the other solvents, γ ′S must be evaluated numerically. As detailed in Suppl. Mat. Section S.II, five different types of methods (labelled I, II, IIP, IIM, III, IV, V+ and V-) were considered for the evaluation of γ ′S . The results for a subset of combinations of solvent models and methods are compared in Suppl. Mat. Table S.2. For the solvent CCL, no calculation is required as γ ′S is analytically and numerically zero in any case (no partial charges). For the solvent H2O, the analytical or numerical estimates of all methods are in excellent agreement. For the solvents DMS, MET and CHL, which are rigid but involve multiple Lennard-Jones interaction sites, the estimates of all numerical methods are in the same ranges, except for one (method IIP). Finally, for the flexible solvent models, where only two methods were tested (methods IV and V), it appears difficult to find a consensus. Tentative estimates of γ ′S considered most reliable and retained for the free-energy calculations, are listed in Table 2. Clearly, the methods proposed here to determine γ ′S for the flexible solvent models are not entirely satisfactory and the estimates reported in Table 2 for these solvents should be 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

considered preliminary at best. Yet, two important points should be kept in mind: (i) the appropriateness of the γ ′S value can be tested explicitly by assessing the methodological independence of the corrected charging free energies (e.g. identical corrected results for the CM, BM and LS schemes); (ii) there is one scheme, the CM scheme, which does not require knowledge of γ ′S for the evaluation of the correction terms. Thus, although the LS and BM schemes are commonly regarded as more accurate than the CM scheme in the context of classical simulations, 209, 210 the latter scheme may actually be advantageous for calculating ionic solvation free energies in non-aqueous solvents, because it does not involve the extra complication of evaluating γ ′S and gives results that are as accurate as those obtained using BM and LS electrostatics after application of the correction scheme. 84, 160, 161 The results of the free-energy calculations in the physical solvents considering the Na− /Na+ ion pair are summarized in Table 6. The corresponding raw results and calculation details can be found in Suppl. Mat. Table S.5. For CCL, the results are analytical (no simulation). For H2O, DMS, MET and CHL, the results obtained using the LS, BM and CM schemes are compared. They are found to be consistent with deviations of at most 5.6 kJ·mol−1 (MET with Na− for LS vs. BM, error amounting to about 1% of the charging free energy). For these solvents, the corrected charging free energies are thus essentially methodology independent, and the LS results will be considered for further discussion. For the other (flexible) solvents, only the CM and BM schemes are compared. The agreement is not as good, with deviations that can reach up to 17.6 kJ·mol−1 (2PL with Na− , error amounting to about 3.5% of the charging free energy). As mentioned above, the BM estimates rely on the somewhat uncertain value of γ ′S whereas the CM estimates do not. For this reason, the CM results will be considered for further discussion. The results concerning the solvation magnitude and asymmetry for the Na− /Na+ ion pair in the 16 physical solvents are illustrated graphically in Figure 6. The successive panels show the charging free energies ∆G∗chg of the anion and of the cation, the intrinsic asymmetry ∆∆G∗chg (Eq. 16), and the quadratic (effective) and linear contributions ∆∆Gchg and ∆∆G∗lin , respectively, to this quantity (Eqs. 21 and 22). The solvents are shown in order of decreasing experimental permittivity ǫS (Table 2). As expected on the basis of the Born model85 (Eq. 2), the magnitude of the solvation free energy tends to decrease with decreasing solvent permittivity. This trend is more systematic for the Na+ cation compared to the Na− anion, where many irregularities are observed. More specifically, the protic solvents H2O, MET, ETL, 2PL and ACA stand out as presenting particularly low solvation free energies 30 ACS Paragon Plus Environment

Page 30 of 74

Page 31 of 74

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 Journal of Physical Chemistry

for Na− . The intrinsic solvation asymmetry ∆∆G∗chg evidences no obvious correlation with the permittivity of the solvent. For example, the two solvents with the highest permittivities, H2O and DMS, exhibit opposite asymmetric solvation properties. This is not surprising considering that in spite of the comparably high polarities, H2O is a protic solvent which is expected to enhance acity, whereas DMS is a non-protic solvent (solely hydrogen-bond acceptor) which is expected to enhance basity. Note that the solvent CCL bears no partial charges in the classical representation employed here. As a result the charging free energy of an ion in this solvent is zero and there cannot be any asymmetry. This special case will not be discussed further in the following. As was observed for the SAMs but not always for the NIMs (Section III.1), the effective contribution ∆∆Gchg dominates the intrinsic asymmetry ∆∆G∗chg , except for the aromatic solvents TOL and BZN. The trends in ∆∆Gchg essentially follow the intuitive picture of Figure 1. Water H2O as well as the alcohols MET, ETL and 2PL, the carboxylic acid ACA, and the solvent chloroform CHL are characterized by positive and large ∆∆Gchg values (range 93.9 to 278.2 kJ·mol−1 ), and are thus acitic. These solvents are all protic, i.e. they present a positive charge concentrated on a relatively exposed hydrogen atom (two atoms for H2O) and negative charges that are more diffuse and less exposed. All the other solvents are characterized by negative and large ∆∆Gchg values (range −224.6 to −77.0 kJ·mol−1 ), and are thus basitic. None of these solvents is protic, given that one does not consider aliphatic, amine (BAN) and aromatic (PYRI, TOL, BZN) hydrogen atoms as protic components. The solvents DMS, DAMD, PPN, PYRI, EAE, BAN and DEE all present a concentrated negative charge on an oxygen or a nitrogen atom, along with positive charges that are more diffuse and less exposed. The solvents BZN and TOL (as well as PYRI) are aromatic and thus capable of favorable interactions between the negatively charged aromatic carbon atoms and a cation, so-called cation-π interactions. 220−223 Generally, the term ∆∆G∗lin represents a marginal contribution to the intrinsic asymmetry ∆∆G∗chg , at most about 20% for BAN, except for the solvents H2O, MET, TOL and BZN. For H2O and MET, ∆∆G∗lin is negative, thereby counteracting the positive ∆∆G∗lin , and amounts to about 40% of ∆∆G∗chg . For the aromatic solvents TOL and BZN, ∆∆G∗lin is positive, and about as large in magnitude as the negative term ∆∆Gchg . As a result, ∆∆G∗chg is nearly zero for these solvents. In contrast to the SAMs and NIMs (Section III.1), ∆∆Gchg and ∆∆G∗lin are not systematically of opposite signs. For the solvents DAMD, PPN, 2PL, 31 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

CHL and DEE, ∆∆Gchg and ∆∆G∗lin actually have the same sign. This suggests that the rules governing the preferential orientation of the solvent molecules at a cavity are more complex in the case of (possibly flexible) physical solvents presenting more than two interaction sites. In particular, the more diffuse and flexible nature of the charge distributions in these physical solvents appears to result in comparatively smaller ∆∆G∗lin contributions. The two exceptions H2O and MET with large and negative ∆∆G∗lin (of opposite sign to ∆∆Gchg ) involve only three sites and are rigid, i.e. arguably closest to the SAMs and NIMs. The two exceptions TOL and BZN, with large and positive ∆∆G∗lin (also of opposite sign to ∆∆Gchg ), probably orient perpendicularly to the cavity interface, as they would likely do around an anion (anion-aromatics hydrogen-bonding 221,224−226 ), whereas they are expected to orient in a tangential fashion around a cation (cation-π interactions 220−223 ).

III.3

Comparison with the Swain parameters

The results of the simulations in terms of the calculated charging free energies for the Na− and Na+ ions (Table 6), for simplicity retaining only one reference electrostatic scheme for the calculation (LS for the rigid models or CM for the flexible models), can be converted to simulation estimates for the Swain acity and basity parameters. Two possible equations were considered to this purpose, resulting in the pair of intrinsic estimates A∗ and B ∗ (Eq. 27) or the pair of effective estimates A and B (Eq. 28). The resulting values can be found in Table 7 along with the experimentally-inferred parameters A and B of the Swain scale, 2, 34 and the comparison is performed graphically in Figure 7. The correlation between the simulation results for A or A∗ and the acity parameter A of Swain is quite remarkable. For A, the linear correlation coefficient is 0.93, with a slope of 0.97 (i.e. very close to unity) and an intercept of 0.09 (i.e. very close to zero). For A∗ , the corresponding values are 0.93, 0.96 and 0.12, respectively. In spite of a solvent that is somewhat off the line (BAN), the excellent agreement suggests that: (i) the A-parameter is indeed an appropriate measure of acity; (ii) the solvent models considered 186, 193−197 are accurate (functional form and parameters), in particular in terms of representing anion-solvent interactions. In contrast, the correlation between the simulation results for B or B ∗ and the basity parameter B of Swain is rather poor. Except for four among the five least polar solvents considered in the simulations, the points essentially fall onto a nearly horizontal line. The four 32 ACS Paragon Plus Environment

Page 32 of 74

Page 33 of 74

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 Journal of Physical Chemistry

remaining solvents (CHL, TOL, BZN and CCL) are characterized by much lower simulation values which, arguably, could be explained by the lack of explicit electronic polarizability in the models (expected to be particularly relevant for low-polarity solvents 227 ). This absence of correspondance between experiment and simulations could have two origins: (i) the Bparameter is not an appropriate measure of basity; (ii) the solvent models considered 186, 193− 197

are inaccurate in terms of representing cation-solvent interactions. When considering the above results, it is interesting to note that the Swain scale, 2 in-

troduced in 1983, has led to some debate. In an article 15 of 1984, Taft et al., authors of an alternative scale of solvation properties based on solvatochromic effects,28−32 have argued that the parameter B is in fact a measure of the dipolarity/polarizability, and not of basity. In a reply 34 published in the same journal issue, Swain provides some very interesting clarifications. First, he observes that the A and B parameters provide a remarkable account of the 1080 experimental observables considered in the parameterization, 2 with an overall correlation coefficient of 0.991. Furthermore, in contrast to other scales including that of Taft et al., the Swain scale involves only two solvent parameters and a single common equation, i.e. no observable-specific conditions are required for deciding about the applicability of a given parameter. In summary, for the 1080 observables considered, the solvent effects appear to be strictly and accurately accounted for by only two solvent parameters. The intuitive interpretation of Swain for this key observation is that one parameter A governs acity, acidity and electrophilicity, whereas the other parameter B governs basity, basicity and nucleophilicity. Based on this interpretation, so-called critical (scale-setting) conditions are introduced to fix the the A and B scales. For a two-parameter scale, four critical conditions are required, which are chosen to be A = B = 1 for water, A = 0 for hexamethyl phosphoric triamide (assumed non-solvating for anions) and B = 0 for trifluoroacetic acid (assumed non-solvating for cations). However, the author also stresses that if the accuracy of the fit is unquestionable, the critical conditions are only fixed by an interpretation and not by the experimental data itself. In fact, considering Eq. 1, it is easily seen that arbitrary affine transformations of the original set of A and B parameters into a new set of parameters A′ and B ′ will leave the quality of the fit invariant, merely requiring a change in the a, b and c coefficients for each observable considered. Based on the above considerations and the observation that, according to our simulations, A appears to correlate with acity but B does not seem to correlate with basity, we suggest in the next section a specific affine transformation of the Swain A,B-scale. This transformation 33 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 34 of 74

results in a new Σ, Π-scale which has the same fitting accuracy with respect to the 1080 observables of Swain, but where the parameters Σ and Π now account for asymmetry and polarity in a way that is optimally compatible with the simulation results.

III.4

Revisiting the Swain scale

The Σ, Π-scale is defined based on the Swain A,B-scale by the affine transformations Σ = −0.77 + 2.17A − 0.40B

(30)

and Π = 0.40 + 0.37A + 0.45B

.

(31)

We also define a Σ∗ , Π-scale variant by replacing Eq. 30 with Σ∗ = −0.82 + 2.72A − 0.90B

.

(32)

Since A = B = 1 for water, one also has Σ = Σ∗ = 1 for this solvent, along with Π = 1.22. The values of Σ, Σ∗ and Π for the 61 solvents and solvent mixtures considered by Swain are reported in Table 1. For the 55 pure solvents, these values are also illustrated graphically in the three lower panels of Figure 8. The three upper panels of this figure show the permittivity ǫS and the original A and B parameters of Swain. Simulation estimates for the 16 solvents considered in the present study are also reported in the panels for Σ, Σ∗ and Π. These were calculated for each solvent S based on the simulation results for the Na− /Na+ pair (Table 6) considering one selected electrostatic scheme (LS for the rigid models or CM for the flexible models), as Σ=

∆∆Gchg,S ∆∆Gchg,R

Σ∗ =

,

∆∆G∗chg,S ∆∆G∗chg,R

(33) (34)

and Π=

∆Gchg,S ∆Gchg,R

,

(35)

where ∆∆G∗chg and ∆∆Gchg are the intrinsic and effective asymmetries (Eqs. 16 and 22), ∆Gchg is the average charging free energy of the anion and the cation (Eq. 19), and the reference solvent R is water. Recall that following from Eq. 19, there is no distinction to be ∗

made between ∆Gchg and ∆Gchg and thus, between a Π and a Π∗ parameter. The numerical values of Σ, Σ∗ and Π calculated based on the simulation results are reported in Table 7, 34 ACS Paragon Plus Environment

Page 35 of 74

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 Journal of Physical Chemistry

along with the corresponding values inferred from the Swain A and B parameters via Eqs. 30-32. Considering Figure 8, the correlation between simulation and experiment in terms of both Σ and Σ∗ is seen to be very good (slightly less for the least polar solvents, especially considering Σ∗ ). This correlation is also illustrated explicitly in Suppl. Mat. Figure S.3 (panels a and b). The correlation is not entirely surprising as the coefficients of the affine transformations of Eqs. 30 and 32 have been selected precisely to enforce optimal agreement (excluding the solvents CHL, TOL, BZN and CCL from the fit) along with the constraint Σ = Σ∗ = 1 for water (i.e. two free fitting parameters for each equation, calibrated against 16 observables). The parameters Σ and Σ∗ are interpreted as providing a measure of the solvent asymmetry, entirely compatible with the original dataset of Swain and optimally compatible with the present simulation results. The coefficients of the affine transformation of Eq. 31 have been selected to optimize linear correlation between Π and 1 − ǫ−1 S for the 55 pure solvents considered by Swain. A comparison between the two quantities can be found in the bottom panel of Figure 8, and the correlation is shown explicitly in Suppl. Mat. Figure S.3 (panel e). As also shown in Figure 8, the agreement between simulation and experiment in terms of Π is seen to be acceptable in terms of trend, but rather noisy for the individual solvents. This correlation is illustrated explicitly in Suppl. Mat. Figure S.3 (panel c; see also panel f for the correlation with 1 − ǫS−1 ). The parameter Π is interpreted as providing a measure of the solvent polarity. It is compatible with the original dataset of Swain and presents an approximate correspondance 85 with the permittivity factor 1 − ǫ−1 S entering into the Born model (Eq. 2). However, we do

not attach too much importance to a perfect concordance between the value of Π from the simulations, that derived from the Swain parameters, and the quantity 1 − ǫ−1 S , for two main reasons: (i) the factor 1−ǫ−1 S is only an appropriate approximate descriptor of polarity within the limit of validity of the Born model (and, in particular, considering the special case of a spherical monoatomic ion); (ii) the solvent models employed here 186, 193−197 are probably not able to provide a quantitatively accurate estimate of the magnitude of the charging free energy (in particular because their permittivities only approximately match the experimental ones). For the 55 solvents, the Π parameter is found to be essentially uncorrelated from the Σ or Σ∗ parameters, which suggests that they indeed measure essentially independent properties of a solvent. The corresponding correlations are shown in Suppl. Mat. Figure S.3 (panels g and h), and have correlation coefficients of 0.52 and 0.40, respectively. Finally, it is interesting 35 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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 observe that the Π parameter presents a high correlation with the π ∗ polarity parameter of the Kamlet-Taft scale 28−32 considering the 45 solvents for which both values are available (Table 1). This correlation is illustrated in Suppl. Mat. Figure S.3 (panel i), and has a correlation coefficient of 0.89. Two alternative measures of asymmetry have been retained, Σ and Σ∗ . The parameter Σ∗ is intrinsic and results from an optimization of the affine transformation (Eq. 33) against simulation results in terms of ∆∆G∗chg (Eq. 16). These include the effects of both the linear and the quadratic components. However, in Sections I and III.1, we have already provided arguments against the inclusion of the linear term ∆∆G∗lin in an asymmetry measure to be compared to experiment. For this reason, we also introduced the effective parameter Σ resulting from an optimization of the affine transformation (Eq. 30) against simulation results in terms of the effective quantity ∆∆Gchg only (Eq. 22). As the ∆∆Gchg term typically dominates over the ∆∆G∗lin term (Section III.2), the two parameters are highly correlated. This correlation is shown explicitly in Suppl. Mat. Figure S.3 (panel d). For the reasons given above, we believe that the parameter Σ is the relevant one. Additional support for this suggestion is provided by the observation that the correlation between simulation and experiment is higher in terms of Σ (correlation coefficient of 0.90) than in terms of Σ∗ (correlation coefficient of 0.84). Therefore, further discussion will consider Σ rather than Σ∗ . Considering again Figure 8 along the series of 55 solvents, the following observations can be made. Although the Π parameter is constructed from comparable contributions of A and B (Eq. 31), the Σ parameter is predominantly determined by A (Eq. 30). As a result, the A and Σ profiles along the solvent series look very similar. There are, however, two important differences. First, A was viewed as a relative measure of the acity of a solvent (affinity for anions relative to that of water). In contrast, Σ is viewed as an absolute measure of the asymmetry of a solvent (affinity for anions relative to that for cations) relative to that of water. As a result, the Σ values can be positive or negative. For example a positive value of one-half indicates an asymmetry favoring anions over cations that is about half as strong as that of water, whereas a negative value of minus one-half indicates an asymmetry favoring cations over anions that is about half as strong in magnitude as the opposite asymmetry of water. We will say that the former solvent is acitic, with an acity of one-half, and that the latter solvent is basitic, with a basity of one-half. Thus, the words acity and basity have somewhat different meanings in the revised context of the Σ parameter, i.e. altered relative to the original definitions of Swain (Section I) for the A and B parameters. 36 ACS Paragon Plus Environment

Page 36 of 74

Page 37 of 74

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 Journal of Physical Chemistry

With these updated definitions, the most striking features in Figure 8 considering the parameter Σ are that: (i) the 23 least polar solvents with ǫS ≤ 6.1 are all basitic, with a basity on the order of -0.9 to -0.5 (except CHL); (ii) most of the more polar solvents are also basitic, with a slightly less pronounced basity on the order of -0.5 to -0.3; (iii) all and only the protic solvents (water, alcohols and carboxylic acids) are acitic, with acities that can differ widely in magnitude. Four solvents worth mentioning specifically are CHL (the least basitic of the non-polar solvents with Σ = −0.15), H2 O (the reference solvent water, which sets the step-size of the asymmetry scale with its acity Σ = 1), CF3 COOH (a protic solvent with particularly high acity, Σ = 2.96), and hexamethylphosphoric triamide (a polar non-protic solvent with a particularly high basity, Σ = −1.20). Recall that the latter two solvents were used by Swain as extremes to anchor his scale. 2 The fact that most of the non-polar solvents as well as most of the polar non-protic solvents are found to be moderately basitic is not entirely surprising. In the absence of strong specific ion-solvent electrostatic interactions (accounted for in the simulations by particularly high or/and exposed solvent atomic partial charges), the asymmetry will be largely determined by the electronic polarizability of the solvent molecules. Since the electrons are intrinsically more exposed than the nuclei, this polarizability is expected to contribute to basity, as exemplified by cation-alkane 228 and aromatic cation-π 220−223 interactions (note, however, that anion-π interactions in tangential geometries have also been acknowledged for electron-poor aromatic systems, 221, 224, 225, 229 as well as anion-aromatics hydrogen-bonding 226 in aligned geometries). Since this type of basity cannot be properly accounted for in the context of a non-polarizable force field, it is not surprising that the calculated Σ values are generally smaller in magnitude (i.e. less negative) than the experimentally-inferred ones for these solvents. This underlines that, owing to the very high polarizability of alkanes,230 it is particularly important to take into account electronic polarizability for low-polarity solvents231 where its influence is particularly relevant. 227 If a basity of -0.9 to -0.3 seems to represent the baseline for most organic solvents in the Swain set, acity appears to be the landmark of protic solvents, with a considerably larger variability, from 0 to 2.96. This is certainly related to the dense positive charge of the exposed hydrogen atoms, which can be modulated over a wide range by inductive effects within the solvent molecule.

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

IV

Conclusion

The goal of the present study was to investigate the asymmetric solvation of ions using MD simulations, and to connect the results to the experimental acity-basity scale of Swain. 2,34 To this purpose, the charging free energies of the ions Na+ , K+ , F− , Rb+ and I− , and of their hypothetical oppositely-charged counterparts Na− , K− , F+ , Rb− and I+ , were calculated in a variety of solvents, taking care of correcting the raw estimates for methodology-dependent errors. In a first set of calculations, we considered artificial solvent models termed SAMs and NIMs (Figure 1) presenting either a charge asymmetry or a shape asymmetry, respectively, at the molecular level. The solvation asymmetry was characterized either by means of an intrinsic difference ∆∆G∗chg in charging free energy (cation minus anion), or of a corresponding effective difference ∆∆Gchg . In addition to the effective component, which is quadratic in the ionic charge, the intrinsic quantity also includes a component ∆∆G∗lin , which is linear in the ionic charge and arises from the preferential orientation of the solvent around the uncharged ion-sized cavity. For these simple solvent models, it was found that the components ∆∆Gchg and ∆∆G∗lin are systematically of opposite signs, an observation that is easily understood. A solvent presenting a more exposed positive charge will orient this charge towards an anion and away from a cation, thereby preferentially solvating the anion in terms of the effective term (∆∆Gchg > 0). However, in contact with a cavity, it will preferentially orient this exposed charge towards the solvent (more polar region) than towards the cavity (non-polar region), thereby preferentially solvating the cation in terms of the linear term (∆∆G∗lin < 0). Because the linear term originates from a potential that is, at least to a good approximation, insensitive to the size and shape of the ionic cavity (in line with the present results considering different ions) as well as homogeneous over the volume of this cavity, 60, 64, 77 it is reasonable to assume that it cannot influence any experimental observable corresponding to a process that preserves electroneutrality. In practice, the experimental data considered when establishing solvent-property scales such as the Swain scale entirely excludes processes involving the absorption or emission of charged subatomic particles, 174−176 or where single ions are transferred between different media. 1,130,177,178 As a result, we suggest that the linear term is in fact irrelevant, and that any connection between the MD simulation results and the experimentally-inferred asymmetric solvation properties of a solvent should rely on effective

38 ACS Paragon Plus Environment

Page 38 of 74

Page 39 of 74

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 Journal of Physical Chemistry

rather than intrinsic quantities. Note that the linear term is expected to be related to the air-liquid interfacial potential of the solvent at a planar interface, the analogy being only limited by the different interface curvatures considered. 161 Unfortunately, the experimental estimates for water are not consistent even in their signs, and there are only very few estimates for other solvents. 1 In a second set of calculations, we considered the Na− /Na+ ion pair in 16 physical solvent models (Table 2) including water as well as non-aqueous solvents. Theoretical estimates for the relative acities and basities were obtained by comparing the charging free energies of the ions in the given solvent to the corresponding values in water. This was done considering both intrinsic and effective charging free energies, although the latter values are considered to be more appropriate (see above). The correlation with the experimentally-inferred values of Swain for the relative acity A and the relative basity B was found to be remarkable for the acity, but very poor for the basity. This led us to question the critical (scale-setting) conditions employed in the derivation of the Swain scale, as these conditions belong to the interpretation of the fit against experimental data without affecting in any way the quality of this fit. Based on the simulation results, the original interpretation of Swain was modified so as to generate an alternative asymmetry-polarity scale by means of affine transformations of the A and B parameters. The resulting (effective) Σ, Π-scale has the same information content as the original scale, but is made to be optimally compatible with the simulation results. In this revised scale, the parameter Σ characterizes the asymmetry in an absolute sense. It can be either negative (basitic solvent) or positive (acitic solvent), the asymmetry magnitude of one being defined by the acity magnitude of water. The parameter Π characterizes independently the polarity of the solvent, i.e. its ability to solvate overall neutral sets of charged species (for which asymmetric effects cancel out). In retrospect, the need for such a reinterpretation of the Swain scale could have been anticipated. Experimental measurements are bound by the electroneutrality constraint (see above), and accessing acities and basities independently is as impossible as trying to define an intrinsic scale of single-ion solvation free energies without invoking an extra-thermodynamic assumption. Considering the simple case of an ion-pair X − /X + , only the sum of the two solvation free energies is experimentally accessible. It follows that the corresponding difference can only be accessed within an unknown constant. By converting the A, B-scale into a Σ, Π-scale, which precisely characterizes the difference and the sum, it is ensured that the 39 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

missing information is only included in the parameter Σ (the parameter Π being formally unambiguous), rather than arbitrarily distributed across A, B-pairs. The second component of the procedure is to set a zero for the Σ scale, which requires an extra-thermodynamic assumption, i.e. a model. In the present case, this model information is supplied by the atomistic MD simulations. After conversion of the A, B-pairs of Swain into alternative Σ, Πpairs, the Σ values are found to be captured remarkably well by the simulations, at least for sufficiently polar solvents when relying on a non-polarizable force field. In comparison, the Π values are more difficult to capture accurately, probably due to force-field inaccuracies and to the very large magnitude of ionic solvation free energies. The omission of the linear term in the comparison with experiment leads to two interesting observations concerning the nature of the extra-thermodynamic assumption anchoring the Σ scale. First, because this linear term corresponds in essence to the effect of the airliquid interfacial potential of the solvent (see above), it suggests that what we called effective quantities here might actually be real quantities 1 (see also Appendix A). In other words, the Σ value may in some way be connected to the real solvation free energy of the proton in the given solvent. Second, this same linear term is expected to be the principal source of asymmetry when comparing the large ions tetraphenylarsonium and tetraphenylborate within the TATB assumption. 43, 53, 55, 76, 94−96 As a result, we suggest that the use of this assumption to define single-ion solvation free energies in different solvents may as well lead to real values rather than intrinsic ones, i.e. that this assumption may in essence be consistent with the Σ scale. In contrast to the original A and B values, which only provided an asymmetry information relative to water, the derived Σ parameter provides via its sign an absolute information on the acity or basity of a solvent. Considering the 55 solvents in the Swain set, it is observed that a moderate basity (Σ between -0.9 and -0.3) represents the baseline for most solvents, while a highly variable acity (Σ between 0.0 and about 3.0) represents a landmark of protic solvents. We interpret these findings by observing that the electrons are more exposed relative to the nuclei, which renders a solvent basitic by default due to electronic polarization. This situation is exemplified by cation-alkane 228 and aromatic cation-π 220−223 interactions. In contrast, a protic hydrogen atom represents a dense and surface-exposed positive charge, and thus induces an acity, the magnitude of which may be heavily influenced by inductive effects within the solvent molecule. This situation is exemplified here by the solvents chloroform and water, the alcohols, and the carboxylic acids, which are all expected to present 40 ACS Paragon Plus Environment

Page 40 of 74

Page 41 of 74

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 Journal of Physical Chemistry

hydrogen-bonding interactions with anions. Clearly, this study only represents a preliminary step towards the connection between the results of atomistic simulations and the experimentally available solvent-property scales. In particular, work is in progress to consider a wider range of solvents as well as polyatomic or/and polyvalent ions. This should lead to a better understanding of the extrathermodynamic assumption anchoring the Σ-scale, as well as its relationship with real solvation free energies and the TATB hypothesis.

Appendix A The goal of this appendix is to restate briefly our main conclusions from previous work 1 , 159

concerning the appropriate definition of the absolute electric potential within a cavity

surrounded by a liquid. This is motivated by the fact that a number of recent (and otherwise excellent) studies 64, 77, 154, 171−173 appear to still overlook or handle incorrectly this issue (see, however, Refs. 96, 165, 232−237 ), and to adopt a convention that is not acceptable in practice, because it introduces an arbitrary constant (the quadrupole-moment trace of a classical solvent model) into the calculation results (air-liquid interfacial potentials, potentials within cavities, and single-ion solvation free energies). Since all the mathematical details as well as additional literature references can be found in Refs. 1, 159 (see also Refs. 84, 160, 161 ), the present discussion is meant to be qualitative and intuitive rather than formal. As illustrated in Figure 9, there are three possible conventions for setting the zero of the electric potential when attempting to define the absolute value of this potential within a cavity surrounded by a liquid, namely : (i) the zero is set outside a macroscopic liquid sample (R-convention); (ii) the zero is set as an average over the bulk liquid (far away from the cavity), averaging over both the interior and the exterior of the solvent molecules (Pconvention); (iii) the zero is set as an average over the bulk liquid, but averaging only over the exterior of the solvent molecules (M-convention). Most experiments only probe processes involving overall neutral sets of ions. 1, 11 When comparing theoretical results to those of such experiments, any of the above conventions is equally acceptable. In very special types of experiments (Voltaic cell measurements1,130,177,178 ), it is possible to obtain (indirectly) information on the transfer free energy of individual ions from the ideal gas phase to the bulk of a liquid. This corresponds to so-called real single-ion solvation free energies. When comparing simulation results to those of such experiments, the

41 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

appropriate convention is the R-convention. The problem with real solvation free energies is that they intermingle two distinct physical effects, 238−240 namely the chemical solvation of the ion by the surrounding solvent molecules and the work required for the ion to penetrate the solvent surface. For this reason, for more than one century, chemists have devoted an incredible amount of effort in trying to access so-called intrinsic single-ion solvation free energies, which are exempt of the interface-crossing component and only account for the effect of ion-solvent interactions. One should realize, however, that: 1, 11 (i) this partitioning is only realizable by adding information from a model which, in this specific context, is often called an extra-thermodynamic assumption; (ii) the resulting intrinsic quantities may be very useful for reasoning within this specific model (and, possibly, across models involving compatible assumptions), but a specific partitioning cannot have a direct bearing on any experimental observable. As long as the extra-thermodynamic assumptions invoked were rooted in continuum electrostatics (Born model 85 and arguments concerning crystallographic ionic radii 48, 86−93 ), there was no reason to distinguish between a P- and a M-convention. This is because the continuum representation of a dielectric medium involves infinitesimal volume elements that carry a dipole but are exempt of internal structure. However, with the advent of atomistic simulations, the granularity of the solvent molecules has been made explicit, and the distinction has become relevant. As discussed elsewhere,1,159 if one eliminates finite-size or/and finite-cutoff artifacts in the configurations generated by atomistic simulation (by taking both parameters to arbitrarily large values), the potential within a cavity surrounded by the liquid will correspond to the M-convention if it is calculated considering Coulombic interactions truncated on a molecular basis at a very large distance from the ion (CM scheme), but to the P-convention if it is calculated using a lattice-sum or reaction-field method instead (LS or BM scheme, the latter in the limit of high reaction-field permittivities). In practical situations, the difference in terms of the calculated potential within a cavity (and of corresponding intrinsic single-ion solvation free energies) can be extremely large. We call this difference the exclusion potential ξS′ of the solvent model. 1 It is defined as the potential in a pure liquid when averaging over the outside of the solvent molecules minus this potential when averaging over the entire sample including the interior of the molecules. Considering Figure 9, for the potential within a cavity, one has thus φP = φM + ξS′ . For a rigid solvent model with a single van der Waals interaction site, it is easily shown1,159 that ξS′ is related to the quadrupole-moment trace γS′ of the solvent model 42 ACS Paragon Plus Environment

Page 42 of 74

Page 43 of 74

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 Journal of Physical Chemistry

relative to this site. Thus, for example, one finds ξS′ values of 0.17, 0.82, 0.85, 0.77, 0.93 and 0.21 V for the ST2, 241 SPC, 186 SPC/E, 242 TIP3P, 243 TIP4P 243 and TIP5P 244 water models, respectively. However, as has been noted earlier, 1, 159, 232 the quadrupole-moment trace of a classical solvent model is an entirely arbitrary quantity that plays no role in the fitting of the force-field parameters against experimental data when this model is developed. For example, by playing with overall-neutral concentric spherical charge layers very close to the oxygen atom (thus well within the van der Waals envelope of the molecule and essentially invisible to the neighboring solvent molecules), one might easily engineer SPC-like water models with exclusion potentials of -1000 V or +1000 V (or any other value), without affecting the reproduction of any measurable thermodynamic or structural property of the liquid. It is also to be noted that a quantum-mechanical representation of the water molecule is expected to have a negative ξS′ instead, 159,172,234,237 because the negative charge (electrons) is, in this case, more exposed than the positive charge (nuclei). Because ξS′ is an essentially meaningless (arbitrary) quantity, it follows from the relation φP = φM + ξS′ that either the M- or the P-convention (or both) is not sustainable, even in the sense of generating intrinsic quantities amenable to theoretical reasoning (if not experimental measurement). Over the last three decades, there has been a very hot debate in the context of the calculation of intrinsic single-ion solvation free energies based on atomistic simulations between the partisans of the M-convention 84,84,159−161, 166−170,245 and those of the P-convention. 49,60,61, 83, 218, 246 −251

A similar debate has taken place in the context of the calculation of the surface potential

of liquids, with the same opposition between the partisans of the M-convention 159, 172, 252−259 and those of the P-convention. 232, 233, 259−268 Building on key considerations by others, 166, 167 we believe to have provided extremely convincing arguments 1, 159 in favor of the M-convention. The main argument concerns the limiting case where a solvent molecule is a perfect isotropic quadrupole (e.g. overall-neutral concentric spherical charge layers or rare-gas atom viewed as a nucleus surrounded by a spherically symmetric electron cloud) or made to rotate so fast (high rotational temperature limit 167 ) that it becomes an isotropic quadrupole on average. Clearly, such a solvent is expected to be non-solvating. In agreement with this expectation, the potential within a cavity using the M-convention will be ΦM = 0. But using the P-convention, one will find ΦP = ξS′ instead i.e., as discussed above, an entirely arbitrary quantity. Similarly, for water at finite temperature, one typically has ΦM < 0 (see below) but, with most of the common classical water models, ΦP > 0. Upon increasing the temperature, both ΦM and ΦP increase. 159 The magnitude of ΦM thus decreases (towards zero), an 43 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

expected consequence of increased thermal motion, i.e. of the increased orientational averaging of the water molecules at the cavity surface. In contrast, the magnitude of ΦP increases (towards ξS′ ), a feature that incorrectly suggests an increased orientational structure of the solvent at the cavity surface (see e.g. Figure 4 in Ref. 173 ; note that the potential displayed there matches −ΦP , as it is measured in the air-to-water direction). For water at finite temperature, the M-based potential within a cavity is found to be negative. 159 This is consistent with the negative ∆∆G∗lin asymmetry term reported in this article, i.e. this term favors cation over anion solvation. This finding is compatible with water molecules preferentially pointing their hydrogen atoms away from the cavity. It is also intuitively reasonable considering that the more polar end of the molecule is expected to preferentially interact with bulk water (polar medium) rather than to be exposed to the empty cavity (non-polar medium). Of course, in practice, this linear term is overcompensated by the quadratic term ∆∆Gchg , which evidences the opposite trend, and water is hence overall acitic and not basitic. In a number of recent studies, however, the P-based convention was used instead, typically linked with the direct application of LS electrostatics. Since the exclusion potentials of the commonly employed atomistic solvent models are positive and relatively large in magnitude (see above), it has been suggested in these studies that the potential within a cavity surrounded by water is positive, 49, 60, 61, 64, 77, 83, 154, 172, 218, 246−251 leading to an opposite trend for ∆∆G∗lin , now apparently acting in the same direction as ∆∆Gchg . Based on this incorrect negative value, some authors have implicitly assumed that water preferentially (slightly) points its dipole towards a cavity (see e.g. Figure 4 in Ref. 64 ). This is structurally incorrect, 237 as exemplified by the orientational-correlation function cIW shown in Figure 5c for the SPC water model 186 around a sodium-sized cavity. Other authors have relied on radialdistribution functions, and observed that the first ion-hydrogen peak is extends closer to the cavity compared to the ion-oxygen peak,77,82 thus apparently justifying a positive potential in the cavity. This observation is correct, as also exemplified by the radial-distribution functions gIO and gIH shown in Figure 5c, but its interpretation is not. This is because one should not consider the gIH function in an absolute sense, but in comparison to a corresponding function generated by isotropically rotating water molecules with their oxygen atoms distributed according to the gIO function. Such an alternative function g˜IH has been calculated and is also displayed in the figure. Clearly, it expands even closer to the cavity compared to the true distribution gIH . This shows that the density of hydrogen atoms close to the cavity is 44 ACS Paragon Plus Environment

Page 44 of 74

Page 45 of 74

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 Journal of Physical Chemistry

actually depleted relative to a bulk-like isotropic situation, in agreement with the positive cIO and with the the negative potential of -0.40 V in the cavity. Generally speaking, at an air-water interface, the majority of the hydrogen atoms will orient parallel to the surface, some will orient towards the bulk, and only very few will orient away from the bulk 261,262,269− 272

(note that the same may not hold for water at interfaces with all hydrophobic media273−276 ).

Note, finally, that in Refs. 158, 237 , the entire issue is circumvented by applying a R-convention instead. This is of course acceptable, but will lead to real electric potentials and single-ion solvation free energies instead of intrinsic ones. As a last comment, one might argue that since intrinsic quantities only have a theoretical (model-bound) value, the choices of a M- or a P-convention for the absolute potential definition merely correspond to alternative and equally acceptable models. This is fine if one accepts that the P-convention results encompass an arbitrary offset which renders them dependent on a specific choice of solvent model (rather than on the sole chemical nature of the solvent considered), and formally incompatible (even in an approximate sense) with all other models such as quantum mechanics, atomistic simulations with M-convention and continuum electrostatics, as well as with most of the extra-thermodynamic assumptions employed to date.

Supporting Information Supplementary material is provided in a separate document. This document includes detailed information concerning the calculations of the dielectric permittivities, exclusion potentials and charging free energies, as well as additional figures illustrating the correlation between different solvent parameters and the solvation structure around ions and uncharged cavities.

Acknowledgements The authors would like to thank Bruno Horta for making the 2016H66 force-field files available prior to publication, Pascal Merz for providing topologies and equilibrated coordinates for the corresponding solvents, and Martin Zacharias for his generosity in terms of computational resources. Last but not least, PHH would like to express his gratitude to Andy McCammon, for the three extraordinary post-doc years spent in his group at UCSD.

45 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

References [1] H¨ unenberger, P.H.; Reif, M.M. Single-ion solvation: Experimental and theoretical approaches to elusive thermodynamic quantities. Edition 1, ISBN: 978-1-84755-187-0.; Royal Society of Chemistry (Theoretical and Computational Chemistry Series), London, UK, 2011. [2] Swain, C.G.; Swain, M.S.; Powell, A.L.; Alunni, S. Solvent effects on chemical reactivity. Evaluation of anionand cation-solvation components. J. Am. Chem. Soc. 1983, 105, 502-513. [3] Reichardt, C. Solvatochromic dyes as solvent polarity indicators. Chem. Rev. 1994, 94, 2319-2358. [4] Reichardt, C.; Welton, T. Solvents and solvent effects in organic chemistry. Edition 4.; Wiley-VCH, Weinheim, Germany, 2011. [5] Katritzky, A.R.; Fara, D.C.; Yang, H.; T¨ amm, K.; Tamm, T.; Karelson, M. Quantitative measures of solvent polarity. Chem. Rev. 2004, 104, 175-198. [6] Nightingale Jr., E.R. Phenomenological theory of ion solvation. Effective radii of hydrated ions. J. Phys. Chem. 1959, 63, 1381-1387. [7] Shannon, R.D.; Prewitt, C.T. Effective ionic radii in oxides and fluorides. Acta Crystallogr. B 1969, 25, 925-946. [8] Johnson, O. Ionic radii for spherical potential ions. I. Inorg. Chem. 1973, 12, 780-785. [9] Holbrook, J.B.; Khaled, F.M.; Smith, B.C. Soft-sphere ionic-radii for Group 1 and Group 2 metal halides and ammonium halides. Dalton Trans. 1978, 12, 1631-1634. [10] Murakami, W.; Yamamoto, M.; Eda, K.; Osakai, T. A non-Bornian analysis of the Gibbs energy of hydration for organic ions. RSC Adv. 2014, 4, 27634-27641. [11] 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. [12] Matyushov, D.V.; Schmid, R.; Ladanyi, B.M. A Thermodynamic Analysis of the π ∗ and ET (30) Polarity Scale. J. Phys. Chem. B 1997, 101, 1035-1050. ¨ [13] Ohrn, A.; Karlstr¨ om, G. Many-body polarization, a cause of asymmetric solvation of ions and quadrupoles. J. Chem. Theory Comput. 2007, 3, 1993-2001. [14] Buyukdagli, S.; Ala-Nissila, T. Alteration of gas phase ion polarizabilities upon hydration in high dielectric liquids. J. Chem. Phys. 2013, 139, 044907/1-044907/11. [15] Taft, R.W.; Abboud, J.-L. M.; Mortimer, J.K. Linear solvation energy relationships. 28. An analysis of Swain’s solvent “acity” and “basity” scales. J. Org. Chem. 1984, 49, 2001-2005. [16] Joerg, S.; Drago, R.S.; Adams, J. Donor-acceptor and polarity parameters for hydrogen bonding. J. Chem. Soc. Perkin Trans. 1997, 2, 2431-2438. [17] Pytela, O. Empirical approach to description of solvent effect on processes in solutions: A review. Coll. Czech. Chem. Commun. 1988, 53, 1333-1423. [18] Catal´ an, J; Toward a generalized treatment of the solvent effect based on four empirical scales: Dipolarity (SdP, a new scale), polarizability (SP), acidity (SA), and basicity (SB) of the medium. J. Phys. Chem. B 2009, 113, 5951-5960. [19] Cabot, R.; Hunter, C.A. Molecular probes of solvation phenomena. Chem. Soc. Rev. 2012, 41, 3485-3492. ¨ [20] Dimroth, K.; Reichardt, C.; Siepmann, T.; Bohlmann, F. Uber Pyridinium-N-Phenol-Betaine und ihre Verwendung zur Charakterisierung der Polarit¨ at von L¨ osungsmitteln. Liebigs Ann. Chem. 1963, 661, 1-37. [21] Dimroth, K.; Reichardt, C. Erweiterung der L¨ osungsmittelpolarit¨ atsskala durch Verwendung Alkylsubstituierter Pyridinium-N-phenol-betaine. Liebigs Ann. Chem. 1969, 727, 93-105. [22] Mayer, U.; Gutmann, V.; Gerger, W. The acceptor number. A quantitative empirical parameter for the electrophilic properties of solvents. Monatsh. Chem. 1975, 106, 1235-1257. [23] Mayer, U. Eine sempiempirische Gleichung zur Beschreibung des L¨ osungsmitteleinflusses auf Statik und Kinetik chemischer Reaktionen. Teil I. Monatsh. Chem. 1978, 109, 421-433. [24] Mayer, U. Eine sempiempirische Gleichung zur Beschreibung des L¨ osungsmitteleinflusses auf Statik und Kinetik chemischer Reaktionen. Teil II. Monatsh. Chem. 1978, 109, 755-790. [25] Gutmann, V. Solvent effects on the reactivities of organometallic compounds. Coord. Chem. Rev. 1976, 18, 225-255.

46 ACS Paragon Plus Environment

Page 46 of 74

Page 47 of 74

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 Journal of Physical Chemistry

[26] Gutmann, V. Empirical parameters for donor and acceptor properties of solvents. Electrochim. Acta 1976, 21, 661-670. [27] Parker, A.J.; Mayer, U.; Schmid, R.; Gutmann, V. Correlation of solvent effects on rates of solvolysis and SN 2 reactions. J. Org. Chem. 1978, 43, 1843-1854. [28] Taft, R.W.; Pienta, N.J.; Kamlet, M.J.; Arnett, E.M. Linear solvation energy relationships. 7. Correlations between the solvent-donicity and acceptor-number scales and the solvatochromic parameters π ∗ , α, and β. J. Org. Chem. 1981, 46, 661-667. [29] Kamlet, M.; Taft, R.W. Linear solvation energy relationships. Part 1. Solvent polarity-polarizability effects on infrared spectra. J. Chem. Soc. Perkin Trans. 2 1979, 1979, 337-341. [30] Kamlet, M.J.; Abboud, J.-L.M.; Abraham, M.H.; Taft, R.W. Linear solvation energy relationships. 23. A comprehensive collection of the solvatochromic parameters, π ∗ , α, and β, and some methods for simplifying the generalized solvatochromic equation. J. Org. Chem. 1983, 48, 2877-2887. [31] Taft, R.W.; Abboud, J.-L.M.; Kamlet, M.J.; Abraham, M.H. Linear solvation energy relations. J. Solut. Chem. 1985, 11, 153-186. [32] Marcus, Y. Linear solvation energy relationships. Standard molar Gibbs free energies and enthalpies of transfer of ions from water into nonaqueous solvents. J. Phys. Chem. 1988, 92, 3613-3622. [33] Abe, T. Improvements of the empirical π ∗ polarity scale. Bull. Chem. Soc. Jpn. 1990, 63, 2328-2338. [34] Swain, C.G. Substituent and solvent effects on chemical reactivity. J. Org. Chem. 1984, 49, 2005-2010. [35] Katritzky, A.R.; Tamm, T.; Wang, Y.; Karelson, M. A unified treatment of solvent properties. J. Chem. Inf. Comput. Sci. 1999, 39, 692-698. [36] Katritzky, A.R.; Fara, D.C.; Kuanar, M.; Hur, E.; Karelson, M. The classification of solvents by combining classical QSPR methodology with principal component analysis. J. Phys. Chem. A 2005, 109, 10323-10341. [37] Mancini, P.M.; Adam, C.G.; Fortunato, G.G.; Votteroa, L.R. A comparison of nonspecific solvent scales. Degree of agreement of microscopic polarity values obtained by different measurement methods. Arkivoc 2007, 2007(16), 266-280. [38] Wrona, P.K. On the correlations between empirical Lewis acid-base solvent parameters and the thermodynamic parameters of ion solvation. J. Electroanal. Chem. 1980, 108, 153-167. [39] Wrona, P.K.; Krygowski, T.M.; Galus, Z. Correlatoins between empirical Lewis acid-base solvent parameters and the thermodynamic parameters of ion solvation. Part II. Acidity parameters of cations and basicity parameters of anions. J. Phys. Org. Chem. 1991, 4, 439-448. [40] Katritzky, A.R.; Tamm, T.; Wang, Y.; Sild, S.; Karelson, M. QSPR treatment of solvent scales. J. Chem. Inf. Comput. Sci. 1999, 39, 684-691. [41] Wicaksono, D.S.; Mhamdi, A.; Marquardt, W. Computer-aided screening of solvents for optimal reaction rates. Chem. Eng. Sci. 2014, 115, 167-176. [42] Zhou, T.; Lyu, Z.; Qi, Z.; Sundmacher, K. Robust design of optimal solvents for chemical reactions. A combined experimental and computational strategy. Chem. Engin. Sci. 2015, 137, 613-625. [43] 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. [44] Jolicoeur, C.; The, N.D.; Cabana, A. Near infrared spectra of water in aqueous solutions of organic salts. A solvation study of Bu4 NBr, Φ4 AsCl, and NaBΦ4 . Can. J. Chem. 1971, 49, 2008-2013. [45] Stangret, J.; Kamie´ nska-Piotrowicz, E. Effect of tetraphenylphosphonium and tetraphenylborate ions on the water structure in aqueous solutions; FTIR studies of HDO spectra. J. Chem. Soc. Farad. Trans 1997, 93, 3463-3466. [46] Stangret, J.; Gampe, T. Ionic hydration behavior derived from infrared spectra in HDO. J. Phys. Chem. A 2002, 106, 5393-5402. [47] Tielrooij, K.J.; van der Post, S.T.; Hunger, J.; Bonn, M.; Bakker, H.J. Anisotropic water reorientation around ions. J. Phys. Chem. B 2011, 115, 12638-12647. [48] Latimer, W.M.; Pitzer, K.S.; Slansky, C.M. The free energy of hydration of gaseous ions, and the absolute potential of the normal calomel electrode. J. Chem. Phys. 1939, 7, 108-111. [49] Lynden-Bell, R.M.; Rasaiah, J.C. From hydrophobic to hydrophilic behaviour: A simulation study of solvation entropy and free energy of simple solutes. J. Chem. Phys. 1997, 107, 1981-1991. [50] Koneshan, S.; Rasaiah, J.C.; Lynden-Bell, R.M.; Lee, S.H. Solvent structure, dynamics, and ion mobility in aqueous solutions at 25◦ C. J. Phys. Chem. B 1998, 102, 4193-4204.

47 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

[51] Roux, B.; Yu, H.-A.; Karplus, M. Molecular basis for the Born model of ion solvation. J. Phys. Chem. 1990, 94, 4683-4688. [52] Hummer, G.; Pratt, L.R.; Garcia, A.E. Free energy of ionic hydration. J. Phys. Chem. 1996, 100, 1206-1215. − [53] 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.

[54] Ashbaugh, H.S. Convergence of molecular and macroscopic continuum descriptions of ion hydration. J. Phys. Chem. B 2000, 104, 7235-7238. [55] Schurhammer, R.; Wipff, G. About the TATB assumption: effect of charge reversal on transfer of large spherical ions from aqueous to non-aqueous solvents and on their interfacial behaviour. J. Mol. Struct. (Theochem) 2000, 500, 139-155. − [56] 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.

[57] Schurhammer, R.; Engler, E.; Wipff, G. Hydrophobic ions in TIP5P water and at water-chloroform interfaces: The effect of sign inversion investigated by MD and FEP simulations. J. Phys. Chem. B 2001, 105, 1070010708. [58] Lynden-Bell, R.M.; Rasaiah, J.C.; Noworyta, J.P. Using simulation to study solvation in water. Pure Appl. Chem. 2001, 73, 1721-1731. [59] Rasaiah, J.C.; Lynden-Bell, R.M. Computer simulation studies of the structure and dynamics of ions and non-polar solutes in water. Phil. Trans. R. Soc. Lond. A 2001, 359, 1545-1574. [60] Rajamani, S.; Ghosh, T.; Garde, S. Size dependent ion hydration, its asymmetry, and convergence to macroscopic behavior. J. Chem. Phys. 2004, 120, 4457-4466. [61] Grossfield, A. Dependence of ion hydration on the sign of the ion’s charge. J. Chem. Phys. 2005, 122, 024506/1024506/10. [62] Mobley, D.L.; Barber II, A.E.; Fennell, C.J.; Dill, K.A. Charge asymmetries in hydration of polar solutes. J. Phys. Chem. B 2008, 112, 2405-2414. [63] Fennell, C.J.; Dill, K.A. Physical modeling of aqueous solvation. J. Stat. Phys. 2011, 145, 209-226. [64] Bardhan, J.P.; Jungwirth, P.; Makowski, L. Affine-response model of molecular solvation of ions: Accurate predictions of asymmetric charging free energies. J. Chem. Phys. 2012, 137, 124101/1-124101/6. [65] Mukhopadhyay, A.; Tolokh, I.S.; Onufriev, A.V. Accurate evaluation of charge asymmetry in aqueous solvation. J. Phys. Chem. B 2015, 119, 6092-6100. [66] Purisima, E.O.; Sulea, T. Restoring charge asymmetry in continuum electrostatics calculations of hydration free energies. J. Phys. Chem. B 2009, 113, 8206-8209. [67] Bardhan, J.P.; Knepley, M.G. Communication: Modeling charge-sign asymmetric solvation free energies with nonlinear boundary conditions. J. Chem. Phys. 2014, 141, 131103/1-131103/6. [68] Ivanishko, I.S.; Borovkov, V.I. Comparison of the mobilities of negative and positive ions in nonpolar solutions. J. Phys. Chem. B 2010, 114, 9812-9819. [69] Gehring, T.; Diffusion of nanoparticles at an air/water interface is not invariant under a reversal of the particle charge. J. Phys. Chem. C 2011, 115, 23677-23681. [70] Chorny, I.; Dill, K.A.; Jacobson, M.P. Surfaces affect ion pairing. J. Phys. Chem. B 2005, 109, 24056-24060. [71] Fennell, C.J.; Bizjak, A.; Vlachy, V.; Dill, K.A. Ion pairing in molecular simulations of aqueous alkali halide solutions. J. Phys. Chem. B 2009, 113, 6782-6791. [72] Perera, L.; Berkowitz, M.L. Structure and dynamics of Cl− (H2 O)20 clusters: The effect of the polarizability and the charge of the ion. J. Chem. Phys. 1992, 96, 8288-8294. [73] Torrie, G.M.; Patey, G.N. Molecular solvent model for an electrical double layer: Asymmetric solvent effects. J. Phys. Chem. 1993, 97, 12909-12918. [74] Kinoshita, M.; Hirata, F. Application of the reference interaction site model theory to analysis on surfaceinduced structure of water. J. Chem. Phys. 1996, 104, 8807-8815. [75] Jungwirth, P.; Tobias, D.J. Ions at the air/water interface. J. Phys. Chem. B 2002, 106, 6361-6373. [76] Scheu, R.; Rankin, B.M.; Chen, Y.; Jena, K.C.; Ben-Amotz, D.; Roke, S. Charge asymmetry at aqueous hydrophobic interfaces and hydration shells. Angew. Chem. Int. Ed. 2014, 53, 9560-9563.

48 ACS Paragon Plus Environment

Page 48 of 74

Page 49 of 74

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 Journal of Physical Chemistry

[77] Cerutti, D.S.; Baker, N.A.; McCammon, J.A. Solvent reaction field potential inside an uncharged globular protein: A bridge between implicit and explicit solvent models. J. Chem. Phys. 2007, 127, 155101/1-155101/12. [78] Gunner, M.R.; Saleh, M.A.; Cross, E.; ud-Doula, A.; Wise, M. Backbone dipoles generate positive potentials in all proteins: origins and implications of the effect. Biophys. J. 2000, 78, 1126-1144. [79] Kim, J.; Mao, J.; Gunner, M.R. Are acidic and basic groups in buried proteins predicted to be ionized. J. Mol. Biol. 2005, 348, 1283-1298. [80] Vanin, A.A.; Brodskaya, E.N. Computer simulation of the surface layer of an ionic micelle with explicit allowance for the contribution of water. Colloid J. 2015, 77, 409-417. [81] Brodskaya, E.N.; Vanin, A.A. Effect of water on the local electric potential of simulated ionic micelles. J. Chem. Phys. 2015, 143, 044707/1-044707/9. [82] Garde, S.; Hummer, G.; Paulaitis, M.E. Free energy of hydration of a molecular ionic solute: Tetramethylammonium ion. J. Chem. Phys. 1998, 108, 1552-1561. [83] Herce, D.H.; Darden, T.; Sagui, C. Calculation of ionic charging free energies in simulation systems with atomic charges, dipoles, and quadrupoles. J. Chem. Phys. 2003, 119, 7621-7632. [84] 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. [85] Born, M. Volumen und Hydrationsw¨ arme der Ionen. Z. Phys. 1920, 1, 45-48. [86] Lange, E.; Mishchenko, K.P. Zur Thermodynamik der Ionensolvatation. Z. Phys. Chem. A 1930, 149, 1-41. [87] Baughan, E.C. The heat of hydration of the proton. J. Chem. Soc. 1940, (Oct.), 1403-1403. [88] Coulter, L.V. The thermochemistry of the alkali and alkaline earth metals and halides in liquid ammonia at -33◦ . J. Phys. Chem. 1953, 57, 553-558. [89] Blandamer, M.J.; Symons, M.C. Significance of new values for ionic radii to solvation phenomena in aqueous solution. J. Phys. Chem. 1963, 67, 1304-1306. [90] Franks, F.; Reid, D.S. Ionic solvation entropies in mixed aqueous solvents. J. Phys. Chem. 1969, 73, 3152-3154. [91] Marcus, Y. The thermodynamics of solvation of ions. Part 2. The enthalpy of hydration at 298.15 K. J. Chem. Soc. Farad. Trans. I 1987, 83, 339-349. [92] Marcus, Y. Thermodynamics of ion hydration and its interpretation in terms of a common model. Pure Appl. Chem. 1987, 59, 1093-1101. [93] 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 1987, 83, 2985-2992. [94] Kim, J.I. Preferential solvation of single ions. A critical study of the Ph4AsPh4B assumption for single ion thermodynamics in amphiprotic and dipolar-aprotic solvents. J. Phys. Chem. 1978, 82, 191-199. [95] Fawcett, W.R. Acidity and basicity scales for polar solvents. J. Phys. Chem. 1993, 97, 9540-9546. [96] Pollard, T.P.; Beck, T.L. The thermodynamics of proton hydration and the electrochemical surface potential of water. J. Chem. Phys. 2014, 141, 18C512/1-18C512/10. [97] Hirata, F.; Redfern, P.; Levy, R.M. Viewing the Born model for ion hydration through a microscope. Int. J. Quant. Chem. 1988, 15, 179-190. [98] Chan, D.Y.C.; Mitchell, D.J.; Ninham, B.W. A model of solvent structure around ions. J. Chem. Phys. 1979, 70, 2946-2957. [99] Blum, L.; Fawcett, W.R. Application of the mean spherical approximation to describe the Gibbs solvation energies of monovalent monoatomic ions in polar solvents. J. Phys. Chem. 1992, 96, 408-414. [100] Fawcett, W.R.; Blum, L. Application of the mean spherical approximation to the estimation of single ion thermodynamic quantities of solvation for monoatomic monovalent ions in aqueous solutions. J. Electroanal. Chem. 1992, 328, 333-340. [101] Babu, C.S.; Lim, C. Theory of ionic hydration: Insights from molecular dynamics simulations and experiment. J. Phys. Chem. B 1999, 103, 7958-7968. [102] Babu, C.S.; Lim, C. A new interpretation of the effective Born radius from simulations and experiment. Chem. Phys. Lett. 1999, 310, 225-228.

49 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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] Babu, C.S.; Lim, C. Incorporating nonlinear solvent response in continuum dielectric models using a two-sphere description of the Born radius. J. Phys. Chem. A 2001, 105, 5030-5036. [104] Sandberg, L.; Edholm, O. Nonlinear response effects in continuum models of the hydration of ions. J. Chem. Phys. 2002, 116, 2936-2944. [105] Laidler, K.J.; Muirhead-Gould, J.S. Continuous dielectric model for hydration of monatomic ions. Trans. Farad. Soc. 1967, 63, 953-957. [106] Goldman, S.; Bates, R.G. Calculation of thermodynamic functions for ionic hydration. J. Am. Chem. Soc. 1972, 94, 1476-1784. [107] Ehrenson, S. Continuum radial dielectric functions for ion and dipole solution systems. J. Comp. Chem. 1989, 10, 77-93. [108] Kumar, A. Surface effects in the Born solvation model. Bull. Chem. Soc. Jpn. 1994, 67, 3150-3152. [109] Basilevsky, M.V.; Parsons, D.F. An advanced continuum medium model for treating solvation effects: Nonlocal electrostatics with a cavity. J. Chem. Phys. 1996, 105, 3734-3746. [110] Hyun, J.-K.; Babu, C.S.; Ichiye, T. Apparent local dielectric response around ions in water: A method for its determination and its applications. J. Phys. Chem. 1995, 99, 5187-5195. [111] Hyun, J.-K.; Ichiye, T. Understanding the Born radius via computer simulations and theory. J. Phys. Chem. B 1997, 101, 3596-3604. [112] Conway, B.E. Factors limiting applications of the historically significant Born equation: A critical review. In: Modern Aspects of Electrochemistry. Volume 35.; R.E. White, B.E. Conway, Ed.; Kluwer Academic/Plenum Publishers, New York, USA, 2002; pp 3596-3604. [113] Mukhopadhyay, A.; Fenley, A.T.; Tolokh, I.S.; Onufriev, A.V. Charge hydration asymmetry: The basic principles and how to use it to test and improve water models. J. Phys. Chem. B 2012, 116, 9776-9783. [114] Xie, D.; Volkmer, H.W. A modified nonlocal continuum electrostatic model for protein in water and its analytical solutions for ionic Born models. Commun. Comput. Phys. 2013, 13, 174-194. [115] Warwicker, J.; Watson, H.C. Calculation of the electric potential in the active site cleft due to α-helix dipoles. J. Mol. Biol. 1982, 157, 671-679. [116] Gersten, J.I.; Sapse, A.M. Generalization of the Born equation to nonspherical solvent cavities. J. Am. Chem. Soc. 1985, 107, 3786-3788. [117] Rashin, A.A.; Namboodiri, K. A simple method for the calculation of hydration enthalpies of polar molecules with arbitrary shapes. J. Phys. Chem. 1987, 91, 6003-6012. [118] Gilson, M.K.; Honig, B. Calculation of the total electrostatic energy of a macromolecular system: solvation energies, binding energies, and conformational analysis. Proteins: Struct. Funct. Genet. 1988, 4, 7-18. [119] Gomez-Jeria, J.S.; Morales-Lagos, D. Free energy of a charge distribution in a spheroidal cavity surrounded by concentric dielectric continua. J. Phys. Chem. 1990, 94, 3790-3795. [120] Madura, J.D.; Davis, M.E.; Gilson, M.K.; Wade, R.C.; Luty, B.A.; McCammon, J.A. Biological applications of electrostatic calculations and Brownian dynamics simulations. In: Reviews in computational chemistry. Volume 4.; Lipkowitz, K.B. & Boyd, D.B., Eds.; VCH Publishers Inc., New York, USA, 1994; pp 3790-3795. [121] Still, W.C.; Tempczyk, A.; Hawley, R.C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127-6129. [122] Onufriev, A.; Case, D.A.; Bashford, D. Effective Born radii in the generalized Born approximation: The importance of being perfect. J. Comput. Chem. 2002, 23, 1297-1304. [123] Onufriev, A.V.; Sigalov, G. A strategy for reducing gross errors in the generalized Born models of implicit solvation. J. Chem. Phys. 2011, 134, 164104/1-164104/15. [124] Sitkoff, D.; Sharp, K.A.; Honig, B. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 1994, 98, 1978-1988. [125] Ginovska, B.; Camaioni, D.M.; Dupuis, M.; Schwerdtfeger, C.A.; Charge-dependent cavity radii for an accurate dielectric continuum model of solvation with emphasis on ions: Aqueous solutes with oxo, hydroxo, amino, methyl, chloro, bromo, and fluoro functionalities. J. Phys. Chem. A 2008, 112, 10604-10613. [126] Hou, G.; Zhu, X.; Cui, Q. An implicit solvent model for SCC-DFTB with charge-dependent radii. J. Chem. Theory Comput. 2010, 6, 2303-2314. [127] Mukhopadhyay, A.; Aguilar, B.H.; Tolokh, I.S.; Onufriev, A.V. Introducing charge hydration asymmetry into the generalized Born model. J. Chem. Theory Comput. 2014, 10, 1788-1794.

50 ACS Paragon Plus Environment

Page 50 of 74

Page 51 of 74

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 Journal of Physical Chemistry

[128] Li, L.; Fennell, C.J.; Dill, K.A. Field-SEA: A model for computing the solvation free energies of nonpolar, polar, and charged solutes in water. J. Phys. Chem. B 2014, 118, 6431-6437. [129] Fawcett, W.R. The solvent dependence of ionic properties in solution in the limit of infinite dilution. Mol. Phys. 1998, 95, 507-514. [130] Fawcett, W.R. Liquids, solutions and interfaces.; Oxford Univ. Press, Oxford, UK, 2004. [131] Dogonadze, R.R.; Kornyshev, A.A. Polar solvent structure in the theory of ionic solvation. J. Chem. Soc. Faraday Trans. II 1974, 70, 1121-1132. [132] Kornyshev, A.A.; Rubinshtein, A.I.; Vorotyntsev, M.A. Model nonlocal electrostatics: . J. Phys. C: Solid State Phys. 1978, 11, 3307-3322. [133] Vorotyntsev, M.A. Model nonlocal electrostatics : II. Spherical interface. J. Phys. C: Solid State Phys. 1978, 11, 3323-3331. [134] Basilevsky, M.V.; Parsons, D.F. Nonlocal continuum solvation model with exponential susceptibility kernels. J. Chem. Phys. 1998, 108, 9107-9113. [135] Hildebrandt, A.; Blossey, R.; Rjasanow, S.; Kohlbacher, O.; Lenhof, H.-P. Novel formulation of nonlocal electrostatics. Phys. Rev. Lett. 2004, 93, 108104/1-108104/4. [136] Weggler, S.; Rutka, V.; Hildebrandt, A. A new numerical method for nonlocal electrostatics in biomolecular simulations. J. Comput. Phys. 2010, 229, 4059-4074. [137] Rottler, J.; Krayenhoff, B. Numerical studies of nonlocal electrostatic effects on the sub-nanoscale. J. Phys. Condens. Matter 2009, 21, 255901/1-255901/7. [138] Basilevsky, M.V.; Odinokov, A.V.; Nikitina, E.A.; Petrov, N.C. The dielectric continuum solvent model adapted for treating preferential solvation effects. J. Electroanal. Chem. 2011, 660, 339-346. [139] Schaaf, C.; Gekle, S. Dielectric response of the water hydration layer around spherical solutes. Phys. Rev. E 2015, 92, 032718/1-032718/9. [140] Dinpajooh, M.; Matyushov, D.V. Free energy of ion hydration: Interface susceptibility and scaling with the ion size. J. Chem. Phys. 2015, 143, 044511/1-044511/9. [141] Lebowitz, J.L.; Percus, J.K. Mean spherical model for lattice gases with extended hard cores and continuum fluids. Phys. Rev. 1966, 144, 251-258. [142] Wertheim, M.S. Exact solution of the mean spherical model for fluids of hard spheres with permanent electric dipole moments. J. Chem. Phys. 1971, 55, 4291-4298. [143] Adelman, S.A.; Deutch, J.M. Exact solution of the mean spherical model for strong electrolytes in polar solvents. J. Chem. Phys. 1974, 60, 3935-3949. [144] Blum, L. Solution of a model for the solvent-electrolyte interactions in the mean spherical approximation. J. Chem. Phys. 1974, 61, 2129-2133. [145] Fawcett, W.R.; Blum, L. The role of dipole-dipole interactions in the solvation of monoatomic monovalent ions in water on the basis of the mean spherical approximation. J. Electroanal. Chem. 1993, 355, 253-263. [146] Blum, L.; Torruella, A.J. Invariant expansion for two-body correlations: Thermodynamic functions, scattering, and the Ornstein-Zernike equation. J. Chem. Phys. 1972, 56, 303-310. [147] Blum, L. Invariant Expansion. II. The Ornstein-Zernike Equation for Nonspherical Molecules and an Extended Solution to the Mean Spherical Model. J. Chem. Phys. 1972, 57, 1862-1869. [148] Chandler, D.; Andersen, H.C. Optimized cluster expansions for classical fluids. II. Theory of molecular liquids. J. Chem. Phys. 1972, 57, 1930-1937. [149] Chandler, D.; Silbey, R. New and proper integral equations for site-site equilibrium correlations in molecular fluids. Mol. Phys. 1982, 46, 1335-1345. [150] Hirata, F.; Pettitt, M.; Rossky, P.J. Application of an extended RISM equation to dipolar and quadrupolar fluids. J. Chem. Phys. 1982, 77, 509-520. [151] Ratkova, E.L.; Palmer, D.S.; Fedorov, M.V. Solvation thermodynamics of organic molecules by the molecular integral equation theory: Approaching chemical accuracy. Chem. Rev. 2015, 115, 6312-6356. [152] Levesque, D.; Weis, J.J.; Patey, G.N. Charged hard spheres in dipolar hard sphere solvents. A model for electrolyte solutions. J. Chem. Phys. 1980, 72, 1887-1899. [153] Yu, H.; Karplus, M. A thermodynamic analysis of solvation. J. Chem. Phys. 1988, 89, 2366-2379.

51 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

[154] Fedorov, M.V.; Kornyshev, A.A. Unravelling the solvent response to neutral and charged solutes. Mol. Phys. 2007, 105, 1-16. [155] Dyer, K.M.; Perkyns, J.S.; Stell, G.; Pettitt, B.M. A site-renormalized molecular fluid theory. J. Chem. Phys. 2007, 129, 194506/1-194506/14. [156] Dyer, K.M.; Perkyns, J.S.; Stell, G.; Pettitt, B.M. A molecular site-site integral equation that yields the dielectric constant. J. Chem. Phys. 2008, 129, 104512/1-104512/9. [157] Dyer, K.M.; Perkyns, J.S.; Stell, G.; Pettitt, B.M. Site-renormalised molecular fluid theory: on the utility of a two-site model of water. Mol. Phys. 2009, 107, 423-431. [158] Misin, M.; Fedorov, M.V.; Palmer, D.S. Hydration free energies of molecular ions from theory and simulation. J. Phys. Chem. B 2016, 120, 975-983. [159] 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. [160] 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. [161] 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. [162] 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. [163] Rocklin, G.J.; Mobley, D.L.; Dill, K.A.; H¨ unenberger, P.H. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: An accurate correction scheme for electrostatic finite-size effects. J. Chem. Phys. 2013, 139, 184103/1-184103/32. [164] Reif, M.M.; Oostenbrink, C. Net charge changes in the calculation of relative ligand-binding free energies via classical atomistic molecular dynamics simulations. J. Comput. Chem. 2014, 35, 227-243. [165] Lin, Y.-L.; Aleksandrov, A.; Simonson, T.; Roux, B. An overview of electrostatic free energy computations for solutions and proteins. J. Chem. Theory Comput. 2014, 10, 2690-2709. [166] ˚ Aqvist, J.; Hansson, T. Analysis of electrostatic potential truncation schemes in simulations of polar solvents. J. Phys. Chem. B 1998, 102, 3837-3840. [167] Vorobjev, Y.N.; Hermans, J. A critical analysis of methods of calculation of a potential in simulated polar liquids : Strong arguments in favor of “Molecule-based” summation and of vacuum boundary conditions in Ewald summations. J. Phys. Chem. B 1999, 103, 10234-10242. [168] Babu, C.S.; Yang, P.-K.; Lim, C. On the charge and molecule based summations of solvent electrostatic potentials and the validity of electrostatic linear response in water. J. Biol. Phys. 2002, 28, 95-113. [169] Yang, P.-K.; Lim, C. Nonconvergence of the solute potential in an infinite solvent and its implications in continuum models. J. Phys. Chem. B 2002, 106, 12093-12096. [170] Heinz, T.N.; H¨ unenberger, P.H. Combining the lattice-sum and reaction-field approaches for evaluating longrange electrostatic interactions in molecular simulations. J. Chem. Phys. 2005, 123, 034107/1-034107/19. [171] 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. [172] Leung, K. Surface potential at the air-water interface computed using density functional theory. J. Phys. Chem. Lett. 2010, 1, 496-499. [173] Sedlmeier, F.; Netz, R.R. Solvation thermodynamics and heat capacity of polar and charged solutes in water. J. Chem. Phys. 2013, 138, 115101/1-115101/12. [174] Rienstra-Kiracofe, J.C.; Tschumper, G.S.; Schaefer, H.F.; Nandi, S.; Ellison, G.B. Atomic and molecular electron affinities: Photoelectron experiments and theoretical computations. Chem. Rev. 2002, 102, 231-282. [175] Rosenstock, H.M. The measurement of ionization and appearance potentials. Int. J. Mass Spectrom. Ion Phys. 1976, 20, 139-190. [176] Riviere, J.C. Work function: Measurements and results. In: Solid state surface science. Volume 1.; Green, M., Ed.; Decker, New York, USA, 1969; pp 139-190. [177] Hush, N.S. The free energies of hydration of gaseous ions. Austral. J. Sci. Res. 1948, 1, 480-493.

52 ACS Paragon Plus Environment

Page 52 of 74

Page 53 of 74

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 Journal of Physical Chemistry

[178] Randles, J.E.B. The real hydration energies of ions. Trans. Farad. Soc. 1956, 52, 1573-1581. [179] Riniker, S.; van Gunsteren, W.F. A simple, efficient polarisable coarse-grained water model for molecular dynamics simulations. J. Chem. Phys. 2011, 134, 084110/1-084110/12. [180] Schmid, N.; Christ, C.D.; Christen, M.; Eichenberger, A.P.; van Gunsteren, W.F. Architecture, implementation and parallelisation of the GROMOS software for biomolecular simulation. Comput. Phys. Commun. 2012, 183, 890-903. [181] van Gunsteren, W.F.; Billeter, S.R.; Eising, A.A.; H¨ unenberger, P.H.; Kr¨ uger, P.; Mark, A.E.; Scott, W.R.P.; Tironi, I.G. Biomolecular simulation: The GROMOS96 manual and user guide.; Verlag der Fachvereine, Z¨ urich, Switzerland, 1996. [182] Scott, W.R.P.; H¨ unenberger, P.H.; Tironi, I.G.; Mark, A.E.; Billeter, S.R.; Fennen, J.; Torda, A.E.; Huber, T.; Kr¨ uger, P.; van Gunsteren, W.F. The GROMOS biomolecular simulation program package. J. Phys. Chem. A 1999, 103, 3596-3607. [183] Hockney, R.W.; Eastwood, J.W. Computer simulation using particles.; McGraw-Hill, New York, USA, 1981. [184] H¨ unenberger, P.H. Lattice-sum methods for computing electrostatic interactions in molecular simulations. In: Simulation and theory of electrostatic interactions in solution: Computational chemistry, biophysics, and aqueous solution.; Hummer, G. & Pratt, L.R., Eds.; American Institute of Physics, New York, USA, 1999; pp 3596-3607. [185] Lennard-Jones, J.E. The equation of state of gases and critical phenomena. Physica 1937, 4, 941-956. [186] Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; Hermans, J. Interaction models for water in relation to protein hydration. In: Intermolecular Forces.; Pullman, B., Ed.; Reidel, Dordrecht, The Netherlands, 1981; pp 941-956. [187] 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 and Gibbs free energy of solvation from cluster-ion solvation data. J. Phys. Chem. A 1998, 102, 7787-7794. [188] Dreyer, I.; Dreyer, R.; Chalkin, V.A. Cations of astatine in aqueous solutions. Preparation and properties. Radiochem. Radioanal. Lett. 1978, 36, 389-398. [189] Dye, J.L. Compounds of alkali metal anions. Angew. Chem. Int. Ed. 1979, 18, 587-598. [190] Milanov, M.; Doberenz, V.; Khalkin, V.A.; Marinov, A. Chemical properties of positive singly-charged astatine ion in aqueous solution. J. Radioanal. Nucl. Chem. 1984, 83, 291-299. [191] Hagler, A.T.; Huler, E.; Lifson, S. Energy functions for peptides and proteins. I. Derivation of a consistent force field including the hydrogen bond from amide crystals. J. Am. Chem. Soc. 1974, 96, 5319-5327. [192] Lifson, S.; Hagler, A.T.; Dauber, P. Consistent force field studies of intermolecular forces in hydrogen-bonded crystals. 1. Carboxylic acids, amides, and the C=O···H hydrogen bonds. J. Am. Chem. Soc. 1979, 101, 51115121. [193] Geerke, D.P.; Oostenbrink, C.; van der Vegt, N.F.A.; van Gunsteren, W.F. An effective force field for molecular dynamics simulations of dimethyl sulfoxide and dimethyl sulfoxide-water mixtures. J. Phys. Chem. B 2004, 108, 1436-1445. [194] Horta, B.A.C.; Merz, P.T.; Fuchs, P.; Dolenc, J.; Riniker, S.; H¨ unenberger, P.H. A GROMOS-compatible force field for small organic molecules in the condensed phase: The 2016H66 parameter set. Submitted to J. Chem. Theory Comput. 2016. [195] Walser, R.; Mark, A.E.; van Gunsteren, W.F.; Lauterbach, M.; Wipff, G. The effect of force-field parameters on properties of liquids: Parametrization of a simple three-site model for methanol. J. Chem. Phys. 2000, 112, 10450-10459. [196] Tironi, I.G.; van Gunsteren, W.F. A molecular dynamics study of chloroform. Mol. Phys. 1994, 83, 381-403. [197] Tironi, I.G.; Fontana, P.; van Gunsteren, W.F. A molecular dynamics simulation study of liquid carbon tetrachloride. Mol. Sim. 1996, 18, 1-11. [198] Hockney, R.W. The potential calculation and some applications. Methods Comput. Phys. 1970, 9, 135-211. [199] Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H.J.C. Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327-341. [200] Berendsen, H.J.C.; Postma, J.P.M.; van Gunsteren, W.F.; di Nola, A.; Haak, J.R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684-3690. [201] Barker, J.A.; Watts, R.O. Monte Carlo studies of the dielectric properties of water-like models. Mol. Phys. 1973, 26, 789-792. [202] Tironi, I.G.; Sperb, R.; Smith, P.E.; van Gunsteren, W.F. A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 1995, 102, 5451-5459.

53 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

[203] 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. [204] Redlack, A.; Grindlay, J. The electrostatic potential in a finite ionic crystal. Can. J. Phys. 1972, 50, 2815-2825. [205] de Leeuw, S.W.; Perram, J.W.; Smith, E.R. Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants. Proc. R. Soc. Lond. A 1980, 373, 27-56. [206] Luty, B.A.; Tironi, I.G.; van Gunsteren, W.F. Lattice-sum methods for calculating electrostatic interactions in molecular simulations. J. Chem. Phys. 1995, 103, 3014-3021. [207] Deserno, M.; Holm, C. How to mesh up Ewald sums. II. An accurate error estimate for the particle-particleparticle-mesh algorithm. J. Chem. Phys. 1998, 109, 7694-7701. [208] H¨ unenberger, P.H. Optimal charge-shaping functions for the particle-particle–particle-mesh (P3 M) method for computing electrostatic interactions in molecular simulations. J. Chem. Phys. 2000, 113, 10464-10476. [209] Kastenholz, M.A.; H¨ unenberger, P.H. Influence of artificial periodicity and ionic strength in molecular dynamics simulations of charged biomolecules employing lattice-sum methods. J. Phys. Chem. B 2004, 108, 774-788. [210] Reif, M.M.; Kr¨ autler, V.; Kastenholz, M.A.; Daura, X.; H¨ unenberger, P.H. Explicit-solvent molecular dynamics simulations of a reversibly-folding β-heptapeptide in methanol: Influence of the treatment of long-range electrostatic interactions. J. Phys. Chem. B 2009, 113, 3112-3128. [211] Berendsen, H.J.C.; van Gunsteren, W.F.; Zwinderman, H.R.J; Geurtsen, R.G. Simulations of proteins in water. Ann. New York Acad. Sci. 1986, 482, 269-285. [212] Kirkwood, J.G. Statistical mechanics of fluid mixtures. J. Chem. Phys. 1935, 3, 300-313. [213] Baker, N.A.; H¨ unenberger, P.H.; McCammon, J.A. Polarization around an ion in a dielectric continuum with truncated electrostatic interactions. J. Chem. Phys. 1999, 110, 10679-10692. [214] Nijboer, B.R.A.; Ruijgrok, T.W. On the energy per particle in three- and two-dimensional Wigner lattices. J. Stat. Phys. 1988, 53, 361-382. [215] Bergdorf, M.; Peter, C.; H¨ unenberger, P.H. Influence of cutoff truncation and artificial periodicity of electrostatic interactions in molecular simulations of solvated ions: A continuum electrostatics study. J. Chem. Phys. 2003, 119, 9129-9144. [216] Peter, C.; van Gunsteren, W.F.; H¨ unenberger, P.H. A fast-Fourier-transform method to solve continuumelectrostatics problems with truncated electrostatic interactions: Algorithm and application to ionic solvation and ion-ion interaction. J. Chem. Phys. 2003, 119, 12205-12223. [217] Goldschmidt, V.M. Geochemische Verteilungsgesetze der Elemente VII: Die Gesetze der Krystallchemie. Skrifter Norske Vidensk. Akad. Oslo 1 Mat.-Nat. Kl. No. 1926, 2, 1-117. [218] Ashbaugh, H.S.; Wood, R.H. Reply to comment on “Electrostatic potentials and free energies of solvation of polar and charged molecules”. J. Phys. Chem. B 1998, 102, 3844-3845. [219] Reif, M.M.; Oostenbrink, C. Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: Restraints on solvent-generated electrostatic potential and solvent polarization. Theor. Chem. Acc. 2015, 134, 2/1-2/19. [220] Ma, J.C.; Dougherty, D.A. The cation-π interaction. Chem. Rev. 1997, 97, 1303-1324. [221] Salonen, L.M.; Ellermann, M.; Diederich, F. Aromatic rings in chemical and biological recognition: Energetics and structures. Angew. Chem. Int. Ed. Engl. 2011, 50, 4808-4842. [222] Magnico, P. Influence of the ion-solvent interactions on ionic transport through ion-exchange-membranes. J. Membrane Sci. 2013, 442, 272-285. [223] Dougherty, D.A. The cation-π interaction. Acc. Chem. Res. 2013, 46, 885-893. [224] Schneider, H.; Vogelhuber, K.M.; Schinle, F.; Weber, J.M. Aromatic molecules in anion recognition: Electrostatics versus H-bonding. J. Am. Chem. Soc. 2007, 129, 13022-13026. [225] Geronimo, I.; Singh, N.J.; Kim, K.S. Can electron-rich π systems bind anions. J. Chem. Theory Comput. 2011, 7, 825-829. [226] Tresca, B.W.; Hansen, R.J.; Chau, C.V.; Hay, B.P.; Zakharov, L.N.; Haley, M.H.; Johnson, D.W. Substituent Effects in CH hydrogen bond interactions: Linear free energy relationships and inuence of anions. J. Am. Chem. Soc. 1992, 137, 14959-14967. [227] Lin, Z.; van Gunsteren, W.F. Effects of polarizable solvent models upon the relative stability of an α-helical and a β-hairpin structure of an alanine decapeptide. J. Chem. Theory Comput. 2015, 11, 1983-1986.

54 ACS Paragon Plus Environment

Page 54 of 74

Page 55 of 74

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 Journal of Physical Chemistry

[228] Premkumar, J.R.; Sastry, G.N. Cation-alkane interaction. J. Phys. Chem. A 2014, 118, 11388-11398. [229] Neela, Y.I.; Sastry, N. Theoretical investigation of anion (F− , Cl− ) and cation (Na+ ) interactions with substituted benzene [C6 H6−n Yn (Y=-F,-CN,-NO2 ;n=1-6). Mol. Phys. 2015, 113, 137-148. [230] Anslyn, E.V.; Dougherty, D.A. Modern physical organic chemistry.; Univ. Science, 2006. [231] Szklarczyk, O.M.; Bachmann, S.J.; van Gunsteren, W.F. A polarizable empirical force field for molecular dynamics simulation of liquid hydrocarbons. J. Comput. Chem. 2014, 35, 789-801. [232] Harder, E.; Roux, B. On the origin of the electrostatic potential difference at a liquid-vacuum interface. J. Chem. Phys. 2008, 129, 234706/1-234706/9. [233] Kathmann, S.M.; I-Feng, W.K.; Mundy, C.J. Electronic effects on the surface potential at the vapor-liquid interface of water. J. Am. Chem. Soc. 2008, 130, 16556-16561. [234] Kathmann, S.M.; Kuo, I.-F.W.; Mundy, C.J.; Schenter, G.K. Understanding the surface potential of water. J. Phys. Chem. B 2011, 115, 4369-4377. [235] Beck, T.L. The influence of water interfacial potentials on ion hydration in bulk water and near interfaces. Chem. Phys. Lett. 2013, 561-562, 1-13. [236] Lu, X.; Cui, Q. Charging free energy calculations using the generalized solvent boundary potential (GSBP) and periodic boundary condition: A comparative analysis using ion solvation and oxidation free energy in proteins. J. Phys. Chem. B 2013, 117, 2005-2018. [237] Remsing, R.C.; Baer, M.D.; 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. [238] Guggenheim, E.A. The conceptions of electrical potential differences between two phases and the individual activities of ions. J. Phys. Chem. 1929, 33, 842-849. [239] Guggenheim, E.A. On the conception of electrical potential difference between two phases. II. J. Phys. Chem. 1930, 34, 1540-1543. [240] Rockwood, A.L. Meaning and measurability of single-ion activities, the thermodynamic foundations of pH, and the Gibbs free energy for the transfer of ions between dissimilar materials. Chem. Phys. Chem. 2015, 16, 1978-1991. [241] Stillinger, F.H.; Rahman, A. Improved simulation of liquid water by molecular dynamics. J. Chem. Phys. 1974, 60, 1545-1557. [242] Berendsen, H.J.C.; Grigera, J.R.; Straatsma, T.P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269-6271. [243] Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926-935. [244] Mahoney, M.W.; Jorgensen, W.L. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J. Chem. Phys. 2000, 112, 8910-8922. [245] Warren, G.L.; Patel, S. Hydration free energies of monovalent ions in transferable intermolecular potential four point fluctuating charge water: An assessment of simulation methodology and force field performance and transferability. J. Chem. Phys. 2007, 127, 064509/1-064509/19. [246] Ashbaugh, H.S.; Wood, R.H. Effects of long-range electrostatic potential truncation on the free energy of ionic hydration. J. Chem. Phys. 1997, 106, 8135-8139. [247] Hummer, G.; Pratt, L.R.; Garcia, A.E.; Berne, B.J.; Rick, S.W. Electrostatic potentials and free energies of solvation of polar and charged molecules. J. Phys. Chem. B 1997, 101, 3017-3020. [248] Hummer, G.; Pratt, L.R.; Garcia, A.E. Ion sizes and finite-size corrections for ionic-solvation free energies. J. Chem. Phys. 1997, 107, 9275-9277. [249] Hummer, G.; Pratt, L.R.; Garcia, A.E.; Garde, S.; Berne, B.J.; Rick, S.W. Reply to comments on “Electrostatic potentials and free energies of solvation of polar and charged molecules”. J. Phys. Chem. B 1998, 102, 38413843. [250] Darden, T.; Pearlman, D.; Pedersen, L.G. Ionic charging free energies : Spherical versus periodic boundary conditions. J. Chem. Phys. 1998, 109, 10921-10935. [251] Sakane, S.; Ashbaugh, H.S.; Wood, R.H. Continuum corrections to the polarization and thermodynamic properties of Ewald sum simulations for ions and ion pairs at infinite dilution. J. Phys. Chem. B 1998, 102, 5673-5682. [252] Christou, N.I.; Whitehouse, J.S.; Nicholson, D.; Parsonage, N.G. Studies of high density water films by computer simulation. Mol. Phys. 1985, 55, 397-410.

55 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

[253] Aloisi, G.; Guidelli, R.; Jackson, R.A.; Clark, S.M.; Barnes, P. The structure of water at a neutral interface. J. Electroanal. Chem 1986, 206, 131-137. [254] Madura, J.D.; Pettitt, B.M. Effect of truncating long-range interactions in aqueous ionic solution simulations. Chem. Phys. Lett. 1988, 150, 105-108. [255] Matsumoto, M.; Kataoka, Y. Study on liquid-vapor interface of water. I. Simulation results of thermodynamic properties and orientational structure. J. Chem. Phys. 1988, 88, 3233-3245. [256] Barraclough, C.G.; McTigue, P.T.; Ng, Y.L. Surface potentials of water, methanol and water+methanol mixtures. J. Electroanal. Chem. 1992, 329, 9-24. [257] Sokhan, V.P.; Tildesley, D.J. The free surface of water: molecular orientation, surface potential and nonlinear susceptibility. Mol. Phys. 1997, 92, 625-640. [258] Parfenyuk, V.I. Surface potential at the gas-aqueous solution interface. Colloid J. 2002, 64, 588-595. [259] Wick, C.D.; Dang, L.X.; Jungwirth, P. Simulated surface potentials at the vapor-water interface for the KCl aqueous electrolyte solution. J. Chem. Phys. 2006, 125, 024706/1-024706/4. [260] Brodskaya, E.N.; The molecular dynamics simulation of water clusters. Mol. Phys. 1987, 62, 251-265. [261] Wilson, M.A.; Pohorille, A.; Pratt, L.R. Molecular dynamics of the water liquid-vapor interface. J. Phys. Chem. 1987, 91, 4873-4878. [262] Wilson, M.A.; Pohorille, A.; Pratt, L.R. Surface potential of the water liquid-vapor interface. J. Chem. Phys. 1988, 88, 3281-3285. [263] Wilson, M.A.; Pohorille, A.; Pratt, L.R. Comment on “Study on the liquid-vapor interface of water. I. Simulation results of thermodynamic properties and orientational structure”. J. Chem. Phys. 1989, 90, 5211-5213. [264] Pohorille, A.; Wilson, M.A. Viewpoint 9 - Molecular structure of aqueous interfaces. J. Mol. Struct. (Theochem) 1993, 284, 271-298. [265] Brodskaya, E.N.; Zakharov, V.V. Computer simulation study of the surface polarization of pure polar liquids. J. Chem. Phys. 1995, 102, 4595-4599. [266] Ashbaugh, H.S. Influence of potential truncation on anisotropic systems. Mol. Phys. 1999, 97, 433-437. [267] Dang, L.X.; Chang, T.-M. Molecular mechanism of ion binding to the liquid/vapor interface of water. J. Phys. Chem. B 2002, 106, 235-238. [268] Patel, S.; Brooks III, C.L. Revisiting the hexane-water interface via molecular dynamics simulations using nonadditive alkane-water potentials. J. Chem. Phys. 2006, 124, 204706/1-204706/14. [269] van Buuren, A.R.; Marrink, S.-J.; Berendsen, H.J.C. Characterisation of aqueous interfaces with different hydrophobicities by molecular dynamics. Colloids Surfaces A: Physicochem. Eng. Aspects 1995, 102, 143-157. [270] Gan, W.; Wu, D.; Zhang, Z.; Feng, R.-R.; Wang, H.-F. Polarization and experimental configuration analyses of sum frequency generation vibrational spectra, structure, and orientational motion of the air/water interface. J. Chem. Phys. 2006, 124, 114705/1-114705/15. [271] Ishiyama, T.; Takahashi, H.; Morita, A. Molecular dynamics simulations of surface-specific bonding of the hydrogen network of water: A solution to the low sum-frequency spectra. Phys. Rev. B 2012, 86, 035408/1035408/10. [272] Kessler, J.; Elgabarty, H.; Spura, T.; Karhan, K.; Partovi-Azar, P.; Hassanall, A.A.; K¨ uhne, T.D. Structure and dynamics of the instantaneous water/vapor interface revisited by path-integral and ab initio molecular dynamics simulations. J. Phys. Chem. B 2015, 119, 10079-10086. [273] Lee, C.Y.; McCammon, J.A.; Rossky, P.J. The structure of liquid water at an extended hydrophobic surface. J. Chem. Phys. 1984, 80, 4448-4455. [274] Alper, H.E.; Bassolino-Klimas, D.; Stouch, T.R. The limiting behavior of water hydrating a phospholipid monolayer: A computer simulation study. J. Chem. Phys. 1993, 99, 5547-5559. [275] Gordillo, M.C.; Mart´ı, J. Molecular dynamics description of a layer of water molecules on a hydrophobic surface. J. Chem. Phys. 2002, 117, 3425-3430. [276] Strazdaite, S.; Versluis, J.; Bakker, H.J. Water orientation at hydrophobic interfaces. J. Chem. Phys. 2015, 143, 084708/1-084708/6. [277] Wohlfahrt, C. Pure liquids: Data. In: Landolt-B¨ ornstein. Numerical data and functional relationships in science and technology. Volume 6.; Madelung, O., Ed.; Springer, Berlin, Germany, 1991; pp 084708/1-084708/6.

56 ACS Paragon Plus Environment

Page 56 of 74

Page 57 of 74

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 Journal of Physical Chemistry

[278] Riniker, S.; Kunz, A.-P.E.; van Gunsteren, W.F. On the calculation of the dielectric permittivity of molecular models in the liquid phase. J. Chem. Theory Comput. 2011, 7, 1469-1475. [279] Smith, P.E.; van Gunsteren, W.F. Consistent dielectric properties of the simple point charge and extended simple point charge water models at 277 and 300 K. J. Chem. Phys. 1994, 100, 3169-3174. [280] Neumann, M. Dipole moment fluctuation formulas in computer simulations of polar systems. Mol. Phys. 1983, 50, 841-858.

57 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 58 of 74

Tables and figures Table 1: Index 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 61

Common name hexane heptane isooctane cyclohexane carbontetrachloride dioxane benzene p-xylene triethylamine toluene m-xylene perchlorethylene o-xylene carbon disulfide dibutylether trichlorethylene diethylether anisole chloroform butanamine bromobenzene chlorobenzene ethylacetate acetic acid methylaniline aniline dimethoxyethane lutidine tetrahydrofuran trifluoroacetic acid dichloromethane dichloroethane tert-butanol pyridine cyclohexanone acetophenone butanol butanone isopropanol propanol acetone ethanol benzonitrile hexamethylphosphoric triamide methanol nitrobenzene acetonitrile dimethylformamide nitromethane dimethylacetamide ethyleneglycol dimethylsulfoxide formic acid water formamide water-methanol (96%) water-ethanol (80%) water-ethanol (60%) water-ethanol (50%) water-acetone (80%) water-acetone (70%)

IUPAC name hexane heptane 2,2,4-trimethylpentane cyclohexane tetrachloromethane 1,4-dioxane benzene 1,4-dimethylbenzene N,N-diethylethanamine methylbenzene 1,3-dimethylbenzene tetrachloroethene 1,2-dimethylbenzene carbon disulfide 1-butoxybutane trichloroethene ethoxyethane methoxybenzene trichloromethane butylamine bromobenzene chlorobenzene ethylethanoate ethanoic acid N-phenylmethylamine phenylamine 1,2-dimethoxyethane 2,6-dimethylpyridine oxolane 2,2,2-trifluoroethanoic acid dichloromethane 1,2-dichloroethane 2-methylpropan-2-ol pyridine cyclohexanone 1-phenylethanone butan-1-ol butan-2-one propan-2-ol propan-1-ol propan-2-one ethanol benzonitrile N-[bis(dimethylamino)phosphoryl]N-methylmethanamine methanol nitrobenzene acetonitrile N,N-dimethylformamide nitromethane N,N-dimethylacetamide ethane-1,2-diol dimethylsulfoxide methanoic acid oxidane methanamide oxidane-methanol (96%) oxidane-ethanol (80%) oxidane-ethanol (60%) oxidane-ethanol (50%) oxidane-propan-2-one (80%) oxidane-propan-2-one (70%)

Abbreviation

CCL BZN

TOL

DEE CHL BAN

EAE ACA

PYRI

2PL PPN ETL

MET

DAMD DMS H2O

ǫS 1.9 1.9 1.9 2.0 2.2 2.2 2.3 2.3 2.4 2.4 2.4 2.4 2.6 2.6 3.1 3.4 4.2 4.3 4.7 5.2 5.4 5.6 6.1 6.2 6.2 7.0 7.2 7.3 8.0 8.3 8.9 10.4 11.5 12.9 15.8 17.4 17.7 18.3 19.1 20.4 20.8 24.3 25.1 29.8

A 0.01 0.00 0.01 0.02 0.09 0.19 0.15 0.06 0.08 0.13 0.04 0.10 0.06 0.10 0.06 0.16 0.12 0.21 0.42 0.15 0.22 0.20 0.21 0.93 0.40 0.36 0.21 0.18 0.17 1.72 0.33 0.30 0.45 0.24 0.25 0.23 0.61 0.23 0.59 0.63 0.25 0.66 0.30 0.00

B -0.01 0.00 -0.03 0.06 0.34 0.67 0.59 0.50 0.19 0.54 0.50 0.25 0.53 0.38 0.28 0.54 0.34 0.74 0.73 0.17 0.66 0.65 0.59 0.13 1.07 1.19 0.50 0.81 0.67 0.00 0.80 0.82 0.50 0.96 0.79 0.90 0.43 0.74 0.44 0.44 0.81 0.45 0.87 1.07

Σ -0.74 -0.77 -0.74 -0.75 -0.71 -0.63 -0.68 -0.84 -0.67 -0.70 -0.88 -0.65 -0.85 -0.71 -0.75 -0.64 -0.65 -0.61 -0.15 -0.51 -0.56 -0.60 -0.55 1.20 -0.33 -0.46 -0.51 -0.70 -0.67 2.96 -0.37 -0.45 0.01 -0.63 -0.54 -0.63 0.38 -0.57 0.33 0.42 -0.55 0.48 -0.47 -1.20

Σ∗ -0.78 -0.82 -0.77 -0.82 -0.88 -0.91 -0.94 -1.11 -0.77 -0.95 -1.16 -0.77 -1.13 -0.89 -0.91 -0.87 -0.80 -0.91 -0.33 -0.56 -0.82 -0.86 -0.78 1.59 -0.69 -0.91 -0.70 -1.06 -0.96 3.86 -0.64 -0.74 -0.05 -1.03 -0.85 -1.00 0.45 -0.86 0.39 0.50 -0.87 0.57 -0.79 -1.78

Π 0.40 0.40 0.39 0.43 0.59 0.77 0.72 0.65 0.52 0.69 0.64 0.55 0.66 0.61 0.55 0.70 0.60 0.81 0.88 0.53 0.78 0.77 0.74 0.80 1.03 1.07 0.70 0.83 0.76 1.04 0.88 0.88 0.79 0.92 0.85 0.89 0.82 0.82 0.82 0.83 0.86 0.85 0.90 0.88

π∗ -0.08 -0.08

33.0 34.7 36.7 36.7 36.7 39.6 40.6 46.4 56.1 78.4 109.6

0.75 0.29 0.37 0.30 0.39 0.27 0.78 0.34 1.18 1.00 0.66 0.76 0.75 0.80 0.82 0.62 0.66

0.50 0.86 0.86 0.93 0.92 0.97 0.84 1.08 0.51 1.00 0.99 0.61 0.65 0.77 0.80 0.70 0.74

0.66 -0.48 -0.31 -0.49 -0.29 -0.57 0.59 -0.46 1.59 1.00 0.27 0.64 0.60 0.66 0.69 0.30 0.37

0.77 -0.81 -0.59 -0.84 -0.59 -0.96 0.55 -0.87 1.93 1.00 0.08 0.70 0.64 0.66 0.69 0.24 0.31

0.90 0.89 0.92 0.93 0.96 0.94 1.07 1.01 1.07 1.22 1.09 0.96 0.97 1.04 1.06 0.94 0.98

0.60 1.01 0.75 0.88 0.85 0.88 0.92 1.00

0.00 0.28 0.55 0.59 0.43 0.14 0.54 0.47

0.24 0.53 0.27 0.73 0.58 0.79 0.71 0.55 0.64

0.53 0.80 0.58 0.50 0.81 0.41 0.87 0.76 0.90 0.47 0.67 0.48 0.52 0.71 0.54 0.90

1.09 0.97

Polarity and asymmetric solvation properties of 55 pure solvents plus 6 solvent mixtures. The solvents are referred to by their common names, IUPAC names, and an abbreviation (only for the solvents considered in the simulations). The static relative dielectric permittivities ǫS from Ref. 277 at 1 bar and 298.2 K (293.2 K for lutidine, isooctane and dibutylether) are also indicated. The relative acities (A) and basities (B) are taken from Swain et al. 2 (see Table III therein; the values implicitly pertain to ambient conditions of pressure and temperature). The asymmetries (Σ or Σ∗ ) and the polarities Π are calculated from A and B using Eqs. 30-32. The polarity parameter π ∗ of the the Kamlet-Taft scale is also reported from Ref. 30 (whenever available) for comparison. The data for the 55 pure solvents is also shown graphically in Figure 8.

58 ACS Paragon Plus Environment

Page 59 of 74

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 Journal of Physical Chemistry

Table 2: Solventa

Ref.

NS

L [nm]

ǫS

ǫ′S

γ ′S [10−3 e·nm2 ]

ρ′S [kg·m−3 ]

Ref. site

SAMs NIMs H2O DMS DAMD MET ETL PPN 2PL PYRI ACA EAE BAN CHL DEE TOL BZN CCL

-

1000 1000 999 1023 512 999 512 512 512 512 512 512 512 999 512 512 512 -m

3.13o 3.13o 3.12 4.95 4.30 4.07 3.68 3.93 4.12 4.10 3.60 4.42 4.37 5.09 4.51 4.52 4.23 -m

78.4i 46.4i 39.6i 33.0i 24.3i 20.8i 19.1i 12.9i 6.2i 6.1i 5.2i 4.7i 4.2i 2.4i 2.3i 2.2i

-b -d 66.7e 39.5e 20.3k 27.8e 17.1k 13.9k 10.3k 6.1k 3.1k 5.0k 3.2k 2.6e 9.4k 1.1k 1.1k -m

-b -d 8.20f 1.45g 5.4h 10.22g 6.3h -0.35h 1.5h 17.4h 2.1h -4.0h 10.4h -7.21g 4.5h -14.4h 5.1h 0f

922 1735 985j 1095j 933j 789j 788j 815j 732j 978j 1096j 869j 747j 1503j 688j 850j 879j 1601n

-c -c O S C(carbonyl) O O C(carbonyl) C(hydroxyl) N C(carbonyl) C(carbonyl) N C O C(aromatic) C(all) -m

186 193 194 195 194 194 194 194 194 194 194 196 194 194 194 197

Rrdf [nm] Na+ Na− -b -b -d -d 0.24 0.22 0.37 0.35 0.35 0.52 0.24 0.22 0.24 0.22 0.35 0.29 0.34 0.30 0.25 0.43 0.34 0.30 0.35 0.50l 0.24 0.23 0.45 0.31 0.24 0.40 0.30 0.33 0.30 0.33 -m -m

(a) Abbreviations as introduced in Section II.1 (see also Table 1 for IUPAC names): spherical anisotropically-charged solvent model (SAM), non-spherical isotropically-charged solvent model (NIM), water (H2O), dimethylsulfoxide (DMS), dimethylacetamide (DAMD), methanol (MET), ethanol (ETL), acetone (PPN), isopropanol (2PL), pyridine (PYRI), acetic acid (ACA), ethylacetate (EAE), butanamine (BAN), chloroform (CHL), diethylether (DEE), toluene (TOL), benzene (BZN) and carbontetrachloride (CCL); (b) the values are reported in Table 4; (c) the central (SAM) or the unscaled (NIM) particle is chosen as the reference site; (d) the values are reported in Table 5; (e) Ref. 278 ; (f) from method I (Suppl. Mat. Table S.2); (g) from method III (Suppl. Mat. Table S.2); (h) from method IV (tentative estimate only; Suppl. Mat. Table S.2); (i) Ref.277 ; (j) calculated from NS and L using Eq. 7, where MS is the molar mass of the solvent (Suppl. Mat. Table S.2); (k) Ref. 194 ; (l) corresponds to the global maximum (height 1.6) of the radial distribution function (the first peak is at about 0.28 nm with height 0.4); (m) no simulations were performed (because all partial charges are zero); (n) Ref.197 , at 293 K; (o) calculated from NS and ρ′S using Eq. 7, where MS is the molar mass of the solvent (17.0 and 32.0 g·mol−1 for the SAM and NIM, respectively).

Properties of the solvent models considered and parameters of the corresponding simulations. The reported entries are the source reference for the model parameters, the number NS of solvent molecules in the cubic computational box, the corresponding box-edge length L, the experimental static relative dielectric permittivity ǫS , the corresponding calculated value ǫ′S for the model, the effective quadrupole-moment trace γ ′S of the model (see Suppl. Mat. Section S.II for the calculation details and Suppl. Mat. Tables S.1 and S.2 for the corresponding exclusion potentials ξS′ ), the density ρ′S (related to NS and L via Eq. 7), the reference site selected for the model (multiple equivalent atoms for TOL and BZN) and the radii Rrdf for the ions Na+ and Na− (first peak in the radial distribution function between ion and reference site of the solvent).

59 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 60 of 74

Table 3:

Na+ , K+ , F+ , Rb+ , I+ ,

Na− K− F− Rb− I−

O(SPC)

C6ii [10−5 kJ·mol−1 ·nm6 ] 7.8840 129.60 205.50 309.60 4445.2 C6ww −5 [10 kJ·mol−1 ·nm6 ] 261.73

ii C12 [10−8 kJ·mol−1 ·nm12 ] 7.2897 186.59 251.13 641.90 56286 ww C12 −8 [10 kJ·mol−1 ·nm12 ] 263.41

σ iw [nm] 0.314 0.326 0.322 0.336 0.391 σ ww [nm] 0.317

σ ˜ iw [nm] 0.264 0.294 0.295 0.307 0.364 σ ˜ ww [nm] 0.293

Rgld [nm] 0.098 0.133 0.133 0.149 0.220

Lennard-Jones interaction parameters of the five ion pairs considered and of the oxygen atom of the SPC water model. The reported entries are the homoionic Lennard-Jones interaction ii (repulsion) for the ions (based on the C iw and C iw parameters C6ii (dispersion) and C12 6 12 parameters of set L in Ref. 84 ; see Tables II and III therein) as well as the corresponding ww for the oxygen atom of the SPC model, 186 the collision diameters parameters C6ww and C12 iw ww σ or σ (zero point of the Lennard-Jones interaction curve with the oxygen atom of SPC) as well as the thermal radii σ ˜ iw and σ ˜ ww (location of the kB T energy in the LennardJones interaction curve with the oxygen atom of SPC, where kB is Boltzmann’s constant and T = 298.15 K), and the Goldschmidt radius 217 (ions only). The parameters of the SPC oxygen atom are also used for the central particle of the SAMs and the unscaled particle of the NIMs. All interaction parameters are connected in this work by application of a geometric mean combination rule. 181, 191, 192

60 ACS Paragon Plus Environment

Page 61 of 74

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 Journal of Physical Chemistry

Table 4: Index

1 5 7 11 13 2 6 10 14 17 3 8 12 16 19 4 9 15 18 20

qs [e] 0.10 0.20 0.30 0.40 0.50 0.10 0.20 0.30 0.40 0.50 0.10 0.20 0.30 0.40 0.50 0.10 0.20 0.30 0.40 0.50

b [nm] 0.06 0.06 0.06 0.06 0.06 0.08 0.08 0.08 0.08 0.08 0.10 0.10 0.10 0.10 0.10 0.12 0.12 0.12 0.12 0.12

ǫ′S

1.3 2.4 4.9 9.5 18.6 1.5 3.8 9.4 22.5 51.9 1.9 6.0 16.7 42.4 109.5 2.4 8.4 25.0 80.8 a

γ ′S e·nm2 ]

[10−3

0.36 0.72 1.08 1.44 1.80 0.64 1.28 1.92 2.56 3.20 1.00 2.00 3.00 4.00 5.00 1.44 2.88 4.32 5.76 7.20

Na+ 0.28 0.27 0.26 0.26 0.25 0.28 0.27 0.25 0.25 0.25 0.27 0.26 0.26 0.25 0.25 0.27 0.26 0.25 0.25 0.24

Na− 0.27 0.26 0.25 0.24 0.23 0.26 0.25 0.24 0.23 0.23 0.25 0.24 0.23 0.22 0.22 0.24 0.22 0.22 0.21 0.20

K+ 0.31 0.31 0.30 0.29 0.29 0.31 0.30 0.29 0.29 0.29 0.31 0.30 0.29 0.29 0.29 0.31 0.30 0.30 0.29 0.28

K− 0.30 0.29 0.29 0.28 0.28 0.30 0.29 0.28 0.27 0.27 0.29 0.28 0.27 0.26 0.25 0.29 0.27 0.26 0.25 0.24

Rrdf [nm] F+ F− 0.31 0.31 0.30 0.30 0.30 0.29 0.30 0.29 0.29 0.28 0.31 0.31 0.30 0.29 0.30 0.28 0.30 0.28 0.29 0.27 0.31 0.30 0.30 0.28 0.30 0.27 0.29 0.26 0.29 0.26 0.31 0.29 0.30 0.28 0.30 0.26 0.30 0.25 0.29 0.25

Rb+ 0.32 0.32 0.31 0.31 0.30 0.32 0.32 0.31 0.30 0.30 0.32 0.32 0.31 0.31 0.31 0.32 0.31 0.31 0.30 0.30

Rb− 0.32 0.31 0.31 0.30 0.29 0.31 0.31 0.30 0.29 0.28 0.31 0.30 0.29 0.28 0.27 0.31 0.29 0.28 0.27 0.26

I+ 0.39 0.38 0.38 0.37 0.37 0.38 0.38 0.37 0.37 0.37 0.38 0.38 0.38 0.37 0.37 0.38 0.38 0.38 0.37 0.37

I− 0.39 0.37 0.37 0.36 0.36 0.38 0.37 0.36 0.35 0.35 0.38 0.36 0.35 0.35 0.34 0.37 0.35 0.34 0.33 0.32

(a) For this extremely polar solvent model, the dielectric permittivity calculation failed to converge (formation of a glassy state).

Properties of the spherical anisotropically-charged solvent models (SAMs). See Figure 1 for the definition of a SAM. For each model, the reported data includes the charge qs , the bond length b, the relative dielectric permittivity ǫ′S (see Suppl. Mat. Section S.I for the calculation details 279, 280 ), the (effective) quadrupole moment trace γ ′S = γS′ (from method I, i.e. calculated analytically as γS′ = qs b2 , see Suppl. Mat. Section S.II; the corresponding exclusion potentials are given by Eq. 15 as ξS′ = Cγ ′S with C = 98.4 V·e−1 ·nm−2 ), and the radii Rrdf of the ions considered in the solvent (first peak in the radial distribution function between ion and solvent, with the central particle as reference solvent site). The molar mass MS of all SAMs is 17.0 g·mol−1 . These solvents were always simulated at a density ρ′S = 922 kg·m−3 (adjusted from the density of SPC water at 298.15 K and 1 bar, and not necessarily equal to the equilibrium density of the SAM under these conditions). Only the properties of the SAMs with qs > 0 are reported as those of the corresponding SAMs with qs < 0 can easily be deduced (same b and ǫ′S , oppositely-signed qs and γ ′S , swapped Rrdf for X + and X − ). The index indicates the ordering according to increasing ǫ′S .

61 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

Page 62 of 74

Table 5: Index

3 4 1 2 7 6 8 5 9 10 11 12 14 13 15 16 18 17 19 20

qs [e]

α

0.10 0.10 0.10 0.10 0.20 0.20 0.20 0.20 0.30 0.30 0.30 0.30 0.40 0.40 0.40 0.40 0.50 0.50 0.50 0.50

0.4 0.5 0.6 0.7 0.4 0.5 0.6 0.7 0.4 0.5 0.6 0.7 0.4 0.5 0.6 0.7 0.4 0.5 0.6 0.7

ǫ′S

4.7 4.7 4.6 4.6 16.7 16.6 17.2 16.5 34.6 35.4 36.3 36.4 59.7 57.0 60.8 65.8 99.6 83.1 a a

γ ′S e·nm2 ]

[10−3

3.88 3.91 3.91 3.91 7.87 7.84 7.81 7.84 11.72 11.73 11.74 11.72 15.97 15.80 15.43 15.64 18.51 19.71 19.69 20.64

Na+ 0.27 0.27 0.27 0.26 0.26 0.26 0.26 0.26 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.24 0.24 0.24

Na− 0.35 0.37 0.39 0.40 0.32 0.35 0.36 0.39 0.30 0.32 0.34 0.36 0.28 0.30 0.33 0.34 0.27 0.29 0.33 0.33

K+ 0.30 0.30 0.30 0.29 0.29 0.30 0.30 0.29 0.29 0.29 0.29 0.29 0.29 0.29 0.29 0.28 0.28 0.28 0.28 0.27

K− 0.36 0.39 0.41 0.43 0.34 0.36 0.38 0.41 0.32 0.34 0.35 0.37 0.31 0.32 0.35 0.35 0.30 0.31 0.34 0.36

Rrdf [nm] F+ F− 0.31 0.37 0.30 0.39 0.30 0.41 0.29 0.43 0.30 0.34 0.30 0.36 0.29 0.39 0.30 0.41 0.29 0.32 0.29 0.34 0.29 0.35 0.29 0.38 0.29 0.31 0.29 0.33 0.30 0.35 0.28 0.35 0.29 0.30 0.29 0.31 0.28 0.34 0.27 0.35

Rb+ 0.32 0.31 0.31 0.31 0.31 0.31 0.31 0.30 0.31 0.30 0.31 0.30 0.31 0.30 0.30 0.30 0.30 0.31 0.30 0.29

Rb− 0.37 0.40 0.42 0.44 0.34 0.37 0.40 0.42 0.33 0.35 0.36 0.38 0.32 0.34 0.35 0.36 0.31 0.33 0.35 0.36

I+ 0.38 0.37 0.37 0.36 0.37 0.37 0.37 0.36 0.37 0.37 0.36 0.36 0.37 0.37 0.37 0.36 0.37 0.37 0.37 0.37

I− 0.41 0.43 0.46 0.49 0.40 0.41 0.42 0.42 0.39 0.40 0.40 0.41 0.38 0.39 0.42 0.40 0.37 0.38 0.40 0.42

(a) For these two extremely polar solvent models, the dielectric permittivity calculation failed to converge (formation of a glassy state).

Properties of the non-spherical isotropically-charged solvent models (NIMs). See Figure 1 for the definition of a NIM, where b is set to 0.2 nm in all cases. For each model, the reported data includes the charge qs , the radius-scaling factor α (Eq. 6), the relative dielectric permittivity ǫ′S (see Suppl. Mat. Section S.I for the calculation details 279, 280 ), the effective quadrupolemoment trace γ ′S (from method III, see Suppl. Mat. Section S.II for the calculation details and Suppl. Mat. Table S.1 for the corresponding exclusion potentials ξS′ ), and the radii Rrdf of the ions considered in the solvent (first peak in the radial distribution function between ion and solvent, with the unscaled particle as a reference solvent site). The molar mass MS of all NIMs is 32.0 g·mol−1 . These solvents were always simulated at a density of ρ′S = 1735 kg·m−3 (adjusted from the density of SPC water at 298.15 K and 1 bar, and not necessarily equal to the equilibrium density of the NIM under these conditions). Only the properties of the NIMs with qs > 0 are reported as those of the corresponding NIMs with qs < 0 can easily be deduced (same α and ǫ′S , oppositely-signed qs and γ ′S , swapped Rrdf for X + and X − ). The index indicates ordering according to increasing ǫ′S .

62 ACS Paragon Plus Environment

Page 63 of 74

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 Journal of Physical Chemistry

Table 6: Solvent

H2O

DMS DAMD MET ETL PPN 2PL PYRI ACA EAE BAN CHL DEE TOL BZN CCL

ǫ′S

Scheme  

LS BM   CM  LS BM 39.5   CM CM 20.3  BM  LS BM 27.8   CM CM 17.1  BM CM 13.9  BM CM 10.3  BM CM 6.1  BM CM 3.1  BM CM 5.0  BM CM 3.2  BM  LS 2.6 BM  CM  CM 9.4  BM CM 1.1  BM CM 1.1 BM 1.0 -a 66.7

∆G∗chg [Na+ ] [kJ·mol−1 ] -440.7 -441.7 -438.8 -439.9 -442.1 -445.3 -448.0 -448.6 -439.4 -440.3 -440.9 -414.5 -421.8 -439.2 -437.4 -393.3 -390.0 -370.9 -376.0 -378.3 -380.8 -412.3 -408.4 -413.6 -414.0 -100.5 -101.2 -97.9 -356.2 -354.1 -122.4 -121.0 -125.8 -125.7 0b

∆G∗chg [Na− ] [kJ·mol−1 ] -639.0 -641.0 -640.3 -271.8 -272.8 -273.5 -228.4 -229.5 -572.0 -577.6 -576.8 -561.3 -567.2 -239.1 -240.0 -502.1 -519.7 -266.4 -267.0 -634.6 -651.6 -190.3 -197.0 -349.2 -355.2 -223.0 -223.8 -222.1 -204.3 -208.7 -116.8 -119.1 -144.9 -144.0 0b

∆∆G∗lin [kJ·mol−1 ] -79.9 -76.9 -74.3 4.2 11.0 7.4 -20.1 -19.9 -51.2 -49.4 -55.8 -0.3 -0.1 -24.5 -24.6 14.9 14.8 10.5 10.4 -13.0 -13.0 2.5 2.5 13.6 13.6 15.0 19.1 5.9 -4.6 -4.5 85.8 85.8 96.0 96.0 0b

∆∆Gchg [kJ·mol−1 ] 278.2 276.2 275.9 -172.3 -180.3 -179.2 -199.5 -199.1 183.9 186.7 191.7 147.2 145.5 -175.6 -172.9 93.9 114.8 -115.0 -119.4 269.4 283.8 -224.6 -214.0 -78.0 -72.4 107.5 103.5 118.2 -147.3 -140.9 -91.3 -87.7 -77.0 -77.7 0b

(a) no simulations were performed; (b) all partial charges are zero. Charging free energies for the Na+ and Na− ions in

∆∆G∗chg [kJ·mol−1 ] 198.3 199.3 201.6 -168.1 -169.3 -171.8 -219.5 -219.1 132.7 137.3 135.9 146.9 145.4 -200.1 -197.5 108.8 129.6 -104.4 -109.0 256.4 270.8 -222.0 -211.5 -64.4 -58.8 122.5 122.6 124.1 -151.9 -145.4 -5.5 -1.9 19.0 18.4 0b

the physical solvents. The solvents considered are water (H2O), dimethylsulfoxide (DMS), dimethylacetamide (DAMD), methanol (MET), ethanol (ETL), acetone (PPN), isopropanol (2PL), pyridine (PYRI), acetic acid (ACA), ethylacetate (EAE), butanamine (BAN), chloroform (CHL), diethylether (DEE), toluene (TOL) and benzene (BZN). For each solvent, the reported entries are the relative dielectric permittivity ǫ′S and, for each electrostatic scheme, the charging free energies ∆G∗chg of Na+ and Na− along with the intrinsic asymmetry ∆∆G∗chg (Eq. 16) and the linear and quadratic contributions ∆∆G∗lin and ∆∆Gchg to this quantity (Eq. 20). The data is reported in more details in Suppl. Mat. Table S.5 and shown graphically in Figure 6.

63 ACS Paragon Plus Environment

The Journal of Physical Chemistry

Table 7:

64

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

Page 64 of 74

Solvent H2O DMS DAMD MET ETL PPN 2PL PYRI ACA EAE BAN CHL DEE TOL BZN CCL

ǫS 78.4 46.4 39.6 33.0 24.3 20.8 19.1 12.9 6.2 6.1 5.2 4.7 4.2 2.4 2.3 2.2

Scheme LS LS CM LS CM CM CM CM CM CM CM LS CM CM CM -a

A (exp.) 1.00 0.34 0.27 0.75 0.66 0.25 0.59 0.24 0.93 0.21 0.15 0.42 0.12 0.13 0.15 0.09

A (sim.) 1.00 0.40 0.35 0.88 0.83 0.37 0.73 0.38 0.94 0.28 0.50 0.32 0.30 0.11 0.14 0b

A∗ (sim.) 1.00 0.43 0.36 0.90 0.88 0.37 0.79 0.42 0.99 0.30 0.55 0.35 0.32 0.18 0.23 0b

B (exp.) 1.00 1.08 0.97 0.50 0.45 0.81 0.44 0.96 0.13 0.59 0.17 0.73 0.34 0.54 0.59 0.34

B (sim.) 1.00 1.10 1.09 1.03 1.03 1.07 1.00 0.94 0.93 1.03 1.05 0.27 0.88 0.41 0.43 0b

B ∗ (sim.) 1.00 1.00 1.02 1.00 0.94 1.00 0.89 0.84 0.86 0.94 0.94 0.23 0.81 0.28 0.29 0b

Σ (exp.) 1.00 -0.46 -0.57 0.66 0.48 -0.55 0.33 -0.63 1.20 -0.55 -0.51 -0.15 -0.65 -0.70 -0.68 -0.71

Σ (sim.) 1.00 -0.62 -0.72 0.66 0.53 -0.63 0.34 -0.41 0.97 -0.81 -0.28 0.39 -0.53 -0.33 -0.28 0b

Σ∗ (exp.) 1.00 -0.87 -0.96 0.77 0.57 -0.87 0.39 -1.03 1.59 -0.78 -0.56 -0.33 -0.80 -0.95 -0.94 -0.88

Σ∗ (sim.) 1.00 -0.85 -1.11 0.67 0.74 -1.01 0.55 -0.53 1.29 -1.12 -0.32 0.62 -0.77 -0.03 0.10 0b

Π (exp.) 1.22 1.01 0.94 0.90 0.85 0.86 0.82 0.92 0.80 0.74 0.53 0.88 0.60 0.69 0.72 0.59

Π (sim.) 1.00 0.66 0.63 0.94 0.90 0.63 0.83 0.59 0.94 0.56 0.71 0.30 0.52 0.22 0.25 0b

(a) no simulations were performed; (b) all partial charges are zero.

Acity and basity vs. asymmetry and polarity measures for the physical solvents considered in the present simulations. The solvents are water (H2O), dimethylsulfoxide (DMS), dimethylacetamide (DAMD), methanol (MET), ethanol (ETL), acetone (PPN), isopropanol (2PL), pyridine (PYRI), acetic acid (ACA), ethylacetate (EAE), butanamine (BAN), chloroform (CHL), diethylether (DEE), toluene (TOL) and benzene (BZN). The reported entries are the experimental permittivity ǫS , the electrostatic scheme used in the simulation, and for both experiment (“exp.”; see Table 1 and Eqs. 30-32) and simulation (“sim.”; see Table 6 and Eqs. 27, 28 and 33-35), the acity (A or A∗ ), the basity (B or B ∗ ), the asymmetry (Σ or Σ∗ ) and the polarity (Π).

ACS Paragon Plus Environment

Page 65 of 74

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 Journal of Physical Chemistry

Figure 1:

Schematic representation of the connection between molecular (charge or shape) asymmetry and solvation asymmetry, along with simplified solvent models evidencing these effects. The spherical anisotropically-charged model (SAM) is exempt of shape asymmetry but characterized by a tunable charge asymmetry (magnitude of the dipolar charges −qs and qs , distance b between the central and off-center sites). The non-spherical isotropically-charged model (NIM) is essentially exempt of charge asymmetry (dipolar charges −qs and qs , where the distance b between the unscaled and scaled sites is held constant) but characterized by a tunable shape asymmetry (scaling factor α applied to the Lennard-Jones radius of the scaled site). Intuitively, one would expect the solvents with a more exposed positive charge to be acitic and those with a more exposed negative charge to be basitic. This is correct in terms of the effective asymmetry, but not necessarily in terms of the intrinsic one (see Sections I and III.1).

65 ACS Paragon Plus Environment

The Journal of Physical Chemistry

Figure 2: +

-1

∆Gchg [kJ mol ]

0

-

+

Na /Na

-

+

K /K

F /F

-

+

+ -

-

I /I

Rb /Rb

-200 b[nm]

-400

0.00 0.06 0.08 0.10 0.12

-600 -800

-0.4 -0.3 -0.2 -0.100.1 0.2 0.3 0.4 -0.5 0.5

-1

∆∆Gchg * [kJ mol ]

600 400 200 0 -200 -400

-1

∆∆Gchg [kJ mol ]

-600 0

1000 750 500 250 0 -250 -500 -750 -1000

-1

0.5

-0.4 -0.3 -0.2 -0.100.1 0.2 0.3 0.4 -0.5 0.5

200 ∆∆G*lin [kJ mol ]

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

Page 66 of 74

100 0 -100 -200 -0.5

0

0.5

0

0 qs [e]

0.5

0.5

0

0.5

0

0.5

Charging free energies of the Na+ /Na− , K+ /K− , F+ /F− , Rb+ /Rb− and I+ /I− ion pairs in the SAMs. See Figure 1 for the definition of a SAM. The successive rows correspond to the average charging free energy (∆Gchg ; over the cation and the anion, Eq. 19), the intrinsic solvation asymmetry (∆∆G∗chg ; Eq. 16), the quadratic contribution to this quantity (∆∆Gchg ; effective asymmetry, Eq. 22), and the linear contribution to this quantity (∆∆G∗lin ; Eq. 21). The different curves correspond to different values of b (excentricity of the off-center site) and are displayed as a function of qs (charge of the off-center site). The simulations were performed for qs > 0 and the graphs are symmetrized with respect to qs = 0. The curves for b = 0.0 nm were not calculated as they are analytical (horizontal line at zero in all graphs). The different ion pairs (columns) are arranged in order of increasing ionic size (radii σ ˜ iw in Table 3). All calculations relied on the LS scheme. The data is reported numerically in Suppl. Mat. Table S.3.

66 ACS Paragon Plus Environment

Page 67 of 74

Figure 3: +

-

+

Na /Na

-

+

K /K

F /F

-

+

-

+ -

Rb /Rb

I /I

-100

-1

∆Gchg [kJ mol ]

0

-200 -300

α 0.4 0.5 0.6 0.7 1.0

-400 -500

00.1 0.2 0.3 0.4 0.5

400

-1

∆∆Gchg * [kJ mol ]

-600

200 0 -200 -400 00.1 0.2 0.3 0.4 0.5

-1

∆∆Gchg [kJ mol ]

400 200 0 -200 -400

00.1 0.2 0.3 0.4 0.5

400 -1

∆∆G*lin [kJ mol ]

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 Journal of Physical Chemistry

200 0 -200 -400 -0.5

0

0.5

0

0 qs [e]

0.5

0.5

0

0.5

0

0.5

Charging free energies of the Na+ /Na− , K+ /K− , F+ /F− , Rb+ /Rb− and I+ /I− ion pairs in the NIMs. See Figure 1 for the definition of a NIM, where b is set to 0.2 nm in all cases. The successive rows correspond to the average charging free energy (∆Gchg ; over the cation and the anion, Eq. 19), the intrinsic solvation asymmetry (∆∆G∗chg ; Eq. 16), the quadratic contribution to this quantity (∆∆Gchg ; effective asymmetry, Eq. 22), and the linear contribution to this quantity (∆∆G∗lin ; Eq. 21). The different curves correspond to different values of α (radius-scaling factor applied to the scaled site) and are displayed as a function of qs (charge of the scaled site). The simulations were performed for qs > 0 and the graphs are symmetrized with respect to qs = 0. The curves for α = 1.0 were not calculated as they are analytical for the three bottom graphs (horizontal line at zero), and would be very similar to the curve for α = 0.7 in the top graph. The different ion pairs (columns) are arranged in order of increasing ionic size (radii σ ˜ iw in Table 3). All calculations relied on the LS scheme. The data is reported numerically in Suppl. Mat. Table S.4.

67 ACS Paragon Plus Environment

The Journal of Physical Chemistry

Figure 4: (a)

(b) Solvent index

20

-

Na 0 Na + Na

15 10 5

SAMs

Solvent index

200 15 10 5 0

NIMs 0.1

0.2

0.3 Rrdf [nm]

0.4

0.5

(c) 20 Solvent index

I

I + I

10 5

-

0

15

SAMs

200 Solvent index

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

Page 68 of 74

15 10 5 0

NIMs 0.1

0.2

0.3 Rrdf [nm]

0.4

0.5

Average distance and radial orientation of the SAMs and NIMs around an anion (blue), an uncharged cavity (green) or a cation (red). Only SAMs and NIMs with qs > 0 are considered. (a) Schematic illustration of the distance between the ionic cavity and the less exposed (here negative, blue) or the more exposed (here positive, red) site of the solvent molecule. (b) Location Rrdf of the first peak in the ion-solvent radial distribution functions for the SAMs and NIMs around a sodium-sized cavity of negative (Na− ), vanishing (Na0 ) or positive (Na+ ) charge as a function of the solvent index (Tables 4 and 5; solvents in order of increasing permittivity). The solvent site is either the less exposed (filled circles) or the more exposed (empty circles) site. Note that for the NIMs with qs = 0.4 or 0.5 e and α = 0.6, it was not possible to identify a peak in the radial distribution function of the more exposed site around the uncharged cavity. (c) As in (b), but for an iodine-sized cavity. Note that for the NIM with qs = 0.4 e and α = 0.6, it was not possible to identify a peak in the radial distribution function of the more exposed site around the uncharged cavity. 68 ACS Paragon Plus Environment

Page 69 of 74

Figure 5:

0.1

2 (a)

0

1

-0.05 0

0.2

0.4

0.6

0.8

1

0.2 c(r)

0.5

gI1(r) gI2(r)

0

gI2(r)

(b)

3

0

2

-0.1

1 0

0.2

0.4

0.6

0.8

1

0.2

0 2

(c)

0.1

1.5

0

1

-0.1 -0.2

~

4

0.1

-0.2

cI1(r)

g(r)

-0.1

g(r)

1.5

g(r)

c(r)

0.05

c(r)

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 Journal of Physical Chemistry

0.5 0

0.2

0.4

0.6 r [nm]

0.8

1

0

Solvation structure around an uncharged sodium-sized cavity for an illustrative SAM, an illustrative NIM and the SPC water model. The curves shown are the orientational correlation function cI1 (r) as well as the radial distribution functions gI1 (r) and gI2 (r), where 1 is the less exposed site (central site for the SAM, unscaled site for the NIM, oxygen atom for SPC) and 2 is the more exposed site (off-center site for the SAM, scaled site for the NIM, hydrogen atoms for water). The function g˜I2 corresponds to the radial distribution function for site 2 in a hypothetical situation where site 1 is distributed according to gI1 , but the molecules are oriented entirely isotropically. (a) SAM with b = 0.12 nm and qs = 0.4 e. (b) NIM with α = 0.4 and qs = 0.2 e. (c) SPC water model. 186 Note the different scales on the left axis for the c function on the right axis for g functions. The curves are based on extended 5 ns simulations with configurations written to file every 0.1 ps. The equations defining the functions can be found in Ref. 160 (see Eqs. 51 and 52 therein).

69 ACS Paragon Plus Environment

The Journal of Physical Chemistry

-

0 -200 -400 -600

X + X

-1

∆∆Gchg * [kJ mol ]

-1

∆Gchg * [kJ mol ]

Figure 6:

200 0 -200 -1

∆∆Gchg [kJ mol ]

200 0

CCL

BZN

TOL

DEE

CHL

BAN

EAE

ACA

PYRI

2PL

PPN

ETL

MET

DAMD

100 50 0 -50 -100

DMS

-1

∆∆G*lin [kJ mol ]

-200

H2O

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

Page 70 of 74

Charging free energies for the Na+ /Na− ions in the physical solvents. The solvents considered are water (H2O), dimethylsulfoxide (DMS), dimethylacetamide (DAMD), methanol (MET), ethanol (ETL), acetone (PPN), isopropanol (2PL), pyridine (PYRI), acetic acid (ACA), ethylacetate (EAE), butanamine (BAN), chloroform (CHL), diethylether (DEE), toluene (TOL) and benzene (BZN). The successive rows correspond to the charging free energies (∆G∗chg ) of the anion (blue) and the cation (red), the intrinsic asymmetry (∆∆G∗chg ; Eq. 16), the quadratic contribution to this quantity (∆∆Gchg ; effective asymmetry, Eq. 22) and the linear contribution to this quantity (∆∆G∗lin ; Eq. 21). The different solvents are arranged in order of decreasing experimental relative dielectric permittivity ǫS . Calculations for rigid and flexible solvent models relied on the LS and CM electrostatic schemes, respectively. The data is reported numerically in Table 6 and Suppl. Mat. Table S.5.

70 ACS Paragon Plus Environment

Page 71 of 74

Figure 7:

B, B* (sim.)

1.5 PPN

ETL 1 ACA

MET EAE 2PL

BAN DEE

DAMD H2O DMS

A*, B* A, B

PYRI

0.5

CHL CCL

0

0

0.2

0.4

TOL BZN 0.6 B (exp.)

0.8

1

1.5 A, A* (sim.)

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 Journal of Physical Chemistry

1

2PL

ACA H2O

BAN PYRI PPN DMS DEE BZN DAMD CHL TOL EAE

0.5

0

ETL MET

0

CCL

0.2

0.4

0.6 A (exp.)

0.8

1

Correlation between calculated and experimentally-inferred acities and basities for the physical solvents considered in the present simulations. The calculated parameters (“sim.”; A∗ and B ∗ from Eq. 27 or A and B from Eq. 28) and the experimental values of Swain2 (“exp.”; A and B) are compared. The solvents are water (H2O), dimethylsulfoxide (DMS), dimethylacetamide (DAMD), methanol (MET), ethanol (ETL), acetone (PPN), isopropanol (2PL), pyridine (PYRI), acetic acid (ACA), ethylacetate (EAE), butanamine (BAN), chloroform (CHL), diethylether (DEE), toluene (TOL) and benzene (BZN). The dashed lines indicate the diagonals. The data is reported numerically in Table 7.

71 ACS Paragon Plus Environment

The Journal of Physical Chemistry

100 75 50 25 0

1 0

B

2 1 0

4

★ ★

★ ★

★ ★★ ★ ★

★ 2 ★ ★★

0 ✩

-2

4 *

A

2

Σ

εS

Figure 8:

Σ

2 0 -2

2 1

Π

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

Page 72 of 74

0 0

10

20 30 Solvent index

40

50

Acity and basity vs. asymmetry and polarity measures for the 55 solvents of the Swain set. The solvents are reported in order of increasing static relative dielectric permittivity ǫS (index and solvent name in Table 1). The six successive panels show (as points connected by lines) the permittivity ǫS , the Swain acity and basity parameters2 A and B, the derived asymmetry parameters Σ (Eqs. 30 and 33) and Σ∗ (Eqs. 32 and 34), and the derived polarity parameter Π (Eqs. 31 and 35). The magenta crosses depict corresponding simulation results for the 16 physical solvents (Table 2). The dashed lines indicate the zero point as a guide for the eye. The dotted line in the bottom panel shows 1 − ǫ−1 S . Filled stars label from left to right the solvents chloroform, acetic acid, trifluoroacetic acid, tert-butanol, butanol, isopropanol, propanol, ethanol, methanol, ethyleneglycol, formic acid, water and formamide. The empty star labels the solvent hexamethylphosphoric triamide. The data is reported numerically in Tables 1 and 7.

72 ACS Paragon Plus Environment

Page 73 of 74

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 Journal of Physical Chemistry

Figure 9:

Schematic illustration of three possible conventions for setting the zero of the electric potential when attempting to determine the absolute value of this potential within a solute cavity (large white circle) surrounded by liquid molecules (small white circles). In the R-convention, the liquid sample is assumed macroscopic and the zero is set outside the sample (blue area). In the P-convention, the zero is set as an average over the bulk liquid, where the averaging is carried out over both the interior and the exterior of the solvent molecules (red area). In the M-convention, the zero is set as an average over the bulk liquid, where the averaging is carried out over the exterior of the solvent molecules (green area). The electric potentials in the solute cavity corresponding to these three conventions are noted φR , φP and φM , respectively. In the context of atomistic simulations, the difference φP − φM is the exclusion potential ξS′ of the solvent model, an essentially arbitrary quantity related to the partial charge distribution within the solvent molecule (and the solvent density).

73 ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

TOC graphic

74 ACS Paragon Plus Environment

Page 74 of 74