1334
Langmuir 1993,9, 1334-1343
Molecular Dynamics Simulations of Long-chain Amphiphilic Molecules in Langmuir Monolayers S.Karaborni KoninklijkelShell-Laboratorium,Amsterdam (Shell Research BV), P.O. Box 3003, 1003 AA Amsterdam, The Netherlands
Received October 13,1992. In Final Form: January 25,1993 Molecular dynamics simulations have been used to investigate phase transitions in a monolayer model of long-chain amphiphilic molecules (20 pseudoatoms) at the air-water interface. Several simulations were performed in the density range near close packing in which a number of phases have been identified experimentally (18.6-25 A2/molecule). Between 24 and 25 A2/moleculethe monolayer undergoes a firstorder phase transition exhibited by a break-up of lattice structure,a strong decrease in trans fraction and end-to-end distance, and an onset of center of mass translational displacement. Between 18.5 and 21 A2 the relaxation time for molecules rotating around their major axes decreases significantly, the average chain conformation in the molecule midsections changes from a mostly trans to a gauche defected conformation, and the frequency of the head group movements normal to the monolayer plane changes noticeably. Simulationsof monolayers of short-chain amphiphilic molecules (16pseudoatoms) indicate that the first-order phase transition does not have the same characteristics as that in the long-chain monolayer.
Introduction Sincethe early surfacetension measurementson oil films by Langmuir' and Harkins: monolayers of amphiphilic molecules at the air-water interface, subsequently known as Langmuir monolayers, have attracted continuous interest from the research ~ommunity.~ In the 1940s the remarkable experiments of Harkins et al.? Dervichian? and Stenhagen and Stenhagen6have strongly questioned the simplicity of Langmuir monolayers. Today, there is mounting evidence that monolayers of amphiphilic molecules at the air-water interface display a very complex phase behavior? but the characteristics of the different phases and the phase transitions, such as the LC-LE or LZ-Ll transitions, are still a subject of debate.at9 The onset of new, powerful experimental techniques such as synchrotron X-ray diffraction has undoubtedly resulted in a molecular understanding of Langmuir monolayers,especiallytheir lattice geometryand molecular tilt and the use of fluorescence microscopy has (1)Langmuir, I. J. Am. Chem. SOC. 1917,39, 1848. (2)Harkins, W. D. J. Am. Chem. SOC. 1917,39,354 and 451. (3)Knobler, C. M. In Advances in Chemical Physics; Rice, S. A., Ed.; Wiley: New York, 1988. (4)Harkins, W. D.; Copeland, L. E. J. Chem. Phys. 1942,10, 262. Copeland, L. E.; Harkins, W. D.; Boyd, G. E. J. Chem. Phys. 1942,10, 357. (5)Dervichian, D. G. J. Chem. Phys. 1939,7 , 931. (6)Stiillberg-Stenhagen, S.;Stenhagen, E. Nature 1945,156,239. (7)Bibo, A. M.; Knobler, C. M.; Peterson, I. R. J. Phys. Chem. 1991, 95,5591;Makromol. Chem. Macromol. Symp. 1991,46,55.Bibo, A. M.; Peterson, I. R. Adu. Mater. 1990,2,309.Peterson, I. R. Ber. Bunsen-Ces. Phys. Chem. 1991,95,1417.Knobler, C . M.J.Phys.: Condens. Matter 1991,3,S17. (8)Hifeda, Y. F.;Rayfield, G.W. Langmuir 1992,8, 197. (9)Kato, T.; Akiyama, H.; Tanaka, T. Chem. Phys. Lett. 1991, 184, 455. Kato, T.; Hirobe, Y.; Kato, M. Langmuir 1991,7, 2208. (10)Dutta, P.; Peng, J. B.; Lin, B.; Ketterson, J. B.; Prakash, M.; Georgopoulos, P.; Ehrlich, S. Phys. Rev. Lett. 1987,58,2228. (11)Barton, S. W.; Thomas, B. N.; Flom, E. B.; Rice, S. A.; Lin, B.; Peng, J. B.; Ketterson, J. B.; Dutta, P. J. Chem. Phys. 1988,89,2257. (12)Kjaer, K.; Ala-Nielsen, J.; Helm, C. A.; Tippmann-Krayer, P.; Mohwald, H. J. Phys. Chem. 1989,93,3200. (13)Jacquemain, D.; Grayer-Wolf, S.; Leveiller, F.; Lahav, M.; Leiserowitz, L.; Deutach, M.; Kjaer, K.; Ab-Nielsen, J. Colloq. Phys. C7 1989,50,C7-29. (14)Lin, B.;Peng, J. B.; Ketterson, J. B.; Dutta, P.; Thomas, B. N.; Buontempo, J.; Rice, S. A. J. Chem. Phys. 1989,90,2393. (15)Lin, B.; Shih, M. C.; Bohanon, T. M.; Ice, G. E.; Dutta, P. Phys. Rev. Lett. 1990,65,191.
led to the direct observationof domain structure and phase Similarly,new microscopytechniquesthat do not use fluorescentprobes are being developed to study domain morphology and the order of phase transitions in the absence of probe i m p u r i t i e ~ . ~Unfortunately, ~-~~ the molecular conformation cannot yet be described quantitativelyusing the abovetechniques,but it has to be inferred from other measurable quantities such as the electron density in the nonpolar section of the monolayer.12 In fact, in all X-ray experiments the chain conformation is believed to play a minor role, thus the moleculesare always assumed to be trans in the interpretation of the data.I2 These assumptions are in contrast with nonresonance Raman studies of Langmuir monolayers at high density which indicate the presence of gauche defects that grow upon monolayer e x p a n ~ i o n .In ~ ~support of the Raman studies, recent X-ray diffraction experiments on fluorinated alkyl acid monolayers show a completely different behavior from fully hydrogenated acid monolayers as a result of the difference in chain flexibility between fluorinated and hydrogenated alkyl acids.31132 (16)Bohanon, T. M.; Lin, B.; Shih, M. C.; Ice, G. E.; Dutta, P. Phys. Rev. B 1990,41,4846. (17)Kenn,R.M.;Bijhm,C.;Bibo,A.M.;Peteraon,I.R.;M6hwald,H.; Ale-Nielsen, J.; Kjaer, K. J. Phys. Chem. 1991,95,2092. (18)Jacquemain, D.; Leveiller, F.; Weinbach, S. P.; Lahav, M.; Leiserowitz, L.; Kjaer, K.; Als-Nielsen, J. J . Am. Chem. SOC. 1991,113, 7684. (19)Schlossman, M. L.; Schwartz, D. K.; Pershan, P. S.; Kawamoto, E. H.; Kellog, G. J.; Lee, S. Phys. Rev. Lett. 1991,66,1599. (20)Tippmann-Krayer, P.; MBhwald, H. Langmuir 1991, 7,2303. (21)Shih, M. C.;Bohanon, T. M.; Mikrut, J. M.; Zschack, P.; Dutta, P. J. Chem. Phys. 1992,96,1556. (22)MBhwald,H. ThinSolid F i l m 1988,159,l.L&che,H.;MBhwald, H. Colloids Surf. 1984,10,217.MBhwald, H.; Kenn, R. M.; Degenhardt, D.; Kjaer, K.; Ale-Nielsen, J. Physica A 1990,168,127. (23)Moore, B. G.;Knobler, C. M.; Akamtau, S.; Rondelez, F. J. Phys. Chem. 1990,94,4588.Moore,B. G.;Knobler,C. M.;Broeeta,D.;Rondelez, F. Faraday Symp. Chem. SOC. 1986,20, 1753. (24)Qiu, X.;Ruiz-Garcia, J.; Stine, K. J.; Knobler, C. M. Phys. Rev. Lett. 1991,67,703. (25)Muller, P.; Gallet, F. Thin Solid F i l m 1992,210/211,138. (26)Ivanov, G.R. Chem. Phys. Lett. 1992, 193, 323. (27)Henon, S.; Meunier, J. Thin Solid F i l m 1992,210/211,121. (28)Takoshima, T.;Masuda, A.; Mukasa, K. Thin Solid Films 1992, 210/211,51. (29)HBnig, D.; NBbius, D. Thin Solid F i l m 1992,210/211,64. HBnig, D.; Overbeck, G.A.; MBbius, D. Adu. Matter 1992,4,419. (30)Kawai, T.; Umemura, J.; Takenaka, T. Chem. Phys. Lett. 1989, 162,243.
0743-7463/93/2409-1334$04.00/00 1993 American Chemical Society
Simulations of Amphiphilic Molecules
Langmuir, Vol. 9,No. 5,1993 1335
In addition to chain conformation and monolayer tilt which the monolayers or bilayers are modeled usingexplicit order, molecular rotation is believed to play an important atoms or groups of atoms (UA),and the solvents are role in the phase behavior of Langmuir m ~ n o l a y e r s . ~ ~ J ~ Jmodeled ~ using external potentials such that the head Monolayers of heneicosanol have been shown to undergo, groups are either fixed in the z-direction or constrained to lie near a plane.56-67 between 5 and 10 OC, a phase transition from a nonrotator and distorted hexagonal phase to a rotator and hexagonal In this study, molecular dynamics simulations are used phase.11J4 Similarly, Kenn et al. have speculated after to investigate the effect of chain conformation,molecular measuring the lattice spacings in different phases that rotation, and lattice order/disorder on the phase behavior molecular rotation is possible in one of the phases that of Langmuir monolayers of molecules having an internal exists at high pressure and temperature, and that restricted flexibility similar to that of alkyl amphiphiles. rotation is possible in other phases that occur at lower pressure or temperature.17 Hautman and Klein also note Model that, in simulations of self-assembled monolayers of The monolayer model used here has been described in alkanethiol on gold, the model monolayer undergoes, with earlier publication^.^^^ The special features of the present rising temperature, a change from a nonrotator to a free model include a provision for the movements, as in real rotator phase via a locked rotator stage.33 monolayers, of head groups and chain methylenes in and As a consequence of the importance of amphiphilic out of the water and the formation of a finite size interface. monolayers from a biological, industrial, and a theoretical In addition,the chain-chain interactionsare modeled using point of view,M there have been over the last decade several an anisotropic united atom model that closely resembles attempts to model monolayers. These models can be methylene-methylene interactions in the full atomic divided into two categories: those that use a rod repremodel. sentation of molecule^^^ and those that stress internal The model amphiphile used in these simulations is molecular conformation^.^^^ Both approacheshave been similar to the one described in ref 70; a skeletal chain partially successful in predicting the complexity of the composed of 19equal diameter pseudoatomsfor the methyl phase behavior of Langmuir monolayers. and methylene segments and a carboxylate-like pseudoIn parallel with the development of new experimental atom for the head The methyl and methylene techniques and mean field and Ising models, molecular pseudoatoms and head group on the same chain interact dynamics and Monte Carlo simulations have been used to via angle bending and torsion forces as well as a Lennardstudy monolayers and bilayers. These can be divided into Jones potential between pseudoatomsthat are separated two classes: those that use full molecular representations by at least three spheres,but the bonds are kept constant.72 to model the monolayers and s ~ l v e n t and a ~ ~those ~ in The angle bending potential is given by (31)Barton, S. W.; Goudot, A.; Boulaesa, 0.; Rondelez, F.; Lin, B.; Novak, F.; Acero, A.; Rice, S. A. J. Chem. Phys. 1992,96,1343. (32)Shin, S.;Collazo, N.; Rice, S. A. J. Chem. Phys. 1992,96,1352. Collazo, N.; Shin, S.; Rice, S. A. J. Chem. Phys. 1992,96,4735. (33)Hautman, J.;Klein, M. L. J. Chem. Phys. 1989,91,4994;J.Chem. Phys. 1990,93,7483.Hautman, J.; Bareman, J. P.; Klein, M. L. J.Chem. SOC.,Faraday Trans. 1991,87,2031. (34)Swalen, J. D.; Allara, D. L.; Andrade, J. D.; Chandrose, E. A.; Garoff, S.; Israelachvili, J.; McCarthy, T. J.; Murray, R.; Pease, R. F.; Rabolt, J. F.; Wynne, K. J.; Yu, H. Langmuir 1987, 3,932. (35)Pearce, P. A.; Scott, H. L. J. Chem. Phys. 1982,77,951. (36)Safran, S. A.; Robbins, M. 0.;Garoff, S. Phys. Rev. A 1986,33, 2186. (37)Hdperin, A.; Alexander, S.; Schechter, 1.J. Chem. Phys. 1989,91, 1383;1987,-86,6550. (38)Chen. Z.-Y.: Talbot,, J.:. Gelbart. W. M., Ben-Shaul, A. Phvs. Rev. Lett. iw, 61,1376. (39)Kreer, M.; Kremer, K.; Binder, K. J. Chem. Phys. 1990,92,6195. (40)Costas, M. E.; Wang, Z.-G.; Gelbert, W. M. J. Chem. Phys. 1992, 96,2228. (41)Somoza, A. M.; Desai, R. C. J. Phys. Chem. 1992,96,1401. (42)Kramer, D.;Ben-Shad, A.; Chen, Z.-Y.; Gelbart, W. M. J. Chem. Phys. 1992,96,2236. (43)Scheringer, M.; Hifler, R.; Binder, K. J. Chem. Phys. 1992,96, 2269. (44)Kaganer, V. M.; Osipov, M. A.; Peterson, I. R. J. Chem. Phys., in press. (45)Marcelja, S. Biochim. Biophys. Acta 1974,367, 165. (46)Gruen, D. W. R. Biochim. Biophys. Acta 1980,595,161. (47)Lemaire, B.; Bathorel, P. Macromolecules 1980,13,311. (48)Mouritsen, 0.G. InMoleculur DescriptionofBiologicalMembrane Components by Computer Aided Conformational Analysk Brasseur, R., Ed.; CRC Press, Inc.: Boca Raton, FL, 1990,Vol. I. Mouritsen, 0. G. Chem. Phys. Lipids 1991,57,179. (49)Cantor,R.S.;McIlroy,P. M. J.Chem.Phys. 1989,90,4431.Cantor, R. S.Mater. Res. SOC.Symp. h o c . 1990,177,345. (50)Szleifer, I.; Ben-Shad, A.; Gelbart, W. M. J. Phys. Chem. 1990, 94,5081. (51)Yamamoto, T. J. Chem. Phys., 1990,93,5990. (52)Egberta, E.; Berendsen, H. J. C. J. Chem. Phys. 1988,89,3718. (53)Smit, B.; Hilbers, P. A. J.; Eseelink, K.; Rupert, L. A.M.; van Os, N. M.; Schlijper, A. G. Nature 1990,348,624;J. Phys. Chem. 1991,95, 6361. Hilbers, P.A. J.; Rupert, L. A. M.; Eseelink, K.; van Os,N. M.; Smit, B. Submitted for publication in Langmuir. (54)Bbcker, J.;Schlenkrich, M.; Nicklas, K.; Brickmann, J.; Bopp, P. J. Chim. Phys. 1991,88,2535. (55)Raghavan, K.; Rami Teddy, M.; Berkowitz, L. Langmuir 1992,8, 233.
60 = 112.15O is the equilibrium bond angle, Bi is the angle between pseudoatoms i, i + 1,and i + 2, and Yb = 1.3 X lo5 J/mol is the bending vibration force constant. The torsion potential used here is similar to the potential used by Ryckaert and belle man^^^
u(ei)= yr (1.116 - 1.462 COS ei - 1.578 cos2+i+ 0.368 cos3 di + 3.156 cos4 &i + 3.788 cos5&) (2) +yr = 8313 J/mol is the constant given by Ryckaert and
and & is the dihedral angle between segments
i - 1,i, i + 1, and i + 2.
(56) Kox, A. J.; Michels, J. P. J.; Wiegel, F. W. Nature (London) 1980, B7, 317. (57)Vacatello, M.; Busico, V.; Corradini, P. Gazz. Chim. Ztal. 1984, 114,117. Vacatello,M.;Busico,V.Mol.Cryst.Liq. Cryst. 1984,107,341. (58)Northrup, S. H.; Curvin, M. S. J. Phys. Chem. 1985, 89, 4707. (59)Khalatur, P. G.; Pavlov, A. S. Makromol. Chem. 1987,188,3029. (60)Cardini, G.; Bareman, J. P.; Klein, M. L. Chem. Phys. Lett. 1988, 145,493. Bareman, J.P.; Cardini, G.; Klein, M. L. Phys. Rev. Lett. 1988, 60,2152. Bareman, J. P.;Klein, M. L. J.Phys. Chem. 1990,94,5202. (61)Harris, J.; Rice, S. A. J. Chem. Phys. 1988, 89,5898. (62)Harris, J.; Rice, S. A. J. Chem. Phys. 1988,88,1298. (63)Millik, M.; Kolinski, A.; Skolnick, J. J. Chem. Phys. 1990,93, 4440. (64)Bishop, M.; Clarke, J. H. R. J. Chem. Phys. 1991,95,540. (65)Moller, M. A.; Tildesley, D. J.; Kim, K. S.; Quirke, N. J. Chem. Phys. 1991,94, 8390. (66) Makovsky, N. N. Mol. Phys. 1991, 72, 235. (67)Siepmann, J. I.; McDonald, I. R. Mol. Phys. 1992,75,255. (68)Karabomi, S.;Toxvaerd, S. J. Chem. Phys. 1992,96,5505. (69)Karaborni, S.;Toxvaerd, S.; Olsen, 0. J.Phys. Chem. 1992,96, 4965. (70)Karaborni, S.;Toxvaerd, S. J . Chem. Phys. 1992, 97, 5876. (71)van der Ploeg, P.; Berendsen, H. J. C. J. Chem. Phys. 1982, 76, 3271;Mol. Phys. 1983,49,233. (72)Edberg, R.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1986,84, 6933. (73)Ryckaert, J. P.;Bellemans, A. Chem. Phys. Lett. 1976,30,123. Discuss. Faraday SOC. 1978,60,95.
Karaborni
1336 Langmuir, Vol. 9, No. 5, 1993
The methyl and methylene spheres interact using a Lennard-Jonespotential and using an anisotropic united atom model (AUA)74
uu = 4 4
(5)"- (;)"I
(3)
where u = 3.527 A and Q = 665 J/mol. The head-head interactions in these simulations are purely repulsive and consist of a strong dipolar repulsion and an excluded volume term (4)
where UH = 4.22 A and QH = 920 J/mol as previously used by van der Ploeg and B e r e n d ~ e n .Head-methyl ~~ and head-methylene spheres also interact via a Lennard-Jones potential
where UHM = 0.5(u + UH) and QHM = (cQH)'/'. The methylene-water interaction is modeled using an external potential based on the free energy of transfer of methylenes from a hydrophobic solvent to water75J6
z 1alu
where 61 = 8.27, a1 = 6, and TI = -32. The head group-water interaction is modeled using a similar potential to the methylene-water interaction. The head group-water potential is based on the free energy of transfer of a carboxylate from water to a hydrophobic solvent76
Uheadgroupwater(Z)
{
= O 1
z Ia26 BZC
+ (z/(azu) -lITZ
z
> azu (7)
where 8 2 = 82.7, 7 2 = -16, and a2 = -6. The external potentials used here are different from the more conventional external potentials that either fix the head group in the z-direction or use an infinite barrier type potential ((1lt-p - (l/r)").Extensive molecular dynamics simulations of the transfer of ions from a nonpolar to a polar solvent support the use of finite energy barriers near a liquid interface.77
Simulation Details Eight molecular dynamics simulations have been performed on monolayers in the density range of 18.5-25 A2/ molecule. The main simulation box consisted of 64 molecules periodically replicated in the x and y directions. The system was initially started by placing the molecules in hexagonal unit cells. For this purpose, the usual square box was not used; instead the box length in the x-direction was (3/4)1/2the length in the y direction. The leapfrog algorithm78was used to solve the equations of motion at (74) Toxvaerd, S. J. Chem.Phys. 1990,93,4290. Padilla, P.;Toxvaerd,
S . J. Chem. Phys. 1991,94,5650; J. Chem. Phys. 1991,95,509. (75) Karaborni, S.; OConnell, J. P. J. Phys. Chem. 1990, 94, 2624;
Langmuir 1990,6, 911. (76) Vilallonga,F. A.; Koftan,R. J.; OConnell, J. P. J.Colloidlnterface Sci. 1982, 90,539. (77) Benjamin, I. J. Chem. Phys. 1992,96, 577. (78) Toxvaerd, S. Mol. Phys. 1991, 72, 159.
1 -10 - 1 18
22
26 30 Area/molecule A*
34
38
Figure 1. Lateral pressure T versus aredmolecule at 300 K (T
units: c/o*, e = 665 J/mol, u = 3.527 A). (a) Monolayer of 20 pseudoatom molecules. (b) Monolayer of 16 pseudoatom molecules using a square simulation box. The isotherm is shifted 5 Azfor clarity. (c) Monolayer of 16 pseudoatommolecules using a simulation box in which the x-length of the box is equal to (3/4)1/2the y-length. The isotherm is shifted 10 Azfor clarity.
a time interval of the order of 1fs, and the temperature was kept constant at 300 K using a Nose-Hoover thermostat. To prevent the initial conditions from affecting the final results, the head groups were first not allowed to move, but they were constrained to lie near a wall using a harmonic potential. Similarly, the chain methylenes were not allowed to cross the z = 0 plane by applying an excluded volume term ( q - 9 After significant equilibration of the constrained monolayer, the constraints were removed and the finite energy potentials were applied (eqs 6 and 7). Then, a second equilibration period followed during which the conditions of trans fraction, interfacial tension, and molecular temperature stability were reached. Then, the actual run at 18.5 A2/molecule was executed. The simulations at higher molecular areas were performed by slowly expanding the monolayer from 18.6 to 25 A2/ molecule using increments no higher than 1 A2. Each simulationat every density was started by an equilibration period of approximately2000&4oooO time steps and was followed by an analysis period of at least 100 OOO time steps. The particle positions are stored on disk every 100 time steps and then analyzed.
Results In addition to reporting the results from the monolayer discussed in the previous section, additional simulations of two monolayers constructed from molecules containing 15 pseudoatoms and a head group will be di~cussed.@~6~ These will be referred to as the short chain monolayers. The present simulation of molecules containing 19 pseudoatoms and a head group will be referred to as the long chain monolayer. The lateral pressure-area isotherms calculated using the virial route79@are shown in Figure 1 for three monolayers. The isotherm from the long chain monolayer is compared to the lateral pressures in the short chain monolayers calculated using a square simulation box and a rectangular box. The isotherms from the two short-chain monolayers show first-order phase transition^.^^ Similarly, the long chain monolayer undergoes a first-order phase transition as indicated by the change in slope sign of the isotherm (79) Buff, F. P. 2.Elektrochem. 1962,56, 311. (80)Toxvaerd, S. J. Chem. Phys. 1981, 74, 1998.
Simulations of Amphiphilic Molecules
at 24 A2/molecule. The transitions in the short chain monolayers occur at a lower area per molecule and at a higher lateral pressure than the transition in the long chain monolayer. The transitions in the short chain monolayers occur at =22 A2/moleculeand a lateral pressure significantly higher than zero while the transition in the long chain monolayer occurs between 24 and 25 A2/molecule and at a pressure of -0. Despite the absence in the simulations of water and counterion molecules and the roughness of the head group potentials, the isotherms calculated have somewhat the same shapes as the isotherms measured in monolayers of C14and higher.6~7~12~18~23~30161-91 However, it is unexpected from the limited number of simulation points to reproduce experimental details such as kinks or plateaus with very small width. Structure: Structure Factor. The following is an attempt to characterize the phase transition in the long chain monolayer between 24 and 25 A2/moleculeusing a mixture of order/ disorder and conformation analysis. A good measure of lattice order/disorder is the structure factor
Langmuir, Vol. 9, No. 5, 1993 1337
200,
. NT
NTis the total number of segments. Rj = (R,,Ry) is the projection of segment j positions in the (x,y) plane, and k = (k,,ky) is a two-dimensional scattering vector whose allowed coordinates are multiples of 2dL, and 27r/Ly,L, and Lybeing the lengths of the simulation box in the x and y directions. The structure factors at 24 and 25 A2/moleculeare shown in Figure 2. The structure factor at 24 A2/molecule is indicative of a crystal structure in which molecules are tilted over their nearest neighbors. Upon expansion of the monolayer to 25 A2/molecule,the peak heights decrease drastically, the peak positions are no longer well-defined, and a ring forms between the various peaks, implying a disappearance of lattice order. Structure: Voronoi Polygons. A second measure of orderldisorder in two dimensional crystals is the calculation of Voronoi polygons. Voronoi polygons are a generalization of the Wigner-Seitz cell to a system not on a lattice.a2 The Voronoi polygons in this case measure the number of disclinations by calculating the fraction of molecules with coordination numbers different from 6 (usually 5 or 7).92193 Between 18.5 and 23 A2/molecule, only occasional disclimination quartets of the 7557 kind are noticed, and the number of molecules that have coordination numbers different from 6 is insignificant (81) Iwahashi, M.; Maehara, N.; Kaneko, Y.; Seimiya, T.; Middleton, 5.; Pallas, N. R.; Pethica, B. A. J. Chem. SOC.,Faraday Trans. 1 1985, 81, 973.
(82) Middleton,S. R.; Iwahashi, M.; Pallas, N. R.; Pethica, B. A. Proc. R. SOC.London, A 1984,396,143. (83) Pallas, N. R.; Pethica, B. A. Langmuir 1985,1, 509.
(84) Fiecher, A.; Sackmann, E. J. Colloid Interface Sci. 1986,112, 1. (85) Benattar, J. J.; Daillant, J.; Belorgey, 0.;Boaio, L. Physica A 1991, 172, 225. (86)Lundquist, M. Chem. Scr. 1971,1,5; Chem. Scr. 1971,1, 197. (87) Benettar, J. J.; Daillant, J.; Bosio, L.; Leger, L. Colloq. C7 . Phys. . 1989,50, C7-39. (88)Gainee. G. L. Insoluble Monolayers at Liquid-Gas Interfaces; Interscience: New York, 1966, p 220. (89) Uyeda, N.;Takenaka, T.;Aoyama, K.; Matsumoto, M.;Fujiyoshi, Y. Nature 1987,327,319. (90)Rettig, W.; Dbrfler, H.-D. Mater. Sci. Forum 1988, 25-26, 577. (91) Kato, T.; Iriyima, K.; Araki, T. Thin Solid Films 1992,210/211, 79. (92) Strandburg, K. J. Reu. Mod. Phys. 1988, 60, 161. (93) McTague, J. P.; Frenkel, D.; Allen, M. P. In Ordering in Two Dimensions; Sinha, S . K., Ed.; North Holland: Amaterdam, 1980.
Figure 2. Structure factor S ( k ) at (a, top) 24 A2/moleculeand (b, bottom) 25 A2/molecule.
= 5 I 0.5
-u.I 18
19
20
21
22
23
24
25
26
Aredmolecule (A’) Figure 3. Voronoi polygon construction: average number of disclinations/molecule in the center of mass plane as a function of area/molecular. The number of disclinations is defined as the number of molecules that have coordination numbers different from six. (Figure 3). At 24 A2/molecule the number of defects becomes higher than zero, and at 25 A2/molecule the
number of defects grows rapidly signaling the loss of triangular order in the monolayer. Although the number of defects at 24 A2/moleculeis higher than zero, the lattice structure is preserved as indicated by the two distinctive peaks in the structure factor. Most of the defecta that appear at 24 A2/moleculeare disclination quartets of the 7557 kind; some of the disclinations are, however, no longer grouped into quartets.92 At 25 A2/molecule,a different defect behavior is present in which some molecules have 8 and 4 nearest neighbors. Average Chain Conformation. Besides the lattice structure, the chain conformation usually plays an im-
Kara borni
1338 Langmuir, Val. 9, No. 5,1993 25
6.9 0.99 6.8 .-6
0.97 6.7
2
0.95
B
0.93
I-
8l
1
6.6
1 2 I
0.91
6.5
20
$
15
cr2 >
&
u 6.4
0
trans fraction
0 8879
5
6.3 6.2
18
19
20
21 22 23 aredmolecule ' A
24
25
26
Figure 4. Trans fraction, and average end-to-end distance as a function of area/molecule. Average end-to-end distance is in units of u.
portant role in the phase behavior of molecules having internal Conformation. To investigate the chain conformation of the amphiphilic molecules, the trans fraction, the end-to-end distance, and the radii of gyration have been calculated. The trans fraction and the end-to-end distance are plotted in Figure 4 as a function of area per molecule. The trans fraction curve is mirrored by the end-to-end distance curve; therefore, future discussions of chain conformation will only focus on the trans fraction. The trans fraction is the number of dihedral angles in trans conformation over the total number of dihedral angles. At 18.5 k/molecule, 99.5% of the dihedral angles are in trans conformation.60~~~94 As the monolayer is expanded, the number of gauche defects increases, in agreement with nonresonance Raman experiments of stearic acid monolayer^.^^ Between 24 and 25 &Vmolecule the fraction of angles in trans conformation decreases from 0.95 to 0.90, a decrease that is more than 5 times the average decrease between 18.5 and 24 AZ/molecule (0.008/A2). Phenomenological Ising models of biological membranes (Pink modelsQ5)also predict a decrease in trans fraction at the LC-LE and gel-liquid phase transitions.48 The molecular conformation behavior of the monolayer prior to the first-order transition is very similar to the well-documented defect formation in long chain rt-alkanes prior to melting.96-102 Translational Self-Diffusion. The translational diffusion in the monolayer plane is calculated for mass centers from mean square displacements using the relation
e--
4t
where
and N is the total number of molecules. xi(t), xi(O), yi(t), (94) Collazo, N.; Rice, S. A. Langmuir 1991, 7 , 3144. (95) Georgallas,A.; Pink, D. A. J. Colloid Interface Sci. 1982,89,107. (96) Ewen, B.; Fischer, E. W.; Piesczek, W.; Strobl, G. J. Chem. Phys. 1974,61,5265. (97) Zerbi, G.; Magni, R.; Gussoni, M.; Moritz, K. M.; Bigotto, A,; Dirlakov, 5.J. Chem. Phys. 1981, 75, 3175. (98) Northrup, S. H. J. Phys. Chem. 1984,88, 3441. (99) Maroncelli, M.; Qi, S. P.; Strauss, H. L.; Snyder, R. G. J. Am. Chem. SOC.1982,104,6237. (100) Craievich,A.; Denicolo, I.; Doucet, J.Phya. Rev. E 1984,30,4782. (101) Maroncelli, M.; Strauss, H. L.; Snyder, R. G. J. Chem. Phys. 1985,82, 2811. (102) Kim, Y.;Strauss, H. L.; Snyder, R. G. J. Phys. Chem. 1989,93, 7520.
0 0
5.0 10'"
1.0 I O L o
1.5 l o L o
2.0 10"
Time (ps)
Figure 5. Center of mass displacement in the monolayer plane at 25 A2/molecule. and yi(0) are the centers of mass coordinates of molecule i at time t and at time 0. Between 18.5 and 24 A2/molecule D is zero. As the monolayer is expanded to 25 A2/moleculethe centers of mass displacements become greater than zero. The average translationaldiffusion calculated at 25 &/molecule (Figure 5) is equal to 2.5 X 1o-Scm2/s. This value should be considered only as approximate since a correct calculation of the translational diffusion requires significantly longer times than the simulation length. Consideringthat no water molecules were explicitly included in the simulations, and the limited duration of most sampling times to 100 ps, the diffusion constant calculated at 25 A2/ molecule is still of the same order of magnitude as the experimental value in a decanol-decanoatebilayer (1.3 X 10-8 cm2/s)lo3and that in a stearic acid monolayer (1.6 X 10-5 cm2/s).lW The molecules at 25 A2 do not appear to be in a true diffusive state. The simulation is divided into stages in which the molecules are displacing and other stages in which molecules do not appear to be diffusing. Voronoi polygons, structure factor plots, chain conformation, and translational diffusion analysis indicate that the phase transition between 24 and 25 A2/molecule is both a lattice break-up as well as a chain melting. On the basis of similar analysis in short chain monolayers, it was previously concluded that the monolayer melting is only a break-up of the lattice s t r u c t ~ r esince , ~ ~ hardly ~ ~ any change in the trans fraction was seen. It appears now that the main transition could be ascribed to either or both chain melting and lattice break-up depending on the molecule chain length. Figure 6 showstwo x-z plane snapshots of the monolayer at 24 and 25 A2/molecule,indicating the phase transition. The monolayer goes from a state where molecules are tilted over their nearest neighbors to a state where molecules are no longer ordered. Other x-y plane snapshots of the monolayer at 25 &/molecule (not shown) show mostly coexisting regions of disordered areas with ordered areas that appear at separate times to be either tilted or nontilted. Molecular Reorientation. Hautman and Klein have discovered in a simulation that self-assembledalkanethiol molecules on gold undergo several transitions in molecular reorientation upon increases in temperature.33 We at(103) Schindler, H., Seelig, J., J. Chem. Phys. 1973,59,1841; 1974,62, 2946. (104) Vassilieff, C. S., Manev, E. D. Ivanov, I. B. Abhaud. Akad. Wiss. DDR Abtl. Math. Naturwiss. Technik. 1986, I , 465.
Langmuir, Vol. 9, No. 5, 1993 1339
Simulations of Amphiphilic Molecules
7,J I'
D 1
I
i
1 I
Figure 6. Snapshots of the x-z plane of the monolayer at (a) 24 A2/moleculeand (b) 25 A2/molecule. Head groups are displayed in yellow and CH2 pseudoatoms are displayed in blue.
tempt here to discern if similar transitions exist in a Langmuir monolayer upon lowering the density. To characterize the molecular reorientation, the approach
taken here is similar to that of Hautman and Klein.33It is briefly described as follows: Hautman and Klein proposed a vector that lies in the C-C-C planes and bisects
Kuru borni
1340 Langmuir, Vol. 9, No.5, 1993 1
0.8
0.6
Ytt) 0.4
t
0.2
0
2.0 Time (sec)
1.0 io-11
io-"
3.0 10-11
Figure 7. Time autocorrelation functions of the reorientation of molecules around their major axes at 18.5-22 &/molecule.
the C-C-C angle. In other words, a vector that lies in the molecular plane, perpendicular to the long axis of the molecule
t 5
18
19
20
21
22
23
24
25
26
Aredmolecule where ria is the vector from pseudoatom i + (i + 1) in molecule CY.dcc is the C-C bond length, BCc is the C-C-C bending angle, and the sum is only performed for those pseudoatoms that are distant from the head group and the end methyl segment (nl = 3 and n2 = 17). The reason is that the usual presence of gauche defects at the chain extremities tends to abnormally affect the value of R,(t). The reorientation of the molecules around their major axes is calculated using a time autocorrelation function y(t)
The relaxation times T~ are then determined by fitting y ( t ) to an exponential function y(t) = A exp(-t/.r,)
(13)
Figure 7 shows the relaxation of y ( t )between 18.5 and 22 AZ/molecule. The relaxation a t 18.5 A2/moleculeis clearly slower than the relaxation at 20,21, and 22 A2/molecule. At 19 A2/molecule y ( t ) decays with an intermediate relaxation rate. In Figure 8, T-, is plotted as a function of molecular area. At 18.5 A2/molecule T~ = 17.5 ps, then it decreases as the molecular area is increased reaching a minimum of =57 ps between 20 and 22 A2/molecule. This value is similar to the relaxation time in the rotator phase of n-alkaneslo5JNand the relaxation time in the rotator phase of a model n-alkanethiol m0nolayer.3~The decrease in relaxation time upon loweringthe density has also been noted in the relaxation time of the rotation of bulky molecules (ORDB) in Langmuir-Blodgett mon01ayers.l~~ According to Anfinrud et al. this decrease is due to a decrease in frictional drag between neighboring molecules as the surface pressure is lowered.107 The increase in relaxation time at area/molecule greater than 21 A2 is due to the presence of molecular tilt.60 (105) Bloor, D.; Bonsor, D. H.; Batchelder, D. N. Mol. Phys. 1977,34, 939. (106) Doucet, J.; Dianoux, A. J. J . Chem. Phys. 1984, 81, 5043. (107)Anfinrud, P. A,; Hart, D. E.; Struve, W. S. J.Phys. Chem. 1488, 92,4067.
Figure 8. Relaxation time T-, of the time autocorrelationfunction y ( t ) as a function of area/molecule.
A rotator phase is usually defined as having a correlation function that decays to zero. Following this definition, the monolayer in the density range of 18.5-25 &/molecule is in a rotator phase. However, the relaxation times are clearly different at different densities. Therefore, they cannot be in the same rotator phase a t all densities. The molecular reorientation a t 20-22 &/molecule has the same relaxation time as that of the molecular rotation of n-alkanes in their free rotator phase. By analogy to crystal alkanes, the monolayer a t 20-22 &/molecule must be also in a free rotator phase. Since the molecular reorientation at 18.5 A2/molecule has a relaxation time approximately 2.5 times higher than that a t 21 &/molecule, the monolayer at 18.5 A2/molecule cannot be in the same rotator state and must be in a restricted rotator or large oscillation phase.lo6 Dihedral Angle Conformation. In the short chain monolayers it was found that, between 18.5 and 20 A2/ molecule, the conformation of the molecules undergoes a change from all-trans to a conformation with kink In order to investigate the presence or absence of such a transition in the long chain monolayer, trans fractions are calculated for each dihedral angle in the amphiphilic molecule. In total, the molecule studied has 17 dihedral angles. These are numbered such that the dihedral angle formed by the head group and the three following methylenes is numbered 2 and the dihedral angle formed by the end methyl segment and the three preceding methylenes is numbered 18. A dihedral angle is considered to be in a trans position whenever cos 4; I -0.5, where 4i is the dihedral angle. Figure 9 shows that at 18.5 &/molecule most of the gauche defects are concentrated at the molecular ends especially near the molecular tails, in agreement with previous simulation^.^^@*^^ The conformation in the central section of the molecules (dihedral angles 4-14) is completely trans. At 19 A2/moleculethe chain changes to a conformation with gauche defects at dihedral angles numbered 7 and 9. A close analysis of these gauche defects shows that they are localized in only one molecule, and they are parts of a kink defect (twogauche angles separated by a trans angle). Figure 10 shows a close-up snapshot of
Langmuir, Vol. 9, No. 5, 1993 1341
Simulations of Amphiphilic Molecules
Ij
0.99
0.98
+18.5A2
i
t .
0.96 0
.
.
' 5
.
" e
.
IO
15
20
Dihedral Angle Number
Figure 9. Probability of finding a particular dihedral angle in the trans conformation at 18.5 and 19 A2/molecule. the monolayer showing the presence of this kink defect. From a packing point of view, kink defects are the most preferred form of chain deformations because the volume created by a kink defected molecule is only slightly higher than the volume of an all-trans and the overall direction of the chain is preserved. The formation of kink defectsis a well-known phenomenon in the phase behavior of long chain n-alkane crystalsW-lo1and biomembranes.lO8 Head Group Movement Normal to the Interface. In addition to the translational diffusion of molecules in the monolayer plane, we looked at the movement of head groups normal to the monolayer plane. A power spectrum has been used here to describe this movement. The power spectrum I&) for the movement of the head groups in the direction normal to the monolayer surface is
I&) = IZ,(P)l2
(14)
where z,(p) is the Fourier transform of the position of the head group in the z direction. The unnormalized spectra at all molecular areas are shown in Figure 11. The spectra displa two regimes. In the first regime (between 18.5and 19 2/molecule)the spectra do not show any clear peaks but indicate that frequencies of the order of 1.5 X loll s-1 or less are most prominent and equally probable within statistical errors. In the second regime (molecularareas greater than or equal to 20 A2/molecule) all spectra show a single dominating frequency at 1.5 X 1011 s-1. The power spectra results show that there is a transition in the monolayer closely related to the movements of the head groups normal to the interface. The relation of this change in vertical movements to other phase transitions is not well understood, but we speculate that this change is related to the decrease in the relaxationtime of molecular rotation around the major molecular axes, or in other words, the decrease in friction between neighboring molecules. The observationof three transitions that occur in the same density range is in agreement with the conclusions of Ginzburg and Manevitch,logwho contend that a phase transition coupled to molecular transition in the direction of the chain axis could occur simultaneously,
l
(108) Pink, D.A.; Green, T. J.; Chapman, D. Biochemistry 1980,19, 349. (109) Ginzburg, V. V.; Manevitch, L. I. Colloid Polym. Sci. 1991,269, 867.
at large chain length, with two transitions that are related to chain rotation and chain conformation.lN The spectra displayed in the short chain monolayers69 are quantitatively similar to the spectra shown in the long chain monolayer. Monolayer Tilt. The collective tilt of amphiphilic molecules in Langmuir monolayers has been the subject of numerous p u b l i c a t i ~ n s ~ ~and J " ~monolayer ~ properties such as tilt angle and tilt direction are still being debated; however, most groups seem to agree on the presence of a transition from a nontilted monolayer at high density to a tilted monolayer at low density. This tilt transition is believed to occur for area coverages between 19 and 21 A2/molecule. From our lateral pressure calculations it is impossible to distinguish any kinks, which are usually indicative of a second-ordertilt transition; however, from the tilt angle results (Figure 12)one can see the presence of a large change in tilt angle between 20 and 21 A2/molecule. The tilt direction between 21 and 22 A2/ molecule is not well-defined. There is a coexistencein the same monolayers of molecules that are tilted toward their nearest neighbors (NN) and others that are tilted toward their next-nearest neighbors (NNN).70For area coverages between 23 and 24A2,all molecules are tilted toward their nearest neighbors.70 The tilt angles calculated in the simulations are in good agreement with experimentally measured tilt angles,12 especiallyat high densities. As the density is lowered, the difference between the simulationand experimentalresults becomes more noticeable. This difference is most probably due to the presence of conformational defects in the model monolayers. It is perhaps more important to note that the relaxation of the monolayer collective tilt is found to be dependent on chain conformation equilibration and on uniaxial expansion or compression. Monolayer simulations at 42 "Cshow that molecules respondvery quicklyto expansion by assuming a large tilt angle. Then, the tilt angle decreases as the chain conformation is equilibrated. Infrared external reflection spectroscopic studies of Buontempo and Rice on a monolayer of heneicosanolshow two relaxation times: one fast relaxation due to uniaxial compression and a slow relaxation due to chain conformation equilibration.ll0 Conclusions Using molecular dynamics simulations the phase behavior of a model Langmuir monolayer has been investigated in the regime near close packing. The simulation results indicate the presence of a first-order phase transition between 24 and 25 &/molecule. The phase transition is characterized by a lattice break-up and a decrease in trans fraction from 95 to 90% . The monolayer goes from a state where molecules are tilted over their nearest neighborsto a state where molecules are no longer ordered. In contrast with the phase transition in the long chain monolayer, the transitions in the short chain monolayers were found to be only a break-up of lattice structure, with no significant decrease in trans fraction. The phase transitions in the short chain monolayers occur between 21 and 23 A2/molecule, and at a significant lateral p r e ~ s u r e , 6 while ~ * ~ ~ in the long chain monolayer the transition occurs between 24 and 25 A2/moleculeand at =O lateral pressure. In addition to the main phase transition, a second continuous transition has been detected between 18.5 and (110)Buontempo, J. T., Rice, S. A. Appl. Spectrosc. 1992, 46, 725. Buontempo, J. T.; Rice, S.A. Two manuscripts in preparation.
Karaborni
1342 Langmuir, Vol. 9, No.5, 1993
I
4
"1 L"
w
c
\
Figure 10. Snapshot of the monolayer at 19 A2/moleculeshowing gauche defects. Head groups are displayed in yellow and CH2 pseudoatoms are displayed in blue. The kink defect in the middle of the monolayer is shown in red. Gauche defects around the monolayer top area are shown in green. 1 8.98
b"
pf,
CA
8.88
1
18
5
19
28
21
22
23
24
25
26
(A') Figure 12. Mean monolayer tilt angles (open circles) and mean fraction of torsional angles in a trans state (filled squares) plotted as a function of molecular area. The continuous line represents tilt angles from X-ray reflectivity measurements in monolayers of C19H&OOH.'2 Area/molecule
5 lo1'
0 1
Frequency (sec- )
Figure 11. Power spectra calculated using equation 14 a t 18.525 A2/molecule.
20 A2/molecule. This transition is a compound of three transformations that occur simultaneously: (1)a significant decrease in the relaxationtime of the autocorrelation function that measures the reorientation of molecules around their major axes, (2) a change in chain conformation from an all-trans to a gauche defected conformation, and
(3)a change in the frequencyof the head group movements normal to the monolayer plane. We are presently concluding an extensive study on the effect of temperature on the phase behavior of amphiphilic molecules. Preliminary results at 0 "C show that solid islands are formed at low densities. The average lattice spacing of 4.0 A (obtained from radially averaged structure factors) inside these islands is close to that found at 22 A2/molecule;supportingthe idea that a heterogeneous island's phase is the preferred structure at low densities
Simulatione of Amphiphilic Molecules
Langmuir, Vol. 9, No. 5, 1993 1343
(>22 A2/mo1ecule).llo This view is further endorsed by Monte Carlo calculationslll which clearly indicate that at sufficiently low densities the monolayers prefer heterogeneous structures over ordered structures with large tilt angles. The orders and locations of phase transitions in molecular simulations are usually subject to finite size and boundary condition effecta.l12 The most obvious effect of finite size is to blur the distinction between first order and
continuous while periodic boundary conditions tend to stabilizethe solid.111 It is therefore necessary in the future to simulate bigger systems to verify the orders and the locations of the phase transitions described here. In summary, the resulta point to a complicated phase behavior that is dominated by several factors including lattice structure, chain conformation, molecular reorientation, tilt, and molecular translation in the monolayer plane.
(111) Siepmann, J. I., McDonald, I. R. Mol. Phys. 1992,75,255; Mol. Phys., in press. (112) Toxvaerd, S. Phys. Rev. A 1981,24,2735; Phys. Rev. Lett. 1983,
Acknowledgment. The author is grateful for the helpful discussionswith 0. Mouritaen, 0. H. Olsen, S. A. Rice, J. I. Siepmann, and S. Toxvaerd.
51, 1971.
