Solubility of Nonelectrolytes: A First-Principles Computational

Apr 28, 2014 - Using a combination of classical molecular dynamics and symmetry adapted intermolecular perturbation theory, we develop a high-accuracy...
0 downloads 3 Views 587KB Size
Article pubs.acs.org/JPCB

Solubility of Nonelectrolytes: A First-Principles Computational Approach Nicholas E. Jackson,† Lin X. Chen,†,‡ and Mark A. Ratner*,† †

Department of Chemistry, Northwestern University, Evanston, Illinois 60208, United States Chemical Sciences and Engineering Division, Argonne National Laboratory, Lemont, Illinois 60439, United States



S Supporting Information *

ABSTRACT: Using a combination of classical molecular dynamics and symmetry adapted intermolecular perturbation theory, we develop a high-accuracy computational method for examining the solubility energetics of nonelectrolytes. This approach is used to accurately compute the cohesive energy density and Hildebrand solubility parameters of 26 molecular liquids. The energy decomposition of symmetry adapted perturbation theory is then utilized to develop multicomponent Hansen-like solubility parameters. These parameters are shown to reproduce the solvent categorizations (nonpolar, polar aprotic, or polar protic) of all molecular liquids studied while lending quantitative rigor to these qualitative categorizations via the introduction of simple, easily computable parameters. Notably, we find that by monitoring the first-order exchange energy contribution to the total interaction energy, one can rigorously determine the hydrogen bonding character of a molecular liquid. Finally, this method is applied to compute explicitly the Flory interaction parameter and the free energy of mixing for two different small molecule mixtures, reproducing the known miscibilities. This methodology represents an important step toward the prediction of molecular solubility from first principles.



which are parametrized via first-principles calculations, have proven successful in predicting solubilities1,20 and free energies of solvation,21,22 but they are also highly dependent on appropriate parametrization. The use of quantum-chemical calculations combined with statistical mechanical approaches possesses great potential for predicting solvation energies and solubilities from first-principles with no dependence upon empirical parametrization. While many methods exist that incorporate quantum-chemical detail,23,24 arguably the most successful of these theories is the Conductor-like Screening Model for Real Solvents (COSMO-RS).25 The COSMO-RS approach has shown great success in determining many parameters relevant to solubility and has become a viable candidate for ab initio solubility determination in a variety of fields.26,27 While the introduction of COSMO-RS represents a huge advance in the field of first-principles solubility determination, the COSMO-RS computation of the interaction energy between two surface patches using sigma profiles does not represent a rigorous solution of the Schrodinger equation, and its accuracy, though typically very high, is also dependent upon extensive parametrization for the radii of cavity construction in the COSMO-RS calculation.28 One particularly interesting and intuitively satisfying approach to the solubility of nonelectrolytes is that of solubility parameters. Originally introduced by Hildebrand7 and modified

INTRODUCTION The solubility of nonionic organic molecules is a complex problem that spans nearly every chemical discipline. While many approaches exist for determining solubility, accurate predictions are still a difficult task; solubility commonly troubles many important fields, including pharmaceuticals development,1−3 protein solubility,4 and the processing of organic semiconductors.5,6 Over the years, the factors that influence molecular solubility have been better understood, with Hildebrand,7 Flory,8 Hansen,9 and others10 proposing quantitative means for assessing the miscibility of two molecular species from rigorously defined thermodynamic principles, with early approaches relying upon regular-solution theories 11 and geometric mean approximations to the interaction energies. In the past 40 years, a variety of new methods for determining the solubilities of nonelectrolytes, particularly in water, have been developed, relying upon statistical regressions against experimental functions,12−14 group contribution UNIQUAC methods,15 and quantitative structure−property relationships,16−18 among others.19 While these methods have proven successful in predicting the solvation free energies and solubilities of a variety of compounds, these predictions are, in many cases, highly dependent upon empirical parametrization and consequently exhibit deficiencies when new compounds are introduced that do not belong to the parametrization groups. To this end, atomistic simulations that are inherently more expensive but capable of describing interactions from first-principles are highly desirable. Molecular dynamics simulations, many of © 2014 American Chemical Society

Received: March 10, 2014 Revised: April 18, 2014 Published: April 28, 2014 5194

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

by Hansen,9 this approach recognizes that a molecule’s selfaffinity, or cohesive attraction to its neighbors, is fundamental to its solubility. Hildebrand formalized the cohesive contribution as a “solubility parameter”,29 which can be quantified as the square root of the cohesive energy density. The cohesive energy density used to determine this solubility parameter, δH, is derived from the molar heat of vaporization (ΔHvap) minus a thermal factor (RT), divided by the molar volume (Vmol) (eq 1). δH =

noncovalent interaction types must be captured to have a successful solubility theory. A variety of modifications to the Hildebrand approach have been developed to incorporate different noncovalent contributions into the molecular cohesion energy. The most common alternative, particularly in the polymer community, was introduced by Hansen.9 In Hansen theory, the sum of the squares of three distinct parameters (polar, δP; dispersion, δD; and hydrogen bonding, δHB) is required to equal the original Hildebrand parameter squared (eq 4).

ΔH vap − RT Vmol

2 δ H2 = δ P2 + δ D2 + δ HB

(1)

This approach accommodates different types of noncovalent interactions while maintaining the intuitive nature of Hildebrand’s approach. The likelihood of two molecules being miscible in the Hansen approach is similar to that in the Hildebrand approach. Instead of comparing single solubility parameters, all contributions are considered by comparing socalled “R” values in Hansen theory (eq 5), relative to a given reference,

The formal theory of solubility in this approach was furthered when it was recognized that one could utilize measured solubility parameters to quantify the free energy of mixing using the regular solution theory formula of Flory and Huggins30 (eq 2), ΔGm = NAkBT (n1 ln ϕ1 + n2 ln ϕ2 + n1ϕ2χ12 )

(2)

where ΔGm is the free energy of mixing, NA is Avogadro’s number, kB is Boltzmann’s constant, T is the temperature, n1 is the number of moles of component 1, n2 is the number of moles of component 2, ϕ1 is the volume fraction of component 1, ϕ2 is the volume fraction of component 2, and χ12 is the Flory interaction parameter, indicating the interaction enthalpy between component 1 and component 2. In formulating eq 2, it is typical to use the geometrical mean approximation to write χ12 as a function of the solubility parameters of the two components (eq 3), χ12 = V1

δ12 + δ22 − 2δ1δ2 NAkBT

(4)

R ∝ (δ1P − δ2P)2 + (δ1D − δ2D)2 + (δ1HB − δ2HB)2

(5)

where a small R value is deemed a good figure-of-merit for two molecules being miscible (we note that it is possible to incorporate Hansen parameters into the original Flory− Huggins theory32). This leads to the conclusion that molecules possessing similar solubility parameters, and thus similar types and magnitudes of noncovalent interactions, will be miscible. Yet, while Hansen parameters are a common and useful method for determining solubility,5,33,34 the tuning of these parameters is highly empirical, and little framework exists for determining these parameters from first principles. It is to be emphasized that while these parametric approaches can be highly effective, there is still no knowledge of the explicit interaction between the two molecular species; it is inferred via the geometric mean of the individual parameters. Given these deficiencies, the utility of a rigorously based first-principles version of these parameter methods could have high impact in the field of solubility determination. A field of chemistry that does not often intersect with solubility parameter approaches is that of ab initio quantum chemistry. While the COSMO-RS model represents the standard for computing accurate measures of solubility (activities, partition coefficients, free energies of solvation), COSMO-RS lacks the simple, qualitative understanding present in the solubility parameters approach that is intuitive to every chemist’s understanding of “like dissolves like”. The benefit of the solubility parameter method is the identification of each molecular entity with a few physical parameters, which can then be compared to the parameters of other molecules to determine solubility with no need for an explicit calculation of their interaction. Classical molecular dynamics has been used with this theme in the past to compute solubility parameters,1,20 but these methods often suffer from the inaccuracies of the force field, as well as the ambiguous nature of differentiating between noncovalent interaction types. Recently, high-accuracy quantum-chemical methods utilizing symmetry-adapted perturbation theory (SAPT) have been developed that are capable of reproducing intermolecular potentials to within chemical accuracy.35,36 Because describing the nature and magnitude of intermolecular interactions is of fundamental importance to any

(3)

where V1 is the segmental volume of component 1, and δ1 and δ2 are the Hildebrand parameters of component 1 and component 2, respectively. These results are fascinating because they allow for a rigorous link between the magnitude of microscopic interactions and the macroscopic miscibility of a mixture. Utilizing eqs 2 and 3, a simple result can be deduced: two molecules with similar solubility parameters (δ1 − δ2 ≈ 0) should be mutually soluble. If the difference of solubility parameters is close to zero, then the solution behaves like an ideal solution (ΔHmixing → 0), where components are always miscible (ΔGm < 0) due to the negative entropy of mixing (assuming ideal mixing). The strength of the solubility parameter approach is clear: solubility parameters provide an experimentally determined value for a given molecule that characterizes its solubility behavior when mixed with another molecular species without any explicit knowledge of the details of their interaction. While Hildebrand solubility parameters are simple and intuitively satisfying, upon inspection it becomes clear that this approach is lacking. Simply speaking, there is no explicit information in Hildebrand solubility parameters regarding the nature of the particular noncovalent interactions present. With this in mind, it is no surprise that Hildebrand parameters fail notably for mixtures involving polar, protic, or anisotropic molecules, where the assumption of purely random dispersive forces as described by Hildebrand and Scatchard31 is no longer valid. From this, it becomes clear that differences in 5195

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

induction terms. The total DFT-SAPT interaction energy, Eij, is represented by eq 7.

solubility parameter approach, it is logical to consider the application of quantum-chemical methods to these issues. In this work, we seek to combine the simple and direct nature of solubility parameters with the high accuracy of firstprinciples quantum-chemical calculations for the purpose of estimating the molecular solubility of nonelectrolytes. These methods will be shown to provide a priori, accurate predictions of cohesive energies (solubility parameters) with no need for empirical parametrization. Further, it will be seen that our approach provides a more rigorous decomposition of the cohesive energy into its different noncovalent contributions, furthering physical intuition with regard to the energetics of molecular solubility and solidifying the approach of Hansen. Finally, we will demonstrate that this approach can be utilized to explicitly compute the interaction energy between two molecular species in solution, providing a complementary approach to COSMO-RS that allows for the rigorous decomposition of interaction types, leading to a more intuitive understanding of solubility through the use of distinct Hansenlike solubility parameters. We will demonstrate these improvements through the analysis of a set of common nonionic organic molecules using high-level theory to analyze solubility energetics, while simultaneously staying grounded in the “like dissolves like” mentality central to every chemist’s intuition.

(1) (2) (2) (2) (2) (1) E ij = Eelec + Eexch + E ind + Eexch −ind + Edisp + Eexch−disp

+ δ HF

While a number of options exist for partitioning the perturbative contributions of eq 7 into physically meaningful terms, we use a grouping similar to that recently employed in the literature.40 The inclusion of E(2) exch−disp in Edisp and of E(2) in E is partially rationalized by the cancellation of exch−ind ind artificial divergences at small intermolecular separations between the paired terms within the DFT-SAPT energy:



all molecules

∑ i,j

1 E ij 2

(1) Eelec = Eelec

(8)

(2) (2) E ind = E ind + Eexch −ind + δ HF

(9)

(2) (2) Edisp = Edisp + Eexch −disp

(10)

(1) Eexch = Eexch

(11)

It becomes immediately apparent from eqs 8−11 that the DFTSAPT decomposition of the total interaction energy bears similarities to the multicomponent Hansen solubility model described earlier. There are differences between the exact terms that appear, but the fact that DFT-SAPT decomposes the energy into different noncovalent interaction contributions makes it particularly applicable to solubility theory, where one desires to reduce any theoretical result back to “like dissolves like”, under some approximation. The relationship between the DFT-SAPT energies and Hansen parameters will be made concrete in a later section. DFT-SAPT Hildebrand Solubility Parameters. To compute the cohesive energy density of a molecular liquid, and subsequently the Hildebrand parameter, a DFT-SAPT calculation must be performed between every pair of molecules in the system (eq 6). This, however, is impractical, particularly for nonpolar organic systems where intermolecular interactions are likely short-range, and not communicated far beyond the first shell of nearest neighbors. Under this assumption, the total interaction energy of a chosen molecule with its nearest neighbors, Eiint, is the sum of individual pairwise SAPT interactions between the center molecule and its m nearest neighbors (eq 12).

THEORY Symmetry Adapted Perturbation Theory. To provide a rigorous grounding for solubility parameter theory, the cohesive energy density of a molecular liquid must be computed with high accuracy, which amounts to determining ΔHvap. In this work, we approximate ΔHvap as ΔEvap, although we are aware that this leaves out the expansion work term in the heat of vaporization. We are motivated by previous work that has employed this approximation successfully for very similar molecular species.20 A pairwise interaction description is used to compute the vaporization energy of a homogeneous molecular liquid (eq 6), where El is the total energy of the interacting molecular liquid, Eg is the total energy of every molecule in the liquid with no intermolecular interactions present, and Eij is the intermolecular energy of interaction between molecule i and molecule j. Note that for nonpolar molecules at their standard liquid densities Eij is always a stabilizing term, and thus, throughout this paper the convention will be that Eij is a negative number. ΔEvap = Eg − E l =

(7)

(6)

m

E iint = −∑ E ij

To compute Eij, it is necessary to choose a method wellsuited for calculating intermolecular interaction energies. SAPT, in particular the density functional theory version37 (DFTSAPT), is appropriate for this application, as it is firstprinciples, high-accuracy, and free from basis set superposition error, requiring no additional calculations to compute a counterpoise correction.38,39 Additionally, DFT-SAPT computes the total interaction energy between closed-shell molecules as a sum of electrostatic, inductive, dispersion, and exchange energies arising from different terms in a perturbative series. This naturally allows for a decomposition of the total interaction energy between two molecules into contributions arising from distinct noncovalent interaction types. In this study, we use second-order DFT-SAPT with an uncorrelated Hartree−Fock correction (δHF) to account for higher-order

(12)

j

To compute the total cohesive energy of the molecular liquid, we then sum Eint i over all molecules in the liquid and divide by two to remove double counting. However, if all molecules in a homogeneous solution can be considered identical, then computing the energy of vaporization amounts to multiplying Eint i by the number of molecules N, and dividing by 2 (eq 13). ΔEvap =

N × E iint 2

(13)

To determine the cohesive energy density of the system, we divide eq 13 by the molar volume, Vm (eq 14). 5196

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B 2 δSAPT =

ΔEvap Vm

=

int 1 N × Ei 2 Vm

Article

this interaction energy explicitly. In Flory−Huggins theory, prior to making the solubility parameter approximation of eq 3, one can write the enthalpic energy change upon mixing using eqs 18 and 19, where ΔHm is the enthalpic change upon mixing component 1 with component 2, z is the number of nearest neighbors, and Δw is a mean-field approximation to the interaction energy when one component 1 contact, w11, and one component 2 contact, w22, are broken, and one component 1−component 2 contact, w12, is formed.

(14)

δ2SAPT

If we wish to compute from a microscopic viewpoint, without knowledge of the molar volume, we can switch to a molecular view and rewrite eq 14 as the more useful eq 15 using the molecular volume Vmolecule (Vm/N). int

2 δSAPT =

1 Ei 2 Vmolecule

(15)

ΔHm = N1ϕ2zΔw

From this formula, one can use the interaction energy of a single molecule with its nearest neighbors, along with its molecular volume, to compute the Hildebrand solubility parameter. In principle the cohesive energy density of a molecular liquid can be computed a priori with no knowledge of the experimental heat of vaporization or molar volume. DFT-SAPT Hansen Solubility Parameters. Hansen solubility parameters are computed analogously to the Hildebrand solubility parameters by substituting eq 7 into eq 15, combining Eelec and Eind into a single term that describes all “polar” interactions, and requiring that the sum of the squares of these parameters be equal to the original DFT-SAPT cohesive energy density, as in Hansen theory. We note that for nonpolar molecules, Eelec, Eind, and Edisp are stabilizing energies (negative number), while Eexch is always a destabilizing term (positive number). E int Vmolecule

=

−Edisp −(Eelec + E ind) Eexch + − Vmolecule Vmolecule Vmolecule

2 2 2 2 δSAPT = δpolar + δdisp − δexch

(18)

1 (w11 + w22) (19) 2 While Flory−Huggins theory utilizes a mean-field approximation to zΔw, it is in principle possible, with some model for the solubility process in mind, to compute this term explicitly. Using the DFT-SAPT formulation, one can compute Eint i for a component 1 molecule interacting in solution with component 2 molecules and use this value along with w11 and w22 (determined from the calculations of Eint for homogeneous i systems) to calculate zΔw. Using these methods and eq 20, the interaction parameter χ12 for use in Flory−Huggins theory can be directly calculated. Δw = w12 −

χ12 =



zΔw kBT

(20)

METHODS To determine the efficacy of our DFT-SAPT solubility parameter theory, we have computed the DFT-SAPT solubility parameters of 26 common organic nonionic small molecules, most of which are commonly used as solvents (Table 1). To calculate the solubility parameters for a given molecular species, a periodic 20 Å cubic box is filled with a number of molecules that matches the experimental room temperature liquid density of that molecule (Supporting Information). A periodic energy minimization is then performed to eliminate bad contacts, followed by a periodic NVT classical molecular dynamics simulation at 298 K using the OPLSAA200941 force field. In the case of water, the TIP3P force field is used. NVT simulations begin with 200 ps for thermal equilibration. After thermal equilibration, snapshots are taken of the trajectory every 10 ps for 300 ps. For each snapshot, a molecule is randomly chosen, and a shell of nearest neighbor molecules is delineated using a 4.5 Å cutoff. This seemingly short cutoff corresponds to a spatial search for atoms within 4.5 Å of atoms on the chosen molecule, and all unique molecules fulfilling this cutoff criteria are delineated. Because this is a quantum-mechanical and not pairwise-classical calculation, the spatial search cutoff leads to an effective interaction cutoff of 4.5 Å + ∼2Rmol, where Rmol is the radius corresponding to the molecular volume of a given molecule. Longer cutoffs were explored, but computed energies (Eint i ) were observed generally to converge to within ∼3 kJ/mol using a 4.5 Å cutoff, with exceptions for strongly polar molecules (notably water), as expected (see Supporting Information for detailed discussion of method and cutoff). On each snapshot, a DFT-SAPT (LPBE0AC/aug-cc-pVDZ) calculation is performed between the chosen molecule and a particular nearest neighbor molecule. This calculation is then iterated through every nearest neighbor molecule surrounding the chosen molecule. For each snapshot, a total interaction

(16) (17)

In this DFT-SAPT formulation of the Hansen parameters (eqs 16 and 17), it is interesting to note the differences from the original definition (eq 4). The dispersion term is essentially identical in meaning to that described by Hansen; however, the polar and hydrogen bonding terms are different, particularly because the DFT-SAPT formulation does not include a specific hydrogen bonding term. In the following analysis we will show that considering the magnitude of the exchange term, in conjunction with the polar term, allows for a quantitative identification of the hydrogen bonding tendency of a molecule. The argument goes as follows, and will be expanded upon further: Hydrogen bonding is an electrostatically dominated process, so when it occurs, we expect the polar interaction energy to be very large. However, when a hydrogen bonding interaction is present, the intermolecular distance between the interacting atoms will be much smaller than for typical weak polar or dispersive interactions. As a result, the increased wave function overlap of the two atoms will result in an increase in the exchange repulsion that is concomitant with the large polar stabilization term, as is true of the Coulombic and exchange energies whenever a covalent bond is formed. Thus, one will be able to distinguish between polar protic and polar aprotic molecules by observing the magnitude of the exchange energy determined by DFT-SAPT and its simultaneous rise with the polar parameter. Interactions between Two Molecular Species. Ultimately, the solubility parameter expression for the interaction parameter χ12 (eq 3) is plagued by the inaccuracies of the geometric mean approximation for the interaction energy and is not valid for polar systems. It would be ideal if one could bypass any approximations to the interaction parameter and compute 5197

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

is superior to that which employs computed molecular volumes. The reason for this is simple: computed molecular volumes using the Voronoi method underestimate the experimental molecular volumes by 50−80%. There is a significant amount of empty space in the condensed molecular liquid that is not described by pure van der Waals radii and must be captured to reproduce experimental values. We note that, for the remainder of this work, molecular volumes will be derived from the experimental density and molar mass of the molecules. The computed DFT-SAPT Hildebrand solubility parameters that employ experimentally derived molecular volumes (Figure 1, red circles) follow a highly linear correlation with their experimental values (R2 = 0.95), demonstrating that DFTSAPT accurately reproduces the relative cohesive energies of a broad set of molecules. Given the experimental error inherent in Hildebrand parameter determination, we consider the majority of our results to be within experimental error, with a consistent slight underestimation due to the neglect of interactions outside of our chosen cutoff radius. Obvious significant deviations occur for molecules such as water, where the nearest neighbor pairwise scheme does not describe the entirety of the liquid’s interactions, and the cohesive energy is underestimated.48 Previous SAPT estimations of many-body interactions in water account for approximately 28% of the total interaction energy49 (in good agreement with our 33% underestimation of the experimental Hildebrand parameters for water, though some of this disagreement is likely the result from neglecting interactions outside of the cutoff radius chosen in this work), suggesting that an extension of our methodology to higher-order many-body interactions is possible for improving accuracy in strongly interacting neutral, and even ionic, systems, albeit at a higher computational cost. Indeed, the fact that our description lacks these many-body interactions is the reason that the data fits in Figure 1 do not go precisely through the origin; molecules where higher order corrections are necessary result in an underestimation of the computed Hildebrand parameters, pulling fits away from linearity through the origin. Regardless of these facts, our current pairwise implementation seems to be an accurate and useful approach for determining Hildebrand solubility parameters and cohesive energies, as it is fairly low cost and requires only the molecular structure, liquid density, and ionization potential as computational inputs, the latter of which can be accurately calculated using a method like CCSD(T) with a large basis set. For novel molecular species in which the experimental density is not known, one could consider performing NPT simulations and utilizing the simulated molar volumes to determine the solubility parameters described in this work in a similar fashion. This further extends the potential of this approach to any case where only the molecular structure (and implicitly, an accurate classical force field for that structure) is known. It is useful to mention that the magnitude of the cohesive energy is reasonably insensitive to the value of the asymptotic shift applied in the DFT-SAPT methodology (∼0.2 eV shift in asymptotic correction leads to ∼0.3 kJ/mol shift in average cohesive energy). This result is encouraging due to the experimental uncertainty present in many reported gas-phase ionization potentials and suggests that this method could be applied accurately to any molecule in situations where highresolution gas-phase measurements are not available, for example, cyclic voltammetry measurements of organic semiconductors.50 Finally, we would like to stress that while

energy between the molecule and all of its nearest neighbors is computed, corresponding to eq 12, which is then decomposed into its DFT-SAPT contributions. The total interaction energy and its components are then averaged over all 30 snapshots of the trajectory to generate a statistically averaged value of the interaction energy. This statistically averaged interaction energy and its components, when divided by molecular volumes, are then used to determine the DFT-SAPT solubility parameters described previously. Regarding the electronic structure, it has been noted that DFT-SAPT calculations are reasonably insensitive to the choice of functional as along as an appropriate asymptotic correction is applied,42 thus our choice of LPBE0AC.37,43 Asymptotic corrections are performed using the monomer highest occupied molecular orbital (HOMO) energy level calculated at LPBE0AC/aug-cc-pVDZ on the OPLSAA2009 minimized geometry and the experimental gas-phase ionization potential (Supporting Information). MOLPRO201244 is used for all density fitting DFT-SAPT calculations. Two methods are employed to compute the molecular volume: The Voronoi method in Chimera 45 is used for the first-principles computation of molecular volumes, and more accurate molecular volumes are derived from the experimental densities and molar masses. The TINKER46 Molecular Modeling package is used for the classical molecular dynamics simulations. The Bussi thermostat is used to maintain constant temperature. A cutoff of 12 Å is used for van der Waals interactions, 8.0 Å for Ewald summation, and 15 Å for all other interactions.



RESULTS Computed DFT-SAPT Hildebrand Parameters. In Figure 1, we have plotted two sets of computed DFT-SAPT

Figure 1. Computed DFT-SAPT Hildebrand parameters versus experimental Hildebrand parameters. Red circles represent Hildebrand parameters computed using experimentally derived molecular volumes; blue triangles represent parameters obtained using computed molecular volumes. Black lines represent linear fits to the data.

Hildebrand solubility parameters (eq 15) versus the experimentally derived Hildebrand parameters from the literature.47 The blue data set represents Hildebrand parameters computed with molecular volumes determined via the Voronoi method, while the red data set represents parameters computed with experimentally derived molecular volumes (Supporting Information). While both approaches capture the qualitative correlation between experimental and computational values, the data set that employs experimentally derived molecular volumes 5198

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

agreeing with common intuition. Specifically, we note that for conventional nonpolar molecules, dispersion interactions contribute ∼60% of the total liquid stabilization energy. Moving further up Table 1, following decreasing fd values, dispersion interactions contribute less to the overall cohesive energy of the system, and polar interactions become more prominent. These molecules include 1,2-dichlorobenzene, acetone, dimethylformamide, propylene carbonate, and acetonitrile, molecules typically classified as polar aprotic solvents. Supporting the characterization of these molecules as polar but aprotic, we note that while δpolar increases relative to those molecules found in the nonpolar section, δexch remains relatively constant. This supports our previous conjecture that δexch can be used as a proxy for describing hydrogen bonding interactions in a system, as conventionally all of these molecules are thought to be predominantly aprotic. We note that thorough characterizations of hydrogen bonding using DFTSAPT exist, and we direct the reader to these.52 It is important to make clear that the solutions studied in this work are homogeneous. It is known that many of the molecules in the polar aprotic section of Table 1 are known to exhibit weak hydrogen bonding when mixed with other molecular species, despite the parameters of Table 1 implying that no hydrogen bonding is present in a homogeneous solution. If one wanted to explicitly compute the hydrogen bonding character of a mixture of two molecules, say propylene carbonate and water, this could be simply done, and δexch and δpolar could be monitored in the same fashion for the mixture to determine hydrogen bonding character. By looking at a variety of mixtures, one could deduce the hydrogen bonding tendencies of one molecule over a range of other molecules, and this is in fact the spirit of the original Hansen parameters. The DFTSAPT approach described here can be used to rigorously determine these parameters over a range of interactions, as in Hansen theory, providing a quantitative foundation for this type of approach. This topic will be discussed further in the following section. Approaching the top of Table 1, representing polar protic species and the smallest values of fd, δpolar increases rapidly, far above that for nonpolar and polar aprotic solvents. We note that all of the molecules appearing in this region of Table 1 are small molecules known to exhibit significant hydrogen bonding (water, formamide, methanol). As we have proposed, concomitant with the increase in δpolar in this region is a rapid increase in δexch, signifying the formation of hydrogen bonds in these molecular liquids. Thus, the DFT-SAPT method is capable of distinguishing between polar protic molecules and polar aprotic molecules by analyzing the magnitudes of δpolar and δexch. By using the DFT-SAPT methods described in this work, it is possible not only to quantitatively predict the magnitude of the cohesive energy of a molecular liquid (Hildebrand parameter), but also to quantitatively describe the noncovalent interaction character of a molecular liquid (fd, δpolar, δdisp, δexch). This is particularly important for molecular solubility, where commonly “like dissolves like” is used to choose solvent types, but there is little to no quantitative formalism to describe this approach. Hansen parameters serve an important role trying to quantify this relationship but are highly empirical and dependent upon extensive parametrization. Here, we present a first-principles method capable of determining the noncovalent interaction character of a particular small molecule quantitatively. We have written the DFT-SAPT decomposition

Hildebrand solubility parameters are not highly useful in regular solution theories of polar molecules, accurate knowledge of the cohesive energy of a cluster of molecules is an important quantity, particularly in many technological applications where solution processability is fundamental51 and side chains are often used to break up cluster aggregation (decrease cohesive energy) and increase solubility. Computed DFT-SAPT Hansen Parameters. In Table 1, we have computed three component DFT-SAPT Hansen solubility parameters (eq 16), as well as fd, the percent contribution of dispersion energy to the overall stabilization of the molecular liquid (eq 21). fd =

100 × Edisp Edisp + Eelec + E ind

(21)

By determining fd, we can simply quantify the nonpolar/polar character of a molecule and see if DFT-SAPT reproduces common trends in solvent classification, which is important for the verification of this approach’s qualitative merits for determining solubility. Particularly, we seek to classify solvent types a priori (nonpolar, polar aprotic, polar protic), and lend quantitative rigor to this traditionally qualitative approach. The bottom of Table 1 represents molecules where dispersion forces are the dominant noncovalent interaction (large fd value) contributing to the cohesion of the molecular liquid. Molecules traditionally classified as nonpolar and dominated by dispersion interactions (hexane, benzene, toluene, diethyl ether) are found in the bottom of the table, Table 1. DFT-SAPT Hansen Parameters and Molecule Classification molecule water formamide ethylene glycol methanol ethanolamine ethanol isopropanol methylacetamide nitromethane acetonitrile propylene carbonate dichloromethane chloroform 1,2-dichlorobenzene dimethylformamide acetone 1,4-dioxane pyridine chlorobenzene benzene tetrahydrofuran 1,2-dimethylbenzene toluene diethyl ether pentane hexane

δpolar

δdisp

Polar Protic 31.8 17.0 24.4 15.6 21.6 14.9 16.9 12.2 20.0 14.9 15.2 12.0 15.0 11.6 15.4 12.9 Polar Aprotic 13.0 12.0 11.3 10.6 12.6 12.4 10.5 11.0 10.4 10.9 11.6 11.9 10.5 11.4 9.1 10.1 10.6 11.8 10.4 12.1 10.2 11.8 Nonpolar 9.4 11.4 9.0 10.9 9.3 11.3 9.0 11.0 7.7 9.5 7.5 9.3 6.8 9.0

δexch

fd (%)

31.5 24.0 21.3 17.1 20.4 15.7 15.3 15.5

22 28 32 33 35 37 37 40

13.0 11.3 12.6 11.8 11.9 11.7 10.8 9.6 11.5 11.4 11.5

45 45 48 50 51 51 52 53 54 56 56

10.7 10.1 10.5 10.2 8.9 8.7 7.9

58 58 58 58 58 58 61 5199

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

in a form akin to traditional Hansen parameters for familiarity, and these parameters can be used in the same fashion as traditional Hansen solubility parameters, but this is not an absolute means of thinking about solubility. In fact, fd may be the most useful parameter involved in this study, as it reproduces the classic solvent categories nearly perfectly and assigns a simple quantitative value to assess the polar or nonpolar character of a molecule. We note that the greatest utility of this methodology could be in the screening of potential solvents for organic semiconductors,53 as well as the prediction of solubility in molecular pharmaceuticals. Computing the Explicit Interaction between Two Molecular Species. While solubility parameters can be a simple and useful method for determining the miscibility of two molecules, utilizing solubility parameters for the prediction of miscibility relies entirely on using the properties of each individual molecule to predict their interaction, without any explicit knowledge of this interaction. It would be useful if the previously described method could be used to compute the energy of interaction between two molecular species, zΔw, as this would allow for the precise determination of interaction parameters used in solution theories, in particular the Flory interaction parameter. Additionally, this would allow for the solubility parameter approach to interaction energies to be verified by explicit calculation of that interaction. To further demonstrate the utility of the DFT-SAPT approach to solubility energetics, we have chosen to compute explicitly the interaction between solute and solvent in two scenarios. The solute to be studied is ethylene glycol (EG), and the two solvent systems it will interact with are ethanol (good solvent) and benzene (bad solvent). Ethylene glycol is experimentally known to be miscible in ethanol and immiscible in benzene. It is well-known in the literature that Flory interaction parameters can exhibit a strong dependence on concentration and temperature.54,55 In principle, the DFTSAPT approach can compute the energy of interaction under any concentration or temperature constraints, as long as the appropriate statistical averaging is performed. However, it is instructive to choose a simple idealized scenario for computing the interaction energy (Figure 2), and we compute zΔw according to eq 22. zΔw = zw12 + zw21 − zw11 − zw22

Figure 2. Schematic solubility process utilized for the computation of zΔw. Green and red circles represent two separate molecular species. In practice, use of periodic boundary conditions and a finite box are employed for the computation of each state, assuring convergence.

Table 2. Computed Two-Component Interaction Energies mixture

fd (%)

ethanol in EG (zw12) EG in ethanol (zw21) benzene in EG (zw12) EG in benzene (zw21) EG in EG (zw11) ethanol in ethanol (zw22) benzene in benzene (zw22)

37 32 56 52 32 37 58

ethanol/EG benzene/EG

Eint i (kJ/mol) −87.6 −100.6 −95.9 −65.0 −118.3 −70.9 −85.6 zΔw (kJ/mol) 1.1 43.1

exhibit fd values (32/37) far closer to those of the respective homogeneous solutions (32/37), whereas the benzene/ethylene glycol mixtures exhibit fd values (52/56) far different from that of the homogeneous ethylene glycol fd. In the case of ethylene glycol/benzene, the necessary dispersion forces exist for solubility, but the polar terms necessary for ethylene glycol’s solubility are not provided by benzene. Consequently, the fd of benzene/ethylene glycol is dispersion dominated. Because ethylene glycol and benzene are immiscible, the use of fd as a parameter to describe “like dissolves like” demonstrates some utility as a quantitative manifestation of this “like dissolves like” rule of thumb. The Table 2 values of zΔw clearly show that the ethanol/ ethylene glycol mixture is far more stable energetically than the benzene/ethylene glycol mixture. While the ethylene glycol/ ethanol mixing energy is slightly endothermic (less than room temperature kT), the benzene/ethylene glycol mixing process is 40 times as endothermic, showing that it is enthalpically impossible for any significant mixing to occur at room temperature. If one assumes that ideal entropies of mixing are valid, uses values of zΔw to determine χ12 via eq 20, and plots the Flory−Huggins expression (eq 2) for the free energy of mixing of ethylene glycol/ethanol and ethylene glycol/benzene mixtures, one finds that the former mixture is miscible over all concentration ranges, while the latter mixture is immiscible over all concentration ranges (Figure 3), in agreement with the known miscibilities. While this is an extremely simple example,

(22)

In the calculation of zΔw, a molecular dynamics trajectory is run of one molecule in a solution of the other molecule, and vice versa. The analysis of the interaction energy parallels the discussion of the previous sections. Because ethanol, benzene, and ethylene glycol are all of similar molecular volumes, we perform the molecular dynamics simulation by eliminating one molecule of solvent from our simulation box and replacing that molecule with a solute molecule. The results of these simulations are presented in Table 2, with computed interaction energies listed, as well as the fd value for a particular interaction. Prior to explicitly computing zΔw, we examine the values of fd in Table 2 pertaining to systems zw11, zw22, zw12, and zw21 for both the EG/ethanol and EG/benzene mixtures. This allows one to evaluate how the “like dissolves like” mentality manifests itself in the explicit evaluation of the interaction of two molecular species. Homogeneous liquids of ethylene glycol, ethanol, and benzene exhibit fd parameters of 32, 37, and 58, respectively (Table 1). Looking at the value of fd for interacting species, it is clear that the ethanol/ethylene glycol mixtures 5200

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B



Article

ASSOCIATED CONTENT

S Supporting Information *

Energy levels used for DFT-SAPT asymptotic corrections; liquid densities used for computation and molecular dynamics simulation details; computed molecular volumes; Figure 1 data and fit details; experimental Hansen parameters; all DFT-SAPT computed values; all molecular dynamics trajectories; discussion of cutoff length for DFT-SAPT + MD methodology. This material is available free of charge via the Internet at http://pubs.acs.org.



Figure 3. Free energy of mixing for ethylene glycol/ethanol (red circles) and ethylene glycol/benzene (blue triangles) mixtures computed using Flory−Huggins regular solution theory.

AUTHOR INFORMATION

Corresponding Author

*Phone: 847-491-5652. E-mail: [email protected]. Notes

The authors declare no competing financial interests.



we provide it as a demonstration of how the DFT-SAPT approach could be used to compute interaction parameters, and then use these values to determine macroscopic miscibilities. As a final note, we stress that we have chosen an idealized scenario (Figure 2) under which to compute zΔw, and thus χ12. As mentioned previously, χ12 is known to exhibit a range of values as a function of concentration and temperature. With the DFT-SAPT methodology described here, it is possible to explicitly compute all concentration and temperature ranges of this parameter, if desired. This would involve doing DFT-SAPT analysis on different volume fractions of solute/solvent, as well as temperature ranges, and would require a large number of simulations. Nevertheless, the interaction parameter could be determined for any constraints on concentration and temperature.

ACKNOWLEDGMENTS N.E.J. thanks the NSF for the award of a Graduate Research Fellowship (NSF DGE-0824162). N.E.J. would like to thank Brett Savoie and Kevin Kohlstedt for useful discussion and Tomekia Simeon for help with the SAPT calculations. M.A.R. thanks the NSF (NSF CHE-1058896) for support.



REFERENCES

(1) Gupta, J.; Nunes, C.; Vyas, S.; Jonnalagadda, S. Prediction of Solubility Parameters and Miscibility of Pharmaceutical Compounds by Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 2014−2023. (2) Tian, Y.; Booth, J.; Meehan, E.; Jones, D. S.; Li, S.; Andrews, G. P. Construction of Drug−Polymer Thermodynamic Phase Diagrams Using Flory−Huggins Interaction Theory: Identifying the Relevance of Temperature and Drug Weight Fraction to Phase Separation within Solid Dispersions. Mol. Pharmacol. 2013, 10, 236−248. (3) Hancock, B. C.; York, P.; Rowe, R. C. The Use of Solubility Parameters in Pharmaceutical Dosage Form Design. Int. J. Pharm. 1997, 148, 1−21. (4) Trevino, S. R.; Scholtz, J. M.; Pace, C. N. Measuring and Increasing Protein Solubility. J. Pharm. Sci. 2008, 97, 4155−4166. (5) Machui, F.; Abbott, S.; Waller, D.; Koppe, M.; Brabec, C. J. Determination of Solubility Parameters for Organic Semiconductor Formulations. Macromol. Chem. Phys. 2011, 212, 2159−2165. (6) Chueh, C.-C. Non-Halogenated Solvents for Environmentally Friendly Processing of High-Performance Bulk-Heterojunction Polymer Solar Cells. Energy Environ. Sci. 2013, 6, 3241−3248. (7) Hildebrand, J. H. The Solubility of Non-Electrolytes; Reinhold: New York, 1950. (8) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NY, 1953. (9) Hansen, C. M. Hansen Solubility Parameters: A User’s Handbook; CRC Press: Boca Raton, FL, 2007. (10) Kumar, S. K.; Szleifer, I.; Panagiotopoulos, A. Z. Determination of the Chemical Potentials of Polymeric Systems from Monte Carlo Simulations. Phys. Rev. Lett. 1991, 66, 2935−2938. (11) Hildebrand, J. H. Solubility. XII. Regular Solutions1. J. Am. Chem. Soc. 1929, 51, 66−80. (12) Abraham, M. H.; Le, J. The Correlation and Prediction of the Solubility of Compounds in Water Using an Amended Solvation Energy Relationship. J. Pharm. Sci. 1999, 88, 868−880. (13) Ruelle, P.; Kesselring, U. W. Prediction of the Aqueous Solubility of Proton-Acceptor Oxygen-Containing Compounds by the Mobile Order Solubility Model. J. Chem. Soc., Faraday Trans. 1997, 93, 2049−2052. (14) Ran, Y.; He, Y.; Yang, G.; Johnson, J. L. H.; Yalkowsky, S. H. Estimation of Aqueous Solubility of Organic Compounds by Using the General Solubility Equation. Chemosphere 2002, 48, 487−509.



CONCLUSION We have developed a computational approach that allows for a rigorous and quantitative determination of the cohesive energy of a molecular liquid and, consequently, the Hildebrand solubility parameter. It is found that solubility parameters computed using the DFT-SAPT method developed herein exhibit good agreement with those determined experimentally. Additionally, we have utilized the natural energy decomposition of DFT-SAPT to develop Hansen-like solubility parameters, representing the various types of noncovalent interactions important to solubility processes. We note explicitly the strong apparent correlation between the first-order exchange energy and the hydrogen bonding character of a liquid. The DFTSAPT method is found to be capable of reproducing general solvent categories (nonpolar, polar aprotic, polar protic) and assigns a variety of quantitative values that help quantify the classic “like dissolves like” approach to solubility. Specifically, we have introduced a parameter fd that quantifies the noncovalent interaction character of a solvent and can be computed for any closed-shell molecular species. Finally, we have utilized the DFT-SAPT methodology to explicitly calculate the interaction between two molecular species and shown it to accurately describe the miscibility of two representative systems, an ethylene glycol/ethanol and an ethylene glycol/benzene mixture. Future work will focus on incorporating many-body corrections and applying this methodology to organic semiconductor and molecular pharmaceutical systems. 5201

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202

The Journal of Physical Chemistry B

Article

(36) Chałasiński, G.; Szczȩsń iak, M. M. State of the Art and Challenges of the Ab Initio Theory of Intermolecular Interactions. Chem. Rev. 2000, 100, 4227−4252. (37) Heßelmann, A.; Jansen, G.; Schütz, M. Density-Functional Theory-Symmetry-Adapted Intermolecular Perturbation Theory with Density Fitting: A New Efficient Method to Study Intermolecular Interaction Energies. J. Chem. Phys. 2004, 122. (38) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of Van Der Waals Complexes. Chem. Rev. 1994, 94, 1887−1930. (39) Szalewicz, K. Symmetry-Adapted Perturbation Theory of Intermolecular Forces. WIREs Comput. Mol. Sci. 2012, 2, 254−272. (40) McDaniel, J. G.; Schmidt, J. R. Physically-Motivated Force Fields from Symmetry-Adapted Perturbation Theory. J. Phys. Chem. A 2013, 117, 2053−2066. (41) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (42) Misquitta, A. J.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. Intermolecular Potentials Based on Symmetry-Adapted Perturbation Theory with Dispersion Energies from Time-Dependent DensityFunctional Calculations. J. Chem. Phys. 2005, 123. (43) Leforestier, C.; Tekin, A.; Jansen, G.; Herman, M. First Principles Potential for the Acetylene Dimer and Refinement by Fitting to Experiments. J. Chem. Phys. 2011, 135. (44) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; et al. MOLPRO, A Package of Ab Initio Programs, version 2012.1; MOLPRO, 2012. (45) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF ChimeraA Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (46) Ponder, J. TINKERSoftware Tools for Molecular Design, version 6.3; Washington University School of Medicine: Saint Louis, MO, 2014. (47) Barton, A. F. M. CRC Handbook of Solubility Parameters and Other Cohesion Parameters; CRC Press: Boca Raton, 1991. (48) Goldman, N.; Saykally, R. J. Elucidating the Role of Many-Body Forces in Liquid Water. I. Simulations of Water Clusters on the VRT(ASP-W) Potential Surfaces. J. Chem. Phys. 2004, 120, 4777− 4789. (49) Milet, A.; Moszynski, R.; Wormer, P. E. S.; van der Avoird, A. Hydrogen Bonding in Water Clusters: Pair and Many-Body Interactions from Symmetry-Adapted Perturbation Theory. J. Phys. Chem. A 1999, 103, 6811−6819. (50) Cardona, C. M.; Li, W.; Kaifer, A. E.; Stockdale, D.; Bazan, G. C. Electrochemical Considerations for Determining Absolute Frontier Orbital Energy Levels of Conjugated Polymers for Solar Cell Applications. Adv. Mater. 2011, 23, 2367−2371. (51) Mei, J.; Bao, Z. Side Chain Engineering in Solution-Processable Conjugated Polymers. Chem. Mater. 2013, 26, 604−615. (52) Hoja, J.; Sax, A. F.; Szalewicz, K. Is Electrostatics Sufficient to Describe Hydrogen-Bonding Interactions? Chem.Eur. J. 2014, 20, 2292−2300. (53) Duan, C. Toward Green Solvent Processable Photovoltaic Materials for Polymer Solar Cells: The Role of Highly Polar Pendant Groups in Charge Carrier Transport and Photovoltaic Behavior. Energy Environ. Sci. 2013, 6, 3022−3034. (54) Schuld, N.; Wolf, B. A. Solvent Quality as Reflected in Concentration- and Temperature-Dependent Flory−Huggins Interaction Parameters. J. Polym. Sci., Part B: Polym. Phys. 2001, 39, 651− 662. (55) Koch, T.; Strobl, G. R. Concentration Dependence of the FloryHuggins Interaction Parameter of a Polymer Blend as Determined by Small-Angle X-ray Scattering Experiments. J. Polym. Sci., Part B: Polym. Phys. 1990, 28, 343−353.

(15) Fredenslund, A.; Jones, R. L.; Prausnitz, J. M. GroupContribution Estimation of Activity Coefficients in Nonideal Liquid Mixtures. AIChE J. 1975, 21, 1086−1099. (16) Duchowicz, P. R.; Castro, E. A. Qspr Studies on Aqueous Solubilities of Drug-Like Compounds. Int. J. Mol. Sci. 2009, 10, 2558− 2577. (17) Balakin, K.; Savchuk, N.; Tetko, I. In Silico Approaches to Prediction of Aqueous and DMSO Solubility of Drug-Like Compounds: Trends, Problems and Solutions. Curr. Med. Chem. 2006, 13, 223−241. (18) Delaney, J. S. Predicting Aqueous Solubility from Structure. Drug Discovery Today 2005, 10, 289−295. (19) Amidon, G. L.; Yalkowsky, S. H.; Anik, S. T.; Valvani, S. C. Solubility of Nonelectrolytes in Polar Solvents. V. Estimation of the Solubility of Aliphatic Monofunctional Compounds in Water Using a Molecular Surface Area Approach. J. Phys. Chem. 1975, 79, 2239− 2246. (20) Belmares, M.; Blanco, M.; Goddard, W. A.; Ross, R. B.; Caldwell, G.; Chou, S. H.; Pham, J.; Olofson, P. M.; Thomas, C. Hildebrand and Hansen Solubility Parameters from Molecular Dynamics with Applications to Electronic Nose Polymer Sensors. J. Comput. Chem. 2004, 25, 1814−1826. (21) Ahmed, A.; Sandler, S. I. Hydration Free Energies of Multifunctional Nitroaromatic Compounds. J. Chem. Theory Comput. 2013, 9, 2774−2785. (22) Mobley, D. L.; Dumont, É.; Chodera, J. D.; Dill, K. A. Comparison of Charge Models for Fixed-Charge Force Fields: SmallMolecule Hydration Free Energies in Explicit Solvent. J. Phys. Chem. B 2007, 111, 2242−2254. (23) Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. New Universal Solvation Model and Comparison of the Accuracy of the Sm5.42r, Sm5.43r, C-Pcm, D-Pcm, and Ief-Pcm Continuum Solvation Models for Aqueous and Organic Solvation Free Energies and for Vapor Pressures. J. Phys. Chem. A 2004, 108, 6532−6542. (24) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161−2200. (25) Klamt, A.; Eckert, F. COSMO-RS: A Novel and Efficient Method for the A Priori Prediction of Thermophysical Data of Liquids. Fluid Phase Equilib. 2000, 172, 43−72. (26) Kolár,̌ P.; Nakata, H.; Shen, J. W.; Tsuboi, A.; Suzuki, H.; Ue, M. Prediction of Gas Solubility in Battery Formulations. Fluid Phase Equilib. 2005, 228−229, 59−66. (27) Klamt, A.; Eckert, F.; Arlt, W. COSMO-RS: An Alternative to Simulation for Calculating Thermodynamic Properties of Liquid Mixtures. Annu. Rev. Chem. Biomol. Eng. 2010, 1, 101−122. (28) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074−5085. (29) Hildebrand, J. H.; Scott, R. L. The Entropy of Solution of Nonelectrolytes. J. Chem. Phys. 1952, 20, 1520−1521. (30) Flory, P. J. Thermodynamics of High Polymer Solutions. J. Chem. Phys. 2004, 10, 51−61. (31) Scatchard, G. Equilibria in Non-Electrolyte Solutions in Relation to the Vapor Pressures and Densities of the Components. Chem. Rev. 1931, 8, 321−333. (32) Lindvig, T.; Michelsen, M. L.; Kontogeorgis, G. M. A Flory− Huggins Model Based on the Hansen Solubility Parameters. Fluid Phase Equilib. 2002, 203, 247−260. (33) Hansen, C. M. 50 Years with Solubility ParametersPast and Future. Prog. Org. Coat. 2004, 51, 77−84. (34) Burgués-Ceballos, I.; Machui, F.; Min, J.; Ameri, T.; Voigt, M. M.; Luponosov, Y. N.; Ponomarenko, S. A.; Lacharmoise, P. D.; Campoy-Quiles, M.; Brabec, C. J. Solubility Based Identification of Green Solvents for Small Molecule Organic Solar Cells. Adv. Funct. Mater. 2013, 24, 1449−1457. (35) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. 5202

dx.doi.org/10.1021/jp5024197 | J. Phys. Chem. B 2014, 118, 5194−5202