The prediction of properties of liquidphase 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.
The prediction of properties of liquidphase 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.
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 nontranslational contribution of the Gibb’s free energy of solvation. This is related by Equations 3.2 and 3.3, as given below:
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.
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 COSMOSAC (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 COSMOSAC 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.
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):
By rearranging the terms, we obtain:
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 twoliquid (NRTL) and UNI QUAsiChemical (UNIQUAC) models. Hence, the final expression takes the form (Equation 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:
Within the COSMOSAC 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 Waallike 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 ensemblebased 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.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/SemiEmpirical) 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.
The segment surface is created by a vdWlike surface, which is determined by the union of spheres drawn around each of the atoms in the geometryoptimized 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}.
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.
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 selfenergy 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
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 selfenergies; this is given by
Here,
is the interaction energy between segments, while is the selfenergy 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):
By reorganizing the terms, we obtain:
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:
The energy of solvation (E ^{solv}) is then minimized with respect to the segments’ charges (q _{α}), as given below:
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 ‘selfconsistent field’ approach. This process of generating and identifying the set of ρ( r ) and q that are selfconsistent is called a ‘conductorlike 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 Gaussian03generated 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 area (Å^{2}) = 371.79 and volume (Å^{3}) = 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.
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 selfconsistent 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 atomcentered 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.
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:
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.Based on statistical thermodynamics, COSMOSAC 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 n_{i} segments. Hence, the total number of segments is
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):
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
Thus, the probability of obtaining a segment pair
will be
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
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:
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:
From statistical thermodynamics, we know:
Hence, Equation 3.25 takes the form:
Here,
are the chemical potentials for the two segments α and β, respectively. Hence, we have:
Equation 3.24 then takes the form:
The nontranslational part of the chemical potential is given as
This implies that the second part is the translational contribution for chemical potential. So, for ideal gases, the nontranslational part is usually zero, or
The final form of Equation 3.29 takes the form:
For obtaining the average number of ensemble pairs, we thus use the relation as below:
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:
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.
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:
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:
By combining the above two expressions, we obtain:
Rearranging the equation provides the following form:
This can be rewritten as
where:
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:
Here, a _{eff} is the area of the effective segment and is an adjustable parameter.
is charge density of segments. Here, c_{es} is a constant and is equal to . The equation above reflects the selfconsistency 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 _{2} … n_{n} . 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 (mol_{i}) 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):
So, by summing up the contributions from all the segments, we obtain:
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:
However, Equation 3.44 shall also be valid when the solvent itself is a perfect conductor. This will be applicable for the following condition:
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:
Thus, we have to subtract the term
from the chemical potential to obtain the nontranslational contribution.
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 threedimensional segment information. Further, no attempt has been made to consider the steric hindrance or the longrange 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 COSMOSAC model, the averaging algorithm is proposed as below:
where:
The averaging of the screening charge densities, therefore, creates a new contribution to the overall Gibb’s free energy. This is provided as below:
Here, the third term on the righthand side is the Gibb’s free energy change occurring due to the averaging process.
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
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.
The threedimensional 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:
where:
For a mixture, the σprofile is determined from the areaweighted average of contributions from all the components present in the mixture. This is given as
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:
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 E_{mf} , which is the energy computed due to the difference in shape and size between a pair of segments. Another contribution, the hydrogen bonding energy E_{hb} , 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.
Where, α’ is a constant for the misfit energy and is calculated from the surface area of a standard segment.
Here, c_{hb} 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 nonzero 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:
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.
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, 1alkyl3methylimidazolium ([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: 1ethyl3methylimidazolium; BMIM: 1butyl3methylimidazolium; OMIM: 1Octyl3methylimidazolium).
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 hydrogenacceptor bonding that dominates the lefthand 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 [PF_{6}], 0.014 e/Å^{2} for [CF_{3}SO_{3}], 0.015 e/Å^{2} for [EtSO_{4}] and 0.011 e/Å^{2} for [BF_{4}], 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.
Finally, we obtain the restoring free energy of molecule in solvent or mixture as
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:
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 meanfree approach, the perturbation approach. For pure component, the expression takes the form:
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:
Where, S_{a} 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. q_{s} 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):
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
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):
This can be reduced as
Here
The dispersion term of the residual chemical potential for species a and b is given as
Hence, we obtain:
Similarly, for bond type b, we have:
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
where:
Here, V_{i} ^{total} is the volume of the molecule obtained from the COSMO file output (Figure 3.3). A _{norm} refers to a normalization constant, as A_{i} ^{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:
The activity coefficient of segment in the mixture and in the pure liquid,
and , can then be rewritten as
The final expression then takes the form:
This equation is referred to as the COSMO Segment Activity Coefficient (COSMOSAC) model. Now, we shall move ahead to compute the activity coefficient using COSMOSAC 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 COSMOSAC 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, c_{hb} = 75006 kcal Å^{4} mol^{–1} e^{–2} (hydrogen bonding energy constant) and σ _{hb} = 0.0084 e Å^{–2} (hydrogen bonding cutoff).
The equilibrium for any ternary liquid–liquid system is defined by the equation:
Here, our aim is to compute the
, the activity coefficient of component I, as per Equation 3.78 from COSMOSAC 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:
Figure 3.9 Modified Rachford–Rice algorithm for COSMOSAC 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 (K_{i}, I = 1,2,3) are then assigned as follows:
Where,
and are predicted using the COSMOSAC model. With the values of , the isothermal flash equation is solved using the conventional flash expression:
subject to,
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.
For the ionic liquid systems described in the book, the inhouse experimental data were used to assume the feed compositions (Equation 3.80). Thus, in a summary, the COSMOSAC 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 higherorder, 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.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 CardaBroch et al., (2003), Berthod and Armstrong (2003) on fluorinated anions such as [PF_{6}].
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 thermophysical 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 [(CF_{3}SO_{2})_{2}N]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 (~10^{18}). 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.
COSMORS 


No. 
Reference 
System 
Single Molecule 
Cation and Anion 
1 
[Emim][(CF_{3}SO_{2})_{2}N] + Acetone 
2.70 
1.81 

2 
[Emim][(CF_{3}SO_{2})_{2}N] + 2Propanol 
3.80 
2.60 

3 
[Emim][(CF_{3}SO_{2})_{2}N] + Water 
9.80 
5.26 

4 
[Bmim][(CF_{3}SO_{2})_{2}N] + Acetone 
3.50 
3.08 

5 
[Bmim][(CF_{3}SO_{2})_{2}N] + 2Propanol 
4.00 
2.89 

6 
[Bmim][(CF_{3}SO_{2})_{2}N] + Water 
9.80 
4.56 

7 
[Mmim][(CH_{3})_{2}PO_{4}] + Acetone 
5.10 
2.50 

8 
[Mmim][(CH_{3})_{2}PO_{4}] + Tetrahydrofuran 
4.60 
2.28 

9 
[Mmim][(CH_{3})_{2}PO_{4}] + Water 
10.60 
6.84 

10 
[Mmim][(CF_{3}SO_{2})_{2}N] + Benzene 
8.92 
3.66 

11 
[Mmim][(CF_{3}SO_{2})_{2}N] + Cyclohexane 
6.26 
4.23 

12 
[Emim][EtSO_{4}] + Benzene 
9.75 
5.62 

13 
[Emim][EtSO_{4}] + 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 COSMORS 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 LLEbased 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 rootmeansquare 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.
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.
RootMeanSquare Deviation 
Experiment 


System No. 
T/K 
System 
Single Molecule 
Cation + Anion 
Error 
1 
298.15 
[Bmim][(CF_{3}SO_{3}] – Ethanolethyltertbutyl ether 
NS 
0.095 
NA 
2 
298.15 
[Bmim][(CF_{3}SO_{3}] – Ethanoltertamyl ethyl ether 
NS 
0.032 
0.004 
3 
298.15 
[Hmim][BF_{4}] – Benzeneheptane 
0.711 
0.043 
0.004 
4 
298.15 
[Hmim][BF_{4}] – Benzenedodecane 
0.087 
0.098 
0.004 
5 
298.15 
[Hmim][BF_{4}] – Benzenehexadecane 
0.481 
0.072 
0.004 
6 
298.15 
[Hmim][BF_{4}] – Ethanolhexene 
NS 
0.147 
0.004 
7 
298.15 
[Hmim][BF_{4}] – Ethanolheptene 
NS 
0.125 
0.004 
8 
298.15 
[Hmim][PF_{6}] – Benzeneheptane 
0.856 
0.042 
0.004 
9 
298.15 
[Hmim][PF_{6}] – Benzenedodecane 
0.385 
0.047 
0.004 
10 
298.15 
[Hmim][PF_{6}] – Benzenehexadecane 
0.457 
0.059 
0.004 
11 
298.15 
[Hmim][PF_{6}] – Ethanolhexene 
NS 
0.315 
0.004 
12 
298.15 
[Hmim][PF_{6}] – Ethanolheptene 
NS 
0.306 
0.004 
13 
298.15 
[Omim][Cl] – Methanolhexadecane 
0.029 
0.005 
NA 
14 
298.15 
[Omim][Cl] – Ethanolhexadecane 
0.03 
0.008 
NA 
15 
298.15 
[Omim][Cl] – Ethanoltertamyl ethyl ether 
0.382 
0.067 
NA 
16 
298.15 
[Omim][Cl] – Benzeneheptane 
0.896 
0.096 
0.006 
17 
298.15 
[Omim][Cl] – Benzenedodecane 
0.058 
0.096 
0.006 
18 
298.15 
[Omim][Cl] – Benzenehexadecane 
0.146 
0.07 
0.006 
19 
298.15 
[Emim][C_{8}H_{17}SO_{4}] – Benzeneheptane 
0.362 
0.07 
0.006 
20 
298.15 
[Emim][C_{8}H_{17}SO_{4}] – Benzenehexadecane 
0.456 
0.14 
0.006 
21 
298.15 
[Omim][MDEG] ^{a} – Benzeneheptane 
NS 
0.103 
0.006 
22 
298.15 
[Omim][MDEG] ^{a} – Benzenehexadecane 
NS 
0.043 
0.006 
23 
313.15 
[Bmpy][BF_{4}] – Xyleneoctane 
NS 
0.03 
0.0025 
24 
348.15 
[Bmpy][BF_{4}] – XyleneOctane 
NS 
0.024 
0.0025 
25 
313.15 
[Bmpy][BF_{4}] – Ethylbenzeneoctane 
NS 
0.041 
0.0025 
26 
348.15 
[Bmpy][BF_{4}] – Ethylbenzeneoctane 
NS 
0.029 
0.0025 
27 
313.15 
[Bmpy][BF_{4}] – Benzenehexane 
NS 
0.048 
0.0025 
28 
333.15 
[Bmpy][BF_{4}] – Benzenehexane 
NS 
0.036 
0.0025 
29 
298.15 
[Bmpy][BF_{4}] – Tolueneheptane 
NS 
0.075 
0.0025 
30 
298.15 
[Mmim][CH_{3}SO_{4}] – Tolueneheptane 
0.455 
0.008 
0.0025 
31 
298.15 
[Emim][C_{2}H_{5}SO_{4}] – Tolueneheptane 
0.116 
0.031 
0.0025 
32 
298.15 
[Bmim][[CH_{3}SO_{4}] – Tolueneheptane 
0.235 
0.032 
0.0025 
33 
313.15 
[Emim][EtSO_{4}] – Ethanolhexene 
NS 
0.124 
0.005 
34 
313.15 
[Emim][EtSO_{4}] – Ethanolheptene 
NS 
0.023 
0.005 
35 
313.15 
[E2,3dmim][EtSO_{4}] – Ethanolhexene 
NS 
0.032 
0.005 
36 
313.15 
[E2,3dmim][EtSO_{4}] – Ethanolheptene 
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 rootmeansquare deviation (RMSD) of 90% ([Omim][Cl]BenzeneHeptane) 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 imidazoliumbased ILs, pyridiniumbased IL (1butyl3methylpyridinium tetrafluoroborate [Bmpy][BF_{4}]) 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][CF_{3}SO_{3}], [Omim][MDeg] and [Bmpy][BF_{4}] systems and also for the alkenebased 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][CF_{3}SO_{3}], [Omim][MDeg] and [Bmpy][BF_{4}] 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%.
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 COSMObased modeling. I am also indebted to his student Dr. Russel I. Burnett for helping me out in the COSMOSAC equations. This visit would not have been possible without the funding from IUSSTF (Indo–U.S. Science and Technology Forum).
The canonical partition function that is Q(N, V, β) is given as
The probability of occurrence of a particular miscrostate i with Energy E _{α} is given as
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:
We now define the probability of occurrence of a particular microstate i with energy E _{α} as
In a similar manner, the probability of occurrence of the energy level E _{α}, p(E _{α}) is
Here Degeneracy = ω(E _{α})
Thus, the partition function can be written in terms of either levels or states, as below:
Now, once the partition function has been defined, we can derive thermodynamics functions such as internal energy U:
However, from the definition, we have:
Thus, the internal energy U takes the form:
Sections 3.2–3.4 reprinted (adapted) from T. Banerjee, K. K. Verma and A. Khanna, Liquid–liquid equilibria for ionic liquidbased systems using COSMORS: Effect of cation and anion combination. AIChEJ 54, 1874–1885, 2008. Copyright WileyVCH Verlag GmbH & Co. KGaA. Reproduced with permission.