Model Performances Evaluated for Infinite Dilution Activity Coefficients

Subscriber access provided by ALBRIGHT COLLEGE

General Research

Model Performances Evaluated for Infinite Dilution Activity Coefficients Prediction at 298.15K Thomas Brouwer, and Boelo Schuur Ind. Eng. Chem. Res., Just Accepted Manuscript • Publication Date (Web): 29 Apr 2019 Downloaded from http://pubs.acs.org on April 29, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 2 3

Industrial & Engineering Chemistry Research

Model Performances Evaluated for Infinite Dilution Activity Coefficients Prediction at 298.15K

4

Thomas Brouwer, Boelo Schuur*

5 6 7 8 9

Sustainable Process Technology group, Process and Catalysis Engineering cluster, Faculty of Science and Technology, University of Twente, PO Box 217, 7500AE, Enschede, The Netherlands.

10 11

*Corresponding author: e: [email protected], p: +31 6 14178235

12 13

Abstract

14

The infinite dilution activity coefficient, 𝛾∞ 𝑖 , is often applied to characterize solvent-solute

15

interactions, and when accurately predicted, it can also serve as an early stage solvent selection

16

tool. Ample data is available on the use of a variety of models, which complicates decision making

17

on which model to apply when. A comparative study was performed for eight predictive models

18

at 298.15K, including the Hildebrand parameter and the Hansen Solubility Parameters. Also, three

19

Group Contribution Methods based on UNIFAC, COSMO-RS, the Abraham model and the

20

MOSCED model were evaluated. The MOSCED model and the Abraham model are overall most

21

accurate for respectively molecular solvents and ionic liquids, with Average Relative Deviations

22

of 16.2±1.35% and 65.1±4.50%. Therefore, cautious decision making based on predicted 𝛾∞ 𝑖 in

23

ionic liquids should always be done, due to the expected significant deviations. A MOSCED model

24

for ionic liquids could be a potential approach for higher accuracy.

25 26

Keywords:

27

Activity coefficients, infinite dilution, thermodynamic modelling, model comparison

28

1 ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 24

29

1. Introduction

30

The global community relies on the chemical industry for the production of goods from complex

31

raw materials, such as oil and biomass. The separation processes required in these production

32

routes account for up to 50 % of the total energy costs in refineries,[1] and improving the efficiency

33

of separations can significantly reduce the environmental impact of the chemical industry.[2] This

34

can only be achieved when the separation processes are understood on the molecular level,

35

which includes a good description of thermodynamic equilibria. An accurate description of these

36

equilibria is possible with models like UNIQUAC and NRTL, but requires labour intensive

37

experimental data. For initial stages of process design, including solvent selection and/or design,

38

less labour intensive approaches for understanding intermolecular interactions are desired. A

39

range of predictive models is available to provide engineers with first estimates for

40

(inter)molecular behaviour.[3-11]

41

An appropriate indicator of intermolecular interactions underlying thermodynamic

42

equilibria is the infinite dilution activity coefficient. [12] Activity coefficients describe the

43

thermodynamic non-ideality between two substances due to intermolecular interactions. These

44

intermolecular interactions are induced by Van der Waals interactions [13-15] and electrostatic

45

interactions, such as inter- or intramolecular hydrogen bonding effects,[16] consequently causing

46

either positive or negative deviation from Raoult’s law. Net attractive interactions result in an

47

activity coefficient, 𝛾, below unity, and net repulsive interactions result in a 𝛾 above unity. Activity

48

coefficients are composition dependent and a limiting case is where a solute is infinitely diluted in

49

a solvent. The non-ideal behaviour of the solute at infinite dilution is solely induced by solvent-

50

solute interaction, i.e. the effect of the molecular properties of the solvent on the activity coefficient

51

of the solute. The activity coefficient in this limiting case is called the infinite dilution activity

52

coefficient, 𝛾∞ 𝑖 .[12]

53

Various methods are in use to estimate activity coefficients. Eight models were evaluated, as

54

these are most commonly used. Seven of them having many similarities and one has a completely

55

separate theoretical framework. Three Solvation Models (SM) are chosen which are, in order of

56

increasing complexity, the Hildebrand parameter, Hansen Solubility Parameter (HSP) and the

57

Modified Separation of Cohesive Energy Density (MOSCED) model. These models attempt to

58

describe the intermolecular interaction strength by increasingly more molecular parameters. The

59

Group Contribution Methods (GCMs), are similar in form to SMs, but differ in origin. GCMs attempt

60

describing the intermolecular interactions by empirical fitting of binary interaction coefficients of

61

segments of the interacting molecules. The three GCMs that are included in this work, are the

62

original Universal Quasichemical Functional-group Activity Coefficients (UNIFAC) method, and

63

two modifications thereof (Lyngby and Dortmund). Despite these differences, all these models

64

still use the same entropic formulations. Next to these, COSMO-RS is also evaluated, which is a 2 ACS Paragon Plus Environment

Page 3 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


65

software package that is called Conductor-like screening model for realistic solvents. COSMO-

66

RS has an entirely different theoretical framework and does not require any input parameters from

67

the user besides the molecular structures. Similar to GCMs, COSMO-RS also divides molecules

68

in segments. COSMO-RS divides the surface of molecules, as calculated by the quantum

69

chemical Density Functional Theory approach in segments, and calculates for each segment the

70

interaction energy. The molecular properties are then calculated by taking the integral over all

71

segments.

72

Group Contribution Methods (GCM) such as variations on UNIFAC, and COSMO-RS are

73

most commonly applied, as can be seen in Figure 1. This, however, does not imply that these

74

models are always most accurate in predicting the 𝛾∞ 𝑖 . Because the existing literature [6, 17-24]

75

is inconclusive, the aim of this contribution is to extensively compare the performance of all these

76

fundamentally different approaches for prediction of the 𝛾∞ 𝑖 of (a)-polar solutes in (a)-polar

77

solvents. The performances of all evaluated models in predictions for a variety of chemical

78

systems were evaluated and explained on the basis of their fundamental assumptions. Examples

79

of such assumptions include that the entropy of the system does not differ from the ideal entropy,

80

that the volume of the molecule does not change in a changing environment, or that there is no

81

distinction between hydrogen bond accepting and donating molecules.

82

The extensive evaluation of model performances yielded insight in the applicability of the

83

models for systems with variations in intermolecular interactions and which models give the most

84

accurate description of the 𝛾∞ 𝑖 at 298.15K. A heuristic approach for the model choice is given for

85

all binary combinations of solute and solvent classes, e.g. apolar compounds, aromatic

86

compounds, halogenated compounds, polar aprotic compounds and polar protic compounds.

87



Page 4 of 24

Figure 1: The number of publications where both the model and activity coefficient is mentioned extracted from Scopus. The search terms where “model” AND “infinite dilution activity coefficient” retrieved at the 3rd of December 2018. 88

2. Theoretical Framework

89

In this section, all the models that were compared for prediction of 𝛾∞ 𝑖 are described. The models

90

can be categorized into solvation models, group contribution methods, linear solvation energy

91

relations, and COSMO-RS predictions that are of statistical thermodynamic nature. Because both

92

the solvation models and the group contribution methods make use of combinatorial and residual

93

contributions that find their origin in the Flory-Huggins model, this model and the variations

94

thereon are first discussed.

95 96

2.1. Flory-Huggins-based models

97

Non-ideal behaviour in mixtures can be induced by intermolecular interactions such as polarity,

98

as well as by shape differences,[16] The simplest cause of non-ideal behaviour is due to shape

99

differences without polarity differences. This situation occurs in alkane mixtures for which activity

100

coefficients can be described by the Flory-Huggins theory (eq. (1)) where the combinatorial

101

contribution to the activity coefficient is formulated solely in terms of molecular volume differences.

102

[25, 26]

103 ln 𝛾𝑐𝑖 = ln

() 𝛷𝑖 𝑥𝑖

+1―

𝛷𝑖 𝑥𝑖

(1)

104 105

Where 𝛷𝑖 and 𝑥𝑖 are the volume fraction and molar fraction of compound i, respectively. 4 ACS Paragon Plus Environment

Page 5 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


106 107

The Flory-Huggins combinatorial approach assumes a very large number of nearest

108

neighbor sites, hence ignores the fact that neighboring sites can be occupied by a segment of the

109

same molecule. Consequently, the Flory-Huggins correlation overestimates the combinatorial

110

contribution. [27] The Stavermann-Guggeheim modification attempts to correct this by

111

incorporating the probability of vacant sites for polymer segments, though still the combinatorial

112

term is overestimated because the coordination number of all molecules is set.[28] The Kikic

113

modifications attempted to correct the correlation by adding an exponent to the number of lattice

114

sites,[29] but this resulted in an underprediction of the combinatorial term. Recently, Krooshof et

115

al. [28] generalized the approach of Guggeheim and set loose the fixed coordination number. For

116

this approach, the coordination number of all molecules in the system has to be additionally

117

specified.

118

When there is also a difference in the polarity of the molecules in the mixture, a residual

119

correction can be added. This can be described by the Flory-Huggins free energy parameter, 𝜒12,

120

as equated in eq. (2). The activity coefficient resulting from both combinatorial and residual terms

121

is given in eq. (3). ln 𝛾𝑅𝑖 = 𝜒12 𝛷2𝑖

(2)

ln 𝛾𝑖 = ln 𝛾𝑐𝑖 + ln 𝛾𝑅𝑖

(3)

122

Starting from eq. (3), a wide range of models have been developed that vary in their formulation

123

of 𝛾𝑐𝑖 and 𝛾𝑅𝑖. In the next subsections, the most relevant solvation models and group contribution

124

methods are described.

125 126

2.2. Solvation Models (SM)

127

In an early attempt to describe and understand the strength of intermolecular interactions,

128

Hildebrand [7] defined the cohesive pressure or cohesive energy density, c, as the net result of

129

the sum of all intermolecular interactions between the molecules. As a measurable quantity for

130

the cohesive energy density in liquids below their boiling point, the molar vaporization energy,

131

ΔUvap or enthalpy, ΔHvap, is considered, [30] and correlated to the Hildebrand parameter, 𝛿,

132

according to eq. (4).

133 𝛿=

𝑐=

―

∆𝑈𝑣𝑎𝑝 𝑖

𝑉𝑚

=

∆𝐻𝑣𝑎𝑝 ― 𝑅𝑇 𝑉𝑚

(4)

134 135

Where iVm is the molar volume of the molecule i.

136 5 ACS Paragon Plus Environment


Page 6 of 24

137

The cohesive pressure is the sum of all attractive interactions which needs to be broken to

138

vaporize the liquid, and large 𝛿 values are therefore obtained for highly polar substances and

139

small 𝛿 values are obtained for weakly interacting substances, e.g. fluorocarbons. [30]

140

Considering mixtures, the difference in 𝛿 of the constituents of that mixture can be interpreted as

141

the difference in nature of these molecules and may be used as a measure for non-ideality. In the

142

Hildebrand-Scatchard equation (eq (5)),[7] the difference in 𝛿 is used in the expression for the

143

Flory-Huggins free energy parameter, 𝜒12.

144 𝑗

𝜒12 =

𝑉𝑚 𝑅𝑇

(𝑖𝛿 ― 𝑗𝛿)2

(5)

145 146

where 𝑖𝛿 and 𝑗𝛿 are resp. the Hildebrand parameters of the solute and solvent.

147 148

Hansen [8] proposed an extension of the Hildebrand parameter by separating the dispersion (𝛿𝐷),

149

polar (𝛿𝑃) and hydrogen bonding (𝛿𝐻𝐵) contribution. These three parameters are called the

150

Hansen Solubility Parameters (HSP) and are linked to the Hildebrand parameter as defined in

151

Equation (6).

152 𝑖 2

2

𝛿 = 𝑖𝛿𝐷 + 𝑖𝛿2𝑃 + 𝑖𝛿2𝐻𝐵

(6)

153 154

Often, the dispersive component was determined by homomorph methods, which estimates the

155

dispersive parameter by evaluating an apolar molecule with nearly the same size and shape of

156

the polar compound. [30] The remainder can be subtracted from the Hildebrand parameter and

157

split into the polar and hydrogen-bonding term. By optimizing the miscibility description an optimal

158

split between both terms was chosen. [30]

159

Common practice in using HSP is to plot the solute and solvents in a 3D Hansen space.

160

The spatial distance between solute and solvent can be correlated to the solvation capability of

161

the solvent, where shorter distances in the Hansen space allow better solubility. Hansen [31] also

162

suggested that the Flory-Huggins parameter can be determined using eq. (7).

163 𝑖

𝜒12 = 𝛼

𝑉𝑚

((𝑖𝛿𝑑 ― 𝑗𝛿𝑑)2 + 0.25(𝑖𝛿𝑝 ― 𝑗𝛿𝑝)2 + 0.25(𝑖𝛿𝐻𝐵 ― 𝑗𝛿𝐻𝐵)2) 𝑅𝑇

(7)

164 165

Where, initially α was taken to be 1, though Lindvig et al.[9] showed an α value of 0.6 increases

166

the average accuracy of the model. 6 ACS Paragon Plus Environment

Page 7 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


167 168

The Modified Separation of Cohesive Energy Density (MOSCED) model may be one of the most

169

extensive Solvation Models.[21] In the MOSCED model, additional contributions to the Flory-

170

Huggins parameter are considered that arise from significant variations in the cybotactic region

171

due to the local organization as a result of electrostatic interactions such as hydrogen bonding.

172

This local organization causes the geometric mean assumption not to be valid anymore for highly

173

polar and associating compounds.[21] To account for hydrogen bonding, the MOSCED model

174

distinguishes acidic (𝛼) and basic (β) contributions to hydrogen bonding. Similar to the HSP

175

model, the summation of the terms results in the Hildebrand parameter, as can be seen in

176

equation (8).[32]

177 𝑖 2

(8)

𝛿 = 𝑖𝜆2 + 𝑖𝜏2 + 𝑖(𝛼β)2

178 179

Where the dispersion constant 𝑖𝜆 and polarity constant 𝑖𝜏 are identical to the HSP parameters 𝑖𝛿𝑑

180

and 𝑖𝛿𝑃, respectively. In addition, the interaction between induced dipoles is accounted for by the

181

induction parameter, 𝑖𝑞. The model furthermore contains two empirical asymmetry factors, 𝜓 and

182

𝜉. More information on the MOSCED parameters can be found eqs. S1-S8 in the Electronic

183

Supplementary Information (ESI).[6] The resulting MOSCED-based equation for the Flory-

184

Huggins parameter is given in eq. (9).

185 𝑖

𝜒12 =

𝑉𝑚

𝑅𝑇

(

(

(

𝑖 2𝑗 2 𝑖 𝑇 2

𝜆 ― 𝜆) +

𝑖

𝑗

𝑞 𝑞

𝜏 ― 𝑗𝜏

𝑖

𝜓

𝑇 2

( 𝑖𝛼𝑇 ―

)

𝑗 𝑇

𝛼

+

)( 𝑖𝛽𝑇 ―

𝑗 𝑇

𝛽

𝑖

𝜉1

)

)

(9)

186 187 188

2.3. Group Contribution Methods (GCM)

189

Also, the Group Contribution Methods of UNIFAC and modifications thereof are the sum of a

190

combinatorial and a residual part. Each GCM model uses a different description for the

191

combinatorial term. The UNIFAC and the mod. UNIFAC (Do) models use the Guggenheim-

192

Stavermann [33] term (eqs. (10) and (11)), while the mod. UNIFAC (Ly) uses the Kikic [29]

193

modification as given in eq. (12). 𝑈𝑁𝐼𝐹𝐴𝐶:

𝑚𝑜𝑑. 𝑈𝑁𝐼𝐹𝐴𝐶 (𝐷𝑜):

( (

ln 𝛾𝑐𝑖 = ln (𝛷𝑖) + 1 ― 𝛷𝑖 ― 5𝑞𝑖 1 ― ln 𝛾𝑐𝑖 = ln (𝛷′𝑖) + 1 ― 𝛷′𝑖 ― 5𝑞𝑖 1 ―

𝛷𝑖 𝜃𝑖 𝛷𝑖 𝜃𝑖

+ ln

+ ln

( )) ( )) 𝛷𝑖 𝜃𝑖

𝛷𝑖

(10) (11)

𝜃𝑖



Page 8 of 24

ln 𝛾𝑐𝑖 = ln (𝛷′′𝑖) + 1 ― 𝛷′′𝑖

𝑚𝑜𝑑. 𝑈𝑁𝐼𝐹𝐴𝐶 (𝐿𝑦):

(12)

194 195

where 𝛷, 𝜃 and 𝑞 are respectively the volume fraction, surface fraction and the coordination

196

number which is often set at 10. The modified UNIFAC (Do) model also has a modified volume

197

fraction, on which more information is given in the ESI.

198 199

The residual contribution of all UNIFAC models is determined via Equations (13) and (14).

200 ln 𝛾𝑅𝑖 =

∑𝜈

(𝑖) 𝑘

(13)

(𝑙𝑛 𝛤𝑘 ― 𝑙𝑛 𝛤(𝑖) 𝑘 )

𝑘

(

𝑙𝑛 𝛤𝑘 = 𝑄𝑘 1 ― ln

(∑ 𝜃

𝑚𝛹𝑚𝑘

𝑚

)

𝜃𝑚𝛹𝑘𝑚

∑∑ 𝜃 𝛹

―

𝑚

𝑛 𝑛

𝑛𝑚

)

(14)

201 202

(𝑖) where 𝛤𝑘, 𝛤(𝑖) 𝑘 ,𝜈𝑘 , 𝑄𝑘 and 𝛹 are respectively the overall activity of moiety k, the activity of moiety

203

k solely surrounded by moiety i, the occurrence of each moiety k in surrounded by moiety i, the

204

Van der Waals surface of group k and the group binary interaction parameter. [4, 5, 11]

205 206

2.4. Linear Solvation Energy Relationship (LSER)

207

Linearly combining solute and solvent descriptors was pioneered by Kamlet, Abboud and Taft

208

[34-38], and later Abraham [3] introduced a now widely used LSER comprised of five solute and

209

five solvent descriptors (eq. (15)). An extension is made for ILs, where both the cation and anion

210

have five unique solvent descriptors (eq.(16)).[39] The Abraham model can be used to obtain

211

either the gas-solvent partition coefficient, KS, or the water-solvent partition coefficient, Ps.

212

(𝑃 ) = 𝑐 + 𝑒𝑬 + 𝑠𝑺 + 𝑎𝑨 + 𝑏𝑩 + 𝑣𝑽 {log log (𝐾 ) = 𝑐 + 𝑒𝑬 + 𝑠𝑺 + 𝑎𝑨 + 𝑏𝑩 + 𝑙𝑳

(15)

𝑆

𝑆

{loglog((𝑃𝐾 ))==𝑐𝑐++(𝑒(𝑒 ++𝑒𝑒 )𝑬)𝑬++(𝑠(𝑠 ++𝑠𝑠 )𝑺)𝑺++(𝑎(𝑎 ++𝑎𝑎 )𝑨)𝑨++(𝑏(𝑏 ++𝑏𝑏 )𝑩)𝑩++(𝑣(𝑙 ++ 𝑣𝑙 )𝑳)𝑽 𝑆

𝑆

𝑐

𝑐

𝑎

𝑎

𝑐

𝑐

𝑎

𝑎

𝑐

𝑐

𝑎

𝑎

𝑐

𝑐

𝑎

𝑎

𝑐

𝑎

𝑐

𝑎

(16)

213 214

The capital variables in eqs. (15) and (14) are defined as the solute descriptors and the lowercase

215

variables are the solvent descriptors. While the c-variable is a fitting constant [-], the latter terms

216

are respectively the excess molar refraction [dm3·mol-1/10], the dipolarity/polarizability [-],

217

hydrogen-bond acidity and basicity [-], the McGowan characteristic volume (dm3·mol-1/100) and

218

the gas-hexadecane partition coefficient at 298.15K [-]. The descriptors describe the tendency to

219

interact via σ- and π-electrons (e/E), the tendency to interact with (induced) multipole moments

220

(s/S), the tendency to accept hydrogen bonds or donate electron (A/b), the tendency to donate 8 ACS Paragon Plus Environment

Page 9 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


221

hydrogen bonds or accept electron (a/B) and the tendency to either form (V/L) or the work

222

required to form (v/l) cavities.

223

The determination of solute and solvent descriptors is done by multilinear regression of

224

experimental 𝛾∞ 𝑖 data of either a set of solutes in a solvent or a solute in various solvents. The

225

gas-solvent partition coefficient can consequently be linked to the 𝛾∞ 𝑖 by eq. (17).[39] log (𝛾∞ 𝑖 ) = log

226

( ) 𝑅𝑇

𝑃𝑜𝑗 𝑗𝑉𝑚

― log (𝐾𝑠)

(17)

where 𝑃𝑜𝑗 and 𝑉𝑗 are resp. the vapor pressure and molar volume of the pure solvent.

227 228

2.5. Conductor like Screening Model for Real Solvents (COSMO-RS)

229

The ab initio method developed by Klamt et al. [40], called COSMO-RS, predicts chemical

230

potentials, which can be used to calculate the 𝛾∞ 𝑖 . The first step is always to perform quantum

231

mechanical calculations to obtain the state of the geometrically optimized molecule. This step

232

needs to be done only once and the result can be stored in a database. The second step estimates

233

the interaction energy of the optimized molecules with other molecules and can henceforth

234

estimate molecular properties such as the 𝛾∞ 𝑖 . [40]

235

The energy of this state can be determined via dielectric continuum solvation methods.

236

Empirical parameters, e.g. atomic radii, are however required to construct the molecular cavity

237

within the conductor exterior. COSMO-RS implements element-specific radii which are ~17 %

238

larger than Van der Waals radii. The state of the molecule is consequently determined using any

239

Self-Consistent Field (SCF) model, e.g. Hartree-Fock (HF) and Density Functional Theory (DFT).

240

Combining the COSMO approach with a SCF model results in not only the total energy of the

241

molecule, but also the polarization charge density, or σ. [40]

242

The -profile is “frozen” into place, while the conductor environment is “squeezed out”. A

243

thin film of a conductor is left in place where different molecules will have an interface. The sum

244

of the σ-value at the interface, (σ + σ’), is then the net polarization charge density. As the

245

conductor is removed, the polarization charge density differences resemble intermolecular

246

interactions with local contact energy, which may be described by; 𝐸𝑐𝑜𝑛𝑡𝑎𝑐𝑡(𝜎,𝜎′) = 𝐸𝑚𝑖𝑠𝑓𝑖𝑡(𝜎,𝜎′) + 𝐸ℎ𝑏(𝜎,𝜎′)

(18)

247

An interaction distinction is made between misfit energy, 𝐸𝑚𝑖𝑠𝑓𝑖𝑡, and due to hydrogen bonding

248

effects, 𝐸ℎ𝑏. Both energy terms dissipate when the conjugant polarization charge densities are

249

equal. If not, the misfit between both σ-values represents the electrostatic interaction energy

250

between those segments. Additional interaction energy can be induced by two very polar surfaces

251

with opposite signed via hydrogen bonding. The capability of COSMO-RS to predict a large



Page 10 of 24

252

number of chemical potentials of solutes in either a pure or mixed solvent enabled fast and

253

versatile predictions regarding various equilibria and also 𝛾∞ 𝑖 .[40]

254 255

3. Methods

256

For all 𝛾∞ 𝑖 predictions a systematic assessment was done at 298.15K and all model specific

257

parameters were imported from literature sources, as can be seen in the ESI. All simulations with

258

COSMO-RS where performed with COSMOthermX C30_1705 in which a pure solvent phase was

259

defined and the activity coefficient of the solute in that phase was estimated. Molecules that were

260

not available in the databases, where created with TurboMole TmoleX 3.4. using the TZVPD-Fine

261

parameterization.

262 263 264

4. Results 4.1. General averaged model predictions

265

The accuracy of all models was determined by comparing the predictions with experimental

266

values taken from literature. An extensive overview of all experimental data is presented in the

267

ESI. The accuracy of the models is evaluated by determining the average relative deviation (ARD)

268

between the predicted and experimental 𝛾∞ 𝑖 values, as given in equation (19).

269

|

∑𝑁 𝑖=1 𝐴𝑅𝐷 =

∞ 𝛾∞ 𝑖 𝑚𝑜𝑑𝑒𝑙 ― 𝛾𝑖 𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙

𝛾∞ 𝑖 𝑒𝑥𝑝𝑒𝑟𝑖𝑚𝑒𝑛𝑡𝑎𝑙

|

(19)

𝑁

270 271

The ARD for prediction of 𝛾∞ 𝑖 for all solutes in molecular solvents and ILs by each model and their

272

95% confidence interval are depicted in Figure 2. For molecular solvents, eight predictive models

273

have been evaluated, for ILs seven predictive models due to lack of MOSCED parameters for ILs.

274

The ARD of the various models differs significantly, as can be seen in Figure 2. For both the

275

molecular solvents and the ILs, the most inaccurate model is the Hildebrand model. Evidently,

276

using only the evaporation enthalpy and the molar volume is insufficient to accurately describe

277

𝛾∞ 𝑖 in systems where intermolecular interactions such as hydrogen bonding take place.[41] This

278

accumulates in a significant ARD of > 105% in the 𝛾∞ 𝑖 prediction.


Page 11 of 24 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: The evaluation of various predictive models for 𝛾∞ 𝑖 for (a) molecular solvents and (b) ionic liquids at 298.15K. On the y-axis the ARD is presented with in the boxes the total amount of comparisons made. The experimental 𝛾∞ 𝑖 is collected in the Supplementary information. The integrated scatter plot depicts similar comparison made in literature for various models, e.g. mod. UNIFAC (Ly)[17-19], UNIFAC[17-22, 42, 43], COSMO-RS[18, 44], mod. UNIFAC (Do) [17, 19, 22, 23, 45-47] and MOSCED[6, 20-22, 24]. 279 280

Using the HSP to calculate the 𝛾∞ 𝑖 with eq (7) is a significant improvement compared to

281

the Hildebrand model (eq (5)), which is due to taking into account hydrogen bonding and polarity

282

effects.[8] Still, an ARD of 66.4±14.4% is observed for molecular solvents. This may be explained

283

by the inability to differentiate between hydrogen acidity and basicity effects.[6] Further refining

284

the model using hydrogen acidity and hydrogen basicity, as well as polarizability effects as taken

285

into account in the MOSCED model, which shows from Figure 2a to be on average the most

286

accurate model with an ARD of 16.2±1.35%. Unfortunately, no parameters for ILs are available,

287

and therefore this model could not be evaluated for ILs. The three group contribution models that

288

were also evaluated showed comparable ARD of 32.2±1.84%, 31.1±1.66% and 24.3±1.63% for

289

respectively the mod. UNIFAC (Ly), UNIFAC and mod. UNIFAC (Do) for molecular solvents. The

290

ARD of COSMO-RS is with 28.3±1.07% similar to the GCM methods. The Abraham model is

291

more accurate than GCM methods and COSMO-RS with an ARD of 21.7±1.19%. The better

292

accuracy of the Abraham model is due to a more elaborate description of the various

293

intermolecular interactions via all descriptors.

294

For ILs the ARD of all models, see Figure 2b, increase due to the presence of not only

295

dispersion, dipole and hydrogen bonding interactions, but also ionic interactions between the ionic

296

species and the solute. The ARD of the HSP-model is determined to be 168±54.5%, while the

297

GCM methods perform better with an ARD of 86.2±14.6%, 86.5±15.7% and 122±55.9% for

298

respectively the mod. UNIFAC (Ly), UNIFAC and mod. UNIFAC (Do). Where COSMO-RS 11 ACS Paragon Plus Environment


Page 12 of 24

299

performed for molecular solvents comparable to the GCMs, for ILs the ARD is larger, 182±16.7%.

300

The larger ARD for ILs is known to be (at least partly) caused by neglecting long range ion-ion

301

interactions and insufficient description of extreme polarization charge densities of ions.[10] The

302

Abraham model performs most accurately for ILs with an ARD of 65.1±4.50%. The accuracy of

303

the Abraham model is interesting for solvent screening purpose, because of the availability of ion-

304

specific Abraham parameters for 60 cations and 17 anions, allowing rapid assessment 1020 ILs

305

as they are binary combinations of these ions. (see ESI)

306

The ARDs reported in Figure 2 appear to be larger than errors described in various literature

307

sources.[6, 17-24, 42, 43, 45-48] Similar comparisons have been made by i.a. Gmehling et al. [17]

308

who assessed the accuracy of the UNIFAC models, and Thomas et al.[20] who assessed the

309

accuracy of the MOSCED model. They reported lower ARD, correspondingly 25.8% (instead of

310

31.1±1.66%) and 9.10% (instead of 16.2±1.35%). Because the error calculation method is the

311

same, and only the used dataset differs, the logical conclusion is that the dataset used in this

312

work includes 𝛾∞ 𝑖 with an average higher error margin. This will be the case for each of the

313

assessed models, hence, these error margins in the data set will not affect the comparison of the

314

relative accuracies between all models. For comparison with Gmehling et al. [17], there is a

315

difference in the selected data, because they state that

316

excluded, whereas in this manuscript these values where included.

𝛾∞ 𝑖 -values larger than 100 where

317 318

4.2. Molecular Solvents

319

Although from the averaged model predictions one general guideline for model selection can be

320

distilled, a more detailed analysis of subgroups of solutes and solvents allows to provide a more

321

sophisticated directive to the use of thermodynamic models for the prediction of 𝛾∞ 𝑖 for various

322

specific chemical families. To this end, all solvents and solutes were classified into five categories,

323

i.e. aliphatic compounds, aromatic compounds, compounds containing a halogen atom, polar

324

aprotic compounds and polar protic compounds, and the model accuracies were evaluated per

325

combination of solvent and solute class.

326

In Table 1 all model evaluations are shown per combination of solvent and solute

327

categories. Comparing the Hildebrand, HSP and MOSCED model, it clearly shows that models

328

with increasing complexity, i.e. taking hydrogen bonding and polarity into account, and describing

329

the hydrogen bonding basicity and acidity separately, as well as including polarizability predict 𝛾∞ 𝑖

330

with increasing accuracy. All of these models show the largest ARD for protic compounds,

331

including amines, alcohols and aldehydes, which is a logic result due to the number and type of

332

intermolecular interactions occurring in these systems.

333 12 ACS Paragon Plus Environment

Page 13 of 24

334 335 336 337

Table 1: The accuracy of eight predictive models differentiated towards aliphatic, aromatic, halogen, aprotic and protic compounds in molecular solvents. All 25 binary solute-solvent combinations have been made at 298.15 K. The colors are indicative: ARD300% >100% 73,7% >100% >100%

Protic polar >105 % >500% >100% >100% 67,7%

Aliphatic

Aromatic

Halogen

22,3% 31,1% 32,9% 47,3% >100%

43,3% 18,2% 25,1% 64,3% 57,7%

44,7% 20,4% 30,5% 47,1% 62,2%

MOSCED Model

Solvent

Aliphatic Aromatic Halogen Aprotic polar Protic polar

Aliphatic

Aromatic

Halogen

7,4% 14,3% 18,1% 15,7% 24,6%

6,7% 19,6% 9,8% 14,8% 16,7%

7,7% 9,7% 18,3% 15,8% 32,7%

Protic polar 89,1% >100% 37,3% 37,3% 19,7%

Aliphatic

Aromatic

Halogen

17,8% 11,8% 34,6% 21,2% 20,6%

11,8% 20,0% 31,8% 21,2% 17,5%

6,8% 13,0% 16,1% 21,7% 22,3%

Solvent

Aprotic polar 46,3% 22,9% 34,9% 22,4% 39,5%

Protic polar 12,5% 27,9% 31,7% 31,7% 34,1%

Aprotic polar 47,6% 32,4% 25,6% 27,9% 36,3%

Protic polar 45,6% 24,9% 38,4% 38,4% 23,4%

UNIFAC

Solute

Solute

Aliphatic

Aromatic

Halogen

Aprotic

8,6% 10,7% 27,6% 32,7% 36,3%

15,4% 13,2% 36,6% 34,4% 32,0%

20,4% 26,7% 15,7% 32,5% 35,0%

35,6% 24,4% 31,4% 19,9% 38,1%

Protic polar 45,2% 56,9% 33,4% 33,4% 23,6%

Aliphatic

Aromatic

Halogen

12,4% 16,6% 34,8% 40,3% 35,8%

11,0% 18,9% 19,4% 23,7% 28,3%

26,6% 20,9% 18,0% 23,0% 32,7%

mod. UNIFAC (Ly)

mod. UNIFAC (Do)

Solute


Protic polar 84,5% 75,5% 62,4% 62,4% 27,8%

Solute Aprotic polar 39,5% 24,2% 23,2% 15,0% 29,8%

COSMO-RS


Aprotic polar 64,2% 28,7% 70,0% 39,8% >300%

Abraham Model

Solute

Solvent

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Solute

Aliphatic

Aromatic

Halogen

Aprotic

9,4% 24,9% 34,0% 41,0% 39,3%

11,9% 19,4% 24,7% 20,8% 41,3%

21,7% 23,9% 19,9% 23,7% 34,2%

47,9% 27,3% 29,1% 23,0% 37,5%

Protic polar 43,6% 43,1% 43,0% 43,0% 16,8%

Aliphatic

Aromatic

Halogen

7,1% 17,6% 32,0% 30,1% 23,2%

5,6% 17,0% 13,5% 17,4% 18,9%

10,6% 17,1% 17,5% 25,9% 28,2%

Aprotic polar 56,8% 26,2% 27,9% 23,6% 43,4%

Protic polar 18,3% 47,8% 35,3% 35,3% 28,5%

338 339

There is no single model predicting 𝛾∞ 𝑖 most accurate for all solvent and solute category

340

combinations. Each of the models, except for the Hildebrand and the HSP, predicts most

341

accurately a category of solute-solvent pairs. COSMO-RS performs best in systems with only

342

(induced) dipole interactions and in absence of hydrogen bonding formation. When polarity and

343

hydrogen bonding systems are concerned, COSMO-RS becomes less accurate. The MOSCED

344

model and Abraham model perform much better in hydrogen bonding systems, due to the multiple

345

parameters that describe these directional interactions. The various UNIFAC models appear to

346

be most accurate for a few categories of chemical interaction systems, which may arise from the

347

empirical nature of the UNIFAC models based on fitting of model parameters to experimental 13 ACS Paragon Plus Environment


Page 14 of 24

348

data. The variation between the UNIFAC models can also arise from the different fitting

349

procedures for the determination of their empirical constants. Finally, the differences in the

350

formulation of their combinatorial term can induce variation in activity coefficient prediction.

351 352

4.3. Ionic liquids

353

Overall, the ARD in predicted 𝛾∞ 𝑖 in ILs is larger than in molecular solvents. Not only the additional

354

electrostatic intermolecular interaction of the charges on the ions with the solutes are responsible

355

for this, but also the additional competition between the solute-related intermolecular interactions

356

and the interactions between the ions in the IL plays an important role. Furthermore, ILs with a

357

cation containing besides the central ionic moiety also a second moiety, e.g. ether, hydroxyl or

358

unsaturated bond, are collectively evaluated as functionalized ILs or FILS.

359 360

4.3.1. Cations

361

A systematic analysis was made for various classes of cations, and the cation class-specific

362

model performances are listed in Table 2. A larger ARD is generally obtained for FILs, due to the

363

fact additional inter- and intramolecular interactions occurs with the moieties present on the cation

364

tails. Also clearly can be seen that COSMO-RS has severe difficulties in predicting accurate 𝛾∞ 𝑖 .

365

Analog to the trends seen in molecular solvents, the ARD increase from apolar to polar solutes

366

due to hydrogen bonding effects. Overall, the Abraham model performs superior over the other

367

models, though some systems can be more accurately described using a variant of UNIFAC. For

368

instance, the predictions for aliphatic, aromatic and halogen solutes in imidazolium cations are

369

more accurately predicted with mod. UNIFAC. (Do). This is most likely due to the large dataset

370

available for imidazolium cations, hence improving the empirical fit of the mod. UNIFAC (Do).

371 372 373 374 375 376

Table 2: The accuracy of five predictive models differentiated towards aliphatic, aromatic, halogen, aprotic polar and protic polar compounds in ionic liquids with various non-functionalized and functionalized cations. All 25 binary solute-solvent combinations have been made at 298.15 K. The colors are indicative: ARD

Model Performances Evaluated for Infinite Dilution Activity Coefficients

Recommend Documents