First-principles prediction of liquid/liquid interfacial tension

: The interfacial tension between two liquids is the free energy per unit surface area required to create that interface. Interfacial tension is a determining factor for two-phase liquid behavior in a wide variety of systems ranging from water ﬂ ooding in oil recovery processes and remediation of groundwater aquifers contaminated by chlorinated solvents to drug delivery and a host of industrial processes. Here, we present a model for predicting interfacial tension from ﬁ rst principles using density functional theory calculations. Our model requires no experimental input and is applicable to liquid/liquid systems of arbitrary compositions. The consistency of the predictions with experimental data is signi ﬁ cant for binary, ternary, and multi-component water/organic compound systems, which o ﬀ ers con ﬁ dence in using the model to predict behavior where no data exists. The method is fast and can be used as a screening technique as well as to extend experimental data into conditions where measurements are technically too di ﬃ cult, time consuming, or impossible.


■ INTRODUCTION
The interfacial tension (IFT) is the free energy per unit surface area that is required for creating an interface between two condensed phases.The IFT is an important factor in a number of applications related to surfactant injection for improved oil recovery, 1 emulsion stability, 2,3 transport of organic contaminants through soil, 4,5 and ensuring water quality in aquifers, which includes the ability to predict the behavior of dense nonaqueous phase liquids such as chlorinated solvents. 5The IFT influences the capillary pressure between two nonmiscible fluids, which in turn impacts how fluids flow in porous media.Other applications where IFT plays a major role include micelle formation and other self-assembly processes of nanoparticles or colloidal particles at liquid/liquid interfaces. 6Self-assembly at liquid/liquid interfaces is an important process in materials science and technology, which is highly relevant to food supplements and cosmetics.The IFT also plays an important role in pharmaceutical applications such as drug delivery, 3,7,8 where surfactants help to form microemulsions, which improves the bioavailability of lipophilic or amphiphilic drugs.
Despite the wide range of applications, no general firstprinciples method for predicting IFT has been developed until now.In the literature, empirical models that are based on input from experimental data 9,10 and models that use molecular dynamics (MD) 11 have been reported.The empirical methods 9 are most commonly extrapolated from mutual solubilitites, 12,13 but this approach makes it difficult to generalize from binary systems to ternary or more complex systems.Classical force field MD has been used to estimate interfacial tension for several fluid systems 14−17 and tends to give reasonable results.On the other hand, MD methods depend on accurate model potentials, which are not always readily available.Furthermore, each MD simulation is time consuming, and a new simulation is required every time the composition or temperature is changed.The MD approach is also less straightforward for very low IFT values, 1−5 mN/m. 16n this paper, we present a method for predicting the IFT between two liquid phases of arbitrary composition using density functional theory (DFT) calculations combined with the COSMO-RS implicit solvent model.This means that no experimental input is required other than what has already been used for parametrizing COSMO-RS.It is important to note that no IFT data are used in the COSMO-RS parametrization.Our method is independent of the number of components in the system, and it is thus immediately applicable to multicomponent systems without modification.To the best of our knowledge, this is the first general model presented for predicting liquid/liquid IFT, which is based on quantum chemistry.In contrast to MD approaches, our method requires only a short series of DFT calculations to capture the most important conformers for each component in the system.No new DFT calculations are required when the composition or temperature is changed, only COSMO-RS calculations are required, and these generally take less than a minute on a standard PC.The method is thus efficient as a screening technique, for example, for testing the effectiveness of surfactants in reducing IFT.

■ METHODS
The COSMO-RS implicit solvent model is fundamentally different from other implicit solvent models used in quantum chemistry.Molecules are placed in a reference state that is a perfect conductor.This places screening charge density, σ, on the molecule surface.The screening surface, that is, the COSMO surface or as it is also called, the σ surface thus reflects the properties of the screened molecule (Figure 1) and is used as the base for further calculations of solvent interactions.By using a universal reference state for all molecules, it is possible to change the solvent composition without making any new DFT calculations.−27 The local solute/solvent interactions of the σ surfaces enable a molecule to interact with two different solvents over different areas on its surface, as shown in Figure 1.
All molecules investigated in this paper were taken from the COSMOtherm databases whenever they were available.All DFT calculations were performed using the Turbomole package, 28,29 Version 6.5 with the Becke−Perdew (BP) exchange−correlation functional 30,31 and Ahlrich-type basis sets. 32,33The f latsurf module of the COSMOtherm program 34 was used together with the BP_TZVP_C30_1301 and BP_TZVPD_FINE_HB2012_C30_1301 parametrizations to calculate the free energies of transferring molecules from one of the phases to an interface between two fluids.
The COSMO-RS f latsurf module calculates the solvation contribution to the change in free energy for transferring a molecule from a bulk phase to a flat and structureless two-phase interface.For each position and orientation, the difference in free energy compared with the bulk phase is calculated as well as the corresponding interfacial area taken up by the molecule.The total change in free energy, denoted G tot , is determined from the partition function of all surface contacting positions and orientations, and the corresponding expectation value of the replaced interfacial area, A av , is calculated.For COSMO-RS interactions, the COSMO surface area for water is scaled down by a factor of ∼1.5. 35For our IFT calculations (eqs 3−5), we used the nonscaled area for water because it is the physical area that is occupied at the interface that is relevant.
The interfacial tension between two phases, A and B, was determined using the following procedure: First, we calculated the partitioning of all components into the two phases using COSMO-RS.We then created a virtual surface phase, S, in the spirit of Guggenheim,36 and placed it between the bulk phases.It thus has two interfaces, both of which contribute to the total IFT of the system.We performed f latsurf calculations to obtain transfer energies from both phases to the surface phase.We calculated the (non-normalized) surface mole fraction, θ′(i), of component i at the interface using where x A (i) represents the mole fraction of component i in phase A and G tot,A→S (i) denotes the f latsurf energy for transferring component i from phase A to the interface between phase A and the surface phase, S.
We then normalized to obtain surface mole fractions, θ(i) Any f latsurf calculation in COSMOtherm requires an IFT value as input because the calculated transfer free energies for bringing molecules to the interface depends on the IFT of the structureless surface.As a starting value, we used 30 mN/m for all calculations.The contribution to the IFT for component i in phase A can be calculated from where A av,A→S (i) represents the expected area of the molecular cross section at the interface between A and S. If the bulk phases and surface phase were in perfect equilibrium, the IFT could equivalently be calculated from Because the composition of the surface phase S is constrained by f latsurf transfer energies from two phases, A and B (eqs 1 and 2), the results using eqs 3 and 4 differ slightly but either expression gives good results for the IFT (Figure S1, Supporting Information).A combined expression, where each contribution has equal weight, turns out to give slightly better agreement with experimental results in order to minimize the mean absolute deviation between the experimental and predicted IFT values.Hence the total contribution to the IFT for phase A is considered to be The contribution for phase B is calculated analogously.The IFT A and IFT B values are generally different from the input IFT value used in the f latsurf calculations.We therefore iterate with the IFT input values for the AS and BS interfaces, which leads to new G tot energies and thus new values for the surface mole fractions.By iterating eqs 1, 2, and 5 with a simple damped linear expression, we obtained self-consistent solutions for IFT A and IFT B .The self-consistency is needed to maintain thermodynamic consistency between the surface phase and the bulk phases in terms of composition, IFT, and transfer free energies.Details are presented in the Supporting Information.
The predicted IFT for the interface, using the self-consistent values for IFT A and IFT B , is thus simply In this study, the IFT value is considered to be converged if the difference of two subsequent iterations is less than 0.02 mN/m.
Error Estimation.Infinite dilution activity coefficients, that is, free energies of transfer of a molecule to another liquid, are typically predicted by COSMO-RS with a root-mean-square deviation of 0.5 ln-units, which corresponds to 1.2 kJ/mol. 37,38he molecules in our evaluation data set have an average surface area of ∼1.5 nm 2 .This results in a value of 1.3 mN/m for COSMO-RS error in free energy per surface area.The predicted IFT in our method results from a sum of IFT A and IFT B , and we should therefore apply a factor of 2 1/2 to the value of 1.3 mN/m to get an error estimate of ∼2 mN/m, which is a lower limit of the expected error of the method.The additional error resulting from the assumptions made in the derivation of our IFT equations cannot be simply estimated.It should be noted that this error can be positive or negative and has to be considered as almost independent of the value of the experimental IFT.Therefore, it can cause large relative errors for the prediction of small IFT values, and it even can result in the prediction of negative IFTs, which obviously should be corrected to an estimate of 0 mN/m afterward.

■ RESULTS
Binary Water/Organic Systems.The predicted IFT values for binary systems are compared with experimental data 5,39 in Figure 2. The agreement is very good for both the TZVPD-fine and TZVP parametrizations, especially considering that no IFT data have been used in the parametrization of COSMO-RS.The points are narrowly scattered around the y = x line, which highlights the significant prediction power of the model with no need for linear regression.Thus, for binary systems, our method gives very good agreement with experimental data.There is no loss of accuracy in absolute terms for low IFT values, which is an improvement over MD methods. 16In particular, the TZVPD-fine performs well for the low IFT values (Figure 2).The results for binary systems are also compiled in Table 1 and show that the TZVPD-fine parametrization performs slightly better than TZVP.The TZVPD-fine parametrization was therefore used for investigating the multicomponent systems.
It is also instructive to compare the errors in predicted IFT for various functional groups, and we have compiled them in Table 2. Some predictions for low IFT of ∼5 mN/m are significantly underestimated (Figure 1), and these come from organic compounds for which the most polar group is the keto group (>CO).Their mean error is −6.6 mN/m.On the other side of the spectrum, there is another class of compounds for which the IFT is overestimated significantly, by 5.9 mN/m.These are molecules for which the most polar group is a phenyl.These have a predicted IFT of about 40 mN/m (Figure 1).In the case where several functional groups are present on a single organic molecule, we found that the most polar group determines the error.This makes sense because the most polar group of the organic compound is the one responsible for interaction with the water phase and the lowering of the IFT.This means that for acetophenone the presence of the keto group determines the sign of the error, not the phenyl group.These systematic errors probably stem from the COSMO-RS parametrization of the interactions.
For the TZVPD-fine level, the COSMO-RS method predicts the following systems to be miscible: water/cyclohexanone, water/2-butanol, and water/methyl acetate.Our method predicts the IFT to be 0 mN/m for these systems (Table 1), and thus, the composition of the surface phase is identical to the bulk phase(s).This shows that our method is thermodynamically consistent in the low IFT limit.
Ternary Systems.A comparison between experimental data 39 and our model for four ternary systems is shown in Figure 3.In all cases, the data are plotted as a function of the mole fraction of the most surface-active component, where partitioning into the aqueous and organic phases is taken into account.The very good agreement between experimental data and our predictions for the IFT as a function of concentration suggests that our method can capture the surface activity and general surfactant behavior for nonionic surfactants.The method is still limited by how well COSMO-RS describes solubilities and partitioning of various components into the two phases.
Few other methods can account for how the partitioning of the surface active compound changes when the concentration is increased and its mole fraction in either of the bulk phases or the surface phase becomes significant.Our method handles all of this automatically.The largest deviation from experimental results, observed in Figure 3, is for the water/propanol/heptane system, where COSMO-RS overestimates the affinity of propanol to water compared to heptane.The partition of acetic acid between water and benzene is also slightly skewed toward the water phase and could result from acetic acid dimerization.This would increase the solubility in the organic phase compared to water and would not be taken into account.In all cases, there is good agreement for a wide range of concentrations, particularly for comparison of IFT with surfactant mole fraction in the organic phase.The surface mole fractions we determine for each system also allow us to quantify surface enrichment of the surfactant molecule (Table S1, Supporting Information) and this gives a picture of the surface composition and how it differs from the bulk phases.We have previously calculated surface enrichment at an oil/water interface for several substituted benzenes using the f latsurf module. 40Our new procedure predicts surface enrichment in a more self-consistent manner thanks to the introduction of a surface phase where the composition is different from the bulk phases.The COSMO-RS interactions between surface-active molecules with a low a All IFT and deviations from experiment are reported in mN/m.The predicted values are plotted as a function of the experimental data in Figure 1.b The single point with negative IFT is not shown in Figure 2. c The TZVPD-fine method predicted these binary systems to be completely miscible and our method predicted the IFT to be 0.These points are excluded from the statistical analysis and from Figure 2. 5.9 a In the case of two functional groups, the most polar one is listed.
bulk phase mole fraction and a high surface mole fraction are taken into account.
Multicomponent System: Mixtures of Dense Nonaqueous Phase Liquids (DNAPL).We obtain a good agreement between experimental data 5 and predicted behavior, even for an organic phase with three or four components (Figure 4).We predict more or less linear dependence on concentration, such as for linear mixing theory.The results again demonstrate the strength and generality of our method.The average relative error is 14% in Figure 4a and 19% in Figure 4b.The accuracy is comparable to what is obtained using empirical methods for these systems. 5The presence of trans-dichloroethylene (trans-DCE) explains one of the largest deviations of the predictions from experimental results (Table 1).Our method can be used to gain understanding and for predicting behavior in organic/ aqueous systems.For example, chlorinated organic contaminants migrate in soil, where the IFT between water and the DNAPL phase is the main control on capillary flow.Changes in composition, partitioning, and temperature can be taken into account in a thermodynamically consistent manner.This would prove useful for being able to predict changes in behavior throughout the contaminated site as the composition of the fluids varies spatially and evolves with time. 41The model would be particularly useful for real groundwater sites, where surfaceactive molecules are also present, to avoid the complications inherent with empirical methods.

■ CONCLUSION AND OUTLOOK
We have developed a method for predicting the interfacial tension between two liquid phases of arbitrary composition using only density functional theory and the implicit solvent model COSMO-RS, which allows partial implicit solvation in two adjacent phases.The implicit solvation method is able to describe liquid/liquid IFT accurately with little computational effort.The method is fast, reliable, and requires no experimental input data.The IFT values predicted by the model agree well with experimental binary, ternary, and multicomponent systems.
The robustness of our method suggests possibilities for predicting interfacial tension in liquid/gas and liquid/solid systems as well.This would allow prediction of the surface tension of a liquid, where a lot of experimental data exists for comparison.It would also open up many new areas of application, including liquid/gas interfacial tension predictions, which are quite important, for example, for CO 2 storage 42,43 and CO 2 -enhanced oil recovery 44 as well as contact angle predictions on solid substrates, which are relevant for separating hydrophobic compounds from mixed phase systems.The possibility of deriving IFT for organic compounds in general makes it possible to predict behavior and design specific surfaces for self-assembly processes that are based on interplay between liquid and solid interfacial tensions.The importance of solid/liquid IFT has also recently been demonstrated for biomineralization rates on polysaccharide templates, where the surface free energy of the polysaccharide/water interface is a determining factor for the nucleation rate of the calcite mineral. 45lthough standard COSMO-RS treatment works well for neutral phases, we are working on generalizing our method to handle ionic species.If ions are present in the system, the surface phase would generally be charged and the degree of charge would depend on the affinity of the ions for the surfaces (eqs 1 and 2).−48 It would make it possible to predict IFT in more complicated systems, including ionic liquids.

■ ASSOCIATED CONTENT
* S Supporting Information Details about the self-consistent calculation for obtaining the IFT and some additional results, a comparison of experimental and predicted temperature dependence, as well as the composition of the surface phase for a few ternary systems.This material is available free of charge via the Internet at http://pubs.acs.org.

Notes
The authors declare the following competing financial interest(s): One of the authors, Andreas Klamt, declares a potential conflict of interest because he is commercially distributing one of the implementations of the COSMO-RS method and COSMOmic from his company COSMOlogic GmbH&CoKG.

■ ACKNOWLEDGMENTS
We are grateful for financial support through the W-EOR Project from the Maersk Oil Research and Technology Centre as well as the Nano-Sand Project in the ExploRe Program with BP Exploration Operating Company, Ltd.The Danish Center for Scientific Computing (DCSC) provided the computing resources.

Figure 1 .
Figure 1.Molecular structure (top), COSMO surfaces (middle), and COSMO images generated with flatsurf (bottom) for water and aniline.In the molecular structure, O is represented by a red sphere, H in white, C in gray, and N in blue.In the COSMO surfaces, red to yellow represents zones where the screening charge is positive, dark to light blue is negative, and green is neutral.The f latsurf COSMO surfaces show the molecule's position with respect to a water/aniline interface.The water phase is represented by dark gray and the aniline phase by light gray.

Figure 2 .
Figure 2. Experimental and predicted interfacial tension for binary water/organic systems.The straight lines in the graphs are the y = x lines.Table1lists the molecules investigated.The mean absolute deviation (MAD) for each method is listed in the graphs.

Figure 3 .
Figure 3. Experimental and predicted IFT for some ternary water/surfactant/organic systems as a function of surfactant concentration.Here, surfactant refers to the most surface active molecule in each ternary mixture.The TZVPD-fine parametrization was used.

Table 2 .
Errors in Binary Organic Phase/Water IFT for Various Functional Groups a