J. Phys. Chem. B 1999, 103, 11191-11195
11191
Simulating Retention in Gas-Liquid Chromatography Marcus G. Martin and J. Ilja Siepmann* Departments of Chemistry and of Chemical Engineering and Materials Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455-0431
Mark R. Schure Theoretical Separation Science Laboratory, Rohm and Haas Company,727 Norristown Road, Spring House, PennsylVania 19477 ReceiVed: September 14, 1999
Accurate predictions of retention times, retention indices, and partition constants are a long sought-after goal for theoretical studies in chromatography. Configurational-bias Monte Carlo (CBMC) simulations in the Gibbs ensemble using the transferable potentials for phase equilibria-united atom (TraPPE-UA) force field have been carried out to obtain a microscopic picture of the partitioning of 10 alkane isomers between a helium vapor phase and a squalane liquid phase, a prototypical gas-liquid chromatography system. The alkane solutes include some topological isomers that differ only in the arrangement of their building blocks (e.g., 2,5dimethylhexane and 3,4-dimethylhexane), for which the prediction of the retention order is particularly difficult. The Kovats retention indices, a measure of the relative retention times, are calculated directly from the partition constants and are in good agreement with experimental values. The calculated Gibbs free energies of transfer for the normal alkanes conform to Martin’s equation which is the basis of linear free energy relationships used in many process modeling packages. Analysis of radial distribution functions and the corresponding energy integrals does not yield evidence for specific retention structures and shows that the internal energy of solvation is not the main driving force for the separation of topological isomers in this system.
1. Introduction Although recent advances in computational chemistry and molecular modeling have improved our understanding of complex chemical systems, little attention has been focused on chromatography, let alone calculations of retention properties.1,2 Gas-liquid chromatography (GLC) is one of the principal methods for the analysis and separation of organic and pharmaceutical molecules. The underlying principles of chromatographic separation are inherently complex, being dictated by the interplay of the sample molecules with the stationary and mobile phases. Predicting the retention characteristics of a solute molecule given only its structure and the experimental chromatographic conditions remains one of the grand challenges in separation science. Our work shows that molecular simulations using advanced techniques (configurational-bias Monte Carlo in the Gibbs ensemble) and transferable force fields (TraPPE-UA force field) can be carried out to predict retention times and to provide microscopic insight into chromatographic separation processes. The application to chromatography poses an extreme challenge on the efficiency of the sampling algorithms and the accuracy of the force fields, because free energies of transfer have to be predicted correctly to within a few tenths of a kJ/mol.3 Here we describe an investigation of the partitioning of linear and branched alkanes with 5-8 carbon atoms between helium and squalane (2,6,10,15,19,23-hexamethyltetracosane), a liquidphase often used in GLC4,5 and the reference material for the Rohrschneider-McReynolds scheme of liquid-phase characterization.6 Two sets of topological isomers, including 2-methylpentane (2-MP), 3-methylpentane (3-MP), 2,5-dimethylhexane
(2,5-DMH), and 3,4-dimethylhexane (3,4-DMH), were selected for this study because the prediction of retention times for topological isomers that are made from the same set of building blocks (methyl, methylene, and methine groups) is very difficult and standard group additivity methods would treat these isomers as identical species. In addition, calculations were also performed for 2,2-dimethylpentane (2,2-DMP, containing a quaternary carbon), 3-ethylpentane (3-EP, star-shaped with ethyl branches), and four linear alkanes: n-pentane (PEN), n-hexane (HEX), n-heptane (HEP), and n-octane (OCT). 2. Simulation Methods The Gibbs ensemble utilizes two separate simulation boxes which are in thermodynamic contact but do not have an explicit interface. In the isobaric version of the Gibbs ensemble, mechanical equilibrium for each phase is established by volume exchanges with an external pressure bath. Coupled-decoupled configurational-bias Monte Carlo (CD-CBMC) is used to swap molecules between the two phases, thereby equalizing the chemical potential of each species.8,9 Thermal equilibrium in each phase is reached via translational, rotational, and CDCBMC conformational moves. As a result, for a given state point the properties of the coexisting phases, such as the partitioning of solute molecules, can be determined directly from a single simulation. In the special case of simulations for GLC systems, particle swap moves have to be performed only for the solutes and the carrier gas. Since the liquid phase in GLC is only used over a temperature range where its vapor pressure is negligible, there is no need to sample the partitioning of squalane.
10.1021/jp9932822 CCC: $18.00 © 1999 American Chemical Society Published on Web 11/24/1999
11192 J. Phys. Chem. B, Vol. 103, No. 50, 1999
Martin et al.
TABLE 1: Summary of Simulation Results for the Helium-Squalane Gas-Liquid Chromatography System (T ) 343 K, p ) 101.3 kPa): Partition Constants, Gibbs Free Energies of Transfer (in kJ/mol), and Kovats Retention Indices system
molecule
K
∆G
PEN/HEX PEN/HEX PEN/HEX PEN/HEX PEN/HEX HEX/HEP HEX/HEP HEX/HEP HEX/HEP HEX/HEP HEP/OCT HEP/OCT HEP/OCT HEP/OCT HEP/OCT
He PEN 2-MP 3-MP HEX He HEX 2,2-DMP 3-EP HEP He HEP 2,5-DMH 3,4-DMH OCT
0.03776 26.68 50.81.8 58.32.0 60.81.9 0.03825 61.22.4 85.74.2 1426 1487 0.03818 1396 19210 31118 31215
+9.355 -9.359 -11.2010 -11.5910 -11.719 +9.315 -11.7311 -12.6914 -14.1311 -14.2414 +9.325 -14.0713 -14.9916 -16.3717 -16.3814
a
I 500 5787 5958 600 600 6388 6969 700 700 7409 79912 800
The subscripts give the statistical uncertainties in the last digit(s).
Figure 2. Gibbs free energies of transfer between helium and squalane phases for normal alkanes as a function of the number of carbon atoms (T ) 343 K, p ) 101.3 kPa). The solid line shows an unweighted linear fit to the data.
methine-methine pairs are 4 and 1, respectively. This ratio of 9:4:1 is reflected in the well depths for the TraPPE-UA model. The Lennard-Jones parameters for helium were obtained from a fit to its critical constants (He/kB ) 4 K, σHe ) 3.11 Å). Lorentz-Berthelot combining rules were used for all unlike interactions.13 A spherical potential truncation at 14 Å and analytical tail corrections were employed for the Lennard-Jones interactions.13 The total lengths of the production periods were 7.8, 4.8, and 4.8 × 105 Monte Carlo cycles for the PEN/HEX, HEX/HEP, and HEP/OCT systems, respectively (one Monte Carlo cycle consists of N trial moves with N being the total number of molecules). 103 Monte Carlo cycles for the PEN/ HEX system take approximately 15 h of CPU time on an Intel Pentium II processor. 3. Results and Discussion
Figure 1. Predicted Kovats retention indices of six branched alkanes (2-MP, square; 3-MP, diamond; 2,2-DMP, circle; 3-EP, cross; 2,5DMH, triangle up; 3,4-DMH, triangle down) for a helium/squalane GLC system (T ) 343 K, p ) 101.3 kPa). The experimental data were taken from ref 4.
Three groups of six independent simulations were carried out where each group sampled the partitioning of a group of four solutes: (i) system PEN/HEX, PEN/2-MP/3-MP/HEX; (ii) system HEX/HEP, HEX/2,2-DMP/3-EP/HEP; and (iii) system HEP/OCT, HEP/2,5-DMH/3,4-DMH/OCT. The simulation systems consisted of 96 squalane, a total of 8 solute (2 each) molecules, and 200 (for PEN/HEX and HEX/HEP) or 500 (for HEP/OCT) helium atoms.10 The TraPPE united-atom alkane force field,9 which has been derived from empirical fitting to vapor-liquid-phase equilibria of pure alkanes, was used for squalane and all solute molecules. This force field employs Lennard-Jones 12-6 potentials11 for the nonbonded interactions of united atoms (that is, methyl, methylene, and methine groups are considered as single interaction centers) and the Lennard-Jones well depth and size σ parameters are as follows: CH3/kB ) 98 K, σCH3 ) 3.75 Å; CH2/kB ) 46 K, σCH2 ) 3.95 Å, CH/kB ) 10 K, σCH ) 4.68 Å; C/kB ) 0.5 K, σC ) 6.4 Å (where kB is Boltzmann’s constant). The large differences in the well depths of the methyl, methylene, and methine pseudo-atoms can be well rationalized because the valence electrons forming the C-H bonds contribute most to the molecular polarizability of alkanes, while the polarizabilities of the C-C bond valence electrons and 1s core electrons on the carbons contribute only a small amount.12 The interaction of a methyl-methyl pair involves nine individual interactions of C-H with C-H bonds, while the corresponding numbers of individual interactions for methylene-methylene and
Gibbs Free Energies of Transfer. Simulations in the Gibbs ensemble7 offer a very elegant and accurate approach to the direct determination of Gibbs free energies of transfer via the corresponding number densities14,15
∆G*S ) -RT ln KS ) -RT ln
() FRS FβS
(1)
where KS is the partition constant, and FRS and FRS are the ensemble averaged number densities of species S in the two phases R and β at equilibrium. The calculated partition constants and Gibbs free energies of transfer for helium and the 10 alkane solutes are listed in Table 1. As expected, the differences in the free energies of transfer between the alkane isomers are relatively small. It is important to note that the free energies of transfer of HEX, HEP, and helium, which were present in more than one system, agree very well between the different simulations (see Table 1). Kovats Retention Indices. Absolute Gibbs free energies of transfer are rarely measured for most chromatographic experiments because the ratio of the mobile to stationary phase volumes must be known for the calculation of equilibrium constants. This phase ratio differs from column to column and often is difficult to measure accurately. Thus it has proven more useful to use the concept of a retention index to give the relative retention time of a given solute.16 The Kovats retention index I of a solute x can be calculated directly from the partitioning using
Ix ) 100n + 100
[
log(Kx/Kn)
log(Kn+1/Kn)
]
(2)
Simulating Retention in Gas-Liquid Chromatography
J. Phys. Chem. B, Vol. 103, No. 50, 1999 11193
Figure 3. Stereoview of the helium-squalane gas-liquid chromatography system HEP/OCT. The helium atoms are depicted as red dots, the backbones of the squalane chains are shown in black, and the alkane solutes are highlighted by thicker lines and coloring: blue for HEP, orange for 2,5-DMH, purple for 3,4-DMH, and green for OCT. To allow for easy viewing, the squalane liquid phase and helium vapor phase are rescaled to similar sizes, but it should be kept in mind that the linear dimension of the periodic cells (outlined in blue) are about 44 and 290 Å, respectively.
where Kx, Kn, and Kn+1 are the partition constants of the solute in question, the highest normal alkane (having n carbon atoms) that elutes prior to the solute, and the lowest normal alkane (having n+1 carbon atoms) that elutes after the solute, respectively. That is, the Kovats index of an n-alkane is 100n, and an index of, for example, 570 means that this solute would elute together with a fictitious normal alkane having 5.70 carbon atoms. The partition constant K and the retention time tr of a solute are related as follows
φK )
tr -1 t0
(3)
where φ and t0 are the phase ratio and the void time (elution time of an unretained compound), respectively. The calculated Kovats retention indices are compared to their experimental counterparts in Figure 1. It is very encouraging that the elution order is predicted correctly for all six branched alkanes and that I values can be predicted with a statistical precision (error of the mean) of about 10 Kovats units (see Table 1). The I values for 2-MP, 3-MP, 2,2-DMP, 3-EP, and 2,5DMH are in satisfactory agreement with experiment, but they are all too high with a root-mean-square error of 10 Kovats units. However, the calculated retention index for 3,4-DMH is
markedly too high. This problem is most likely caused by shortcomings of the TraPPE-UA force field. In particular, the calculated standard boiling points for 3,4-DMH and n-octane are too close (382 and 386 K) compared to their experimental counterparts (391 and 399 K).9 With regard to the retention order, it might be surprising that 3,4-DMH with its bulky center elutes after 2,5-DMH, but analysis of the microscopic structures points to greater flexibility of 3,4-DMH which allows it to pack more efficiently (see below). Martin Equation. The Martin equation predicts that the Gibbs free energies of transfer (and also the net retention times) are a linear function of the number of carbon atoms n in any homologous series
∆G* n ) A + Bn
(4)
It is widely accepted and based on many experimental observations. However, using theoretical arguments it has been pointed out that the Martin equation should only be approximately correct because of size dependent entropic contributions to the free energy.19-23 Our simulation results for ∆G of the n-alkanes indicate linearity at least within the accuracy of the simulations (see Figure 2) and yield a methylene group increment of 2.39 ( 0.10 kJ/mol.
11194 J. Phys. Chem. B, Vol. 103, No. 50, 1999
Martin et al.
Figure 5. Energy integrals for the intermolecular interactions of methyl, methylene, and methine united atoms averaged over all squalane united atoms. Solid, dashed, dotted, and dashed-dotted lines depict the integrals calculated for united atoms belonging to HEP, 2,5-DMH, 3,4DMH, and OCT, respectively.
Figure 4. Liquid-phase radial distribution functions for squalanes solute pairs for the 2,5-DMH (solid line) and 3,4-DMH (dotted line) solutes in system HEP/OCT. Radial distribution functions are shown for all nine combinations of united-atom pairs and are shifted with respect to each other by 2 units.
Microscopic Analysis. A snapshot of the coexisting vapor and liquid phases for the HEP/OCT system is shown in Figure 3. This snapshot is “typical” in two important ways. First, there are equal numbers of solute molecules in the vapor and liquid phases (see ref 10). Second, in neither the vapor nor liquid phase is clustering of the solute molecules apparent. Ensemble averaged over the four types of solute molecules in system HEP/ OCT, the mean number of solute molecules within 14 Å (measured between centers of mass) of another solute molecule is 0.0011 ( 0.0007 and 0.17 ( 0.01 in the vapor and liquid phases, respectively. Furthermore, in the vapor phase the mean number of helium atoms surrounding a solute molecule within 14 Å is only 0.22. To obtain additional microscopic information about the liquid phase, radial distribution functions were calculated for the different combinations of united atoms belonging to the solute and solvent molecules. The liquid-phase radial distribution functions for the two topological octane isomers are shown in Figure 4. Some of the pair distributions are very similar for 2,5-DMH and 3,4-DMH. However, differences are more apparent for the distributions of solute methylene and methine groups with the interior segments of squalane (methylene and methine). Here the first peaks belonging to 2,5-DMH are shifted to shorter separations than is the case for 3,4-DMH. Evidently, the bulky center of 3,4-DMH prevents a close approach. Nevertheless, we should keep in mind that 3,4-DMH partitions more strongly into the liquid phase than 2,5-DMH. Since it is not straightforward to relate the information obtained from the radial distribution functions to differences in the retention (partitioning) process, we have also evaluated the corresponding energy integrals24
〈Uij(r)〉 )
4πNiNj 2V
∫0r uij(r′)gij(r′) r′2 dr′
(5)
where Uij(r) is the internal energy arising from interactions of species i with all neighbors of type j within a distance of less than r. Ni, Nj, and V are the total numbers of species i and j in
the simulation box of volume V, and uij(r′) and gij(r′) are the potential energy and the radial distribution function of pair i and j at separation r′. The energy integrals for the three different types of united atoms present in the solute molecules (system HEP/OCT) and summed over all squalane united atoms are shown in Figure 5. The first observation is that the long-range values of the integrals are clearly separated between the three types of solute groups (methyl, methylene, and methine) and that the difference between methyl and methylene groups is similar to the difference between methylene and methine groups. While there appears to be a noticeable difference between methyl groups belonging to linear or branched isomers, there is little difference based on the chain length (for HEP and OCT) and on the position of the methyl groups (for 2,5-DMH and 3,4-DMH). The difference between methyl groups in linear or in branched isomers is caused by the more compact structure of the branched alkanes which reduces the space available for intermolecular contacts and leads to smaller numbers of intermolecular neighbors. Summing the energy integrals over all united atoms of a given molecule (and adding the tail corrections for contributions beyond the potential truncation) yields the total interaction energy for this molecule with the entire solvent phase. The difference of the summed interaction energies for the liquid and vapor phases should be a good estimate for the internal energy change of transfer (solvation) assuming that the internal energy of the liquid solvent changes insignificantly upon addition of the solute.15 The calculated values for the internal energy change of transfer for the HEP/OCT system are -16.0, -17.1, -17.2, and -17.9 kJ/mol for HEP, 2,5-DMH, 3,4-DMH, and OCT, respectively. While the difference between HEP and OCT is of similar magnitude as the corresponding difference in the Gibbs free energies of transfer (see Table 1), there is little difference in the internal energy change of transfer between 2,5-DMH and 3,4-DMH (similar observations can be made for 2-MP and 3-MP in system PEN/HEX). At the standard pressure used in our simulations, the differences in the volume work of solvation should be very small,15 thus one might expect similarly small differences between internal energy changes and enthalpies of transfer. Thus our results indicate that the internal energy change (or enthalpy) of transfer plays a major role for the normal alkanes, but this internal energy change is not the driving force for the separation of topological isomers. Analysis of the conformational distributions of the topological octane isomers points to greater flexibility of 3,4-DMH. The central dihedral angle of 3,4-DMH is found with roughly equal probability in its trans state or any gauche state, whereas the central dihedral angle of 2,5-DMH is found with a probability
Simulating Retention in Gas-Liquid Chromatography of about 95% in its trans state. Thus the two bulky ends in 2,5DMH might enforce extended conformations that do not pack efficiently with the squalane phase or itself.9 4. Conclusions Using configurational-bias Monte Carlo simulations in the Gibbs ensemble and the TraPPE force field we are able to quantitatively predict the relative retention times of 10 alkane isomers in a helium-squalane gas-liquid chromatographic system. The simulations are able to distinguish isomers solely on the topology of otherwise identical building blocks. The Gibbs free energies of transfer for the normal alkanes are well described by a linear regression, supporting Martin’s empirical law and showing that size dependent contributions to the free energy are not significant for the molecules studied here. The analysis of radial distribution functions and the corresponding energy integrals shows only relatively small differences for individual united atoms belonging to different solutes. However, while the internal energy change (or enthalpy) of transfer contributes substantially to the differences in Gibbs free energies of transfer of the normal alkanes, its role in the separation of topological isomers is much less evident. Finally, the simulation techniques employed here are not limited to gas-liquid chromatography but can be applied to many other forms of chromatography and the estimation of partition constants in general. Acknowledgment. We thank Pete Carr and Ray Mountain for many stimulating discussions. Financial support from the National Science Foundation (Grant CHE-9816328), a Camille and Henry Dreyfus New Faculty Award, a McKnight LandGrant Assistant Professorship, an Alfred P. Sloan Research Fellowship, and a Department of Energy Computational Science Graduate Fellowship (M.G.M.) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. References and Notes (1) Klatte, S. J.; Beck, T. L. J. Phys. Chem. 1996, 100, 5931-5934. (2) Slusher, J. T.; Mountain, R. D. J. Phys. Chem. B 1999, 103, 13541362. (3) Schure, M. R. In AdVances in Chromatography; Marcel Dekker: New York, 1998; Vol. 39, pp 139-200.
J. Phys. Chem. B, Vol. 103, No. 50, 1999 11195 (4) Tourres, D. A. J. Chromatogr. 1967, 30, 357-377. (5) Jennings, W. A. Analytical Gas Chromatography; Academic Press: Orlando, 1997. Harold, M.; McNair, J. M. Basic Gas Chromatography; Wiley: New York, 1997. (6) (a) Rohrschneider, L. J. Chromatogr. 1966, 22, 6-22. (b) McReynolds, W. O. J. Chromatogr. Sci. 1970, 8, 685-691. (7) (a) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813-826. (b) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527-545. (8) (a) Siepmann, J. I. Mol. Phys. 1990, 70, 1145-1158. (b) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59-70. (c) Mooij, G. C. A. M., Frenkel, D.; Smit, B. J. Phys: Condens. Matt. 1992, 4, L255-L259. (9) (a) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569-2577. (b) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508-4517. (10) Increasing the number of He atoms allows us to change the phase ratio so that the solutes are found with roughly equal probabilities in the vapor and liquid phases. This improves the precision of the simulations with little additional cost. (11) Lennard-Jones, J. E. Proc. R. Soc. London, Ser. A 1924, 106, 441462. (12) Kaplan, I. G. Theory of Molecular Interactions; Fraga, S., Klobukowski, M., Translation Eds.; Elsevier Science Publishers: Amsterdam, 1986 (English transl.). (13) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (14) (a) Martin, M. G.; Siepmann, J. I. J. Am. Chem. Soc. 1997, 119, 8921-8924. (b) Martin, M. G.; Siepmann, J. I. Theor. Chem. Acc. 1998, 99, 347-350. (15) Ben Naim, A. Statistical Thermodynamics for Chemists and Biochemists; Plenum Press: New York, 1992. If vapor-liquid equilibria are considered, then the standard convention is to take the liquid phase as phase R and the vapor phase as phase β. (16) Budahegyi, M. V.; Lombosi, E. R.; Lombosi, T. S.; Me´sza´ros, S. Y.; Nyiredy, Sz.; Tarja´n, G.; Tima´r, I.; Taka´cs, J. M. J. Chromatogr. 1983, 271, 213-307. (17) (a) Kovats, E. HelV. Chim. Acta 1958, 41, 1915-1932. (b) Kovats, E. In AdVances in Chromatography; Marcel Dekker: New York; 1965, Vol. 1; pp 229-247. (18) Giddings, J. C. Unified Separation Science; Wiley: New York, 1991. (19) Matire, D. E. J. Chromatogr. 1989, 471, 71-80. (20) (a) De Young, L. R.; Dill, K. A. J. Phys. Chem. 1990, 94, 801809. (b) Krukowski, A. E.; Chan, H. S.; Dill, K. A. J. Chem. Phys. 1995, 103, 10675-10688. (21) (a) Sharp, K. A.; Nicholls, A.; Fine, R. F.; Honig, B.Science 1991, 252, 106-109. (b) Kumar, S. K.; Szleifer, I.; Sharp, K. A.; Rossky, P. J.; Friedman, R.; Honig, B. J. Phys. Chem. 1995, 99, 8382-8391. (22) Hildebrand, J. H. J. Chem. Phys. 1947, 15, 225-28. (23) Ben-Naim, A.; Mazo, R. J. Phys. Chem. 1993, 97, 10829-10834 1993. (24) (a) McQuarrie, D. A. Statistical Mechanics; Harper Collins: New York, 1976. (b) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academics Press: London, 1986.
