Constrained Molecular Dynamic Simulation of the Potential Mean Force of Lithium Bromide Ion Pairs in Acetonitrile

Molecular dynamic simulations of Li+, and Br− ions in acetonitrile were carried out. The simulated structural properties were compared to experimental data. The solvent potentials of Li+-Br−, Li+-Li+, and Br−-Br− were evaluated using constrained molecular dynamics (CMD) simulations, to determine the solvent contribution to the total force acting on the solute and estimate the liquid arrangements according to the potential of mean force (PMF) values. The PMF of friction kernels and the passage across the Li+-Br− barrier was studied using the Grote–Hynes theory. The union-separation development happens in the polarization confining system.


Introduction
Understanding ionic association in polar solvents at the molecular level is vital to interpreting many chemical and biological processes [1]. Much research is being carried out based on molecular dynamics simulations of ionic liquids in polar diluents, such as methanol and water [2][3][4]. Constrained molecular dynamics (CMD) [5] simulations are performed while maintaining a fixed interior separation. In this way, the ion-ion potential of mean force (PMF) and coefficients of relative friction can be obtained to determine the association constants. In addition, such simulations can reveal the behavior of ions and adjacent molecules at the microscopic level. A characteristic of the liquids mentioned is that, unlike acetonitrile (MeCN), they form hydrogen bridges.
Equations for deriving PMF values for lithium bromide in acetonitrile through simulations have not yet been established, although theoretical calculations can be performed via the Ornstein-Zernike equation using a hypernetted-chain (HNC) approach [6]. Acetonitrile is a widely used organic solvent whose molecules do not form hydrogen bonds, unlike water (H 2 O) and methanol (MeOH). In addition, it has a relative permittivity (ε MeCN r = 36) significantly lower than that of water (ε H 2 O r = 78) and similar to that of methanol (ε MeOH r = 34). It can be beneficial to compare the results obtained with different solvents. The ion pair (Li + -Br − ) was chosen in this study because data on neutron diffraction [7] and dielectric relaxation [8] are available for comparison. In the last 20 years, a great deal of progress has been made in the study of ions in solution (both aqueous and non-, and both experimental and theoretical). For example, the importance of explicit polarizability in models to describe solvation of anions such as Br − is well established, and certainly new experimental data using neutron scattering or other methods are available. Metal salt solutions in polar aprotic solvents have a wide range of uses, including electrolytes in electrochemical devices (such as lithium batteries with high energy density, capacitors, and so on) and reagents in organometallic chemistry. The details of interactions between dissolved salt ions and organic solvent molecules influence the density, viscosity, volatility, ionic conductivity, and other essential properties of ionic solutions. As a result, ion-solvent interactions appear to be critical for furthering our understanding of ionic solutions. The interaction of solvents with positive and negative ions differs significantly, making it difficult to obtain a thorough knowledge of the anion solvation mechanism. Because of the huge size of anions, their intricate geometry (for polyatomic anions) and charge distribution, anion solvent interactions are weaker than cation solvent interactions. Second, in aprotic solvents, anion solvation is mostly accomplished through hydrogen bonding, which is absent in anionic solvents. For these reasons, it is commonly assumed that metal salt dissolving in polar aprotic solvents occurs solely due to the solvation of metal cations, whereas anions remain uncoordinated by the solvent molecules. Nevertheless, a thorough investigation of monoatomic halide anions using both theoretical and experimental methods has revealed that these notions are oversimplified, and that interactions of the halide anion with polar aprotic solvent molecules is possible. The weak halide-solvent interactions, on the other hand, make theoretical analyses of halide solvation difficult, necessitating a fine balance between ion-solvent and solvent-solvent interactions. For bromide, as well as aprotic and non-polar solvents, this trend becomes very pronounced. When acetonitrile is utilized as the solvent, solvation numbers in the range of 1-9 are discovered in the literature for this anion. Furthermore, the solvent molecules' molecular organization is not well understood [9][10][11].
The PMF values (W (r)) between particles with electrical charge in liquid perform a crucial part of the analysis of liquefied mixtures. CMD, as a computer simulation technique, is very useful for calculating W (r). CMD simulations can be used to calculate the mean force (MF) exerted on dissolved substance by solvent, whereas the separation between the solute is constrained. We used CMD and non-constrained (non-CMD) simulations to calculate the W(r) of Li + -Br − , Li + -Li + , and Br − -Br − and the dynamics of these particles with electrical charge.
This article is structured as follows. Details of our simulations and interaction potential models are given in Section 2. The ion-solvent structural properties obtained from the non-CMD simulations, the ion-ion MF along with PMF results obtained for all the simulations and the Li + -Br − pair friction kernel and the interconversion process are analyzed in Section 3. The conclusions of this paper are summarized in Section 4.

Materials and Methods
The PMF values for Li + -Br − , Li + -Li + , and Br − -Br − in acetonitrile have been evaluated. We carried out CMD simulations of systems built of two ions and 214 solvent molecules located in a cubic box with periodic boundary conditions. The volume of the cubic box (L = 26.58 Å) was selected to give a solvent density of 0.777 g/cm 3 and with a temperature of 298 K. Additionally, non-CMD simulations were performed for two systems with free ions: one with Li + and another with Br − . All code/programs-software used to perform the simulations were developed in-house. Mainly the simulations codes were developed using Fortran 77. The code included routines which evaluated the interaction potential, periodic conditions, integration algorithms, implementation of constraints (SHAKE procedure) and the mean force. See the Appendices A and B for simulation details.
We used an "acetonitrile rigid molecule pairwise additive potential model" [12,13] and performed calculations for the nitrogen (N), carbon (C) and methyl group (Me) interaction sites. Me is considered to be a unit atom model. The methyl-carbon and carbon-nitrogen distances were 1.458 Å and 1.157 Å, respectively, as shown in Figure 1. The dipole moment of each acetonitrile molecule is equal to 3.96 D. The ion-ion, ion-solvent and solvent-solvent interactions were evaluated using the following mathematical expression: in terms of the Lennard-Jones parameters, i and q i is the charge allocated to the place i.
We used the parameters for Li + and Br − derived by Chandrasekhar et al. [14] and Lybrand [15] from first principles computations of free ions in H 2 O. The values of the parameters are summarized in Table 1. By using the same ion's Lennard-Jones parameters for acetonitrile as for water, we are assuming that ε and σ transfer from one solvent to another.  [12]. b Values from [14]. c Values from [15]. Table 2 summarizes the range of the interionic separation for each pair and the total number of simulations performed. Each simulation consisted of a stabilization interval of 25 ps followed by a generation time of 50 ps. For certain interionic separations, CMD simulations of 100 ps were carried out to calculate statistical errors and friction kernels, and for detailed study of the configuration of the solvent over the ion pairs. In the case of acetonitrile, an integration step of 0.005 ps was used. We used the integration algorithm of Berendsen et al. [16] with the SHAKE technique [17] to fix the molecular distances. The Lennard-Jones interactions were cut off at L/2 and the long-range exchanges were calculated by the Ewald sum technique [18]. Radial distribution functions (g(r)). The number of particles in a volume dr centered at r, provided that one is at the origin is N(r + dr) = ρ bulk ·g(r)·dr; we know that (g(r)) depends solely on the modulus r, so we may better consider the number of particles contained in a spherical shell of radius r and r + ∆r, given by N(r, r + ∆r) = 4π·ρ bulk ·r 2 ·g(r)·∆r, from this formula g(r) can be calculated as, g(r) = lim ∆r→0 N(r,r+∆r) 4π·ρ bulk ·r 2 ·∆r . In a practical simulation, a finite ∆r is considered. At each time step of the numerical integration, all the pairs separated by a distance (r, r + ∆r) would be counted and added to a histogram, then a coarse-grained g(r) is obtained.
Finally, we carried out two unconstrained simulations of Li + and Br − , with 215 solvent molecules. The information obtained from these ion-solvent structure simulations is compared to the neutron diffraction data in Section 3. In addition, these simulations are representative of an infinitely diluted solution (zero concentration limit). According to Chandrasekhar and Jorgensen [19], the enthalpy of dissolution (∆H sol ) can be determined from these simulations by means of the following expression: where the last two terms are negligible under normal environmental conditions, and ∆U sol is the difference between the interaction energy of the infinite dilution system (U ∞ sol ) and the pure solvent (U * ): although these expressions are commonly used to determine ∆H sol from the simulations, we will use them in reverse to obtain an experimental estimate of U ∞ sol from the experimental results of ∆H sol and U*. Table 3 shows the U ∞ sol values corresponding to our simulations, which are compared against the experimental values obtained from the procedure described above. Table 3. Comparison between the experimental and simulation results regarding the interaction energy of a solution at infinite dilution U ∞ sol (kcal mol −1 ).

Ion-Acetonitrile Structural Properties
The structural properties obtained from the non-CMD simulations of Li + and Br − ions in acetonitrile are shown in Figure 2 (left axis) in the form of ion-acetonitrile radial distribution functions [g IMeCN (r) = g IN (r), g IC (r), g IMe (r)]. We similarly calculated the running coordination numbers (n(r)) which is giving by n(r) = 4πρ r 0 r 2 g(r )dr where ρ is the numeric solvent density. The function n(r) gives the mean number of molecules within a sphere of radius r centered on the ion. The coordination number is defined as the plateau value of n(r) at distances close to the first g(r) minimum ( Figure 2, right axis). Table 4 summarizes the positions of the first peaks of g IMeCN (r), and the coordination numbers obtained from our simulations and neutron diffraction experiments. The most probable distances, Li + -N and Li + -C, are well-reproduced. For Li + in acetonitrile, we obtain a coordination number slightly higher than the experimental one (5.5 instead of 3). Elastic neutron scattering measurements obtained by applying the first-order difference approach [22,23] yield experimental information about the arrangement of solvent molecules surrounding ions. This procedure establishes a difference function ∆G I , which is a combination of the partial radial distribution functions corresponding to each solvent atom: g IMe (r), g IC (r), g IN (r). The experimental solubility of LiBr in acetonitrile is 0.58 M. We calculated the ∆G MeCN Li + corresponding to this concentration from the g IMeCN (r) for Li + as follows: where the constants A-F depend on the concentration and scattering lengths. The contribution of the first two terms of (4) pertains to distances beyond the first hydration shell, A = B ∼ 0. The shape of the ∆G MeCN Li + function resulting from the simulations is in agreement with the experimental results ( Figure 3). We obtained the contribution to ∆G MeCN Li + for the last three terms of (4). The results are consistent with the interpretation regarding the principal ∆G MeCN Li + peaks. The discrepancies in the height of the ∆G MeCN Li + peaks can be attributed to the deficiency in the interaction potentials used, as well as to experimental inaccuracies and the hydrogen atoms of the methyl group not being taken into account. Nevertheless, nearly all of the computer simulation studies described for acetonitrile and acetonitrile solutions are based on the hypothesis of rigid molecular models that do not explicitly consider the hydrogen atoms of the methyl group. With this assumption, the computing time is substantially decreased, and longer simulations of larger systems can be carried out.  Table 4 summarizes the positions of the first peaks of g IMeCN (r), and the coordination numbers obtained from our simulations and neutron diffraction experiments. The most probable distances, Li + -N and Li + -C, are well-reproduced. For Li + in acetonitrile, we obtain a coordination number slightly higher than the experimental one (5.5 instead of 3). This failure of simulations can be attributed to the deficiency in the interaction potentials used and the hydrogen atoms of the methyl group not being considered explicit.

Mean Force Potentials
The force (∆F(t;r)) applied by N solvent molecules along the intermolecular axis of two A and B ions is indicated by: the resultant forces F AS (t; r) and F BS (t; r) acting on the solute caused by the molecular liquid, m A and m B are the masses of the solute,r is the unit vector alongside the AB path and µ = m A m B /(m A + m B ) is the reduced mass. This expression is calculated at every single time interval and averaged throughout the entire simulation. Where F d (r) is the direct solute force, and the resultant MF between ions is: where ∆F(r) ≡ ∆F(t; r) . W(r) was calculated by integrating the mean force: W(r 0 ) was calculated so the PMF values matched the macroscopic Coulomb potential at large distances (r 0 = 9.0 Å) and used the experimental acetonitrile relative permittivity (ε MeCN r = 36). Table 2 summarizes the range of interionic separation for each pair and the total number of simulations performed.
The statistical errors for W(r) and ∆F were analyzed as stated by the method outlined in the ref. [24]. Table 5 shows the statistical inefficiency (s ∆F(r) ), standard deviation (σ ∆F(r) ) and estimated statistical errors ( ∆F(r) ) values according to several interionic distances and simulation lengths. Simulations beyond 50 ps provided no further benefit. The estimated errors for the mean force potentials for the various ion pairs range between 0.4 k B T and 0.7 k B T.
The MF acting on every single ion is shown in Figure 4, jointly with the direct (F d ) and solvent (∆F(r)) contributions to the Li + -Br − , Li + -Li + , and Br − -Br − pairs. The contributions of the solvent (∆F(r)) are opposite to the direct ion-ion force in all the systems. This is constant with the trend of polar liquids to produce stable complexes in the case of like-ion pairs, and to separate ions with opposite signs [2][3][4][5]. Equilibrium between the direct and solvent contributions results in a repulsive force acting on the Li + -Li + pair at all distances. For the Li + -Br − pair, the resultant force is repulsive at extremely small distances and attractive at middle separation. The total force on the Br − -Br − pair has similar tendencies, but the attractive component is considerably less significant.   Figure 5 shows the PMF values for the various ion pairs. Oscillations in the potential corresponding to Li + -Br − occur in association with the structure of the solvent in the vicinity of the ions. W(r) shows a noticeable minimum at 2.2 Å and lower one at 8.0 Å. The contact ion pair (CIP) region corresponding to the first minimum, and for the solvent separated ion pair (SSIP) region the second one. W(r) shows a maximum at r = 4.4 Å. In the CIP region, the ion pairs are enclosed by acetonitrile particles, with the corresponding nitrogen atoms geared to the Li + and their methyl groups in the direction of Br − , as illustrated in Figure 6. The SSIP region values are attributed to attractive forces produced by solvent molecules between ion pairs. The stability of the CIP and SSIP regions is indicated by the association constants [25]. For the CIP and SSIP state we have the following expressions where N A is the Avogadro's number and r b is a cut-off parameter that we have chosen as the position of the first potential of mean force W(r) maximum, and r m is some value of r for which it is considered that the diffusive motion is a good approximation. Moreover, the equilibrium constant is giving by we calculate the association constants K CIP and K SSIP and the equilibrium constant, as shown in Table 6. We found that K CIP K SSIP . The value obtained from the equilibrium constant (K eq = 7.1 × 10 −12 ) indicates much less stable configurations in the SSIP region compared to the CIP region.  W(r) for Li + -Li + is repulsive and decreases monotonically to zero. Even though the acetonitrile molecules located near the ion pairs have their nitrogen atoms toward Li + , as shown in Figure 7, the electrostatic interaction N-Li + is less important than F d . The Br − -Br − pair conformation is steadied at a separation corresponding to the minimum (r = 7.0 Å) of W(r) by the acetonitrile molecules enveloping it. Their methyl atoms are oriented toward the Br − ions due to an attractive interaction. The structure is stabilized by the methyl atoms of the liquid molecule in the central area.
To carry out a more quantitative analysis of the distribution of the liquid molecules around the ionic particles, we determined the radial ion-solvent distribution functions (left axis) and their corresponding coordination numbers (right axis) for different interionic separations. Figures 8 and 9 show the radial distribution functions, Ion-Me, Ion-C and Ion-N, obtained in the situation of the Li + -Br − ion pair for the interionic separations corresponding to r CIP , r SIP and the maximum value of W(r). Table 7 summarizes the positions of the maximum and minimum radial distribution functions for the various cases analyzed. The corresponding coordination numbers are given in Table 8. A coordination number (cn N ) of 23 for a free bromide ion in acetonitrile is large. This value may be attributed to the deficiency in the interaction potentials used. Nevertheless, other studies indicate that the solvation shell of the Br-anion can have a maximum number of molecules of acetonitrile bigger than 10 [9]. It should also be pointed out that in aqueous solutions coordination numbers increase when the concentration decrease [26].

Li + -Br − Interconversion and Friction Kernels
The friction kernels ξ(t) depend on the ionic species and interionic distances [27]. We determined the transmission coefficient (k) of the barrier at 4.4 Å according to the Grote-Hynes theory of chemical reactions in a solution. In this model, the reaction coordinate develops by a global Langevin equation. The transmission coefficient (subscript GH or Kr stand for Grote-Hynes or Kramers) can be obtained as: where ω b is the barrier frequency, obtained by extrapolating an inverted parabola with an apex equivalent to the maximum PMF value. λ r (reactive frequency) is the result of the next expression: (12) Figure 9. Radial distribution functions and running coordination numbers for Br − at several interionic separations for the Li + -Br − ion pair.
It is follows that λ r is the characteristic frequency for the unstable motion at the transition state. For the coefficients involved in the previous theory regarding the Li + -Br − , Li + -Li + , and Br − -Br − systems, the initial values of the friction kernels (ξ(0)) are summarized in Table 9. The initial time values of the ξ(t) kernels are not only dependent on the ionic species but also show noticeable changes with the interionic separation. This interionic separation dependence of ξ(0) is strongly related to the magnitude of the changes in the solvent structure around the ion pairs as a function of these separations. The normalized friction kernels (ξ N (t) = ξ(t)/ξ(0)) related through all cases are displayed in Figures 10 and 11. There are features common to the functions of both systems, all of which show a rapid and very similar initial decay before 0.2 ps. A lengthy decay follows, characterized for each individual system. Table 7. Maximum (R ion−atom ) and minimum (r ion−atom ) radial distribution functions for the various cases.  Table 9. Initial values of the friction kernels for all ion pairs and interionic separations. Br − -Br − 7.0 0.14 The coefficients involved in the Li + -Br − system are summarized in the Table 10. λ −1 r , gives us a temporal measure of the transition state in the presence of the solvent. We determined that, for acetonitrile, λ −1 r = 2.17 ps. The correlation time of the friction of the solvent acting on the ion pair (τ c ) is defined by the value of the integral of the normalized kernel function:   Table 10. Li + -Br − association-dissociation process parameters and reaction constants.

Magnitude Value
In the case of acetonitrile, the correlation time is shorter than the reaction time.
A polarization caging regime is verified because a negative amount for the nonadiabatic barrier frequency (ω 2 N A = ω 2 b − ξ(0)) was found [28]. In this regime, the solvent immediately catches the solute in a well or "polarization cage" of solvent molecules. Hence, the motion of solvent molecules is extremely crucial. A limiting case of the polarization caging regime is when τ c λ −1 r and designated as adiabatic. In such circumstances, the predictions of the Grote-Hynes theory are reduced to Kramers doublets. Our results confirm that for the Li + -Br − ion pair, k GH k kr .
Lastly, we computed the full rate constants for the dissociation (k f ) and association (k b ) processes, given by [27,29]: whereby µ I is the reduced mass of the ion pair, k GH is the transmission coefficient, r = is the interionic separation at the transition state (i.e., the top barrier), β = 1/k B T and r m is an arbitrary interionic distance beyond the outer boundary of the SSIP region, where there is no further interactivity. The results are summarized in Table 10. We have obtained a k f rate constant that is roughly smaller than k b . This is due to the greater stability of the CIP complex for the present model system since a value of 7.1 × 10 −12 was obtained for K eq . Our overall constant rate k f + k b

Conclusions
The results reported in this paper on structural properties are in agreement with the experimental data; for instance, we carried out two unconstrained simulations of Li + and Br − , with 215 solvent molecules. The information obtained from these ion-solvent structure simulations is compared to the neutron diffraction data (Table 3). We also obtain a coordination number for Li + in acetonitrile, slightly higher than the experimental one (5.5 instead of 3- Table 4). This failure of simulations can be attributed to the deficiency in the interaction potentials used and the hydrogen atoms of the methyl group not being considered. We obtained PMF values for Li + -Br − , Li + -Li + , and Br − -Br − . For the Li + -Br − ion pair, the diluter allows sets two stable interionic separations, a CIP and a SSIP; where the CIP is more stable. The Li + -Li + system is repulsive. The Br − -Br − pair composition is equilibrated at a distance corresponding to the lowest value (r = 7.0 Å) of W(r) by the adjacent acetonitrile molecules. The friction kernels and interconversion process of the Li + -Br − ion pair were evaluated. The passage of the barrier happens in the polarization caging regime. The correlation time of the solvent τ c is much shorter than the reaction time λ −1 r and the Kramers regime of the Grote-Hynes theory was found for the system. In the Grote-Hynes theory, the reaction coordinate obeys a generalized Langevin equation [30,31].
where the friction coefficient is replaced by a memory function of friction kernel ξ(t), µ is a reduced mass of particles that conform the solute and β = 1 K B T . The first term of the above equation is the gradient of the potential of the mean force regarding the reaction coordinate. R(t) is a stochastic force with zero mean and is uncorrelated with the initial velocity and modeled by a Gaussian process. In this case, If the dynamics are studied in the vicinity of the transition state, the barrier can be approximated as an inverted parabola of frequency ω b . Grote and Hynes have deduced the following expression for the transmission coefficient, From the analysis of the model, it follows that λ r is the characteristic frequency for the unstable motion at the transition state. λ r is determined through the equation, ξ(λ r ) is the Laplace transform of the time dependent friction coefficient, The basic assumptions of this approach are: (a) The generalized Langevin equation is a valid description; (b) The friction kernel ξ '(t) does not depend on the reaction coordinate; (c) The mean force barrier is parabolic.
If we approach the friction coefficient dependent on time to a delta function; namely, (t) = ξδ(t), with the friction coefficient constant given by, Then, the Kramers theory can be recovered from this model and the transmission coefficient is given by, Another case is the non-adiabatic limit, where the solvent remains essentially static during the time scale of the motion at the barrier. The effect of the friction is almost constant and equal to its time zero value ξ(0). The transmission coefficient in the extreme case of constant friction is, where ω N A is called the non-adiabatic barrier frequency. In this case it can be interpreted as the solute moving under the influence of an instantaneous barrier of frequency ω N A instead of ω b , One particular situation is that in which ω 2 N A is negative; it is called the polarization caging regime. The reaction coordinate is trapped in a well rather than moving on an inverted parabola. Here, it is clear that the movement of the solvent is key to allowing this instantaneous configuration to relax.

Appendix A.2. Friction Kernels from Molecular Dynamics
The Grote-Hynes description of the reaction coordinate assumes previous knowledge of the potential of mean force and the friction kernel. There is a method for evaluating the friction kernel in the limit case called the frozen bond approximation [31,32].
The method is based on the following idea: when we assume the generalized Langevin equation, one of the hypotheses is that the whole dependence on the reaction coordinate is contained in the potential of the mean force. On this basis, it seems reasonable that ξ (t) should not depend on the interparticle potential.
When there is a harmonic potential of interaction for the coordinate of interest, it is possible to obtain the corresponding memory function as, Following the Mori formalism, it follows that with The renormalized frequency,ω where ξ (t) can be obtained through the relations .
Following this method, we can obtain the memory function parameterized as a function of x 0 and ω ξ (t) = ξ (t; x 0 , ω).
It has been demonstrated that the ω → ∞ limit corresponds to the aforementioned frozen bond approximation. This method is implemented in a molecular dynamics simulation program; the reaction coordinate is constrained to the x 0 during the simulation and the friction kernel is calculated as the autocorrelation function, ξ (r 0 , t) ≡ β µ δ∆F(t, r 0 ) δ∆F(0, r 0 ) , where the instantaneous values ∆F(t, r 0 ) are given by Equation (5) and, δ∆F(t, r 0 ) ≡ ∆F(t, r 0 ) − ∆F(r 0 ) .
The particles of interest lie in a central cell, which is surrounded by exact replicas. When one of the particles leaves the central shell through one face, its periodic image will enter it through the opposite face of the cube, so that we always have the same number of particles. The particles in the central cube are considered to interact with those of the neighboring images. The usual approximation is to fix a cut-off for the range of the interaction potentials. Each particle in the central box is considered to interact only with those that are at a distance lower than some value R c . In this work, R c was chosen equal to L/2, L being the length of the cubic central shell. This procedure is fine for the short-range interactions since they show a rapid decay with distance. However, the Coulomb interactions decay as r −1 . To overcome this problem, we used the Ewald summation technique. The basic idea is that the potential energy associated with the Coulomb interactions should be computed as a sum over all pairs for which at least one of the atoms lies in the central box [31,32].
(A22) → n = n x , n y , n z with n i integers. The prime indicates that the case i = j should not be considered for → n = 0, and that the given pair i-j must be excluded when → n = 0 if the charges belong to the same molecule.
If we take as additional parameters for A the positions of N solute particles (A (N, V, T, r 1 , . . . , r N )), the work done by the system in an ideal reversible process in which the N particles are infinitesimally moved is: We are interested in the free energy as a function of the interionic vector where → F r is the force on the interionic axis. For symmetry reasons, A will only depend on the modulus of r and multiplying (A29) with a unitary vectorr, then dA dr ≡ −F(r), the equality comes from the fact that < → F r > will be directed alongr. From this equation, we can calculate the free energy with a simple integration: where F(r) = F d (r) + ∆F(r), F d is just the direct force that can be calculated from the interparticle potential and ∆F(r) is the contribution from the solvent.