Solvation thermodynamics: two formulations and some misunderstandings

Abstract

Conventional (Fowler–Guggenheim) and Ben-Naim's formulations of solvation thermodynamics are analyzed in parallel, emphasizing their differences and stressing their interconnections. The pivotal equations relating the thermodynamic functions in both theories are derived. Connections with Pierotti–Abraham's cavity-interaction partition model are also contemplated in detail. Evidence is presented that misinterpretations of some of the derived equations lead to wrong estimates of solvation thermodynamic quantities that have been detected in the output of some of the most popular quantum chemistry software packages (Gaussian 09).

---

1. Introduction

Prediction of solvation free energies, ΔG^sol, is a subject of great interest in different fields like chemistry (thermodynamics¹ and kinetics²), biophysics and biochemistry,³ or drug design,⁴ among others. Particularly relevant and fruitful has been the use of solvation models in the study of systems of paramount importance in molecular biology by means of molecular simulations.⁵

In their classical text,⁶ first published in 1939, Fowler and Guggenheim presented the theory of liquids and solutions of non-electrolytes. In particular, these authors derived the properties of ideal dilute solutions by regarding them as special class of regular solutions.

On the other hand, Ben-Naim published in 1978 a groundbreaking contribution,⁷ where solvation thermodynamics was developed at the molecular level in a nonconventional manner by redefining the concept of solvation and by introducing an auxiliary quantity, the pseudochemical potential (PCP, μ*), which in the case of a two-component solution represents the change in the Gibbs energy caused by the addition, at a fixed position, of one solute molecule to the solution while keeping T, P and the number of solvent molecules unchanged. The conventional definition of the chemical potential (CP, μ), did not include any restriction to the location of the additional solute molecule. Ben-Naim has emphasized that this apparently subtle difference allowed for a much more genuine definition of the solvation Gibbs energy and its derived thermodynamic quantities.

Indeed, Ben-Naim has shown^1,7 that classical statistical thermodynamics leads to a very powerful relationship between both chemical potentials CP and PCP, namely, μ = μ* + kT log{ρΛ³} (see next sections), which is the pivotal expression of his novel formulation of solvation thermodynamics.

Ben-Naim has also stressed that while the conventional expressions for CP, like those found in the Fowler–Guggenheim text,⁶ only apply to very dilute solutions of solutes, the above mentioned equation is valid for any kind of molecule in any fluid mixture, and for any concentration. It is this general character that has contributed greatly to the present popularity of Ben-Naim's formulation which is being widely used in the recent literature.

While both formulations lead to correct results, they must be handled and interpreted within their own context because, as shown in this paper, both the meaning and the numerical values of conventional (Fowler–Guggenheim's) and Ben-Naim's thermodynamic magnitudes are different from each other.

Unfortunately, as usually happens, the coexistence of different formulations brings along the possibility of some degree of confusion related to the fact that different magnitudes can receive the same denomination, thus causing misunderstandings and misinterpretations. Indeed, as will be shown in this article, it is not uncommon to find wrong values of solvation thermodynamic quantities in the outputs of some of the most popular quantum chemical software packages. We do believe that a substantial part of the existing confusion can be avoided by carrying out an elemental parallel analysis of the conventional (Fowler–Guggenheim) and Ben-Naim's formulations, making clear their differences and emphasizing their interconnections.

The lack in the literature of any in-depth study where the equations connecting both formulations are derived in a systematic manner, prompted us to tackle the present research. As will be shown in this paper, our analysis helps to prevent potential confusions arising from the simultaneous references in the literature to both formulations. In particular, we will focus on the solvation thermodynamic functions of a solute in solution as computed by means of continuum solvation medium models. We will present evidence that some packages of software, which are employed as a reference software in quantum calculations by a large number of researchers, do not employ the correct expression for computing neither the partial enthalpies (H̄^sol_s) nor the partial entropies (S̄^sol_s) of a solute in solution. Our aim is to highlight, with the help of a chosen specific example, the fact that a parallel analysis of the conventional and Ben-Naim's theories of solvation can be very helpful in computing and interpreting, in a proper manner, the solvation thermodynamic functions.

The article is organized as follows: we start in Section 2 from Fowler–Guggenheim's formulas for the Helmholtz free energies of an ideal gas and an ideal dilute solution, and we derive mathematical expressions for the rest of thermodynamic functions. In Section 3, we start from Ben-Naim's equation for the chemical potential of a molecule in gas-phase and in solution, and we obtain mathematical expressions to compute the thermodynamic functions in this alternative formulation. Appropriate formulas to compute conventional and Ben-Naim's solvation thermodynamic functions are deduced in Section 4. Their relations to Ben-Naim's solvation ρ-process and x-process as well as to Pierotti–Abraham's solvation functions (widely employed in experimental studies for the prediction of thermodynamic solution parameters) are analyzed in some detail, thus allowing for a much deeper understanding about the meaning of solvation magnitudes in both formulations. The pivotal equations connecting the conventional and Ben-Naim's solvation thermodynamic functions are also obtained. Finally, in Section 5 we present a numerical application of the equations derived in previous sections, showing how the analysis carried out throughout this paper, helps to uncover errors detected in the output of some of the most popular quantum chemistry software.

2. Thermodynamic functions in Fowler–Guggenheim's formulation

Throughout this article, both molecular and molar units will be employed, as indicated. The notation used is the traditional one found in the literature. A slightly more elaborated notation is chosen for denoting the different standard states in order to stress their meaning.

2.1 Ideal gases

The Helmholtz free energy of an ideal gas can be written as⁶where N is the number of gas molecules of mass m occupying the volume V^g at temperature T, k is Boltzmann's constant, h Planck's constant, and j^g(T) denotes the partition function for all the internal degrees of freedom of the solute in gas-phase. E^g₀ represents the solute electronic energy. This term shifts the origin of energies to allow for a direct comparison with the quantum packages' outputs (see Tables 2 and 3). It is not included in Fowler–Guggenheim's original formulation. Superscripts “g” and “sol” indicate gas-phase and solution, respectively. The symbol “log” stands for natural logarithm.

From eqn (1), one gets

On the other hand, from eqn (3)

where Λ = h/(2πmkT)^1/2, ρ^g = N/V^g is the number density and kT(∂ log j^g(T)/∂T)_P = kT(∂ log j^g(T)/∂T) = E^g_int/T. In classical systems the momentum partition function Λ³ is independent of the environment, whether it is a gas or a liquid phase.¹

Eqn (3) can be written as

μ^g = E^g₀ − kT log kT + kT log{Λ³/j^g(T)} + kT log P (5)

It should be recalled at this point that the adoption of a specific standard state represents the choice of a given value for the molar volume (or number density) appearing in the equations for the calculations of chemical potential and partial molecular (molar) entropies (their translational part). We can define the standard partial Gibbs free energy and entropy of solute as

μ^o,g = E^g₀ + kT log{Λ³/[kTj^g(T)]} (6)

andwhere the standard state “o,g” corresponds to an hypothetical ideal gas at P = 1 bar (1 atm). Thus, using molar units and standard conditions (T = 298.15 K), we have RT/1 bar = V^g_m = 24.5 L mol⁻¹.^8,9 Let us symbolize such an standard state by: {g-ideal, P = 1 bar; V_m = 24.5 L mol⁻¹}. In some occasions (see Section 5), a concentration scale standard state is considered. Under such circumstances eqn (6) and (7) become

μ^⊕,g = E^g₀ + kT log{ρ^⊕Λ³/j^g(T)} (8)

where ρ^⊕,g is the number density defining the standard state “⊕,g”. For an ideal gas at P = 1 bar and 298.15 K we have ρ^⊕,g = 1 bar/RT = 1/24.5 mol L⁻¹ = 0.041 M. Let us symbolize such an standard state by {g-ideal, ρ^g = 0.041 M; V_m = 24.5 L mol⁻¹} or in general {g-ideal, ρ^g; V_m}, when the value of ρ^g is different from 0.041 M (P ≠ 1 bar).

It should be noted that standard states “o,g” and “⊕,g” are fully equivalent when ρ^⊕,g = ρ^g = 0.041 M [i.e. μ^⊕,g = μ^o,g and S̄^⊕,g = S̄^o,g; see eqn (6)–(9)].

It is straightforward to show that eqn (3) and (2) lead, respectively, to

andwhere, for an ideal gas, the thermal expansion coefficient α^g_p = (∂V^g/∂T)_P/V^g = 1/T.

2.2 Ideal dilute solutions

According to Fowler and Guggenheim,⁶ the Helmholtz free energy of an ideal dilute solution of a solute, “s”, in a solvent, “l”, (solute–solvent assembly) can be given by

F^sol = F_l(T, V^sol − N_sV̄^sol_s, N_l) + N_s[E^sol_0,s − χ_ls + kT log{Λ_s³/V^sol} − kT + kT log N_s − kT log j^sol_s(T)] (12)

where N_s solute molecules, having a constant potential energy −χ_ls (representing the solute–solvent interaction), move freely in a region of volume V^sol. The N_l solvent molecules are occupying the V^sol_l = V^sol − N_sV̄^sol_s volume available and make a contribution F_l(T, V^sol − N_sV̄^sol_s, N_l) to the total Helmholtz free energy. j^sol_s(T) represents the contribution of the internal motions of the solute molecules and E^sol_0,s the solute electronic energy, the latter not included in Fowler–Guggenheim's original formulation.

From eqn (12), one obtains for the solute CP

It is easy to show (see Appendix) that (∂F_l/∂N_s)_T,V,Nl ≈ PV̄^sol_s, and then

μ^sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + kT log{Λ_s³/V^sol} + kT log N_s − kT log j^sol_s(T) (14)

According to Pierotti,⁸ −χ_ls is the molar potential energy of the solute in the solution relative to infinite separation, P the hydrostatic pressure, and −χ_ls + PV̄^sol_s represents the reversible work required to introduce one solute molecule into a solution of concentration N_s/V^sol (see in the next section its relation to Ben-Naim's W(s|l) coupling work of solute to the system).^1,7 As remarked by Fowler and Guggenheim,⁶ PV̄^sol_s is usually negligible at ordinary pressures (see below) and it is not considered when calculating the derivative with respect to T. On the other hand, Fowler and Guggenheim also mention⁶ that the term −χ_ls may depend on T, but in this work (see Section 5), −χ_ls [or W(s|l); see eqn (29) and (30) in the next section] has been estimated through the Gaussian 09 (G09) package of programs,¹⁰ by adopting the generalized Born approximation and the formalism of atomic surface tensions and the solvent-accessible surface areas as a model for the solvent (including cavity, dispersion and solvent reorganization energy contributions), where no dependence on T is considered (the T dependence inherent to the dielectric constant⁶ is not taken into account). From the operational viewpoint, the above implies that the term [∂(−χ_ls + PV̄^sol_s)/∂T]_P,Nl,Ns, which contributes to H̄^sol_s through and to S̄^sol_s through , becomes zero. It is straightforward to confirm that the terms [∂(−χ_ls + PV̄^sol_s)/∂T]_P,Nl,Ns appearing in the expressions for H̄^sol_s and S̄^sol_s when the T dependence of −χ_ls and PV̄^sol_s terms is taken into consideration, cancel each other out when the corresponding μ^sol_s magnitude is computed from H̄^sol_s − TS̄^sol_s. In other words, the expression found for μ^sol_s [i.e. eqn (14)] does not depend on whether or not the T dependence of −χ_ls + PV̄^sol_s has been assumed. As will be remarked in Section 4, the assumption that −χ_ls + PV̄^sol_s does not exhibit T dependence is not, in general, a very accurate approach.

Consequently, from eqn (14) we find

where α^sol_p = (∂V^sol/∂T)_P,Ns,Nl/V^sol ≈ (∂V^l/∂T)_P/V^l = α^l_p is the thermal expansion coefficient of the solvent and kT(∂ log j^sol_s(T)/∂T)_P,Ns,Nl = kT(∂ log j^sol_s(T)/∂T) = E^sol_int,s/T.

On the other hand, eqn (14) can be written as

μ^sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + kT log{Λ_s³/j^sol_s(T)} + kT log N_s/N_T − kT log(V^sol/N_T) (16)

where N_T = N_s + N_l is the total number of molecules in solution (solute plus solvent).

Now, for a ideal dilute solution where N_l/N_T → 1 and V^sol → V_l (V^sol/N_T → V_l/N_l = V_l,Nl) one gets

μ^sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + kT log{Λ_s³/[j^sol_s(T)V_l,Nl]} + kT log x_s (17)

where V_l,Nl represents the pure solvent molecular volume (inverse of the pure solvent number density). We can now define the standard partial Gibbs free energy of the solute (conventional standard CP) as

μ^o,sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + kT log{Λ_s³/[j^sol_s(T)V_l,Nl]} (18)

Similar arguments applied to eqn (13) lead to the following equation for the standard partial entropy of the solute

The standard state “o, sol” (mole fraction scale) corresponds to an hypothetical ideal solution (fulfilling Henry's law) at T and P of the solution, where the solute mole fraction is the unity (x_s = 1). Thus, using molar units and standard conditions for V_l,Nl, we have V_l,m = 0.018 L mol⁻¹.^8,9 Let us symbolize such an standard state by: {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹}. Fig. 1 collects, in a p_solute − x_solute graphic, the different standard states adopted for the solute throughout the present work.

Fig. 1 The red point represents an ideally dilute solution fulfilling Henry's law. The orange: {sol-Henry, ρ^sol_s = 0.041 M; V_m = 24.5 L mol⁻¹}, green: {sol-Henry, ρ^sol_s = 1 M; V_m = 1 L mol⁻¹} and blue: {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹} points represent fictitious states of the solute in which each solute molecule experiences the same intermolecular forces it experiences in the ideally dilute solution, where it is surrounded by solvent molecules. They are the hypothetical states arising from an extrapolation of the properties of solute in the very dilute solution to the limit of 0.041 M, 1 M and x_solute = 1 concentrations, corresponding to the different standard states mentioned throughout the present work. The orange point corresponds to the standard state adopted by the Gaussian 09 software in solvation calculations.

In order to emphasize that eqn (14) and (15) correspond to a concentration scale standard state, they can be rewritten as

μ^⊕,sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + kT log{ρ^⊕,sol_sΛ_s³/j^sol_s(T)} (20)

where ρ^⊕,sol_s is the number density defining the standard state “⊕, sol”. It corresponds to an hypothetical ideal solution (fulfilling Henry's law) at T and P of the solution, where the solute concentration is ρ^⊕,sol_s. Let us symbolize such a standard state: {sol-Henry, ρ^sol_s; V_m}, for a solution of concentration ρ^⊕,sol_s = ρ^sol_s as mentioned in Section 2.1, we need to introduce this additional standard state here to analyze the G09 output in Section 5.

On the other hand, from eqn (14)

It is important to stress at this point that while PV̄^sol_s arises from the first term contribution (F_l) to the solution Helmholtz free energy (F^sol) in eqn (12), kT²α^l_p comes from the term −N_skT log V^sol. Both contributions are in most cases negligible (for instance, using molar units and standard conditions, in the case of solute water in water solution: RT²α^l,water_p = 4.5 × 10⁻² kcal mol⁻¹ and PV̄^sol_water = 4.4 × 10⁻⁴ kcal mol⁻¹), but according to Fowler and Guggenheim,⁶ they should not be omitted until the final stage to avoid apparent inconsistencies.

On the other hand, from eqn (22) we obtain

3. Thermodynamic functions in Ben-Naim's formulation

Ben-Naim starts from the relation between the CP of a molecule s, μ^φ_s, in gas-phase (φ = g) or in solution (φ = sol), and the corresponding PCP, μ^*,φ_s^1,11

μ^φ_s = μ^*,φ_s + kT log{ρ^φ_sΛ_s³} (24)

Both μ_s and μ^*_s represent the change in Gibbs energy caused by the addition of one molecule s to the system keeping T, P, and the number of other type of molecules present in the system (if any) unchanged. While for μ^*_s the center of mass of the added molecule is placed at a fixed position, in the case of μ_s the new molecule is released to wander in the entire volume.

According to Ben-Naim¹

μ*^,g = E^g₀ + W(s|g) − kT log j^g(T) = E^g₀ − kT log j^g(T) (25)

μ^*,sol_s = E^sol_0,s + W(s|l) − kT log j^sol_s(T) (26)

where the contribution of the electronic energy of solute, E^φ_0,s, has been included.

W(s|x), Ben-Naim's coupling work, is the average Gibbs energy of interaction of s with its entire surroundings (other gas molecules, x = g, or solvent molecules, x = l). Of course, for an ideal gas, the coupling work, W(s|g), is zero. Note that in the gas phase, where only one type of molecules is present, the subscript “s” becomes unnecessary and has been eliminated, as we also did in the case of Fowler–Guggenheim's formulation (see Section 2).

Eqn (24)–(26) lead to

μ^g = E^g₀ + kT log{ρ^gΛ³/j^g(T)} (27)

which, of course, matches eqn (5) (ρ^g = N_s/V^g = P/kT), and

μ^sol_s = E^sol_0,s + W(s|l) + kT log{ρ^sol_sΛ_s³/j^sol_s(T)} (28)

Comparison of eqn (28) with (14) in Fowler–Guggenheim's formulation (ρ^sol_s = N_s/V^sol) allows one to conclude that

W(s|l) = −χ_ls + PV̄^sol_s (29)

Thus, Ben-Naim's coupling work of the solute s to the system, corresponds to Fowler–Guggenheim's reversible work required to introduce one solute molecule into the solution (at one fixed position).^12,13 Indeed, Ben-Naim describes in detail how the solvation of any solute in any solvent can always be decomposed into two parts: creation of a suitable cavity (PV̄^sol_s) and then turning on the other parts of the solute–solvent interaction (−χ_ls).^12,13 As PV̄^sol_s is usually negligible (see Section 2.2), we have

W(s|l) ≈ −χ_ls (30)

Ben-Naim's definition of μ^*_s can be generalized to any partial molecular thermodynamic function, X̄^*_s: it will represent the change in the particular function considered due to the addition of one molecule of solute at a fixed position in the system. Let us obtain explicit expressions for Ben-Naim's partial molecular entropies, enthalpies and internal energies.

Bearing in mind that S̄_s = −(∂μ_s/∂T)_P, one gets from eqn (25)–(28)

where, as already mentioned above, for an ideal gas, α^g_p = 1/T, and for an ideal dilute solution, α^sol_p ≈ α^l_p (thermal expansion coefficient of the solvent). Combination of eqn (31) and (4), and eqn (32) and (15) leads to the following explicit expression for Ben-Naim's partial molecular entropies of the solute in gas-phase and in solution

S̄^*,φ_s = k log j^φ_s(T) + E^φ_int,s/T (φ = g, sol) (33)

The solute partial molecular enthalpies H̄^g and H̄^sol_s are obtained from H̄^φ_s = μ^φ_s + TS̄^φ_s, using eqn (24), (31) and (32)

Although rather obvious, it is important to stress at this point that kT²α^g_p ≠ kT²α^l_p. Indeed, in the case of water, for instance, using molar units and standard conditions, RT²α^g,water_p = 0.59 kcal mol⁻¹ and RT²α^l,water_p = 4.5 × 10⁻² kcal mol⁻¹. While the latter can be neglected without appreciable lost of accuracy for most of applications, the former becomes of the order of standard state corrections and should be taken into consideration.

A second important remark is that for an ideal gas, the term kT that must be added to Ū^g_s to get H̄^g_s (see eqn (10)–(11)), arises from the contributing term kT²α^g_p. Since kT = PV̄^g, one can write H̄^g = Ū^g + PV̄^g. In other words, for an ideal gas kT²α^g_p = PV̄^g_s. In sharp contrast, for the case of a solute in solution, the terms kT²α^l_p and PV̄^sol_s are different contributions to H̄^sol_s [see eqn (22)], as stressed in the previous section.

The above points, although trivial, should be taken into account in order to avoid misunderstandings leading to some serious mistakes detected in the literature. In particular, we will show in Section 5 that some rather popular packages of software (for instance, the G09 package of programs),¹⁰ with a huge number of users, do contribute to create confusion in this regard.

It is easy to show that eqn (25) and (34) with S̄*^,g = −(∂μ*^,g/∂T)_P,Nl,Ns, within Ben-Naim's formulation, lead to eqn (10) in Fowler–Guggenheim's formulation. Similarly, eqn (26), (29) and (35) with S̄^*,sol_s = −(∂μ^*,sol_s/∂T)_P,Nl,Ns, result in eqn (22). Consistency between both formulations is thus, once again, confirmed.

Combination of eqn (34) and (10) on one hand, and eqn (22), (29) and (35) on the other hand, leads to explicit expressions for Ben-Naim's partial molecular enthalpies of the solute in gas-phase and in solution,

H̄*^,g = E^g₀ + E^g_int (36)

H̄^*,sol_s = E^sol_0,s − χ_ls + PV̄^sol_s + E^sol_int,s (37)

In order to get the expressions for Ben-Naim's partial molecular internal energies, we start from Ben-Naim's expression for the partial molecular volume¹

V̄^φ_s = V̄^*,φ_s + kTβ^φ_T (φ = g, sol) (38)

where β^φ_T = (∂V^φ/∂P)_T/V^φ is the isothermal compressibility coefficient, that for an ideal gas reduces to 1/P and for a liquid phase becomes so small that can be safely neglected.¹

Bearing in mind that H̄*^,φ = Ū*^,φ + PV̄*^,φ (φ = g, sol), eqn (36)–(38) lead to

Ū*^,g = E^g₀ + E^g_int (39)

Ū^*,sol_s = E^sol_0,s − χ_ls + E^sol_int,s (40)

4. Solvation thermodynamic functions

In the previous section we have derived explicit expressions for Ben-Naim's solute partial molecular properties X̄^*,φ_s. The relevance of such magnitudes lies in the fact that the solvation thermodynamic functions can be defined in a very convenient way in terms of Ben-Naim's X̄^*,φ_s (see below).

In Fowler–Guggenheim's formulation, the solvation Gibbs free energy is obtained from eqn (6) and (18)

ΔG^o,sol = μ^o,sol_s − μ^o,g_s = E^sol_0,s − E^g₀ − χ_ls + PV̄^sol_s + kT log(V^g_m/V_l,m) (41)

where we have used the fact that kT/(1 bar V_l,Nl) = RT/(1 bar V_l,m) = V^g_m/V_l,m. As stressed in Section 2, according to the standard states employed in this formulation (see Section 2),⁸ V^g_m = RT/1 bar is the molar volume of an ideal gas at P = 1 bar and T = 298.15 K (V^g_m = 24.5 L mol⁻¹), and V_l,m is the pure solvent molar volume (V_water,m = 0.018 L mol⁻¹ for water solutions under standard conditions). As usual, we assume throughout this work that the thermal contributions from internal motions of the system (i.e. rotation and vibration) are very similar in gas-phase and in solution, thus cancelling each other out [j^g(T) = j^sol_s(T) and E^g_int = E^sol_int,s].

Eqn (41) is fully consistent with the works by Pierotti⁸ and Abraham and Nasehzadeh¹⁴ where solvation molar Gibbs free energies are computed as

ΔG^o,sol = Ḡ_CAV + Ḡ_INT + RT log(V^g_m/V_l,m) (42)

with Ḡ_CAV and Ḡ_INT being the partial molar Gibbs energies associated with the creation of a cavity (PV̄^sol_s) and with the solute–solvent interactions when the solute is in the cavity (−χ_ls), respectively. Solute electronic contribution E^sol_0,s − E^g₀ must be added to eqn (42) to make it compatible with eqn (41).

Eqn (41) represents, within the conventional formulation of thermodynamics of solution, the energy involved in the transference of a solute between gas-phase in the standard state: (hypothetical gas ideal with P = 1 bar) whose chemical potential is μ^o,g_s [see eqn (6)], and solution in the standard state: (hypothetical ideal dilute solution fulfilling Henry law with x_solute → 1) whose chemical potential is μ^o,sol_s [see eqn (18)].¹⁴ It is the so-called “solvation x-process” in Ben-Naim's formulation. That is to say, ΔG(x-process) = ΔG^o,sol. In other words, Ben-Naim's x-process corresponds to Fowler–Guggenheim's solvation process.

Ben-Naim also defines what he denotes as the “solvation ρ-process”, in which one molecule of solute is transferred from an ideal gas phase into an ideal solution (Henry's law) at fixed temperature and pressure and such that ρ^sol_s = ρ^g. Application of eqn (5) and (14) leads to

Combination of eqn (25), (26), (29) and (43) leads to

ΔG*^,sol = μ^*,sol_s − μ^*,g_s = ΔG(ρ-process) (44)

where ΔG*^,sol represents Ben-Naim's solvation Gibbs free energy.

It is important to stress, as Ben-Naim does,¹ that while ΔG*^,sol and ΔG(ρ-process) are equal in magnitude, they correspond to two different processes. The process involved in the definition of ΔG*^,sol represents the change in Gibbs free energy when moving one solute molecule of solute from a fixed position in an ideal gas phase into a fixed position in an ideal dilute solution at constant temperature and pressure, which obviously differs from the ρ-process as defined above.

The claimed superiority of Ben-Naim's solvation magnitudes ΔX*^,sol when compared with their related Fowler–Guggenheim's ΔX^o,sol values have been summarized by that author as follows:¹⁵ (a) while the conventional standard state free energy is not a bona fide measure of the average free energy of interaction of a solute with its surroundings, ΔG*^,sol is a direct measure of it. Indeed, it can be shown that ΔG*^,sol = −kT〈exp[−B_s/kT]〉, where B_s is the total interaction energy of a solute s at a fixed position, (b) while ΔG*^,sol = μ^*,sol_s − μ*^,g applies to solutions of any concentration, the conventional standard free energies ΔG^o,sol = μ^o,sol_s − μ^o,g are only valid for very dilute solutions, (c) while the magnitudes ΔX*^,sol pertain to the same process of solvation, when using the conventional standard magnitudes of solution, one usually applies different processes to different thermodynamic magnitudes, and (d) while the specification of a standard state is mandatory when using conventional standard free energies, no standard state definition is required in Ben-Naim's formulation.

Combination of eqn (41), (43) and (44) leads finally to the pivotal equation which establishes the relationship between solvation Gibbs free energies in Fowler–Guggenheim's and Ben-Naim's formulations¹¹

ΔG(ρ-process) = ΔG*^,sol = ΔG^o,sol − kT log(V^g_m/V_l,m) (45)

On the other hand, eqn (43) and (44) lead to

ΔG*^,sol = E^sol_0,s − E^g₀ − χ_ls + PV̄^sol_s (46)

which is the explicit expression for Ben-Naim's Gibbs free energy of solution.

The solvation molecular enthalpy is obtained from eqn (10) and (22)

ΔH^o,sol = H̄^sol_s − H̄^g = E^sol_0,s − E^g₀ − χ_ls + PV̄^sol_s − kT + kT²α^l_p (47)

Eqn (47) is fully consistent with the results of Pierotti⁸ and Abraham and Nasehzadeh¹⁴ who compute the solvation molar enthalpies according to

ΔH^o,sol = H̄_CAV + H̄_INT + α^l_pRT² − RT (48)

with H̄_CAV and H̄_INT being the partial molar enthalpies associated with the creation of a cavity and to the solute–solvent interactions when the solute is in the cavity, respectively. Solute electronic contribution E^sol_0,s − E^g₀ must be added to eqn (48) to make it compatible with eqn (47).

It should be recalled that Pierotti⁸ and Abraham and Nasehzadeh¹⁴ assume that Ḡ_INT = H̄_INT but Ḡ_CAV = H̄_CAV − TS̄_CAV [with S̄_CAV = −(∂PV̄^sol_s/∂T)_P,Nl,Ns], thus implying that while −χ_ls (interaction term) is independent of T, PV̄^sol_s (cavity term) is allowed to depend on T. As mentioned in Section 2.2, in the present work, neither −χ_ls nor PV̄^sol_s are considered to be T dependent and therefore, not only Ḡ_INT = H̄_INT = −χ_ls but also Ḡ_CAV = H̄_CAV = PV̄^sol_s [see eqn (42)]. In this regard, it is easy to show from eqn (7) and (19) that

TΔS^o,sol = TS̄^o,sol_s − TS̄^o,g_s = −kT + kT²α^l_p− kT log(V^g_m/V_l,m) (49)

and also from eqn (7), (19) and (33) we find

TΔS^o,sol = TΔS̄*^,sol − kT + kT²α^l_p − kT log(V^g_m/V_l,m) (50)

which represents the pivotal equation connecting the solvation entropies in Fowler–Guggenheim's and Ben-Naim's formulations.

From eqn (41), (47) and (49) it is seen that, as expected, ΔG^o,sol = ΔH^o,sol − TΔS^o,sol. Of course, eqn (47) could also be derived from this latter equation, using ΔS^o,sol = −(∂ΔG^o,sol/∂T)_P,Nl,Ns together with eqn (41), and assuming that −χ_ls + PV̄^sol_s is independent of T (E^sol_0,s − E^g₀, involving the electronic energies of the molecule in solution and gas-phase, does not depend on T).

It is clear either from eqn (33) or from (49) and (50) that ΔS*^,sol = 0. Of course, this is again a direct consequence of the simplification adopted in this work that W(s|l) (or equivalently −χ_ls + PV̄^g_s) is independent of T. As stressed in Section 2.2, such a restriction has been imposed for consistency with G09,¹⁰ for we will analyze the thermodynamic solvation G09 output in the last section of this article. However, we must insist that it is a rather crude approximation. Indeed, in the case of water solute in water solution at 298.15 K, ΔG*^,sol = −6.324 kcal mol⁻¹, ΔH*^,sol = −9.974 kcal mol⁻¹, and ΔS*^,sol = −12.24 cal K⁻¹ mol⁻¹.¹¹

It is easy to show¹ that, unlike solvation Gibbs free energies, solvation enthalpies for Ben-Naim's ρ- and x-processes are identical to each other. Bearing in mind that Fowler–Guggenheim's solvation enthalpy, ΔH^o,sol, does correspond to the change in enthalpy of Ben-Naim's x-process (i.e. the change in enthalpy involved in the x-process is H̄^sol_s − H̄^g), we can finally write

ΔH^o,sol = ΔH(ρ-process) = ΔH(x-process) (51)

From eqn (34) and (35) we obtain

ΔH^o,sol = H̄^sol_s − H̄^g = H̄^*,sol_s − H̄*^,g − kT + kT²α^l_p = ΔH*^,sol− kT + kT²α^l_p (52)

which represents the pivotal equation connecting the solvation enthalpies in Fowler–Guggenheim's and Ben-Naim's formulations.¹¹

On the other hand, comparison of eqn (47) and (52) leads to

ΔH*^,sol = E^sol_0,s − E^g₀ − χ_ls + PV̄^sol_s (53)

which is the explicit expression for Ben-Naim's enthalpy of solution.

It is important to emphasize that unlike solvation Gibbs free energy [see eqn (45)], Ben-Naim's solvation enthalpy, ΔH*^,sol, differs from the change in enthalpy for the ρ-process, ΔH(ρ-process). Eqn (47) and (53) imply that, likewise for Gibbs energies, when computing solvation enthalpies, two different estimates can be done, namely, Fowler–Guggenheim's ΔH^o,sol, or Ben-Naim's ΔH*^,sol. Thus, while authors like Wilhelm et al.^16,17 or Abraham et al.¹⁴ report ΔH^o,sol (as well as ΔG^o,sol) values, Ben-Naim^1,10 reports ΔH*^,sol (as well as ΔG*^,sol) values. In this context, it is interesting to note that in Cabani et al.'s compilation,⁹ one finds Fowler–Guggenheim's ΔH^o,sol enthalpies and Ben-Naim's ΔG*^,sol Gibbs free energies.

Regarding solvation internal energies, eqn (11) and (23) lead to

ΔU^o,sol = Ū^sol_s − Ū^g = E^sol_0,s − E^g₀ − χ_ls + kT²α^l_p (54)

and from eqn (39) and (40) we get

ΔU*^,sol = Ū^*,sol_s − Ū*^,g = E^sol_0,s − E^g₀ − χ_ls (55)

Combination of eqn (54) and (55) gives

ΔU*^,sol = ΔU^o,sol − kT²α^l_p (56)

which represents the pivotal relationship between solvation internal energies in Fowler–Guggenheim's and Ben-Naim's formulations.

Before ending this brief parallel analysis on Fowler–Guggenheim's and Ben-Naim's formulations of solvation thermodynamics, we would like to mention that Jorgensen's group¹⁸ computed solvation enthalpies by using the expression (in molar units)

ΔH^o = ΔE^o + PΔV^o − RT (57)

where¹⁸ ΔE^o is the energy change on transferring the solute from the gas-phase to solution [ΔU^o,sol in our notation: see eqn (54)], ΔV^o is the partial molar volume of the solute (V̄^sol_s in our notation). Eqn (47) and (54) lead to

ΔH^o,sol = ΔU^o,sol + PV̄^sol_s − RT= ΔU^o,sol + P(V̄^sol_s − V̄^g) = ΔU^o,sol + PΔV^o,sol (58)

with the solvation volume ΔV^o,sol = V̄^sol_s − V̄^g. Eqn (58) is equivalent to eqn (57) (note that PΔV^o in Jorgensen's notation is PV̄^sol_s), thus showing that the equations derived in this work are consistent with Jorgensen's formula.

In the next section, we will show that a certain confusion has been detected in the literature regarding the correct application of the equations derived in this article. In particular, we will present evidence that the G09 package of programs,¹⁰ which is employed as a reference software in quantum calculations by a large number of researchers, does not employ the correct expression for computing neither the partial enthalpies (H̄^sol_s) nor the partial entropies (S̄^sol_s) of a solute in solution.

5. Practical application

Table 1 gathers the most important equations derived in the previous sections. Tables 2 and 3 collect the computed values of the thermodynamic magnitudes considered in the present work for the particular case of water solute in water solution, which has been chosen as a representative system where we focus our analysis.

Table 1 Thermodynamic functions for a solute in gas-phase, (X̄^g, X̄*^,g), in solution, (X̄^sol_s, X̄^*,sol_s), and solvation thermodynamic functions (ΔX^o,sol, ΔX*^,sol). Contribution −χ_ls + PV̄^sol_s has been assumed to be independent of T (see the text for more details). Molar units are employed

Solute in gas-phase
	Ū*^,g = E^g0 + E^gint
	H̄*^,g = E^g0 + E^gint
	TS̄*^,g = RT log j^g(T) + E^gint
Ḡ^o,g = μ^o,g = E^g0 + RT log Λ^3/[RTj^g(T)]	Ḡ*^,g = μ*^,g = E^g0 − RT log j^g(T)
Solute in ideal dilute solution
	Ū^*,sols = E^sol0,s − χls + E^solint,s
	H̄^*,sols = E^sol0,s − χls + PV̄^sols + E^solint,s
	TS̄^*,sols = RT log j^sols(T) + E^solint,s
Ḡ^o,sols = μ^o,sols = E^sol0,s − χls + PV̄^sols + RT log Λs^3/[j^sols(T)Vl,m]	Ḡ^*,sols = μ^*,sols = E^sol0,s − χls + PV̄^sols − RT log j^sols(T)
Solvation functions
ΔU^o,sol = E^sol0,s − E^g0 − χls + RT^2α^lp	ΔU*^,sol = E^sol0,s − E^g0 − χls
ΔH^o,sol = E^sol0,s − E^g0 − χls + PV̄^sols − RT + RT^2α^lp	ΔH*^,sol = E^sol0,s − E^g0 − χls + PV̄^sols
TΔS^o,sol = −RT + RT^2α^lp − RT log[RT/Vl,m]	TΔS*^,sol = 0
ΔG^o,sol = E^sol0,s − E^g0 − χls + PV̄^sols + RT log[RT/Vl,m]	ΔG*^,sol = E^sol0,s − E^g0 − χls + PV̄^sols

---

Table 2 Solute electronic energies (Schrödinger equation) (E_o), energy contributions of the internal molecular degrees of freedom (E_int), molar potential energy of the solute in the solution (−χ_ls), partial molar internal energies (Ū), partial molar enthalpies (H̄) and partial molar Gibbs free energies (Ḡ = μ) for water in gas phase (g) and water in solution (sol). All numbers in hartree (T = 298.15 K)^a

E^go	E^gint	Ū^g	H̄^g	TS̄^o,g	Ḡ^o,g = μ^o,g	Ū*^,g	H̄*^,g	TS̄*^,g	Ḡ*^,g = μ*^,g
a The Gaussian 09 (G09) calculations leading to the output energies collected in this table were obtained by running the G09 route cards: #b3lyp/6-31G(d,p) opt(calchffc) freq and #b3lyp/6-31G(d,p) opt(calchffc) scrf=(solvent=water) freq for the gas-phase and solution calculations, respectively. When single values appear in the table, they correspond to the magnitude as computed with the formulae in Table 1. Values in parentheses, are the ones in the G09 output that differ from the values computed with formulas in Table 1.b Assuming PV̄^sols ≈ 0.c Computed from , as proposed in this work (see eqn (22) with PV̄^sols ≈ 0 and RT^2α^lp ≈ 0).d Computed from , with Vl,m = 0.018 L mol^−1, as proposed in this work [see eqn (19) with RT^2α^lp ≈ 0].e Computed from μ^o,sols = E^sol0,s − χls + RT log{Λs^3/[j^sols(T)Vl,m]}, with Vl,m = 0.018 L mol^−1, as proposed in this work [see eqn (18) with PV̄^sols ≈ 0].f Computed from , as proposed in Gaussian 09 package of programs [see eqn (62)].g Computed from , with V^gm = 24.5 L mol^−1, as proposed in Gaussian 09 package of programs [see eqn (60)].h Computed from μ^⊕,sols(G09) = H̄^sols(G09) − TS̄^⊕,sols(G09) = E^sol0,s − χls + RT log{Λs^3/[j^sols(T)V^gm]}, with V^gm = 24.5 L mol^−1, as proposed in Gaussian 09 package of programs.
−76.419737	0.022784	−76.395537	−76.394593	0.021436	−76.416029	−76.396953	−76.396953	0.004997	−76.401950
E^solo,s − χls^b	E^solint,s	Ū^sols	H̄^sols	TS̄^o,sols	Ḡ^o,sols = μ^o,sols	Ū^*,sols	H̄^*,sols	TS̄^*,sols	Ḡ^*,sols = μ^*,sols
−76.426730	0.022736	−76.402578	−76.402578^c	0.013689^d	−76.416267^e	−76.403994	−76.403994	0.005002	−76.408996
			(−76.401634)^f	(0.021441)^g	(−76.423076)^h

---

Table 3 Solvation thermodynamic functions for water solute in water solution, as computed by means of the conventional (Fowler–Guggenheim) and Ben-Naim's formulations together with the corresponding values as estimated from the Gaussian 09 outputs for gas-phase and solution runs. The approximations PV̄^sol_s ≈ 0 and RT²α^l_p ≈ 0 have been adopted in all cases. All numbers in kcal mol⁻¹

	ΔU^sol	ΔH^sol	TΔS^sol	ΔG^sol	ΔG^sol(exp.)^d	Gas-phase standard state^e	Solute in solution standard state^f
a As computed from ΔX^o,sol = X̄^sols − X̄^g (X = U, H, S^o, G^o), using the values collected in Table 2.b As computed from ΔX*^,sol = X̄^*,sols − X̄*^,g (X = U, H, S, G), using the values collected in Table 2.c As computed from ΔX^o,sol = X̄^sols(G09) − X̄^g(G09) (X = U, H, S^⊕, G^⊕) with X̄^sols(G09) and X̄^g(G09) taken from the Gaussian 09 outputs for the gas-phase and solution calculations, respectively.d From ref. 20.e The adopted gas-phase standard state corresponds to an ideal gas at standard conditions (P = 1 bar, T = 298.15 K), occupying the molar volume Vm.f The adopted standard state for the solute in solution corresponds to an hypothetic ideal solution (fulfilling Henry's law) at P and T of the solution (we assume in this work standard conditions), where: (i) the solute mole fraction tends to unity (xs → 1) or (ii) RT/1 bar = 24.5 L mol^−1 substitutes Vl,m in eqn (19) (i.e. concentration scale is employed to define the standard state, with ρ^⊕,sols = ρ^sols = 0.041 M; see eqn (21) and the text for details).g Ben-Naim's solvation magnitudes, ΔX*^,sol, do not involve the definition of any standard state. The values in the table were computed from data in columns 7–10 in Table 2.h G09 uses the same value of V^φm [see eqn (59)] for the calculations in both phases (24.5 L mol^−1). As discussed in the text, it means using concentration scale standard states with identical concentration for gas-phase and solution (ρ^⊕ = ρ^g = ρ^sols = 0.041 M, in molar units).
Fowler–Guggenheim^a	−4.4	−5.0	−4.9	−0.1	−2.05	{g-ideal; Vm = 24.5 L mol^−1}	{sol-Henry; Vm = 0.018 L mol^−1}
Ben-Naim^b	−4.4	−4.4	0.0	−4.4	−6.32	—^g	—^g
Gaussian 09^c	−4.4	−4.4	0.0	−4.4		{g-ideal; ρ^g = 0.041 M}^h	{sol-Henry; ρ^sols = 0.041 M}^h

---

Calculations were carried out by means of the G09 package of programs.¹⁰ The equations implemented in G09 for computing thermochemical data have been described in detail in a technical report available at the official Gaussian Website.¹⁹ Tables 2 and 3 compile the values appearing in the G09 output file together with the results obtained from the application of the equations derived in this work (see Table 1). Since the conclusions of the present work do not depend at all on the level of theory employed in the calculations, we have chosen a rather standard low-level (see Table 2 for details) to avoid diverting attention of the reader from the nuclear point analyzed in this work, namely, the correct use of solvation thermodynamic formulas in practical cases.

The values in Table 2 corresponding to the gas-phase magnitudes (row 2) in columns 1–6 have been taken from the G09 output and match those computed from the equations derived in this work and collected in Table 1. However, in the case of solution calculations (row 4), we detected discrepancies between the H̄^sol_s, TS̄^o,sol_s and μ^o,sol_s values in the G09 output (columns 4–6, row 5; values in parentheses) and the ones computed with the formulae in Table 1 (columns 4–6, row 4). Thermodynamic magnitudes corresponding to Ben-Naim's formulation (X̄*^,φ with φ = g, sol; columns 7–10) have been computed using formulas in Table 1 and the appropriate data from the G09 output (Ben-Naim's magnitudes are not provided explicitly in the G09 output). Let us analyze the origin of such discrepancies.

The expression employed in G09 to compute the translational component of the partial molar entropy of the solute in both gas-phase and solution calculations is¹⁹

with V^φ_m = 24.5 L mol⁻¹. That is to say, G09 uses the same value for V^φ_m in gas-phase and solution calculations. It is straightforward to show that it is equivalent to state that G09 employs concentration scale standard states for gas-phase and solution with identical concentrations, namely, ρ^g = ρ^sol_s = 1 bar/RT = 1/24.5 mol L⁻¹ = 0.041 M, in molar units [i.e. {g-ideal, ρ^g = 0.041 M; V_m = 24.5 L mol⁻¹} standard state for gas-phase and {sol-Henry, ρ^sol_s = 0.041 M; V_m = 24.5 L mol⁻¹} standard state for solution]. After adding the contributions from the internal degrees of freedom, the resulting equation for the partial molar entropy, S̄^⊕,φ_s, agrees with eqn (7) and (9), namely, the ones derived in this work for the gas-phase calculations. Then, the G09 value for TS̄^⊕,sol_s is computed in molar units as

As a consequence, the G09 values for TS̄^⊕,g (0.021436 hartree) and TS̄^⊕,sol_s (0.021441 hartree) are virtually identical (the negligible discrepancy observed, less than 2 × 10⁻³ kcal mol⁻¹, arises from the difference in the contributions from the internal degrees of freedom after geometry optimizations in each phase). It should be noted that the difference between the two mentioned entropy contributions as computed with eqn (19) (the equation derived in this work: TS̄^o,sol_s) and (60) (the equation employed by G09: TS̄^⊕,sol_s), respectively, is (molar units), RT + RT log(RT/V_l,m). The first term is associated with the error involved in the G09 incorrect use of the gas-phase eqn (59) and (60) for calculations in solution. The proper equation to be used in solution calculations is eqn (19) that involves (molar units) (with RT²α^l_p ≈ 0) instead of . The second term arises from the use of different standard states and consequently different values for the volume in eqn (19) [standard state: {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹}] and (60) [concentration scale standard state: {sol-Henry, ρ^sol_s = 0.041 M; V_m = 24.5 L mol⁻¹}].

Thus we can conclude that the TS̄^⊕,sol_s(G09) value overestimates by RT the correct magnitude. When properly computed with eqn (19), which is referred to the {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹} standard state, TS̄^o,sol_s becomes 0.013689 (see Table 2). If the {sol-Henry, ρ^sol_s = 0.041 M; V_m = 24.5 L mol⁻¹} standard state is chosen for the solution calculations (as G09 does), the correct TS̄^⊕,sol_s(G09) value [as computed with eqn (21)] would be: 0.021441 − RT = 0.020497 hartree. It must be emphasized that while the overestimation by RT is a mistake, the choice of any of the two mentioned standard states is free. Consequently, both values (0.013689 and 0.020497 hartree) are correct, provided the appropriate standard state is specified.

It is noteworthy that the TΔS^⊕,sol(G09) solvation magnitude is computed as [see eqn (60) and (7)]

which agrees with TΔS*^,sol = 0 [see eqn (49) and (50)]. Consequently, the G09 output provides gas-phase and solution partial molar entropies that, although do not correspond to S̄*^,φ (φ = g, sol) Ben-Naim's magnitudes (i.e. the values in column 9, rows 2 and 4, in Table 2 do not match the G09 output: column 5, rows 2 and 5), their difference, by a lucky fluke, leads to Ben-Naim's solvation entropies (i.e. TΔS^⊕,sol(G09) = TΔS*^,sol; see Table 3).

Regarding partial molar enthalpies, results collected in Table 2 show that G09 computes H̄^sol_s as

However, according to the equations derived in this work, H̄^sol_s should be computed from eqn (22) that, after neglecting the kT²α^l_p + PV̄^sol_s contribution (see Section 2), differs from eqn (62) in (molar units) RT (0.000944 hartree; see column 4, rows 4 and 5, in Table 2). That is, G09 overestimates H̄^sol_s in this amount. The mistake arises from the wrong tacit assumption in G09 that (molar units) kT²α^l_p + PV̄^sol_s = RT. As mentioned, in Sections 2 and 3, RT²α^l,water_p = 4.5 × 10⁻² kcal mol⁻¹, PV̄^sol_water = 4.4 × 10⁻⁴ kcal mol⁻¹ and RT = 0.59 kcal mol⁻¹, and hence kT²α^l_p + PV̄^sol_s cannot be identified with RT.

Thus, the G09 solvation enthalpy is computed as

which agrees with ΔH*^,sol, once the PV̄^sol_s term (4.4 × 10⁻⁴ kcal mol⁻¹) is neglected [see eqn (53)]. Therefore, like in the case of entropy, the G09 output provides partial molar enthalpies for gas-phase and solution calculations such that, although they are not H̄*^,φ (φ = g, sol) Ben-Naim's partial molar enthalpies (i.e. the values in column 8, rows 2 and 4 in Table 2 do not match the G09 output: column 4, rows 2 and 5), their difference leads to Ben-Naim's solvation enthalpies (i.e. ΔH^sol(G09) = ΔH*^,sol; see Table 3).

Let us now comment on G09 chemical potentials. The μ^⊕,sol_s(G09) value provided by G09 is “correct”, but it is again a fluke. Indeed, μ^⊕,sol_s(G09) is computed as H̄^sol_s(G09) − TS̄^⊕,sol_s(G09). According to the previous paragraphs, the RT incorrect contributions to H̄^sol_s(G09) and TS̄^sol_s(G09) values, will cancel each other out, thus leading to a correct μ^⊕,sol_s(G09) value in the G09 output (−76.423076 hartree) that, according to what has been mentioned above, will correspond to the {sol-Henry, ρ^sol_s = 0.041 M; V_m = 24.5 L mol⁻¹} concentration scale standard state. The value μ^o,sol_s = −76.416267 hartree (see Table 2) computed from eqn (18), neglecting the PV̄^sol_s term, is indeed also a correct value for the solute (water) chemical potential in solution (water), but this time, it corresponds to the {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹} standard state.

It has been stressed in Section 4 that Ben-Naim's solvation Gibbs free energy ΔG*^,sol is identical in magnitude (although conceptually different!) to the ΔG(ρ-process) value. As noted there, this latter, is the energy involved in the transference of a solute from the gas-phase {g-ideal, ρ^g; V_m} standard state to the solution {sol-Henry, ρ^sol_s; V_m} standard state under the condition that number densities in both phases are identical (ρ^g = ρ^sol_s). This is exactly the situation considered in G09 where the number densities employed in both phases are 0.041 M (molar units). Thus, according to what has been established in Section 4, ΔG(ρ-process) coincides in magnitude with ΔG*^,sol, Ben-Naim's solvation Gibbs free energy. Consequently, the solvation Gibbs free energy computed in G09, ΔG^⊕,sol(G09) = μ^⊕,sol_s(G09) − μ^⊕,g(G09) = −76.423076 + 76.416029 = −0.007047 hartree = −4.4 kcal mol⁻¹, coincides with Ben-Naim's solvation Gibbs free energy ΔG*^,sol = μ^*,sol_s − μ*^,g = −76.408996 + 76.401950 = −0.007046 hartree. We have used the fact, mentioned in Section 2.1, that since in G09, ρ^g = 0.041 M, then μ^⊕,g(G09) = μ^o,g. Thus, likewise for entropies and enthalpies, the G09 output provides chemical potentials for gas-phase and solution calculations such that, although they are not μ*^,φ (φ = g, sol) Ben-Naim's pseudochemical potentials¹ (i.e. the values in column 10, rows 2 and 4 in Table 2 do not match the G09 output: column 6, rows 2 and 5), their difference leads to Ben-Naim's solvation Gibbs free energies (i.e. ΔG^⊕,sol(G09) = ΔG*^,sol; see Table 3).

It is very instructive and conceptually clarifying to emphasize that the ΔG^o,sol solvation Gibbs free energy computed within Fowler–Guggenheim's formalism (transference of a solute from the gas-phase {g-ideal, P = 1 bar; V_m = 24.5 L mol⁻¹} standard state to the solution {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹} standard state (see Table 3), i.e. (see Table 2) ΔG^o,sol = μ^o,sol_s − μ^o,g = −76.416267 + 76.416029 = −0.000238 hartree = −0.1 kcal mol⁻¹, differs, in 4.3 kcal mol⁻¹ from the corresponding value in Ben-Naim's formulation, i.e. (see Table 2) ΔG*^,sol = μ^*,sol_s − μ*^,g = −76.408996 + 76.401950 = −0.007046 hartree = −4.4 kcal mol⁻¹. These 4.3 kcal mol⁻¹, that according to the previous paragraph correspond to the difference ΔG^o,sol − ΔG(ρ-process), must be ascribed to the difference in Gibbs free energy between the two following transferences: (a) ΔG^o,sol: {g-ideal, P = 1 bar; V_m = 24.5 L mol⁻¹} → {sol-Henry, x_s → 1; V_m = 0.018 L mol⁻¹} and (b) ΔG(ρ-process): {g-ideal, ρ^g; V_m} → {sol-Henry, ρ^sol_s; V_m}, with ρ^g = ρ^sol_s. Indeed, Cabani et al. mention in their compilation⁹ that the above difference (4.3 kcal mol⁻¹) holds for the particular case where both number densities (gas-phase and solution) become 1 M in the ρ-process (as stressed above, in the case of G09, both number densities are also identical, but this time their respective values are 0.041 M).

Before ending, we think it is important to remark that G09 does not compute any Ben-Naim's thermodynamic quantity X̄* (the values in the columns 7–10 in Table 2 have been computed by us from data in the G09 output, but they are not found explicitly in it). Furthermore, as concluded in this article, G09 provides wrong values for the partial molar enthalpies and entropies of a solute in solution. It is hoped that the present contribution can be helpful in order to reorganize the solvation thermodynamics section of the Gaussian packages of programs as well as of any other software where the solvation code may call for corrections. Despite the apparent present superiority of Ben-Naim's formulation over the conventional approach, we do suggest that an optimum output should provide detailed information regarding the two formulations considered in this article.

6. Conclusion

A parallel derivation of conventional (Fowler–Guggenheim) and Ben-Naim's formulations of solvation thermodynamics, emphasizing their differences and stressing their interconnections, is presented. Explicit equations for conventional (ΔX^o,sol = X̄^sol_s − X̄^g) and Ben-Naim's (ΔX*^,sol = X̄^*,sol_s − X̄*^,g) solvation magnitudes, with X being internal energies, enthalpies, entropies, and Gibbs free energies, have been obtained. Pivotal equations connecting both formulations have also been inferred.

Both theoretical frameworks were used to compute the solvation thermodynamic magnitudes for a standard system: water solute in water solution. Comparison with the Gaussian 09 output shows that this software does not estimate properly neither the partial molar entropies nor the partial molar enthalpies of a solute in solution.

The present article intends to be of help to researchers who need to apply solvation models in order to tackle the study of chemical, biophysical or biochemical systems which they may be interested in, as well as to software developers who implement the solvation thermodynamic equations in the available codes.

Appendix

The first term contribution to the Helmholtz free energy of solution (F^sol) in eqn (12)⁶ is

F_l = F_l(T, V^sol − N_sV̄^sol_s, N_l) = F_l[T, V^sol_l(V^sol, N_s, V̄^sol_s), N_l] (A1)

where V^sol_l (= N_lV̄^sol_l) is the volume available for the N_l solvent molecules

V^sol_l = V^sol − N_sV̄^sol_s (A2)

Let us use in the following V, V_l and V̄_s, instead of V^sol, V^sol_l and V̄^sol_s, for the sake of simplicity.

From (A1) and (A2) we get

andBut

On the other hand, for an ideal dilute solution N_s → 0 and V_l → V. Then,

where, in order to write the last equality of (A6), we have taken into account that −(∂F_l/∂V_l)_T,Nl ≈ −(∂F_l/∂V)_T,Nl is the pressure at which the pure solvent has the volume V − N_sV̄_s. As stressed by Fowler and Guggenheim,⁶ such a pressure will be quite close to the pressure P on the actual solute–solvent assembly.

Consequently, bearing in mind that according to (A2), (∂V_l/∂N_s)_V,V̄s = −V̄_s, (A4) leads to

in agreement with Fowler and Guggenheim [see eqn (823, 3 in p. 373 of ref. 6].

Acknowledgements

Profs J. Fernández-Rico, J. M. García de la Vega, A. Largo, J. M. Lluch and E. Ortí were kept abreast of the present research at different stages of its development. I wish to thank their stimulating comments. The author is indebted to Ms María Viña García-Bericua for her help with the artwork.

References

1. A. Ben-Naim, Solvation Thermodynamics, Plenum, NY, 1987.
2. See for example, A. Varela-Álvarez, D. Markovic, P. Vogel and J. A. Sordo, J. Am. Chem. Soc., 2009, 131, 9547–9561.
3. See for example, Solvation Effects in Molecules and Biomolecules. Computational Methods and Applications, S. Canuto, Springer, 2009.
4. See for example, Drug Design. Structure and Ligand-Based Approaches, ed. K. M. Merz jr, D. Ringe and C. H. Reynolds, Cambridge Univ. Press, 2010.
5. J. A. Sordo, Anal. Bioanal. Chem., 2014, 406, 1825–1828.
6. R. Fowler and E. A. Guggenheim, Statistical Thermodynamics, Cambridge Univ. Press, Cambridge, 1956.
7. A. J. Ben-Naim, Phys. Chem., 1978, 82, 792–803.
8. R. A. Pierotti, Chem. Rev., 1976, 76, 717–726.
9. S. Cabani, P. Gianni, V. Mollica and L. Lepari, J. Solution Chem., 1981, 10, 563–595.
10. M. J. Frisch, et al., Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford CT, 2009.
11. A. Ben-Naim and Y. J. Marcus, Chem. Phys., 1984, 81, 2016–2027.
12. A. Ben-Naim, Water and Aqueous Solutions. Introduction to a Molecular Theory, Plenum, NY, 1974.
13. A. Ben-Naim, Statistical Thermodynamics for Chemists and Biochemists, Plenum, NY, 1992.
14. M. H. Abraham and A. J. Nasehzadeh, J. Chem. Soc., Faraday Trans. 1, 1981, 77, 321–339.
15. A. Ben-Naim, J. Solution Chem., 2001, 30, 475–487.
16. E. Wilhelm and R. Battino, Chem. Rev., 1973, 73, 1–9.
17. E. Wilhelm, R. Battino and R. J. Wilcock, Chem. Rev., 1977, 77, 219–262.
18. W. L. Jorgensen, J. Gao and C. Ravimohan, J. Phys. Chem., 1985, 89, 3470–3473.
19. J. W. Ochterski, Thermochemistry in Gaussian, http://www.gaussian.com, 2000.
20. M. H. Abraham, G. S. Whiting, R. Fuchs and E. J. Chambers, J. Chem. Soc., Perkin Trans. 1, 1990, 2, 291–330.

---