Langmuir 1995,11,237-246
237
Molecular Dynamics Simulations of Dioctadecyldimethylammonium Chloride Monolayers D. B. Adolf Department of Physics, University of Leeds, Leeds L S 2 9JT, United Kingdom
D. J. Tildesley* Department of Chemistry, University of Southampton, Southampton SO17 lBJ, United Kingdom
M.R. S. Pinches and J. B. Kingdon Parallel Applications Center, University of Southampton, Southampton SO16 7NP, United Kingdom
T. Madden and A. Clark Unilever Research Port Sunlight Laboratory, Wirral L63 3JW, United Kingdom Received June 20, 1994. I n Final Form: September 8, 1994@ Molecular dynamics are performed on systems of two-dimensional periodicity composed of 64 ionic dioctadecyldimethylammonium chloride amphiphiles arranged in a monolayer at 298 K with surface coverages ranging from 57 k per amphiphile to 150 A2 per amphiphile. Bond lengths are constrained
whereas valence and torsional angles interactions are described by conventional expressions. Nonbonded interactions are introduced through an anisotropic united atom method due to Toxvaerd. The amphiphiles adhere to the surface through electrostatic and nonbonded interactions between the amphiphiles and an underlying substrate. Results are presented for the density normal to the monolayerplane, order parameters, and conformational transition rates. The structure of the two-dimensional layer at the lowest headgroup area studied is quite disordered. Small islands of empty surface surrounded by amphiphiles are formed. As the headgroup area increases, we observed further translational disordering and an increase in the number of amphiphiles aligned with the surface. The g+g+tg+g+defect observed in one of the chains is stable on a timescale of 250 ps. The equilibrium structures of the two chains are quite different but the conformational dynamics as observed by the transition rates are almost indistinguishable.
Introduction An understanding of the properties of self-assembling systems is a t the root of much research in biomolecular organization and cellular biophysics. Phospholipids spontaneously form aggregates when dispersed in water and are a principal component ofbiological membranes. Their structure is characterized by an ionic headgroup and two hydrocarbon tails. Singer and Nicolsonl review a fluid mosaic model where mitochondria and chloroplast membranes are represented as globular proteins intermixed within a phospholipid bilayer. Depending on the degree of hydration, phospholipids exhibit a range of lyotropic liquid crystalline structures. Heating the crystalline phase (L,) leads to a tilted (Lp)or untilted (Lp) lamellar gel phase. Further heating eventually affords the fluid phase (La).The hydrocarbon tails are tightly packed in the crystalline phase. Within the gel phase, the tails remain closely packed but concerted rotational motions of the chains increase. The L, phase has completely disordered hydrocarbon tails which exhibit fluid-like configurational dynamics. Experimental studies have addressed the assembly structures and their molecular constituents. The structure and dynamics within these assemblies can be monitored using fluorescence depolarization2 and electron spin resonance techniques3 which are sensitive to the motion of probe molecules. Similar @
Abstract published inAdvanceACSAbstracts, October 15,1994.
(1)Singer, S.J.; Nicholson, G. L. Science 1972,175,720. (2) (a) Zannoni, C. Mol. Phys. 1981,42,1303. (b) Johansson, L. B.-
A.; Lindblom, G. Q.Rev. Bwphys. 1980,13,63.
information is also available from nonprobe techniques such as X-ray diffraction4 and deuterium (D)NMR measurement^.^ The lamellar gel to liquid crystalline phase transition has been studied using models with varying degrees of torsional flexibility describing the hydrocarbon chains.6 The freedom to choose input potentials within molecular dynamics simulations affords an opportunity to observe the dependence of aggregate dynamics and structure on the details of molecular interactions. The combination of experiments, theory, and simulation serves to clarify the structure-property relationships governing self-assembly. To understand the chemical detail within self-assembly, considerable attention has been focused on the cationic dialkylamine salts represented by the general chemical formula ( C H ~ ) ~ N + [ ( C H Z ) , - ~ C H ~ ~ [ ( C H Z ) ~ -The ~CH~~C~-. two long alkyl chains and ionic headgroup make these molecules structurally similar to the phospholipids. Kunitake' shows that the length of each alkyl chain plays a crucial role in determining the aggregate's structure. (3)van Ginkel, G.; Korstanje, L. J.; van Langen, H.; Levine, Y. K. Faraday Discuss. Chem. SOC. 1988,81,49. (4)(a) Luzzati, V. Biological Membranes; Chapman, D., Eds.; Academic Press: New York, 1968;Vol. 1,p 71. (b) Cevec, G.; Marsh, D. Phospholipid Bilayers: Physical Principles and Models; Wiley-Interscience: New York, 1987. (5) (a)Seelig, J. &.Rev. Biophys. 1977,10,353.(b) Seelig, J.;Seelig, A. Q.Rev. Biophys. 1980,13,19. (c) Bocian, D.F.; Chan, S. I. Annu. Reu. Phys. Chem. 1978,29,307. (6)(a) Caill6, A,; Pink, D.; de Verteuil, F.; Zuckermann, M. J. Can. J.Phys. 1980,58,581. (b) Nagle, J. F. Annu. Rev. Phys. Chem. 1980, 31,157. (7)Kunitake, T.Angew. Chem., Int.Ed. Engl. 1992,31,709.
0743-746319512411-0237$09.00/0 0 1995 American Chemical Society
Adolf et al.
238 Langmuir, Vol. 11, No. 1, 1995
Dihedral Angle / '
a.)
38 39 Figure 1. Simulated form of dioctadecyldimethylammonium chloride. The circles represent united methylene (2-18 and 20-36) or united methyl (1and 37-39) atoms. The darkened circles possess acharge of+0.251el. Beads l-18comprise chain 1and beads 20-37 comprise chain 2. The backbone torsional angles are numbered from 4 to 37 with angle 4 involvingatoms 1through 4 and angle 37 involving beads 34 through 37.
This behavior has been studied for hydrocarbon chain lengths where n and m range from 1 to 22 and where the chains can be of different lengths (n * m). Our aim in studying these systems is to understand how the friction between surfaces can be modified using surfactant coatings of the amphiphile dioctadecyldimethylammonium chloride (DOAC1, m = n = 18). This molecule is one of the principal active constituents within commercial cationic fabric softeners. A necessary prerequisite to the study of these systems under shear is an understanding of the equilibrium behavior of adsorbed surfactant layers. In this paper, we present the results of molecular dynamics simulations of DOACl monolayers adsorbed on a dielectric continuum. The model is initially described in terms of the input potentials and starting configurations. The simulations are used to calculate the structural and dynamical properties of a monolayer. We present results for the molecular tilts, the density profiles normal to the surface plane, the orientational order parameters, and the conformational properties of the amphiphiles. The final section contains our conclusions.
Model The simulated structure ofthe DOACl molecule is shown in Figure 1. Each circle represents either a methylene (bead numbers 2-18 and 20-36) or a methyl united atom grouping (bead numbers 1 and 37-39). The numbering system begins at the terminal methyl united atom of alkyl chain 1 and proceeds down this chain through the headgroup and up alkyl chain 2. There are 39 simulated sites per amphiphiles plus 1 chloride counterion. The modeled amphiphiles interact through a variety of potentials which are now presented in detail; Table l summarizes the potential parameters. Intramolecular Potentials. Bond lengths are fixed at their equilibrium values for the duration of the simulation using constraint dynamics implemented through the SHAKE algorithm.8 Three-atom valence angles interactions are described by the potential
(8) (a) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J . Comput.
Phys. 1977,23,327.(b)Ciccotti,G.;Ryckaert,J.-P. Comput. Phys.Rep. 1986,4,345.
Dihedral Angle /
b.) Figure 2. (a) The torsional potential plots for the CCCC torsional angles within the simulation. The solid line is generated from an expression developed by Steele whereas the dashed line originates from the equation proposed by Ryckaert and Bellemans. (b) The torsional potential plots for the simulated torsional angles involving nitrogen. The solid line refers t o the CCCN torsional angles whereas the dashed line refers t o the CCNC types. These plots were generated using molecular mechanics on small molecules (see text).
where ke and 80 are 520 kJ mol-l and 113.3' for all the simulated valence angle^.^ Distinctions between CCN-, CNC-, and CCC-type valence angles are not made within the simulations reported here. The potential expression for four-atom torsional interactions is 6
where the parameters ai depend on the type of torsional angle. For quartets composed solely of methylene andor methyl beads, the Ryckaert and Bellemanslo (RBI or Steelell coefficients are used. Both parameter sets have been used in modeling n-butane. The RB coefficients are based on experimental datal2 while the Steele set is generated from ab initio calculations using a 4-31G basis set. The RB and Steele butane torsional potentials are shown in Figure 2a. For both these potentials the trans (t)conformation is at 0'. The local minimum at negative torsional angle values corresponds to the gauche minus (g-) conformation and the local minimum a t positive torsional angle values corresponds to the gauche plus (g+) conformation. The shapes of the potentials are qualitatively similar; however several differences are apparent. The local minima of the Steele potential are at f112.8', whereas the local minima of the RB potential are f120". The Steele potential has a t to g* barrier of 15.1 k J mol-l and the energy of the gauche states is 4.41 kJ mol-'. The (9) Clarke, J. H. R.; Brown, D. Mol. Simul. 1989,3, 27. (10)Ryckaert,J.-P.;Bellemans, A. Chem. Phys. Lett. 1975,30,123. (11)Steele, D.J . Chem. SOC.,Faraday Trans. 2 1985,81, 1077. (12)(a) Scott, R. A.; Scheraga, H. A. J . Chem. Phys. 1966,44,3054. (b) Woller, P. B.;Garbisch, E. W., Jr. J . A m . Chem. SOC.1972,94,5310.
MD of Surfactant Monolayers
Langmuir, Vol. 11, No. 1, 1995 239
Table 1. Potential Parameters (a) Bond Lengths
SHAKE convergencetolerance 1.0
10-4
constrained bond lengths (A) C-N c-c 1.45
1.53
(b) Valence Angle Interactions 520
113.3
torsional plot has a t to g* barrier height of 3.18 k J mol-' where the gauche states are located at f86"with an energy of9.12 kJ mol-l. The CCNC torsional plot has one barrier height of 18.9 kJ mol-' and three equivalent states of 0 k J mol-' located a t 0" and f120". Nonbonded Potentials. This class of potential accounts for the interaction between sites i andj which are separated by at least three atoms along the contour of the amphiphile and between sites on different molecules. This interaction is of the Lennard-Jones form
(c) Dihedral Angle Interactions (J mol-')
RB Steele CCCN CCNC
9287 8832 9250 5780
12 155 -13 119 -3059 26239 -31493 0 18087 4880 -31800 0 0 0 -3809 28570 10584 -53244 -27371 36316 28070 31710 -36540 -79300 -1 160 51410
For the interaction between unlike atoms the LennardJones length and energy parameters are given by
(d) Anisotropic Nonbonded Interactions CHdMe) mass (kg mol-') 15.0 x c7 (A) 3.527 E (J/mol) 997.8 Lennard-Jones cutoff
CHZ 14.0 x 3.527 665.2 8.818 A
N 14.0 x 3.250 711.3
c135.4 x 3.400 992.4
(e) Continuum Surface Interactions GSS (A) 3.40 EQS (J/mol) 232.8 e (atoms/&) 0.142
N Ymaxa
&axb
nmm del
dJL
(0 Electrostatic Interactions 320 real charges and 320 dielectric image charges 0 (for real charges interacting with real charges) 1 (for real charges interacting with dielectric image charges) 15 (for real charges interacting with real charges) 15 (for real charges interacting with dielectric image charges) 1 (eqs 6 and 7) -l(chloron) and +0.25 (CH2 or CH3 a to nitrogen) 0.08
aDenotes magnitude of lattice vector index in real space. Denotes magnitude of reciprocal lattice vector index in reciprocal space.
corresponding g* to t barrier height is therefore 10.7 k J mol-'. The RB potential has a t to g* barrier of 12.4 k J mol-l and the energy ofthe gauche states is 2.93 k J mol-l. The corresponding g* to t barrier height is therefore 9.47 k J mol-'. Furthermore, the energy ofthe syn conformation a t 4 = f180" is 27.4 k J mol-l using the Steele potential and 44.8 k J mol-' as calculated using the RB potential. Ab initio studies ofAllinger etaZ.13report the syn rotational barrier of butane as 21.0 k J mol-l which supports the lower syn barrier energy found in the Steele potential. At room temperature, the RB potential precludes direct transitions between the two gauche states. These transitions are observed infrequently when using the Steele potential. For torsional angles involving the headgroup nitrogen atom, molecular mechanics minimi~ationsl~ were performed on small molecules resembling fragments of the headgroup and fitted to eq 2 to afford parameterized ai values. n-Propyltrimethylammonium was the model compound for the CCCN torsional angle parameterization and n-ethyltrimethylammonium for the CCNC torsional angle parameterization. The best fit ai values are reported in Table l a and the corresponding torsional potential energy plots are illustrated in Figure 2b. The CCCN (13) Allinger, N. L.; Grev, R. S.;Yates, B. F.;Schaefer, H.F., I11 J . Am. Chem. SOC.1990,112, 114. (14) The results published were generated using the programs CHARMm and QUANTA. These programs have been developed by Molecular Simulations, Inc.
0..=
w
Eli
+
oii oj 2
=
For each bead i , the link-list method is used to search for all nonbonded interactions, j , with rg -= Rcut,where RCut is the cutoff distance. For these simulations, Rcutis 8.82 A, which corresponds to 2.5aMeMe. In eq 3, ru is normally measured from the center of one united atom, i , to its neighbor j indicating that the force site and mass site coincide. Simulation work on alkane crystals15and Langmuir-Blodgett film@ reveals that at close packing densities the shapes of the CH2 and CHB moieties are poorly represented as spherical united atoms. We use a technique due to Toxvaerdl' to model the anisotrophy of the methylene and methyl units while still retaining the computational savings inherent in the united atom approximation. This approach has been shown to reproduce the simulated equation of state behavior of liquid alkanes" better than the conventional isotropic united atom model. This approach moves the force sites to the geometric centers of the CH2 and CH3 groups and leaves mass sites at the nuclear positions. For a methylene united atom, the force site is moved toward its hydrogen atoms by a distance of 0.320 A along avector which bisects the HCH valence angle. The force site for a methyl united atom is shifted a similar distance toward its hydrogen atoms in a direction collinear to the bond vector between the methyl grouping and its neighboring carbon or nitrogen atom. The nitrogen atoms and chloride counterion are not modified in this way. Expressions for the corresponding forces on a particular mass site must be consistent with the change in potential and must compensate for the internal torques which are produced within this model. These modified force expressions can be found in ref 17. The 40 sites which comprise each simulated amphiphilel counterion pair interact with a surface through an integrated Lennard-Jones potential. For a particle at a distance z from the surface the potential is
The surface number density, e = 0.142 atoms A-3, is derived from crystal structures of cellulose18since cotton (15) Ryckaert, J.-P.; Klein, M. L.J . Chem. Phys. 1988, 85, 1613. (16) Moller, M. A.; Tildesley, D. J.; Kim, K. S.; Quirke, N. J.Chem. Phys. 1991,94, 8390. (17)Toxvaerd, S. J. Chem. Phys. 1990,93,4290. (18) Marsh, J. T.; Wood, F. C. An Introduction to the Chemisty of Cellulose; Chapman Hall: London, 1938; Chapter 21.
240 Langmuir, Vol. 11, No. 1, 1995
is one substrate of interest. Cellulose has a complicated molecular structure and we have used the E,, and a,, correspondingto a carbon atom.lg The cross parameters in eq 4 are calculated using the normal Lorentz-Berthelot mixing rules. ElectrostaticPotentials. The conventional Coulomb potential is used to represent the interaction between partial charges in different molecules. These interactions are long-ranged and are normally evaluated using Ewald summation techniques20 for systems which are periodic in three dimensions. In this paper, the simulation cell is periodic in the monolayer plane but not orthogonal to it. Within the all-image convention,the electrostatic potential for a system of N charges is given as21
T
* x-
a.)
I
0- - - -
with rg denoting the distance between charged species i and j . T represents a two-dimensionallattice vector, ( Y ~ , Y,), of the periodic system where r ~ ,= ,, 1 7 i - 7j - 31. The prime indicates that the i =j terms are skipped when 1 31 = 0. Two efficient methods have been proposed to compute electrostatic interactions within systems of two-dimensional periodicity. The first technique due to Hautman and Klein22(HK) adopts a modified Ewald approach. The potential is decomposed into a short-range interaction
and a long-range interaction
The long-range potential is computed in reciprocal space and the short-range potential in real space. sG is the inplane separation between charges i and j and ZG is the out-of-plane separation. h,(sg,a) represents an Ewald convergence function with convergence distance a. Hautman and Klein use the error function and associated derivatives to represent h,(sc,a). Effectively, an nth order Taylor series in powers of z ~ / s Cwith binomial coefficients a, has been added and subtracted to eq 5. The corresponding forces are determined by analytical differentiation of eqs 6 and 7. An alternative approach due to L e k n e ~ starts - ~ ~ with summations for the forces over all periodic cells. These summations are manipulated using the Poisson- Jacobi identity to obtain convergent results. The potential is obtained by integrating the equations for the force. The constant of integration is calculated using well-known results for perfect two-dimensional ionic lattices. We have developed algorithms for both approaches and a detailed comparison of their performance will be presented in a future publication. The runs within this paper use the (19) Kim, K. S. Ph.D. Thesis, Department of Chemistry, University of Southampton, SO9 5NH, cites e = 0.114 atoms A-3 for graphite. (20) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (21) Tildesley,D. J. Computer Simulation in Chemical Physics;Allen, M. P., Tildesley, .D. J., Eds.; Kluwer: Boston, MA, 1993; Chapter 2. (22) Hautman, J.; Klein, M. L. Mol. Phys. 1992,75, 379. (23) Lekner, J. Physica A 1991,176,485.
.
Adolf et al.
L b.)
-Electrostatic Surface - - - - - -
Figure 3. (a) Schematic representation of a starting configuration looking normal to the monolayer (xy)plane. The stick figures represent the cationic amphiphiles with g+g+tg+g+ defects. The small dots represent the chloride counterions. For a given number of amphiphiles, the surface coverage is altered by modifying the values of X and Y. (b) A portion of a typical starting configuration viewed in the monolayer (xy) plane. The circles above (+z) the physical surface represent united atoms where darkened circles are charged. Darkened circles within the physical surface represent dielectric image charges used to account for the polarization of the continuum substrate by the amphiphiles and counterions.
HK method which is significantly faster than the Lekner approach for charges in a monolayer. The ionic charges within the amphiphiles also polarize the 9-3 surface. To account for this polarization, dielectric image charges are created within the surface. If a charge is at a position x, y , z above the image plane, its corresponding dielectric image is at x, y , -z below the image plane (see Figure 3). In this paper, the image plane is located ( a ~ e ~ J=3 1.175 A) above the centers of the top atoms of the surface plane. The charge magnitude of a dielectric image, qimage, is related to the magnitude of the corresponding charge, q, by
where cr is the relative permittivity of cotton within this For cr = 7, qimage = -0.75q.25 The surface polarization is easily incorporated into the HK formalism. For an interaction between a charge and its dielectric image, half of the potential afforded by eqs 6 and 7 is apportioned to the dielectric image and half to the charge. Interactions between a charge and its own dielectric image are not skipped when Y = 0 in eq 5. All charges interact with all dielectric images in the central cell and with all periodic images of the dielectricimages. Dielectric images do not interact with each other and electrostatic interactions between charges on the same cationic amphiphile (24) Torrie, G. M.; Valleau, J. P.; Patey, G. N. J . Chem. Phys. 1982, 76, 4615. (25) Brandrup, J.;Immergut, E. H. Polymer Handbook, 2nd ed.; Wiley Interscience: New York, 1975.
MD of Surfactant Monolayers
Langmuir, Vol. 11, No. 1, 1995 241
are zero. To obtain a reasonable starting configuration, a balance must be found between the repulsion of the 9-3 surface interaction and the attraction between the charges and their dielectric images. If the image plane is located at 1.175 a suitable balance is obtained. If the height of the image plane is increased, the attraction between the headgroups and the surface decreases causing the molecules to leave the surface. If the height of the image plane is decreased, the chargeldielectricimage interactions overcome the surface repulsion and the molecules penetrate the surface in an unphysical manner. The HK method is most accurate in the limit of a thin film of charges where z&"
