J. Phys. Chem. 1994,98, 6675-6678
6675
Monte Carlo Simulation of the Liquid-Vapor Coexistence in a Langmuir Monolayer of Pentadecanoic Acid J. I. Siepmann,'J*@ S. Karabomi,* and M. L. Klein? Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania I91 04-6323, and KoninklijkelSheN- Laboratorium, Amsterdam, P.O. Box 3003, I003 A A Amsterdam, The Netherlands Received: April 1, 1994; In Final Form: May 30, 1994'
Configurational-bias Monte Carlo simulations in the Gibbs ensemble have been used to perform direct simulations of fluid-phase equilibria in a Langmuir monolayer. The calculations were carried out for a pentadecanoic acid monolayer using the model of Karaborni and Toxvaerd. The results show that the model monolayer exhibits a vapor-liquid coexistence with a critical point a t Tc = 104 f 5 O C and pc = 1.2 0.2 molecules/nm2, but agreement with experimental data is only qualitative. However, the results are of sufficient quality to establish the dimensionality of the model system and to judge between different experimental data sets.
Introduction The phase behavior of monolayers formed by amphiphilic molecules at the air/water interface has been the topic of intense research during the past decades.' Interest in these films has been largely motivated by their technological importance, and in addition, they can serve as model structures to gain a better understanding of one of Nature's most advanced designs, the lipid membrane. Langmuir films of pure fatty acids or phospholipids display a rich phase diagram with a variety of lowtemperature,high-density structures and three distinct fluid phases (for excellent reviews, see refs 2-5). The phases can be classified as f0llows:5-~At sufficiently low densities, monolayers exist in a gaseous state (G). A two-phase coexistence regime marks the onset of condensation to the liquid-expanded (LE) phase. At higher surface densities, many mesomorphous phases (liquidcondensed, LC) analogous to those of a smectic liquid crystal are found. Further compression results in transition to crystalline polymorphs. The experimental determination of the phase boundaries and of, if present, the critical points for the coexistence regimes is a very difficult tasks2 Experimentsto determine the 'K-Aisotherms of monolayers of pentadecanoic acid (PDA)-the system which has been studied in most detail-have been carried out as early as 1926.8 More extensive measurements were later performed by Harkins and co-workers.gJ0 In the 1970s, Kim and CanneWlJ2 showed that the LE and G phases are connected by horizontal isotherms inferring a first-order transition. They also carried out surfacepotential measurements to support the isotherm data.13 From the coexistence curve they concluded that the two-phase regime ends in a critical point at around 26 OC. In the 198% Pallas and Pethical4undertookvery careful isotherm and surface potential experiments using extremely pure systems. They confirmed the first-order nature of the L E 4 transition but found markedly different coexistingdensities and argued that the critical temperature should be higher than 50 OC.14 In addition, Pallas and Pethicals were the first to obtain horizontal isotherms for the L S L C phase regime. The controversy between the twoisotherm studies could only be resolved using a new technique,fluorescence microscopy, which allows observation of the phase coexistence directly due to different solubilities of a fluorescence impurity in each phase. 16sL7 Comparing isotherm and fluorescence data, Knobler and co-workers18 confirmed the coexistence densities t University of Pennsylvania. t
Koninklijke/Shell-Laboratorium.
I Permanent address: Department of Chemistry, University of Minnesota, 207 Pleasant St. SE, Minneapolis, MN 55455. Abstract published in Advance ACS Abstracts, July 1. 1994.
0022-3654/94/2098-6675sO4.5QIQ
reported by Pallas and Pethica and showed that the L E 4 and L E L C coexistence lines merge in a triple point. Using the fluorescence technique, phase coexistence was observed at temperatures approaching 70 OC.18 However, this value should be viewed with caution, since the effect of the fluorescenceprobe on the stabilities of the coexisting phases is unknown. In principle, computer simulations are an ideal tool to investigate the transitions of the Langmuir monolayers and the coexistence regions in microscopic detail. However, until recently it was considered impossible to calculate the phase diagram of such complex systems using conventional techniques.19 Even if reliable potential functions of the aqueous monolayer were available, direct simulations in the canonical or microcanonicalensemble cannot yield reliable data on the phase coexistence,since phase segregation is dramaticallyperturbed by the finite sizeof the simulated systems (usually N < lQQ)20-23 Recent results on the critical behavior of n-alkanes,2k27however, have indicated that the limitations of the conventional techniques can be successfully overcome using a combination of the Gibbs-ensemble approach and the configurational-bias MonteCarlomethod. In this Letter wedemonstrate that the same approach allows us to perform direct simulations for the phase equilibria in Langmuir monolayers.
Simulation Details Ih the Gibbs-ensemble Monte Carlo (GEMC) scheme,28.29 simulations of the liquid and vapor phases are carried out in parallel, utilizing two coupled simulation boxes. Monte Carlo rules which allow for changes in the number of particles and the volume ensure that the two boxes (liquid and vapor) are in thermodynamic equilibrium (same chemical potential, pressure, and temperature) with each other. The values of the excess chemical potential and the pressure are not predetermined but are an outcome of the simulation. Because the two boxes are not in 'physical" contact, there is no interface, and the bulk properties can be determined with a surprisinglysmall number of particles. This makes the GEMC method extremely efficient for phase equilibriumcalculations. However, there is one major limitation of the GEMC technique: one of the steps involves the exchange of particles between the two boxes. For liquids consisting of small molecules this does not cause a serious problem. For chain molecules, like alkanes or amphiphilic molecules, however, the probability of successfulexchanges can become vanishingly small. The configurational-bias Monte Carlo (CBMC)method was designed to avoid this problem by biasing the insertion of chain molecules.3G32 During this biased insertion, the chain is "grown' atqm-by-atom in such a way that favorable regions of space are found. CBMC insertions were first used to calculate the excess chemical potential of amphiphiles in simple monolayers.3'J The 0 1994 American Chemical Society
6676
The Journal of Physical Chemistry, Vol. 98, No. 27, 1994
Letters
r '
380.0
360.0
E CI H1
340.0
nm
Figure 1. Snapshots of the coexisting phases at 65 OC (top) and 95 OC (bottom). Vapor phases are on the left and liquid phases on the right. The bold horizontal lines represent the water/air interface (I = 0), and their lengths indicate the dimensions of the periodic cells.
bias in the GEMC particle exchange step can be removed by adjusting the Metropolis acceptance Besides allowing simulations of fluids with articulated structure in the Gibbs' ensemble, the CBMC approach has been shownto be very efficient to relax the conformational degrees of freedom which are believed to be the rate-limiting step for equilibration of monolayer^.^^^^^ The model of Karaborni and c o - w ~ r k e r s has ~ ~ ,been ~ ~ used throughout the present simulations. A PDA molecule consists of 15 pseudoatoms, each representing either the carboxylate headgroup, an internal methylene group, or the terminal methyl group. The interactions of tail segments are modeled by Toxvaerd's anisotropic united atom model which accounts implicitly for the hydrogen atoms bonded to the carbon atoms.37 Amphiphile-water interactions are represented by external potentials which define a nonzero width interface as obtained by X-ray reflectance measurement^.^^ The interactions of the headgroups are governed by a combination of excluded volume and dipolar repulsion.39 This potential model has previously been used in many molecular dynamics simulations of the high-density phases of fatty acid m o n o l a y e r ~ . ~ ~ J ~ ~ ~ Nine GEMC simulations for 60 PDA molecules were carried out covering the temperature range from 25 to 100 OC. Most of the simulations were started with the chains equally distributed between the two simulation boxes, standing normal to the surface in all-trans conformation and arranged in a herringbone pattern. However, some runs commenced from an equilibrated configuration of a run with similar temperature. Five different types of Monte Carlo moves have been used: displacement or rotation of a molecule, change of the conformation of a molecule (using CBMC), exchange of area or molecules (using CBMC) between the boxes. Further details of the simulation algorithm are given in refs 26, 34, and 41. The number of trial orientations in the CBMC moves was set to 10, resulting in acceptance rates for particle exchanges of around 0.5% for the higher temperatures and close to zero for the lower temperatures (see later). Between 30 OOO and 60 000 MC cycles were employed for equilibration of the samples, and production runs were performed over at least 30 000 cycles. The simulations were carried out on a HewlettPackard 735 workstation, and 1000 cycles require approximately 115 min of CPU time; i.e., a typical run requires around 150 h of CPU time. The coexistence densities-the primary results of the GEMC simulations-were obtained using density histog r a m ~ . ~The * errors in the densities were estimated from the histograms and from block averages. Results and Discussioas In Figure 1 snapshots of the PDA monolayer at 65 and 95 OC areshown. Theseparation intocoexisting phasesisclearly evident from the configurations. The liquid phase at 65 OC is very homogeneous, most of the chains possess conformational defects, and there is no preferred tilt direction, as expected for the LE phase. Throughout the run, the population of the vapor box fluctuates between 0 and 4 molecules with an average close
~
280.0'. ' 0.0 0.5
1.0
1.5
' 2.0
'
I
2.5
.
"
3.0
' "
3.5
4.0
p [mo~ecdes/nm '1
Figure 2. Vapor-liquid coexistenccs curves for pentadecanoic acid monolayers. The results of the isotherm measurements of Kim and Cannelll*and Pallas and Pethi~a~~arcshown as open trianglesand squares, (only given for the raspectively. The fluorescencedata of Moore et dL8 liquid branch) are represented by open diamonds. The simulation data are depicted by open circles. The solid line is the scaling law fit through the simulation data using an exponent of 0.32,the dashed line is obtained from the rectilinear law, and the filled circle is the correspondingcritical point.
to 1. The distribution of molecules in the liquid phase at 95 OC is more heterogeneous, a sign that this temperature is close to the critical point. Also, the molecules in the vapor phase begin to form loose clusters consisting of up to 10 molecules. The phase diagram of the model PDA monolayer is shown together with three sets of experimental results in Figure 2. The vapor densities for the three runs with the lowest temperatures have been omitted in the plot and in the analysis of the coexistence curve since the number of successful particle exchanges was too low to ensure well-equilibrated samples. As is evident from the figure, there is no quantitative agreement between the coexistence densities obtained from our simulations and the experimental data. The densities of the liquid phase of the model monolayer are approximately 25% higher than those obtained experimentaiiy.14.18 Using the law of rectilinear diameters and the scaling law and following the procedure described by Smit and we have tried to estimate the coexistencecurve and the critical point. The shape of the coexistence curve cannot be accurately described by fits using the Ising exponent for a two-dimensional system (j3 = 1/#445or the mean-field exponent (j3 = l/2) (see, for example, refs 46 and 47). Therefore, we have calculated numerous coexistencecurves by varying the exponent between 0.12 and 0.5 (in 0.01 intervals) to obtain the critical exponent which describes our coexistence densities. It turned out that a value close to the one expected for three-dimensional systems, i.e. j3 = 0.32,2s*46.47 yields the best agreement. This is perhaps not surprising, considering that the fluctuations perpendicular to the interface which arise from the flexibility of the amphiphilic tails together with the finite width provided by the external potentials approach themagnitudeof thefluctuationsallowed by theperiodic boundary used for the interface plane. Thus, a simulation using a much larger system size might yield a different exponent than that obtained in this work. Results of fits for the simulation and experimental densities are given in Table 1. As was already pointed out by Kim and Cannell,11J2their data can be best described using themean-field exponent. Note that our procedure of estimating the critical temperature yields an almost identical value to that obtained by Kim and Cannel1 from isothermal compressibilities.l* Estimation of the scaling exponent for the results of Pallas and Pethica14 is difficult since the data are not in the vicinity of the critical point.18 It is, however, evident that the two-dimensional Ising exponent would yield a critical
The Journal of Physical Chemistry, Vol. 98, No. 27, I994 66V
Letters
TABLE 1: Summary of the Fits of the Simulation and Experimental Data to the Law of Rectilinear Diameters and the scplina Law. ref ‘J: B ‘Jr Tc (K) pc (molcculcs/nmz) 0.125 0.32
12 12 12 14 14 14
0.5 0.125 0.32 0.5 0.125 0.32 0.5
this work this work this work
298.9 298.8 299.4 318 33 1 345 373 378 390
0.438 0.441 0.419 0.982 0.845 0.703 1.24 1.21 1.14
ur and us are the corresponding standard
1.02
I
I
0.03 1 0.03 1 0.031 0.036 0.036 0.036
0.008 0.008 0.008
0.156 0.788 0.090 0.250 0.093 0.063 0.932 0.067 0.153
deviation^.^^
I
0.90
‘
0.82 0.0
1 .o
2.4
3.0
I 4.0
Augmentation of the model to include this part of the potential should substantially lower the heat of vaporization and thereby decrease the critical temperature. In the liquid phase, this competing interaction would enhance the formation of loop and mushroom-like chain conformationsand hence reduce the density of this phase. Finally, it will also be of interest to extend the calculations and determine whether the model used throughout this work yields a L E L C coexistence.
Conclusions Five conclusions emerge from the present study. First, with a combination of Gibbs-ensemble Monte Carlo and configurational-bias Monte Carlo, it is feasible to obtain reliable estimates of the coexistence densities for a given monolayer model. Next, the monolayer model of Karabomi and Toxvaerd yields coexistence between a liquid and a vapor phase. However, the liquid phase of the model PDA monolayer is substantially denser than the liquid-expanded phase of the “real” monolayer, and thecritical point seems to be shifted to higher temperatures. Third, the coexistence curve obtained from the simulations can best be fit with a scaling exponent of approximately0.32, supporting a threedimensional character of the finite model system. Scaling of the coexistence data by the reduced temperatures and densities yields good agreement with the results of Pallas and Pethica, but not with those of Kim and Cannell. Finally, new experimental data closer to the critical point (perhaps, only feasible for shorter fatty acids or alcohols) and simulations with improved potentials and larger system sizes are highly desirable to further explore the nature of the scaling exponent for the coexistence curve.
P Figure 3. Reduced vapor-liquid cocxistences curves for pentadecanoic acid monolayers (see text). The same symbols as in Figure 2 are used.
temperature much lower than the 70 OC at which coexistence was still observed by fluorescence microscopy.’* Phase diagrams are sensitive to the choice of interaction potentials and are therefore a very good test on the overall correctness of a given potential model. Assuming that our model monolayer falls into the same universality class as the “real” PDA monolayer, we should be able to fit all data sets with the same critical exponent. Following from the law of corresponding states, a plot of the coexistencecurves in reduced coordinates (P = TIT, and p* = p / p c ) should help us to determine whether our simulation results are in agreement with one of the experimental data sets. Figure 3 shows the curves obtained from fits with /3 = 0.32. (A similar picture is obtained for /3 = 0.5.) In reduced units, the simulation densities are in very satisfactory agreement with the data of Pallas and Pethica, but those of Kim and Cannell differ significantly. This observation further supports the conclusion by Moore et a1.I8 based on their fluorescence microscopy data that the coexistencedensities reported by Pallas and Pethica are correct. Since the simulation model gives better than qualitative agreement with experiment, it should be profitable to modify our model to obtain quantitative agreement. From calculations of the phase diagram of normal we know that the anisotropic potential of Toxvaerd, which governs the tail group interactions, is not an ideal choice with the original parametrization. For example, it could be replaced by theisotropic model recently suggested by usz5-z6or with a reparametrized version of the Toxvaerd potential. Furthermore, the model of Karaborni and co-workers was developed mainly to reproduce the behavior of monolayers at high densities. For this reason the attractive forces due to dispersion interactions between tail segments and the water surface were not considered. Omission of these interactions in combination with the external potential (the energy penalty for ”solvating”methylene beadszz)manifests itself by missing the tendency of the amphiphilic molecules to lie flat on thewater surface in the vapor phase (see Figure 1).
Acknowledgment. We acknowledge many stimulating discussions with D. Frenkel, J. Hautman, K. Laasonen, I. McDonald, H. Mbhwald, C. Mundy, N. van Os, R. Pastor, and B.Smit. This work was in part supported by a National Institute of Health grant to M.L.K. (R01 GM40712). References and Notes (1) Gaines Jr., G. L. Insoluble Monolayers at Liquid Gas Interfaces; Wiley: New York, 1966. (2) Knobler, C. M. Adu. Chem. Phys. 1990, 77, 397. (3) MBhwald, H. Annu. Reu. Phys. Chem. 1990, 41,441. (4) McConnell, H.M. Annu. Rev. Phys. Chem. 1991, 42, 171. (5) Knobler, C. M.; Dcsai, R. C. Annu. Rev. Phys. Chem. 1992,43,207. (6) Harkins, W. D. The Physical Chemistryofsurface Films; Reinhold: New York, 1952. (7) Bibo, A. M.; Knobler, C. M.; Peterson, I. R. Macromol. Chem.: Macromol. Symp. 1991, 46, 55. ( 8 ) Adam, N. K.; Jessop, G. Proc. R. Soc. London A 1926,110,423. (9) Harkins, W. D.; Young, T. F.; Boyd, E. J . Chem. Phys. 1940,8,954. (IO) Harkins, W. D.; Boyd, E. J. Phys. Chem. 1941,45,20. (11) Kim, M. W.; Cannell, D. S.Phys. Rev. Lett. 1975, 35, 889. (12) Kim, M. W.; Cannell, D. S.Phys. Reu. A 1976, 13, 411. (13) Kim, M. W.; Cannell, D. S.Phys. Reu. A 1976, 14, 1299. (14) Pallas, N. R.; Pethica, B. A. J. Chem. Soc.,Faraday Trans. 1 1987,
83, 585. (15) Pallas, N. R.; Pethica, B. A. Lungmuir 1985, I, 509. (16) von Tscharner, V.;McConnell, H.M. Biophys. J. 1981, 36, 409. (17) LQche, M.; Sackmann, E.; Mbhwald, H.Ber. Bunsm-Ges. Phys. Chem. 1983,87, 848. (18) Moore, B. G.; Knobler, C. M.; Akamatsu, S.;Rondelez, F. J. Phys. Chem. 1990.94, 4588. (19) Allen, M. P.; Tildesley,D. J. Computer Simulationofliquids; Oxford University Press: Oxford, 1987. (20) Bareman, J. P.; Cardini, G.; Klein, M. L. Phys. Rev. Lett. 1988,60, 2152. (21) Bareman, J. P.; Klein, M. L. J . Phys. Chem. 1990, 94, 5202. (22) Karaborni, S.;Toxvaerd, S.;Olsen, 0. H.J. Phys. Chem. 1992,96, 4965. (23) Karabomi, S.;Siepmann, J. I. Mol. Phys., in press. (24) Siepmann, J. I.; Karaborni, S.; Smit, B. J. Am. Chem. Soc. 1993, I IS, 6454. (25) Siepmann, J. I.; Karaborni, S.;Smit, B. Nature 1993, 365, 330. (26) Smit, B.; Karaborni, S.;Siepmann, J. I. J. Chem. Phys., submitted
for publication. (27) Laso,M.;dePablo, J. J.;Suter,U. W.J. Chem.Phys. 1992.97.2817. (28) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813.
6678 The Journal of Physical Chemistry, Vol. 98, No. 27, 1994 (29) Panagiotapouloo. A. Z.; Quirke, N.;Stapleton, M.;Tildesley, D. J.
Mol. Phys. 1988.63, 527. (301 Siepmann, J. I. Mol. Phys. 1990, 70, 1145. (31) Siepmann, J. I,; Frenkel. D. Mol. Phys. 1992, 75, 59.
(32) Frenkel, D.; Mooi,G.C. A. M.;Smit, B.J. Phys.: Ondens. Matter 1992.4, 3053. (33) Mooij, G. C. A. M.;Frenkel, D.;Smit, B. 1.Phys.: Condens. Matter 1992, 4, L255. (34) Sicpmann, J. I.; McDonald,I. R. Mol. Phys. 1993, 79,457. (35) Buontmpo, J.T.;Rice,S. A.;Karaborni,S.;Siepmann,J. I. Lungmuir 1993,9, 1604. (36) Karabomi, S.;Toxvaerd, S . J. Chem. Phys. 1992,96, 5505. (37) Toxvacrd, S . J. Chem. Phys. 1990, 93,4290. (38) Braslau, A.; Pershan, P. S.;Swislow, G.;Ocko, B. M.;Als-Nielsen, J. Phys. Rev. A 1988,38,2457.
Letters (39) vander Ploeg, P.; Berenhn, € J.IC. . J. Chem. Phys. 1982,76,3271. (40) Karabomi, S.Lungmuir 1993, 9, 1334. (41) Siepmann, J. I. In Computer Simulation ofBiomolecular Systems: Theoretical and Experimental Applfcations;van Gunsteren, W. F., Weiner, P. K., Wilkinson. A. J., Eds.; Escom Sciencc Publisher: Leiden, 1993; p 249. (42) Smit, B.; de Smedt, P.; Frenkel, D. Mol. Phys. 1989,68,931. (43) Smit. B.; Williams, C. P. J. Phys.: Condens. Matter 1990,2,4281. (44)Onsager, L. Nuovo Cimenro (Suppl.) 1949,6261. (45) Kim,H. K.;Chan, M.H. W. Phys. Rev.k i t . 1%4,53, 170. (46) Stanley, H. E. Introduction to Phase Transitions and Critical Phenomena; Oxford University Press: New York, 1971. (47) Binney, J. J.; Dowrick, N. J.; Fisher, A. J.; Newman, M.E. J. The Theory of Critical Phenomena; Clarendon Press: Oxford, 1992.
