Prediction of Hydrogen-Bonding Interaction Free-Energies with Two New Molecular Descriptors

1. Introduction

The role of hydrogen bonding (HB) in numerous physicochemical processes in biology and life itself, in drug and xenobiotic interaction with biota, in aquatic environments and in chemical industry cannot be overemphasized [2,3,4,5,6,7,8,9,10,11,12]. HB is studied experimentally and theoretically for more than one century now [13]. It is one of the principal causes of mixture non-idealities and is of particular interest in modern molecular thermodynamic developments. The current status in the open literature regarding the estimation of the strength of hydrogen bonding interactions and the open problems are extensively discussed in numerous monographs, reviews and relevant articles [1,2,3,4,9,10,11,14,15,16,17].

At the outset it should be stressed that HB interaction energies and free-energies are obtained indirectly via modelling their contribution to measurable thermodynamic properties and, thus, they are attributed a quasi-thermodynamic character. One of the most exploited thermodynamic properties for extracting HB information is solvation free energy and the associated solute partitioning / transfer and the limited Henry’s law. The solvation free-energy, $\Delta G_{12}^S$, and its components enthalpy, $\Delta H_{12}^S$ and entropy, $\Delta S_{12}^S$, are directly connected with phase equilibrium studies involving activity or fugacity coefficients through the following most useful classical working equation for the equilibrium solvation constant $K_G^S$:

$$2.303Log{K}_G^S={\displaystyle \frac{\Delta G_{12}^S}{RT}}={\displaystyle \frac{\Delta H_{12}^S-T\Delta S_{12}^S}{RT}}=\ln {\displaystyle \frac{{\phi }_1^0{P}_1^0{V}_{m2}{\gamma }_{1/2}^{\infty }}{RT}}$$

(1)

V_m2 in Equation (1) is the molar volume of the solvent (component 2) and ${\gamma }_{1/2}^{\infty }$is the activity coefficient of solute 1 at infinite dilution in solvent 2. $P_1^0$is the vapor pressure of pure solute at temperature, T, and ${\phi }_1^0$ its fugacity coefficient (typically, set equal to 1 at ambient conditions).

It is from the assumed contribution to quantities like the above solvation free energy and enthalpy that HB interaction energies (enthalpies), $\Delta H_{12}^{hb}$, and free energies, $\Delta G_{12}^{hb}$, are very often obtained. In this regard, many attempts have been made to model and quantify HB interactions for use in quantitative structure-property and structure - activity relationships (QSPR / QSAR) and in the development of molecular thermodynamic models of highly non-ideal mixtures [11,12,14,15,16,17,18,19,20,21,22,23,24,25,26,27].

One of the most successful QSPR-type approaches is Abraham’s LSER (Linear Solvation Energy Relationship) model [28,29,30,31,32,33] and has led, ultimately, to the valuable and freely available comprehensive and rich in thermodynamic information LSER Database [34].

The LSER model uses simple linearity equations for the quantitative description of solute transfer between two phases. These important linear relationships, in Abraham’s LSER approach, take the form of Equations (2) and (3) for the equilibrium constant, $K_G^S$, of solute partitioning between gas and liquid phases and for the corresponding solvation energy (enthalpy) constant, $K_E^S$, respectively:

$$Log{K}_G^S={\displaystyle \frac{-\Delta G_{12}}{2.303RT}}=c_{g2}+e_{g2}{E}_1+s_{g2}{S}_1+a_{g2}{A}_1+b_{g2}{B}_1+l_{g2}{L}_1\hspace{1em}\left(\mathrm{LSER}\right)$$

(2)

and

$$Log{K}_E^S={\displaystyle \frac{-\Delta H_{12}}{2.303RT}}=c_{e2}+e_{e2}{E}_1+s_{e2}{S}_1+a_{e2}{A}_1+b_{e2}{B}_1+l_{e2}{L}_1\hspace{1em}\left(\mathrm{LSER}\right)$$

(3)

The upper-case letters in these equations indicate solute molecular LSER descriptors while the lower-case letters indicate the corresponding complementary but solvent-phase-specific LFER coefficients. The solute molecular LSER descriptors V_x, L, E, S, A, and B, are corresponding to the McGowan’s characteristic volume, the equilibrium constant for the gas–liquid partition in n-hexadecane at 298 K, the excess molar refraction, the dipolarity / polarizability, the hydrogen bonding (HB) acidity A, and hydrogen-bonding basicity B, respectively [28,29,30,31,32,33]. Thus, in Abraham’s LSER model, the HB contribution to solvation energy (enthalpy) constant of a solute 1 by solvent 2 is modelled as the sum, a_e2A₁ + b_e2B₁. The corresponding HB contribution to the equilibrium constant for the solvation free-energy is modelled as the sum a_g2A₁ + b_g2B₁. Analogous to the above equations apply for the solute transfer / partitioning between two condensed phases.

There are three main limitations of this otherwise highly successful LSER model. The descriptors A and B and the coefficients a and b are essentially obtained from extensive experimental data correlations and depend on their availability. In fact, for just a few decades of solvents are a and b coefficients available [17,34,35]. The second limitation arises from the very way this LSER model treats hydrogen bonding. On self-solvation, where solute and solvent become identical, one would expect the acid – base (aA) interaction to be identical to the base – acid (bB) interaction between the very same pair of donor – acceptor sites. However, in LSER the product aA is in general different from bB [40,41,42]. In fact, this is the case for other popular QSPR – type models, including Raevski’s HYBOT model [36]. This restricts rather severely the transfer of this HB information into other models, especially in Molecular Thermodynamics modeling of hydrogen bonding. The third limitation arises from the fact that the LFER coefficients, a and b, are determined simultaneously with all other LFER coefficients, including the constant c, by multilinear regression of available experimental solvation data. This implies that some caution must be exercised when attributing to the above sums, a_e2A₁ + b_e2B₁ and a_g2A₁ + b_g2B₁ sums the (exclusive or total) HB contribution to solvation energy and free energy, respectively.

Recently [1], a predictive method was developed for the HB contribution to solvation energy. The method is using new molecular descriptors [37] based on quantum – chemical (QC) calculations and has a simple structure analogous to Abraham’s LSER model [28,29,30,31,32,33]. The novel set of molecular descriptors, referred to as QC-LSER descriptors, is based on the molecular surface charge densities or σ – profiles widely used by the quantum-mechanics based COSMO-RS model [38,4]. These σ- profiles are available, free of charge, for thousands of molecules in the open literature, as an example in ref. [45]. They may, of course, be obtained also by using appropriate quantum-chemical calculation suites, such as the TURBOMOLE, DMol³ of BIOVIA’s MATERIALS STUDIO suite, or the SCM suite [46,47,48]. In our previous work [1], the σ-profiles from COSMObase [44] at the DFT / TZPVD-Fine level of quantum chemical calculations are used and the corresponding QC- LSER descriptors for a number of common hydrogen bonded solutes were reported.

Two of the QC-LSER descriptors, the HB acidity, A_h, and the HB basicity, B_h, were particularly explored for the development of the predictive scheme due to their sound basis and insightful character [37]. In fact, the effective HB acidity descriptor of a solute was given by the product α = f_AA_h and the effective HB basicity descriptor by the product β = f_BB_h. It was observed that the “availability fractions” f_A and f_B are characteristics of homologous series (have the same value for all solutes of the homologous series). The factors f_A and f_B and the descriptors α and β were reported for a number of common hydrogen-bonded solutes [1]. With these descriptors, the HB interaction energy for a solute 1 – solvent 2 pair may be obtained at 25 °C by the simple equation:

$$\begin{array}{l}-\Delta E_{12}^{hb}=2.303RT\left({\alpha }_1{\beta }_2+{\beta }_1{\alpha }_2\right)=2.303\times 8.314\times 298.15\left({\alpha }_1{\beta }_2+{\beta }_1{\alpha }_2\right)=\\ \hspace{1em}\hspace{1em}\hspace{1em}=5.71\left({\alpha }_1{\beta }_2+{\beta }_1{\alpha }_2\right) \hspace{0.33em}kJ/mol\end{array}$$

(4)

The method was tested against LSER data (a_e2A₁ + b_e2B₁) but also against corresponding estimations of $\Delta E_{12}^{hb}$from COSMO-RS model [44]. In nearly all cases the predictions were close to one or both of these sets of data.

The present work is attempting to extend this previous method [1] to a simple and versatile method for the prediction of the HB contribution to solvation free-energy, $\Delta G_{12}^{hb}$, of solute 1 in solvent 2. This task, however, is not trivial. In contrast to $\Delta E_{12}^{hb}$, COSMO-RS does not provide estimations of $\Delta G_{12}^{hb}$and there are good reasons for this: Besides the energy (enthalpy) component, $\Delta G_{12}^{hb}$has also an entropy component which is not easily amenable to a priori estimation in complex hydrogen-bonded solvent environments. Thus, the developments of the present work will be confined and exclusively based on the corresponding Abraham’s LSER estimations, a_g2A₁ + b_g2B₁. It is pointed out again that LSER is currently one of the best predictive models for solvation free energies. The goal of the present work is to attribute effective acidity and basicity descriptors, α_G1 and β_G1 to each solute, which will enable us to predict their HB contribution to solvation free energy in a solvent with effective acidity and basicity descriptors, α_G2 and β_G2, by an equation analogous to Equation (4), or:

$$-\Delta G_{12}^{hb}=5.71\left({\alpha }_{G1}{\beta }_{G2}+{\beta }_{G1}{\alpha }_{G2}\right) \hspace{0.33em}kJ/mol$$

(5)

In this effort, one of the major issues that must be addressed, is the fact that in LSER calculations the quantity $\Delta G_{12}^{hb}$is often solvent – specific, especially in multi-sited solvents (molecules with many and/or distant HB sites) or $\Delta G_{12}^{hb}\ne \Delta G_{21}^{hb}$. In contrast to $\Delta E_{12}^{hb}$ the LSER model may provide estimations of $\Delta G_{12}^{hb}$for many more hydrogen-bonded mixtures, including multi-sited solvent systems. This will help reveal the solvent-specific character of$\Delta G_{12}^{hb}$ in these systems.

The fact that our calculations are heavily based on the corresponding LSER calculations reorients also the central scope of the present work: The LSER model reproduces successfully the experimental solvation data for some decades of hydrogen-bonded solvent systems. The developments in the present work aim, ultimately, at complementing and extending the LSER method to many more solvent systems by removing its above-mentioned drawbacks. The rationale and the structure of the new predictive method are presented in the next section.

2. The Predictive Method for HB Interaction Free Energies

As mentioned above, the developments in the present work will be in line with the corresponding developments in our previous work on HB interaction energies [1]. The effective acidity and basicity descriptors will be obtained by, first, determining the corresponding “availability fractions” f_AG and f_BG, which will give:

$$\begin{array}{l}{\alpha }_G=f_{AG}{A}_h\\ {\beta }_G=f_{BG}{B}_h\end{array}$$

(6)

An updated list of QC-LSER descriptors, including A_h and B_h descriptors, are reported in Table SI1 of the Supplementary Information SI1 file. Thus, the first step is the determination of fractions f_AG and f_BG, mainly, from the corresponding LSER estimation, a_g2A₁ + b_g2B₁ [28,29,30,31,32,33] and along the lines of our previous work [1]. In this first step we will focus on relatively simple systems. Complex multi-sited solvent systems, exhibiting solvent – specific peculiarities, will be examined later.

In order to make the presentation of the method clear, let us consider two self-associated compounds, say ethanol (1) and 2-pyrrolidone (2) with corresponding acidity descriptors ${\alpha }_{G1}$ and ${\alpha }_{G2}$, and basicity descriptors β_G1 and β_G2, respectively. By adapting the rationale of the previous work on HB interaction energies [1] to the present case of HB interaction free-energies, we may write for the pair 1 – 2 HB interaction free energy at 298.15 K the following equation:

$$-\Delta G_{12}^{hb}=2.303RT\left({\alpha }_{G1}{\beta }_{G2}+{\beta }_{G1}{\alpha }_{G2}\right)=5.71\left({\alpha }_{G1}{\beta }_{G2}+{\beta }_{G1}{\alpha }_{G2}\right) \hspace{0.33em}kJ/mol$$

(7)

The corresponding equations for their self-association interaction in a self-solvation process are:

$$-\Delta G_{11}^{hb}=11.42\left({\alpha }_{G1}{\beta }_{G1}\right) \hspace{0.33em}$$

(8)

$$-\Delta G_{22}^{hb}=11.42\left({\alpha }_{G2}{\beta }_{G2}\right) \hspace{0.33em}$$

(9)

Let us consider also their interaction with a proton – donor compound 3 (say, chloroform), which cannot self-associate but it can form hydrogen bonds with both compounds 1 and 2. The corresponding equations for their HB interaction free energy are:

$$-\Delta G_{13}^{hb}=5.71\left({\alpha }_{G3}{\beta }_{G1}\right) \hspace{0.33em}$$

(10)

$$-\Delta G_{23}^{hb}=5.71\left({\alpha }_{G3}{\beta }_{G2}\right) \hspace{0.33em}$$

(11)

Equations (7) to (11) are five independent equations with five unknowns. If all interaction free energies are known, they may be solved, in principle, and give ${\alpha }_{G1}$, ${\alpha }_{G2}$, β_G1, β_G2, and ${\alpha }_{G3}$.

Alternatively, or in addition, one may use a proton-acceptor compound 4 (say, ethyl acetate), which cannot self-associate but it can form hydrogen bonds with both compounds 1 and 2. The corresponding equations for their HB interaction free energy are now:

$$-\Delta G_{14}^{hb}=5.71\left({\alpha }_{G1}{\beta }_{G4}\right) \hspace{0.33em}$$

(12)

$$-\Delta G_{24}^{hb}=5.71\left({\alpha }_{G2}{\beta }_{G4}\right) \hspace{0.33em}$$

(13)

Equations (7)–(9), (12) and (13) are five independent equations with five unknowns. If all interaction free energies are known, they may be solved, in principle, and give ${\alpha }_{G1}$, ${\alpha }_{G2}$, β_G1, β_G2, and β_G4. In principle, the first four descriptors would be identical to the ones obtained above with compound 3. Alternatively, or in addition, one could use another proton acceptor compound 5 (say, tetrahydrofuran) and obtain ${\alpha }_{G1}$, ${\alpha }_{G2}$, β_G1, β_G2, and β_G5.

The above procedure could be repeated by substituting compound 1 with another compound 6 from the same homologous series (say, 1-butanol in place of ethanol). Solution of the set of the new independent equations would give the descriptors ${\alpha }_{G6}$and β_G6. What is important now is to closely examine the descriptors of compounds 1 and 6 of the same homologous series (1-alkanols, in our case).

Following our previous practice [1], we may write the following equations for compound 1:

$$\begin{array}{l}{\alpha }_{G1}=f_{AG1}{A}_{h1}\hspace{0.33em}\\ {\beta }_{G1}=f_{BG1}{B}_{h1}\end{array}$$

(6a)

and, similarly, for compound 6:

$$\begin{array}{l}{\alpha }_{G6}=f_{GA6}{A}_{h6}\hspace{0.33em}\\ {\beta }_{G6}=f_{GB6}{B}_{h6}\end{array}$$

(6b)

As mentioned already, the coefficients f_AGi and f_BGi in the above equations are considered to reflect the effective fractions of the acidity and basicity sites of compound i which are available for hydrogen bonding. For a given HB functional group (-OH, in our case).it is rather reasonable to expect that these effective fractions are close, if not identical, for all members of the homologous series. Thus, the main focus of the above calculation scheme is on the determination of the common average values of fractions f_AGi and f_BGi in the homologous series of compounds i. For this purpose, many representative interacting compounds were used for each examined compound of the various homologous series and the determined descriptors are reported in Table 1. It is recalled again that the above calculation scheme was heavily based on the estimation of $\Delta G_{ij}^{hb}$by Abraham’s LSER method [28,29,30,31,32,33], or:

$$-\Delta G_{ij}^{hb}=2.303RT\left(A_i{a}_j+B_i{b}_j\right)=5.71\left(A_i{a}_j+B_i{b}_j\right)\hspace{0.33em}kJ/mol\hspace{1em}(\mathrm{at} \mathrm{298.15} K)$$

(2)

3. Results

In this section we will compare the calculations by Equation (7), using the molecular descriptors of Table 1, with the LSER calculations [17] by Equation (14). The detailed results are assembled in six tables, SII1 to SII6, and are reported in Supplementary Information file SI2. The full list of results with all examined hydrogen-bonded solutes are reported in each of these tables.

Table SII1 contains the results for the solvation in methanol, ethanol and 1-octanol solvents. All three solvents, members of the homologous series of 1-alkanols, are characterized by one single set of availability fractions, namely, f_AGi = 0.54 and f_BGi = 0.65. The results for solvation in 1-octanol solvent might be improved by using a different set with lower availability fractions but, following our previous practice [1], we have refrained from introducing such refinements at this stage of development of the predictive method. What is essential to notice is that, in this case of relatively simple HB functional groups, the equality $\Delta G_{12}^{hb}=\Delta G_{21}^{hb}$ holds true.

As observed in Table SII1, in the overwhelming majority of cases, especially in methanol and ethanol solvents, the discrepancies between predicted and LSER estimations of $\Delta G_{12}^{hb}$are less than 0.5 kcal/mol (or 2 kJ/mol). The peculiarities with phenols, which were noticed also in our previous work [1], seem to persist in HB contributions to solvation free energy as well. As previously [1], we have refrained from changing their availability fractions until it is clarified that their deviations are not due to rather excessively high HB LSER molecular descriptors.

For the bulk of solutes, only their most stable conformers C₀ are reported in Table SII1 and in the rest of tables in SI2 file. For some solutes, however, prone to intramolecular hydrogen – bonding, such as ethoxy-alkanols, more extended conformers are also reported. It is clear from Table SII1 that these conformers, when used with the availability fractions of conformers C₀, lead to higher deviations from LSER estimations. This means that they must be attributed their own set of availability fractions if selected to replace conformers C₀ in the $\Delta G_{12}^{hb}$ predictions. Alternatively, as discussed previously [1], relatively large deviations of predictions with C₀ conformers may imply that additional solute conformers must be used in order to properly represent its solvation behavior in the studied solvent.

Table SII2 contains the results for solvation in acetic acid, formamide and 2-pyrrolidone solvents. There are two points to be made regarding acetic acid solvent. The selected values for the availability fractions give rather satisfactory results for solvent or solute acetic acid but not as good for its self-solvation. It is not quite clear how the acid dimerization could explain this discrepancy [49]. The second point regards the solvation of 1-alkanols in acetic acid. As observed, the predicted results are higher than the corresponding LSER estimations. But, tt is also observed that the LSER estimations for the case of solvation of acetic acid in 1-alkanols (cf. Table S2) are significantly higher than the corresponding LSER estimations for the solvation of 1-alkanols in acetic acid. It is not quite clear if the tendency for esterification could explain this peculiarity. It should be mentioned, however, that this discrepancy between LSER estimations was also observed previously [1] in the case of HB contributions to solvation enthalpies but it was not in agreement with corresponding COSMO-RS estimations [44]. It should be clarified that the present predictive method accounts for plain HB interactions. Peculiarities of cooperativities, reactions, ionizations or cooperativities [12] are beyond the scope of the method at this stage of its development.

The predicted results for the case of solvation in formamide are rather satisfactory for the bulk of the studied solutes. The selected values for the availability fractions give also satisfactory predictions for the solvation of formamide in various solvents as well as for its self-solvation.

The results for the solvation in 2-pyrrolidone deserve some discussion. The availability fractions for 2-pyrrolidone are (cf. Table 1): f_AG = 0.35 and f_BG = 0.38. As observed in Table SII2, for nearly all heterosolvated (non-self-associated) proton acceptor solutes, like ethers, esters, ketones, etc. the predicted results are higher than the corresponding LSER estimations. One could consider this discrepancy as an indication of a rather exaggerated high proton-donor capacity of 2-pyrrolidone or an assumed high value for the availability fraction f_AG. By lowering its value from 0.35 to a value below 0.15, indeed, satisfactory results could be obtained for these heterosolvated compounds. However, this lowering of f_AG would imply a significant increase in f_BG above 0.85 in order to retain self-solvation result in satisfactory agreement with LSER estimations. But this alternative pair of availability fractions would lead to rather very high discrepancies for nearly all other studied solutes. Thus, the source of this discrepancy for the heterosolvated solutes is not clear. As in the case of 1-alkanol solvents, the solvation of phenols and carboxylic acids in 2-pyrrolidone exhibits the same type of discrepancies, again, of not clear origin. The case of phenolic / aromatic solutes will be discussed in the next section.

The results for solvation in aniline are reported in Table SII3 along with the results for solvation in chloroform and ethyl acetate. Besides acidity and basicity character, aniline has also an aromatic character. The availability fractions, heavily dictated by its self-solvation, are: f_AG = 0.45 and f_BG = 0.95. The rather high value of f_BG underlines its basic character. For the bulk of solutes, the predictions are in rather satisfactory agreement with LSER estimations. The agreement for carboxylic acids is also satisfactory. Phenols are again an exception. Their discrepancy could be removed by increasing their f_BG, but, as argued above, we have refrained from doing so in this work. The solvation results for the solute aniline in the studied solvents are also satisfactory (cf. Tables SII1–SII4).

Chloroform is a proton donor and cannot self-associate (heterosolvated compound). As seen in Table SII3, the availability fraction f_AG = 0.37, leads to predictions which are nearly always higher (lower in absolute terms) than the corresponding LSER estimations. One would be tempted, then, to increase f_AG in order to improve results. This, however, would deteriorate the results for the solvation of the solute chloroform in the various other solvents and, thus, we have refrained from doing so in this work. In fact, the source of the observed discrepancies may be due to the fact that part of the HB contribution to solvation free energy may be present in the constant LFER coefficient of the LSER model or in some of the other terms of the LSER Equation (1). If this is the case, better agreement with LSER estimations could be reached by simply adding a constant term of the order of 1.2 in the predictions by Equation (7).

Ethyl acetate is a proton acceptor and, like chloroform, cannot self-associate. As seen in Table SII3, the availability fraction f_BG = 0.49 leads to predictions in rather satisfactory agreement with LSER estimations for nearly all studied solutes. Similar is the picture with the other proton acceptor heterosolvated solutes like butyl acetate, tetrahydrofuran, diethyl ether, butanone, etc., as one can easily verify by just using Equation (7) with the α_G and β_G descriptors of Table 1.

Of much interest, however, are the results reported in Table SII4 for the solvation in diethylene and triethylene glycols. As seen, for each of these solvents, two alternative predictions are reported based on the two alternative sets of their α_G and β_G descriptors reported in Table 1. The very molecular structure of these solvents differentiates their solvent behavior from their solute one, as will be discussed in the next section. For the time being let us focus on the results in Table SII4.

As observed in Table SII4, the predictions based on the “solute” descriptors of the solvent are significantly lower compared to the predictions based on the “solvent-specific” descriptors of the solvent. This is more pronounced for triethylene glycol having the bigger multi-sited molecular structure. As seen in Table SII4, the average deviation in the case of triethylene glycol, using “solute” descriptors, is 8.99 kJ/mol and, when using the “solvent – specific” descriptors, the average deviation is reduced to 0.08 kJ/mol. In the case of diethylene glycol, the average deviation, using the “solute” descriptors, is 5.08 kJ/mol and is reduced to 0.39 kJ/mol when using the “solvent-specific” α_G and β_G descriptors reported in Table 1. For both glycols, the predictions with the solvent – specific α_G and β_G descriptors are in rather satisfactory agreement with the corresponding LSER estimations for nearly all studied solutes. Phenols are, again, exceptions with the predictions being lower (in absolute terms) than the LSER estimations, as with all other studied solvents.

Similar comments apply to the results reported in Table SII5 for solvation in solvents ethylene glycol and 1,2 – propylene glycol. In the case of ethylene glycol, the above average deviations are 1.63 kJ/mol with the “solute” descriptors and 0.23 with the solvent – specific descriptors. In the case of 1,2 – propylene glycol, these average deviations are 3.08 kJ/mol and 0.96 kJ/mol, respectively.

In contrast to the organic solvents discussed so far, the solvation in water presents numerous peculiarities and satisfactory agreement of predictions with LSER estimations is rather difficult to obtain with just one solvent – specific set of α_G and β_G descriptors. By using the “solute” descriptors of water, the deviations are often large, as seen in Table SII6. In nearly all cases, the predictions are lower (in absolute terms) than the corresponding LSER estimations. Part of these deviations might be attributed to the relatively high constant term of Equation (1) in the case of solvation in water [17]. But still the observed discrepancies elude such a simple explanation. The picture that emerges from Table SII6 is pointing to a rather solute – specific source of discrepancy in these aqueous systems. Thus, in Table SII6 are also reported the changes in f_BG of solutes required to bring agreement between predictions and LSER estimations. As seen, these changes must often be well above 100%. An attempt to explain these peculiar discrepancies is made in the next section.

4. Discussion

The results presented in the previous section indicate that the prediction of $\Delta G_{12}^{hb}$ cannot be made with one single set of α_G and β_G descriptors for all solvents. One single set of α_G and β_G descriptors (the same set for the compound as solute and as solvent) seems to be sufficient for the organic solvents of relatively simple molecular HB structure, like the solvents in Tables SII1 to SII3 in SI2 file. Organic solvents with more complex molecular HB structure, like the solvents of Tables SII4 and SII5 of SI2 file require two separate sets of α_G and β_G descriptors, one to account for their HB behavior as solvents and one for their HB capacity as solutes. Water is a stand-alone case in which the above predictive scheme does not seem to apply. The key question is why the prediction of$\Delta G_{12}^{hb}$ in these solvents require more than one set of α_G and β_G descriptors. The answer to this question is not trivial but some hints could be drawn from an inspection of the molecular structure of these solvents.

In Figure 1 are shown the geometries and COSMO-surfaces [44] of three representative conformers of triethylene glycol. The open structure of conformer C0 exhibits its two proton donor (acidic) sites as well as its four proton acceptor (basicity) sites. The other two relatively closed structures of conformers C3 and C9 exhibit the basicity sites but nearly hides the acidity sites because they are involved in intramolecular hydrogen bonding. When triethylene glycol is a solute in a relatively simple hydrogen-bonded solvent like ethanol, a number of solvent molecules may hydrogen - bond with one solute molecule and this brings up much of its hydrogen bonding capacity translated to relatively high availability fractions for the solute molecule. In contrast, if triethylene glycol is used as solvent, then, one simple solute molecule like ethanol may form, at most, two hydrogen bonds with either one or two solvent molecules. In either case this solute brings up part only of the hydrogen bonding capacity of the solvent triethylene glycol. Thus, when solvent, the availability fractions of triethylene glycol are lower than the ones corresponding to solute triethylene glycol.

Similar arguments apply to the other HB multi-sited solvents and could explain the use of separate sets of availability fractions, one as solvents and one as solutes. In the particular case of intramolecularly hydrogen-bonded solutes, like triethylene glycol (cf. Figure 1), their prevailing conformation may change with composition if this lowers the overall system free energy implying that, at the two ends of mixture composition, molecules like triethylene glycol may adopt different prevailing conformations and, thus, different hydrogen bonding capacities.

The above arguments, however, do not seem sufficient to explain the LSER estimations of$\Delta G_{12}^{hb}$ in solvent water. The peculiarity of water is its small molecular size exposing two proton donor sites and its two basic lone-pair sites. Its small size permits water molecule to probe most of the hydrogen bonding capacity of the interacting molecule. In fact, the results in Table SII6 indicate that nearly all solute molecules expose much more of their HB capacity in water solvent than in any other organic solvent. In order to get an idea of this additional exposure of HB solute capacity, in Table SII6 are reported the availability fractions f_BG which make predictions equal to the LSER $\Delta G_{12}^{hb}$estimations. The required percent increase in f_BG is shown in the last column of Table SII6. As shown, even simple alkanols exhibit about 8% more hydrogen bonding capacity when solvated by water. Apparently, water molecule may reach both lone pairs of oxygen of the hydroxyl groups and two water molecules may simultaneously interact with hydroxylic oxygen of one alkanol molecule. In other more complex molecules the percent increase in f_BG is often well above 100%, as shown in Table SII6. Such large increase in f_BG is also required for phenols and this may indicate that there is much more basicity capacity in these phenolic solutes than one would expect from just their hydroxylic oxygen and this requires now some attention.

In contrast to COSMO-RS model [38,39,40,41,42,43], the LSER model [28,29,30,31,32,33] attributes some basic character (proton acceptor capacity) to aromatic hydrocarbons, like benzene or toluene. Obviously, this basicity is a feature of the aromatic ring and is manifested in the HB interactions of all aromatic solutes. In fact, by increasing the β_G descriptor of aromatic solutes by 0.5, the predictions of $\Delta G_{12}^{hb}$are brought in rather satisfactory agreement with the corresponding LSER estimations, as shown in Table 2.

If LSER estimations reflect more accurately the basicity character of aromatic solutes than the COSMO-RS model, the question is why the latter fails. It is recalled that the cutoff surface charge density for a surface site to be considered as basic is 0.01 e/Ǻ². Screen charge densities lower than this cutoff value are not considered sufficient to attribute a proton acceptor capacity to their sites. In the case of benzene, there are no surface sites with surface charge densities higher than 0.008 e/Ǻ². COSMO-RS could attribute basicity to benzene if the above cutoff value were shifted to somewhere below 0.008 e/Ǻ². Such a change, however, would affect the whole structure of the COSMO-RS model and its HB calculations. In our case, it would affect not only the HB descriptors but also the polarizability molecular descriptors [37] (cf. Table SI1 of SI1 file).

A rather more severe drawback of COSMO-RS model is its inability to distinguish distant HB sites of a solute molecule from the neighboring ones. As mentioned already, COSMO-RS calculations are done through the sigma – profiles of the molecules [38,39,40,41,42,43]. But the three-dimensional specificities of the molecule as lost to a significant degree in its sigma – profile. This is a more general problem in COSMO calculations and has also an impact in polarity related calculations. All these COSMO-RS drawbacks are inherited in the present predictive method for $\Delta G_{12}^{hb}$. Luckily, the use of two sets of descriptors for each HB multi-sited solvent seems to alleviate these problems in organic solvents. The problems persist, however, in aqueous systems and require more information on the stereo-chemical structure of the solute molecules.

The results reported in SI2 file constitute, in essence, an extensive validation of the present predictive method. Thus, the development of the method was heavily based on LSER estimations of$\Delta G_{12}^{hb}$. An alternative source of information on HB interactions is Hansen’s solubility parameters (HSP) [50]. It would be interesting, then, to compare the above predictions with the HB information from HSPs. In order to translate the above predictions to HB HSP, the bridging concept of partial solvation parameters (PSP) [51,52] may be used. The hydrogen-bonding PSP of a self-associated compound i is given by

$${\sigma }_{hb}^2={\displaystyle \frac{-N_{hb}\Delta E_{ii}^{hb}}{V_m}}$$

(3)

where, V_m is the molar volume, $\Delta E_{ii}^{hb}$ is the HB contribution to self-solvation energy obtained in our previous work [1], and the number of hydrogen bonds, N_hb, is given by

$$N_{hb}=1+1/2K-\sqrt{1/K+1/4K^2}$$

(4)

with

$$-RT\ln K=\Delta G_{ii}^{hb}$$

(5)

At 25 °C, the HB parameter σ_hb is identical to Hansen’s HB solubility parameter, δ_hb. The HB contribution to self-solvation free energy may be obtained from Equation (5) or Equation (8).

The predictions, σ_hb,pred, of Equation (15) are compared with the corresponding experimental Hansen’s HB solubility parameter, δ_hb.[50] in Table 3 for a number of common self-associated solutes. As observed, for the bulk of solutes, the agreement is rather satisfactory. It is worth observing the discrepancies in the case of phenols, which also indicate that the aromatic ring is a source of basicity beyond what COSMO-RS model estimates. In such a case, an analogous to the above correction of the basicity descriptor, β_G, of aromatic solutes by a constant term should be applied also to the corresponding descriptor, β_E, for the interaction energy. However, as mentioned above, this type of corrections should account for their impact on the rest of intermolecular interactions. This is a major issue and will be dealt with in a forthcoming publication.

As mentioned in our previous work [1], the valuable information on the spatial or 3D distribution of HB sites in a molecule (cf. figure 1) is lost to a significant extent in the sigma - profile of the molecule. This information on the special distribution of molecular surface charges as well as the information on conformational changes of solute upon solvent changes are crucial for understanding the HB behavior of complex multi-sited molecules. Accounting for this information is not trivial but the needed tools are probably already available [7,8,9,10,53,54,55]. Hopefully, the present work will stir a broader interest in the scientific community with expertise on these tools.

5. Conclusions

A simple predictive method for the HB contribution to solvation free energy was presented and tested against corresponding LSER estimations for a variety of solute – solvent systems. The method is, in essence, a direct extension of our previous work on the prediction of HB interaction energies [1]. With their combination, an integral approach is being proposed for the prediction of the hydrogen – bonding contribution to solvation energy (enthalpy), free energy and, consequently to solvation entropy as well. In general, the agreement of predictions with LSER estimations is rather satisfactory. The predictive method, besides being simple, it is also free from the main limitation of LSER method which is heavily dependent on the availability of extensive experimental data. In this respect, the present method may significantly expand the range of applicability of LSER method. In particular, it may turn it into a rich source of valuable information on hydrogen bonding interactions directly transferable to Molecular Thermodynamic models, to other QSPR – type methods or to Hansen solubility parameters. In spite of its limitations, it is hoped that the present work will contribute in the establishment of a widely accepted reference scale for the strength of hydrogen bonding interactions. Such a broader collaboration could bring valuable feedback for further improvement of te predictive method.