Activity coefficients of binary methanol alcohol mixtures from cluster weighting.

The hydrogen bond network of different small alcohols is investigated via cluster analysis. Methanol/alcohol mixtures are studied with increasing chain length and branching of the molecule. Those changes can play an important role in different fields, including solvent and metal extraction. The extended tight binding method GFN2-xTB allows the evaluation and geometry optimization of thousands of clusters built via a genetic algorithm. Interaction energies and geometries are evaluated and discussed for the neat systems. Thermodynamic properties, such as vaporization enthalpies and activity coefficients, are calculated with the binary quantum cluster equilibrium (bQCE) approach using our in-house code peacemaker 2.8. Combined distribution functions of the distances against the angles of the hydrogen bonds are evaluated for neat and mixed clusters and weighted by the equilibrium populations achieved from bQCE calculations.

Introduction

Aliphatic alcohols are readily accessible low cost solvents, which are easy to extract and recover. Due to those properties their applications are considered and investigated in different fields such as alternative fuel production, solvometallurgy, hydrometallurgy, and solvent extraction. A potential extraction solvent must meet many criteria for a successful implementation, such as extraction performance, chemical stability, solvent regeneration, safety, and low environmental risk.1, 2, 3, 4, 5 Since they well match those criteria, for a long time, alcohols have been used to extract different metal ions. For example, Co(II) and Tl(I) chelat complexes are easily extracted by aliphatic alcohols in aqueous solvents.3 Indium forms a chelat with pyridylazonaphthol that can be extracted by butyl and pentyl alcohols.4 Gold(I) can be extracted from cyanide solutions such as Au(CN)2 by various alcohols.4 Vanadium and Niobium can be extracted from n‐octanol.6 Aliphatic alcohols are also known for their applications in the field of liquid/liquid solvent extractions.2, 5 There is evidence in literature, that the branching of an alcohol as well as the location of the hydroxyl group within the molecule does affect the extraction's capacity.5, 7, 8 If the branching of an alcohol is increased while the molecular weight stays constant, it was found that the separation factor increases as well.5

Offeman et al.5 proved that the position of the hydroxyl group, as well as the branching and chain length are important parameters that affect the ethanol extraction performance. From these considerations, it can be seen that the size and branching of alcohols affect their properties as hydrogen bond donors/acceptors. It is fundamental to understand how the hydrogen donor effect is related with the molecular configuration of alcohols.

Infrared and Raman spectroscopy were used to investigate the hydrogen network in alcohols and how it is affected by the alcohol's branching. Both experimental and quantum mechanical techniques were employed.9, 10, 11

Fourier transform microwave spectroscopy has been used in the past to investigate clusters of simple alcohols. In particular, many works investigated chiral dimers of methanol, ethanol, propanol, and butanol.12, 13, 14, 15, 16 These works found that the possible conformations of alcohol dimers involve significant dispersion interactions.

Classical molecular dynamics simulations17 and Monte Carlo simulations18 have been performed in the past, in order to investigate the hydrogen bond networks formed in different linear and branched alcohols. The results have shown systematic differences in their hydrogen‐bonded structures, depending both on hydroxyl group position and the molecular weight.

The hydrogen bond network of alcohols was investigated in the past with quantum mechanical approaches. Many works on this topic are present in literature, for instance, the conformation of 1‐butanol in the liquid phase was already studied in 1994 by Ohno et al.,19 wherein they demonstrated the importance of taking into account different conformations. The donor/acceptor configuration was investigated by Finneran et al.20 for the ethanol/methanol dimer. Rowley et al.21 analyzed the potential surface of many small alcohols, and Vargas et al.22 showed how the global minimum of the ethanol dimer is stabilized by the hydrogen bond.

An alternative approach to face the challenge is based on the binary quantum cluster equilibrium (bQCE) theory.23, 24, 25 bQCE is an extension of Weinhold's quantum cluster method for pure liquids26, 27, 28, 29, 30 and has been successfully applied to predict the miscibility of binary mixtures, the ionic product of water, activity coefficients, and mole fraction dependent dissociation for weak acids.25, 31, 32, 33, 34

By applying models of statistical thermodynamics to quantum chemically calculated clusters, the thermodynamic description of neat liquids and their mixtures at non‐zero temperature and pressure is possible in the condensed and gaseous phase. Self‐consistent‐field calculations lead to equilibrium populations of these clusters and thus an ensemble of different structural states is generated similar to molecular dynamics simulations.23, 24, 28 A first step in order to study hydrogen bond donor/acceptor systems was done by Brüssel et al.23 investigating the dimethyl sulfoxide/water system. Later, Matisz et al. were the first to study the binary methanol/water system.35 They found that cubic and spiro clusters are the dominant motifs in the mixed phase. Studies on methanol found that the liquid phase is formed mainly by cyclic ring structures.36, 37, 38 Liquid ethanol was found to be comprised mainly of the monomer, cyclic tetramer, and cyclic pentamer.39

From the quantum cluster approach we are able to evaluate the activity coefficients of binary mixtures.25 Those values are needed to determine phase equilibria,40, 41 and they are directly related to different phenomena, such as vapor pressure lowering and freezing point depression.40, 41 Activity coefficients are a convenient indicator for the deviation from ideal behavior42 and their theoretical determination is desirable, since in many cases they are not easily accessible experimentally. In particular, activity coefficients can be an important tool in the investigation and design of novel solvent mixtures. One example are deep eutectic solvents (DES),43 which since the beginning of last decade generated great interest44 and find a wide range of applications, such as metal extraction processes.45

In this article, we apply the quantum cluster approach to binary mixtures of methanol with different alcohols. In particular, we investigate the effect of molecule size and branching on the deviation from ideal behavior for small size alcohols (one to four carbon atoms). Our methodology can be found in the appendix, including details on the theoretical derivation of the binary quantum cluster equilibrium approach and the properties obtained from it, the computational details, and the generation of cluster sets.

Results and Discussion

Here, a large range of alcohols and their binary mixtures with methanol are investigated. The alcohols are chosen considering two factors: the number of carbon atoms (ranging from two to four) and the branching. Hence, both propanol isomers, n‐propanol and iso‐propanol, are considered as well as three isomers of butanol, namely n‐, iso‐, and tert‐butanol. Figure 1 shows a selection of clusters used in this work. Displayed are methanol clusters at different cluster sizes, dimers of all neat alcohols, and a set of mixed methanol/iso‐propanol clusters with different compositions. In total, 5760 cluster structures where quantum chemically optimized, and subsequently 1144 of them where selected by geometric and vibrational criteria (see the appendix for further explanations of the selection methodology). The binary mixture methanol/ethanol was already studied with the quantum cluster approach in an earlier work.25

Figure 1

Ball‐and‐stick models of some selected clusters. Top: methanol clusters at the size of 1–6 molecules. Center: Dimer geometries of the different alcohols. Bottom: Mixed methanol/iso‐propanol pentamers at different compositions. Please note that this selection shows only a small excerpt of all 1144 clusters.

Our results are mainly obtained via the extended tight binding method GFN2‐xTB46 (henceforth called xTB, see the appendix for details), which includes the D4 dispersion correction47, 48 accounting for the London dispersion energy and is an improved revision of the GFN‐xTB method,49 which we successfully employed for the calculation of activities and vaporization enthalpies in the past.25 xTB is a highly efficient method optimized for the calculation of geometries, vibrational frequencies, and noncovalent interactions, allowing the evaluation of thousands of cluster conformations which would not be feasible on DFT level. Additionally, xTB was found to perform well at computing the interaction energies of hydrogen‐bonded water clusters, outperforming even some GGA and hybrid DFT functionals such as BLYP and PBE0.46 Hence, we find this method is optimally suited for our approach.

Hydrogen bond analysis

Cluster analysis

Here, we will consider the interaction energy per monomer , where is the total adiabatic interaction energy in a cluster of the size of n molecules. Figure 2 shows the averaged interaction energies plotted against the cluster size n for the neat alcohols, as obtained from xTB. This average is taken from up to ten clusters per cluster size. The exact numbers of clusters per cluster size are included in the supporting information.

Figure 2

Interaction energy per monomer for the neat alcohols, averaged over multiple geometries at each cluster size. The left panel shows the effect of increasing chain length, the center and right panels display the effect of branching. Lines are meant to guide the eye.

In the case of the linear systems, depicted in the left panel, an increase in the cluster size leads to stronger (i. e. lower) interaction energies per monomer. In the middle panel n‐propanol and iso‐propanol are compared, but no particular differences are present. The right panel shows the different isomers of butanol investigated in this article. Whereas less stable in the case of the dimer, increasing the cluster size, the unbranched alcohol has slightly lower interaction energies compared to its more branched isomers. Table 1 shows the interaction energies of the global minimum methanol and n‐butanol clusters in more detail, along with the dispersion energy and remaining energy where . The set of global minimum structures of methanol is displayed in Figure 1. Again, we can observe that increases with the cluster size. With each additional molecule there's a gain in interaction strength that can be attributed to cooperativity effects, calculated as where is the average interaction energy per monomer in a cluster of size n. This cooperative gain decreases rapidly and seems to be mostly saturated at a cluster size of four molecules. In comparison, whereas is of similar size in methanol and n‐butanol, dispersion forces are considerably stronger in the latter. While smaller in magnitude, levels out less rapidly than and benefits from cooperative effects even in larger clusters.

Table 1

Average interaction energies , remaining energies , and dispersion energies per molecule in kJ/mol for the global minimum of methanol and n‐butanol clusters of size n, as obtained on the xTB level of theory, as well as the relative cooperative gain in %.

	Methanol				n‐Butanol
n	Δ_intE‾	Δ_dispE‾	Δ_remE‾	coop.	Δ_intE‾	Δ_dispE‾	Δ_remE‾	coop.
2	−10.2	−1.7	−8.8		−12.2	−5.4	−6.8
3	−20.0	−3.4	−16.9	97.5	−30.6	−8.6	−22.0	150.7
4	−26.8	−4.2	−22.6	33.4	−38.8	−12.8	−26.0	26.7
5	−28.5	−5.5	−23.0	6.1	−40.6	−13.9	−26.8	4.8
6	−29.5	−5.9	−23.6	3.7	−41.0	−15.8	−25.2	1.0
7	−30.2	−6.2	−23.9	2.3	−42.8	−15.4	−27.4	4.3
8	−31.4	−6.4	−25.0	4.0	−41.8	−16.4	−25.4	−2.4
9	−31.3	−7.0	−24.3	−0.3	−42.3	−18.2	−24.1	1.3

Table 2 lists the interaction energies of the global minimum structures of (ROH)2 dimers, in order to compare the performance of xTB against GGA methods. The complete list of interaction energies of all the 1144 clusters are given in the supporting information. The interaction energies obtained via the xTB method are compared to those obtained from single point calculations on the same geometries employing DFT methods, namely the GGA functional BP86 and the hybrid functional B3LYP. Overall, xTB interaction energies are weaker for unbranched alcohols and stronger for branched ones with respect to the DFT methods energies. Nevertheless, the trends are similar. Regardless of the method, the lowest and highest interaction energies are found for i‐BuOH and MeOH, respectively.

Table 2

Interaction energies (kJ/mol) of alcohol dimer geometries optimized at xTB level, as obtained from single point calculations at different levels of theory. The energies correspond to the global minimum geometries of the pure ROH‐ROH and mixed MeOH‐ROH dimers.

	Δint E(ROH‐ROH)			Δint E(MeOH‐ROH)
ROH	xTB	BP86	B3LYP	xTB	BP86	B3LYP
MeOH	−20.4	−22.7	−23.8	–	–	–
EtOH	−23.8	−26.1	−27.4	−23.9	−25.8	−27.1
n‐PrOH	−24.0	−26.6	−27.7	−23.9	−26.1	−27.4
i‐PrOH	−29.4	−27.3	−28.4	−27.5	−27.4	−28.4
n‐BuOH	−24.4	−27.0	−28.1	−24.2	−26.5	−27.7
i‐BuOH	−33.8	−28.1	−28.8	−24.7	−26.0	−27.3
t‐BuOH	−30.4	−27.5	−28.7	−27.8	−27.8	−28.8

–

–

–

The differences in interaction energies become less pronounced for the mixed dimers formed by methanol and an additional alcohol, listed in Table 2. Increasing the branching of the molecule, the differences between the methods become smaller. The intermolecular hydrogen‐oxygen distances and the complementary O−H⋯O angle within the dimers are listed in Table 3, for both xTB and BP86 optimized geometries. Overall, the distances are in good agreement, whereas the angles are slightly different. In the following, we will focus on our xTB results, exclusively, based on the 1144 calculated clusters optimized on that level. We observe that xTB can reproduce energetic and geometric features with sufficient accuracy and find that the ability to quantitatively evaluate a wide range of potential cluster geometries justifies the use of this method. Results obtained via the GGA method BP86 are available in the supporting information.

Table 3

Hydrogen bond distances r (pm) and the complementary O−H⋯O angles α (°) of geometry optimized dimer structures at xTB and BP86 level of theory. The distances and angles correspond to the global minimum geometries of the pure ROH‐ROH and mixed MeOH‐ROH dimers.

	Pure dimers ROH‐ROH				Mixed dimers MeOH‐ROH
	r(xTB)	r(BP86)	α(xTB)	α(BP86)	r(xTB)	r(BP86)	α(xTB)	α(BP86)
MeOH	189.6	190.1	9.3	7.6	–	–	–	–
EtOH	188.7	189.2	2.4	10.7	188.4	189.2	2.3	10.2
n‐PrOH	188.3	188.8	2.4	10.4	182.2	189.1	1.8	10.0
i‐PrOH	188.0	190.1	6.1	13.2	186.2	188.8	7.8	12.3
n‐BuOH	188.7	189.0	3.0	10.7	188.4	189.2	1.5	9.7
i‐BuOH	190.0	189.6	1.9	6.7	188.0	189.1	0.3	10.3
t‐BuOH	187.9	190.1	4.5	14.2	185.5	188.0	6.5	12.4

–

–

–

–

Increasing either chain length or branching of the molecule, the interaction energy is decreasing, both for the mixed methanol/alcohols dimers and the pure systems. No similar trend is observed for the hydrogen bond distance and angle. Nevertheless, the methanol dimer shows both a larger distance and wider angle as compared to the other dimers. The iso‐propanol containing dimers show a wider hydrogen bond angle compared to the other alcohols except methanol. The tert‐butanol containing dimers show the widest angle of the investigated butanol isomers. The hydrogen bond distance of 188.8 pm in the neat ethanol dimer is lower than the literature value of 191.0 pm found by Vargas et al.22 This difference can be imputed to the different level of theory. Nevertheless, both results are in acceptable agreement. The distance of 189.6 pm in the neat methanol dimer is in good agreement with the MP2 value of 187.2 pm reported by Provencal et al.50

Population‐weighted analysis

In earlier works, we presented sophisticated methods for detecting and quantifying hydrogen bonds.51, 52 From a geometrical perspective, hydrogen bonds are often characterized by their length and angle.50, 53 Different bond lengths and angles can bring to light distinct behaviors of the investigated species forming the hydrogen bonds. For this reason, we show combined distribution functions (CDF) of the different alcohols, constructed from the intermolecular hydrogen‐oxygen distances and the angular distribution of the complementary O−H⋯O angle. Since our cluster sets not only include global minimum structures but also those more distant from the enthalpically optimal binding situation, by combining the collected data of all investigated clusters, in total 1144, and weighting them by their bQCE populations (see appendix for method), we obtain CDFs similar in appearance to those of a MD simulation. Through the weighting by population, these CDFs are accessible for any temperature and pressure investigated in the bQCE calculation. Here, we investigate methanol, ethanol, n, and tert‐butanol in order to include both linear and branched alcohols. For all of them, the complete cluster set is analyzed with our in‐house trajectory analysis code travis;54 then, the data of each cluster is collected and weighted by the cluster population at 298.15 K, as obtained by bQCE calculations.

In Figure 3 the CDFs of the neat systems are reported. The color scale is relative and referenced to the maximum value of all systems. The average of the hydrogen bond distance is in the range of 170–180 pm, which is in good agreement with the literature values of methanol of both ∼180 pm53 (obtained by a combined experimental and molecular dynamics investigation). Several ab initio molecular dynamics studies employing the BLYP functional find the first peak of the radial distribution function of the O−H⋯O distance around 190 pm,55, 56 in good agreement with our results of 189.6 pm for methanol dimer (Table 3). Our CDFs show a range lower than this value, but are still in good agreement.

Figure 3

Combined distribution fuctions of the intermolecular hydrogen‐oxygen distance against the complementary of the O−H⋅⋅⋅O angle in (left to right) methanol, ethanol, n‐butanol, and (bottom) tert‐butanol.

Comparing the different alcohols in Figure 3, it can be seen that the CDFs become more localized with increasing size of the alcohol, from methanol to n‐butanol and tert‐butanol, with an increasingly distinct maximum observable in the distribution. Comparing the other alcohols to tert‐butanol, a shift of this maximum from 1–2° to a wider angle of 9–10° can be observed, which can be attributed to the different steric restrictions in the branched alcohol.

Figure 4 shows the CDFs of n‐butanol at a temperature of 248.15 K, 298.15 K, and 348.15 K, comprising a large span within the liquid range of n‐butanol. While the overall structure of the distributions remains the same, as expected, with increasing temperature, the hydrogen bonds become less localized in the lower range of angles and an increased distribution over the range of 8–12° is observed. Thus, with rising temperature, the average hydrogen bond angle increases.

Figure 4

Combined distribution functions of the intermolecular hydrogen‐ oxygen distance against the complementary of the O−H⋅⋅⋅O angle of n‐butanol, at the temperature (from left to right) 248.15 K, 298.15 K, 348.15 K.

Figure 5 shows the CDFs of methanol/ethanol (a), methanol/n‐butanol (b), and methanol/tert‐butanol (c) for increasing mole fractions of methanol of 0.2, 0.5, and 0.8. In general, for small mole fractions of methanol the systems show a more localized maximum of the hydrogen bond distribution. For the methanol/tert‐butanol system, it is visible that the preferred angle is shifting from larger to smaller values with an increasing mole fraction of methanol. Even more, compared to the other binary systems, this mixture is less localized due to the branching of tert‐butanol.

Figure 5

Combined distribution fuctions of the intermolecular hydrogen‐oxygen distance against the complementary of the O−H⋅⋅⋅O angle for the mixtures of methanol with ethanol (a), n‐butanol (b), and tert‐butanol (c) at mole fractions of methanol of (from left to right) 0.2, 0.5, and 0.8.

Thermodynamic properties of neat systems

In earlier works, we established our approach of calculating from QCE, by performing so‐called QCE0 calculations wherein a mf is set to 0, as reference for the gas phase.25, 57 For QCE0 calculations, the cluster sets are reduced to clusters up to the size of three molecules. can then be obtained for any temperature simply as difference of the total enthalpies in the liquid and gaseous phase respectively. In Table 4 our calculated enthalpies of vaporization at 298.15 K are listed next to their experimental reference value for every neat system investigated in this work. A good agreement with experimental data can be seen for most systems. For ethanol, an improved result with respect to the previous work25 (44.09 kJ mol−1) is obtained, due to the increased size and number of clusters. Our approach appears more accurate for the non‐branched systems methanol, ethanol, n‐propanol and n‐butanol. The system that deviates most from experiment is also the most branched, namely tert‐butanol. In general a larger deviation from experimental results can be seen with an increased branching of the molecule.

Table 4

Calculated and experimental enthalpies of vaporization and in kJ mol−1 for the neat systems at standard conditions. Experimental enthalpies of vaporization are taken from the NIST Chemistry WebBook.58

	Δ_vapH	Δ_vapH^exp
methanol	39.33	37.60
ethanol	43.36	42.30
n‐propanol	47.17	47.00
iso‐propanol	47.49	45.00
n‐butanol	51.37	52.00
iso‐butanol	49.81	51.00
tert‐butanol	46.71	39.70

Thermodynamic properties of binary mixtures

As shown in earlier works, via the quantum cluster approach we are able to reproduce quantitatively the experimental Gibbs energies of mixing , using the density and phase transition temperature as only experimental input data.24, 25, 59 The Gibbs energy of mixing at 328.15 K is depicted in Figure 6 for the binary mixtures of methanol with ethanol, n‐propanol, iso‐propanol, n‐butanol, iso‐propanol, tert‐butanol. Of all the investigated systems, methanol/ethanol is the one that most resembles an ideal mixture. In contrast, methanol/n‐butanol is the system deviating the strongest from ideality. In general, an increase in the deviation from ideality can be seen with an increasing size of the molecule, from ethanol to butanol. In order to investigate the deviation from the ideal mixture in more depth, activity coefficients are calculated. Activity coefficients are directly connected to the excess Gibbs energy of mixing as shown in Equation 9. In Table 5, activity coefficients of all mixtures are shown over the complete mixing range. As described before and in our previous work,25 the mixture methanol/ethanol is the most ideal system; the activity coefficients and of both methanol in ethanol and the opposite respectively are near to one for every mole fraction. An increase in the activity coefficient at infinite dilution can be observed with increasing size of the alcohol; methanol at infinite dilution in n‐propanol and n‐butanol shows activity coefficients of 4.39 and 148.21 respectively as compared to 1.03 in ethanol. Increasing the branching of the molecule, the activity coefficients are decreasing to values closer to one; this is in good agreement with the Gibbs energies of mixing in Figure 6, where the most branched system, methanol/tert‐butanol, is the closest to ideality second to only methanol/ethanol. The activity coefficients of methanol in iso‐propanol compared to methanol in n‐propanol confirm this behavior. For a better visualization, the activity coefficients are shown in Figure 7 and Figure 8 with respect to the methanol mole fraction. From these graphs it can be seen that for low mole fractions of methanol its activity coefficient in alcohols is increasing with the size of the solvent and decreasing with its branching. The same behavior can be found for the activities of alcohols diluted in methanol. We can conclude that in mixtures of methanol with an alcohol, increasing the size of the alcohol leads to a larger deviation from ideality. On the other hand, an increase in the branching of the alcohol leads to a more ideal mixture. The same behavior can be found in the experimentally calculated excess Gibbs energy of mixing from Polak et al.59 for isomeric butanol, where n‐butanol shows the largest deviation from ideality, and tert‐butanol the lowest.

Figure 6

Calculated excess Gibbs energies of mixing ▵mix G e for binary mixtures of methanol with an alcohol ROH at 298.15 K. x indicates the mole fraction of methanol.

Table 5

Activity coefficients f ROH of both components in mixtures of methanol (MeOH) with an alcohol ROH, where xm is the molar fraction of methanol.

	ethanol		n‐propanol		iso‐propanol		n‐butanol		iso‐butanol		tert‐butanol
x m	f MeOH	f EtOH	f MeOH	f nPrOH	f MeOH	f iPrOH	f MeOH	f nBuOH	f MeOH	f iBuOH	f MeOH	f tBuOH
0.00	1.03	1.00	4.39	1.00	2.03	1.00	148.21	1.00	21.09	1.00	2.49	1.00
0.10	0.99	1.00	2.58	1.03	1.51	1.01	20.01	1.10	6.99	1.06	1.87	1.01
0.20	0.99	1.00	1.83	1.09	1.31	1.04	5.82	1.36	3.40	1.20	1.53	1.05
0.30	0.99	1.00	1.45	1.18	1.21	1.07	2.77	1.74	2.12	1.40	1.33	1.10
0.40	0.99	1.00	1.24	1.28	1.14	1.10	1.78	2.20	1.56	1.64	1.21	1.16
0.50	1.00	1.00	1.13	1.39	1.09	1.14	1.37	2.72	1.29	1.92	1.13	1.23
0.60	1.00	0.99	1.06	1.49	1.05	1.19	1.18	3.27	1.15	2.21	1.07	1.30
0.70	1.01	0.98	1.03	1.57	1.03	1.24	1.08	3.85	1.07	2.51	1.04	1.38
0.80	1.01	0.97	1.02	1.64	1.02	1.30	1.03	4.44	1.03	2.81	1.02	1.47
0.90	1.01	1.00	1.01	1.73	1.01	1.37	1.00	4.99	1.01	3.19	1.01	1.58
1.00	1.00	1.16	1.00	1.98	1.00	1.55	1.00	5.43	1.00	3.82	1.00	1.76

Figure 7

Activity coefficients of methanol in binary mixtures with an alcohol ROH at 298.15 K. x indicates the mole fraction of methanol.

Figure 8

Activity coefficients f ROH of alcohols ROH in a binary mixture with methanol. x indicates the mole fraction of methanol.

Conclusions

In order to understand the liquid behavior we optimized 5760 cluster geometries of which 1144 were further analyzed. The average interaction energies per monomer have been evaluated for all neat dimers and stronger interactions with increasing chain length and branching of the molecule are found. Combined distribution functions of distances and angular distributions of hydrogen bonds are calculated for several neat and mixed systems, demonstrating how the different size and branching of the alcohol lead to different geometrical conditions of the hydrogen bonds. In this article we applied the quantum cluster approach to calculate the activity coefficient of different methanol/alcohol mixtures. This approach relies on the binary quantum cluster equilibrium theory.23, 24, 25, 60 With respect to our previous work,25 we increased the maximal size of the clusters from six to nine molecules, ensuring a better accuracy in the quantum cluster approach simulations. Furthermore, we studied a much wider range of MeOH‐alcohol mixtures. The vaporization enthalpies of all pure substances were calculated at room temperature; while overall a good agreement with experimental data can be observed, the deviation increases with the alcohol's branching, with tert‐butanol showing the largest deviation. Using Redlich−Kister polynomials and calculating their derivatives allows access to the activity coefficients, further establishing the bQCE approach as a novel and conceptually outstanding method of computing such quantities. In this article we demonstrated that in mixtures of methanol with an alcohol, increasing the size of the alcohol leads to a larger deviation from ideality. On the other hand, an increase in the branching of the alcohol leads to a more ideal mixture. This case study will help to move our approach to complex solvent media, adding to the tools used in application driven solvent design.

Computational Details

The bQCE method

The bQCE method has been established and its underlying theory detailed in depth in many earlier works.23, 24, 25, 27 Here, we will present only a short overview of the key equations of bQCE.

First, we consider a system of non‐interacting clusters in thermodynamic equilibrium, built up from one (neat substances) or two (binary systems) monomers. The equilibrium reaction between clusters of a binary system can be written as

where and are the number of monomers of each component and that form the cluster . The system's total partition function at volume V and temperature T is given by

where is the partition function corresponding to the single cluster and is the full set of total cluster populations . From calculating all the thermodynamic properties of the system are accessible. Each cluster partition function can be evaluated as product of partition functions corresponding to the cluster's different degrees of freedom:

where , , and are the translational, rotational, and vibrational partition function. They can be calculated from standard equations for the particle in a box, rigid rotator, and harmonic oscillator respectively.27, 62 The electronic partition function is calculated from the adiabatic binding energy of the cluster.63

In order to describe , the phase volume V must account for an exclusion volume V ex which attributes a volume to the non‐punctiform clusters. The exclusion volume is calculated as

wherein is the cluster volume. Since cluster volumes are sensitive to the choice of atomic radii, a scaling parameter b xv must be introduced. Additionally, is extended by a correction term to account for the interactions between clusters in form of a volume and cluster size dependent mean‐field energy. The electronic partition function then reads

where k B is the Boltzmann constant and the mean‐field parameter a mf (Jm3/mol2) is a second empirical parameter, that is scaling the strength of inter cluster interactions. In an optimized bQCE calculation, the parameters are chosen such that the deviation of the bQCE results from a given experimental input, such as densities and phase transition temperatures, becomes minimal.

To calculate , all independent variables need to be known. The temperature must be set by the user and the volume is restricted with respect to an externally applied pressure

If several combinations of V and exist that fulfill this condition, then the solution with the lowest Gibbs energy

is chosen. Using this approach, good performance has been demonstrated for mixed systems in several studies.25, 57, 64

Activity coefficients from bQCE

An accurate and detailed description of the calculation of activity coefficients via the quantum cluster approach can be found in a previous work.25

For a binary mixture, the excess Gibbs energy of mixing can then be calculated as

where and are the Gibbs energy and the ideal Gibbs energy of mixing respectively. Activity coefficients f are directly related to the excess Gibbs energy of mixing (labeled as ) by

where R is the ideal gas constant, T is the temperature, and N is the particle number of component i. Since no analytical expression for is available, we calculate its derivative numerically through a Redlich−Kister (RK) style polynomial fit which smooths out all local inconsistencies.65 We denote the Gibbs energy with

wherein g are the Redlich–Kister parameters.65 In this work we used up to five parameters . With we can write as function of N and N:

Equation 11 gives an analytical expression for that can be differentiated with respect to the number of particles N:

Analogously, we can evaluate

wherein and . Inserting these expressions into Equation 9 allows the evaluation of the activity coefficients.

Computational details and cluster search

The selection of the cluster set is a crucial step of the quantum cluster equilibrium approach.23, 28 In this work, we increased the cluster size (up to nine molecules) compared to our previous work,25 in order to increase the accuracy of the calculations. In Figure 9 the cluster generation procedure is explained. In the first step, the global minimum structure of each cluster composition is searched for by running a genetic structure optimization algorithm at the classical force field level of theory. For this purpose the ogolem framework,66, 67 interfaced with the amber 2016 molecular dynamics package68 and the general amber force filed (GAFF),69 is used. During the optimization, the number of individuals in each generation as well as the total number of individuals are set accordingly with the cluster size. For each cluster a number between 2000 and 6000 structures are evaluated, with each generation consisting of 100 to 300 individual structures. For each respective cluster composition, fifteen clusters are chosen from the final generation in order to sample a set of energetically and geometrically diverse individuals, which represent the global and local minima of that structure. These clusters are geometrical optimized at semi‐empirical level of theory, using the extended tight binding method GFN2‐xTB 6.0.1,46, 49 which includes the D4 dispersion correction47, 48 accounting for the London dispersion energy. Frequency calculations are performed with the same method. All the clusters with a first normal mode below a threshold (in this case, 10 cm–1) are removed from the cluster pool in order to avoid imaginary or low frequencies that could affect the simulations. Likewise, structural duplicates of already existing clusters are removed from the cluster set. The conformational similarity of two clusters is quantified by their geometrical distance d:70

Figure 9

Cluster generation procedure.

wherein I and are principal moments of inertia of the clusters and respectively. Clusters with a geometrical distance of were removed from the cluster set. For the geometrical investigation of hydrogen bonds, the cluster sets were analyzed with our in‐house trajectory analysis code travis.54 Note, that the all cluster geometries can be obtained from the authors upon request.

Additionally, a select set of systems were optimized on the DFT level of theory. Here, we used the orca 4.0.071 quantum chemical code employing the GGA functional BP86 with the 6‐31G* basis set, D3 dispersion correction, and geometrical counter‐poise correction.72 Due to the increased computational cost, we reduced the cluster set size to a maximum of six molecules. We observed no improvement in the results, hence, we excluded the analysis of those clusters in this article. However, interaction energies and thermodynamic results of those systems are presented in the supporting information. Additionally, single point simulations where also performed employing the hybrid functional B3LYP, using the same 6‐31G* basis set.

bQCE calculations were performed with the peacemaker 2.8 program package61 which has successfully been used to describe binary mixtures previously.24, 25, 35, 59 All calculations were performed at a fixed pressure of 1.01325 bar and temperature ranging from 273 to 500 K. Cluster volumes were calculated employing van der Waals volumes. The parameters a mf and b xv were fitted so that the deviation of bQCE results to experimental boiling points and densities taken from literature60, 73, 74, 75, 76, 77 become minimal. Different to our earlier works, here, we employed the Nelder–Mead optimization algorithm.77

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.