A Novel Approach to Phase Equilibria Predictions Using Ab Initio

Alexander A. Novitskiy , Jie Ke , Victor N. Bagratashvili , and Martyn Poliakoff. The Journal ... Eli Ruckenstein, Ivan L. Shulgin, and Jeffrey L. Til...
0 downloads 0 Views 134KB Size
Download PDF
Ind. Eng. Chem. Res. 1999, 38, 2849-2855

2849

A Novel Approach to Phase Equilibria Predictions Using Ab Initio Methods Amadeu K. Sum and Stanley I. Sandler* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716

Molecular orbital ab initio calculations have been used to compute interaction energies between pairs of molecules in a large molecular cluster. These energies are then used as the interaction energy parameters in the widely used Wilson and UNIQUAC activity coefficient models. Lowpressure vapor-liquid equilibria predictions based on the calculated parameters have been computed for binary systems of water with methanol, ethanol, 1-propanol, 2-propanol, formic acid, acetic acid, acetone, acetonitrile, acetaldehyde, and m-methylformamide. Excellent predictions are obtained with the UNIQUAC model, whereas poor results are found with the Wilson model. In several cases, our predictions are also superior to those obtained from UNIFAC. In addition, using the same parameters and the UNIQUAC model, high-pressure vapor-liquid equilibria predictions were made using the Peng-Robinson-Stryjek-Vera equation of state and the Wong-Sandler mixing rule for methanol, ethanol, 2-propanol, and acetone separately with water. The low- and high-pressure results demonstrate that this unique approach can lead to accurate vapor-liquid equilibrium predictions for hydrogen-bonding mixtures based only on pure-component properties and the structure of the molecules. Introduction Activity coefficients play a central role in phase equilibria calculations. Currently, UNIFAC1 is the best predictive model for activity coefficients; however, it is limited by the availability of parameters that have been regressed from experimental data and by the uncertain accuracy of the method. In addition, as a result of technological advances and environmental concerns, thermodynamicists are faced with the need to predict the phase behavior of more complex and/or hazardous systems for which experimental data may not exist. Parameters for activity coefficient models such as Wilson2 and UNIQUAC3 are usually determined by regression of experimental data. However, in the development of these models, the adjustable parameters are interpreted as interaction energies between molecular species. If these models are theoretically sound and their parameters are indeed true interaction energies, one should be able to directly calculate the energies in these models. Ab initio quantum mechanics methods provide a way to calculate the interaction energies between molecules. Quantum mechanics describes the behavior of systems at the most fundamental level of a particle through electronic and nuclear interactions as governed by the Schro¨dinger equation. Here we present a new and unique approach to phase equilibrium predictions based solely on the use of ab initio quantum mechanics to determine the parameters in activity coefficient models, without the use of experimental data. This is done by using ab initio methods to determine a minimum-energy configuration of a cluster of molecules. Then, from this cluster, the interaction energies of like and unlike pairs of molecules are calculated and used as the energy parameters in the Wilson and UNIQUAC activity coef* To whom correspondence should be addressed. Phone: 302-831-2945. Fax: 302-831-4466. Email: [email protected].

ficient models. Vapor-liquid equilibria (VLE) predictions using this method are presented for several binary aqueous systems at both low and high pressures. Computational Procedure The focus of this work is the use of ab initio methods for the calculation of the interaction energy parameters used in activity coefficient models. A detailed description of the mathematical formalism of quantum mechanics is beyond the scope of this paper. (Excellent texts are available with detailed derivations and discussion of the methods used here.4,5) Instead, we will focus on the development of the approach used for the calculations. Our goal was to develop a systematic calculational procedure so that the interaction energies calculated from ab initio methods were consistent and unambiguous. The following is an outline of the steps used for the computation of the interaction energies: (1) A cluster of molecules is constructed to represent a dense fluid. A cluster size composed of eight molecules (four of each kind) has generally been found to give the best results. (2) The cluster energy is initially minimized using the PM3 semiempirical method.6 This initial energy minimization of the cluster is computationally inexpensive and is a good way to obtain a starting geometry for the following step. (3) Another energy minimization is performed on the cluster thus obtained using the Hartree-Fock (HF) method with a 6-31G** basis set (HF/6-31G**). This minimization is more easily done by gradually increasing the basis set, that is, going from 6-31G to 6-31G* and then 6-31G**. (4) Directly interacting molecular pairs are selected from this optimized cluster, that is, all like and unlike pairs that are in close proximity, and their separation distances and relative orientations recorded. Pairs that

10.1021/ie9900263 CCC: $18.00 © 1999 American Chemical Society Published on Web 05/25/1999

2850 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999

are directly interacting are usually hydrogen bonded or in close proximity. (5) The interaction energy of each molecular pair is computed using the HF method with the larger 6-311++G(3d,2p) basis set (HF/6-311++G(3d,2p)) at the separation and orientation obtained in the previous step. The pair interaction energy is then computed from

Eint AB ) EAB{AB} - EA{AB} - EB{AB}

(1)

where E is the energy of molecule(s) A and/or B and {AB} refers to the combined basis set for both A and B molecules (dimer basis set). (6) Finally, as there will usually be several of the same interacting molecular pair (e.g., H2O-H2O), energies of the sets of the same molecular pairs are linearly averaged to obtain the energy parameters to be used in the activity coefficient models. The steps outlined above have proven to give consistent and good results, as will be shown in the next section. However, it is also important to mention some of the issues associated with the procedure just described. The cluster size composed of eight molecules (four of each kind) was selected after several attempts were made with smaller clusters (e.g., two and four). It was found that with the smaller clusters we were required to repeat the above procedure for several different initial configurations because as the results were strongly dependent on the initial arrangement of the molecules. However, with eight molecules, the results obtained were more consistent and less dependent on the initial configuration of the cluster. Also, a cluster size of eight molecules was chosen as a compromise between computational cost and a reasonable representation of the phase space of a dense, liquidlike fluid. The PM3 semiempirical method was chosen for the initial optimization of the cluster because it qualitatively describes hydrogen bonding, and all of the systems studied thus far include water. As mentioned above, this initial energy minimization was done to obtain a better starting point for the subsequent, more rigorous minimization. One could also initially minimize the cluster energy with some other method, such as molecular mechanics or another appropriate semiempirical method. The energy minimization step at the HF level is the most time-consuming portion of the calculations. For systems with large molecules, the optimization may require several days to over a week on a multiprocessor supercomputer. However, the better the initial geometry obtained, the quicker the optimization. Hartree-Fock theory is known to result in reasonably accurate geometries of molecules, and it is the least expensive method except for density-functional theory7 (DFT), which was not used for our calculations because that method is best suited for studying strongly interacting systems, not weak intermolecular interactions (i.e., hydrogen bonding and van der Waals forces) as in the systems under study here. The use of more sophisticated methods (e.g., Møller-Plesset perturbation theory4,5) could not be justified by the enormous additional computational cost for slight improvement that would be obtained in the final optimized geometry. The basis set 6-31G** for the final geometry optimization was selected in order to properly model the intermolecular interactions, both quantitatively and qualitatively; smaller basis sets did not give reasonable results and a larger basis set would

make the calculations too expensive computationally. The cluster optimized at HF/6-31G** is taken to represent a sample of the phase space of a dense fluid; however, we cannot be assured that the final configuration obtained is the absolute minimum-energy configuration for that particular cluster (for systems as complex as the ones under study, many local minima can be found. Once the energy of the cluster of molecules was minimized, the coordinates of the molecules were frozen in space. Interaction energies were calculated only for those pairs of molecules within close proximity. The interaction energies were calculated by the so-called supermolecular approach.8 To obtain the interaction energy between a pair of molecules, three energy calculations were required: (i) total energy of the AB pair with the {AB} basis set, (ii) total energy of molecule A with the {AB} basis set, and (iii) total energy of molecule B with the {AB} basis set. Steps ii and iii are calculated with ghost atoms of the other molecule, thus the use of the {AB} basis set. The interaction energy was then computed as given by eq 1. The {AB} basis set was used for all calculations in order to eliminate the basis set superposition error9 (BSSE), that is, the nonphysical error in the calculated energy associated with the limited size of the basis set describing the wave function of the electrons. The energies were calculated using the HF method with the 6-311++G(3d,2p) basis set. Hartree-Fock theory does not account for electron correlation and therefore weak van der Waals forces. However, because the final quantity of interest is the difference of interaction energies, we assumed electron correlation would largely cancel out so that it need not explicitly be considered. Methods accounting for electron correlation, such as DFT, were not used for the reason discussed earlier. Also, Møller-Plesset perturbation theory was not used because of the computational cost for the large systems. The large basis set was chosen in order to minimize the effect of the BSSE (other, smaller basis sets gave similar total interaction energies, but the BSSE was a large fraction of these energies). For the 6-311++G(3d,2p) basis set, the BSSE was always just a few percent of the interaction energy. All of the ab initio computations were performed using the Gaussian 9410 program running on multiprocessor supercomputers (Cray J-90 and SGI Origin2000). The two activity coefficient models considered in this study are the Wilson and UNIQUAC. Both models have two adjustable parameters that are interpreted as the difference of interaction energies between the like and unlike molecular pairs. For the Wilson model these int parameters are ∆λij ) Eint ij - Eii , and for UNIQUAC int int they are ∆uij ) Eij - Ejj . In some forms of the UNIQUAC model, there is a factor of z/2 multiplying the interaction energies, and in other cases, this factor is unity; in this work it was taken to be unity. With these energy parameters and the pure-component vapor pressures, we were able to predict the phase behavior of several binary systems of water separately with methanol, ethanol, 1-propanol, 2-propanol, formic acid, acetic acid, acetone, acetonitrile, acetaldehyde, and m-methylformamide. The low-pressure VLE predictions were performed using the usual γ-φ method assuming an ideal vapor phase, except for formic and acetic acid systems where dimerization was taken into account, as described in

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2851 Table 1. Interaction Energies from Ab Initio Calculations Using HF/6-311++G(3d,2p)//HF/6-31G** system

E11 (kJ/mol)

E22 (kJ/mol)

E12 (kJ/mol)

∆u12a (J/mol)

∆u21a (J/mol)

methanol (1) + water (2) ethanol (1) + water (2) 1-propanol (1) + water (2) 2-propanol (1) + water (2) formic acid (1) + water (2) acetic acid (1) + water (2) acetone (1) + water (2) acetonitrile (1) + water (2) acetaldehyde (1) + water (2) m-methylformamide (1) + water (2)

-12.300 -12.426 -12.678 -10.709 -22.813 -13.253 -8.678 -10.827 -6.884 -17.509

-11.852 -11.294 -12.110 -12.784 -11.582 -9.578 -12.094 -11.490 -12.692 -10.862

-12.247 -11.332 -11.462 -10.944 -15.573 -11.001 -9.138 -9.900 -8.231 -13.987

-395.05 -37.36 647.26 1840.21 -3990.53 -1423.44 2956.54 1589.88 4461.15 -3124.90

52.97 1093.91 1215.45 -234.51 7239.87 2252.00 -460.11 926.97 -1346.24 3522.72

a

Parameters for the UNIQUAC model.

Figure 1. VLE diagram for methanol (1) + water (2).

Figure 2. VLE diagram for ethanol (1) + water (2).

Gmehling et al.11 Because our interest is also in the description of VLE over large ranges of temperature and pressure, we have also used the results obtained to make VLE predictions at high pressure using the PengRobinson-Stryjek-Vera (PRSV) equation of state12 and the Wong-Sandler mixing rule.13 Results Table 1 summarizes the energies calculated from the ab initio computations as well as the interaction parameters used for the UNIQUAC activity coefficient model. Other relevant parameters (Vm for the Wilson model and ri and qi for the UNIQUAC model) and the constants in the Antoine equation were obtained from Gmehling et al.11 In addition, for the high-pressure predictions, critical properties were obtained from Reid et al.,14 and the kij’s were determined by matching the excess free energy calculated from the equation of state with that of UNIQUAC using the quantum-mechanically calculated interaction parameters. Figures 1-10 show the low-pressure VLE diagrams for the binary aqueous mixtures of methanol, ethanol, 1-propanol, 2-propanol, formic acid, acetic acid, acetone, acetonitrile, acetaldehyde, and m-methylformamide. All of the figures are P-x-y diagrams (except Figure 10 for mmethylformamide + water which is a T-x-y diagram) plotting the results of the VLE calculations using the interaction parameters calculated from this work. Our results are compared with experimental data (from compilations of Gmehling et al.11) and predictions based on the UNIFAC group contribution model. The results

Figure 3. VLE diagram for 1-propanol (1) + water (2).

for the binary mixtures of methanol, ethanol, formic acid, and acetic acid, separately with water, were reported previously at different temperatures.15 Highpressure predictions were only done for those systems for which experimental data were available for comparison. Data for the methanol + water and acetone + water systems were obtained from Griswold and Wong16 and those for the ethanol + water and 2-propanol + water systems from Barr-David and Dodge.17

2852 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999

Figure 4. VLE diagram for 2-propanol (1) + water (2).

Figure 5. VLE diagram for formic acid (1) + water (2).

Discussion The results presented in this paper only deal with binary aqueous systems. We first explored the calculational methodology with aqueous systems because water is a small molecule (has few electrons) and its interactions with the other molecules is mainly through hydrogen bonding. In addition, to use the Hartree-Fock method to calculate interactions energies, we needed to study systems for which electron correlation resulting in van der Waals or dispersion energies could be neglected (Hartree-Fock theory does not implicitly account for electron correlation). This latter restriction is a reasonable approximation for hydrogen-bonding systems, which are predominantly characterized by exchange energies as the main contribution to the interaction energy between molecules (van der Waals or dispersion energies are negligible). We found that as the nonaqueous molecule increased in size and flexibility, it was more difficult to obtain a minimum-energy cluster, mainly because of the large number of geometric degrees of freedom. Another observation from our computations is that it is important to sample the phase space of many possible interactions

Figure 6. VLE diagram for acetic acid (1) + water (2).

Figure 7. VLE diagram for acetone (1) + water (2).

between the molecules in order to obtain the correct interaction energies, especially for the systems with large nonaqueous molecules. As shown in Figures 1-10, the low-pressure predictions based on the UNIQUAC model using the calculated interaction parameters are very good and in some cases superior to predictions based on UNIFAC. However, predictions with the Wilson model, using the same interaction energies, did not give good results for any of the systems studied. In the development of both the UNIQUAC and Wilson models, the adjustable parameters are interpreted as interactions energies. Given the results presented here, it appears the UNIQUAC model has a reasonable theoretical basis because the energies that appear in this model are the same as those calculated from quantum mechanics. In contrast, the parameters in the Wilson model must be treated as completely adjustable to obtain good agreement with experimental data. However, we should point out that the Wilson model is frequently a good one when its parameters are fitted to experimental data. Even though the results we have presented are encouraging, it is important to note the approximations and assumptions made in the calculations:

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2853

Figure 8. VLE diagram for acetonitrile (1) + water (2).

Figure 9. VLE diagram for acetaldehyde (1) + water (2).

1. Electron correlation has been neglected by using HF theory in the ab initio calculations, for both the energy minimization and interaction energy calculations. As mentioned above, by using a more sophisticated method (e.g., Møller-Plesset perturbation theory) the energy minimization calculation would probably add little to the final configuration for hydrogen-bonding fluids at a large additional computational expense. Perhaps the same would not be true for the interaction energies, but because of cancellation of errors when taking the difference of interaction energies, the results may not be much affected by the contribution of electron correlation. Also, because hydrogen-bonding and Coulombic interactions are the predominant interactions for the systems studied, the contribution of electron correlation resulting from the van der Waals or dispersion energy contributing to the interaction energy is ignored in our calculational scheme and assumed to be negligible. 2. A cluster size of eight molecules optimized at HF/ 6-31G**, at 0 K, and in vacuo is a good representation of the system. This assumption implies that the state of a fluid at temperatures other than 0 K (as long as it is below the critical point of the fluids involved) can be

Figure 10. VLE diagram for m-methylformamide (1) + water (2).

approximated by the equilibrium configuration of a sample cluster at 0 K. 3. The calculated interaction energies at 0 K are applicable as parameters in activity coefficient models to predict the phase behavior of systems at temperatures other than 0 K. This, along with the previous assumption, is perhaps the most difficult to justify. It is known that the intermolecular structure of a fluid is a function of the temperature and density; however, nowhere in the ab initio calculations is a temperature or density specified. 4. The energy parameters in the UNIQUAC and Wilson models are the quantum-mechanically calculated interaction energies obtained from the minimum-energy configurations. This assumption presumes that the sample clusters describe the most likely states in the dense fluid and also that the parameters in these semiempirical models are true interaction energies. The high-pressure predictions were based on the UNIQUAC model with the previously determined parameters and the Wong-Sandler mixing rule. No further ab initio calculations were performed, and the interaction parameters were the same as those used in the low-pressure predictions. As seen from Figures 1114, the predictions are remarkably good. Both the lowand high-pressure results are true predictions in the sense that the only input required to determine the interaction parameters was the atomic connectivity of the molecules. Another interesting and perhaps surprising result was the range of conditions for which the calculated parameters were applicable. On the basis of a single set of parameters for the UNIQUAC model, when used in an equation of state with an appropriate mixing rule, accurate VLE predictions were obtained over a large range of temperatures and pressures. It is known that the parameters for the UNIQUAC model are correlated, and a nonunique set is obtained when correlating a particular set of data. However, the method here results in a unique set of parameters for this model. So far a limited subset of polar/polar systems was studied. It is desirable to extend this approach to other systems including nonpolar/nonpolar, polar/polar, and nonpolar/polar interactions in order to test the general applicability of this method. Currently, our computa-

2854 Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999

Figure 11. High-pressure VLE prediction for the methanol (1) + water (2) system using the quantum-mechanically determined parameters, PRSV EOS, and Wong-Sandler mixing rule.

Figure 13. High-pressure VLE prediction for the 2-propanol (1) + water (2) system using the quantum-mechanically determined parameters, PRSV EOS, and Wong-Sandler mixing rule.

Figure 12. High-pressure VLE prediction for the ethanol (1) + water (2) system using the quantum-mechanically determined parameters, PRSV EOS, and Wong-Sandler mixing rule.

Figure 14. High-pressure VLE prediction for the acetone (1) + water (2) system using the quantum-mechanically determined parameters, PRSV EOS, and Wong-Sandler mixing rule.

tional methodology is being tested on other larger systems, both aqueous and nonaqueous. However, so far there has been limited success in obtaining good VLE equilibria predictions, one reason being that, for larger molecules and weakly interacting systems, calculation of interaction energies at the HF method level does not account for weaker interactions (i.e., van der Waals and dispersion energies). Using methods that properly account for electron correlation in large systems is very demanding of currently available computational resources, which at present is a limitation. For example, in the next least computationally intensive method to calculate the interaction energies that includes electron

correlations, MP2, the CPU and storage requirements scale with approximately the fourth power of the total number of basis functions and number of electrons. Conclusions Interaction energies between molecules have been calculated from ab initio methods. These energies, in turn, have been used as the interaction energy parameters in the Wilson and UNIQUAC activity coefficient models. VLE predictions, both at low and high pressures, with the UNIQUAC model gave excellent results

Ind. Eng. Chem. Res., Vol. 38, No. 7, 1999 2855

for the binary aqueous systems considered; however, the same was not true with the Wilson model. This result suggests that the UNIQUAC model has a sounder theoretical basis than the Wilson model. Acknowledgment We are very grateful to the University of Delaware for providing a Presidential Fellowship to A.K.S. and the computational resources required for this work. This research also used computational resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Energy Research of the U.S. Department of Energy under Contract No. DEAC03-76SF00098. Also, we acknowledge the funding of this project by the Department of Energy (DE-FG0285ER13436) and National Science Foundation (CTS9521406). Literature Cited (1) Fredenslund, Aa.; Gmehling, J.; Rasmussen, P. VaporLiquid Equilibria Using UNIFAC; Elsevier: Amsterdam, The Netherlands, 1977. (2) Wilson, G. M. Vapor-Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing. J. Am. Chem. Soc. 1964, 86, 168. (3) Abrams, D. S.; Prausnitz, J. M. Statistical Thermodynamics of Liquid Mixtures: A New Expression for the Excess Gibbs Energy of Partly or Completely Miscible Systems. AIChE J. 1975, 21, 116. (4) Levine, I. N. Quantum Chemistry, 4th ed.; Prentice-Hall: Englewood Cliffs, NJ, 1991. (5) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry. Introduction to Advanced Electronic Structure Theory; Dover Publications: New York, 1996. (6) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods. I. Method. J. Comput. Chem. 1989, 10, 209. (7) Parr, R. G.; Yang, W. Density-Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989.

(8) Chalasinski, G.; Gutowski, M. Weak Interactions between Small Systems. Models for Studying the Nature of Intermolecular Forces and Challenging Problems for Ab Initio Calculations. Chem. Rev. 1988, 88, 943. (9) van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Lenthe, J. H. State of the Art in Counterpoise Theory. Chem. Rev. 1994, 94, 1873. (10) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T. A.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Peng, C. Y.; Ayala, P. A.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revisions D.1, D.3, D.4; Gaussian, Inc.: Pittsburgh, PA, 1995. (11) Gmehling, J.; Onken, U.; Arlt, W. Vapor-Liquid Equilibrium Data Collection; DECHEMA: Frankfurt, Germany, 1977; and onward. (12) Stryjek, R.; Vera, J. H. An Improved Peng-Robinson Equation of State for Accurate Vapor-Liquid Equilibrium Calculations. Can. J. Chem. Eng. 1986, 64, 334. (13) Wong, D. S. H.; Sandler, S. I. A Theoretically Correct Mixing Rule for Cubic Equations of State. AIChE J. 1992, 38, 671. (14) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1987. (15) Sum, A. K.; Sandler, S. I. Use of Ab Initio Methods to Make Phase Equilibria Predictions Using Activity Coefficient Models. Eighth International Conference on Properties and Phase Equilibria for Product and Process Design, Noordwijkerhout, The Netherlands, April 1998. (16) Griswold, J.; Wong, S. Y. Phase-Equilibria of the AcetoneMethanol-Water System from 100 °C into the Critical Region. Chem. Eng. Prog. Symp. Ser. 1952, 48, 18. (17) Barr-David, F.; Dodge, B. G. Vapor-Liquid Equilibrium at High Pressures. The Systems Ethanol-Water and 2-PropanolWater. J. Chem. Eng. Data 1959, 4, 107.

Received for review January 12, 1999 Revised manuscript received April 15, 1999 Accepted April 21, 1999 IE9900263