J. Phys. Chem. 1995, 99, 1393-1402
1393
ARTICLES Molecular Dynamics Simulation of Tetradecyltrimethylamonium Bromide Monolayers at the AirNater Interface Mounir Tarek, Douglas J. Tobias, and Michael L. Klein* Department of Chemistry, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 Received: August 18, 1994; In Final Form: November 9, 1994@
Molecular dynamics simulations of tetradecyltrimethylaonium bromide (C14TAB) monolayers at the air1 water interface were performed at surface areas of 45 and 67 A2/molecule. The computed density profiles normal to the interface for the headgroups and hydrocarbon tails are in good agreement with those derived from neutron reflectivity data. The simulations predict a decrease in the widths of the headgroup and the hydrocarbon chain distributions with increasing surface aredmolecule. Analysis of the distributions of the nitrogen atoms within the monolayers by radial and three-dimensional distribution functions showed that the headgroups form a disordered layer. The large widths of the headgroup density profiles seem to be related to an intrinsic roughness of the monolayers and not to a long-wavelength periodic structure. The headgroup radial distribution functions display split first peaks whose inner portion corresponds to about one neighbor at a distance of 6.5 A. The hydrocarbon chain structure was analyzed by taking into account the headgroup out-of-plane position. This analysis showed that the chains are highly disordered (2-3 gauche conformations/ chain) and that their internal structure depends on the vertical position with respect to the interface. A slight difference observed between the chain conformations in the two systems is related to the different widths of the headgroup distributions. The analysis of the distribution functions of water molecules and bromide counterions with respect to the headgroups indicated that the coordination shell structures around the headgroups are similar in the two systems. The in-plane diffusion constants for the centers-of-mass of the c14TAB molecules were of the same order of magnitude as those in other lipidwater systems.
Introduction In contrast to the situation with dense, solid-like (near closepacked) states of monolayers of amphiphilic molecules at the airlwater interface where there is a fair amount of detailed experimental structural data available,' relatively little information is available about the structure of the liquid-like states. This is unfortunate since it is the liquid state that most closely corresponds to the biologically relevant phase of lipid membranes2 Aside from the study of the surface pressure versus the averaged area/molecule isotherms, the main experimental insight into the properties of expanded films has been achieved by specular reflection of neutrons3 and X-raysS4 These techniques yield information on the density profiles of the surfactant perpendicular to the interface. In the case of neutrons, isotopic labeling of the molecules also allows one to obtain profiles for different fragments of the molecules comprising the monola~er.~,~fj Systems containing C,-trimethylammonium bromide (C,,TAB) surfactants in water have been the subject of a series of experimental structural investigations. Indeed, Penfold, Thomas, and colleague^^,^-^ have used neutron reflectivity measurements to determine the profiles of these molecules adsorbed at low densities at the water/air interface. Two types of structural parameters were estimated by fitting model profiles to the reflectivity data: the separations between the centers of the distributions of the distinct molecular fragments and the widths of these distributions. While the separations are determined @
Abstract published in Advance ACS Abstracts, January 1, 1995.
0022-365419512099-1393$09.00/0
rather accurately, the widths are usually influenced by the roughness of the interface, which arises from both static disorder (static roughness) and thermal fluctuations (capillary waves). After accounting for the surface roughness due to capillary waves, the results of Simister et ala3v7suggest that the thickness of the c14TAB monolayer at 44 P/molecule is less than the fully extended chain length (19.2 A). Furthermore, this thickness and the width of the headgroup density profile decrease with increasing surface area. In principle, molecular dynamics (MD) simulations can complement the experimental studies, furthering our understanding of the surfactant structure and behavior at the molecular level. In this paper, we report MD simulations of tetradecyltrimethylammonium bromide (C14TAB) monolayers at the air1 water interface. The aim of the present work is, first, to demonstrate the extent to which simulations employing current generation force fields can reproduce the experimental results and, second, to refme the low-resolution picture that has emerged from experimental density profile analyses by studying details of the structure and dynamics of the monolayer that are not yet available experimentally. We performed simulations of two expanded liquid monolayers of C14TAB with surface areas of 45 and 67 A21molecule, the areas estimated by Simister et aL5 for solutions with concentrations of 4.5 x and 1 x M, respectively. In the next section we will describe the setup of the systems and the simulation methods employed. In the following section the simulation and experimental results for the density profiles will be compared in order to assess the reliability of the simulations. Then, we will proceed to develop a microscopic picture of the 0 1995 American Chemical Society
1394 J. Phys. Chem., Vol. 99, No. 5, 1995
Tarek et al.
TABLE 1: Details of the Two Simulations of C14TAB Monolayer9 area/ molecule concentration box dimensions system (8.2) (10-3 ~ ) b (8.3) I
I1
45 67
4.5 1.0
37.5 x 37.5 x 80 46.3 x 46.3 x 80
no. Nn of water length molecules (ps) 1117 1685
428 377
Each simulation contained 64 C14TAT.3 molecules (2 monolayers, 32 moleculesflayer). Concentrations studied by Simister et al. (ref 5). The surface aredmolecule values used in the simulations are those estimated by Simister et al. for these concentrations. surfactant behavior as a function of surface density at the air/ water interface. In particular, we will focus on the association of the surfactant molecules with each other, the interactions of the headgroups with the solvent and the counterions, the hydrocarbon chain conformations, and lateral diffusion constants.
Simulation Setup and Methodology There are several choices that need to be made when setting up the system and specifying the boundary conditions in simulations of monolayers of surfactant/water systems.25 One possibility is to simulate a single monolayer in contact with the aqeuous solution and confine the solution by placing a hard wall below the interface.17 But this requires using a relatively large bath of water molecules for a given number of surfactant molecules to avoid artifacts due to the long-ranged structuring of water by the Another possibility is to simulate two monolayers on opposite sides of a slab of aqueous solution, with the thickness of the slab chosen sufficiently large that the monolayers are effectively isolated (Le., there is a region of bulk solution between the monolayers). In either case, there is the additional question of how the long-ranged electrostatic interactions will be treated. The simplest approach is to use spherical truncation,17but this seems inappropriate for the highly polar systems of interest here. In principle, an Ewald summation method suitable for roughly two-dimensional systems with a finite extent in the third dimension could be employed, but in order to be computationally feasible, the dimension in the nonperiodic direction should be much less than in the periodic directions, and this criterion is not met by the system sizes chosen in this Given these considerations, we have opted to use the two monolayers on a slab arrangement and use three-dimensional periodic boundary conditions and an Ewald summation with a large value for the dimension of the simulation cell in the direction normal to the interface (z). The assumption here is that the approximate mirror symmetry combined with the large z dimension results in weak (insignificant) interactions between periodic replicas in the z direct i ~ n The . ~ ~agreement between experimental and simulation results presented below for the density profiles suggests that this assumption is reasonable. The simulation cell dimensions and other details of the simulations are given in Table 1. The initial configurations were set up by arranging 32 tetradecyltrimethylammonium ions in a monolayer with the headgroup nitrogen atoms on a 4 x 4 body-centered square lattice in the xy plane and the hydrocarbons chains extending perpendicular to the lattice in all-trans conformations. The lattice constants we;e chosen to give surface aredmolecule values of either 45 A2 (system I) or 67 A2 (system 11). Then, for each system, two monolayers were placed such that their headgroups were in contact with a roughly 25 8, thick slab of water with appropriate x and y dimensions. The thickness of the water layer was just large enough to give a small region of bulk solution in the middle of the system. Finally, 64 water molecules were selected at random and replaced with bromide ions.
The potential energy function used in the calculations was of the CHARMM form.8 We employed the “polar-H’ model where only hydrogen atoms attached to oxygen are explicitly represented and the nonpolar hydrogens are implicitly represented through the heavy atoms to which they are attached, Le., the methylene and methyl groups are considered a single “extended atom”. This approximation should be adequate for the expanded (disordered) films studied here.g Except for the charges and van der Waals parameters for the tetraalkylammonium groups, which were taken from Jorgensen and Gao,l0 the ammonium ions were modeled using the CHARMM PARAM19 parameters” (the PARAM19 nonbonded interactions are essentially those of the TIPS potential28). We used the SPCEE model for water12 and the parameters of Lybrand et al.13 for the bromide ions. The Ewald method14 was used to compute the electrostatic interactions. The minimum image conventionI4 was employed to calculate the van der Waals interactions and the real-space part of the Ewald sum with simple truncation at 10 A. In both systems I and I1 the energy of the system was first minimized with respect to the bromide and water positions using steepest descent minimization. The position of all but the hydrocarbon chain atoms were then fixed, and 2 ps of dynamics were performed at 1000 K to randomize the chain conformations. The resulting configurations were then used to initiate MD simulations performed at a constant temperature of 293 K by using the NosC-Hoover chains method15 with separate thermostats coupled to the ammonium ions, bromide ions, and water molecules. The equations of motion were integrated using a Verlet-like algorithmI6 with a time step of 2 fs. The SHAKE algorithm’6was used to constrain the lengths of bonds involving hydrogen atoms. The fictitious masses of the thermostat variables were chosen according to the prescription of Martyna et al.,15 with time scales of 0.5 ps and a NosC-Hoover chain length of five. The lengths of the MD simulation runs are given in Table 1; 2089 and 1828 configurations were saved for the analysis of systems I and 11, respectively. Although most of the various components of the force field have been tested independently and exhibited reasonable agreement with experimental results for structural and thermodynamic properties of water,12 alkylammonium chloride solutions,lo and polar-alkyl liquids,28there is no guarantee that they will work equally well when combined in the system under study here. One way to evaluate the parameter set is to analyze the average components of the pressure tensor. If the potential model is appropriate, the average component of the pressure tensor in the direction normal to the interface (z) should be essentially zero (1 atm), the components in the plane of the interface (xy) should be related to the experimentally derived surface tension, and the off-diagonal components should be zero. It is wellknown that simulation estimates of the surface tension for systems containing ions and hydrogen-bonded molecules are subject to large errors. At the surface aredmolecule values studied here, it suffices to say that the in-plane pressure components should be “small”. The average components of the pressure tensor in the x, y, and z directions were -14 f 282, -44 f 280, and 7 f 333 atm, respectively, for system I, and -83 f 207, -84 & 223, and -14 f 232 atm, respectively, for system 11. In both systems the average magnitudes of the off-diagonal components were less than 7 atm. Thus the magnitudes of the pressure components are indeed “small”. We conclude that, insofar as the pressures are concemed, the potential is reasonable for the systems studied here. The agreement between the experimental and simulation results presented below for the structural parameters of the monolayers will provide further verification of the potential model.
J. Phys. Chem., Vol. 99, No. 5, 1995 1395
Tetradecyltrimethylammonium Bromide Monolayers Results and Discussion
1. Overall Monolayer Organization. The number density profiles normal to the interface of the water molecules, the bromide counterions, the polar (CH3)3N headgroups, and the hydrocarbon chain carbons are shown in Figure 1. The main features are in good agreement with the well-known picture of the surfactant organization at the &/water interface: the headgroups are localized in the interface, the hydrocarbon tails are mostly excluded from the water layer, and the bromide counterions are largely associated with the cationic headgroups but form a broad distribution which does not vanish in the bulk solution in the middle of each system. These features can be observed in the snapshots of instantaneous configurations of the two systems at the ends of the simulation runs displayed in Figure 2. Given the experimental concentrations corresponding to our systems and the numbers of water molecules in our simulations (Table l), we expect to have about 0.1 and 0.03 ammonium ion in the bulk solution in systems I and n, respectively. Indeed, we occasionally see an ammonium ion venturing out into the bulk in system I but not in system I1 (Figure 2). Finally, we note that the density profiles from both simulations are quite symmetric, indicating that the systems are well-equilibrated in the sense that the same overall structures are observed in the two (independent) monolayers comprising each system. 2. Comparison with Experiment. In order to compare our results with the density profiles obtained by fitting to the specular reflectivity data,5 the headgroup density was assumed to include the bromide ions, and all of the density profiles were averaged over both of the monolayers in each system. These resulting profiles are shown in Figure 3. The water distributions are essentially identical, while the molecular fragment profiles are narrower and shifted toward the interface for the less-dense system 11. The profiles are well represented by the functional forms assumed by the experimentalist^,^ namely, a tanh profile for the water (s) distribution and Gaussian profiles for the headgroup (h) and chain (c) distributions:
The parameters 6 and ai in the above equations are the widths of the distributions, and zp their centers ( i = h, c). The monolayer structure is described in terms of these widths as well as the separations, d c h (chain-headgroup) and d h s (headgroup-solvent), between the centers of these distributions. The parameters determined from the best fits to the simulation data are given in Table 2. Before comparing our results to the experimentally derived parameters, we note that, in the neutron reflectivity experiments, the observed widths of the density profiles include contributions from both the static and the thermal (capillary wave) roughness. The contribution from the latter depends to some extent on the resolution cutoff of the experiment and has been shown to be important for the reported d a h 3 Bocker et al." recently carried out two simulations of a C16TAC1 monolayer at 45 8,2/molecule with simulation cell xy dimensions of 29 x 25 A2 (32 molecules/ layer) and 58 x 50 A2 (64moleculesflayer), and they noted a wave-like structure with a 29 8, periodicity in the large system but none in the small system. It is not possible to tell from their account whether the observed wave is static or dynamic as they only presented an analysis of a single configuration. In
n
m
0.03
" , 0.02 n N
W
Q
0.01
0
0.03
0.02 n
N Q
W
0.01
0
Figure 1. Number density profiles in the direction ( z ) noma1 to the interface for the water molecules (solid lines), headgroups (dotted lines), hydrocarbon chains (long dashed lines), and bromide ions (short dashed lines): (a) system I (45 A2/molecule) and (b) system I1 (67 molecule).
Az/
the present study, perhaps because of the sizes of the systems simulated (cf. Table l), no wave-like arrangement of the monolayer was observed in either the high- or low-density systems. Although the wavelengths of capillary fluctuations in the real monolayer systems are not known, it appears that they are appreciably longer than our simulation cells (-50 A). X-ray scattering results sug est that capillary waves with wavelengths as small as 400 contribute to the roughness of the pure water-vapor interface, but whether or not shorter waves exist is not known due to experimental limitation^.^^ It seems reasonable that the simulated systems can be viewed as locally flat regions of a long-wavelength, periodic structure. Due to the absence of capillary waves in our simulations, the widths of our density profiles should not be compared directly to the experimental ones. Rather, the static and thermal contributions to the surface roughness should be separated. Such a separation has become part of the experimental data reduction and anal~sis.~ The thermal roughness contribution w to the experimentally observed density profiles Gobs (for both the headgroups and the chains) is defined by
1
UOb?
= uco;
+ w2
(3)
where a,,, is the corrected profile width.3 Lu et aL7 have determined w for the C16TAB monolayer at different surface areas, assuming that w equals the capillary wave value of the pure liquid at the corresponding surface tension. They estimated w to be 9.0 and 6.7 8, at 43 and 60 A2/molecule, respectively. For the C14TAB monolayer, w cannot be calculated directly since no surface tension measurements are available, but we assume that its value is similar to that of C16TAB at a given surface density. With this assumption, the corrected experimental full widths at half-maximum for the chain distributions are 10.9 & 0.8 and 8.3 f 1.5 A for C14TAB at 45 and 67 A2/ molecule, respectively, which are in good agreement with the results of the present study (Table 2). The separations, dv, do not contain any roughness contributions. As shown in Table 2, the simulation correctly reproduces
Figure 2. Final monolayer configurations from the simulations: (left) system I (45 A*/molecule) and (right) system I1 (67 A2/molecule).The atom coloring scheme is N, green; C, black; Br, yellow; 0, red; H, gray. The differences in the layer thickness and number of solvated bromide ions between the two systems should be noted.
Tetradecyltrimethylammonium Bromide Monolayers
J. Phys. Chem., Vol. 99, No. 5, 1995 1397
n
n
&
0.02
n
0 0.01 Q
n
h M
!-
Y
0
10
20
30
40
n I
04 I
2
v
n
0
N
4) Figure 4. Headgrou N N radial distribution functions, gm(i-), (solid line) system I (45 dP2hn:lecule) and (dotted line) system 11 (67 Az/ molecule).
W
Q
n n
04
N
2
v
n
N
v
Q
Figure 3. Number density profiles averaged over the two monolayers in each simulation (dotted lines) and fits (solid lines) of (a) water (solid line fit to tanh), (b) headgroups (solid lines fit to Gaussians), and (c) hydrocarbon chains (solid lines fit to Gaussians). The headgroup and chain distributions are shifted toward the water bulk (z = 0) for system II (recall Figure 2). the experimental values for system I at 45 8,*/molecule. For system I1 at 67 8,z/molecule, the neutron reflectivity profiles did not allow a determination of the structural parameters of the monolayer. We note, however, that the decrease in the separation between the headgroup and solvent distributions observed in our simulations implies, as expected, an increase in the water penetration of the monolayer with increasing surface aredmolecule . 3. Molecular Picture. The averaged structural parameters of the monolayer obtained from the present simulations correlate and are in good agreement, in terms of solvent, headgroup, and chain density profiles, with the experimental results. It is, however, interesting to develop the molecular picture of the structure of these surfactants afforded by the simulations. In the following sections we will investigate the headgroup three-
dimensional structure and the chain conformations, comparing and contrasting the two systems. Then, we will proceed to an analysis of the distributions of the water molecules and the bromide counterions around the headgroups. a. Headgroups. The nitrogen-nitrogen radial distribution functions, gm(r),18 are shown in Figure 4. The overall shapes are similar for both systems: each displays a bimodal first peak with maxima around 6.5 and 8.5 8, and a second peak arround 15 8,. As the concentration of the surfactant increases, the coordination number at 6.5 8, increases, and the second shell maximum is shifted to 1 8, shorter distance. A high degree of disorder exists in both systems as can be seen by the width of the first (bimodal) peak and the height of the first minimum, the latter indicating an exchange between the two coordination shells. Integration of the g m ( r ) up to the first minimum (11.5 and 12.5 8,) gives approximatly 7 neighbors for the system I and 6 for system 11. These values are similar to, but smaller than, that reported for the C16TAC1 at 45 A2/molecule.17 The spatial disposition of the headgroup nitrogen atoms was analyzed in terms of a three-dimensional distribution function, gm(s,z), which depends on the difference between the nitrogen atom positions along the normal (z) and the distance between their projections on the xy plane (s) and accounts for the inhomogeneity in the direction normal to the interface. The derivation and details of the calculation of gm(s,z) are given in the Appendix. The g m ( s , z ) functions for the two systems are displayed Figure 5 as both 3-d mesh and contour plots. Very broad distributions are observed for both systems. The two coordination shells of g(r) correspond here to the regions beyond and above s = 11- 12 8,. For the low-density system, the first coordination shell has a maximum around (s,z) = (8,0), and extends to z = 5 8, (Figure 5b). The shape of this shell (s =
TABLE 2: Parameters from the Fitted Density Profiles for the N(CH3)JBr Headgroups (h), Hydrocarbon Chains (c), and Water (sy system (aredmolecule) a h (A) a c (A) 6 (4 d c h (A) dhs (A) 6,s (A) 1(45 '42)
this workb
8.6 f 0.3 10.2 f 0.3 5.6 f 0.3 7.2 f 1.0 0.3 f 1.0 7.5 f 1.0 expf 16i1 5 .O 6.0 f 1.5 l i l 7.0f 0.75 expt 11 f 1 II (67 Az) this workb 6.7 f 0.3 8.4 f 0.3 5.1 i0.3 6.5 i 1.0 -1.5 f 1.0 4.8 f 1.0 expf 12 f 2 expt ( ~ o r r ) ~ 8.0 i 1.5 The widths a (full width at half-maximum) and 6 are defined in eqs 1 and 2 in the text. The 6, are the separations between the distribution centers. Obtained from fits to averaged density profiles from MD simulations (Figure 3). Obtained from fits to neutron reflectivity data in ref 6 (lle widths for the headgroups and chains). Full widths at half-maximum, a,,,, after correction of the results from ref 6 by the capillary wave contribution (cf. eq 3 in the text).
Tarek et al.
1398 J. Phys. Chem., Vol. 99, No. 5, 1995
2.5
4 35 3 25 2 15 1
2 1.5 1
8O
05 z(A) 0
z (A)
05 0
228 182 137 0456
-
-
3.83 3.06 2.3 - . - - . 153 .
I lo
2
4
6
8
1
0
1
2
s(A)
2
4
6
8
1
0
1
OB5
-
2
s (A)
Figure 5. Distribution functions, g”(s,z), for the nitrogen atoms (a, left) system I(45 A2/molecule)and (b, right) system I1 (67 A*/molecule).z is the out-of-plane distance between the nitrogen atoms, and s is the in-plane distance.
constant) implies an out-of-plane fluctuation of the molecules. When the concentration increases, the shell extends toward larger z values (up to 7 A) and, as seen in the g m ( r ) distribution function (Figure 4),is shifted toward a shorter radial neighbor distance (6.5 A). The integration of g m ( r ) in Figure 4 up to 6.5 8, (the first half of the split first peak) gives approximately 1 for the system at 45 A*/molecule. To gain insight into the spatial distribution of this neighbor, local density i s o s u r f a ~ e s ~have ~-~~ been calculated up to this distance. Following Shelley,19the instantaneous positions of the atomic nuclei (which generally give too noisy results) were replaced by normalized Gaussian distributions, yielding a smooth effective local density with little loss of detail. A surface of the nuclear distributions was then obtained via linear interpolation for contours of a specified value. The result, displayed graphically in Figure 6, shows, as expected from the gW(s,z) plots, that the maximum density is located around z = 0, with a spread of about 1.5 A, which implies that the closest neighbor pairs are approximately in the same plane. To conclude the analysis of the headgroup structure, the lack of a peak at positive z values, the broadness of the gm(s,z) profiles, and the shapes of the local density isosurfaces of the nearest neighbor indicate that there is no evidence for a static “ordered staggering” of the headgroups within the concentration range studied. Rather, the headgroups form a “disordered” layer which spreads in the out-of-plane direction at high concentration. The increase of the headgroup-headgroup mean distance in the out-of-plane direction takes place simultaneously with a decrease of the in-plane mean distance which leads principally to an increase of the density profiles normal to the interface. The shapes of the radial distribution functions, g m ( r ) , show that the spreading process affects weakly the nitrogen-nitrogen mean distance which remains roughly the same (8.5 A) in the two systems. At the higher concentration (system I), there is about one neighbor at the shortest separation (6.5 A), but this neighbor does not appear to account for the increase of the
interface width as it is mainly located in the same xy plane as the reference headgroup. b. Hydrocarbon Chains. The chain density profiles shown in Figures 1 and 3 represent mean distributions over the monolayers. The mean out-of-plane length of the hydrocarbon chains cannot be directly estimated from these profiles as previously suggested? since there is a contribution from the headgroup distribution which should be accounted for in a proper analysis of the chain structure. The “corrected” chain profiles, obtained by taking the difference (for each molecule) between the positions of the headgroup nitrogen and hydrocarbon chain carbon atoms, are shown in Figure 7 . Furthermore, since the width of the headgroup region is large, one expects that the intemal structure of the chains will depend on their vertical position in the interface. To underline this, the density profiles were calculated for the “outer” and the “inner” parts of the layer, Le., for molecules whose headgroups are located away from or toward the bulk with respect to the center of the headgroup distribution, respectively. From Figure 7 it appears that, indeed, for each system, the intemal configuration of a hydrocarbon chain depends on its position. The molecules in the outer part of the interface are more compact but display wider distributions for all carbons. The increased width of the distributions of the first methylene group in the outer part shows a wider orientational distribution of the N-CHz vector. The distributions of the terminal methyl groups in the outer part are nearly symmetric and span almost up to the region of the first carbon. The mean out-of-plane length of the molecules belonging to outer part (taken as the region where 90% of the terminal group is located) is 13.5 and 12 8, for the systems I and 11, respectively. The distributions of the terminal fragments for the molecules belonging to the inner part of the monolayer are shifted to larger distances from the nitrogen atom. The resulting confoqation leads to mean lengths of 15.5 A for system I and 14.5 A for system 11. The average gauche populations plotted in Figure 8 show that there is a noticeable difference between chain conformations
J, Phys. Chem., Vol. 99, No. 5, I995 1399
Tetradecyltrimethylaonium Bromide Monolayers
n c)
0.02
s n
N
w
Q
0.01
0
0.025 0.02 n c3
5
0.015
n
N
w
Q
0.01
0.005 Figure 6. Two views of density isosurfaces for the nearest nitrogen neighbor to a nitrogen at the origin, i.e., at a distance less than 6.5 A, in system I: (top) in-plane view looking along the x-axis and (bottom) view looking along the z-axis. Two isosurfaces with a ratio of 1:2 are shown. The darker region is the higher density isosurface and corresponds to the first relative maximum in the split first peak of g m (r) in Figure 4.
in the two parts of the interface. A larger gauche population in the first bond is observed for the chains in the outer part, and this appears to be related to the orientation of the N-CH:! bond. The distributions in Figure 7 imply that this bond tends to be oriented parallel to the interface. The kink in the chain at the first bond compensates for this orientation and points the tail away from the interface. In summary, we describe the hydrocarbon chain region of the monolayer as consisting of two distinct regions: an outer region where there is more space and, hence, the chains are more disordered (have more gauche conformations)and an inner region where the local density is high and the chains are, therefore, more extended. The present analysis performed on the systems at low and high concentrations shows that the fluidlike monolayer state of the C14TAB surfactant is characterized, as expected, by a high degree of disorder in the hydrocarbon chains (2- 3 gauche conformationskhaio overall). We found only small differences in the structures of the hydrocarbon tails between the two concentrations and conclude that the observed
0
0
5
10
15
20
Figure 7. Number density profiles of the hydrocarbon chain individual methyl(ene) groups normal to the interface with respect to the headgroup-nitrogen position. The solid (dashed) lines are for the chains bonded to a nitro en in the inner (outer) part of the interfacial region: (a) system I (45 I2/molecule) and (b) system I1 (67 A2/molecule).
difference in the chain density profiles primarily reflects differences in the headgroup out-of-plane distributions. c. Headgroup Association with Water and Counterions. We begin our analysis of the association of the headgroups with the water molecules and the bromide ions by discussing the the N-0 and N-Br g(r) functions plotted in Figure 9. The bromide distributions in both systems are characterized by a large first peak in the g(r) with a maximum at 4.5 A, in agreement with the simulation of C16TAC1,17 and a smeared out second coordination shell in the region above 8 A. The N - 0 g(r) functions also display large first peaks centered around 4.5 A, but the second coordination shells are more pronounced and are at shorter distances (7.5 A) than in the N-Br functions. When the two systems are compared, the bromide association to the headgroup decreases and, as expected
Tarek et al.
1400 J. Phys. Chem., Vol. 99, No. 5, I995
40.3-t
I \
-I
0
5
0
10
3
....... inner
0
5
10
0
5
10
15
Bond
Figure 9. N - 0 (a) and N-Br (b) radial distribution functions: (solid lines) system I (45 A2/molecule) and (dotted lines) system 11 (67 A2/
0.3
molecule). Note that system 11 is characterized by a slightly larger hydration and slightly lower counterion coordination number.
Bond
Figure 8. Fractions of gauche conformations in the hydrocarbon chains as functions of carbon-carbon bond number (counting away from the headgroup). The dashed (solid) lines are for the chains bonded to a nitrogen in the inner (outer) part of the layer: (a) system I (45 A2/ molecule) and (b) system 11 (67 A2/molecule).
from the density profiles (Figure 3), the water molecule association increases with increasing concentration. The f i s t coordination shells of the N-Br and N-0 g ( r ) functions are mostly contained within half the mean N-N distance (8.5 A). The coordination numbers obtained by integrating g(r) to the first minimum are 17 and 18 water molecules and 2.5 and 1.6 bromide ions in systems I and 11, respectively. Considering together the numbers of headgroups and bromide ions in the first coordination shell, we conclude that the counterions are shared by the headgroups. The water coordination numbers found in our simulations agree with the number 18 found in the ClaTACl simu1ationl7 and are slightly less (due to steric restriction) than the number 20 found in a simulation of a tetramethylammonium ion in water.1° To gain further insight into the association of water molecules and counterions with the headgroups, the local density isosurfaces of the water oxygens and the bromide ions around the (CH&NCHz headgroup were calculated and are shown in Figure 10 for system I. The C3" symmetry of the headgroup was considered, and the data were averaged over the entire simulation range. For clarity only the coordination shells up to the f i s t minima in the radial distribution functions are displayed. Both isosurfaces are characterized by rings around the CH3
groups that fuse along the C3 axis of the headgroup. The density isosurfaces of the bromide anions have holes in the N-CH3 bond directions, reflecting the preference for a bromide ion to be shared between positively charged methyl groups rather than associated with a single methyl group. The closer approach of the water oxygen atoms reflects their smaller size with respect to the bromide ions. The main differences between the two distributions appear to be the extension above the headgroup for the bromide distribution and the location of the density maxima: the water distribution is centered around the CH3 groups, Le., slightly below the headgroup N, while the bromide ion distribution is centered slightly above the nitrogen atom around the charged methylene group. The overlap of the water and the bromide density distributions around the headgroup indicates that the water molecules within the first coordination shell are simultanously solvating the trimethylammonium groups and the ions. The present results are in broad agreement with the C16TACl simulation results of Backer et al.17 Using a three-dimensional distribution function, they have found analogous solvation shells for the water and the counterions around the headgroup. In the present work we have employed the density isosurfaces to give further details on the N - 0 and N-Br association. 4. Lateral Diffusion. To characterize the in-plane dynamics of the ammonium ions in the monolayer, we have calculated the center-of-mass mean-squared displacements (MSDs) in the xy plane as functions of time. If the in-plane (lateral) motion is diffusive (with diffusion constant D),the MSD obeys the Einstein re1ation:l8
-
lim (s2(t)) = 4Dt
(4)
The center-of-mass lateral MSDs shown in Figure 11 for the tetradecyltrimethylammonium ions in the two systems both display linear behavior, indicative of diffusive motion in the plane of the monolayer, at long times. The slopes of these functions indicate, not surprisingly, that the lateral mobility of the molecules increases with increasing surface area per molecule. The diffusions constants, D, calculated from eq 4 are (1.9 f 0.1) x cm2 s-l and (4.8 f 0.1) x cm2 s-l for systems I (45 A*/molecule) and I1 (67 A2/molecule), respectively. There have been no experimental measurements of the diffusion constants in C,TAB systems, but the calculated
J. Phys. Chem., Vol. 99, No. 5, 1995 1401
Tetradecyltrimethylaonium Bromide Monolayers
Figure 10. Density isosurfaces around the headgroup (CH2N(CH3)3) of the water oxygens (green) and the bromide ions (red). Two isosurfaces with a ratio of 1:2 are shown, with the darker region corresponding to the higher density. The headgroup atoms are colored black. Note the holes in the density isosurfaces in the directions of the three N-CH3 bonds.
diffusion constants are of the same order of magnitude as those found in simulations of other amphiphilic molecules in monolayers,21bilayers,22and a Newton black film.23 Note that (s2(200 ps))lI2 = 4.5 and 6.5 A, respectively, for systems I and 11. Given that the radius of the first coordination shells of the headgroups is about 10 A, the ammonium ions have, on the average, hardly changed place at the ends of the simulations (428 and 377 ps, respectively, for systems I and 11).
Summary and Conclusions MD simulations of C14TAB monolayers at the aidwater interface were performed at two surface areas: 45 and 67 A2/ molecule. The density profiles normal to the interface are in good agreement with those fit to neutron reflectivity data. The profiles predict an increase of the headgroup and the hydrocarbon chain widths with increasing surface area. Due to the small system sizes, no long-wavelength capillary wave fluctuations were allowed, and thus the large widths of the headgroup density profiles are related to an intrinsic roughness of the monolayers. The distributions of the nitrogen atoms within the monolayers were analyzed by radial distribution functions and threedimensional distribution functions that accounted for the inhomogeneity of the system in the direction normal to the interface. The analysis showed that there is no static roughness within the layer and that the headgroups form a disordered layer whose thickness increases with decreasing surface area per molecule. The headgroup radial distribution function for the system at 45 A2/molecule displays a split first peak whose inner portion corresponds to about one neighbor at a distance of 6.5 A. This may be related to a tendency of the molecules to form dimers at high (near close-packing of the headgroups) densities which probably contributes to increasing disorder at the monolayer interface.
(u n
5. n A 3 (v W
v)
V
?* 0
50
100
150
200
Figure 11. Lateral mean-squared displacements of the centers-of-mass of the tetradecyltrimethylammonium ions as functions of time: (solid line) system I (45 W2/molecule) and (dotted line) system I1 (67 A*/ molecule).
The hydrocarbon chain structure was analyzed by taking into account the headgroup out-of-plane position. This analysis showed that the chains are highly disordered (2-3 gauche conformations/molecule) and that their internal structure depends on the molecular position with respect to the interface. The slight difference between the chain conformations in the two systems is related to the variations in the width of the headgroup distribution. The analysis of the distribution functions of the water molecules and the bromide counterions with the headgroups
Tarek et al.
1402 J. Phys. Chem., Vol. 99, No. 5, 1995
have shown that the coordination shell structures around the headgroup are similar in the two systems. As the concentration increases, the “hydration” decreases and the association with the counterions increases slightly. Finally, the in-plane diffusion constants for the centers-ofmass of the C14TAB molecules were found to be of the same order of magnitude as those in other lipidwater systems and were found to be in a ratio of 1:2.5 for the C14TAB monolayers with surface areas of 45 and 67 Az/molecule, respectively. In conclusion, the present simulations, carried out at two concentrations, enabled us to predict many structural properties of the C14TAB monolayers and to investigate their concentration dependence. Given the lack of high-resolution experimental data, the findings presented herein provide new insights into the fluid-like states of monolayers at the air water interface and give support to the assumptions made in fitting model density profiles to experimental neutron reflectivity data.
Acknowledgment. This work was supported by generous grants from the National Institutes of Health (R01 GM40712 and F32 GM14463) and the Procter and Gamble Company. We thank Bob Thomas for sending us experimental results prior to publication, Bob Laughlin for his interest and for discussions regarding this work, and John Shelley for invaluable advice concerning the graphics. Appendix In this Appendix we derive the gm(s,z) function used in this paper to describe the headgroup structure. The starting point is a distribution function, g(sij,zi,z,), similar to a function introduced by Snook and Henderson24to describe the structure of an inhomogeneous system, e.g., a hard sphere fluid near a wall. These distribution functions depend on the out-of-plane positions zi and z, (z is the direction of the inhomogeniety, normal to the interface) and the in-plane separation si, (xu2 = (xi (yi - Y,)~)of two atoms (in our case, headgroup nitrogen atoms) i and j . Following Snook and Henderson,24 we define the function g(sij,zi,z,) as
+
where (N(sij,zi,z,)) is the average number of pairs with out-ofplane positions in the ranges zi f dz and z, f dz, and in plane separation in the range so f ds. The normalization constant C is determined by
(N(sij,zi,zj))=
e(zi)e(zj) g(sij,zi,zj)dsi ds,
(A.2)
where e(zi)is the density of the nitrogen atoms at z,, and the integration volumes are the slabs of width dz (centered at zi and zj) and the cylindrical shell of radius sij and thickness ds (centered on atomj). If ds and dz are sufficiently small such that the integrand is constant, eq A.2 can be rewritten as
As we are interested in seeing if there is a static roughness of the headgroups within the interface, we define a restricted correlation function, gm(s,z), for pairs of nitrogen atoms with a separation z = Izi - zjl in the direction normal to the interface and an in-plane separation s gw(s,z) =
(A.4)
where the integration is over all atoms i and the atoms j at distances z = Izi - zjl. Combining eqs A.3 and A.4, we obtain a computationally convenient equation for g(s,z) in terms of (N(sy,zi,zj)) and e(&functions which are readily computed from simulation data:
Note that care must be taken to consider the density e(zj) at
z, satisfying (zi - zjl = z. References and Notes (1) Kenn, R. M.; Bohm, C.; Bibo, A. M.; Peterson, I. R.; Mohwald, H.; Als-Nielsen, J.; Kjaer, K. J. Phys. Chem. 1991, 95, 2092 and references therein. (2) Small, D. M. Handbook of Lipid Research. 4. The Physical Chemistry ofLipids; Plenum Press: New York, 1986. (3) Lu, J. R.; Simister, E. A,; Thomas, R. K.; Penfold, J. J. Phys. Condens. Matter 1994, 6, A403. (4) Btlorgey, 0.;Benattar, J. J. Phys. Rev. Lett. 1991, 66, 313. (5) Simister, E. A,; Lee, E. M.; Thomas, R. K.; Penfold, J. J. Phys. Chem. 1992, 96, 1373. (6) Lu. J. R.: Simister. E. A,: Thomas. R. K.: Penfold. J. J. Phvs. Chem. 1993,97, 6024. 17) Lu. J. R.: Hromadova. M.: Simister, E. A,: Thomas, R. K.; Penfold, J. J. ‘phys. Chem. 1994, 98, 11519. (8) Brooks, B. R.; Bruccoleri, R.; Olafson, B.; States, D. J.; Swaminathan, s.;Karplus, M. J. Comput. Chem. 1983, 4, 187. (9) Bareman, J. P.; Klein, M. L. J. Phys. Chem. 1990, 94, 5202. (10) Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986, 90, 2174. (11) Reiher, W. E. Ph.D. Thesis, Harvard University: Cambridge, 1985. (12) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (13) Lybrand, T. P.; Ghosh, I.; McCammon, J. A. J. Am. Chem. SOC. 1985, 107, 7793. (14) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1989. (15) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (16) Ciccotti, G.; Ryckaert, J. P. Comput. Phys. Rep. 1986, 4, 345. (17) Blicker, J.; Schlenkrich, M.; Bopp, P.; B r i c k ” , J. J. Phys. Chem. 1992, 96, 9915. (18) Friedman, H. L. A Course in Statistical Mechanics; Prentice-Hall: Englewood Cliffs, 1985. (19) Shelley, J. C. Ph.D. Thesis, University of Pennsylvania, Philadelphia, 1993. (20) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1993, 9, 916. (21) Ahlstrom, P.; Berendsen, H. J. C. J. Phys. Chem. 1993, 97, 13691. (22) Egberts, E.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 3718. (23) Gamba, Z.; Hautman, J.; Shelley, J. C.; Klein, M. L. Langmuir 1992, 8, 3155. (24) Snook, I. K.; Henderson, D. J. Chem. Phys. 1978, 68, 2134. (25) Bareman, J. P. Ph.D. Thesis, University of Pennsylvania, Philadelphia, 1992. (26) Lee. C. Y.: McCammon. J. A.: Rosskv, P. J. J. Chem. Phvs. 1984, 80, 4448. (27) Hautman, J.; Klein, M. L. Mol. Phys. 1992, 75, 379. (28) Jorgensen, W. L.; Ibrahim, M. J. Am. Chem. SOC.1981,103,3976. (29) Schwartz, D. K.; Schlossman, M. L.; Kawamoto, E. H.; Kellogg, G. J.; Pershan, P. S.; Ocko, B. M. Phys. Rev. A 1990, 41, 5687. ~
where A is the surface area of the system.
g(sij,zi,zj)@(zi>d(z - Izi - z,l)A dz
JP9422320
