511

# COSMO-SAC: A Predictive Model for Calculating Thermodynamic Properties on a-priori Basis

Authored by: Anand Bharti , Debashis Kundu , Dharamashi Rabari , Tamal Banerjee

# Phase Equilibria in Ionic Liquid Facilitated Liquid–Liquid Extractions

Print publication date:  March  2017
Online publication date:  March  2017

Print ISBN: 9781498769488
eBook ISBN: 9781315367163

10.1201/9781315367163-3

#### Abstract

The prediction of properties of liquid-phase thermodynamics is conventionally performed using the Gibbs energy of solvation (ΔG solv), which is merely the energy required to transfer a solute molecule from vacuum to a solvent phase. Two of the important properties computed from this route are the activity coefficient and the vapor pressures.

#### 3.1  Liquid-Phase Thermodynamics

The prediction of properties of liquid-phase thermodynamics is conventionally performed using the Gibbs energy of solvation (ΔG solv), which is merely the energy required to transfer a solute molecule from vacuum to a solvent phase. Two of the important properties computed from this route are the activity coefficient and the vapor pressures.

3.1 ()

In the equation,

is the Gibb’s free energy of solvation of molecule i in solvent j, while refers to the term for molecule i in itself. Here, the term refers to the non-translational contribution of the Gibb’s free energy of solvation. This is related by Equations 3.2 and 3.3, as given below:

3.2 ()

3.3 ()

Here, Equation 3.3 refers to the De Broglie’s wavelength originating from the dual particle theory of material, where m i is the mass of the particle, h is the Planck’s constant (J.s) and k is the Boltzmann constant (J/K). In a similar manner, the vapor pressures can be predicted by Equation 3.4, where an ideal phase has been assumed.

3.4 ()

It should be noted that the thermodynamic properties, namely partition coefficient, heat of vaporization and phase equilibria, are computed from

and . Hence, the objective of COSMO-SAC (Lin and Sandler, 2002) is to determine these properties from quantum chemical and statistical mechanical framework. The model is based on the original pioneering work of Klamt (Klamt, 1995; Klamt, Jonas, Bürger, & Lohrenz, 1998; Klamt & Schuurmann, 1993). The primary contribution from COSMO-SAC is that . is mainly bifurcated from two contributors, namely electrostatic (ES) and van der Waals (vdW) force (Lin, Chang, Wang, Goddard, & Sandler, 2004; Wang, Lin, Chang, Goddard, & Sandler, 2006; Wang, Sandler, & Chen, 2007; Wang, Stubbs, Siepmann, & Sandler, 2005). Contributions such as those from vibrational, rotational and dipole–dipole moment are ignored. Let us write out the contribution for each part in Equation 3.5.

3.5 ()

is the Helmholtz energy of molecule i in the solvent. Assuming the fact that > , we further write out the Helmholtz energy as a sum of two parts, namely and . The former refers to the London Dispersion term. This is usually referred to as temporary attractive force, which occurs when electrons in two nearby atoms occupy positions such that they form temporary dipoles. This force is also called ‘induced dipole–dipole’ interaction. refers to the amount of energy required to construct a cavity of component i within a solvent S. This primarily depends on the shape and size of the molecule.

The predictions of

take a different pathway each for activity coefficients and vapor pressure. The activity coefficients are usually predicted from the difference of their properties in solvent and individual phase, as given below (Burnett, 2012):

3.6 ()

By rearranging the terms, we obtain:

3.7 ()

Both the phases considered are liquid; hence, the dispersion energies are almost similar and equal. In other words, their difference is fractionally small as compared with other terms. So, we neglect this contribution, that is, the second term of Equation 3.7. On similar lines, the last two terms reflect the differences in cavity formation and density. These are used through a combinatorial contribution, which is similar to that used in nonrandom two-liquid (NRTL) and UNI QUAsiChemical (UNIQUAC) models. Hence, the final expression takes the form (Equation 3.8):

3.8 ()

For the estimation of vapor pressures, we require the Gibb’s free energy of solvation for a single phase. Hence, the dispersion terms do not cancel out; this invariably implies that the cavity formation term in Equation 5.8 (second term) is not straightforward. The Equation 3.8 takes the form:

3.9 ()

#### 3.2  Estimation of Electrostatic Contribution *

Within the COSMO-SAC model, the largest contribution, as evident from Equation 3.8, is from the electrostatic contribution, which is due the inherent charges between molecules. The description of

is a tedious job and is usually computed with a continuum solvation model, as defined in quantum mechanics. Here, the solute molecule is first modeled to be surrounded by a van der Waal-like segmented surface. Thereafter, it is immersed in a conductor with infinite dielectric constant, such that it totally screens the molecule’s inherent electronic density. Density Functional Theory (DFT) is then applied on it to solve for the Hamiltonian, so as to determine the properties of each segment, including the charge and area associated with each segment. Thereafter, an ensemble-based theory along with a statistical mechanical framework is used, which relates the segment interactions to useful thermodynamic properties such as chemical potential. In a nutshell, this incorporates steps, namely geometry optimization of the molecule, construction of the vdW segmented surface, immersion of the molecule with the surface into a conductor and, finally, the computation of the segment properties via statistical mechanical framework. We shall now discuss all the steps in detail.

#### 3.2.1  Geometry Optimization of the Molecule

The geometry for any given molecule is obtained by minimizing the total energy with respect to the atomic coordinates. This is usually performed by solving the Schrödinger equation, so as to predict the energies. Initially, the structures are drawn in some suitable freeware, such as Avogrado, Molden or GaussView, to name a view. Thereafter, the optimization is done in vacuum by giving a suitable level of theory (HF/DFT/Semi-Empirical) and a basis set in commercial packages such as Gaussian and TURBOMOLE. This is done in vacuum by using where after the end of optimization, an absence of negative frequencies indicates a global minimum. This step also gives us the optimized geometry, its partial charges and the energies of the molecule in vacuum. This optimized geometry is then used for all subsequent calculations. In so doing, the effect of any changes in geometry between vacuum and solvent is assumed negligible.

An interesting question may be: Should the geometry optimization be carried out in a known solvent? However, the choice of solvent may be then biased either towards high dielectric constants such as water or towards low dielectric constants such as hexane. Since neither of them is applicable, universal phenomena such as vacuum are chosen. This also reduces the computational time, as optimization in vacuum is faster. Moreover, reported experimental data (Pye, Ziegler, van Lenthe, & Louwen, 2009) suggest that optimization in vacuum provides better reproducibility.

#### 3.2.2  Segmented Surface

The segment surface is created by a vdW-like surface, which is determined by the union of spheres drawn around each of the atoms in the geometry-optimized molecule. Here, each element corresponding to the particular atom is given a specific radius, with which the spheres are created. These radii, which are determined empirically, are called the ‘COSMO radii’ and these are roughly 20% more than their vdW radii. The COSMO radii for 11 elements, namely H, C, N, O, F, Si, P, S, Cl and Br, have been optimized by Klamt et al. (1998). For unknown elements, Klamt and Schuurmann (1993) have recommended a value 20% larger than the vdW radius. The optimized radii of the elements were obtained from 15,000 DFT calculations (Klamt et al., 1998). Given this extensive optimization, the commercial packages directly use these data in their database (Banerjee, Singh, & Khanna, 2006; Grensemann & Gmehling, 2005).

Once the molecular surface or the molecule is defined, it is then divided into tiny polyhedra (segments), each with its own position and area. Figure 3.1 shows a sample segmented sphere. Each segment can be either a pentagon or a hexagon, but for computational affordability, each is assumed circular, with a radius,

.

Figure 3.1   Sample surface segmentation into pentagons (left). The COSMO model approximates each of these by circles, with area equal to those of the original polygons (right).

In the initial DMOL3 implementation of COSMO, which was Klamt’s original DFT software, the average segment size is 0.24 Å2.

#### 3.2.3  Determination of Segment Properties

After the segmented surface is obtained by the union of all the vdW spheres outlined around each atom, where each atom possesses a ‘COSMO radii’, the molecule is then immersed in a homogeneous conductor having infinite dielectric constant (Klamt, 2005). In such a scenario, the screening charges from the conductor migrate to the segmented surface to negate the underlying molecular charge distribution. This is physically modeled by placing a screening charge from the conductor on each segment. The total energy of solvation (E solv) for the entire molecule is then determined from the interaction of the screening charges with the solute nuclei, the solute electronic density and within themselves. The segment properties are now calculated by defining the energy of solvation (E solv), which solves the interaction of the segments with the solute nuclei (Z) (first term of Equation 3.1), solute electron density (second term of Equation 3.1) and interaction within the segments (third term of Equation 3.1), as given below. q α represents the screening charge of αth segment.

3.10 ()

where:

Here B, C, D and r represent vector quantities and represent the contributions due to Coulomb’s law. Z α and r α denote the molecular charge and position of nucleus A, respectively. V α, possessing a unit of inverse of length, signifies the distance between the segment and the inherent charge density

of the molecule. D αα represents the energy required to provide a charge on each segment.

The expression for the self-energy of a segment was made by assuming the molecule as a sphere of radius R . In such a situation, the total charge Q is distributed homogenously across its surface. Hence, the electrostatic energy of the sphere is given by

3.11 ()

As per the COSMO protocol, the surface has to divide into n homogeneously similar segments. This implies that the energy then becomes the total of the sum of segment interaction energies and the segment self-energies; this is given by

3.12 ()

Here,

is the interaction energy between segments, while is the self-energy of segment. The latter reflects the energy required to maintain a charge of a given density on a particular segment. denotes the charge of segment α, while provides the distance between segments α and β. It is assumed that the area and charge of each segment are similar and the summations as above are exact. However, we still do not have an exact expression for , since the energy of the sphere must be the same and will not depend on the distribution of charge on the sphere. An explicit way to determine the same is to an expression for can be obtained by equating Equations 3.11 and 3.12. Therefore, by rearranging the above equations, we obtain (Burnett, 2012):

3.13 ()

By reorganizing the terms, we obtain:

3.14 ()

Here,

and is equal to , A is the surface area of the sphere and is the area of each segment. Thus, in a nutshell, Equation 3.14 is exact only when the charges are homogenously placed on top of the sphere. However, as the charge and area are known from a COSMO calculation, this dependence is limited to only the last part of the equation, that is, . The expression for this term is the highest (i.e. 1.250) at n = 4, while it steadily drops down to 1.078 at n = 20 and finally to 1.069 at n = 60. Hence, it has been recommended to keep this constant at n = 20 (i.e. 1.070) such that the final expression takes the form:

3.15 ()

The energy of solvation (E solv) is then minimized with respect to the segments’ charges (q α), as given below:

3.16 ()

As the molecule is placed on the continuum with infinite dielectric constant, the charges in the nearby region tend to cancel the electron density ρ( r ) of the solute molecule. This causes the molecule to polarize. Hence, this calculation needs to be performed in an iterative manner, so as to achieve consistency. At each step, the screening charge q is found for the respective ρ( r ), after which ρ( r ) is updated with new values of q . This procedure, which utilizes an iterative method, is called ‘self-consistent field’ approach. This process of generating and identifying the set of ρ( r ) and q that are self-consistent is called a ‘conductor-like screening model’ (COSMO).

In addition to q and area of the segments a and r , it also provides the potential for each segment φ. The calculation also provides us with the total energy of the molecule in a conductor, which we will denote as ( E COSMO).Sample screening charges for different molecules are depicted in Figure 3.2. It should be noted that the positive screening charges depict the negative regions of the molecule. Likewise, the negative screening charges represent the positive regions of the molecule. The screening charge distribution of the molecules is drawn by using the COSMO builder within TURBOMOLE. The package can also evaluate the segment potential for the segments, which we shall discuss in the ensuing section.

Figure 3.2   Gaussian03-generated COSMO surfaces for three molecules: Acetone (739 segments), water (292 segments) and toluene (1053 segments). Surfaces colored by screening charge density, where blue regions have negative screening charge values, green regions are neutral and red regions have positive values.

For the acetone COSMO file, ‘nps’ represents the number of segments, which is 739. The values of area2) = 371.79 and volume3) = 561.92 represent the total cavity area and volume, respectively, within the conductor, where acetone molecule is placed. The next section denotes the position of all the 10 atoms, along with the total screening charge on each atom ‘COSMOCharge’ and the last column represents the division of this ‘COSMOCharge’ with the area corresponding to each atom, called ‘area’. The last column is the most desirable quantity, as it reflects the screening charge density. It is interesting to note that except for oxygen atom ‘O3’ (COSMOCharge = 0.24046), all other atoms posses a negative charge. This follows the obvious trend, as negatively charged atoms possess positive screening charges, and vice versa. Further, the positive charge on the oxygen atom is then divided into a number of segments for which information is available in the concluding part of the COSMO file. For example, information for the segments of ‘C1’ and ‘H9’ atoms of acetone is depicted in Figure 3.3. The screening charges q of the solvent medium are usually calculated by a scaling q*, so that q = f(ε) q*, where f(ε) is the scaling factor given by

Figure 3.3   An excerpt of a COSMO file, as generated from Gaussian03 version C.02 for acetone molecule.

3.17 ()

Thus, a simpler boundary condition of a conductor appears in the above equation and also in Figure 3.3 (fepsi = 1.00). It should be noted that a conductor has an infinite supply of charge, so as to screen the entire solute molecule. Hence, there is a point of thought that the selection of solvent as a conductor can alter the charges, as evident in Figure 3.2. Therefore, a decision has to be arrived for the selection of a solvent in such a manner so that it is able to screen the entire solute component. This brings out to the point that it should possess sufficient charge, and thus, for obvious reason, it has to be a perfect conductor. Here, the COSMO approximation is exact in the limit of ε = ∞ and is within 0.5% accuracy for strong dielectrics such as water with a permeability of (ε = 80). Even for a lower dielectric limit of solvents such as ε = 2, COSMO coincides with the exact dielectric model within 10%. Hence, the COSMO approach is becoming a standard Continuum Solvation Model (CSM) in quantum chemical codes.

It should be noted that the self-consistent field (SCF) calculation or the COSMO scheme has to be done once for each component. Once these charges are obtained, then a framework has to be created so as to restore these charges when placed in an actual solvent of a known dielectric constant. This is where the derivation of chemical potentials and activity coefficients is performed. Our aim will then be to create a repository of COSMO files, which will include ions, solutes, solvents and inorganic compounds, subject to their dissociation constants. A sample description of an input for generating the COSMO file is given in Figure 3.4.

Figure 3.4   Input file for COSMO file generation in Gaussian03.

The different terms and their meaning are as follows:

‘% mem = 540 MW’ refers to the total amount of internal memory required to perform this calculation. The route section with ‘# P BVP86/SVP/DGA1’ indicates the following form ‘DFT Theory/Basis Set/optional’. So, the form of functional used for DFT is the P BVP86, while the basis set used is SVP or Split Valence Polarized. The orbital coefficients for SVP can be obtained by the GFprint and GFInput commands in Gaussian03. The density gradient approximation or ‘optional’ used is DGA1. It expands the electron density of the atoms in the form of atom-centered functions for saving computational time. It only integrates the coulomb interaction in place of all the electron integrals. The ‘SCF = tight’ indicates a full convergence of energy.

#### 3.2.4  Computation of the Restoring Free Energy ( Δ G i ​ / ​ s ∗ e s )

Here, the contribution to the restoring free energy can be divided into two parts, namely

and , where the former implies the electrostatic energy due to the component in the solution and the latter explains the change in energy when the molecule is taken from a conductor to a solvent. These are mathematically defined as below:

3.18 ()

Both the energies in COSMO conductor and vacuum optimization can be obtained as per the procedures mentioned in Sections 3.2.2 and 3.2.1, respectively. For example, the total energy in the conductor, as per Figure 3.3, is ‘Total energy corrected [a.u.] = −193.0450861463’.

represents the restoring free energy for the movement of the solute from the conductor to real solvent. For a real conductor, , which implies that . Hence, the major portion from now on will be the Segment Activity Coefficient (SAC) approach, which will involve the derivation of segment potential, as discussed in the next section.

#### 3.2.5  Segment Chemical Potential

Based on statistical thermodynamics, COSMO-SAC assumes the ensemble of molecules by an ensemble of independently interacting segments (Lin and Sandler, 2002; Burnett, 2012). Hence, for any mixture, segments from each component are mixed in proportion with the desired component mole fractions. The resulting interaction is then captured by the energy required to bring them from an infinite point in space. Let us consider an ensemble of k different types of molecules, each having ni segments. Hence, the total number of segments is

3.19 ()

Here, b represents the total number of molecules of type i and

represents the total number of segments of type α. For a single component, this reduces to . Assuming a constant geometry, where all the segments are allowed to move freely, implies that the movements of the segments are not dependent on neighboring segments. This ensemble then takes a form of mixture of independently interacting segments. We can thus propose any such formation as below (Burnett, 2012):

3.20 ()

Each row in the above matrix notation denotes the interaction of the first with the remaining segments (same component or of other component). Hence, the total neighbors will be of type ‘1’ and will be

such that . So, from the same logic, the number of pairs in the mixture will be

3.21 ()

Thus, the probability of obtaining a segment pair

will be

3.22 ()

Here,

is an ensemble averaged value approximated from a canonical ensemble. Hence, the whole mixture can be written in terms of a partition function Q(N, V, T) with a Boltzmann exponential term (see Appendix I for canonical partition function). This is given as below

3.23 ()

Here, Q and Q represent the partition functions of the mixture with and without including the segment pair, respectively. Here,

represents the interaction between the segment pair α and β. Thus, when α = β, then = 1, since both the segments become indistinguishable. The value of the degeneracy or becomes 2 when α#β. Equation 3.23 takes the following form when average probability is defined:

3.24 ()

The difference between Equations 3.23 and 3.24 is that they denote a segment pair and an individual segment, respectively. It is a known fact the Helmholtz energy for the ensemble having the two segments

takes the form:

3.25 ()

From statistical thermodynamics, we know:

3.26 ()

Hence, Equation 3.25 takes the form:

3.27 ()

Here,

are the chemical potentials for the two segments α and β, respectively. Hence, we have:

3.28 ()

Equation 3.24 then takes the form:

3.29 ()

The non-translational part of the chemical potential is given as

3.30 ()

This implies that the second part is the translational contribution for chemical potential. So, for ideal gases, the non-translational part is usually zero, or

3.31 ()

The final form of Equation 3.29 takes the form:

3.32 ()

For obtaining the average number of ensemble pairs, we thus use the relation as below:

3.33 ()

The total number of degrees of freedom for such a system can be shown to be equal to

; hence, we require the same number of equations so as to solve them and get a unique solution. Rewriting Equation 3.22 by diving it with ωαβ, we obtain:

3.34 ()

Here, the total probability is divided with the degeneracy and finally with total number of segment pairs. In a way, the above expression represents all possible probabilities of segment α with remaining segments across all the components c.

3.35 ()

Equation 3.35 points out to the fact that the total interaction with the remaining segments merely points out to the total number of pairs with segment α. Hence, we obtain n:

3.36 ()

Further, Equation 3.36 points out to the fact that a probability of locating a particular segment pair in the ensemble is merely locating a single pair or the segment mole fraction

.When applied to Equation 3.29, it takes the form:

3.37 ()

3.38 ()

By combining the above two expressions, we obtain:

3.39 ()

Rearranging the equation provides the following form:

3.40 ()

This can be rewritten as

3.41 ()

where:

3.42 ()

Thus, the equation is implicit in nature, as μα * appears on both sides of the equation. Hence, this needs to be solved iteratively. The E αβ possesses the form as given below:

3.43 ()

Here, a eff is the area of the effective segment and is an adjustable parameter.

is charge density of segments. Here, ces is a constant and is equal to . The equation above reflects the self-consistency approximation, as discussed earlier. Here, all the segments are assumed to be of the same size as per ‘area of the effective segment’; hence, the segment pair itself can be assumed to be an individual identity, with a total charge as the sum of the two charge segments. As the two molecules interact with each other, their electronic densities and subsequently, the wave function get disturbed. To capture this effect, the factor f pol is defined.

Now, we shall move to the most important aspect of thermodynamics, that is, the derivation of the chemical potential from the information available from the segment charge densities. In such a scenario, we shall assume the entire mixture of N molecules to be divided into number of segments, n 1, n 2nn . Assuming that we have large number of molecules implies that removal of a single molecule or segment will define another system. An analogy can be devised from an NPT ensemble, where a deduction of a particular molecule (moli) or a segment (segα) gives a pathway for calculating the chemical potential of the particular molecule or a segment, respectively.

Here,

refers to the chemical potential of the ith molecule from the mixture. Now, if we consider the same ensemble but in terms of segments (segα), we obtain (Burnett, 2012):

3.44 ()

So, by summing up the contributions from all the segments, we obtain:

3.45 ()

So, the total energy or the restoring free energy required will be equal to the sum of all segments in a mixture S. This is given as below:

3.46 ()

However, Equation 3.44 shall also be valid when the solvent itself is a perfect conductor. This will be applicable for the following condition:

3.47 ()

Here,

denotes the effective number of segments of component I, while gives us the effective number of segments of molecule I of type α. So, when it becomes a perfect conductor, then Equation 3.47 takes the form:

3.48 ()

Thus, we have to subtract the term

from the chemical potential to obtain the non-translational contribution.

3.49 ()

This merely reflects the fact that for a perfect conductor, the restoring free energy is zero and the chemical potential is due to its translational contribution. Equation 3.47 is based on the premises that all segment interactions are possible. In order to assume such information, we need to neglect the three-dimensional segment information. Further, no attempt has been made to consider the steric hindrance or the long-range liquid structure. This may be termed the most limiting assumption in the model. On the contrary, this assumption only makes the model simplistic in nature, as it eliminates the need to explicitly sample the full configuration space of the entire ensemble.

In the stated COSMO calculation, the surface of the molecule is dissected in small segments and screening charges are determined for each segment, such that the net potential everywhere at the surface is zero (perfect screening). These charges are averaged in order to make the predictions uniform and also to capture the different shapes and sizes of the segments uniformly. Thus, to obtain apparent screening charges to be used in COSMO-SAC model, the averaging algorithm is proposed as below:

3.50 ()

where:

• is the radius of the standard surface segment
• is the radius of the segment n
• dmn is the distance between segments m and n

The averaging of the screening charge densities, therefore, creates a new contribution to the overall Gibb’s free energy. This is provided as below:

3.51 ()

Here, the third term on the right-hand side is the Gibb’s free energy change occurring due to the averaging process.

3.52 ()

Thus, the first term, that is, the restoring free energy, is a function of the solvent, while the remaining terms are dependent on the molecule configuration. Here, E diel is given by

3.53 ()

is directly obtained from COSMO output file and is the energy needed for dielectric polarization, that is, the screening charge distribution. The factor has a constant value of f pol = 0.6917. These two free energies, namely the second and third terms of Equation 3.51, are pure species properties. So, they have the same values in ideal gas phase (i.e. in vacuum) and in the solvent. Mixture properties are thus obtained by subtracting free energy values of mixture and vacuum phase respectively. This finally leads to the cancellation for mixture property calculation. This is the reason that they are not included in the final expression in solvation free energy.

#### 3.2.6  Sigma Profiles of Simple Compounds

The three-dimensional screening charge density distributions, which are obtained as per the procedure given in Figure 3.3, are first averaged using Equation 3.50. This gives us charges ranging from –0.03 e/Å2 to +0.03 e/Å2. These segment charge densities are grouped using a histogram, which is commonly known as σ-profile (p(σ)). This essentially is the probability of locating a surface segment with screening charge density σ by the expression:

3.54 ()

where:

• Ai (σ) is the surface area with a charge density of value σ
• ai is the total surface area of species i

For a mixture, the σ-profile is determined from the area-weighted average of contributions from all the components present in the mixture. This is given as

3.55 ()

Since the current theme of this book is on ionic liquids (ILs), we shall introduce sigma profiles for simple components and then move towards ILs. It should be noted that computation of segment potential is quite a cumbersome process, concerning the number of segments a molecule possesses. A different and easier approach would be to use the histogram known as σ-profile or p(σ). In such a scenario, Equation 3.50 takes the form:

3.56 ()

3.57 ()

Here,

refers to the charge density within the histogram. Now, it is known that the histograms lie between –0.03 e/Å2 and +0.03 e/Å2. So, if we allocate each histogram with a value of 0.001 e/Å2, we shall then reduce the whole sigma profile to 61 histograms. Hence, N bins takes the number 61. In the above altered equations, xα/i and x α/S are invariably replaced by and , respectively. The latter two determine the sigma profiles of the component in itself (Equation 3.54) and in the mixture (Equation 3.55), respectively. However, x α is segment mole fraction, while is based on the segment area fraction. Notwithstanding as Equation 3.47 is not derived, the difference is usually very small when computing from .

It is interesting to note that the exchange energy (E pair) contains contributions from electrostatic interactions. One of them is referred as misfit energy Emf , which is the energy computed due to the difference in shape and size between a pair of segments. Another contribution, the hydrogen bonding energy Ehb , occurs when an electropositive atom (hydrogen) forms additional bonds with electronegative atoms such as oxygen, nitrogen and fluorine. The expressions of the two terms are given by Equations 3.58 and 3.59, respectively.

3.58 ()

3.59 ()

3.60 ()

Where, α is a constant for the misfit energy and is calculated from the surface area of a standard segment.

Here, chb is a constant for the hydrogen bonding interactions, while σ hb is cutoff value for hydrogen bonding interactions. When Equation 3.59 is used, all possible pairs are first grouped. Now, each pair is dealt with separately, having screening charge densities as σα and σβ. Then, σacc and σdon are allocated the larger and smaller values of σ m and σ n . Here max and min indicate the larger and smaller values of their arguments, respectively. Under this definition, the hydrogen bonding contribution is non-zero only if one segment has a negative charge density less than

and the other has a positive charge density greater than . In such a case, the contribution is limited to segment pairs of opposite charge possessing larger magnitudes. This contribution will also increase higher difference in charges densities between the pair of segments. The SAC then takes up the form:

3.61 ()

This equation is then solved iteratively. In the next section, we will focus our attention on the sigma profile and potential of the commonly used cations and anions.

#### 3.2.7  Sigma Profiles and Potentials of Ionic Liquids

The sigma profile and potential of cations and anions are obtained using Equations 3.55 and 3.56, respectively (Figures 3.5 and 3.6). For such an exercise, we have chosen a common cation, that is, 1-alkyl-3-methylimidazolium ([RMIM]; R=O, H, B etc.; O:Octyl; H:Hexyl; B:Butyl) (Figure 3.5) and then paired it up with three common anions, as given in Figure 3.7.

Figure 3.5   Sigma profiles of individual cations (EMIM: 1-ethyl-3-methylimidazolium; BMIM: 1-butyl-3-methylimidazolium; OMIM: 1-Octyl-3-methylimidazolium).

Figure 3.6   Sigma potential of individual cations.

Figure 3.7   Sigma profile of individual anions.

The sigma profiles and potentials for cations are shown in Figures 3.5 and 3.6, respectively. In all the normalized profiles, a small part of the sigma profile lies to the left of the hydrogen bonding cutoff radius (i.e. σ < –0.0082 e/Å2), indicating a hydrogen bonding donor capacity. The sigma profiles for [BMIM] and [OMIM] are of the same nature. Most of the prominent peaks lie in the negative side of the sigma profile. The negative screening charges are due to the positive charge residing inside the imidazolium ring. [EMIM], [BMIM], [OMIM], all show peaks at the outmost position in the negative direction. They also show up in the sigma potential (Figure 3.6), that is, negative values of chemical potentials are encountered on the positive side of the screening charge densities (SCDs). The negative values are encountered since extra free energy is gained by forming hydrogen bonds. It should be noted that extra energy is gained with the formation of hydrogen bonds; this energy is negative in value when compared with the free energy required (which is positive in nature) for removing the SCDs.

The sigma profiles and potentials for anions are shown in Figures 3.7 and 3.8, respectively. For the anions, it is evident that a negative contribution to the sigma potential is due to positive screening charges, as obtained from their sigma profiles. It should be noted that the contribution of chemical or segment potential is a sum of two values, that is, the misfit term and the hydrogen bond term. For the anions, a negative value on the profile implies that it is the hydrogen-acceptor bonding that dominates the left-hand side of Figure 3.7. However, this is the opposite in the donor region, where the peak rises sharply owing to the difference in the shape and size of the segments or the misfit contribution is known to dominate the total interaction.

Figure 3.8   Sigma potentials of individual anions.

For all the anions, peaks are lying on the right of the cutoff zone for hydrogen bonding, that is, σ hb = +0.0082 e/Å2. All the prominent peaks for the anions, that is, 0.0085 e/Å2 for [PF6], 0.014 e/Å2 for [CF3SO3], 0.015 e/Å2 for [EtSO4] and 0.011 e/Å2 for [BF4], lie to the right of the cutoff SCD, that is, σ hb = +0.0082. This is due to the inherent negative charges of the anions. The sigma potential is a similar one as observed for cations, except that negative values of chemical potentials are encountered on the negative side of the SCDs.

#### 3.2.8  Restoring Free Energy

Finally, we obtain the restoring free energy of molecule in solvent or mixture as

3.62 ()

This original definition was given by Klamt et al. (1998). Although this gave significant improvement for in hydrogen bonded systems, it was not fully consistent with the experimental observations. It should be noted that such policy made the hydrogen bonding contribution biased towards acceptors, which are limited to small electronegative atoms such as nitrogen, oxygen and fluorine. On the other hand, the donors originates from the hydrogen atoms which are in turn bonded to these electronegative atom(s). Hence, the process needs some course correction, as it leads to occasional meaningful errors. These weaknesses leave the process prone to occasional errors.

To incorporate a more meaningful hydrogen bond model, the definitions of acceptor and donor were modified. We assumed positively charged segments from the surfaces of all nitrogen, oxygen and fluorine atoms to be considered as acceptors. Similarly, the negatively charged segments from the surfaces of hydrogen atoms bound to any one of the acceptor atoms were chosen as HB donors. The modification will be taken up with test cases in Chapter 5 of this book. Now let us revert to Equation 3.7, which is given below:

3.63 ()

Here, we have obtained relation for the electrostatic contribution (first part of Equation 3.7, as above).

Now, we move our attention to the second term, that is, the dispersion or the van der Waals expression. This is the contribution due to the Helmholtz free function, which is modeled using a first order mean-free approach, the perturbation approach. For pure component, the expression takes the form:

3.64 ()

Here,

(K Å3) is the pair interaction energy between atom type j and atom type k, while is the effective number of atoms of type j within a species i molecule, which is obtained from:

3.65 ()

Where, Sa is the solvent accessible area of atom a and S a0 is the bare surface area calculated using the set of atomic radii in the quantum chemical calculations. qs is the scaling factor in the range from 0 to unity. k is the Boltzmann constant in J/(molecule K), while V is the volume in Å3/molecule. Here, r denotes the atomic type: C, H, N, O, F, Cl, S, P, Br or I. For mixtures, it takes the following form (the interaction here is considered using all possible element pairs):

3.66 ()

Here,

is the dispersion free energy for the mixture in joules. N is the number of molecules, where and the volume of the ideal mixture given by

3.67 ()

is the mole fraction given by: . Here, i and j sum over all the components bond types, namely and . a refers to the number of components in the mixture. This contribution is used to account for all possible interactions between the atoms in different ensembles. For a binary mixture (Wang et al., 2007):

3.68 ()

This can be reduced as

3.69 ()

Here

3.70 ()

The dispersion term of the residual chemical potential for species a and b is given as

3.71 ()

Hence, we obtain:

3.72 ()

Similarly, for bond type b, we have:

3.73 ()

For large complex molecules, additional corrections are required, which depend on the molecular size. It should be noted that the dispersion term is primarily responsible for the pure component calculations. However, it is the reverse for a mixture, where the contribution to the activity coefficient is negligible because of its cancellation between the solvent and pure solute phases. Hence, we shall neglect this contribution when we compute the activity coefficient.

Now again referring back to Equation 3.7, we will now concentrate on the remaining two terms, that is,

and . These are calculated using the Staverman−Guggenheim contribution (Smith, Van Ness, & Abbott, 2001), where the liquid volumes can be accessed from the DIPPR database (Design Institute for Physical Properties, Sponsored by AIChE DIPPR Project 801). The contribution was derived based on a lattice model where the activity coefficient is caused by the difference in the shape and size of the molecule. It is given as

3.74 ()

where:

3.75 ()

Here, Vi total is the volume of the molecule obtained from the COSMO file output (Figure 3.3). A norm refers to a normalization constant, as Ai total is not dimensionless. For molecules having similar shape and sizes,

, and hence, this term can be neglected. It should also be noted that for pure component, so it is used to model the cavity formation between a solvent and a solute only. We now are in a position to write the total expression for the activity coefficient of a component in the mixture. Thus, from Equations 3.62 through 3.65, we have the following equation:

By writing out each contribution, we have

Thus, the activity coefficient of a segment that is

and sigma potential are connected through the following relation:

3.76 ()

The activity coefficient of segment in the mixture and in the pure liquid,

and , can then be rewritten as

3.77 ()

The final expression then takes the form:

3.78 ()

This equation is referred to as the COSMO Segment Activity Coefficient (COSMO-SAC) model. Now, we shall move ahead to compute the activity coefficient using COSMO-SAC model (Equation 3.78), as these values will be required to predict the tie line composition in case of ternary or binary system. The entire model is characterized by our earlier COSMO-SAC parameters, a eff = 6.32 Å2 (surface area of a standard segment), α = 8419 kcal Å4 mol–1 e–2 (misfit energy constant) for misfit energy interaction, chb = 75006 kcal Å4 mol–1 e–2 (hydrogen bonding energy constant) and σ hb = 0.0084 e Å–2 (hydrogen bonding cutoff).

#### 3.3  Predictions of Tie Lines

The equilibrium for any ternary liquid–liquid system is defined by the equation:

3.79 ()

Here, our aim is to compute the

, the activity coefficient of component I, as per Equation 3.78 from COSMO-SAC model. These values will be different for different phases (I and II). The corresponding mole fractions in phases I and II will be denoted by and , respectively. Assuming an adiabatic operation and that both feed F and solvent S enter a theoretical stage at same temperature, we are then left to account for the heat of mixing. This value for most thermodynamics calculations involving liquid−liquid equilibria (LLE) is sufficiently small, as we observe that a small temperature change occurs. For calculating the LLE split, we adopt the modified Rachford−Rice algorithm (Figure 3.9). Initially, a random feed composition to any known LLE data is used. In any case, we can use the reported data to get an initial feed point. Thus, the feed concentration ( ) is then calculated using the following equation:

3.80 ()

Figure 3.9   Modified Rachford–Rice algorithm for COSMO-SAC model.

Assuming a unity feed rate, that is, F = 1, we visualize the process of one feed flow, namely F, entering a flash column and then flashing into one liquid phase (L) and one vapor phase (V), such that F = L + V. In the next step, the values of distribution coefficient (Ki, I = 1,2,3) are then assigned as follows:

3.81 ()

Where,

and are predicted using the COSMO-SAC model. With the values of , the isothermal flash equation is solved using the conventional flash expression:

3.82 ()

subject to,

3.83 ()

Here, we have assumed L =

and V = , as it contains both liquid phases. Both the liquid phases represent the flow rate of the extract and the raffinate phase . Equation 3.83 is nonlinear in nature. This is first iteratively solved for , which is a fraction and varies between zero and unity. Thereafter, the mole fractions in both phases are calculated via Equations 3.84 and 3.85, respectively.

3.84 ()

3.85 ()

For the ionic liquid systems described in the book, the in-house experimental data were used to assume the feed compositions (Equation 3.80). Thus, in a summary, the COSMO-SAC model can predict the tie lines based on the initial input of

. Here, the effect of pressure on LLE computations has been assumed to be negligible. Thus, not only ionic liquids but also any unknown system can be predicted once its COSMO files (Figure 3.3) are generated. This also makes use of the fact that a successful computation of by Equation 3.82 indicates a heterogeneous region. A feed point outside the bimodal cure, as observed in a ternary LLE diagram, will usually fail to converge Equation 3.82. Thus, based on any random feed composition ( ), one can generate a tie line data. These predictions can be further used for other higher-order, namely quaternary and quinary, systems. It has been described in detail for the computation of VLLE (vapor−liquid−liquid equilibria) diagrams by making necessary adjustments. This has been discussed in Chapter 4.

#### 3.4  Predictions of LLE of Ionic Liquid Systems

Recent experimental work suggests that ILs can be considered as fully dissociated cations and anions. Further, the mechanism of hydrogen bonding (Scheiner, 1997) in the ILs is also not fully available. This fact was used to study the coordination chemistry extraction mechanisms of metal ions, particularly actinide (Cocalia, Gutowski, & Rogers, 2006; Jensen, Neuefeind, Beitz, Skanthakumar, & Soderholm, 2003; MacFarlane & Seddon, 2007). Experimentally, these prove that ILs are completely dissociated into cations and anions when in solution. Also, the definition of IL states it to be a liquid with melting point below 100°C and containing ions exhibiting ionic conductivity. The fact that the ionic conductivity can be measured experimentally proves that ILs consist of ions in solutions. Ionic conductivity studies have been carried out by Carda-Broch et al., (2003), Berthod and Armstrong (2003) on fluorinated anions such as [PF6].

Further, in molecular dynamic simulations, a cubical box is taken, which contains equal number of cations and anions, instead the whole IL molecule. For example, in the work of Morrow and Maginn (2002), the authors studied the thermo-physical properties of ionic liquids by using 300 cations and 300 anions. Similarly, in the work of Lee, Jung and Han (2005), the ionic conductivity of ILs were found using a cubical box of 100 cations and 100 anions. On the contrary, Earle et al. (2006) recently distilled IL without decomposition, at very low pressures. In their pioneering work, they were successful in distilling the [(CF3SO2)2N]-based ILs at 100 Pa and 573 K, without decomposition. However, in LLE experiments, ILs are assumed to be in isobaric condition of 1 atm, whereas the distillation is observed with pressures ~0.000009 atm. Since the operating pressures of LLE are 1 atm, our approach holds true.

Dissociation was also observed during vapor−liquid equilibria (VLE) predictions on reported systems (Kato & Gmehling, 2005; Kato, Krummen, & Gmehling, 2004) in ILs, where we obtained better results than obtained using a single composite molecule (Diedenhofen, Eckert, & Klamt, 2003; Table 3.1). This prompted us to propose that the IL be considered as a dissociated pair of cation and anion for both the VLE and LLE predictions. This novel approach can be applied to ILs containing a limitless combination of cation and anion (~1018). Therefore, one has to do separate quantum mechanical calculations for cations and anions. Thereafter, one can test an IL for a particular application by addition of sigma profiles and sigma potential of cation and anion only.

### Table 3.1   VLE Comparison of RMSD with Different Models

COSMO-RS

No.

Reference

System

Single Molecule

Cation and Anion

1

[Emim][(CF3SO2)2N] + Acetone

2.70

1.81

2

[Emim][(CF3SO2)2N] + 2-Propanol

3.80

2.60

3

[Emim][(CF3SO2)2N] + Water

9.80

5.26

4

[Bmim][(CF3SO2)2N] + Acetone

3.50

3.08

5

[Bmim][(CF3SO2)2N] + 2-Propanol

4.00

2.89

6

[Bmim][(CF3SO2)2N] + Water

9.80

4.56

7

[Mmim][(CH3)2PO4] + Acetone

5.10

2.50

8

[Mmim][(CH3)2PO4] + Tetrahydrofuran

4.60

2.28

9

[Mmim][(CH3)2PO4] + Water

10.60

6.84

10

[Mmim][(CF3SO2)2N] + Benzene

8.92

3.66

11

[Mmim][(CF3SO2)2N] + Cyclohexane

6.26

4.23

12

[Emim][EtSO4] + Benzene

9.75

5.62

13

[Emim][EtSO4] + Cyclohexane

6.94

2.36

Average

6.59

3.66

Moreover, till date, an application concerned with the assumption of complete dissociation of ILs into cations and anions with equimolar concentrations has been carried out earlier by Klamt and Schuurmann (1993) has assumed a complete dissociation of cations and anions for their prediction of infinite dilution activity coefficients (IDAC) of solutes in ILs by COSMO-RS implementation. In their work, the ILs have been described by an equimolar mixture of two distinct ions, the cation and the anion, which finally contribute to the sigma profile of the mixture as two different compounds.

Hence, for the prediction of LLE-based system involving ILs, two modes of approach were assumed: (1) IL as a pure solvent and (2) IL consisting of a pair of cations and anions (Table 3.2). The latter approach assumes a solvent of a mixture of two components: cation and anion. The latter approach has proved to give better phase split with lower root-mean-square deviation. Hence, we proposed a complete dissociation of ILs into cations and anions with equimolar quantities, which is in line with the linear addition of the sigma profiles of the cation and anion.

3.86 ()

Where,

and are the sigma profiles for cation and anion, respectively. This is based on the assumption that the IL in a solution is fully dissociated into its respective cations and anions with equimolar concentrations. It should be noted that the entropy of mixing will not exist, since we are dealing with a sigma profile of a single compound/mixture. Thereafter, the linear additions of COSMO area and volume are done, so as to obtain the profile of single mixture/compound.

### Table 3.2   LLE Comparison of Single and Cation–Anion Pair by Using COSMO-RS

Root-Mean-Square Deviation

Experiment

System No.

T/K

System

Single Molecule

Cation + Anion

Error

1

298.15

[Bmim][(CF3SO3] – Ethanol-ethyl-tert-butyl ether

NS

0.095

NA

2

298.15

[Bmim][(CF3SO3] – Ethanol-tert-amyl ethyl ether

NS

0.032

0.004

3

298.15

[Hmim][BF4] – Benzene-heptane

0.711

0.043

0.004

4

298.15

[Hmim][BF4] – Benzene-dodecane

0.087

0.098

0.004

5

298.15

0.481

0.072

0.004

6

298.15

[Hmim][BF4] – Ethanol-hexene

NS

0.147

0.004

7

298.15

[Hmim][BF4] – Ethanol-heptene

NS

0.125

0.004

8

298.15

[Hmim][PF6] – Benzene-heptane

0.856

0.042

0.004

9

298.15

[Hmim][PF6] – Benzene-dodecane

0.385

0.047

0.004

10

298.15

0.457

0.059

0.004

11

298.15

[Hmim][PF6] – Ethanol-hexene

NS

0.315

0.004

12

298.15

[Hmim][PF6] – Ethanol-heptene

NS

0.306

0.004

13

298.15

0.029

0.005

NA

14

298.15

0.03

0.008

NA

15

298.15

[Omim][Cl] – Ethanol-tert-amyl ethyl ether

0.382

0.067

NA

16

298.15

[Omim][Cl] – Benzene-heptane

0.896

0.096

0.006

17

298.15

[Omim][Cl] – Benzene-dodecane

0.058

0.096

0.006

18

298.15

0.146

0.07

0.006

19

298.15

[Emim][C8H17SO4] – Benzene-heptane

0.362

0.07

0.006

20

298.15

0.456

0.14

0.006

21

298.15

[Omim][MDEG] a – Benzene-heptane

NS

0.103

0.006

22

298.15

NS

0.043

0.006

23

313.15

[Bmpy][BF4] – Xylene-octane

NS

0.03

0.0025

24

348.15

[Bmpy][BF4] – Xylene-Octane

NS

0.024

0.0025

25

313.15

[Bmpy][BF4] – Ethylbenzene-octane

NS

0.041

0.0025

26

348.15

[Bmpy][BF4] – Ethylbenzene-octane

NS

0.029

0.0025

27

313.15

[Bmpy][BF4] – Benzene-hexane

NS

0.048

0.0025

28

333.15

[Bmpy][BF4] – Benzene-hexane

NS

0.036

0.0025

29

298.15

[Bmpy][BF4] – Toluene-heptane

NS

0.075

0.0025

30

298.15

[Mmim][CH3SO4] – Toluene-heptane

0.455

0.008

0.0025

31

298.15

[Emim][C2H5SO4] – Toluene-heptane

0.116

0.031

0.0025

32

298.15

[Bmim][[CH3SO4] – Toluene-heptane

0.235

0.032

0.0025

33

313.15

[Emim][EtSO4] – Ethanol-hexene

NS

0.124

0.005

34

313.15

[Emim][EtSO4] – Ethanol-heptene

NS

0.023

0.005

35

313.15

[E-2,3-dmim][EtSO4] – Ethanol-hexene

NS

0.032

0.005

36

313.15

[E-2,3-dmim][EtSO4] – Ethanol-heptene

NS

0.065

0.005

Note: NS: No splitting of the phase; NA: Not available

Notes:

a   MDEG: monomethyldiethyleneglycol

The results of all the reported IL ternary systems are given in Table 3.2. Using the sigma profile for the single composite molecule, it is observed that for nearly half of the systems, the values of the activity coefficients were not able to predict the split between the extract (IL phase) and raffinate phase for the given tie line. Predictions giving an average RMSD of 36.5% and a maximum root-mean-square deviation (RMSD) of 90% ([Omim][Cl]-Benzene-Heptane) were observed. NO SPLIT (NS) condition for 15 systems was encountered. For 13 systems, an RMSD much greater than 10% was observed.

The data sets contain 12 different ILs comprising seven different cations and eight anions. Apart from imidazolium-based ILs, pyridinium-based IL (1-butyl-3-methylpyridinium tetrafluoroborate [Bmpy][BF4]) has also been studied. Switching from composite to additive sigma profile, we notice that a ‘NO SPLIT’ situation has been converted to a ‘SPLIT’ between the two phases for all the 15 cases. It is clear that the improvement is drastic for the IL containing system namely [Bmim][CF3SO3], [Omim][MDeg] and [Bmpy][BF4] systems and also for the alkene-based systems [System No. 6, 7, 11, 12], for which no prediction had been possible when considering the IL as a single molecule. For the [Bmim][CF3SO3], [Omim][MDeg] and [Bmpy][BF4] systems (System No. 1, 2, 21, 22 and 23–29), the RMSD is less than 10%. Out of the 36 systems studied, most of the predictions are excellent, considering the fact that only six systems gave RMSD greater than 10%.

#### Acknowledgments

A part of this chapter has been adapted from my earlier National Program for Technology Enhanced Learning (NPTEL) web course entitled ‘Molecular Simulation in Chemical Engineering’. I am deeply indebted to the Ministry of Human Resource and Development, Government of India, for providing me financial support for the same. Further, I would like to thank Prof. Stanley I. Sandler, University of Delaware, for helping me in learning the basics of COSMO-based modeling. I am also indebted to his student Dr. Russel I. Burnett for helping me out in the COSMO-SAC equations. This visit would not have been possible without the funding from IUSSTF (Indo–U.S. Science and Technology Forum).

#### Canonical Partition Functions

The canonical partition function that is Q(N, V, β) is given as

A.1 ()

The probability of occurrence of a particular miscrostate i with Energy E α is given as

A.2 ()

The significance of β can be made from the fact that as β > 0, the state of higher energy is less likely as the state of lower energy. We now shall see the process of evaluating the thermodynamic properties based on canonical partition functions. It should be noted that the probability of finding a system in any microstate of energy E α would be the product of the probability that the system is in a particular microstate with energy E α and the number of states having that energy (also known as degeneracy). Thus, we have:

A.3 ()

We now define the probability of occurrence of a particular microstate i with energy E α as

A.4 ()

In a similar manner, the probability of occurrence of the energy level E α, p(E α) is

A.5 ()

Here Degeneracy = ω(E α)

Thus, the partition function can be written in terms of either levels or states, as below:

A.6 ()

Now, once the partition function has been defined, we can derive thermodynamics functions such as internal energy U:

A.7 ()

However, from the definition, we have:

A.8 ()

Thus, the internal energy U takes the form:

A.9 ()

#### References

Banerjee, T. , Singh, M. K. , & Khanna, A. (2006). Prediction of binary VLE for imidazolium based ionic liquid systems using COSMO-RS. Industrial & Engineering Chemistry Research, 45(9), 3207–3219. doi:10.1021/ie051116c.
Burnett, R. I. (2012). Predicting liquid-phase thermodynamic properties using COSMO-SAC. PhD Thesis, University of Delaware, Newark, DE.
Carda-Broch, S. , Berthod, A. , & Armstrong, D. W. (2003). Solvent properties of the 1-butyl-3-methylimidazolium hexafluorophosphate ionic liquid. Analytical and Bioanalytical Chemistry, 375(2), 191–199. doi:10.1007/s00216-002-1684-1.
Cocalia, V. A. , Gutowski, K. E. , & Rogers, R. D. (2006). The coordination chemistry of actinides in ionic liquids: A review of experiment and simulation. Coordination Chemistry Reviews, 250(7–8), 755–764. doi:10.1016/j.ccr.2005.09.019.
Diedenhofen, M. , Eckert, F. , & Klamt, A. (2003). Prediction of infinite dilution activity coefficients of organic compounds in ionic liquids using COSMO-RS. Journal of Chemical & Engineering Data, 48(3), 475–479. doi:10.1021/je025626e.
Döker, M. , & Gmehling, J. (2005). Measurement and prediction of vapor–liquid equilibria of ternary systems containing ionic liquids. Fluid Phase Equilibria, 227(2), 255–266. doi:10.1016/j.fluid.2004.11.010.
Earle, M. J. , Esperanca, J. M. S. S. , Gilea, M. A. , Canongia Lopes, J. N. , Rebelo, L. P. N. , Magee, J. W. , … Widegren, J. A. (2006). The distillation and volatility of ionic liquids. Nature, 439(7078), 831–834. doi:10.1038/nature04451.
Grensemann, H. , & Gmehling, J. (2005). Performance of a conductor-like screening model for real solvents model in comparison to classical group contribution methods. Industrial & Engineering Chemistry Research, 44(5), 1610–1624. doi:10.1021/ie049139z.
Jensen, M. P. , Neuefeind, J. , Beitz, J. V. , Skanthakumar, S. , & Soderholm, L. (2003). Mechanisms of metal ion transfer into room-temperature ionic liquids: The role of anion exchange. Journal of the American Chemical Society, 125(50), 15466–15473. doi:10.1021/ja037577b.
Kato, R. , & Gmehling, J. (2005). Measurement and correlation of vapor–liquid equilibria of binary systems containing the ionic liquids [EMIM][(CF3SO2)2N], [BMIM][(CF3SO2)2N], [MMIM][(CH3)2PO4] and oxygenated organic compounds respectively water. Fluid Phase Equilibria, 231(1), 38–43. doi:10.1016/​j.fluid.2005.01.002.
Kato, R. , Krummen, M. , & Gmehling, J. (2004). Measurement and correlation of vapor–liquid equilibria and excess enthalpies of binary systems containing ionic liquids and hydrocarbons. Fluid Phase Equilibria, 224(1), 47–54. doi:10.1016/​j.fluid.2004.05.009.
Klamt, A. (1995). Conductor-like screening model for real solvents: A new approach to the quantitative calculation of solvation phenomena. The Journal of Physical Chemistry, 99(7), 2224–2235. doi:10.1021/j100007a062.
Klamt, A. , Jonas, V. , Bürger, T. , & Lohrenz, J. C. W. (1998). Refinement and parametrization of COSMO-RS. The Journal of Physical Chemistry A, 102(26), 5074–5085. doi:10.1021/jp980017s.
Klamt, A. , & Schuurmann, G. (1993). COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. Journal of the Chemical Society, Perkin Transactions 2(5), 799–805. doi:10.1039/P29930000799.
Klamt, A. (2005). COSMO-RS: From quantum chemistry to fluid phase thermodynamics and drug design. Amsterdam, the Netherlands: Elsevier.
Lee, S. U. , Jung, J. , & Han, Y.-K. (2005). Molecular dynamics study of the ionic conductivity of 1-n-butyl-3-methylimidazolium salts as ionic liquids. Chemical Physics Letters, 406(4–6), 332–340. doi:10.1016/j.cplett.2005.02.109.
Lin, S.-T. , Chang, J. , Wang, S. , Goddard, W. A. , & Sandler, S. I. (2004). Prediction of vapor pressures and enthalpies of vaporization using a COSMO solvation model. The Journal of Physical Chemistry A, 108(36), 7429–7439. doi:10.1021/jp048813n.
Lin, S.-T. , & Sandler, S. I. (2002). A priori phase equilibrium prediction from a segment contribution solvation model. Industrial & Engineering Chemistry Research, 41(5), 899–913. doi:10.1021/ie001047w.
MacFarlane, D. R. , & Seddon, K. R. (2007). Ionic liquids–Progress on the fundamental issues. Australian Journal of Chemistry, 60(1), 3–5. doi:10.1071/CH06478.
Morrow, T. I. , & Maginn, E. J. (2002). Molecular dynamics study of the ionic liquid 1-n-Butyl-3-methylimidazolium hexafluorophosphate. The Journal of Physical Chemistry B, 106(49), 12807–12813. doi:10.1021/jp0267003.
Pye, C. C. , Ziegler, T. , van Lenthe, E. , & Louwen, J. N. (2009). An implementation of the conductor-like screening model of solvation within the Amsterdam density functional package–Part II. COSMO for real solvents. Canadian Journal of Chemistry, 87(7), 790–797. doi:10.1139/V09-008.
Scheiner, S. (1997). Hydrogen bonding: A theoretical perspective. New York: Oxford University Press.
Smith, J. M. , Van Ness, H. C. , Abbott, M. M. (2001). Introduction to chemical engineering Thermodynamics (6th ed.). Boston, MA: McGraw-Hill.
Wang, S. , Lin, S.-T. , Chang, J. , Goddard, W. A. , & Sandler, S. I. (2006). Application of the COSMO−SAC−BP solvation model to predictions of normal boiling temperatures for environmentally significant substances. Industrial & Engineering Chemistry Research, 45(16), 5426–5434. doi:10.1021/ie050352k.
Wang, S. , Sandler, S. I. , & Chen, C.-C. (2007). Refinement of COSMO−SAC and the Applications. Industrial & Engineering Chemistry Research, 46(22), 7275–7288. doi:10.1021/ie070465z.
Wang, S. , Stubbs, J. M. , Siepmann, J. I. , & Sandler, S. I. (2005). Effects of conformational distributions on sigma profiles in COSMO theories. The Journal of Physical Chemistry A, 109(49), 11285–11294. doi:10.1021/jp053859h.

Sections 3.23.4 reprinted (adapted) from T. Banerjee, K. K. Verma and A. Khanna, Liquid–liquid equilibria for ionic liquid-based systems using COSMO-RS: Effect of cation and anion combination. AIChEJ 54, 1874–1885, 2008. Copyright Wiley-VCH Verlag GmbH & Co. KGaA. Reproduced with permission.

## Use of cookies on this website

We are using cookies to provide statistics that help us give you the best experience of our site. You can find out more in our Privacy Policy. By continuing to use the site you are agreeing to our use of cookies.