Calculation of the intrinsic solvation free energy profile of an ionic penetrant across a liquid-liquid interface with computer simulations.

We introduce the novel concept of an intrinsic free energy profile, allowing one to remove the artificial smearing caused by thermal capillary waves, which renders difficulties for the calculation of free energy profiles across fluid interfaces in computer simulations. We apply this concept to the problem of a chloride ion crossing the interface between water and 1,2-dichloroethane and show that the present approach is able to reveal several important features of the free energy profile which are not detected with the usual, nonintrinsic calculations. Thus, in contrast to the nonintrinsic profile, a free energy barrier is found at the aqueous side of the (intrinsic) interface, which is attributed to the formation of a water "finger" the ion pulls with itself upon approaching the organic phase. Further, by the presence of a nonsampled region, the intrinsic free energy profile clearly indicates the coextraction of the first hydration shell water molecules of the ion when entering the organic phase.

Introduction

The transport of ions as well as of uncharged penetrants across fluid (i.e., liquid–liquid or liquid–vapor) interfaces, in particular, between an aqueous and an apolar phase is of fundamental importance in various fields of physics, chemistry, and biology, from heterogeneous catalysis to ion extraction and from electrochemical applications to drug delivery.[1−5] In understanding the thermodynamic background of such transport processes, the determination of the solvation free energy profile of the penetrant across the interfaces is of great importance, because the free energy profile represents the thermodynamic driving force for the transport. In computer simulations, the density profile ρ(Z) of the penetrant molecule along the interface normal axis, Z, can, in principle, be easily converted to its solvation free energy profile A(Z) because A(Z) = −RT ln ρ(Z) (R and T being the gas constant and absolute temperature, respectively). Unfortunately, in the majority of systems of scientific interest the penetrant molecule is present in rather low concentrations, leading to a negligibly small probability of visiting certain regions of high free energy. This insufficient sampling can then make the above determination of the solvation free energy profile statistically unreliable. To circumvent this problem, several methods providing enhanced sampling of rare events,[6] such as harmonic[7] and adaptive umbrella sampling,[8] thermodynamic integration,[9] metadynamics,[10,11] the Widom test particle insertion method[12] and its cavity insertion variant,[13] or potential of mean force (PMF) calculation by, for example, constrained molecular dynamics,[14,15] have been proposed over the years. These methods of solvation free energy calculation have been applied for a number of ionic[16−33] and nonionic penetrants[13,34−43] at various fluid interfaces in the past decades, complementing theoretical predictions based on various continuum dielectric calculations.[44,45]

In calculating such solvation free energy profiles, however, one has to face the additional difficulty of taking into account the roughness of the interface. Fluid interfaces are corrugated by capillary waves which have the particularly harmful side effect of growing with the second power of the simulation box edge length, effectively smearing (in a system-size dependent way) any local property related to the interface.[46] Thus, any physically meaningful calculation of the solvation free energy profile of a penetrant, like any other kind of profile across such interfaces, requires the profile to be determined relatively to this molecularly rugged interface (usually referred to as the intrinsic interface) rather than to an external reference frame (or, equivalently, relative to an ideally flat planar surface, also called the nonintrinsic interface), as is usually done.

The first step toward the calculation of intrinsic solvation free energy profiles is clearly the determination of the intrinsic surface of a liquid phase. In other words, instead of simply determining the position of the Gibbs dividing surface along the surface normal axis of the basic simulation box, Z, the precise location of the surface along this axis Z has to be given at every {X,Y} point of the macroscopic plane of the surface. The determination of the intrinsic surface of a given condensed phase is usually coupled to the determination of the list of molecules that are located right at the boundary of this phase.

In early simulations of fluid interfaces, the intrinsic surface was approximated by dividing the system into several slabs along the macroscopic surface normal axis and determining the surface in each slab separately.[36,47−52] In further elaborating this approach, Jorge and Cordeiro proposed the use of a considerably finer grid than what is prescribed by the capillary wave theory and also determined the number of slabs required for convergence.[53] This method of grid-based intrinsic profiles (GIP) has been applied for a number of fluid interfaces.[53−55] The identification of the truly interfacial molecules (ITIM) approach,[56] also applied for a number of liquid–vapor[56−64] and liquid–liquid interfaces,[65−68] made one further step in this respect. In an ITIM analysis, a probe sphere of given radius, Rp, is moved along test lines from the bulk of the opposite phase to the surface of the phase to be analyzed. Once the probe sphere touches the first molecule of the phase of interest, it is stopped, and the touched molecule is marked as an interfacial one. The intrinsic surface itself is then approximated by the position of the interfacial molecules.[69] In the particular case when the distance between two neighboring test lines is equal to the probe sphere diameter, ITIM becomes equivalent to the GIP method of Jorge and Cordeiro. However, the ITIM method also allows to use a considerably denser grid, leading to a more precise detection of the intrinsic surface.[55] Another method, based on the relative distance between molecules of opposing phases, has been proposed by Chowdhary and Ladanyi for liquid–liquid interfaces.[70] This method was later extended by introducing an adjustable free parameter, and it was shown that this extended variant leads to results comparable with those of ITIM and GIP, but without applicability to vapor–liquid systems.[55] In their pioneering work, Chacón and Tarazona proposed a completely different approach.[46] Their intrinsic sampling method (ISM), applied later to a number of systems,[71−76] is based on a self-consistent algorithm in which the intrinsic surface is determined as the surface of minimum area that covers a set of appropriately chosen pivot atoms.[46,72] However elegant, this approach is too computationally demanding for direct application in PMF calculations. Several methods have also been proposed to identify the intrinsic surface for macroscopically nonplanar interfaces, such as the circular variance method of Mezei,[77] the coarse-graining approach of Wilard and Chandler,[78] or the generalized ITIM (GITIM) method based on α-shapes.[79] A recent comparison of various intrinsic surface determining methods revealed that ITIM provides an excellent compromise between the accuracy and the computational cost.[55]

Having determined the intrinsic surface location, the profile of any physical quantity relative to this surface can, in principle, be calculated. Several possible ways of calculating the density profiles of the components of the two phases were proposed and compared in terms of accuracy and computational cost in our recent study.[69] However, to the best of our knowledge, such an intrinsic analysis has never been applied to the calculation of free energy profiles. The reason for the lack of such studies is probably related to the computational cost required both by the intrinsic analysis and by the free energy calculation. Clearly, the calculation of solvation free energy profiles in simulations is a computationally far more demanding task than simply the generation of an equilibrium ensemble. The reason for this originates from the fact that in the canonical ensemble the Helmholtz free energy, A, is directly related to the full partition function, Q, through the equation A = −RT ln Q, and hence, in free energy calculations the entire configurational space rather than only its lowest energy domains has to be sampled. Indeed, the methods for free energy profile calculation require always considerably more computational effort than simply calculating density profiles. In the calculation of intrinsic free energy profiles, however, the computational burden is further exacerbated by the calculation of the intrinsic surface for—in principle—every sampled configuration during the simulation, which poses significant problems for the practical implementation of such methods.

In the present paper we propose a computationally feasible way of calculating the intrinsic solvation Helmholtz free energy profile of a single penetrant particle across fluid interfaces and apply it to the calculation of the intrinsic solvation free energy profile of a Cl– ion across the water–1,2-dichloroethane (DCE) liquid–liquid interface. The choice of this system is dictated by the fact that it is relatively well-studied in the literature: studies targeting the nonintrinsic free energy profile have been reported several times.[17,80,81] The method proposed here is based on the idea that the calculation of the intrinsic profile does not require the determination of the entire intrinsic surface, but just of a small portion of this surface that lies right behind the penetrant molecule. In the context of the ITIM method, this means that, instead of checking the full set of test lines in the entire basic box, only a few of them, located in the close vicinity of the penetrant, have to be considered. We “pull” the ion from the aqueous to the organic phase and use the intrinsic surface of the water phase as a reference, but the method is equally applicable to the reverse scenario (i.e., moving the ion from the organic to the water phase and/or using the DCE interface as a reference).

The paper is organized as follows. In section 2 details of the computer simulations, intrinsic surface analysis, and free energy calculations are given. In section 3 the resulting intrinsic solvation free energy profile is presented and compared to the corresponding nonintrinsic one, and the physical meaning of the obtained results is discussed in detail. Finally, in section 4 the main conclusions of this study are summarized.

Computational Methods

Molecular Dynamics Simulations

Molecular dynamics simulations of the water–DCE liquid–liquid interfacial system containing one single chloride ion at different positions have been performed in the canonical (N,V,T) ensemble at 298 K using the GROMACS 3.3.2 simulation program package.[82] The lengths of the X, Y, and Z edges of the rectangular basic simulation box (Z being perpendicular to the macroscopic plane of the interface) have been set to 50, 50, and 104 Å, respectively. The system consisted of 4000 water and 1014 DCE molecules.

Water molecules have been modeled using the TIP4P potential,[83] whereas standard OPLS force field parameters have been used for the Cl– ion and for the DCE molecules (see Table 1).[84] All bond lengths and bond angles have been kept fixed in the simulations, while torsional flexibility of the DCE molecule around its C–C bond has been allowed (Table 2). The CH2 groups of the DCE molecules have been treated as united atoms. The total potential energy of the system is calculated as the sum of the pair interaction energies of all molecule pairs. The interaction energy of two molecules has been calculated as the sum of the Lennard–Jones interactions acting between the sites and of the Coulomb interactions acting between the fractional charges. Bond lengths and bond angles of the DCE and water molecules have been kept fixed by means of the LINCS[85] and SETTLE[86] algorithms, respectively. The temperature of the system has been controlled using the Nosé–Hoover thermostat.[87,88] All dispersion interactions have been truncated to zero beyond the center–center cutoff distance of 9.0 Å. The electrostatic interaction has been calculated using the smooth particle mesh Ewald method[89] with a real space cutoff of 9.0 Å, a mesh grid with a spacing of 1.2 Å, and a spline order of 4. The charge of the Cl– ion has been compensated by a uniformly distributed positive background charge. An integration time step of 1 fs has been used in all simulations.

Table 1

Interaction Parameters of the Molecular Models Useda

molecule	interaction site	σ (Å)	ε (kJ/mol)	q (e)
waterb	Ow	3.154	0.649	0.000
	Hw	0.000	0.000	0.520
	Mwc	0.000	0.000	–1.040
Cld	Cl–	3.550	1.046	–1.000
DCEd	CH2	3.800	0.494	0.227
	Cl	3.400	1.255	–0.227

Lennard–Jones parameters correspond to the OPLS model, ref (84), fractional charges are taken from ref (48).

TIP4P model, ref (83).

Nonatomic interaction site, placed along the H–O–H bisector 0.15 Å away from the O atom toward the hydrogens.

OPLS model, ref (84).

Table 2

Geometry Parameters of the Molecular Models Used in the Simulations

molecule	bond	length (Å)	bond angle	angle (deg)
watera	O–H	0.9572	H–O–H	104.52
DCEb	CH2–CH2	1.53	Cl–CH2–CH2	108.2
	CH2–Cl	1.79

TIP4P model, ref (83).

OPLS model, ref (84).

Calculation of the Distance from the Intrinsic Surface

Since the main aim of this study is to calculate a free energy profile which is free from the artifacts due to the presence of capillary waves, a proper, local reference frame has to be chosen. In other words, the solvation free energy has to be calculated as a function of the position of the penetrant with respect to the local, molecularly rugged, intrinsic interface. ITIM[56] has been proven to be one of the most efficient methods for this task.[55] However, the full ITIM method has turned out to be computationally too costly to analyze the large number of trajectories required for the PMF calculation in a reasonable time. Thus, before applying it to the calculation of the free energy profile, the algorithm had to be simplified to reduce its computational cost. Since it is only the exact position of the penetrant, relative to the intrinsic surface, that has to be determined, it is sufficient to know the position of the surface only in the neighborhood of the point in the XY plane where the penetrant is located, and not in the entire simulation box. Therefore, instead of using a uniformly distributed fine grid of test lines across the entire XY plane as in the original ITIM method, we used only a reduced set of test lines, determined on the fly for every frame during the ITIM analysis from their lateral distance (i.e., the distance in the XY plane of the simulation box) from the ion. More precisely, only test lines whose lateral distance from the ion is smaller than N × dgrid are considered, whereas all the other test lines are disregarded from the analysis. Here, N is a conveniently chosen small integer leading to the use of N2 test lines, and dgrid is the grid spacing. Performing the ITIM analysis only on this set of test lines leads to significant savings in computer time. In the present study we have used a probe sphere radius of 1.25 Å, dgrid = 0.5 Å (in accordance with the suggestion of Jorge et al.[55]) and N = 4.

The last step needed to compute the position of the penetrant relative to the intrinsic surface is to interpolate the position of the surface (which is so far known only at a discrete number of points, corresponding to the positions of the detected interfacial atoms) to the XY projection of the penetrant. To do this task, several methods can be used,[69] out of which we have tested two, namely, the method of lifted Voronoi polygons[90] and the triangular interpolation.[69] These algorithms have been described and compared thoroughly in the study of Jorge et al.;[69] thus, here we provide only a brief explanation of the differences between them (see Figure 1c and e of ref (69)). In the lifted Voronoi polygons method, the Z position of the intrinsic surface is approximated by that of the surface molecule with smallest lateral distance from the ion.[90] In the triangular interpolation, on the other hand, the Z position of the intrinsic surface is obtained by linear interpolation between those of the three nearest interfacial molecules enclosing the XY projection of the ion.[69]

Reconstruction of the Intrinsic Potential of Mean Force

To compute the intrinsic Helmholtz free energy profile we have applied the constrained molecular dynamics approach.[14,15] The free energy profile is obtained by integrating the mean force required to keep the penetrant at predefined positions along a reaction coordinate. In our present intrinsic variant the reaction coordinate is the distance zintr of the penetrant from the intrinsic surface along the macroscopic surface normal axis. As a starting point for the series of constrained simulations, an equilibrated configuration of the ion-free water–DCE interfacial system has been taken from our previous work,[67] and a randomly chosen water molecule, situated well inside the bulk aqueous phase, has been replaced by the Cl– ion. The energy of the system has been minimized, and a 1 ns long relaxation run has been performed.

To create the starting configurations required for each subsequent run, the Cl– ion has been moved manually, in increments of 1 Å, from the position it occupied at the end of the previous run toward the bulk organic phase along the macroscopic interface normal axis Z. After each of such moves, the energy of the system has been minimized, and the system has been equilibrated for 100 ps. After the equilibration, an additional 500 ps production run has been carried out, during which configurations have been saved every 5 ps for the subsequent analysis.

For each saved configuration, the force FZ(t) required to keep the Z component of the ion position fixed with respect to the global coordinate frame has been recorded. Then, with the intrinsic ion–interface distance, zintr(t), provided by the ITIM analysis, one can readily obtain the constraint force as a function of the intrinsic distance, F(zintr). The values F(zintr) are then rebinned using the same bin width as the distance used to displace the ions, namely, 1 Å, therefore providing a conditional-constraint averaging ⟨F(zintr)⟩. The integration of ⟨F(zintr)⟩ along the reaction coordinate zintr provides the intrinsic Helmholtz free energy profile. Since zintr depends linearly on the absolute position of the ion along the macroscopic Z axis, no Fixman correction[91] is required.

Cluster Analysis

When applying the above-described procedure one has to be aware that the first hydration shell is usually coextracted with the ion as it moves into the organic phase. The coextraction of water by ions has been observed both by experimental and by theoretical methods. The group of Kusakabe[92,93] have determined the number of waters coextracted with halide anions, among them Cl–, in different organic solvents including 1,2-DCE. They have shown that each ion brings a certain number of water molecules to the organic phase regardless of the chemical nature of the organic phase. A wide range of theoretical works in the field also demonstrates the presence of coextracted waters in the organic phase. Benjamin and co-workers mention the problem of coextraction and water finger formation for instance in their work from 1995.[49] In a later paper Benjamin and Rose address the question of how the number of coextracted waters influence the calculated free energy of solvation.[94] Their results showed that indeed it is only possible to reproduce quantitatively the experimental free energies using a hydration number equal to the experimental one, whereas for the pure ion they obtain results that largely overestimate the experimental free energy difference.

In the present set of simulations the number of coextracted water molecules in the hydration shell of the Cl– ion turned out to be 20. To check how reproducible this value is, we have also performed a set of 100 out-of-equilibrium runs in which the ion was pulled by a constant force through the interface. (Unfortunately, repeating the full set of the original, quasi-equilibrium simulations several times turned out to be a computationally far too extensive task.) The number of the coextracted water molecules has been found to be 19 ± 3, in an excellent agreement with the 20 waters coextracted in our quasi-equilibrium set of runs.

The coextracted water molecules in the hydration shell of the Cl– ion should clearly not be considered for the intrinsic surface calculation once the hydrated ion exits from the aqueous phase. Otherwise, the ITIM method will misidentify these molecules as interfacial ones as they will stop the probe sphere before it reaches the region of the real interface. In this case, the ion will seemingly never cross the interface. This situation is depicted in the last panel of Figure 1. To overcome this problem, the two situations with the hydration shell attached to and detached from the aqueous phase (see panels b and c of Figure 1) have to be clearly distinguished.

Figure 1

Snaphosts of the system as taken from our simulations, with the penetrant at different intrinsic distances. Water O and H atoms are marked by red and white colors, respectively, while the penetrant Cl– ion is represented by a green sphere. (a) The Cl– ion is in the aqueous phase, close to the surface. The start of the formation of a water finger can be observed. (b) The Cl– ion is at the intrinsic water surface and well beyond the Gibbs dividing surface, situated at the tip of a clear water finger it pulled out of the surface. The hydration shell of the ion is still attached to the aqueous phase. (c) The hydrated Cl– ion is already in the organic phase, located farther from the intrinsic water surface than its radius. The hydration shell of the ion is already detached from the aqueous surface.

To make this distinction, we have used the method recently proposed by Pártay et al.[68] Thus, we define the aqueous phase to be the largest hydrogen bonding cluster of the water molecules in the system, and the hydration shell molecules as those water molecules the O atom of which is closer to the ion than 3.9 Å, i.e., the position of the first minimum of the corresponding radial distribution function. Two water molecules are defined to belong to the same cluster if they are connected via a chain of hydrogen bonding pairs, and two water molecules are regarded to be hydrogen bonded to each other if the distance of their O atoms and closest O–H pair are less than 3.35 Å and 2.45 Å, i.e., the first minimum position of the corresponding radial distribution functions, respectively. Finally, the hydrated ion is considered to be in the aqueous phase only if at least one of its hydration shell molecules belongs to the largest water cluster.

It turned out that in those cases when the size of the cluster to which a hydration shell water molecule belongs was at least three times the size of the hydration shell itself, this cluster was always also the largest cluster in the system. This allows reducing the computational costs of the procedure by stopping the algorithm when this critical cluster size is reached. However, the cluster analysis still turned out to be computationally as expensive as the MD simulation itself, and nearly as expensive as the ITIM analysis, which is the slowest step of the entire procedure. Hence, performing cluster analysis added an extra 30% to the computational cost. The computational costs of various steps of the algorithm, expressed in percentage of the total time needed to perform the entire procedure, are listed in Table 3.

Table 3

Computational Cost of the Subprocedures Constituting the Total Protocol of Reconstructing the Intrinsic Free Energy Profile, Expressed in Percent of the Total Time Needed To Run the Subroutine for One Single Configuration

procedure	time demand
MD simulation	28%
cluster analysis	27%
ITIM analysis	42%
triangular interpolation	1%
Voronoi method	1%
merging	1%
reslabbing and integration	1%

Solvation Free Energy Profiles

Nonintrinsic Solvation Free Energy Profile

The free energy profile of the chloride ion across the water–DCE interface is shown in Figure 2. For comparison, the density profiles of the entire aqueous and organic phases as well as of their first three ITIM layers are also shown as obtained in the ion-free reference system.[67] The complete free energy profile, interpreted thus together with the density profiles of the system, provides us with more detailed information. As seen in Figure 2, the free energy profile can be divided into three main regions, namely, a more or less constant part in the bulk aqueous phase (region I), a smoothly increasing part starting from the subsurface region of the aqueous phase (region II), which eventually turns smoothly into a plateau in the subsurface and bulk regions of the organic phase (region III). Region I corresponds to the ion being solvated in the bulk aqueous phase, where the opposite phase is far enough to exert negligible forces on it. Region II, i.e., the monotonically increasing part of the profile spans through the interfacial region, i.e., where the density of both water and DCE changes substantially. The free energy increase in this region is due to the progressively increasing net force the ion is subjected to upon approaching the organic phase. However, this monotonic increase of the profile reflects also the effect of the averaging over the fluctuating water surface, and hence no clear conclusion on the microscopic details of the energetics of the ion in the transition region can be drawn from this profile. This part of the profile turns smoothly, without exhibiting any other feature, into the plateau of region III, where the average free energy value of 55 kJ/mol corresponds to the work required to bring the hydrated ion from water to 1,2-dichloroethane. This value turns out to be comparable with the results of previous computational studies, which report values ranging between 50 and 60 kJ/mol for similar systems.[49,81,94] Furthermore, it is also in good agreement with the corresponding experimental value of 51.9 kJ/mol.[95]

Figure 2

Nonintrinsic free energy profile of the Cl– ion across the water–DCE interface (bottom panel). For reference, the (nonintrinsic) mass density profiles of the water and DCE molecules (filled and open squares, respectively) together with those of the first three molecular layers of both phases, detected by the ITIM method (solid lines), are also shown (top panel). The density profiles were obtained in the ion-free system.[67]

Intrinsic Free Energy Profile

We have reconstructed the intrinsic free energy profile of the chloride ion through the water/DCE interface using both the method of lifted Voronoi polygons and triangular interpolation (see section 2.2). The obtained profiles are shown in Figure 3 and compared to each other and to the nonintrinsic curve. The two intrinsic free energy profiles look very similar, apart from a small shift along the Z axis of about 2–3 Å, thus confirming the robustness of the approach. Similar shifts have been seen previously in various intrinsic density profiles calculated using these methods.[69] Since the method of lifted Voronoi polygons resulted in a slightly less noisy curve, in the following we only discuss the intrinsic free energy profile obtained this way. In Figure 3 it is seen that the free energy difference between the two bulk phases calculated from the intrinsic profile is in good agreement with that obtained in the nonintrinsic case. This fact shows that the bulk regions are not affected significantly by the instantaneous fluctuations of the interface caused by the capillary wave effect. However, upon approaching the interfacial region, the effect of the fluctuations becomes clearly visible in the intrinsic profile. Here, in contrast to the smooth and monotonic character of the nonintrinsic profile, the intrinsic curve exhibits several additional features, such as a local maximum, a local minimum, and an under-sampled region.

Figure 3

Comparison of the nonintrinsic (top panel) and intrinsic (bottom panel) free energy profiles of the penetrant Cl– ion at the water–DCE liquid–liquid interface. The intrinsic profiles obtained by the method of lifted Voronoi polygons and by using triangular interpolation are shown by full and open circles, respectively.

To have a better understanding of the features of the intrinsic free energy profile, we present it together with the intrinsic density profiles of the two phases (relative to the intrinsic water surface) in Figure 4. However, one has to keep in mind that the density profiles correspond to the ion-free system,[67] and hence any comparison can only be regarded as semiquantitative. The most apparent feature of the intrinsic profile is the presence of a clear maximum of about 63 kJ/mol around −5 Å. To understand the origin of this maximum, it has to be reminded that as the ion crosses the Gibbs dividing surface it pulls out a water “finger” with itself, being situated on the top of it.[17,80,81] In terms of an intrinsic analysis, however, the chloride ion inside such a water finger is still at the aqueous side of the (intrinsic) interface (see snapshot in Figure 1b). The free energy maximum can thus be attributed to the increased number of water–organic contacts (i.e., increase of the area of the interface between the two phases) and also that the chloride ion, located at the tip of the finger, has more unfavorable contacts with DCE than when it is located either in the bulk aqueous phase or, being surrounded by its hydration shell, in the bulk DCE phase. To quantify this effect we have calculated the average hydration number of the chloride ion at different values of its (signed) distance from the intrinsic surface, zintr, and found that around this free energy maximum it has, on average, 15% fewer contact water neighbors than in either of the two bulk phases.

Figure 4

Intrinsic free energy profile of the Cl– ion across the water/DCE interface, obtained by the method of lifted Voronoi polygons (bottom panel). For reference, the intrinsic mass density profiles of the water and DCE molecules (filled and open squares, respectively) relative to the water surface are also shown as obtained in the ion-free system[67] (top panel). The delta-like peak at the precise position of the water surface has been removed for clarity of visualization (for details, see ref (67)). The arrows show the profile at the intrinsic distance values corresponding to the three snapshots shown in Figure 1.

This maximum is followed by a descending part of the free energy profile, still at the aqueous side of the intrinsic surface. This reflects the progressive narrowing of the bottleneck of the water finger through which it is attached to the aqueous phase and the consequent reorganization toward the formation of a quasi-spherical isolated hydration shell around the chloride ion. This process corresponds to an increased shielding effect of water molecules around the ion and the decrease of the unfavorable water–DCE contacts, and eventually to a decrease of the free energy.

This descending part of the profile is followed by a nonsampled region between the zintr values of about −0.3 and 8 Å. The presence of this region indicates that no configuration characterized by such intrinsic distance values occurs during the course of the simulation. The reason for this can be understood considering the fact that the chloride ion is coextracted to the organic phase together with its first hydration shell,[96] and hence it cannot be closer to the intrinsic surface at its organic side than the radius of this shell. On the other hand, when the ion is located at the tip of a water finger, it can be arbitrarily close to the surface, and the value of −0.5 Å simply reflects the resolution of our profile. Stated differently, before its detachment the hydrated ion still belongs to the aqueous phase, although the chloride itself is rather close to the surface (being situated at the tip of the water finger), whereas just after the detachment the entire hydrated ion is regarded to be in the organic phase, evidently being farther from the intrinsic surface than its radius. It should be emphasized that the presence of this nonsampled region does not indicate any sudden displacement of the ion, rather, the breaking of the hydrogen bonds connecting its hydration shell to the aqueous phase. This means that two infinitesimally close configurations (i.e., when the hydrated ion is just before and after the detachment) correspond to two markedly different values of our reaction coordinate, zintr. As a consequence, the integration of the mean force cannot be used inside the nonsampled region, and hence the free energy profile beyond this region is, in principle, defined only up to an additive constant. However, since the configurations corresponding to the boundaries of this region are infinitesimally close, the continuity of the free energy allows us to find the physically meaningful value of this constant, requiring that the two boundaries of the nonsampled region correspond to the same free energy value.

Once the hydrated ion is detached from the aqueous phase but is still close enough to it, the attraction between the hydrated ion and the first molecular layer of the aqueous surface is still significant. Hence this layer experiences the competing attraction of the bulk aqueous phase and the hydrated ion. As a consequence, the first molecular layer is less strongly attached to the aqueous phase than either when the ion is in bulk water or when it is far enough from the interface at its organic side (see snapshots of Figure 1a and c). Thus, the progressive displacement of the chloride ion farther from the interface at its organic side corresponds to a decrease in the free energy (beyond the zintr value of about 10 Å), reflecting the relaxation of the first molecular layer of the aqueous phase.

Summary and Conclusions

The calculation of free energy profiles has become widespread during the last decades, thanks to the development of a series of methods which introduce an enhanced sampling in regions of the phase space that are otherwise rarely visited by standard simulation algorithms.[6−15] However, when calculating the free energy profile of a penetrant through an interface, using a simulation box-fixed reference frame to define the reaction coordinate, capillary waves can smear the obtained profile to such an extent that many of its important features are obscured. Here, we have introduced the concept of the intrinsic free energy profile, which is based on the adoption of a local reference frame to properly describe the position of the penetrant with respect to the interface. The intrinsic surface has been computed with the ITIM approach, and two different variants of the calculation of the penetrant position relative to the interface, namely, the methods of lifted Voronoi polygons[90] and triangular interpolation,[69] have been tested. We showed that a combination of constrained dynamics and conditional sampling can efficiently give access to the true free energy profile and reveal those features that are concealed by capillary waves in conventional calculations. Our approach requires only the constraint force and trajectory of a conventional potential of mean force calculation and can therefore be used to extract the intrinsic free energy profile from existing simulation data of the conventional, nonintrinsic free energy profile.

The importance of the intrinsic treatment in free energy profile calculations has clearly been demonstrated by the particular case of a chloride ion at the water–DCE liquid–liquid interface. Thus, the intrinsic treatment revealed several peculiar features of the free energy profile that were not seen in the conventional, nonintrinsic one. Among others, it turned out that the transfer process of the ion from one phase to the other is a process requiring activation energy, as the free energy profile exhibits a clear maximum at the aqueous side of the interface. Further, the fact that the ion coextracts its first hydration shell into the organic phase is also clearly indicated by the presence of a nonsampled region in the intrinsic profile, the width of which roughly corresponds to the radius of the hydrated ion. Our results suggest that, in the case of an apolar penetrant, which is not expected to coextract hydrating water molecules, no such nonsampled region of the intrinsic free energy profile should occur, and hence the calculation of such a profile would provide further support to the present results. Work in this direction is currently in progress.