Langmuir 1999, 15, 3589-3594
3589
Molecular Dynamics Simulation of the Vibrational Spectra of Stearic Acid Monolayers at the Air/Water Interface Emiko Okamura,*,† Nobuhiro Fukushima,‡ and Soichi Hayashi‡ Institute for Chemical Research, Kyoto University, Uji, Kyoto 611-0011, Japan, and Nihon Silicon Graphics K.K., Yebisu Garden Place 4-20-3 Ebisu, Shibuya-ku, Tokyo 150-6031, Japan Received July 13, 1998. In Final Form: February 26, 1999 Infrared (IR) spectra of the Langmuir monolayers of stearic acid were simulated by the molecular dynamics (MD) method, for the first time, where the translational, rotational, and vibrational degrees of freedom of all atoms in the monolayer were considered. The IR spectra simulated by this method completely reproduced surface area dependence of the monolayer structure and the chain axis orientation. The characteristics of the spectra were also consistent with the chain orientation evaluated by in situ experimental techniques. In addition, the simulated spectra gave evidence of no area dependence of the head group orientation in the monolayer.
Introduction Recently, molecular dynamics (MD) simulation has been widely attempted with Langmuir monolayers, which are amphiphilic monomolecular films spread at an air/water interface.1-6 The structure and molecular orientation of these monolayers can be predicted theoretically by this method. The results can also be compared to the numbers of experimental observations. The purpose of the present work is to simulate the IR spectra of the Langmuir monolayers, for the first time, by the MD calculation. In this study the MD method is applied to evaluate the structure and molecular orientation of the amphiphiles theoretically from the simulated IR spectra, and to compare the results with those of the IR and some other in situ experiments. We deal with the MD simulation of a Langmuir monolayer of stearic acid, simple but forming various phases in a two-dimensional system. Here, we use the potential model including the effective bond charges and their fluxes with respect to certain internal coordinates, which are used also as the infrared intensity parameters. We consider the translational, rotational, and vibrational degrees of freedom of all atoms of molecules in the monolayer. These have not been considered yet in this system. The molecular orientation of the long-chain fatty acid monolayer at the water interface has also been experimentally determined by different kinds of physicochemical techniques.7,8 Therefore, we finally compare the simulated IR spectra with these results. Molecular Dynamics Simulation Method of Simulation. Here, we consider the translational, rotational, and vibrational degrees of freedom of * To whom correspondence should be addressed. † Kyoto University. ‡ Nihon Silicon Graphics K.K. (1) Harris, J.; Rice, S. A. J. Chem. Phys. 1988, 89, 5898. (2) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 96, 5505. (3) Karaborni, S.; Toxvaerd, S. J. Chem. Phys. 1992, 97, 5876. (4) Karaborni, S.; Toxvaerd, S.; Olsen, O. H. J. Phys. Chem. 1992, 96, 4965. (5) Alper, H. E.; Bassolino, D.; Stouch, T. R. J. Chem. Phys. 1993, 98, 1. (6) van Buuren, A. R.; de Vlieg, J.; Berendsen, H. J. C. Langmuir 1995, 11, 2957. (7) Kjaer, K.; Als-Nielsen, J.; Helm, C. A.; Tippman-Krayer, P.; Mo¨hwald, H. J. Phys. Chem. 1989, 93, 3200. (8) Sakai, H.; Umemura, J. Bull. Chem. Soc. Jpn. 1997, 70, 1027.
all molecules in the system. The force acting on an atom i ()Fi) is assumed to be the sum of the intra- and intermolecular force, Fiintra and Fiinter, respectively. Thus, for a system consisting of N′ atoms,
Fi ) Fiintra + Fiinter ) -∂Vintra/∂ri - ∂Vinter/∂ri, i ) 1, 2, ..., N′ (1) where V is the corresponding potential energy. We use the intramolecular potential function expressed by a set of internal coordinates St as9
Vintra )
1
1
∑ftt′StSt′ + 6tt′t′′ ∑gtt′t′′′StSt′Stt′t′′ + 2 tt′ 1
∑ htt′t′′t′′′StSt′Stt′t′′Stt′t′′t′′′,
24tt′t′′t′′′
t,t′,t′′,t′′′ ) 1, 2, ..., 3N′-6 (2)
where St’s represent the bond stretching, the bend angles, the out-of-plane angles, and the torsion angles. Potentials. As the quadratic valence force constants ftt′ in eq 2 for stearic acids, we adopted the values used by the molecular mechanics simulation of fatty acids.10 Further, we used the combination of 2-, 3-, and 6-fold symmetrical torsional potentials, considering the conformational change.10 For the bond stretching, we used the anharmonic Morse potential as given by
VMorse )
∑Di exp{-ai (ri - r0i )} ×
[exp{-ai(ri - r0i )} - 2] (3)
where Di is the dissociation energy, ai the dissociation parameter, ri the distance between atoms, and r0i the equilibrium value of the bond length. The force constant f ′ is given by f ′ ) 2ai2Di in this case. According to the parameters ai and Di by Yokoyama et al.,10 the force constants were obtained as Table 1. For the water, we adopted quartic and cubic force constants in addition to the quadratic ones. We used the (9) Ishioka, T.; Murotani, S.; Kanesaka, I.; Hayashi, S. J. Chem. Phys. 1995, 103, 1999. (10) Yokoyama, I.; Miwa, Y.; Machida, K.; Umemura, J.; Hayashi, S. Bull. Chem. Soc. Jpn. 1993, 66, 400.
10.1021/la9808723 CCC: $18.00 © 1999 American Chemical Society Published on Web 04/15/1999
3590 Langmuir, Vol. 15, No. 10, 1999
Okamura et al.
Table 1. Force Constants Used for the Bond Stretchinga H-O C-O CdO H-CR H-Cβ H-C C-CR CR-Cβ (CH2CH2) C-C (CH3CH2) C-C (CH2CH2) a
7.103623 6.043547 11.726251 4.702529 4.940328 4.839276 5.017210 4.595265 4.479989 4.574709
Force constants are in the unit of aJ‚Å-2.
Table 2. Equilibrium Atomic Charges (Electron Unit) of Stearic Acids H(OH) O(OH) O(CdO) C(CdO) C(CRH2) C(CβH2) C(CH2) C(CH3) H(CRH2) H(CH2) H(CH3)
0.3346 -0.3539 -0.3187 0.2000 -0.0300 -0.0520 -0.1000 -0.1200 0.0600 0.0500 0.0400
anharmonic force constants of water,11 fitting well with the experimental vibrational spectra. The intermolecular potential energy Vinter is assumed as the sum of the Coulomb interaction and the van der Waals interaction. The Coulombic potential VC is given by
VC )
1
qiqj
∑ 2 i,j r
(4)
ij
Here, qi and qj are the instantaneous point charge on the atoms i and j, where the equilibrium value and the charge fluxes are taken into account and rij is the distance between the ith and jth atoms. The bond charge parameters and the bond charge fluxes were transferred from the values of fatty acids.10 The equilibrium atomic charges of stearic acids are listed in Table 2. The van der Waals potential VvdW is assumed to be the Buckingham-type of eq 5 having the attractive constant A and the repulsive constants B and C as
VvdW )
Figure 1. A model of a stearic acid Langmuir monolayer used in the simulation.
{(-Ar-6 ∑ ij + B exp(-Crij)} i,j
(5)
The cutoff distance was 14 Å for all atom-atom pairs. We adopted potentials reported by Stillinger and Rahman12 for O‚‚‚O, O‚‚‚H, and H‚‚‚H. For the other atomatom pairs, we used the values of fatty acids.10 Model of Simulation. For convenience, we use a symmetric system consisting of five layers (air/stearic acid/ water/stearic acid/air) shown in Figure 1, as a model of a stearic acid Langmuir monolayer. We consider two layers of stearic acids with the thickness of 25 Å lying upon and under the water layer (the thickness, 30 Å). As starting conditions, we introduce the following two points: (1) the density of the water corresponds to that of the bulk liquid, while that of the stearic acids is equal to the surface density obtained by the surface pressure (π)-area (A) isotherm (11) Smith, D. F., Jr.; Overend, J. Spectrochim. Acta 1972, 28A, 471. (12) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978, 68, 666.
Figure 2. Initial configurations of three runs of the simulation. Runs I, II, and III correspond to the surface area of the stearic acid monolayer of 22.5, 25.3, and 33.7 Å2 molecule-1, respectively.
of the monolayer, and (2) stearic acids do not interact with each other in the z-direction, which is normal to the monolayer surface. To study the surface pressure dependence, we examine three runs I, II, and III, each unit cell with its dimension 14.22 × 14.22 × 25 Å3 containing 9, 8, and 6 molecules of stearic acids, respectively. The runs I, II, and III correspond to the surface area of 22.5, 25.3, and 33.7 Å2 molecule-1 in the π-A isotherm.13 Initial configurations of the runs are depicted in Figure 2. For water, 200 molecules exist in the unit cell (14.22 × 14.22 × 30 Å3). Conditions of the MD simulations are summarized in Table 3. The water phase with the initial condition (14.22 × 14.22 × 30 Å3), where no monolayers are spreading, was also simulated separately, to characterize the interfacial water near the monolayer as compared to the bulk water. For given initial values of the position of atom i ()ri, the positional vector represented by a set of Cartesian coordinates) and the velocity ()vi), the Newton’s equation of motion was solved stepwise for a given time step ∆t (13) Sakai, H.; Umemura, J. Chem. Lett. 1993, 2167.
Stearic Acid Monolayers at the Air/Water Interface
Langmuir, Vol. 15, No. 10, 1999 3591
Table 3. Conditions of the MD Simulation of a Stearic Acid Monolayer on the Water Surfacea surface area/Å2 molecule-1 no. of molecules stearic acid water a
I
II
III
22.5
25.3
33.7
18 200
16 200
12 200
Temperature, 293 K; step time, 0.244 fs; N, V, T const.
(2-12 ps ) 0.244 fs) at the given energy and density by a modified Verlet method.14 The time step was chosen for the convenience of the computational program of the Fourier transformation in the IR spectral simulation described below. The initial values of the translational velocities, the angular momenta of the molecules, and the amplitudes of each normal vibrations were chosen from normal random numbers based on the Maxwell-Boltzmann distribution.9 The number of the molecules, N, the volume of the unit cell, V, and the temperature, T ) 293 K, were always constant in each simulation. An equilibration run was performed within 100 ps, from the judgment of the total energy and the temperature of the system. Data analyses were based on the runs of the last 93, 86, and 139 ps in equilibrium for I, II, and III, respectively. The simulation time of the bulk water was 290 ps. The calculation was made by a CRAY Y-MP2E computer at the Supercomputer Laboratory, Institute for Chemical Research, Kyoto University and a CRAY Y-MP8E computer of SGI/Cray Research Co. Ltd., U.S. Infrared Spectra. We can use the linear response theory15 to calculate infrared absolute intensities from the ensemble average of the dipole moment time correlation function 〈µd(0)‚µd(t)〉 where d ) x, y, or z. We obtained the absorption cross section, Rd(ω), corresponding to the infrared absolute intensities of both intra- and intermolecular vibrations, as
Rd(ω) ) Id(ω) )
4π2ω[1 - exp(-βpω)] Id(ω) 3pcn
(6)
∫-∞∞dt exp(-iωt)〈µd(0)‚µd(t)〉
(7)
1 2π
where ω is the angular frequency, β ) (kBT)-1 in which kB is the Boltzmann constant and T the temperature, p ) h/2π in h is Planck’s constant, c the velocity of light, n the averaged refractive index of the Langmuir monolayer ()1.5), and Id(ω) the absorption line shape. Here, we consider polarization effects in terms of n in eq 6. Results and Discussion Monolayer Structure. The structure of stearic acid Langmuir monolayers can be easily visualized by the simulation. Snapshots of the monolayers on the water surfaces are given in Figure 3, for three runs, I, II, and III, corresponding to the surface area of the monolayer 22.5, 25.3, and 33.7 Å2 molecule-1, respectively. Each figure clearly illustrates the stable and well-organized monolayers forming at the air-water interface. Distribution profiles of atom positions projected to the surface normal z also provide a lot of information about the area-dependent monolayer structure. The result is presented in Figure 4 for stearic acid carbonyl O and chain methyl C atoms, together with the position of the water O atom, where the abscissa is in accord with the five-layered model in Figure (14) Beeman, D. J. Comput. Phys. 1976, 20, 130. (15) Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1981, 74, 4872.
1. The water layer exists from 60 to 90 Å, indicating the appropriateness of the model used. It is obvious that there is no significant area dependence of the position of the carbonyl O (head group of the stearic acid), always lying at the air-water interface. In contrast, the position of the terminal methyl C remarkably depends on the surface area; its position is far away from the water surface with a decrease of the surface area. This result means that the monolayer thickness is obviously increased upon the monolayer compression. Since the positions of the head groups are the same at each area, as mentioned above, it is demonstrated that the increase in the monolayer thickness is directly combined with the structural change in the chain region. Therefore, the variation of the thickness in Figure 4 predicts the different chain conformation and/or the different orientation of the chain axis upon the monolayer compression. These are discussed in the following sections. It should be noted that the distribution profiles are not strictly similar between two layers of stearic acid, presumably due to the statistical error. Hydrogen Bonding and the Water Structure. The head group of the stearic acid monolayer can interact with water and/or the neighboring stearic acid molecules via the hydrogen bonding. It is therefore considered that features concerning the hydrogen bonding of the monolayer head group, as well as the interfacial water, should be changed upon the compression of the monolayer. How many numbers of hydrogen bonds (H-bonds) can be formed in the respective surface area of the monolayer? What is the surface area dependence of the interfacial water structure? To make these clear, we estimate the probability of the O‚‚‚H distance smaller than 2.5 Å which is accounted for as an H-bond, for two kinds of O molecules of stearic acid and one water. In Figure 5A-C the number of H-bonds per molecule formed in the stearic acid hydroxyl O, carbonyl O, and the water O, respectively, are shown. It is found that the head group of the stearic acid molecules forms H-bonds irrespective of the surface area (compare the head group position of stearic acid in Figure 4 with Figure 5A,B). This means that the stearic acid head groups interact with each other and/or with the interfacial water via H-bonds. More H-bonds are formed in the stearic acid carbonyl dO (Figure 5B) than the hydroxyl -O- (Figure 5A). This is due to the large space remaining free around the carbonyl dO, in which the intermolecular H-bond can be easily formed. In contrast, the formation of the intermolecular H-bond is restricted at the hydroxyl -Osite, because of the intramolecular chemical bond of O-H already existing. One interesting and suggestive thing in Figure 5A,B is that at 22.5 Å2 molecule-1 H-bonds are formed over the wide range, especially z of 42-50 and 100-108 Å. This is also the case in water at the same area, as shown in Figure 5C. In these regions of z, it is certain that very small amounts of the head groups of stearic acid are actually distributed, as represented by the two wings in the distribution profile (see the arrows in Figure 4). This is also the case for the profile of water (Figure 4). One can therefore find that water does not penetrate into the hydrocarbon chain region of the monolayer, but that H-bonds between the head groups and water are also formed in this region. Besides, these H-bonds can be at first explained by considering the zigzag orientation of the stearic acid molecules such as in the solid phase of the monolayer. Although the liquid-condensed phase is dominant at the surface area of 22.5 Å2 molecule-1, it is probable that the solid phase locally appears during the monolayer compression.
3592 Langmuir, Vol. 15, No. 10, 1999
Okamura et al.
Figure 3. Snapshots of a stearic acid Langmuir monolayer for the three MD runs I, II, and III.
On the other hand, each water molecule forms two H-bonds in the bulk phase (z of 60-90 Å), as shown in Figure 5C, irrespective of the presence of the monolayer. The interfacial water can be distinguished at the surface area of 22.5 Å2 molecule-1, also forming almost two H-bonds per molecule. Chain Conformation. Stearic acid molecules have 15 twist angles (φ). A twist angle is considered to be trans whenever φ > 120°, otherwise, to be a gauche conformation. Figure 6 gives the fraction of the gauche position against the carbons C1-C15. Here, the carbon number 1 corresponds to the twist angle between the bond angles ∠(C1-C2-C3) and ∠(C2-C3-C4), where C1 is the carbonyl carbon. The number 15 corresponds to the angle at the terminal methyl side. It is shown that the fraction of the gauche conformation decreases with decreasing the surface area at all methylene and methyl sites except for the methylene site near the head group. There are two sites having the large fraction of the gauche position; C1 near the head group region and C15 at the terminal methyl side. Further, no surface area dependence at the former site is recognized. These facts mean the trans position of the twist angle at this site is excluded at all areas, since the steric hindrance of the trans position is relatively large. On the other hand, the large gauche percentage at the
latter site is due to the large degrees of freedom at the terminal methyl group. This is also shown by the MD simulation of the Langmuir monolayers of single-chain amphiphiles having 15 methylene segments and a dipolar head group.2 Bond Orientation. The bond orientation function S is estimated for each bond of the hydrocarbon chain:
S ) (3 cos2 θ - 1)/2
(8)
where θ, usually called the orientation angle, is the one formed by the bond vector connecting the adjacent segment centers and the vector normal to the monolayer surface. In Figure 7 the surface area dependence of the S value at the respective bond of the chain is shown. The orientation function evidently increases with the surface area decrease (the monolayer compression) at all methylene and methyl sites of the chain. Although the area dependence of the orientation function is quite roughly obtained (only three points of the surface area) in this study, this behavior is completely similar to the previous MD simulation of the Langmuir monolayers.2 According to eq 8, the obtained change of S from 0.2 to 0.8 during the compression (Figure 7) corresponds to the change of an average orientation
Stearic Acid Monolayers at the Air/Water Interface
Figure 4. Distribution profiles of atom positions projected to the surface normal z for Langmuir monolayers at the surface areas of 22.5 (solid line), 25.3 (dashed line), and 33.7 (dasheddotted line) Å2 molecule-1. Thin and bold lines demonstrate the positions of chain methyl C and carbonyl O atoms of stearic acid at the corresponding surface area, respectively. The position of the O atom of water at the respective area is also given in the bottom.
Langmuir, Vol. 15, No. 10, 1999 3593
Figure 6. Gauche percentages as a function of a carbon atom number of a stearic acid at the surface areas of 22.5 (+), 25.3 (O), and 33.7 (b) Å2 molecule-1. The carbon atom numbers 1 and 15 correspond to the carbonyl and the terminal methyl sides of the stearic acid, respectively.
Figure 7. Orientation function S as a function of the respective bond of the hydrocarbon chain at the surface areas of 22.5 (+), 25.3 (O), and 33.7 (b) Å2 molecule-1. The atom numbers 1 and 15 correspond to the carbonyl and the terminal methyl sides of the stearic acid, respectively.
Figure 5. Number of H-bonds per molecule formed in the stearic acid hydroxyl O (A), the carbonyl O (B), and the water O(C) atoms against the position projected to the surface normal z at the surface areas of 22.5 (solid line), 25.3 (dashed line), and 33.7 (dashed-dotted line) Å2 molecule-1. In (C) the number of H-bonds per molecule formed in the bulk water phase is also given by the dotted line.
angle of the chain θ h from 47° to 21° against the surface normal. Thus, the θ h ’s obtained in this study are also close to the experimental values of arachidic acid Langmuir monolayers estimated by the X-ray reflectivity at the corresponding areas.7 Polarized Infrared Spectra. The polarized infrared spectra of stearic acid Langmuir monolayers simulated
by the MD method clearly show the orientational changes of the chain region of stearic acid caused by the change in surface area. Most typical is the antisymmetric (νas) and symmetric (νs) CH2 stretching vibrations in the 28003000 cm-1 region. The simulated spectra in this region are shown in Figure 8 as a function of the surface area. The xy-averaged components of both CH2 stretching bands gradually increase in the intensity with decreasing the surface area (the monolayer compression) (Figure 8A). The z components of the same bands, however, decrease in intensity (Figure 8B). The dichroic ratio, determined by Ixy/Iz, is increased from 1.5 ( 0.3 to 6.0 ( 1.0 for both bands during the compression. This behavior is ascribed to the two reasons overlapped: (i) the molecular chain axis of the stearic acid is oriented normally against the film surface and (ii) the fraction of the gauche conformers is decreased with increasing the surface density during the monolayer compression. Remember that the dipole moments of both CH2 stretching vibrations lie vertical to the chain axis of stearic acid. It is noteworthy that the area dependence of the simulated polarized spectra strongly supports the orientational change of the hydrocarbon chain axis of stearic acid. The simulated spectra completely reproduce the surface area dependence of the atom distribution profiles (Figure 4), the chain conformation (Figure 6), and the orientation function of the chain axis (Figure 7). Since the difference in the force constant between gauche and alltrans conformations10 is not considered in this study, the
3594 Langmuir, Vol. 15, No. 10, 1999
Okamura et al.
Figure 9. Simulated IR spectra (0-2000 cm-1 region) of a stearic acid Langmuir monolayer at the surface areas of 22.5 (solid line), 25.3 (dashed line), and 33.7 (dashed-dotted line) Å2 molecule-1. (A) xy-averaged; (B) z.
Conclusions
Figure 8. Simulated IR spectra (2600-3200 cm-1 region) of a stearic acid Langmuir monolayer at the surface areas of 22.5 (solid line), 25.3 (dashed line), and 33.7 (dashed-dotted line) Å2 molecule-1. (A) xy-averaged; (B) z.
frequency shift in both CH2 stretching modes is unfortunately not noticed with the change in the gauche/trans population. However, all results obtained verify the vertical orientation of the chain with the all-trans conformation upon the monolayer compression. The CdO stretching (νCdO) band around 1700 cm-1, on the other hand, exhibits almost no significant spectral change during the monolayer compression, as shown in Figures 9A,B. This means that the carbonyl group in the hydrophilic region always hydrates via the hydrogen bonding at the water interface, and this hydration prevents the direct orientational change of the carbonyl site during the compression. It is still probable that the degree of the hydration of the carbonyl group varies with the surface density change. However, almost no change in the νCdO region is a striking contrast to the orientational change of the CH2 chain region where the vertical orientation is stabilized by the interchain hydrophobic interaction. Finally, it should be noted that our theoretical prediction is especially valuable for the subject of the orientational change in the head group region of stearic acids during the monolayer compression. This is because there have been no IR experimental evidence of the head group orientation in the fatty acid Langmuir monolayer because of the strong IR absorption of water in the νCdO region.
Infrared spectra of stearic acid Langmuir monolayers are simulated for the first time, by considering the translational, rotational, and vibrational degrees of freedom of all atoms in the monolayer, as well as considering effective bond charges and their fluxes with respect to certain internal coordinates. The surface area dependences of the polarized spectra in the chain methylene stretching region and carbonyl stretching region clarify the area dependence of the chain axis orientation and no significant dependence of the orientation of the head group. Moreover, the surface area dependence of the simulated spectra in the methylene stretching region is verified by the area dependences of the distribution profiles of the atom positions, the gauche fraction, and the orientation function of the hydrocarbon chains. The limitation of the present simulation is the smaller dimension size. Much larger time and length scales are needed to solve several problems especially with longrange interaction. For example, the island formation at a low surface pressure region is considered to be one of the most important phenomena. We are now preparing the simulation of the island formation with larger system sizes, although further progress in the computational power is expected. Acknowledgment. Computation time was provided by the Supercomputer Laboratory, Institute for Chemical Research, Kyoto University and SGI/Cray Research Co. Ltd., U.S. The authors also thank Professor Masaru Nakahara of this laboratory for his constant support during this work. LA9808723
