Quantum Chemical (QC) Calculations and Linear Solvation Energy Relationships (LSER)-Hydrogen-Bonding Calculations with New QC-LSER Molecular Descriptors

1. Introduction

The new technological advances and their ever changing landscape brings up new challenging problems and opens new research horizons for efficient and insightful solutions in new material or new process designs. Current developments of modern Molecular Thermodynamics exploit the widely available high computational speed, even at the laptop level, and derive valuable information from quantum chemical calculations and molecular simulations. Quantum chemical suites at various levels are now widely and freely available. Atomistically-detailed force fields developed with information from quantum chemical calculations are extensively used allowing the computer-based reliable prediction of thermodynamic / thermophysical properties, phase equilibria or solvation quantities [1,2,3,4,5,6].

On the other hand, in numerous sectors, including the pharmaceutical and biomedical sector, the need often arises for an efficient fast solvent screening or reactant selection with specific properties. Quantitative structure – property or structure – activity relationships (QSPR / QSAR) are the preferred approach in most practical applications. One of the most successful QSPR-type approaches is Abraham’s LSER (Linear Solvation Energy Relationship) model [7,8,9,10,11] and has led, ultimately, to the valuable and freely available comprehensive and rich in thermodynamic information LSER Database [12]. Abraham’s legacy is continued by his collaborators and other developers worldwide. Bill Acree, having been one of Abraham’s closest collaborators, is conti nu ing the development and expansion of the LSER model with particular focus on solvation enthalpies [10,13,14]. The LSER model is currently one of the simplest, robust and widely used tools for thermodynamic calculations in a remarkably broad range of applications and often more successful than other more involved models [15].

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

$$Log{K}_G={\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 \text{(LSER)}$$

(1)

and

$$Log{K}_E=-{\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 \text{(LSER)}$$

(2)

Analogous equations apply for the solute transfer / partitioning between two condensed phases. 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 hydro gen bonding (HB) acidity A, and hydrogen-bonding basicity B, respectively [7,8,9,10,11].

Besides the traditional solute-partitioning / transfer studies and the limited Henry’s law, the solvation free-energy, ΔG₁₂, and its components enthalpy, ΔH₁₂ and entropy, ΔS₁₂, are directly connected with phase equilibrium studies involving activity or fugacity coefficients through the following most useful classical working equation:

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

(3)

V_m2 in Equation (3) 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). As seen, by knowing ΔG₁₂ and ΔH₁₂ from the LSER equations 1 and 2, one has a direct information on ${\gamma }_{1/2}^{\infty }$over a range of external conditions, which is crucial, among others, for phase equilibrium calculations [16,17,18]. This explains the interest of thermodynamicists, especially chemical engineers, in LSER-type models and their solvation equations.

There are, however, two main drawbacks of this otherwise very useful LSER model [7,8,9,10,11,12]. The LSER descriptors and the corresponding LFER coefficients are, typically, determined by multi linear regression of experimental data. The model expansion is, thus, restricted by the availability of experimental data. The second drawback is the thermodynamic inconsistency of the way the linear equations 1 – 2 are currently used. Although there is a sound thermodynamic basis of the LSER equations, their current use leads to peculiar results when applied to self-solvation of hydrogen-bonded solutes and very far from the expected equality of the comple mentary hydro gen-bonding interaction energies, since solute and solvent become identical on self-solvation. This problem was extensively discussed in the recent literature along with the significant advantages an eventual consistent reformulation would have [16,17,18].

New reliable experimental data on solvation or hydrogen-bonding (HB) quantities are becoming more and more scarce following the related scarcity of research groups on experimental thermodynamics worldwide. There is, in fact, an enormous amount of experimental information and calculations on these quantities already available [4,12,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34] but, unfortunately, the scatter of data often reach very high levels, of several thermal energy (RT) units, even for well-studied systems such as water or alkanols and their mixtures. It seems also that little care for thermodynamic consistency was applied to the proposed equations and scale correlations, which seriously restricts or rather forbids the reliable exchange of information between the various compilations and scales [16,17,18]. In this respect, the decisions for use of data exclusively from ab initio or high level DFT quantum chemical calculations for the development of predictive tools of solvation quantities [6] or the correlation of exclusively their own experimental solvation data for testing new descriptors [35] appear well plausible.

Strong specific interactions, such as hydrogen-bonding or Lewis acid / base interactions are ubiquitous and very often of primordial importance in chemical and biochemical processes and systems of practical or industrial interest [20,21,22,23,24,25,36,37]. Developments of modern Molecular Thermodynamics do focus on these specific interactions of directional character. As an example, the core of one of the most advanced equation-of-state models today, the family of SAFT (statistical associating fluid thermodynamics) models [38,39,40], consists of Wertheim’s perturbation approximation of these association interactions [41]. However, the above mentioned uncertainty regarding hydrogen bonding quantities pervades also SAFT models with often very large discrepancies in reported hydrogen bonding strengths [42]. An alternative equation-of-state model [43,44,45], which handles hydrogen-bonding with the more versatile Veytsman statistics [46], was also developed and applied to numerous hydrogen-bonded systems exhibiting intermolecular, intramolecular or cooperative association, forming oligomeric associates up to three – dimensional network structures [47]. However, this NRHB (non-randomness with hydrogen-bonding) model, the SAFT models and, essentially, all similar models of current Molecular Thermodynamics require external information on hydrogen bonding interactions as well as information on conformational changes upon solvation or mixing, typically, obtained from quantum chemical calculations.

Nearly thirty years ago, Klamt has published his pioneering work on Conductor-like Screening Model for Real Solvents (COSMO-RS) [48], which became one of the central directions in the development of modern Molecular Thermodynamics. Sandler, Lin and collaborators [49,50] have further contributed in familiarizing thermodynamicists, especially chemical engineers, with quantum – chemical / COSMO – calculations. A key feature of COSMO-RS is the adoption of the simple nearest-neighbor pairwise additive interaction approach for the calculation of potential and solvation energies [48,49,50,51,52,53,54]. The combination of this approach with detailed quantum-chemical information on the distribution of molecular charge densities has turned COSMO-RS into a powerful a-priori predictive and popular tool in Molecular Thermodynamics, especially, in solute solvation / partitioning calculations [55,56]. Due to the structure of the model, calculation of the separate hydrogen-bonding contribution to the solvation free energy cannot be made but the corresponding contribution to the solvation enthalpy can be made. Thus, COSMO-RS may be used as a predictive tool for this contribution and compared with the corresponding LSER contribution. Above all, however, COSMO-RS may be used for a well based development and expansion of a consistent LSER-type model, not only for solvent screening applications but also as a reliable source of information on intermolecular interactions and conformational changes needed in activity-coefficient or equation-of-state models [18].

Thus, the main task of the present work is to explore the potential of using COSMO-RS for the development of simple enough but thermodynamically consistent linear solvation energy relationships. These simple equations will be used to correlate experimental solvation data from recent compilations and extract their hydrogen bonding information. For this purpose, new molecular descriptors will be derived in the next section based on the molecular charge-density distributions or sigma-profiles of the COSMO-RS model [48,49,50,51,52,53,54]. The new linear relationships will be tested against experimental solvation data, including self-solvation data, and compared with LSER calculations [12,13,14] and, of course, COSMO-RS estimations [57]. The significance and perspectives of these calculations will be critically discussed.

2. The New Quantum-Chemistry Based LSER (QC-LSER) Molecular Descriptors

In Abraham’s LSER approach only the molecular descriptors V₁ and E₁ are well defined by the molecular structure and the refractivity of the solute [7,8,9,10,11,12]. The descriptors S₁, A₁, B₁, although attributed a more or less clear physical meaning, are determined, in essence, via global optimization of data correlations or error minimization in the prediction of various solute / solvent properties in a variety of solvents. Thus, our first task in this work is to derive new molecular descriptors obtained solely from the molecular structure of the solute and the COSMO-based molecular surface charge distributions or sigma-profiles [48,49,50,51,52,53,54,57], which could replace the LSER descriptors S, A, B. This is not a novel idea. Abraham, Klamt et al. [58] have already used the moments of the σ-profiles (cosmoments) in order to derive expressions for the LSER descriptors. These cosmoments are used extensively in the literature for the calculation of numerous quantities in a variety of applications [51]. More recently, Dohnal [35] has also used the σ-profiles and derived new molecular descriptors. However, in none of these applications the concern was the consi stent reformulation of the linear relationships that would permit the extraction of reliable information on intermolecular interactions. This is what will be attempted in this and the next sections.

From its very definition, a sigma-profile (σ-profile) is the distribution of probabilities P_i = A_i /A to find a surface charge density σ_i on the (COSMO solvation / screening) cavity of the studied molecule [51]. The sum of the partial surface areas A_i (not to be confused with the LSER descriptor A_i) gives the total surface area, A, of the molecular cavity, while the sum of the partial charges σ_iA_i gives the total charge of the molecule or ion. Figure 1 is an example of a sigma-profile.

Klamt has derived from the σ-profile A^X(σ) of a solute X, a number of molecular descriptors or cosmoments M_i^X of type i, defined by the integral ${\displaystyle \int A^X\left(\sigma \right)}{f}_i\left(\sigma \right)d\sigma$, for various functional forms f_i(σ) [51]. Dohnal’s hydrogen-bonding descriptors are derived by setting f_i(σ) = σ_i integrating in the hydrogen-bonding ranges of the charge density, σ .and using the “centroid” concept [35].

Of interest to us is the derivation of descriptors reflecting the strength of intermolecular forces or interactions. A proper handling of the σ-profiles could give us information on these forces of electrostatic character. In this respect, the product of the surface charge density σ_i by the charge σ_iA_i or the product σ_i²A_i (being charge by charge over unit area, with units C²m^-2) appears more appropriate for reflecting an electrostatic force. Such a distribution function turns Figure 1 into Figure 2. A key difference of the two distributions is in the central dome above zero charge density, which becomes a valley in Figure 2 touching the ground zero of electrostatic interactions.

According to Klamt’s suggestion [51], charge densities to the left of the blue line in Figure 2 correspond to acidic or donor sites able to participate in strong specific or hydrogen-bonding interactions (HBD sites). Similarly, charge densities to the right of the red line in Figure 2 correspond to basic or acceptor (HBA) sites. The molecular sites with charge densities between the two lines interact with much weaker forces, primarily, dispersion / induction forces. Of course, they may also interact with sites beyond the two cutoff lines when the latter are not involved in hydrogen bonding interactions. Thus, interactions due to polarizability do not exclusively involve sites from the two central sub-sections of Figure 2. Representative distributions for other molecules are reported in Figures S1 and S2 in the Supplementary Information (SI) file.

The two cutoff lines and the zero charge density divide the total area under the distribution curve of Figure 2 in four sub-areas, two to the left towards the acidic sites and two to the right towards the basic sites. The areas pertaining to hydrogen bonding sites are marked by the subscript h (A_h the HBD sites and B_h the HBA sites). The two central “polarizability” sub-areas are marked by the subscript p (A_p the sites to the left with screening charge sign negative (as acidic ones) and B_p the sites towards the basic sites of positive screening charge sign). The four new QC-LSER descriptors are just the areas of these four sub-sections in Figure 2 and they may be obtained from simple integration as follows (the factor 10⁶ accounts for the units of σ and A in σ-profiles):

$$A_h={10}^6{\displaystyle \underset{-\infty }{\overset{-1}{\int }}{A}_i{\sigma }_i^{2}d{\sigma }_i}$$

(4)

$$A_p={10}^6{\displaystyle \underset{-1}{\overset{0}{\int }}{A}_i{\sigma }_i^{2}d{\sigma }_i}$$

(5)

$$B_p={10}^6{\displaystyle \underset{0}{\overset{1}{\int }}{A}_i{\sigma }_i^{2}d{\sigma }_i}$$

(6)

$$B_h={10}^6{\displaystyle \underset{1}{\overset{\infty }{\int }}{A}_i{\sigma }_i^{2}d{\sigma }_i}$$

(7)

σ- Profiles are available, free of charge, for thousands of molecules in the open literature, as an example in ref. [50]. 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 [59,60,61]. In this work, the σ-profiles from COSMObase [57] at the DFT / TZPVD-Fine level of quantum chemical calculations are used and the corresponding QC- LSER descriptors for common solutes are reported in Table 1.

The descriptors A_h and B_h are considered reflecting the HBD and the HBA capacity, respectively, of the molecule. They are the analogues of the A and B descriptors of the Abraham LSER model [7,8,9,10,11,12]. The descriptors A_p and B_p do not have direct analogues in the Abraham’s LSER model. Abraham’s polarity / polarizability descriptor, S, is much different from descriptors A_p and B_p. Dipole or quadrupole moments and dipolar interactions cannot be extracted from σ-profiles since there is no information in them on the spatial (3D) charge distribution. Conformer distribution provides with valuable information but not direct polarity information. Yet, σ-profiles may provide with valu able information for the electrostatic interactions between sites of the A_h or B_h subareas with the B_p or A_p subareas of Figure 2. In fact, COSMO-RS model has very successfully used this information via the “misfit” energy contribution to the potential energy [48,49,50,51,52,53,54]. Following this idea, the sum A_hB_p + B_hA_p could be considered a substitute to the LSER S descriptor. This sum is considered to reflect the intermediate between hydrogen-bonding and dispersion-interaction contributions to solvation energy, that is, to play the role the polarity/ polarizability interactions are playing in the Abraham’s LSER model. It is in this sense that the descriptors A_p and B_p are used in this series of works not as measures of the actual polarizability of the molecule.

An important advantage of the new QC-LSER descriptors of Table 1 over the Abraham’s LSER descriptors is their capacity to account for conformational changes of the solute upon transfer to the solvent. As seen in Table 1, two conformers are reported for ethylene glycol, 2-methoxy etha nol and 1,3-propylene glycol and 3 conformers for 2-ethoxyethanol (many more conformers are reported in ref. [57]). The σ-profiles and sigma-surfaces of these conformers are shown in Figure 3 and in the Supplementary Information (SI) file. As shown in Figure 3, the conformer c0 of 1,3-propylene glycol has a significantly reduced HBD capacity compared to conformer c3 (acidity peaks on the left-hand side of σ-profiles). Conformer c0 is hydrogen-bonded intra-molecularly and, thus, the proton donor is not much available for intermolecular association.

In Table 1 are also reported the conformational energies of the less stable conformers, namely, their excess energy over the stable conformer c0. As shown also in the SI file, the c0 conformers favor intramolecular association while the less stable expose much more their hydrogen-bonding sites and favor intermolecular association. In the case of 2-ethoxyethanol, both conformers c1 and c9 favor intermolecular hydrogen-bonding. However, the conformational energy of conformer c9 is much higher than that of conformer c1 and this will diminish the gain from hydrogen-bonding. Most likely, then, conformer c1 will be a more favorable structure of 2-ethoxyethanol for intermolecular hydrogen-bonding We will come back to these conformers in a later section.

Before using these new QC-LSER descriptors for calculations in real systems, it is worth comparing them with corresponding descriptors in the literature. As mentioned already, Klamt has proposed various cosmoments of the σ-profiles as molecular descriptors [51]. In Figure 4 are compared the cosmoment sig2 [51,57] with the square root of the product of the two QC-LSER polarizability descriptors of alkanes. As shown, the ratio of the two is nearly constant as both seem to be nearly proportional to the molecular size.

The third - order hydrogen-bonding cosmoments [51] are often considered in the literature as adequate hydrogen-bonding descriptors [14,15]. In Figure 5, the new hydrogen-bonding descrι ptors B_h are compared with the corresponding B LSER descriptor [12] and the cosmoment HB_acc3 [57] for a number of common hydrogen-bonded hydroxyl (OH-) containing solutes. As shown, although the descriptors are not equal, their trends are rather similar. A corresponding graph with the HBD descriptors shows even more scatter and, thus, each set of descriptors reflects the hydrogen-bonding interaction capacity of the molecules in its own way. More information on the new QC-LSER descriptors is reported in the SI file.

3. Linear Solvation Energy Relationships with the New QC-LSER Descriptors

Hydrogen bonding and strong polar interactions belong to the so called “complementary” interactions in the sense that a donor site interacts with an acceptor site (not with a “similar” donor site), a positively charged site with a negatively charged one or a basic with an acidic one. When two different molecules 1 and 2 having, both, donor and acceptor sites interact, two types of intermo lecular hydrogen bonds form, namely, one HBD₁ X HBA₂ and one HBA₁ X HBD₂. The strength of these two bonds are considered here to be proportional to the complementary products A_h1B_h2 and B_h1A_h2, respectively. This is in conformity with the thermodynamic basis of linearity of all these LSER – type models [16,17,18]. The form, then, of thermodynamically consistent linear relationships for the enthalpy of solvation of solute 1 by solvent 2 cannot be arbitrary and, in the limit of self-solvation, the sum of the above two products should be equal to 2A_h1B_h1 = 2A_h2B_h2 since solute and solvent become identical on self-solvation (subscript 2 becomes subscript 1). This is guaranteed by the following simple expression for the hydrogen-bonding term of the linearity equation:

$$HB-term=c_{h12}\left(A_{h1}{B}_{h2}+B_{h1}{A}_{h2}\right)$$

(8)

where, c_h12 is a characteristic constant of the 1 – 2 pair.

In a similar manner, the contribution of the intermediate “polarity / polarizability” interactions is considered to be given by the expression:

$$Polar-term=c_{p12}\left(A_{h1}{B}_{p2}+B_{h1}{A}_{p2}\right)+c_{p12}\left(A_{h2}{B}_{p1}+B_{h2}{A}_{p1}\right)$$

(9)

where, c_p12 is a characteristic constant of the 1 – 2 pair. In conformity, however, with the main stream of relevant linear relationships in the literature, the polar contribution is not activated when the solute is non-polar. In this way, the solute must possess at least one hydrogen bonding site and only the first of the two terms on the right-hand side of equation (9) must be kept. We will further expand on these two alternative views in a forthcoming publication.

As mentioned above, one of the deficiencies of the sigma profiles is the fact that they do not quite reflect the shape or the three – dimensional character of the molecule. The contribution to solvation energy of the dispersion interactions are affected by this 3D character and, thus, in this series of works the polarizability descriptors are not used for the estimation of dispersion interactions except for a rather minor correction. Instead, the well-defined and easily obtained Abraham’s descriptors V and E for the McGowan volume and the excess refractivity, respectively, will be used for the dispersion - induction interactions with the following simple expression:

$$Dispersion-term=c_v{V}_1^{\alpha }+c_e{E}_1+c_2$$

(1)

where, c_v, c_e and α are (expected to be) universal constants and c₂ is a characteristic constant of solvent 2. The sum of the above three terms gives the solvation enthalpy, ΔH₁₂, of the 1 – 2 pair of molecules. Thus, the logarithm of the solvation enthalpy (energy) constant, K_E, is given by the following working equation:

$$Log{K}_E=c_2+c_v{V}_1^{\alpha }+c_e{E}_1+c_{p12}\left(A_{h1}{B}_{p2}+B_{h1}{A}_{p2}\right)+c_{h12}\left(A_{h1}{B}_{h2}+B_{h1}{A}_{h2}\right)$$

(11)

or a similar equation containing both terms on the right-hand-side of equation 9. It is this form of equation 9 which will facilitate the direct use of LSER descriptors V and E and, also, will facilitate the direct comparison of the new molecular descriptors with the corresponding Abraham’s LSER descriptors.

In addition to the above sum, the solvation free energy, ΔG₁₂, will require a solvation – entropy contribution, ΔS₁₂, so that at any temperature, T, the left-hand side of equation 3 will hold true. In order to facilitate comparison and handle solvation entropy data, the final working equation of the present linear solvation energy relationship is cast in direct analogy with Abraham’s LSER model in the following form:

$$Log{K}_G={\displaystyle \frac{-\Delta G_{12}}{2.303RT}}={\displaystyle \frac{-\Delta H_{12}}{2.303RT}}+{\displaystyle \frac{\Delta S_{12}}{2.303R}}=Log{K}_H-Log{K}_S$$

(12)

LogK_G may be cast also in a form similar to equation 11 but the coefficients cannot be arbitrary since always equation 12 should be valid. When adopting the LSER equations 1 – 2, thermodynamic consistency implies that the following equation should hold true:

$$Log{K}_S=\left(c_{e2}-c_{g2}\right)+\left(e_{e2}-e_{g2}\right)E_1+\left(s_{e2}-s_{g2}\right)S_1+\left(a_{e2}-a_{g2}\right)A_1+\left(b_{e2}-b_{g2}\right)B_1+\left(l_{e2}-l_{g2}\right)L_1 \text{(LSER)}$$

(1)

4. Applications

As mentioned already, one of the main tasks of this work is the extraction of information on hydrogen-bonding quantities from the corresponding solvation quantities. COSMO-RS [57] gives information on hydrogen-bonding enthalpy but no information on hydrogen-bonding free-energy or entropy. Thus, this section is divided in two parts. The focus in the first part is on self-solvation calculations in a variety of common solutes and the extraction of hydrogen-bonding quantities of self-associated solutes. The focus on the second part is on the extraction of consistent hydrogen-bonding energies in representative mixtures including aqueous systems. The calculations make extensive use of relevant experimental data on solvation quantities from the CompSol compilation [19] and COSMO-RS estimations [57].

4.1. Self-Solvation Calculations

The CompSol database [19] contains information on self-solvation free-energy, enthalpy and entropy for a number of common solutes and has systematically used the thermodynamic equation 3. This is a very useful information for model developments and verification of calculations. In order to directly apply the equations of the previous section, the experimental data are, first, transformed into their logarithmic equivalents (Log₁₀) as is the practice of Abraham’s LSER [7,8,9,10,11,12] model (cf. eqs. 1–2]. In order to recover solvation free-energies and enthalpies in kJ/mol units, one has to multiply LogK_G and LogK_E by the transformation factor -2.303RT/1000 = -5.71 kJ/mol. Similarly, the solvation entropy in JK^-1mol^-1 is recovered by multiplying LogK_S by the transformation factor -2.303R = -19.15 JK^-1mol^-1.

For the calculation of these self-solvation quantities, we must, first, determine the constants and coefficients in equations 10 and 11. On the basis of equation 12, if we will be able to calculate the logarithms for energy (enthalpy) and entropy, the logarithm of self-solvation free energy will be obtained by simple subtraction. On the other hand, since there are plenty of experimental data on heats of vaporization (the negative of self-solvation energy) it would be useful to focus also on one of the other two quantities (their logarithms) and, by subtraction, get also the third quantity.

In order to simplify things, we will start with the solvation quantities of normal alkanes. For alkanes only the dispersion – induction term will contribute to the solvation energy since there are no polarity or hydrogen-bonding interactions present. As mentioned above, the general dispersion – induction term is considered here to be given by equation 10. Since, by definition, the excess refractivity E is equal to zero for alkanes, their working equation will be:

$$Log{K}_E=c_v{V}_1^{\alpha }+c_1 \text{(n-alkanes)}$$

(14)

In fact, as shown in Figure 6, the experimental self-solvation enthalpies of n-alkanes are described well by the following simple equation:

$$Log{K}_E=5.55V_1^{1.07} \text{(n-alkanes)}$$

(15)

In a similar manner, the self-solvation entropy of n-alkanes is described well by the following simple equation:

$$Log{K}_S=2.5V_1^{1.07} \text{(n-alkanes)}$$

(16)

On the basis of the above equations 15-16, the self-solvation free-energy of n-alkanes is given by:

$$Log{K}_G=3.05V_1^{1.07} \text{(n-alkanes)}$$

(17)

Calculations with equations 15, 16 and 17 are compared with experimental data in Figure 6.

In view of the simplicity of the equation for LogK_S, it is interesting to examine the corresponding expression for non-hydrogen-bonded solutes. Obviously, some correction must be applied to equation 16, which would take into account at least the size difference from the corresponding “homomorph” alkane with carbon atoms equal to the number of non-hydrogen atoms of the solute. In Table 2 are reported the experimental LogK_s and the calculated ones for some representative solutes. The calculations were made with the following simple equation:

$$\begin{array}{l}Log{K}_S=2.5V_1^{1.07}+2\ln {\displaystyle \frac{V_{\mathrm{hom}}}{V_1}}\\ or\\ Log{K}_S=2.5V_1^{1.07}+4.6Log{\displaystyle \frac{V_{\mathrm{hom}}}{V_1}}\end{array}$$

(18)

V_hom is the McGowan volume of the homomorph alkane of the solute with volume V₁ [7,8,9,10,11,12].

This definition of the homomorph alkane is operational for solvents with more than two non-hydrogen atoms. Our focus here is on the self-solvation of liquid solvents and their homomorph counterparts should also be in a condensed / liquid state. Thus, instead of methane, the homomorph volume for water is taken equal to 1/4 of the McGowan volume of butane and, instead of ethane, the homomorph volume for methanol is taken equal to 2/4 of the volume of butane. It is worth mentioning that this choice of homomorph volume for water turns its size correction term in equation 18 negligible and, thus, the simple equation (16) applies also to water for the non-hydrogen-bonding contribution to LogK_S. It is also worth mentioning that the general validity of equation 18 may be used to define the homomorph in the case of halogenated compounds. Thus, for carbon tetrachloride, equation 18 reproduces its known LogK_s by setting its homomorph volume equal to 0.835 of the volume of n-heptane or V_hom = 0.835 V_n-C7.

Equation 18, though empirical, is an important result. The observed discrepancies in Figure 6 and Table 2 between experimental and calculated LogK_S are rather tolerable, in view of the scatter of the experimental data [19]. In fact, a deviation of the order of 0.2 in LogK_S implies a deviation of 0.2 X 19.15 = 3.8 JK^-1mol^-1 in solvation entropy and this, in turn, implies an error of 3.8 X 0.298 = 1.13 kJ/mol in solvation free energy or enthalpy, which is well within the experimental uncertainty of these quantities. Thus, equation 18 will be adopted, in the present work, for the estimation of the non-hydrogen-bonding contribution, LogK_S,nhb , to the solvation entropy. If the overall solvation entropy is known, the corresponding hydrogen-bonding contribution is obtained by subtraction, or:

$$Log{K}_{S,hb}=Log{K}_S-Log{K}_{S,nhb}$$

(19)

Equations 18 and 19 may, now, be coupled with information on hydrogen bonding enthalpy from COSMO-RS [57] in order to obtain the hydrogen-bonding free-energy of the studied solute through the classical equation:

$$\Delta G_{HB}=\Delta H_{HB}-T\Delta S_{HB}$$

(20)

This comprehensive hydrogen – bonding information for a number of common self-associated solutes is reported in Table 3. The corresponding LSER estimations of ΔG_HB are available for some of the solutes [62] and are reported in the last column of Table 3 for comparison .

The results reported in Table 3 are particularly interesting. The hydrogen-bonding enthalpies are, essentially, very close to the values more often used in current equation of state or activity coefficient models for the strength of hydrogen-bonding interactions [42,47]. Quite often in practice, these self-association strengths are used as a basis and the corresponding cross-association strengths are obtained by simple combining / mixing rules [42]. However, as we will see in the next sub-section, there is an alternative route for direct information on these cross-association strengths. On the other hand, the calculated hydrogen-bonding free-energies are also close to the estimations by Abraham’s LSER model [62], as shown by a comparison of the data in the last two columns of Table 3. The observed relative large discrepancy in the case of acetic acid is, most likely, due to the fact that lower carboxylic acids are dimerized in their pure state. It is not certain that all reported estimations have properly taken into account this dimerization for acetic acid. As shown in the table, the impact of this dimerization is significantly reduced after pentanoic acid and the entropy contribution assumes the expected larger absolute values of strongly hydrogen-bonded compounds. This is one example of the importance of accounting for self-solvation, especially for self-associated solutes, by the various linear relationships.

4.2. Hydrogen-Bonding Calculations in Solute (1) – Solvent (2) Systems

More often needed is the strength of cross-association interaction. If an alcohol is interacting with an ester or an ether, the required information for alcohol is its donor capacity. If it is interacting with chloroform, the required information is its acceptor capacity. And if it interacts with phenol or alkylamine, both, its donor capacity and its acceptor capacity must be known. In view of this, it is crucial to extract information on hydrogen-bonding interactions from solvation quantities in a consistent manner, caring for the separate donor and acceptor contributions.

Let us consider again a binary mixture of cross-associated solute (1) with solvent (2) and let their hydrogen-bonding donor and acceptor capacities be HBD₁ and HBA₁ for the solute and HBD₂ and HBA₂ for the solvent. The crucial information for Molecular Thermodynamic models is the strength of the two cross-association interactions HBD₁ X HBA₂ and HBA₁ X HBD₂ and of course the two self-association interactions HBD₁ X HBA₁ and HBD₂ X HBA₂. Obviously, when the solute becomes identical to the solvent the above cross-association strengths should become identical to the corresponding self – association strengths. This is precisely what is guaranteed with equation 8 or 11. Unfortunately, the existing LSER models and acidity / basicity scales in the literature have not paid much attention to this fact and their interaction strengths cannot be used for thermodynamic calculations. Luckily, the overall hydrogen-bonding enthalpies calculated by Abraham’s LSER model [7,8,9,10,11,12] are, very often, close to the corresponding COSMO-RS and equation-of-state calculations [18] and, thus, still are a useful information.

In Table 4 are reported the HB interaction strengths for representative solute – solvent systems. The LSER calculations [62,63] for the HBD₁ X HBA₂ and HBA₁ X HBD₂ interactions are reported in the second and third column, respectively. The fifth and sixth column of the Table reports the corresponding calculations from the COSMO-RS estimations [57]. COSMO-RS estimates the overall HB interaction HBD₁ X HBA₂ + HBA₁ X HBD₂. This overall strength is split in its separate contributions through equation 8. The equality of the two results specifies the binary coefficient c_h12 which is reported in the fourth column of Table 4. A number of comments regarding these calculations are in order.

As observed in Table 4, the overall LSER estimations of the cross-association interactions are often not very different from COSMO-RS estimations. It is their splitting in the separate interactions which matters for thermodynamic calculations since these values, extracted from infinitely dilute solutions, are typically used throughout the full composition range of the mixture. As seen, in all self-solvation cases the LSER calculations give two significantly different estimations for the same quantity HBD₁ X HBA₁. When the solvents are ethanol or 1-octanol and the solute is an alkanol, the hydrogen-bonding interactions are of the OH---OH type and one would expect some symmetry in the HBD₁ X HBA₂ and HBA₁ X HBD₂ interaction values. As seen in Table 4, the two LSER estimations are significantly different and this is rather difficult to rationalize.

In contrast to LSER calculations, the splitting of the COSMO-RS estimations gives two identical values on self-solvation and symmetric values for the OH---OH interactions. In fact, the reported values in the last two columns of Table 4 are quite close to the typically adopted values in equation-of-state calculations [42,47]. This is important in practice and in model developments since there are no limitations in COSMO-RS calculations for both the hydrogen-bonding contribution to the solvation enthalpy as well as to the QC-LSER descriptors even for non-synthesized compounds. If not already available, the required quantum chemical calculations can be done with available software suites [59,60,61].

It is worth observing in Table 4 that the binary coefficient c_h12 takes nearly a common value for each class of solute molecules or, rather, for each type of donor – acceptor pair. Thus, in the case of aqueous systems, this coefficient is nearly equal to 0.555 for all alkanols except for methanol. In the case of esters, it is nearly equal to 0.385 and in the case of alkanals it is equal to 0.35. This indicates that the new QC-LSER descriptors discriminate properly the type of interacting HB sites and the above coefficients may be used for the prediction of HB interaction strengths with other members of the homologous series. Of course, many more and diverse systems have to be studied in order to derive safe rules for these predictive calculations.

Of much interest is the availability of QC-LSER descriptors for conformers, when conformer-specific calculations are needed. Conformational changes may occur upon transfer of the solute from the gaseous to the liquid state. Of particular interest are conformers with intramolecular association. Missopolinou et al. [64,65] have studied the thermodynamics of alcoxyethanol mixtures with inert solvents and have shown that their skewed heat-of-mixing diagrams can be explained by accounting for the conformational changes of alcoxyethanols upon change of mixture composition. At infinite dilution of ethoxyethanol in solvent n-octane, the prevailing conformer is the c0 shown in Figure 7 and Figure S5 of SI file. Intramolecular hydrogen-bonding is possible with this conformation. The practically isolated ethoxyethanol molecules in the dilute concentration region prefer this conformation since the intramolecular association reduces the system potential energy and free energy. However, as the concentration increases, the probability for intermolecular association increases and it is the preferred association since it is stronger than the intramolecular one. The more open conformers c1 or c9 of Figure 7 favor intermolecular interaction much more than the conformer c0. There is of course a penalty of conformation energy in going from conformer c0 to these more extended conformers but this is offset by the gain from the stronger intermolecular interactions. These considerations explain the dispersive, intermolecular and intramolecular contributions to the experimental heat-of-mixing, especially the last one which is negative [64].

As mentioned already, information on hydrogen-bonding strengths and conformational changes are essential but extraneous to activity-coefficient or equation-of-state thermodynamic models. It is on account of this that the above COSMO-RS information, incorporated in the new QC-LSER descriptors, is particularly useful. As seen in Table 4, the solute – solvent HB interactions change significantly with conformation change and this must be accounted for by the molecular thermodynamic models, like the above mentioned SAFT or NRHB models.

COSMO-RS provides also with values for the total solvation energy [57]. In principle. then, one could combine it with the more easily available [19] total solvation free-energy and obtain the corresponding total solvation entropy. A scheme analogous to the one discussed above for self-solvation (cf. equations 16 and 18) could be devised for the mixtures. Again, the starting point could be n-alkane solvents. In this case, a correction should be made to equation 18 accounting at least for the size difference of solute and solvent. For all other solvents the additional correction should again account at least for the size difference of the solvent from its homomorph n-alkane. In this way, both, the hydrogen-bonding and the non-hydrogen-bonding contributions to the solvation entropy may be obtained. In a forthcoming publication we will see how this scheme may be applied, not only to systems of inert solvents like n-hexane or n-hexadecane but also to self-associated solvents such as ethanol or water and aqueous systems.

5. Discussion and Conclusions

As mentioned in the Introduction, hydrogen-bonding (HB) interactions are ubiquitous, playing a vital role in biological systems and processes but also in systems of practical or industrial interest. As such, they must be attributed a concrete value as any other properties and parameters in thermodynamic calculations. Unfortunately, in spite the enormous experimental and theoretical work done for nearly one century now, there is not even one single system whose strength of HB interaction is unanimously accepted as precisely known quantity. Various research schools define and measure differently in various environments the presumably very same interaction. The various acidity / basicity scales have their own essentially arbitrary scale – origin or zero interaction strength, which prevents the reliable exchange of information. Often, no due attention is paid to the thermodynamic consistency of the proposed correlations and predictive equations. One example is the inconsistency of the current linear relationships of solvation quantities when applied to self-solvation. This is really unfortunate because self-solvation is probably more appropriate for devising a rational widely acceptable scheme for the precise determination of HB quantities.

The present work is not solving this problem but it indicates a potential route towards this end. Quantum chemical calculations or quantum-mechanics based molecular simulations may be used to determine the strength of HB interaction in a pure self-associated fluid at well-defined external conditions and express it in understandable and replicable terms. Ethanol or propanol are examples of molecules that might be carefully studied. They are small enough with no extensive conformer populations and large enough to sterically hinder the simultaneous attack of both acceptor sites (lone pairs) of oxygen. Well-founded averaging schemes could be adopted for the determination of the OH--O interaction. Various attempts from various research schools would, sooner or later, coincide on the value of it within a tolerance less than the thermal energy (RT) unit. This is the information needed by statistical thermodynamic equation-of-state models empowered with well-founded hydrogen-bonding statistics.

The solvation quantities are clearly defined from the thermodynamics standpoint and are measurable precisely. Consistent self-solvation entropy values appear highly informative and the scheme leading to equations analogous to equations 16 and 18 deserve further exploration as a potential route for the sufficiently precise determination of the non-hydrogen-bonding contribution to the self-solvation entropy and, thus, to hydrogen-bonding contribution as well. A comprehensive table analogous to the above Table 3 could be produced by a well concerted international effort which could report unanimously accepted values for the hydrogen-bonding contributions to self-solvation entropy, ΔS_HB, enthalpy, ΔH_HB and free-energy, ΔG_HB, of well-studied solutes. Such a table could act as a reference or standard measure of hydrogen-bonding interactions.

In a self-associated pure solvent at equilibrium, at given external conditions, hydrogen bonds are continuously formed and disrupted but, on average, there is one degree of association and, thus, one average number of hydrogen bonds which share the overall hydrogen bonding contribution to the solvation or potential energy of this well-defined pure system. It is this share of hydrogen-bonding energy which may be determined with good and eventually widely accepted precision. Statistical thermodynamic HB modeling might, then, be used to decouple multi-site interactions in solutes like water, glycols, alkoxyalkanols, etc. The fact that ethylene glycol has two hydroxyls does not mean that its overall averaged HB interaction strength should be twice the strength of each separate OH--O interaction. The latter strength, when intermolecular, would be close to the one determined in self-associated propanol or butanol but it would be different when it is an intramolecular interaction. Statistical thermodynamics may account for this with the proper information from quantum chemical calculations. This is a first key step for turning the vast available relevant information into a thermodynamically usable quantity with very important implications.

The second crucial step is the use of thermodynamically consistent solvation energy relationships, such as the above equation 11. There is a key prerequisite for this step: the existence of consistent and coherent molecular descriptors, especially the descriptors of the HBD and HBA scales or of the acidity and basicity scales. Quantum mechanics seems to be very useful in this respect. Schemes analogous to the above equations 4 to 7 could be implemented as scripts in the existing relevant databases and produce consistent descriptors for thousands of compounds. At the same time, quantum chemical calculations of the strength of hydrogen-bonding energies in solute–solvent systems could be split into their HBD₁ X HBA₂ and HBA₁ X HBD₂ interactions as shown in Table 4. This may become a very useful information for testing hydrogen-bonding theories and for the development of advanced molecular thermodynamic models of mixtures. The above two steps are essential for the systematic thermodynamic study of hydrogen-bonding. Besides the descriptors used in the present work, analogous descriptors like the above mentioned Klamt cosmoments [57] or Dohnal descriptors [35] or any other analogous descriptors may be explored for analogous consistent HB calculations.

The HB enthalpies from COSMO-RS [57] or from LSER model [62,63] are close to the widely used values by activity–coefficient and equation-of-state models [42,47]. In contrast, the corresponding HB entropies differ significantly. As an example, the HB entropies used in the NRHB equation-of-state model [47] for alkanols, are in the range from -25 to -30 JK^-1mol^-1, which differs by more than 50% from the reported ones in Table 3. Care must be exercised before judging this difference. Entropy changes are sensitive to density changes and depend on molecular sizes / shapes. Density changes with external conditions, as well as molecular sizes / shapes, are taken explicitly into account by equation-of-state models. The above type of linear relationships do not account explicitly for the impact of density changes on HB interactions, although their cavitation terms account indirectly for solvent densities. Thus, the transfer of HB entropies in an equation-of-state model must be adapted to the specifics of the system and its treatment of hydrogen-bonding. As discussed before [16,17,18], partial solvation parameters or partial solubility parameters are a useful bridge for this transfer as they do take explicitly into account the molar volumes or densities of solutes and solvents. However, a systematic and conclusive study on this HB entropy transfer has not been done so far.

The focus in the present work was on hydrogen-bonding interactions. If these interactions were sufficiently precisely determined, the determination of the intermediate electrostatic or polarity / polarizability interactions could be very much facilitated. A forthcoming article is focusing on this and on calculations with equations 9 and 11.

Summing up, this work has shown how quantum-mechanics based calculations may be combined with experimental information on solvation quantities in order to extract information on hydrogen-bonding interactions. A simple scheme was used for the derivation of two pairs of molecular descriptors, one for the hydrogen-bonding interactions and one for the electrostatic interactions of the intermediate strength of polarity / polarizability interactions from the σ-profiles of COSMO calculations [57]. These descriptors were used to split in a consistent manner the COSMO-RS estimations of the HB interaction energies into their HBD₁ X HBA₂ and HBA₁ X HBD₂ component interactions. These COSMO-RS calculations were combined with experimental solvation entropies which permitted the extraction of information on hydrogen-bonding free-energy, enthalpy and entropy for a number of common solutes. The potential for extension and thermodynamic use of these calculations is significant for practical applications and molecular model developments.