Using Grand Canonical Monte Carlo Simulations to Understand the

Subscriber access provided by Northern Illinois University

Article

Using Grand Canonical Monte Carlo Simulations to Understand the Role of Interfacial Fluctuations on Solvation at the Water-Vapor Interface Kaustubh Rane, and Nico F. A. van der Vegt J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b05237 • Publication Date (Web): 17 Aug 2016 Downloaded from http://pubs.acs.org on August 23, 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 36

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

Using Grand Canonical Monte Carlo Simulations to Understand the Role of Interfacial Fluctuations on Solvation at the Water-Vapor Interface Kaustubh Rane∗,†,‡ and Nico F. A. van der Vegt† †Eduard-Zintl-Institut für Anorganische and Physikalische Chemie and Center of Smart Interfaces, Technische Universität Darmstadt, Alarich-Weiss-Strasse 10, 64287 Darmstadt, Germany ‡Current address: Chemical Engineering, Indian Institute of Technology, Gandhinagar, Palaj, Gujarat, India - 382355 E-mail: [email protected] Phone: 0091 8433596760

1

ACS Paragon Plus Environment


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

Abstract The present work investigates the effect of interfacial fluctuations (predominantly capillary wave-like fluctuations) on the solvation free energy (∆µ) of a monoatomic solute at the water-vapor interface. We introduce a grand-canonical-ensemble-based simulation approach that quantifies the contribution of interfacial fluctuations to ∆µ. This approach is used to understand how the above contribution depends on the strength of dispersive and electrostatic solute-water interactions at the temperature of 400 K. At this temperature, we observe that interfacial fluctuations do play a role in the variation of ∆µ with the strength of the electrostatic solute-water interaction. We also use grand canonical simulations to further investigate how interfacial fluctuations affect the propensity of the solute towards the water-vapor interface. To this end, we track a quantity called the interface potential (surface excess free energy) with the number of water molecules. With increasing number of water molecules, the liquid-vapor interface moves across a solute, which is kept at a fixed position in the simulation. Hence, the dependence of the interface potential on the number of waters models the process of moving the solute through the water-vapor interface. We analyze the change of the interface potential with the number of water molecules to explain that solute-induced changes in the interfacial fluctuations, like the pinning of capillary-wave-like undulations, do not play any role in the propensity of solutes towards water-vapor interfaces. The above analysis also shows that the dampening of interfacial fluctuations accompanies the adsorption of any solute at the liquid-vapor interface, irrespective of the chemical nature of the solute and solvent. However, such a correlation does not imply that dampening of fluctuations causes adsorption.

1. Introduction Understanding the propensity of ionic and non-ionic solutes towards liquid-vapor interfaces is important in many fields. The attraction or repulsion of ions to/from the air-water interface has important implications in atmospheric chemistry as well as in biochemistry where ions, 2


Page 2 of 36

Page 3 of 36

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


present in the hydration shell of proteins, affect their stability and function. 1 The water-vapor interface shares many features with the interface between water and extended hydrophobic surfaces, 2,3 which may be studied as models for nonpolar surface patches on proteins and for solvent-exposed areas of macromolecules in aqueous solution. 2,3 Therefore, studies on the affinity of small molecules or ions towards water-vapor interfaces can provide useful insights in the role of those molecules in biochemical processes and in affecting the water solubility of macromolecules. A detailed understanding of interfacial solvation thermodynamics is important in various other fields as well. For example, the prediction of the stability of thin liquid films requires to account for interfacial solvation. The adsorption of small molecules at the interfaces of such films influences their stability, and therefore, determines the properties of foams that contain such films. 4 The questions that we want to address in the present manuscript are the following: 1) How can we use molecular simulations to quantify the role of interfacial fluctuations to the free energy of solvation? This question is important from the computational perspective as well as from the perspective of conceptually isolating a specific molecular-scale phenomenon. 2) Do the solute-induced changes in the interfacial fluctuations affect the adsorption? This question emerges from the fact that strong correlations between adsorption and the dampening of interfacial fluctuations are frequently observed. 5–7 3) Are the above mentioned correlations specific to the chemical nature of water and some ions? Can they be predicted from the thermodynamics of a generic liquid-vapor interface? In order to address the above questions, we have divided the discussion of results into two parts. The first part addresses the first question. Here, the objective is to show whether interfacial fluctuations have a measurable effect on the free energy of solvation of a solute. The term ’interfacial fluctuations’ is used to denote the fluctuations in solvent density near the liquid-vapor interface. These mainly include the capillary-wave-like undulations that are characterized by long-ranged transverse correlations between density-fluctuations in the plane of the interface. 8,9 In the Appendix of our recent manuscript we discussed the rigorous

3



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

calculation of their contribution to the free energy of solvation from molecular simulations. 10 We concluded that those methods are impractical because of the large computational expenses and proposed an indirect method that uses an external field to dampen the interfacial fluctuations in canonical ensemble (CE) simulations. Interfacial solvation free energies were calculated based on Monte Carlo simulations with and without external field, and their difference was shown to quantify the effect of interfacial fluctations. The use of the canonical ensemble makes that approach suitable also for the molecular dynamics simulations. However, it is necessary to select the external field such that the structure of the solvent near the solute is not altered significantly. In the present manuscript, we propose a grand-canonicalensemble-based approach to quantify the effect of interfacial fluctuations on solvation. The main advantage of this approach is that the use of a grand canonical (GC) ensemble allows us to rigorously relate the computed quantity to the fluctuations in solvent densities. Secondly, the new method is more sensitive to the effect of interfacial fluctuations as compared to the previous approach. This enhanced sensitivity allows us to compare the role of interfacial fluctuations in the solvation of different solutes. Specifically, we compare the influence of electrostatic and dispersive interactions between solute and solvent. The second part addresses the last two questions. It discusses the adsorption of a single solute molecule at the interface between two saturated (liquid and vapor) phases studied in the grand canonical ensemble. We describe the analogy between the relevant free energy (interface potential) that is collected using grand canonical simulations and the potential of mean force (PMF) that is frequently used to study adsorption in the canonical ensemble simulations. It is shown that the profile of the interface potential vs the number of solvent molecules provides important insights about solvent reorganization induced by solute adsorption discussed in earlier work. 5 The analogy between the interface potential and the PMF therefore enables us to conclusively comment on the importance of solvent reorganization in the process of adsorption. In our previous work, 10 we explained that solute-induced changes (including those in solvent fluctuations) do not affect the free energy of solvation,

4


Page 4 of 36

Page 5 of 36

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


and consequently, have no effect on solute adsorption. This was done by referring to the potential distribution theorem (PDT) 11,12 and by describing the mutual cancellation of enthalpic and entropic components associated with solvent reorganization. 13–15 We refer to the first explanation as "informative" and to the latter one as "interpretative". In more general terms, the informative explanation can be summarized as follows: Given any two equilibrium states of a system in a particular ensemble, the difference between their free energies can be computed, in principle, from the information gathered at one of the states. This argument is related to the fact that free energies are state variables and the difference is independent of the path taken to go from one equilibrium state to another. In the process of solvation, the two states are 1) the solute molecule located outside the solvent and 2) the solute molecule located inside the solvent. The information gathered when the solute is located outside the solvent is sufficient for computing the free energy of solvation and therefore, the solute-induced changes are unimportant. The analogy between interface potential and PMF, discussed in section 6.3, helps us to relate the informative and interpretative explanations. Particularly, we show how thermodynamic quantities associated with solvent reorganization (related to the interpretative approach) depend on the path between equilibrium states, and following the informative explanation, should not influence the free energy difference. We further show that the strong correlation that is observed between the solute adsorption and solvent reorganization (like dampening of capillary waves) is independent of the nature of solute and solvent. The paper is organized as follows: Section 2 discusses the theory behind the grandcanonical-based strategy to quantify the contribution of interfacial fluctuations to the free energy of solvation. Section 3 explains the simulation strategy. Section 4 discusses the models that are used for solute and solvent. Section 5 provides the details about the simulations. Section 6 presents the results and discusses the important observations and section 7 concludes the paper.

5



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 36

2. Theory In this section we introduce a quantity that can be used to compare the contribution of interfacial fluctuations to the solvation free energy of different systems. In the grand canonical ensemble, interfacial fluctuations can be characterized by the response of the fluid density near the interface to a perturbation in the chemical potential. A higher magnitude of this response (χ) indicates stronger interfacial fluctuations. Below, we show that the response of the solvation free energy of a solute to a perturbation in the chemical potential of the solvent (see

d∆µ df

defined in equation (4)) is related to χ. Specifically,

d∆µ df

quantifies the coupling

between the solute and the interfacial fluctuations. Since the above quantity can be easily computed using the grand canonical simulations, it can serve as a useful measure to study the role of interfacial fluctuations in solvation thermodynamics. Note that such an indirect procedure becomes necessary because the direct calculation of the part of ∆µ that comes from interfacial fluctuations is computationally impractical. 10 We consider a system containing rigid solvent molecules and a solute molecule fixed at a particular location. The interaction energy between a solute and a solvent molecule is denoted by VΨ (R1 , R2 ; λ). Here, λ is the appropriate coupling parameter such that VΨ is a linear function of λ. The vector Ri contains translational and orientational coordinates of the molecule (i = 1 for the solute and i = 2 for the solvent). The interaction of an applied field with a solvent molecule at R is denoted by VF (R; f ), where f denotes the strength of the applied field. Note that the field does not act on the solute molecule. The grand potential of the system containing a solute fixed at R1 is denoted by Ω(λ; R1 ). Let ∆µ(λ; R1 ) denote the free energy difference associated with changing the coupling parameter to λ starting from a reference value λref . For a rigid solute, ∆µ is equal to the change in grand potential of the system. 16

∆µ(λ; R1 ) = Ω(λ; R1 ) − Ω(λref ; R1 )

6


(1)

Page 7 of 36

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


If the solvent molecules are also assumed to be rigid, then one can express the difference in grand potentials as follows:

Ω(λ; R1 ) − Ω(λref ; R1 ) =

Z

λ

dλ λref

Z

dR2 ρ(R2 ; R1 , λ)

∂VΨ (R1 , R2 ; λ) ∂λ

(2)

Here, ρ(R2 ; R1 , λ) denotes the density of solvent molecules when the coordinates of the solute are fixed at R1 and the coupling parameter is λ. In this work, VΨ is a linear function of λ. Let vΨ (R1 , R2 ) =

∂VΨ (R1 ,R2 ;λ) . ∂λ

Equations (1) and (2) can then be used to derive the

functional derivative of ∆µ with respect to the applied field VF : δ∆µ(λ; R1 ) = δVF (R3 ; f )

Z

λ

dλ λref

Z

dR2

δρ(R2 ; R1 , λ) vΨ (R1 , R2 ) δVF (R3 ; f )

(3)

2 ;R1 ,λ) Since the solute is fixed in our system, we can denote G(R2 , R3 ; R1 , λ) = −kT δρ(R δVF (R3 ;f )

where G is the solvent density-density correlation function. If we assume that Vf is a linear function of a parameter f , then Equation (3) can be integrated to derive the following expression:

d∆µ(λ; R1 ) = −β df Here, vF (R3 ) =

∂VF (R3 ;f ) ∂f

Z

λ

dλ λref

Z Z

dR2 dR3 G(R2 , R3 ; R1 , λ)vΨ (R1 , R2 )vF (R3 )

(4)

and β = 1/kT with k the Boltzmann constant and T the tempera-

ture. Notice that the right hand side of Equation (4) is only influenced by the density-density correlations of the solvent. Therefore, d∆µ/df quantifies the contribution of solvent fluctuations to ∆µ. We stress that the interpretation of the Equation (4) depends on the definition of λref because G(R2 , R3 ) can change significantly with λ. The definition of the external field vF in equation 4 governs the physical interpretation. It is known that fields directed perpendicular to the liquid-vapor interface dampen the capillary waves, and therefore, the interfacial fluctuations. 8,9 If vF corresponds to such a field, then d∆µ/df denotes the effect of dampening of interfacial fluctuations by vF . 10 Let fdamp be 7



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 36

the field-strength where capillary waves are almost completely absent with negligible effect on the other aspects of the interface. Equation (4) can then be integrated between f = 0 and fdamp to approximately obtain the contribution of interfacial fluctuations to ∆µ. Notice that this quantity is ∆∆µ that was introduced in our previous work 10 and is also discussed in the supporting document of the present manuscript. Herein, we are interested in a constant field vF which is independent of the solvent positions. In this case, the physical interpretation of

d∆µ df

is different as will be explained

below. First, such a field is meaningful only in the grand canonical ensemble. For example, the field vF (R3 ) = −1 can be interpreted as a shift in the chemical potential of the solvent from its value at liquid-vapor coexistence as µ = µsat + f . Notice that the negative sign can be understood from the definition of the grand potential. Then, the integral over G can be expressed in terms of the local compressibility 17 or susceptibility 9 χ as follows:

χ(R2 ; R1 , λ) =

Z

dR3 G(R2 , R3 ; R1 , λ)

(5)

From Equations (4) and (5) we get d∆µ(λ; R1 ) =β df We now explain how

d∆µ df

Z

λ

dλ λref

Z

dR2 χ(R2 ; R1 , λ)vΨ (R1 , R2 )

(6)

quantifies the contribution of interfacial fluctuations to ∆µ.

Solvent-density fluctuations contribute to ∆µ via the fluctuations in solute-solvent interaction energy. The instantaneous interaction energy (say, Ψ(R2 )) between a fixed solute and the solvent molecules located in an infinitesimally small region at R2 is given by ρ(R2 )vΨ (R2 ). The overline on the variables indicate that they are instantaneous quantities and the position of the solute is not included for clarity. We see that the fluctuations in Ψ are proportional to those in ρ. Since χ(R2 ) characterizes the density-fluctuations at R2 , R the integral dR2 χ(R2 )vΨ (R2 ) quantifies the solvent-density fluctuations to the free energy of solvation. This is exactly the right hand side of Equation (6). It is also known that density

8


Page 9 of 36

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


fluctuations are much greater near the liquid-vapor interface due to long-ranged transverse correlations. 9 Therefore, Equation (6) predicts that the magnitude of

d∆µ df

will be greater

for a solute positioned near the liquid-vapor interface than in the bulk liquid phase. Hence, the above derivative quantifies the role of interfacial fluctuations to ∆µ. This quantity can be numerically calculated from the grand canonical simulations via histogram reweighting techniques. 18 A positive shift (f ) in solvent chemical potential corresponds to an increase of the vapor pressure and a resulting dampening of density fluctuations, which are larger at the watervapor interface than in the bulk liquid. Therefore, the sign of d∆µ/df provides information on the direction in which the interfacial solvation free energy changes if surface capillary waves are dampened. The sign of d∆µ/df can be easily related to the microscopic phenomena. If λ > 0, then the sign is determined by the functional form of vΨ because χ > 0. If the integral is dominated by terms where vΨ > 0, then the right hand side of Equation (6) is positive. On the other hand, if terms containing vΨ < 0 dominate, then it is negative. Let us now consider different cases based on vΨ . Case1: vΨ > 0 d∆µ/df will be always positive. Case2: vΨ < 0 d∆µ/df will be always negative. Case3: vΨ > 0 for a certain range [Ra ,Rb ] and vΨ < 0 for [Rc ,Rd ]. d∆µ/df > 0 implies that solvent fluctuations corresponding to [Ra ,Rb ] are important. d∆µ/df < 0 implies that solvent fluctuations corresponding to [Rc ,Rd ] are important. Herein, we study different solutes in SPC/E water. 19 The SPC/E model contains a static electric dipole created by two positively charged hydrogen sites and a negatively charged oxygen site. Only the oxygen site interacts via the Lennard Jones (6-12) potential. Let us first consider the case of a Lennard Lones (LJ) solute where vΨ only depends on the translational coordinates of the water molecule. If the range of λ is such that the density

9



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

of solvent molecules at R2 where vΨ > 0 is negligible, then the system belongs to Case2. The system will always belong to Case 2 when λ < λref because the effective hard-core diameter of the solute decreases with decrease in λ. For λ > λref , the system will belong to Case2 if the effective hard-core diameter of the solute is not significantly greater than that corresponding to λref . Therefore, we expect that d∆µ/df < 0 for a LJ solute. The ionic solutes considered in this work interact with the SPC/E water molecules via a Lennard Jones potential as well as via electrostatic interactions. We use the coupling parameter λ to manipulate the electrostatic interaction only (see Equation (11)). We will use λref = 0 and λ > 0. Then, vΨ will be the electrostatic potential energy and the orientational coordinates of water molecules are important. Water molecules near an ion predominantly orient with their oxygen sites facing the cation and hydrogen sites facing the anion. Hence, vΨ < 0 for solvent molecules near an ion and d∆µ/df < 0 indicates that the solvent fluctuations near the ion are important. Moreover, this analysis indicates that dampening of surface capillary waves leads to a more favorable electrostatic ion-water free energy at the air-water interface. This may lead to enhanced ion adsorption provided that electrostatic interactions overcompensate repulsive, excluded volume interactions (Case1).

3. Simulation strategy All the simulations are performed in a parallelepiped box with the longest dimension along the z-axis (See Figure 1). Periodic boundary conditions are applied along the x- and ydirections. The box is bounded by a repulsive upper wall and a short-ranged attractive (sticky) lower wall along the z-axis. The details of these surfaces are provided in Section 4. This setup enables us to maintain a liquid-vapor interface that is parallel to the x-y plane and that is not affected by the two walls. Solvation is studied by positioning the solute in the center of the x-y plane at a desired distance (z) from the lower wall. Solvation at the liquid-vapor interface is sensitive to the saturation conditions. Therefore,

10


Page 10 of 36

Page 11 of 36

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


Figure 1: Molecular simulation setup that is used in this work. The filled violet circle represents the solute. The interaction potentials for the hard wall and the sticky wall are explained in the text. GC simulations use the chemical potential of the solvent µsat that corresponds to liquid-vapor coexistence. We determine µsat by tracking the interface potential V (N ) as a function of number of solvent molecules N for a system without solute. 20 The interface potential V is related to the Helmholtz free energy F by the following relation:

V (N ) = F (N ) − µsim N

(7)

Here, µsim is the chemical potential that is used to perform the GC simulations. The interface potential is always discussed relative to that at reference number of molecules Nref . In the grand canonical simulations V (N ) is related to the probability distribution over the number of molecules P (N ) as V (N ) = −kT ln[P (N )/P (Nref )]. If the liquid-vapor interface is not influenced by the two walls, V is expected to be constant at liquid-vapor coexistence.

11



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 36

This is because, at liquid-vapor coexistence, the relative amounts of liquid and vapor (and, therefore, the location of the liquid-vapor interface) can be changed without any change in free energy. We perform a GC simulation by using any initial guess of the chemical potential µo . This results in a monotonically increasing or decreasing profile of V (N ). We then employ histogram reweighting 18 to determine µsat that results in a constant value of V (N ). In order to determine ∆µ(λ; R1 ) in Equation (1) we use grand canonical expanded ensemble (GCEE) simulations. 21 The desired range of λ is discretized. GCEE simulations are then performed with moves that attempt to change λ along with insertion/deletion and displacement of solvent molecules. The number of molecules are allowed to vary over a small range in the GCEE simulations. This is done to ensure that the position of the macroscopic liquid-vapor interface is negligibly affected. This is important for studying the solvation of solute at a fixed position relative to the interface. If the number of molecules is not bounded, then the position of liquid-vapor interface can vary significantly at coexistence conditions. The probability distribution ΠEE (λ) is collected using the transition matrix scheme. 22 We also collect the probability distribution P (N ; λi , R1 ) for each value of λ. If λref denotes the coupling parameter of the reference, then

ΠEE (λi ; R1 ) ∆µ(λi ; R1 ) = −kT ln ΠEE (λref ; R1 )

(8)

The effect of perturbing the solvent chemical potential on ∆µ is evaluated by using the histogram reweighting technique. 18 The quantity d∆µ/df is expressed with respect to a suitable coupling parameter λref and is numerically estimated as follows: P P (N ; λi , R1 )eβN f d∆µ(λi ; R1 ) kT N ≈ − ln P βN f df f N P (N ; λref , R1 )e

(9)

Here, the distribution is reweighted to the new chemical potential µsat + f of the solvent, where f is chosen small compared to µsat . We also perform grand canonical simulations where the solute position and the solute-

12


Page 13 of 36

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


solvent coupling parameter are fixed. These simulations are used to track the interface potential V (N ) as the liquid-vapor interface shifts relative to the position of the fixed solute. The shift of the liquid-vapor interface in z-direction mimics the transfer of the solute from the vapor to the liquid. Here, the number of molecules are allowed to vary over a wider range than GCEE simulations. Results are also discussed for the expanded ensemble simulations performed in the canonical ensemble mode. The strategy employed for these simulations was described in our previous work. 10 Here, we apply an external field that acts in the direction perpendicular to the liquid-vapor interface. The quantity ∆∆µ is then used to denote the difference between the canonical version of ∆µ calculated with and without the external field.

4. Models Most calculations are performed with the SPC/E water model 19 as solvent. We also show results for a Lennard Jones (LJ) solvent. Two types of solutes are studied combined with the SPC/E water model. The first type (LJ) interacts with the oxygen site of the water molecule via the following potential:

VΨ (R1 , R2 ; λi ) = 4λi εO

"

σO | R 1 − R2 |

12

−

σO | R1 − R 2 |

6 #

(10)

Where R1 and R2 are the positions of the solute and the oxygen site of a water molecule, respectively. The coupling parameter for a particular bin i in the GCEE simulations is denoted by λi . εO and σO denote the LJ energy and distance parameters, respectively, for the oxygen site of the SPC/E model. In the rest of this manuscript, we will refer to the coupling parameter λi in Equation (10) by εslt−sol . The second type of solute (ionic) interacts with the sites of the water molecule via the following potential:

13



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

λ i q 1 q2 + 4εO VΨ (R1 , R2 ; λi ) = 4πǫo | R1 − R2 |

"

σO | R1 − R2 |

12

Page 14 of 36

−

σO | R 1 − R2 |

6 #

(11)

Where q1 and q2 denote the point charges on the solute and the interaction site of the water molecule, respectively. The permitivity of vacuum is denoted by ǫo and R2 denotes the position of the oxygen or hydrogen site of the SPC/E model. The parameter λi is used to manipulate the electrostatic interaction. In the rest of this manuscript, we will refer to λq1 as the effective charge q of the solute. The second term (LJ part) is absent for the interaction with the hydrogen site of the SPC/E model. The symbols are the same as described for Equation (10). We also show results for the chloride and iodide ions using the interaction parameters from the work of Horinek et al. 23 The studies employing the LJ solvent are performed using the "LJ", "HC" and "SC" solutes that are described by the Equations (7), (8) and (9) of our previous work. 10 All simulations with an ionic solute are performed with a uniform background charge to maintain the electroneutrality of the system. The appropriate corrections for the interaction with the background charge 24 are used in GCEE simulations. We expect that the presence of a background charge does not influence the main conclusions of this manuscript. The "sticky wall" shown in the Figure 1 interacts with the oxygen site of the water molecule via the following potential:

Vsw

    ∞     = −εsw (lsw −z) lsw       0

z≤0 0 < z ≤ lsw

(12)

z > lsw

Where z is the distance of the oxygen site of water molecule from the surface. We use εsw = 2.67kT and 2kT for simulations performed at T = 300 K and 400 K, respectively, and lsw = 16Å. The external field that is applied on the liquid-vapor interface in the canonical 14


Page 15 of 36

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


ensemble simulations has the same functional form as Equation (12). For the field, we use εsw−field = 0.1kT and 0.075kT for simulations performed at T = 300 K and 400 K, respectively, and lsw−field = 65.4Å. The "hard wall" in Figure 1 interacts with a large step potential that restricts the molecules in the simulation cell. For the simulations performed with the LJ solvent, the details of the above two walls are identical to those used in our earlier work. 10

5. Simulation details Simulations with water are performed in a simulation box with a square cross section of 33 × 33 Å2 and a height of 65 Å. The GC simulations with water are conducted at 400 K to avoid sampling difficulties associated with the insertion/deletion of molecules. Results from the canonical ensemble simulations are discussed for 300 K and 400 K. The number of water molecules range between 700 and 1450 in the GC simulations, and between 1050 and 1100 in the GCEE simulations. We do not show results for GCEE simulations with the LJ solvent. The GCEE simulations with the LJ solute are performed with εslt−sol ranging between 0.01 and 2.5. The difference between adjacent bins is 0.01. The GCEE simulations with the ionic solute are performed with the effective charge q of the solute ranging between −1.0 and 1.0 and the difference between adjacent bins is 0.01. In all simulations, the microstates were sampled using displacement and rotation of the solvent molecules. The insertion or deletion of molecules and the change in εslt−sol or q were used to sample macrostates in GC and GCEE simulations. The Lennard Jones and the electrostatic interactions were calculated using the Ewald summation with the real space cut-off of 16 Å. The uncertainties are calculated by performing four independent simulations. The GC simulations with the LJ solvent are performed at T = 0.9 in scaled units. Results are shown for the number of molecules ranging between 1800 and 2250. The remaining details of the simulations with the LJ solvent are described our previous work. 10

15



6. Results and discussion 6.1 Variation of ∆µ with εslt−sol and q at T = 400K

∆µ (kT)

0 -1 interface

a

-2

bulk liquid

-3 -4

εslt-sol 0

1

0.5

1.5

2

2.5

3

0

∆µ (kT)

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 16 of 36

-50

b

-100

-150 -1

0 q

-0.5

0.5

1

Figure 2: Relative solvation free energies in water for two positions of the solute at T = 400 K. (a) The variation of ∆µ with εslt−sol (q = 0). The values are plotted relative to ∆µ(εslt−sol = 1.0). (b) The variation of ∆µ with q (εslt−sol = 1.0). The values are plotted relative to ∆µ(q = 0.0). Black and red lines denote the results for solute fixed at z = 15 Å(bulk liquid) and 29 Å(interface), respectively. The straight and dashed lines denote results obtained using the GCEE and the canonical EE simulations, respectively. Uncertainties are shown at selected points. ∆µ is in units of kT and q is in units of the electronic charge. We first discuss the variation of ∆µ with ǫslt−sol and q for a solute positioned in the bulk liquid and near the liquid-vapor interface. Figure 2 shows the results for LJ and ionic solutes calculated using the GCEE and EE simulations. The results from the canonical and the grand canonical ensemble are expected to converge in the macroscopic limit. The slight disagreement between the data computed in these two ensembles is due to the limited size of the simulation cell and the slight difference in the position of the liquid-vapor interface for the two ensembles (the solvent density profile across the liquid-vapor interface is shown in 16


Page 17 of 36

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


Fig. S2 for the canonical ensemble calculations and in Fig.3 for the grand canonical ensemble simulations). For the LJ solute, ∆µ increases for small values of ǫslt−sol due to the increasing effective size of the repulsive cavity of the solute. This increase is less pronounced for the LJ solute at the liquid-vapor interface due to the lower density of surrounding water molecules compared to the bulk. For higher values of ǫslt−sol , ∆µ decreases monotonically because of the increase in the solute-solvent attractive interaction. For the ionic solute, increasing the magnitude of the charge q results in smaller ∆µ-values. For the anion, ∆µ decreases more rapidly with q than for the cation. This is expected because of the relative size of oxygen and hydrogen sites of the SPC/E water model. Due to the smaller size of the hydrogen site, the water molecules can approach an anion more closely than a cation, leading to a stronger attractive interaction. 24,25 We also notice that the change in ∆µ with q is much larger than that with εslt−sol . This is expected because the electrostatic interaction is longer ranged than the LJ interaction.

6.2 Variation of d∆µ/df with ǫslt−sol and q Figure 3 shows the change in ∆µ produced by perturbing the chemical potential of water by a small amount µ = µsat + f . The results are shown for f = 0.01kT , which is very small in magnitude as compared to µsat = −10.84kT . The density profiles in Figure 3c show that the effect of this shift on the position of the liquid-vapor interface is negligible. Therefore, the observed response of the solvation free energy, d∆µ/df , is due to perturbation of the local water density as explained earlier. For a LJ solute positioned in the bulk liquid, the magnitude of d∆µ/df is statistically insignificant for most values of εslt−sol . In order to understand this observation, we recollect that the value of χ, as given by Equation (5), is negligible in the bulk liquid as compared to the liquid-vapor interface. Therefore, the contribution coming from solvent fluctuations to ∆µ is negligible for the investigated range of εslt−sol . For the LJ solute positioned at interface, the magnitude of d∆µ/df is statistically significant for many values of εslt−sol . This shows that the interfacial fluctuations do have a 17



d∆µ/df

5 0

a

-5 0

0.5

1

2

1.5

2.5

d∆µ/df

εslt-sol 0

b

-10 -20

-1.0

0.0 q

-0.5

1.0

0.5

80

-3

ρ (nm )

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 36

60

c

40

40 20

20 0 0.0

0 24

10.0

20.0

26

30.0

28

40.0

30

32

50.0

z (Å)

Figure 3: Change in solvation free energies due to a small change in the chemical potential of water at T = 400 K calculated using GCEE simulations. The chemical potential of water is increased by 0.01kT . (a) The variation of d∆µ/df with εslt−sol (q = 0). The values are plotted relative to d∆µ/df at εslt−sol = 1.0. (b) The variation of d∆µ/df with q (εslt−sol = 1.0). The values are plotted relative to d∆µ/df at q = 0.0. Black circles and red squares denote the results for the solute fixed at z = 15 Å(bulk liquid) and 29 Å(interface), respectively. Uncertainties are shown at selected points. (c) Number density of oxygen sites of the water model. Green and blue straight lines denote the densities at the saturation chemical potential and after increasing it by 0.01kT , respectively. The inset shows the magnified profile near the liquid-vapor interface. small effect on ∆µ. Also, we find that d∆µ/df < 0 as predicted from the theory (Section 2). For the ionic solute in the bulk liquid, the magnitude of d∆µ/df is statistically insignificant. Hence, solvent density fluctuations have a negligible effect on the free energy of charging the solute in the bulk liquid. For the ionic solute at the liquid-vapor interface, the effect of q is stronger. The magnitude of d∆µ/df increases monotonically with that of

18


Page 19 of 36

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


q, which indicates that the contribution of interfacial fluctuations to ∆µ increases with the charge on the solute. Following the discussion in the theory section, d∆µ/df < 0 shows that the fluctuations of solvent molecules near the ion are important. The slightly lower magnitude of d∆µ/df for the anion as compared to the cation may be due to the closer approach of hydrogen sites of water towards the anion than that of oxygen sites towards the cation. However, we do not comment on this difference further because of the low precision of our data. The variation of d∆µ/df with εslt−sol or q provides additional physical insight. The positive perturbation in the chemical potential of water results in the increase in the average pressure of the system. Consequently, the fluctuations at the interface are damped to a greater extent than those in the bulk water because of the larger compressibility of the water-vapor interface. An observed dependence of d∆µ/df on a particular order parameter (εslt−sol or q) therefore indicates that the interfacial solvation free energy is affected by interfacial fluctuations. Coming to the LJ solute, the size of the repulsive core of the solute increases with the increase in εslt−sol between 0.01 and 0.25 (Also see the discussion of Figure 2). Hence, the slight increase of d∆µ/df for the above range of εslt−sol (Figure 3) implies that the interfacial solvation of the repulsive core of the LJ solute is more favored in presence than in the absence of interfacial fluctuations. For the ionic solute, the decrease in d∆µ/df with q indicates that interfacial fluctuations make solvation less favorable for stronger electrostatic interactions. In other words, dampening of those fluctuations causes electrostatic stabilization of the ion at the water-vapor interface. We note that the above interpretation is sensitive to the position of solute with respect to the interface and further investigation is necessary before deducing general trends. The objective of the present manuscript is to quantitatively compare the contribution of interfacial fluctuations to ∆µ of different solutes. We reserve the detailed discussion about whether the fluctuations favor or disfavor the solvation for a future manuscript. The above results are qualitatively consistent with those obtained using the canonical-

19



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

ensemble-based approach. 10 In Figure S2 of the supporting document, we show how ∆µ for the above solutes change with the application of an external field in the canonical ensemble simulations. We note that in these simulations the location of the liquid-vapor interface is slightly different than in the GC simulations. The observed change (∆∆µ) is negligible for the LJ solute positioned in bulk water and near the interface. For the ionic solute positioned near the interface, the magnitude of ∆∆µ is statistically significant for some values of q. Thus, both the canonical as well as the grand canonical ensemble-based approaches indicate that the interfacial fluctuations play an important role in the solvation of ions via electrostatic interactions. In the supporting document we also show results from the canonical ensemble simulations performed at 300 K (Figures S1 and S3). The calculation of d∆µ/df was not performed at this temperature due to sampling difficulties encountered during the grand canonical simulations. The values of ∆∆µ indicate that the interfacial fluctuations have negligible influence on the ∆µ for all the solutes.

6.3 Variation of V with N The earlier discussion focussed on the variation in the free energy of a solute that is fixed at a particular position. It showed that interfacial fluctuations affect the free energy of solvation predominantly via electrostatic interactions between solute and solvent. We now discuss whether the process of adsorption is indirectly affected by the solute-solvent interaction due to the solute-induced changes in the solvent. In order to understand the propensity of a solute towards the water-vapor interface, we compute the interface potential V (N ) as a function of number of water molecules N . The position of the solute and its interaction with water molecules is kept constant, and N is allowed to vary over a larger range than in GCEE simulations. The process is depicted in Figure 4. The solute is on the vapor side of the interface at small N , and on the liquid side for large N . The order parameter N can be approximately expressed either in terms 20


Page 20 of 36

Page 21 of 36

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


Figure 4: The analogy between the calculation of the interface potential V (N ) via grand canonical ensemble (GC) simulations and the calculation of the potential of mean force (PMF) via canonical ensemble (CE) simulations. The dark blue and light blue regions denote the liquid and vapor phases, respectively, and red circle denotes the solute. The upper and lower rows depict configurations from the GC and CE simulations, respectively. The two configurations in a particular column have solute at the same distance from the liquid-vapor interface. Here, li denotes the distance from the interface. The number of solvent molecules in the grand canonical ensemble is denoted by Ni . of the position of the liquid-vapor interface for a fixed solute, or in terms of the position of the solute with the fixed interface. Therefore, the interface potential V (N ) is analogous to the potential of mean force (PMF) that is computed with respect to the interface in canonical ensemble simulations. 5,6,23 (Figure 8). Notice that this analogy is possible only at liquid-vapor coexistence because, like the PMF, V should approach a constant value as the solute moves to the bulk phase. One advantage of tracking the interface potential is that the shape of the profile can be rigorously related to the molecular-scale phenomena. In the following discussion, we will use this aspect to explain the role of interfacial fluctuations in the adsorption of solutes at the liquid-vapor interface. Figure 5 shows V (N ) for different solutes. The upper horizontal axis shows Gibbs dividing surface, which is calculated by averaging the densities over the entire plane of the interface. We observe that the LJ solute is slightly more stable in the vapor than in the bulk liquid. We 21



l (Å) 20.9

27

33.2

39.4

0 LJ anion cation ClI-

-10

-20

V(kT)

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 36

-30

-40

-50

-60

800

1000

N

1200

1400

Figure 5: Effect of different fixed solutes on the interface potential of the water-vapor interface at T = 400 K. Black, red, green, blue and violet straight lines denote the results for the LJ solute, a model anion, a model cation, a chloride ion and an iodide ion, respectively. The LJ parameters for the LJ solute and the model ions are the same as that of the oxygen site of the SPC/E water model. The parameters of the chloride and iodide ions are taken from the work of Horinek et al. 23 The locations of the Gibbs dividing surface corresponding to different N are indicated on the upper axis. They are calculated by averaging the density over the entire cross section of the simulation cell and are statistically the same for different solutes. Each solute is fixed at z = 30 Å. The dashed maroon line indicates the value of N where the Gibbs dividing surface passes through the solute. Uncertainties are shown at selected points. V is in units of kT . also observe a a slight drop in the free energy on the vapor side of the water-vapor interface. This behavior has been observed and explained in previous studies. 26–28 All the ions are more stable in liquid water due to favorable energetic contributions from the electrostatic interactions. The behavior of iodide and chloride ions is consistent with that observed at high temperatures in an earlier simulation study. 6 We also observe a small drop in the free energy on the liquid side of the interface for the model cation. This indicates the propensity 22


Page 23 of 36

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


of the above type of ion to adsorb at the interface. The effect of the nature of ion-water interactions on such propensity has been discussed before. 5,6,23,26 We focus on the role that is attributed to the interfacial fluctuations by some of the above works. 5,7 Figure 3 shows that interfacial fluctuations do have a small effect on ∆µ. Since the propensity towards the interface depends on ∆µ, interfacial fluctuations will affect such propensity. However, as explained earlier, this is not due to the solute-induced changes in the fluctuations. The interpretation of the shape of V (N ) in Figure 5 helps us to explain the above point. In the following explanation, the location of the solute is taken to be its position with respect to the liquid-vapor interface. We only consider equilibrium states. In Figure 5, different locations of the solute (relative to the liquid-vapor interface) correspond to different N . The propensity of the solute to be at one location with respect to another is given by the difference between the interface potentials corresponding to those two locations. Notice that there may be numerous other ways of computing the difference. For example, via PMF calculations for different locations of the solute, 5,6,23,26 or by calculating ∆µ at different locations. 10 Thus, the propensity of a solute to favor one location with respect to the another is independent of the process of transferring the solute between those two locations. One can, in priciple, compute the free-energy difference between two locations of solute by using the information gathered at one of them (recollect the informative approach mentioned in the introduction). Let us consider an order parameter X that is used to track the relevant free energy F . Let X(R1 ) and X(R1′ ) denote the values of the order parameter at the two locations of the solute. Then, according to the above argument, the shape of the free energy profile between X(R1 ) and X(R1′ ) is unimportant for the propensity of the solute to be at a particular location. The same can be said for the derivatives (of any order) of F with respect to X at X(R1 ) and X(R1′ ). The effect of the molecular-scale phenomena associated with those derivatives should therefore cancel out. In Figure 5, F ≡ V and X ≡ N . Equation (7) can be used to show that

23



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

dV = µ(N ) − µsat dN

Page 24 of 36

(13)

where µsat is the chemical potential of solvent that corresponds to the liquid-vapor saturation conditions and µ(N ) is the free energy change of inserting a solvent molecule to a system containing N solvent molecules. Notice that dV /dN is the free energy change associated with transferring a molecule of solvent from the bulk liquid phase (at coexistence) to a system containing N molecules of solvent and a solute. This process is explained in Figure 6.

Figure 6: The physical interpretation of dV /dN . The arrow is used to denote the transfer of a single molecule of solvent from the bulk liquid at coexistence to a system containing the solute. The red circle denotes the solute. The solvent molecule favors the bulk saturated liquid when dV /dN > 0; it favors the vicinity of the solute when dV /dN < 0 . Since our system is at liquid-vapor coexistence and the liquid-vapor interface is not influenced by any external field, dV /dN = 0 in the absence of the solute. In the presence of the solute, addition of a solvent molecule will be disfavored if dV /dN > 0. That is, the solute causes unfavorable changes in the solvent. If dV /dN < 0, then the solute causes favorable changes in the solvent. It follows from the discussion in the previous paragraph that dV /dN , and consequently, the above solute-induced changes in the solvent, are unimportant for the propensity of a solute towards a particular location. For example, if we focus on the profile of the model cation in the Figure 5, there is a minimum on the liquid side of the interface. The model cation has a greater propensity to be at this location than another (say, one corresponding to N = 1440). At N = 1440, we observe that dV /dN > 0 and therefore, the solute causes unfavorable changes in the solvent. According to our argument, these changes 24


Page 25 of 36

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


should not influence the abovementioned propensity. The second derivative of V (N ) relates to fluctuations in the number of molecules: 29 d2 V kT = dN 2 hδN 2 i

(14)

where hδN 2 i is the fluctuation in the number of molecules of the solvent. Since the system is at liquid-vapor coexistence and the liquid-vapor interface is not influenced by any external field, hδN 2 i → ∞ in absence of the solute. This is due to the divergence of the local compressibility that is associated with the interfacial fluctuations. Therefore, cates that the solute decreases the interfacial fluctuations. However, as

d2 V dN 2

d2 V dN 2

> 0 indi-

is unimportant

for the propensity of solute towards the interface, solute-induced dampening of interfacial fluctuations should not play any role in interfacial adsorption. Taking the example of the model cation in the Figure 5, the minimum in the profile (the favored location of the solute) implies that

d2 V dN 2

> 0. Thus, the model cation dampens the interfacial fluctuations near its

most preferred location. We do not discuss cases where to unstable states in real systems. Notice that

dV dN

and

d2 V dN 2

d2 V dN 2

< 0 because these correspond

are related to the enthalpic and

entropic changes associated with the solute-induced reorganization of solvent (recollect the interpretative approach mentioned in the introduction). The above discussion thus relates two approaches of understanding solvation thermodynamics.

6.4 V (N ) and its derivatives with respect to N for the LJ solvent The noise that is present in V (N ) of water does not allow the calculation of derivatives with suitable precision. In order to explain their physical significance, we show V (N ), and

d2 V dN 2

dV dN

for the LJ solvent in Figure 7. The location of the Gibbs dividing surface is also

indicated. The types of solutes and their interaction parameters have been discussed in our previous work. 10 In the following discussion we will refer to the solutes corresponding to the black, red, green and blue profiles in Figure 7 as LJ2.5, HC2.5, SC2.5 and HC1, respectively.

25



l 23.3 0

24.8

26.1

27.4

28.7

a V

-10 -20 -30

dV/dN

0.00 -0.05

b 6 4

2

4

d V/dN (x 10 )

-0.10

2

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

Page 26 of 36

2

c

0 1800

1900

2000

2100

2200

N Figure 7: Interface potential of a LJ fluid in the presence of LJ, HC and SC solutes, and its derivatives with respect to N . Black and green symbols denote results for LJ and SC solutes, respectively. Red and blue symbols denote results for the HC solute. The interaction potential of the LJ, SC and HC solutes is given by the expressions (7), (8) and (9), respectively, of our previous work. 10 Each solute is fixed at z = 26σLJ . Black, red and green symbols correspond to εslt−sol = 2.5εLJ , whereas blue symbols correspond to εslt−sol = 1.0εLJ . σLJ and εLJ denote the LJ distance and energy parameters of the solvent, respectively. The locations of the Gibbs dividing surface corresponding to different N are shown on the upper axis. They are calculated by averaging the density over the entire cross section of the simulation cell and are statistically identical for different solutes. The dashed maroon line indicates the N -value where the Gibbs dividing surface passes through the solute. (a) The variation of V with N . Uncertainties are shown at selected points. NOTE: The results shown with dV with blue symbols are magnified by a factor of 10. (b) The variation of first derivative dN N . These are obtained by the numerical differentiation of V (N ). Results are shown at few d2 V points for clarity. (c) The variation of the second derivative dN 2 with N . These are obtained by first smoothening the profiles in "b" using the Savitzky Golay filter 30 and then performing d2 V the numerical differentiation. Results are shown for the range of N -values where dN 2 > 0. Few points are shown for clarity.

26


Page 27 of 36

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


Note that we have multiplied V of HC1 by a factor of 10 for clarity. The profiles of the interface potentials V (N ) show that HC1 will prefer a location close to the liquid-vapor interface, while the other solutes prefer to be in the bulk liquid phase. The profiles of the first derivative,

dV , dN

show that the LJ2.5, HC2.5 and SC2.5 solutes cause favorable changes in

the LJ solvent for all locations near the interface. That is, the addition of a solvent molecule is favorable in their presence. The profile of the interface potential (V (N )) demonstrates that this will always be the case when a solute prefers to be in the bulk liquid. The HC1 solute brings about favorable changes at some locations, and unfavorable changes at others. Such a behavior is always expected when a solute adsorbs at the liquid-vapor interface (minimum dV = 0) is also the one in V (N )). It is interesting to note that the most preferred position ( dN

where the solute has no net effect on the solvent structure (µ(N ) = µsat ), i.e. the solute is invisible from the perspective of the solvent. Coming to the second order derivatives, only positive values of

d2 V dN 2

are plotted because

negative values correspond to the configurations that are unstable in real systems. We observe that all the solutes decrease the interfacial fluctuations when positioned at certain locations on the liquid side of the interface. From the profile of the interface potential, one may notice that such a behavior is always expected when the solute prefers to be on the liquid side of the interface. The above discussion shows that the analysis of the derivatives of V (N ) with respect to N provides quantitative information about the solute-induced changes in the structure and the fluctuations of the solvent. The above analysis also shows that any type of correlation between the observed propensity of a solute towards a liquid-vapor interface and the observed changes in the solvent can be explained by referring to the thermodynamics of a generic liquid-vapor interface. Such a correlation is therefore independent of the chemical nature of the solute and the solvent. Some recent observations 5,6 regarding the adsorption of ions at the water-vapor interface can be interpreted in terms of the derivatives of V (N ). Otten et al. 5 observed, using canonical-ensemble simulations, that the free energy of a system containing a fixed iodide ion

27



V(N) (Fixed solute)

N (~ Position of interface) Nads N

PMF (Fixed interface)

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

Vapor

Liquid

Liquid

Vapor R

Rads

Position of solute Figure 8: The analogy between the PMF from canonical ensemble simulations and the interface potential V from grand canonical simulations for a generic system where the solute adsorbs near the interface. Note that the profiles are representative. The upper panel shows V as a function of the number of solvent molecules. This scenario is depicted in the upper panel of Figure 4. The change in position of the liquid-vapor interface is proportional to N . Here, Nads denotes the number of solvent molecules where V is minimum. The lower panel shows the PMF as a function of the position of the solute with respect to the fixed liquid-vapor interface. This scenario is depicted in the lower panel of Figure 4. Here, Rads denotes the position of the solute where the PMF is minimum. The red arrow denotes the shift in the interface (upper panel) or the transfer of the solute (lower panel) such that the solute is initially on the liquid side and finally near the interface. is at a minimum when the ion is located near the water-vapor interface. They also observed a minimum in the enthalpic component of the free energy. It was proposed that, when the iodide ion is transferred from the bulk liquid water to the interface, the water molecules surrounding the ion move to the bulk liquid. This process results in the net reduction of the energy of the system. Let us consider the same system in the grand canonical ensemble and the simulation setup similar to the one used for Figures 5 and 7. There exists a direct 28


Page 28 of 36

Page 29 of 36

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


analogy between the PMF and V as depicted in Figure 8. Then V will have a minimum at Nads corresponding to the location Rads of the ion where Otten et al. observe the minimum in free energy. Since

dV dN

= 0 at Nads , Equation (13) shows that µ(N ) = µsat . This implies

that the water molecules will equally favor to be near the ion at Rads and in the bulk liquid at liquid-vapor coexistence. However, when N is slightly greater than Nads (or the ion is away from the interface on the liquid side), it is expected that

dV dN

> 0. This implies that

the water molecules will prefer to be in the bulk liquid than near the ion. The process of dV dN

tending to 0 as N → Nads can also be interpreted as the transfer of water molecules

surrounding the ion to the bulk liquid. Thus, our analysis is consistent with the observation of the reference 5 and shows that similar reorganization is expected for any solute-solvent system where the solute adsorbs at the liquid-vapor interface. While this reorganization of the solvent does not contribute to the free energy of adsorption, it affects its enthalpic and entropic components. Previous canonical ensemble simulation results 5,6 also observe a strong correlation between the adsorption of iodide at the water-vapor interface and pinning of capillary-wave-like fluctuations. According to our analysis, this is expected because the minimum in V (N ) implies

d2 V dN 2

> 0. Ou and Patel 6 also observe an enhancement in the above fluctuations for

some locations of the iodide ion below the observed minimum in the potential of mean force (PMF) (See the Figures 1 and 7 of the reference 6). The authors also highlight a slight barrier in the PMF at about the same positions. Since such a barrier in PMF will correspond to the negative curvature of V (N ) in the grand canonical ensemble, we think that the enhanced fluctuations may characterize unstable states.

7. Conclusions The grand canonical ensemble simulation is a useful tool to study solvation at the liquidvapor interface. We employed the grand canonical version of expanded ensemble simulations

29



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

and histogram reweighting techniques to quantify the effect of interfacial fluctuations on the solvation free energy of different (ionic and non-ionic) solutes in water. The quantity d∆µ/df that is computed here is more rigorously related to the characteristics of the interfacial fluctuations than ∆∆µ obtained from the previously proposed canonical-ensemble-based approach. 10 Using this technique, we conclude that, at T = 400 K, the interfacial fluctuations play a non-negligible role in the variation of the solvation free energy with the strength of the electrostatic interaction between ionic solutes and the SPC/E model of water. We also used grand canonical simulations to study the affinity of different solutes towards the water-vapor interface. We studied model solutes and two real solutes (chloride and iodide ions). These results, along with those mentioned in the previous paragraph, indicate that the propensity of an ion towards a particular location near the water-vapor interface may be influenced by interfacial fluctuations. The interface potentials computed from the grand canonical simulations show that the above influence is due to direct ion-water interactions. The indirect effect of the ion, like pinning or dampening of fluctuations, does not affect the propensity of the ion towards the water-vapor interface. We explained the physical significance of the first and second-order derivatives of the interface potential with respect to the number of solvent molecules. These quantities provide important insights into the solute-induced changes in the solvent structure and solvent density fluctuations. Based on the analysis of these quantities, we conclude that solute-induced dampening of interfacial fluctuations will be always observed when solutes adsorb at the interface, irrespective of the chemical nature of solute and solvent. The general nature of such correlations can help in verifying the phenomena studied using molecular simulations. Finally, we note that all the results were presented for a solute fixed at a particular position relative to the interface. This was necessary because we wanted to study the interfacial solvation of solutes that need not adsorb at the interface. In real systems containing the adsorbed solute, the position of the solute is coupled to the interfacial fluctuations. 15 Therefore, fixing the position of solute is similar to applying an additional external field. This external

30


Page 30 of 36

Page 31 of 36

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


field will influence the interfacial fluctuations and therefore, the estimated role of interfacial fluctuations will be different from the case when the solute is free to move. It should be possible to employ our method with a ’free’ solute, provided that the solute adsorbs at the interface. Our calculations with a fixed solute show that an increase of the pressure (and a corresponding dampening of interfacial density fluctuations) causes electrostatic stabilization of the ion at the liquid-vapor interface. This conclusion may be tested experimentally by raising the pressure of an inert gas (e.g. air) while monitoring the vibrational spectrum of the air-water interface containing adsorbed ions.

Acknowledgement KR acknowledges the fellowship from the Alexander von Humboldt foundation. We have benefited from the discussion with Prof. Dor Ben-Amotz and the comments of the anonymous reviewers. We are thankful to Prof. Jeffrey R. Errington for sharing the computer programs developed in his group. The computational resources were partly provided by the high performance computing center of Technische Universität Darmstadt.

Supporting Information Available The supporting document shows the results obtained from the canonical ensemble simulations using the approach described in our earlier work. 10 Results are shown for the temperatures 300 K and 400 K. This material is available free of charge via the Internet at http://pubs.acs.org/.

31



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) Jungwirth, P.; Winter, B. Ions at aqueous interfaces: from water surface to hydrated proteins. Annu. Rev. Phys. Chem. 2008, 59, 343–366. (2) Patel, A. J.; Varilly, P.; Jamadagni, S. N.; Hagan, M. F.; Chandler, D.; Garde, S. Sitting at the edge: How biomolecules use hydrophobicity to tune their interactions and function. The Journal of Physical Chemistry B 2012, 116, 2498–2503. (3) Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of water near extended hydrophobic and hydrophilic surfaces. The Journal of Physical Chemistry B 2010, 114, 1632–1637. (4) Angarska, J.; Dimitrova, B.; Danov, K.; Kralchevsky, P.; Ananthapadmanabhan, K.; Lips, A. Detection of the hydrophobic surface force in foam films by measurements of the critical thickness of the film rupture. Langmuir 2004, 20, 1799–1806. (5) Otten, D. E.; Shaffer, P. R.; Geissler, P. L.; Saykally, R. J. Elucidating the mechanism of selective ion adsorption to the liquid water surface. Proceedings of the National Academy of Sciences 2012, 109, 701–705. (6) Ou, S.; Patel, S. Temperature Dependence and Energetics of Single Ions at the Aqueous Liquid–Vapor Interface. The Journal of Physical Chemistry B 2013, 117, 6512–6523. (7) Ou, S.; Cui, D.; Patel, S. Liquid–Vapor Interfacial Properties of Aqueous Solutions of Guanidinium and Methyl Guanidinium Chloride: Influence of Molecular Orientation on Interface Fluctuations. The Journal of Physical Chemistry B 2013, 117, 11719–11731. (8) Wertheim, M. Correlations in the liquid–vapor interface. The Journal of Chemical Physics 1976, 65, 2377–2381. (9) Stecki, J.; Toxvaerd, S. Susceptibilities of liquid–vapor interface of simple liquids. The Journal of chemical physics 1996, 105, 4191–4196.

32


Page 32 of 36

Page 33 of 36

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


(10) Rane, K.; van der Vegt, N. F. A. Understanding the influence of capillary waves on solvation at the liquid-vapor interface. The Journal of Chemical Physics 2016, 144, 114111. (11) Widom, B. Some topics in the theory of fluids. The Journal of Chemical Physics 1963, 39, 2808–2812. (12) Ben-Naim, A. A simple model for demonstrating the relation between solubility, hydrophobic interaction, and structural changes in the solvent. The Journal of Physical Chemistry 1978, 82, 874–885. (13) Yu, H.-A.; Karplus, M. A thermodynamic analysis of solvation. The Journal of chemical physics 1988, 89, 2366–2379. (14) Ben-Amotz, D.; Raineri, F. O.; Stell, G. Solvation thermodynamics: theory and applications. The Journal of Physical Chemistry B 2005, 109, 6866–6878. (15) Ben-Amotz, D. Interfacial solvation thermodynamics. Journal of Physics: Condensed Matter 2016, in press. (16) Henderson, J. In Fundamentals of inhomogeneous fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; Chapter 2, pp 23–84. (17) Evans, R.; Stewart, M. C. The local compressibility of liquids near non-adsorbing substrates: a useful measure of solvophobicity and hydrophobicity? Journal of Physics: Condensed Matter 2015, 27, 194111. (18) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo technique for studying phase transitions. Physical review letters 1988, 61, 2635. (19) Berendsen, H.; Grigera, J.; Straatsma, T. The missing term in effective pair potentials. Journal of Physical Chemistry 1987, 91, 6269–6271.

33



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

(20) Grzelak, E. M.; Errington, J. R. Calculation of interfacial properties via free-energybased molecular simulation: The influence of system size. The Journal of chemical physics 2010, 132, 224702. (21) Rane, K. S.; Murali, S.; Errington, J. R. Monte Carlo simulation methods for computing liquid–vapor saturation properties of model systems. Journal of Chemical Theory and Computation 2013, 9, 2552–2566. (22) Wang, J.-S.; Swendsen, R. H. Transition matrix Monte Carlo method. Journal of statistical physics 2002, 106, 245–285. (23) Horinek, D.; Herz, A.; Vrbka, L.; Sedlmeier, F.; Mamatkulov, S. I.; Netz, R. R. Specific ion adsorption at the air/water interface: The role of hydrophobic solvation. Chemical Physics Letters 2009, 479, 173 – 183. (24) Hummer, G.; Pratt, L. R.; Garcia, A. E. Free energy of ionic hydration. The Journal of Physical Chemistry 1996, 100, 1206–1215. (25) Gruziel, M.; Rudnicki, W. R.; Lesyng, B. Hydration free energy of a Model LennardJones solute particle: Microscopic Monte Carlo simulation studies, and interpretation based on mesoscopic models. The Journal of chemical physics 2008, 128, 064503. (26) Venkateshwaran, V.; Vembanur, S.; Garde, S. Water-mediated ion–ion interactions are enhanced at the water vapor–liquid interface. Proceedings of the National Academy of Sciences 2014, 111, 8729–8734. (27) Benjamin, I. Theory and computer simulations of solvation and chemical reactions at liquid interfaces. Accounts of chemical research 1995, 28, 233–239. (28) Vaikuntanathan, S.; Rotskoff, G.; Hudson, A.; Geissler, P. L. Necessity of capillary modes in a minimal model of nanoscale hydrophobic solvation. Proceedings of the National Academy of Sciences 2016, 201513659. 34


Page 34 of 36

Page 35 of 36

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


(29) Hansen, J.-P.; McDonald, I. R. Theory of simple liquids; Elsevier, 1990. (30) Savitzky, A.; Golay, M. J. Smoothing and differentiation of data by simplified least squares procedures. Analytical chemistry 1964, 36, 1627–1639.

35



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

36


Page 36 of 36

Using Grand Canonical Monte Carlo Simulations to Understand the

Recommend Documents