UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
THERMODYNAMIC MODELING AND MOLECULAR DYNAMICS SIMULATION OF SURFACTANT MICELLES By ROBERT ANTHONY FARRELL A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1988 To Peggy, for her support, encouragement, and patience TABLE OF CONTENTS ABSTRACT . . . CHAPTERS 1 INTRODUCTION . . 2 A MOLECULAR THERMODYNAMIC MODEL OF MICELLE FORMATION . . 2.1 Background . . 2.2 Stoichiometry and Reaction Equilibria for Multicomponent Micelles . 2.3 Estimation of Free Energy Changes . 2.4 The Calculational Technique . 3 RESULTS OF THERMODYNAMIC MODELING . 3.1 Parameter Estimation . 3.2 Behavior of the Model for SingleComponent Nonionic Systems 3.3 Aggregate Size Distribution and Concentration Behavior . 3.4 Chain Length and Temperature Effects. 4 MOLECULAR DYNAMICS SIMULATION OF SURFACTANT MICELLES . . 4.1 Background . 4.2 The Molecular Dynamics Method . 4.3 The Model Surfactant Molecule . 4.4 The Model Micelle ......... 4.5 Summary of Computer Simulations . 5 RESULTS OF MOLECULAR DYNAMICS SIMULATION . 5.1 Mean Radial Positions of Groups . 5.2 Probability Distributions of Group Positions. 5.3 Conformations of Chain Molecules . 5.4 Shape Fluctuations . . 5.5 Pair Correlations of Groups . . 6 CONCLUSIONS . . . iii Page v 4 4 8 11 17 23 23 30 40 43 52 52 55 60 65 71 78 78 90 107 126 133 137 I * o . APPENDICES A MICELLE SIZE AND SHAPE . . B HEAD GROUP INTERACTION IN A BINARY MICELLE C PROGRAM LISTING FOR SINGLECOMPONENT NONIONIC MICELLE CALCULATION . D MOLECULAR DYNAMICS PROGRAM LISTING FOR SIMULATIONS 1 AND 2 . E DERIVATION OF HEAD GROUP INTRAMOLECULAR INTERACTIONS . F MOLECULAR DYNAMICS PROGRAM LISTING FOR SIMULATION 3 . . REFERENCES . . . BIOGRAPHICAL SKETCH . . . 143 . 150 . 160 . 174 * 218 . 224 . 268 . 276 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THERMODYNAMIC MODELING AND MOLECULAR DYNAMICS SIMULATION OF SURFACTANT MICELLES By ROBERT ANTHONY FARRELL April, 1988 Chairman: John P. O'Connell Major Department: Chemical Engineering The association of surfactant molecules into aggregates known as micelles gives them a broad range of applications. In spite of the widespread use of this class of chemicals, there is not yet sufficient scientific understanding to predict their behavior in solution. A molecular thermodynamic model has been developed to describe the formation of micelles in multicomponent surfactant solutions. Using a hypothetical, reversible sevenstep process, the total free energy of micellization is calculated by summing the contributions due to solvophobic interaction, mixing, surface formation, confor mational change, head group interactions and electrostatics. Distributions of micelle sizes and compositions can then be generated through a set of reaction equilibria. Where possible, the free energy contributions are related to comparable processes on which experimental measurements have been made. Aggregate size distributions have been generated from the model for singlecomponent solutions of nonionic surfactants of different chain lengths and at different temperatures and solution concentrations. It has been found that a detailed description of micelle structure and entropy effects on chain conformation is necessary to fully describe the thermodynamics of micelle formation without empirical parameterization. To this end, computer simulations of model micelles have been conducted by the molecular dynamics method. Micelles of three different head group characteristics and a comparable hydrocarbon droplet have been simulated. A spherical shell is used to contain the aggregate, providing estimated solvophobic interactions with the molecules. The simulation results reveal that internal structure of the aggregates is relatively insensitive to the head group characteristics with the greatest effect resulting from head group size. While the micelles all showed chain ordering and the hydrocarbon droplet did not, the bond conformations averaged approximately 71 percent trans in all cases. Results of other simulations and experimental studies are generally similar to those of the present work for the effects of chain length, aggregate size, and simulation technique on static properties. vii CHAPTER 1 INTRODUCTION There are few classes of chemical species which have received more attention in the scientific literature, or have exhibited a more ubiquitous presence in everyday life, than the amphiphilic molecules known as surfactants. Their unique solution properties of association and adsorption at interfaces give rise to a broad range of applications. Long the essential component in detergency applications, in more recent times surfactants have assumed roles of importance in applications as diverse as enhanced oil recovery and pharmacy. The association of surfactant molecules in solution into aggregates known as micelles is the primary attribute which has garnered interest among the scientific community. Since the discovery of micelles in solution (McBain and Salmon, 1920), many experimental and theoretical studies have been conducted, yet full understanding of these systems of molecules has not been achieved. True predictive capability is not yet a reality. This investigation takes a twofold approach toward the ability to make quantitative predictions of the behavior of systems of surfactants in solution. Since the behavior of the surfactant solution is a consequence of its thermodynam ics, a model set in the framework of molecular thermodynamics is developed to describe the formation of micelles in a multicomponent surfactant solution. A thermodynamic description of a micellar system is limited by a precise knowledge of the structure of micelles. To this end, computer simulations of model micelles by the molecular dynamics method are conducted. By accurately modeling the forces present, pertinent descriptions of micelle structure are obtained. In Chapter 2 of this work, further background on micellar systems is given and the development of a model for the free energy change upon formation of micelles in a multicomponent surfactant solution is described. A multi step reversible process is employed to generate contribu tions to the total free energy change due to hydrophobic interaction, mixing effects, conformational change, head group interaction, and electrostatic interaction. Some brief results obtained with the thermodynamic model are presented in Chapter 3. Distributions of free energy change of micellization and aggregate size are calculated for surfactants of different chain lengths and for solutions of different temperature. Due to limitations in the scope of this portion of the investigation, calculations were carried out only for cases of singlecomponent solutions of nonionic surfactants. The computer simulation of surfactant micelles is described in Chapter 4. Following a concise description of the molecular dynamics method, its application to the simulation of micelles is discussed. A summary of the computer simulations of the four modelsthree micelles and one hydrocarbon dropletis given. In Chapter 5, the results of analyses of the computer simulations are presented. Elements of aggregate internal structure, chain conformations, aggregate shapes, and changes in the aggregate with time are investigated for the four model aggregates. Chapter 6 provides a summary of the significant conclusions of this work. Recommendations are made toward the future progress of both of the projects described in the previous chapters. CHAPTER 2 A MOLECULAR THERMODYNAMIC MODEL OF MICELLE FORMATION 2.1 Background The forces present in liquid solutions dictate that solution of a polar solute in a polar solvent is more favored than is a polar solute in a nonpolar solvent. Similarly, the nonpolar solvent is more accommodating toward a nonpolar solute than a polar one. Therefore, molecules which contain both polar and nonpolar groups exhibit a unique behavior when present in a polar or nonpolar solvent. Such molecules, known as surfactants, will tend to minimize the unfavored contact (i.e., polarnonpolar) while maximiz ing the favored contact (i.e., polarpolar). At an interface between polar and nonpolar liquids, the surfactants will penetrate the interface to achieve "like" interactions on both sides, reducing the "unlike" interactions encountered in the bulk liquid. The polar groups of a large number of surfactant molecules packed closely at an interface will repel each other, producing a spreading pressure which reduces the interfacial tension. In the bulk liquid, surfactant molecules will aggregate into structures known as micelles which can afford much the same benefits as the interface. In the typical case, a surfactant with a polar "head group" and a nonpolar "tail," when present in sufficient quantity in a polar solvent such as water, will form micelles having an interior consisting of tails and possibly some nonpolar solubilizate and a surface comprised mostly of head groups. The head groups remain in contact with the watera favored interac tionwhile the tails reduce their contact with the wateran unfavored interaction. When the solvent is nonpolar, inverted micelles can form. The ability of surfactants to adsorb at interfaces and aggregate into micelles makes them a very useful class of compounds. The reduction of interfacial tension has many applications, ranging from oil recovery to biological processes. In addition, micelles can solubilize other solutes in their interior, as in drug delivery processes, and reactions can even take place there, as in emulsion polymerization. The earliest and best known use of surfactants, detergency, uses both aspects of their behavior. As useful as these phenomena are, they are not understood to a degree that would allow their full potential to be realized. The ability to predict behavior rather than just explain it is the goal of this undertaking. This requires a knowledge of the thermodynamics of surfactant phenomena. The formation of surfactant aggregates in solution instills a certain ambiguity in the description of the system by a traditional thermodynamic formalism. The aggregation of surfactant monomers into micelle structures has been treated as the formation of a "phase" (Blankschtein et al., 1985; Kamrath and Franses, 1983; Matijevic and Pethica, 1958) or as a stepwise association "reaction" (Tanford, 1974; Mukerjee, 1972; Murray and Hartley, 1935). Although the former description may aid in visualizing certain aspects of micellar solutions, the thermodynamic idea of a phase cannot be used in a rigorous manner. Its requirements of continuity and homogeneity are not met by a collection of micelles in solution and a single micelle cannot be treated as a phase since its properties are sizedependent. The treatment of micelle formation as reaction equilibrium is plagued by the lack of a single stoichiometry. Since a distribution of products is formed (Vold, 1950), one must consider each micelle to be in reaction equilibria with the dispersed monomers. The determination of the many equilibrium constants by experimental methods is impossible. Hall and Pethica (1967) proposed using a smallsystems thermodynamics approach to avoid the difficulties of these two treatments. But their approach cannot be used with ionic systems and is mainly formal, not lending itself to practical use. The thermodynamics of micelle formation remains an interesting problem. The literature is abundant with studies. In addition to the quantities of temperature, pressure, and composition which typically define the thermodynamic state of a typical solution, the thermodynamic behavior of solutions containing surfactant species can depend on the size, shape, and structure of the aggregates which are formed. Thermodynamic properties have been measured and correlated (Burchfield and Woolley, 1984; Woolley and Burchfield, 1984, 1985). Aggregate formation has been investigated from the points of view of classical thermodynamics (Moroi et al., 1984; Muller, 1973) and statistical thermodynamics (Hoeve and Benson, 1957; Owenson and Pratt, 1984). Investigations have focused on size distributions (Ruckenstein and Nagarajan, 1975; BenNaim and Stillinger, 1980), the role of micelle shape (Tanford, 1974; Israelachvili et al., 1976; Ljunggren and Eriksson, 1984, 1986; Eriksson and Ljunggren, 1985; Vold, 1985), and shape transitions (Van de Sande and Persoons, 1985; Ikeda, 1984; Missel et al., 1983; Mukerjee, 1977). While contributing to our understanding of the complex nature of micelle formation, none of these works produced a practical model with predictive capabilities. A semiempir ical model for the thermodynamic properties of surfactant aggregate formation based on molecular thermodynamic processes was developed by Hourani (1984) and was successful at predicting thermodynamic quantities and aggregate size distributions for systems of a single surfactant species in solution. Benedek (1985) developed a different model in the framework of molecular thermodynamics. While it was demonstrated successfully, the use of empirical parameters was more extensive than in Hourani's work. The model of Hourani showed promise of being extendable to multicomponent systems and of being more closely related to other observable molecular phenomena. The beginnings of such an extension are given in this chapter. 2.2 Stoichiometry and Reaction Equilibria for Multicomponent Micelles Since the free energy of a process is independent of the path chosen between the initial and final states, Hourani proposed a process consisting of a series of steps for which the free energy change can be modeled. In this process, the monomers were removed from the solvent and placed in a vacuum at their original densitya gaseous state. Cavities of excluded volume remained in the solvent, to be coalesced in a subsequent step. The monomer gas, considered ideal, was compressed to micellar density, and placed into larger cavities which had been formed in the solvent. Essential to the modeling was the elimination and creation of the solvent cavities. The counterions were handled in the same fashion, with the addition of the necessary electrostatic calcula tions. Extending Hourani's molecular thermodynamic model to solutions containing two or more surfactant species required the addition of new steps and modification of others. The solvent cavity steps were eliminated and the monomers are removed to a liquid state rather than the gaseous. These changes facilitate the handling of multiple components. The development of the multicomponent model is detailed below. Micelle formation in solution yields a distribution of micelle sizes. In addition, a multicomponent surfactant solution has a distribution of compositions (Warr et al., 1983; Scamehorn et al., 1982; Birdi, 1975; Moroi et al., 1974, 1975a, 1975b; Rubingh, 1979; Clint, 1975). To describe this, a set of equilibrium reactions can be written for the formation of J micelles of distinct sizes and/or compositions. For I different surfactant species, Zi, and K counterion species, Bk, aggregating to form J micelles, ZJ, KI NI,Z + N,21Z2+...+ N,,Z+ MllB, + M2lB2 ++A...+ MK ,BK < Z K2 N 2Zl + N22Z2+... + N2Z,+ M12B + M22B2+...+ MK2 BK Z2 (2.1) KJ NjZ ,+ N2jZ 2+...+N jZ ,+M1jB1 +M 2B2+...+ M j B = ZJ where Nij is the number of monomers of species i present in the jth micelle and Mkj is the number of k counterions bound to it. The equilibrium constant Kj for the formation of the micelle is given by K = (2.2) [zilN[z2 ] 21..[ZlJ [B ij M[B21\ 1.[Bkl These are related to the standard state free energy of micellization of the jth micelle by G, =RTlnK, (2.3) The total number of monomers in this micelle is NjA= N (2.4) so that N ,AG J=RTInKj (2.5) where Gj is on a per monomer basis. The mole fraction of monomers in the jth micelle is found by combining equations (2.2) and (2.5). For the dilute concentration range of micellar solutions, ideality of the monomer solution can be assumed, and N exp + xT iln[zJ,])+MkJln[B,]lnCj (2.6)  IxRT i k where Co is the total solution concentration [ ] indicates concentration of species XJ is the monomer mole fraction in jth micelle Nj is the aggregation number of the jth micelle Within the material balance constraint for each species, xi= x, (2.7) equation (2.6) describes the distributions of micelle size and composition in solution when provided with the free energy change as a function of the size and composition of the micelle formed. 2.3 Estimation of Free Energy Changes The standard state free energy of micellization is calculated via the sevenstep process shown in Figure 21. The standard state free energy of formation for a single micelle in solution is the sum of the free energies of the seven steps. Each step is modeled as closely as possible from an observable phenomenon of a similar nature. Certain free energy terms are dependent on the shape of the micelle. The micelle grows as a roughly spherical aggregate until the additional volume of one more monomer would cause the radius of the sphere to exceed the length of the longest alltrans chain. With the addition of the next C )= C1) !4 01 Q)O 4ri ri 4 00 4JH 000 au) ri 0 ri r3 Z 0 4J r.4 0) 4) 0 r1 41 m 4J M a)o i $ C 4 0 'd ] 0 o r '0 4r c Q CO0m 0 0 , P4 a 4) Q p (z4 0 0 (0 < N P > U "0 m Cl C)i Gr U C i 0 0 0 (U (1 w r P p 0 (1) r4 z 0 0 $ 0 0 uM 'i Wn u t z 0 u C W 4J 0 U) 1 ] Q4) ) H U r0 4a ur a I 0 rn it >4 p p 0) 0 *H m r1 0 S4+J o uM 0 tn U2 0 44 U] V0 FO 4 L o p0 m .0 0 4 z w I 4C I 0 1 3 ) 4 Z 4 4) 41M (1) 0) H U a) 0) U U pl 0 pM> p 0 0 P zc pl r ul ro p H 0 Q 0) H >4 S fZ U 0 ] m 4V4 U monomer the micelle must grow with a nonspherical geometry to avoid the formation of a materialfree core. Several geometries have been proposed: spherical dumbells, oblate ellipsoids, prolate ellipsoids, and spherocylindrical rods. While Ljunggren and Eriksson (1984, 1986; Eriksson and Ljunggren, 1985) have proposed that the shape fluctuates between spherical, rodshaped, and even discshaped, Vold (1985) has found little effect of the particular geometric model on the thermodynamics of micelle formation. In this work nonspherical micelles are modeled as prolate spherocylindrical rods. The derivation of micelle dimen sions and surface area based on this geometric model is given in Appendix A. In step 1, a standard state solution of surfactants is transformed into a solution of hydrocarbon chains by reversibly removing the head groups and counterions. Since these are reversibly replaced in steps 6 and 7, the net free energy change for the mere removal and replacement of the head groups and counterions is zero. If no free energy contributions due to the replacement of head groups and counterions are contained in steps 6 and 7, then G I O= 0 (2.8) In step 2 the hydrocarbon solution is separated into I pure hydrocarbon liquids and the pure solvent. This is the reverse process of hydrocarbon solubility, so GRT =,xlnCr (2.9) In step 3 the I pure hydrocarbons are placed into J ideal hydrocarbon mixtures of different compositions and amounts. For this ideal mixing, = J= x,,lnx,, (2.10) In step 4 the J hydrocarbon mixtures are formed into J droplets and placed into the solvent. The free energy of this step is the free energy of forming the hydrocarbonsol vent interface. There are both surface area and curvature contributions to this step. The expression for AG4, based on Buff (1955) and Stillinger (1973), is the surface area of the droplet, S, times the planar interfacial tension, y, of the hydrocarbon mixture, corrected for curvature: ART4 NkTSY (1 (2.1 la) The curvature effect is dependent on the parameter, g. The surface area depends on the size and shape of the micelle. For the spherocylindrical micelle, a second curvature parameter, gc, is used for the cylindrical portion: 2R( I 2R 1 L )](2.11b) RT NkT I XIR)Lk IR) This approximates the cylindrical curvature effect, whose uncertainty has been discussed by Henderson and Rowlinson (1984). In step 5 conformational changes in the hydrocarbon chains are made so that one end of each chain is at the surface of the droplet. The contribution from this step is entirely entropic and may only be estimated. The expression used in this model is AG  S c (2.12) RT ( _R where SC is a parameter of the model. The squared ratio of chain length, 1c, to micelle radius, R, takes into account the very severe conformational restrictions present when the micelle radius is much smaller than the chain length. As indicated in the discussion of step 1, step 6 contains no contribution due to the reattachment of head groups. The quantity Gg6 is the free energy change due to the interaction between the head groups in their positions at the micelle surface. The head groups are modeled as dipoles. For the ionic surfactant species, a charged head group paired with a counter ion forms a strong dipole. Nonionic head groups exhibit weaktomoderate dipole moments. The dipoledipole interactions of adjacent pairs are summed for the free energy contribution of this step. The separation and orientation between two dipoles are dependent on the size and shape of the micelle, with the head groups evenly spaced over the surface of the micelle. The derivation of the head group interactions is carried out in Appendix B. The potential between a pair of adjacent head groups is E P2 ( r )2 E 12 = 3 [ 1 (2.13) where p is the dipole moment r is the pair separation R is the micelle radius D is the dielectric constant of the solvent As indicated in the derivation, this form of the potential takes into account the angle between adjacent dipoles as a function of micelle radius. The total energy contribution from this step is found by summing the contributions from the different pairs in the manner described by equations (B.25) through (B.28) in Appendix B. Step 7 contains no contribution for the replacement of the counterions back into solution. The free energy of this step is due to the difference between the original random distribution of counterions in the micellefree solution and the final PoissonBoltzmann distribution of counterions around the surface of the micelle with a fraction bound in a GuoyChapman electrical double layer. The derivation of Hourani (1984) for the numerical solution of this charge distribution model is applicable here. The entire step is actually the process of discharging the counterions in their original distribution, compressing them into the bound layer and final distribution, and then recharging the counterions. Therefore AG7 contains an entropy contribution from the compression and an enthalpy contribution from the distribution. The bound layer will be populated with dipoles formed by head group/counterion bound pairs and unbound head groups. The charge interactions in the bound layer are included in the dipole and point charge pairings of the head group term, G65. 2.4 The Calculational Technique The free energy of formation for a single mixed micelle with nonionic head groups is obtained by summing the contributions from steps 1 through 6: AG j=(xi, ln Ceq + xi, Inx,,)+ 2nR L 1 +2R,  RT NkT R1) \. R) j SC 2+ 1 (2.14) SR Rja R heads To include the presence of ionic head groups and counterions, the contributions of step 7 must be added to equation (2.14). Such calculations were not carried out here, so this section will pertain only to mixtures of nonionic surfactants. Table 1 summarizes the model's variables, parameters, and required data. Variables Data Parameters Table 1 Arguments of the Thermodynamic Model xil Compositions of micelles Nj Aggregate sizes of micelles T Temperature of system, Kelvin xt Composition of system nc Number of carbons in surfactant chain Ctot Overall system concentration, moles/liter Ceq Solubility concentrations of hydrocarbons, moles/liter yi Interfacial tensions of hydrocarbon mixtures, dynes/cm P Dipole moments of head groups, Debyes D Dielectric constant of water Ic Surfactant chain length, Angstroms Vc Surfactant chain volume, Angstroms3 g Spherical curvature parameter, Angstroms gc Cylindrical curvature parameter, Angstroms Sc Entropy of conformation parameter The values of chain length (Ang.) and chain volume (Ang.3) used in the determination of micelle size and shape are calculated from Tanford's correlations (Tanford, 1972): = 1.265nc + 1.5 (2.15) v, =26.9nc+ 27.4 (2.16) To facilitate the use of computer programs in carrying out the calculations, correlations are used for the required physical data. The aqueous solubility of hydrocarbons used in the calculation of AG20 is obtained from a correlation due to Leinonen et al. (1971): eq = exp[(I 2x q xeq)]} (2.17) The parameter K is fit to the solubility data of McAulliffe (1966), Polak and Lu (1973), and Sutton and Calder (1974). It is found to be a linear function of the hydrocarbon chain length for the nalkanes, but since hydrocarbon solubility in water exhibits a break at decane, two linear relationships for K are used, one for the longer chains and one for the shorter chains. Equation (2.17) is solved iteratively for the hydrocarbon solubility, xeq. For the hydrocarbonwater interfacial tension required in the calculation of AG40, a correlation based on the works of Aveyard et al. (1972), on the surface tension of hydrocarbons, and Jasper (1972), on the surface tension of water, is used: 57.868nc 117.99 (.059nc 1768T (2.18) y=1.381 (2.18) nc+2.4 where the temperature is in Kelvin and the surface tension is in dynes/cm. The value of the dielectric constant of water used in the evaluation of iG6 is given by D =252.422.806329T+.0007469T2 (2.19) which is a polynomial fit of the data of Owen et al. (1961) at atmospheric pressure with temperature in Kelvin. The values of dipole moments used in this calculation are estimated as the dipole moments of molecules of similar structure to the head groups. The three parameters, g, gc, and SC, are fit to the measurable data on the micellar system. These are the mean aggregate size of the micelles and the set of mixture critical micelle concentrations (CMCs) at the system temperature. Equation (2.6) generates an Idimensional surface of the aggregate size distribution of the multicomponent micelles. This is accomplished by first choosing values of the temperature, T, the total concentration, Ctot, and the system composition, xi. Then for each possible micelle composition, xij, equation (2.6) is evaluated at values of Nj from two to infinity (in practice, until Xj/Nj becomes negligible). This generates the distribution of aggregate sizes with micelle composi tion. The parameters must be chosen such that the distribution meets the material balance constraints of equation (2.7). To find the proper values of the parameters, the mean size is calculated as N= ((2.20) ZXJ/N, and the CMC is taken to be the value for which 4(C I+C ) lim 0.5 (2.21) C(01CMC 3Ctot as put forth by Hall and Pethica (1967). Here, C1 is the free monomer concentration and CM is the concentration of micelles, calculated as Ceot Cl CM=, (2.22) Once a set of parameters is determined for a system, the response of the size distribution to changes in any of the input variables can be investigated. The program listing in Appendix C gives the FORTRAN source for the interactive fitting of parameters and calculation of aggregate size distribution on systems of a 22 single nonionic surfactant species in water. This is the simplest application of the model and was used to generate the results presented in the next chapter. CHAPTER 3 RESULTS OF THERMODYNAMIC MODELING 3.1 Parameter Estimation Three parameters of the model, g, gc, and SC, must be found by fitting the mean aggregate size generated by the model to the experimental value at the critical micelle concentration (CMC). The model output is considered to be at the CMC when the total surfactant concentration is equal to the experimental value of the CMC and the derivative of the monomer concentration with respect to the total concentration satisfies equation (2.21), indicating that one out of two surfactant monomers added to the solution at this concentration would join a micelle. The model convergence is quite sensitive to the values of the parameters, so in fitting them to the data, one must either use good initial guesses or approach the values conservatively. Failure to do either of these can result in an aggregate size distribution which has infinite values of Ctot and N, providing no useful information. Each of the parameters has a physical significance which can aid in the choice of the initial guesses. The parameters g and gc are used in the curvature corrections to the planar interfacial tension for the spherical and cylindrical geometries, respectively. Physically, g is the thickness of the spherical interface in Angstroms (Buff, 1955). Although experimental values are not available, it is expected to be positive and small relative to the micelle radius. The parameter gc does not have the same physical connection to the interface (Henderson and Rowlinson, 1984), but by comparison it is expected to be positive and less than g. The parameter SC is the dimensionlesss) conforma tional entropy change for a monomer joining a large micelle. Conformational contributions of the aggregated monomers were studied via a statistical thermodynamic theory by BenShaul and coworkers (BenShaul, Szleifer and Gelbart, 1985; Szleifer, BenShaul and Gelbart, 1985) and values of the conformational entropy were found to be in the range of 8 to 7 for a chain with seven bonds. While these values are based strictly on theoretical considerations with no experimental corroboration, the parameter SC is expected to be of the same sign and magnitude. In approaching the parameter fitting conservatively, the initial values of the parameters are chosen such that the aggregate size distribution will definitely converge. That is, the condition lim must be satisfied. Equation (2.6) defines the aggregate size distribution. In this equation, AjG is the only term influenced by the parameters. In order for the distribution to converge as N becomes large, this free energy change must be greater than the logarithm of the monomer concentration. Overestimating the parameters toward a less negative free energy change will assure that the size distribution converges. The parameters can then be adjusted toward a more negative free energy change as they are fit to the data. The effect of adjusting the parameters can be foreseen by analyzing the corresponding terms. The parameters g and gc affect the behavior of the surface free energy, AG4. For g < R and gc < ic, the usual cases, increasing the parameter decreases the free energy, consistent with equation (2.11). In equation (2.12) it can be seen that increasing the parameter Sc decreases the conformational free energy, AG5. The aggregate size distribution can be divided into two regions. By designating the value of N for which R = lc as Ntrans, the transition of the micelle geometry from spherical to spherocylindrical defines the two regions of the size distribution. For N < Ntrans, the micelle is spherical with radius R. For N > Ntrans, the micelle is spherocylindrical with radius lc and length L. For sizes above Ntrans, the conformational free energy is constant with a value of Sc. At Ntrans the surface free energy makes the transition from being equal to the spherical contribution to becoming dominated by the cylindrical contribution for large N. Figures 31 and 32 show the effect of aggregate size and parameter values on the surface and conformational free energies. A typical aggregate size distribution generated for a nonionic surfactant is shown in Figure 33. A peak occurs at Ntrans, beyond which the fraction of aggregate decreases with increasing N. The total surfactant concentration is proportional to the area under the distribution and is thus influenced by both the height of the peak and the slope of the distribution above Ntrans. Because the distribution below Ntrans rises so rapidly, the mean aggregate size, N, is dependent only on the slope of the distribution above Ntrans The parameter SC, since it contributes over the entire range of N, affects both Ctot and 9. The parameter gc, contributing only above Ntrans, has a significant effect on N but only a slight effect on Ctot, while g influences only Ctot and not N. This causes such interplay between the parameters that a range of parameter sets can produce identical values of Ctot and N. The derivative constraint,  rc' N 0o~ oi o r N o_ o*1 n nr N . 0 O D Uc .0 a Hri c4 a Ir 0 (1 C0 1 4 M 11  0 4 F= r 4) " 0 U,) 1M 14 P4 It 4a J r OZ )  N a m 41 a, U ( 1) o (0 Z ) 44 04 U a)4n 0) 4 ) U ) r r. 4Ju ro0 OC 40 a) D E (1) 00 C > > o 0 +J a ) 0 0 n 4) H o 4x a) Ca) 0C = 4 0 0U 0 r. E 0 < E 4> U)0 a) co a) 4 H 1 H $4 toa po (0 a W 1 0 44 1 1i X 044 0 b M h cw 0 tP b 01X S0 I 4,40 S; 4 )I V 0I () (S N 4J ( 0 0 m) r Q4 10 & 0 c 0 , 4 4 4 E 43 44 C 0U 01 Q) C% C14 0 E Z Q) 1 40 3 I cI c 0) CO r CD I t K) 00 t C E 0  el . Ir4 Q) 40 r 4) U) < I C /~ ~ / )0 MQ) / ~ ~ ( / J=1 4 / / E~iu~ i M< __^'^ .^ ^^ I tM 3 + n + ) (1)u hC    1     [         [   1  [  0^ 1 r f Q 1 i i i i i i i t i 3 (i ( , '4* ^ C T O C) OO P^ O iO *^ r/} ^ r 0C )l g *) r 4 ( ^ *^ klla, 0 C ) 0 >4 l H EI ,J o 0 C0 C0 C 0' H (o ( oo r. 0 0 44 9 Q H 03 Cr z a) gIn H C: ) H 0) N4J 4) H 0 0 4 0 o oo o0 ,r S4I VQ) NU 5 (2.21), can also be satisfied by these parameter sets if they are fit at a value of the free monomer concentration, C1, just below the critical micelle concentration. The method used to fit parameters consisted of first choosing low values of the parameters, resulting in a Ctot essentially equal to C1 and an N insignificantly larger than Ntrans. The parameter gc was increased until N reached the desired value, and then g was adjusted to result in the proper Ctot. This is repeated for different values of Sc to generate the range of parameter sets fitting the data. The experimental data used to generate the parameters is given in Table 2. 3.2 Behavior of the Model for SingleComponent Nonionic Systems Parameter fitting was carried out on systems of different surfactant chain lengths and at different temperatures. For each system, linear relationships were found to exist between each pairing of the three parameters. Corresponding to the manner in which the fitting was accomplished, the curvature parameters g and gc can each be expressed as a straight line plotted against the Sc abscissa. Such a plot is given in Figure 34. These expressions are of the form g= mSc+b (3.2) gc=nSc+ d (3.3) Table 2 Data Used in Fitting Model Parameters for Aqueous Solutions of Nonionic Surfactants Surfactant T (Kelvin) CMC (M) Mean Size C0lH21(OC2H4)60H 298.15 9.0E04 a 73 a 308.15 6.6E04 260 318.15 6.4E04 640 C12H25(OC2H4)60H 298.15 8.7E05 b 400 c C14H29(OC2H4)60H 298.15 1.1E05 a 3100 a a: Balmbra et al. (1964) b: Corkill et al. (1961) c: Balmbra et al. (1962) CN r, L (Q) ED rj O_ 0 (9 CL E ri 0 0 ci cO l to O 4' r) Nc 0  JaLGWDJDCd aJnfDAJnl. u U 4C 0 a) 4 n 0 0 C 0) r! e w u z SC~ 04 r4. 0 4i CO U) cco rC iO to II C) ''r4 U 0 fflw 13 a(N U a) C C) H 4. Co tc) r3 i C0 ) rl C) SiH Ei r a)a) 0 0i cl, 0 0 a) n :s i (a1 i *H C r p I 4 0 0 (C 11 o 0) The values of the slopes and intercepts found for the systems modeled are given in Table 3. Many aspects of surfactant behavior have been found to be linearly dependent on chain length and/or temperature (Rosen, 1978). The information in Table 3 indicates that such relationships are also possible for the interaction of the parameters of this model. Though the three chain lengths and three temperatures investigated cannot conclu sively indicate linear behavior, they can indicate whether this behavior is likely. Linear regression of the slopes and intercepts of Table 3 versus chain length and temperature resulted in the following equations and correlation coefficients: m=.0089T+.8817 b= .0549T20.20 n=.0137T+ 1.421 d =.1 125T45.09 for the 10 carbon nonionic, and m=.1464n, .3170 b= 1.418n, 18.10 n=.2201n,.4676 d= 1.553n, 27.29 at a temperature of 298.15 Kelvin. R=.9998 R= .9778 R = .9999 R=.9810 R =.9999 R= .9988 R =1.000 R =.9991 (3.4) (3.5) (3.6) (3.7) (3.8) (3.9) (3.10) (3.11) Table 3 Slopes and Intercepts of Curvature Parameter Dependence on Conformational Parameter T (Kelvin) g slope g int. gc slope gc int. 318.15 1.96163 2.80775 2.94330 9.43866 308.15 1.86885 3.15286 2.80330 10.1776 298.15 1.78296 3.90518 2.66898 11.6881 298.15 2.07048 1.32359 3.10874 8.80775 298.15 2.36862 1.74162 3.54946 5.47762 Substituting (3.4) through (3.11) into (3.2) and (3.3), the temperature relations for C10 nonionic are g =(.0549.0089Sc)T 20.20 +.8817Sc (3.12) gc=(. 1125 .01372c)T45.09+ 1 .4214Sc (3.13) and the chain length relations at 298.15 Kelvin are g=(1.418.1464S,)n, 18.10.3170Sc (3.14) gc =( 1.553 .2201S,)nc27.29.4676Sc (3.15) These relations are only valid for values of SC which, for each chain length, produce values of the curvature parameters which are neither negative nor larger than the micelle radius. In this sense, SC is dependent on chain length, but there remains a lack of uniqueness in the parameter fitting for the single component case. This results from the finding that the CMC and the concentration derivative (2.21) are not independent of each other. The free energy of micellization and its various contributions, as given by the model, are shown in Figures 35 through 38. These plots were made for two chain lengths at the same temperature, two temperatures with the same chain length, and two parameter sets at the same chain length and temperature. Figures 35 and 38 show the effect of two different parameter sets for the same system. While the surface and conformational contributions are different due to the different parameters, the total free energy of micellization is the same in the two figures. This follows c o C N C) C CN o ) '1) o 0 o I o 040 (1) a U4 in UC> C0 faQ o > 0) .4^ U C 3 0 ) 44; i r.? Si~ U 1 C)( 4c) 4 0 0c U 0 *C 03 O J S 4Jo 0 i ol 0 (a 0S.u 0 CL a) 4 0 0) 0 1 0 4 05 CC Q) 44S >(0 04C : 41 1 (Z 1 Q) 4C 4 P 0r H 43 U) (0 4C (a  rcu nU a la)a 0 *P CP ) ^r * 4 i fE (0 3 c a0 So f (l E=  u L I CO I 00 II a  I I ( 0l I 0 I o ,:0 i O 1. LI C I u M O 0 I 0o N 0 0 U 1*o f,4 o 1 >I < 0 I o 1.4 (0 C 0 0 *Q)4 ao) to r= 4 ; l.p 0 S0 i 41 I( ) 44 4 0 O O O U 0 (, II r I i t o r P~ac 0 la) I 0 * SU I "I 0 I i o3 4 C) i' i 4 r. 0 4 C 1) 4) C if 0o O 0 ( i i I IN 4C4 S _0 0 4 0 /r / CO ) S/0 I I I I I o a 0 n to 0 LO 0 to 0 : 1 a / Q, 4) UlC U!' 0 C 0 CO N 0'3 C cD CM 0 LO 0 'O 0 L" 0 to (1 "  r " ,,NI I (1) 0 0 Uo I)  1rr1 CE0 0 C m U r4 0 ca, U 4c) 04 cra) 0 *P ( 0 0 4t 0( Qp a, a P c co) ( r( C0 U E cd from the fact that the parameter sets were fit to the same concentration and mean size data and generate identical size distributions. 3.3 Aggregate Size Distribution and Concentration Behavior With the parameter sets determined by the methods described above, the aggregate size distributions and the effects of surfactant concentration were investigated for the nonionic alkyl surfactants. The size distributions for this class of surfactants exhibit a rapid increase in the spherical region of N, peaking at Ntrans In the nonspherical region, the distribution decreases very slowly through large values of N, producing the large mean aggregate sizes which have been measured for the nonionic surfactants. The size derivative of the distribution in the nonspherical region is negative. With increasing free surfactant monomer concentration, it becomes less negative, resulting in a higher mean aggregate size and total surfactant concentration. The effect of surfactant concen tration on the modeled size distribution and mean aggregate size of the C12 nonionic surfactant is shown in Figures 39 and 310. Below the CMC, the model calculates a mean aggregate size approximately equal to Ntrans, present as the initial plateau in figure 310. This value is meaningless, however, since the concentration of micelles in this range 41 o co CV 4 IZ 4 Z: I II I I/ / C 1J /43 c cr 01 rr4 _0 (N 0 o .J N Co 0 0) cc o o O*c 4 C ( C 4C 0 0 //08 N 0 r c00 04 I I n IIII I  44 ) Co r C 0 40 C ) H / tN r i T T ,I N)501 Q) / / N) 42 r 4 O E1 \ . 0 ,. o r C) \j 0 i OO 0 \LJED o L oC 0 \ 0 L 4 \C C) \ . I 0 m (NO \ CM ( O CO U') (C( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 U 0s w w "0 w u w w .   o~~~~~a o ooooo ^r CM~~~O O Q t NC O O O 'C 0 0 3 t( 11( r^ n cM M CN CM M ' ' *r 'r QZIS9^D~j66vUDQ^ of surfactant concentration is insignificant. The result is merely the ratio of two extremely small numbers, since the model always produces a finite value for the concentration of micelles (equation (2.6)). Measurement of properties such as surface tension, which depend on the free monomer concentration, show that beyond the CMC the free monomer concentration is virtually unaffected by an increase in the total surfactant concentration. The relationship between the free and total surfactant concentrations, as predicted by the model, is exhibited in Figure 311. At total concentrations lower than the CMC, the majority of any surfactant added to the solution exists as free monomers, while at concentrations higher than the CMC the bulk of any additional surfactant micellizes. The definition of the CMC used in this work, as expressed in equation (2.21), quantifies the break between these two behaviors of the solution. Figure 312 shows the concentration derivative of equation (2.21) for the data of Figure 311 with the CMC identified. 3.4 Chain Length and Temperature Effects It is known that the hydrophobic nature of the alkyl chain increases with its chain length, as evidenced by the decreasing solubility in water of the increasingly long 44 I 4 LJ 0 0> E 0  c. r 4 j 0 o ro I (U 0 Li D '7 o L. U C o i U 0 03 OZ 0 o 0 U iC 40 * O L J ir 0 r 0 m i < ,i0 O n o oi 01 0 0 4 PH 40 4 C o Co o "I Iuu I III I I Li L O O Li N LI UL LO Li 'J Lr ) J UJ (0 UJ Ud UJ UrJ *r1 0S 04 03 4 0 (0 04 0 0 444 F. *J 'UOIpDJjUGOUOQ JGWOUO0j GGJj 45 cG UJ 0 0, 0 I N ~ L. 0) 0 C .C 0 0 (a C a) W o 4 4 u 3 O r N 4 1 0 0 C C0 4; (o L. o o 4J U) a) n ^ 4J 4> 4 > E 43 4 0 3 0) rO ri 1 vC o I 0 H 0 .rH 0 I (a0) > ((o + o c____ "t n 06 M O 0 00 l 0 N > (11 + +e chains (McAulliffe, 1966). With alkyl surfactants, this increasing hydrophobicity promotes micellization. While no significant aggregate formation is observed for chains of six carbons or less, increasing chain length brings about the formation of larger micelles at lower critical micelle concentrations (Table 2). The elements of the model which depend on chain length are the hydrophobic free energy, taken directly from the aqueous solubility, the curvature parameters, accounting for the difference in micelle radius, and the conformational free energy. Equations (3.14) and (3.15) show the relationship of the parameters to the chain length. The aqueous solubility of the alkane, which is chainlength dependent, is used by the model in the fitting of the parameters. The dependence of the aggregate size distribu tion and free energy change on the surfactant chain length is shown in Figures 313 and 314. The size distribution broadens with increasing chain length, yielding a much larger mean aggregate size with a less significant increase in the CMC. The free energy of micellization becomes more negative with increasing chain length, showing that micelle formation by surfactants with more hydrophobic area is more thermodynamically favored. 47 C c'4 o o 04 N ra 'I NW Q44J rI 0o 0 C E / C) I0 0 c: c o (0 C 0 > r4 L   o o Pr4 ro N4 I i r S1 4JI 1 iI I I I I (,),4 0", 1 1 1 1 1 1 1.^ ^ y co rc LJ 44 1 2 .) u U I I " 0 I I " II  4.4 r  ) 0 Q)2) c 0 a) 0u O C Q4 A C0 N Q . 4 ~ c4 i1 I 0 U < j 0 N .I1 4 ONC S i co C J 1 0 ",' II I"  rl  / ___/ nM Just as the aqueous solubility of an alkane decreases with increasing temperature, the CMC of the corresponding surfactant decreases and the mean aggregate size increases. Again, by incorporating the solubility data into the model, the proper temperature dependence is given to the parameters. The temperature dependence of the parameters is in the form of equations (3.12) and (3.13). The dependence of the aggregate size distribution and free energy change on temperature is shown in Figures 315 and 316. The effect of increased temperature is similar to that of increasing the surfactant chain length: a broadened size distribution and a more negative free energy of micellization. But while the first increase of ten degrees Kelvin has a much more pronounced effect on the free energy change (and hence the size distribution) than the next increase of ten degrees, the mean aggregate size is more affected by the latter. O O O C C 0 Ca co o 0 CD CM  "' _/ 0 C3cn o'NX) m c0 S4 0 03 4 30 C 4 U 0 J .c 3 E4 cO a) m ,. '0. Q)r c1 g 3I =) 0c Fp ic L 0 C C c C C c0 C C oa o 0 ,C 0 L0 1 N 0(D k_ en to lO ( 0 ) rN T c4 n rD D (D ( I rp r,. r, 4 cl OrJ is 4, C) a) Q 4Q E) C) 41O U4 C .) .U. Cl ) C 1X1 CQ CHAPTER 4 MOLECULAR DYNAMICS SIMULATION OF SURFACTANT MICELLES 4.1 Background The development of a predictive model of micelle behavior, such as that of Chapter 2, requires a knowledge of the structure of the aggregate in order to describe its thermodynamic properties. The calculation of a surface free energy contribution necessitates certain assumptions about micelle size, shape, and surface roughness. The solvent interface must be modeled and solvent penetration must be estimated. Calculation of the conformational contribution to the free energy calls for details of the molecular conformations in the free and aggregated states. The head group contribution depends on the average positions of the head groups. While experimental results are the most desirable basis for thermodynamic models, measurements of the structure of micelles are limited in the information they can provide. Therefore, current models of micelle behavior, including that derived in Chapter 2, are based primarily on theoretical considerations. The most informative techniques for studying micelles in solution are spectroscopic. Light scattering, NMR, and smallangle neutron scattering (SANS) have been used extensively to study micelles. However, results of these experiments often conflict, as in the question of whether there is solvent penetration into the micelle core, and their description of the micelle surface and internal structure is incomplete (Tabony, 1984). Quasielastic light scattering has been used effectively to determine micelle size and shape (Corti and Degiorgio, 1981a, 1981b) and Raman light scattering has been used to investigate chain conformations (Kalyanasundaram and Thomas, 1976). The absence of solvent from the micelle core was determined with NMR techniques (Mitra et al., 1984), though the opposite had been reported earlier (Menger, 1979). The most promising experimental means available of studying micelle structure is SANS. It has been used to study the micelle surface (Hayter and Zemb, 1982), solvent penetration (Tabony, 1984; Cabane et al., 1985), and chain conformations in the micelle interior (Bendedouch et al., 1983a, 1983b). The spectra generated by the above methods must be carefully interpreted to yield correct information about the micelle surface and interior. It is in this interpretation of the results, particularly in accounting for intermicellar effects, that inconsistencies arise (Cabane et al., 1985). Computer simulation is an alternative means to answer questions raised by the experimental measurements. For a given model of surfactant and micelle behavior, simulation can provide exact quantitative descriptions of micelle structure. Simulations of micelles have been conducted by both Monte Carlo (Haan and Pratt, 1981a, 1981b; Owenson and Pratt, 1984) and molecular dynamics (Haile and O'Connell, 1984; Woods et al., 1986; Jonsson et al., 1986) methods. Thus far, results of the simulation studies have been promising. Their contribution to a better understanding of micelle structure has been limited only by the uncertainties about the model micelle. Monte Carlo calculations are carried out on a lattice of sites. Movement of a molecular segment from one site to another is allowed or disallowed based on the energy change associated with the move. The restriction of segment positions to the lattice sites limits the possible conformations available to the molecules. To allow more possible conformations, the number of lattice points, and hence the computing effort, must be increased. This restricts how far such a simulation can go in accurately representing reality; in addition, the Monte Carlo method precludes the opportunity to study the dynamics of a system by simulation. The molecular dynamics method provides a continuum of sites for molecular positions and is limited only by the accuracy of the forces built into the model. In the present work, molecular dynamics models were developed based on the work of Haile and coworkers (Haile and O'Connell, 1984; Woods et al., 1986) to investigate the effects of surfactant chain length, head group mass, and head group size. 4.2 The Molecular Dynamics Method The molecular dynamics approach to computer simulation consists of summing the forces on each particle in the system and solving Newton's equations of motion for the resultant positions and moment (Haile, 1980). The problem of constructing a realistic model system then becomes one of supplying appropriate models for the intermolecular and intramolecular potentials. It is assumed that a given particle's interactions with its neighbors are pairwise additive. Thus, the potential between two particles at time t is of the form U (r,,(t))= U( ,(t),,c(t)) (4.1) The force between the two particles is given by F(r i(t)) = VU(r, (t)) (4.2) and the total force on the particle of interest is the sum of the forces resulting from interactions with all other particles: () = (r (t)) (4.3) j The resulting acceleration and velocity of the particle are obtained from Newton's Second Law of Motion, relating force, mass, and acceleration: F,(t)= m (4.4) dt2 The set of secondorder differential equations for all of the particles in the system is solved at each time step of the simulation. The means employed in this work for determining the new positions of all of the particles from the forces at the previous time step is a fifthorder predictorcorrector algorithm (Gear, 1971). This algorithm consists of three procedures which are repeated at each time step. The position and its first five time derivatives at time (t+At) are predicted for all of the particles by Taylor's series expansion of their values at time t: dr1(t) d2ri(t)(At)2 d3r1(t)(At)3 ri(t + Lt)= r,(t)+ A t + dj + d t dt dt2 2! dt3 3! d4r(t)(Azt)4 d5r(t)(t)5 (4 + + (4.5) dt4 4! dt5 5! d d dr(t d )ri(t) dari(t)(t)2 r ,(t+dt)= + (t) + dt l dt dt2 dt3 2! d4ri(t)(ot)3 dsri(t)(At)4 + + (4.6) dt4 3! dt5 4! 2r (t+ t)= + ) + + (4.7) dt2 dt2 dt3 dt4 2! dt5 3! d3 d3r1(t) d4rr(t) dri(t)(t)2 (4.8) dtr.(t+ tL) = + )+* v ; (4.8) dt3 dt3 dt4 dt5 2! d4 d4r (t) d5ri(t) dt4 dt4 dt5 d5 d5r1(t) Ar(t+dt)= (4.10) dts dts The forces Fi(t+At) are then evaluated at the predicted positions. Finally, the predicted values of the position and its derivatives are corrected according to the error between the acceleration predicted by the series expansion and the acceleration calculated from Fi(t+At). These calculations are carried out by the subroutines PREDCT, EVAL, and CORR in Appendix D. The simulations described in this work were carried out at constant temperature. The temperature of the system is related to the velocities of the particles by T=3kZmj4 (4.11) The velocities of all particles are scaled at each time step so that the temperature is maintained at a constant value. Subroutine EQBRAT in Appendix D accomplishes this task. To start a simulation, values of the positions are required to evaluate the initial forces. Initial values of the five derivatives of position are required by the fifth order predictorcorrector algorithm. In these simulations, initial positions were assigned according to a lattice of sites and random velocities were assigned the particles (subroutines INTPOS and INTVEL). The moment were scaled so that there was no net linear momentum of the system. The initial accelerations are calculated by equation (4.4) and the third and higher derivatives of position are assigned initial values of zero, since there is no way of evaluating them. The algorithm recovers rapidly from this initial condition, arriving at proper values of all derivatives within 1020 time steps. Since intermolecular potentials are significant over a pair separation which is usually small relative to the dimensions of the molecular system, computing effort can be reduced by employing truncated potentials (Haile, 1980) and maintaining a neighbor list (Verlet, 1967). In a simulation where the majority of possible pair separations are so large that the corresponding pair potentials will be very small, it is more efficient not to calculate the potential for pair separations which are beyond a certain value and do not contribute significantly to the system properties. Instead, these contributions are neglected during the simulation and a correction is made to the system properties. A cutoff distance, rc, is set such that the potential can be neglected for rij > rc. The intermolecular forces are truncated at rc and shifted vertically, going to zero smoothly at rc (Nicolas et al., 1979). The truncated potential is obtained from this shifted force by equation (4.2). This eliminates a step change in the potential and force at rij = rc, which could have a disruptive effect on the energy conservation of the simulated system. The general equations of a shifted force, Fs(r), and its potential, Us(r), are given below: dU3(r) [d U(ry Fd,(r) + d dr r r, r U,(r) = U(r) U(r,) (r r By checking the value of rij and bypassing the calculation of F(rij) for rij > rc, a significant savings in computing effort can be realized. Further savings can be made by not checking pairs which have a low probability of rij < rc. This is accomplished through the use of a neighbor list. A particle (i) in the system is surrounded by a sphere of radius rc. Other particles (j) inside this sphere will interact with it, since rij < rc. From one time step to the next, particles will enter and exit this sphere. Over some short period of time, all of the particles which have interacted with the central particle will lie inside a larger sphere whose radius is rlist. If, over this brief period of time, only those particles inside the larger sphere were checked, all of the particles interacting with the central particle would be found without checking all of the other particles in the system. Through a neighbor list a record is maintained of all the "neighbors" within a distance rlist of each particle in the system. Only these neighbors are checked for rij < rc. The neighbor list is updated periodically based on the value of rlist and the dynamics of the system. The molecular dynamics method represents a massive computing effort for any system large enough to be of interest. However, using the techniques described above to improve efficiency, and with the development of supercomput ers and parallel processors, it has become a practical tool for quantitatively studying molecular behavior. 4.3 The Model Surfactant Molecule The surfactant molecule studied in this work has a linear alkyl chain of eight carbon groups. All of the groups are considered to be methylenes. A polar head group is attached to one chain end; head groups of different sizes and masses were used in the simulations. The bond lengths and angles along the chain are those of a normal alkane. These are given in Figure 41. Three types of intramolecular interaction take place in the model surfactant molecule: 1. Bond vibration 2. Bond bending 3. Bond rotation 112.150 1.539 A Head Group Figure 41. The model surfactant molecule, consisting of an eightmember alkane chain attached to a polar head group, shown in the alltrans conformation. Bond angles and lengths of the alkane chain are indicated. The bond to the head group is head group dependent in the different simulations conducted. The bond vibration potential used is that of Weber (1978), taken from a simulation of nbutane. It is a harmonic potential about the equilibrium bond length. The potential and the forces on the two groups are given by U, b. Iy. b.b, (4.14) F: b, =yo b,bo (4.b1 ., bo =y. bb , (41) where bi is the bond length between groups i and i+l, bo is the equilibrium bond length, and y. is the force constant. The bond bending potential is also taken from the simulation of Weber (1978). It is a harmonic potential about the cosine of the equilibrium bond angle with forces on the three groups forming the bond angle. For the angle  formed by groups i, i+1l, and i+2, UbG.=21 yb.cos0 cos,12 (4.17) 1 [b,S 1 b, \ 1e, =yb cosO,cosQ b, b(6 cos (4.18) F ,.yb cosO6cosB ( :: b 1 (bb I b ~ +l [b b S(4. 1 .) ib2(Oi)=Yb(COS oCOS ) b,+ 1 bi l b ,ybtcocos Jb1k+b b,+ cose1, (4.20) The model surfactant molecule is comprised of eight bonds connecting its nine groups. Rotation about any of the inner six bonds results in a change in the conformation of the molecule and a change in its potential energy. A set of four groups in sequence along the chain molecule contains a dihedral angle formed by the three bonds connecting the groups. This angle is a measure of rotation about the middle bond, and the potential energy is a function of the dihedral angle. The bond rotation potential chosen for these simulations is that of Ryckaert and Bellemans (1975): U(O)=y,(1.1161.462cos 1.578 cos2+ 0.368 cos3 +3. 156cos 4 +3.788cos 5 ) (4.21) where 0t is the dihedral angle in radians. The bond rotation potential is illustrated in Figure 42. The force on each of the four groups resulting from this potential and equation (4.2) has been derived by Woods (1985) and the resulting calculations are carried out by subroutine RYFOR in Appendix D. The intramolecular potentials account for all of the interactions among adjacent groups on a molecule as they are taken two (vibration), three (bending), or four (rotation) at a time. Two groups on the same molecule which are separated by more than three intervening groups do not o C) r, Q) CC Q}. 4' H 4 0 4 VC) 43 / 'I C C < o ,. ,4 / 0 0 4 0 4J (a 0 0 fu Q m*4 m > ON \ ov0 0) a) a: ) S 4J J  4z U 0 TS a) C)/\ () F r  a, 0 4J >0 I S I 00 WI O IC" 0 ,4" 0 C 0 ( 0 M \0 E ) rC 0 4) r \  1 C o P 0 E rmrlP r. (3 0 014) p O (a 4 0 3 I C *"4 n .t .e) 14 ia ) (044 ) P O L i C C1 . *H*H P Hr >1 1 :zP > FX4rOM r interact through any of these three potentials, but may experience the intermolecular potential if they become sufficiently close to one another. 4.4 The Model Micelle The micelle used in the simulations contains 24 surfactant molecules. The individual groups of the molecules interact through intermolecular potentials, and the micelle is surrounded by a spherical shell which models the solvophilic effect on the head groups and the solvophobic effect on the chains. The intermolecular interactions which take place in the micelle are 1. Head groupshell 2. Alkane groupshell 3. Head grouphead group 4. Alkane groupalkane group 5. Head groupalkane group The shell does not attempt to explicitly model the solvent molecules. To do so greatly increases the computing requirements of the simulation and a previous attempt at such explicit modeling (Jonsson et al., 1986) did not show it to be a clear improvement over the present poten tialfield model. Rather, the solvophilic and solvophobic natures of the head group and alkane group are modeled in a simple manner intended to keep the head groups close to the surface of the micelle and prevent the chains from extending out from the micelle. This simulates the minimum chainsolvent contact of accepted experimental results (Tabony, 1984). The head group interaction with the shell is through a harmonic potential about an equilibrium radial position: U(r,) = Yh(riro)2 (4.22) F(ri)=2yh(rir)r (4.23) where ri is the position of the head group and ro is the equilibrium position. The force constant is y,. The alkane group interaction with the shell is through a repulsive potential: U(r,) = c ) 12 (4.24) T(r,) )= 2c(r.)12 (4.25) where riw is the quantity rwri, rw being the radial position of the spherical repulsive shell, and rm and E are the radius and energy of minimum potential for the alkanealkane intermolecular potential. The interaction between two head groups is described by a potential comprised of both dipole (r3) and softsphere (r12) repulsions: r= i r + r h)12(4.26) \ 'i/ / r'i/ This potential is designed to spread the head groups apart on the exterior of the micelle. The force between a pair of head groups is obtained by applying equation (4.2) to equation (4.26). Fr + )1 ]+ (4.27) Applying equations (4.12) and (4.13) to use the shift edforce potential for rij rhh 3 hh 4hr.h 3 rrhh )r S+  13 ruj rij ( r, ) r, + 3r hr4y +12r 1Y3J (4.28) rhh ,\ r r, F(rij) = E 3 + + 12 (4.29) r i. r, rI I J The intermolecular potential between two alkane groups is a (6,9) form of the Mie (m,n) potential (Reed and Gubbins, 1973): U(r) = F[2( 3 (4.30) (rFj) =18E r,,, (4.31) ) r _) r) r _) 6 \r r, J =rj 0r,j 2 rJ 18 7 (4.32) rmq [ r, r, S9 6 F 7( f'f '\1 = 18E r Fm r] (4.33) r10 r, rL, r, rj (3 This potential is also used for the interaction between a head group and an alkane group, except that the radius of minimum potential, rm, is adjusted to account for the difference between the diameter of the head group and that of the alkane group: headalkane 14. m 2(rm +rhh) (4.34) Figure 43 illustrates the intermolecular potentials used in the model micelle, and Table 4 gives the values of the parameters used in the intramolecular and intermolecular potentials. The model micelle is assembled by initially placing the 24 molecules with their head groups evenly spaced over the surface of a sphere of twice the expected micelle radius. Each chain is directed radially inward, in the alltrans conformation. With the bond rotation force constant reduced to onetenth of its normal value to facilitate chain packing during startup, the simulation is begun. The radius of the confining sphere is reduced every ten time steps until the appropriate radius is achieved. This is the value that would give the density of the analogous liquid alkane, adjusted to a pressure which fluctuates about zero with an amplitude of less than ten atmospheres. After this, the 2,:93  1800 I 1600 1400 " 1200 (6,9) 00 (12) \ 1000 ............ (3,12) \ U(r) 800 \ J/mol 600 \ " 400 200 \ ..... o " . ... .. .. .. .. .. .. .. . 200 400 600 ,i  0 2 4 6 8 Pair Separation, A Figure 43. The intermolecular potentials for segmentseg ment (6,9), segmentshell (12), and headhead (3,12). A value of 419 Joule/mole was used for c; both rm and rhh were assigned values of 4 Angstroms. Table 4 Parameter Values for Potentials Used in Potential Parameter Value Bond vibration + Bond bending + Bond rotation + Headshell Segmentshell Headhead Segmentsegment Headsegment YV bo Yb 0o Yr Yh lo E rm E rfhh rm ha r 9.25x105 1.539 1.3x105 112.15 8313 785 4 419 4 419 419 4 419 (rm+rhh) 2 Simulations Units Joule/A2/mole Angstrom Joule/mole degree Joule/mole Joule/A2/mole Angstrom Joule/mole Angstrom Joule/mole Angstrom Joule/mole Angstrom Joule/mole Angstrom a Weber (1978) b Ryckaert and Bellemans (1975) c Woods (1985) * Parameter varies with head group size. + See also Table 6. rotational force constant is restored to its desired value and the system is allowed to equilibrate. When there is no net drift in the energy of the system or in the fraction of trans bonds, equilibrium is considered to have been reached. The simulation is continued, periodically saving the positions and velocities of all of the groups for subsequent analysis. The time step used in simulation must be small enough to track the motions of the molecules and maintain stability, yet large enough to avoid unnecessary computing effort. The ideal time step can only be found by trial, which was conducted on a simulated liquid alkane system by Weber (1978). That result, At = .002 psec., is the basis for the time steps used in the present work. 4.5 Summary of Computer Simulations The following simulations were carried out on the model micelle: 1. Head group mass and size same as chain segment. 2. Head group mass greater than segment, size the same. 3. Head group mass and size greater than segment. In addition, a fourth simulation was carried out on a hydrocarbon droplet of 24 ninesegment molecules by setting the head group to the size and mass of a methylene and replacing headhead and headshell interactions with segment interactions. From these simulations the independent effects of head group size, mass, and solvophilic nature can be observed, and by comparison with comparable simulations of dodecyl surfactant micelles (Woods et al., 1986) the effect of chain length can be evaluated. A brief summary of the simulations is given in Table 5. To carry out the first two simulations, 24 model surfactant molecules were placed with their head groups on a sphere of radius 24 Angstroms and their chains directed radially inward in the alltrans configuration. The groups were given initial velocities, the rotational force constant was reduced to onetenth of its normal value, and the simulation was begun. Every ten time steps, the radius of the confining sphere was decreased by .01 Angstrom. This scaling down of the model micelle was continued until the radius of the spherical shell reached a value of 12.00 Angstroms, slightly higher than the radius of an analogous hydrocarbon droplet. Tenthousand time steps were run at this radius to allow the micelle to recover from the scaling down procedure and observe the pressure at the shell. The pressure fluctuated about zero with an amplitude of approximately two atmospheres and analysis of the positions of the groups revealed no effect of the initial conditions. This model micelle was chosen as the starting point for simulation 1. ) I k.0 col co co U > 0I a E4C VQ 1 r4 H c4) 5 I 0ay U I I 1 * S CO co C Co C 4 tn 4) o H ( to H W u 4J Q i rl (1S Z U) 0 'U U 8 i (1 ) o oi 0 00 (a o m U2 U) 41 C J AI c ~ ' o (S ( ) oQQ u u i] >1 14Cl to cp ri (a 14 C.CW a ( 00 >n c m Q 4 o o co (a ^ n ^O *~lQQ o j\ _i ^ ^ '1 ' (C ^i fa m 0l~ ,I y U o m U) ^ (n (o CO U)U3 '0 U) S 0 040 :; I4 (0 Q LO 3 aj c 0 ou ^ n *c Q) H N P4  H0 P>1 U) 4' 0) 6i < rl g *rt .H0 o (a 0 o &i c~ *Q t M2 CO 0a C(O Q) Simulation 2 was started from this same model micelle, after increasing the mass of all head groups by a factor of seven and allowing 10,000 time steps for equilibration. This mass was chosen to represent a common ionic head group, the sulfate ion. Only the mass of the head group was changed; all of the intermolecular and intramolecular interactions remained the same as in the previous simulation. The expected result of this would be simply a change in the dynamics of the head group. The program used to compute simulations 1 and 2 is listed in Appendix D. In simulation 3, the change in the head group was more extensive. An attempt was made to represent the sulfate head group in all of its intermolecular and intramolecular interactions. The mass of the group remained at seven times that of a methylene, while the diameter was increased to 2.45 times that of a methylene. The equilibrium bond length and bond angle for the head group were changed, as were the force constants for bond vibration, bondangle bending, and bond rotation involving the head group. These values are given in Table 6 and compared with the methylene values. The derivation of the head group intramolecular interaction parameters based on the work of Muller and coworkers (Muller and Nagarajan, 1967; Muller et al., 1968), Cahill et al. (1968), and Blukis et al. (1963) is given in Appendix E. 75 Table 6 Bond Parameters of "Sulfate" and "Methylene" Groups Parameter "Methylene" Value "Sulfate" Value Units 9.25x105 1.539 1.3x105 112.15 8313 2.7x104 2.6 9.1x105 140 20000 Joule/A2/mole Angstrom Joule/mole degree Joule/mole Due to the increase in the size of the head group, the range of head group interactions was significant in comparison to the size of the micelle and truncation of the potentials could no longer be justified. Therefore, neighborlisting and truncation of intermolecular potentials were not employed in this simulation. It was also found to be necessary to decrease the time step to maintain stability in the energetic of the system. It is for these two reasons that the computing efficiency (simulation time per CPU time) of simulation 3 is dramatically lower than the other simulations (Table 5). The program used to compute simulation 3 is listed in Appendix F. Simulation 4, the hydrocarbon droplet, was conducted in the same manner as simulation 1, except that the head group was replaced by a chain segment (all segments have the same properties to simplify computation) in its interactions with the shell and other groups. The radius was scaled down to 11.928 Angstroms to achieve the liquid hydrocarbon density of .7176 g/cc. Since there was no sorting out of the head groups for separate interactions, the efficiency of this simulation was 18 percent greater than that of the first two simulations. The simulation computations were conducted on a Control Data Corporation (CDC) Cyber 205 computer at the Supercomputer Research Institute (SCRI) in Tallahassee, Florida. Access to the SCRI facilities was provided through a grant from the United States Department of Energy. All of the simulation programs were modified for the CDC FORTRAN 200 Vectoroptimizing compiler. By employing parallel processing techniques when possible, this compiler provides very efficient operation from code which is highly compatible with ANSI standard FORTRAN. One important difference lies in the default word length of the Cyber 205. Real numbers on this computer are eight bytes in length, compared to four bytes on most other computers. To achieve eightbyte precision in standard FORTRAN, the REAL*8 variable type is used. The programs in Appendices D and F achieve this level of precision on the Cyber 205 with default variabletype coding. Further detail about the simulations, along with the results of their analysis, is given in the next chapter. CHAPTER 5 RESULTS OF MOLECULAR DYNAMICS SIMULATION In this chapter the results of structural and dynamic analyses performed on the output of the simulations are reported. Certain of these results are reported as timeaveraged properties; others are given as the change in a property with time. The former are averaged over the entire length of the equilibrium run (Table 5); the latter appear in the figures over a common time period for ease of comparison. The simulations are referenced by number, Run 1 being the case with head groups of the size and mass of a chain segment (methylene); Run 2 the case with head groups of the same size as a chain segment, but a mass seven times greater; Run 3 the case with head groups of 2.45 times the size and seven times the mass of a chain segment, and a revised head group intramolecular potential; and in Run 4 the head groups replaced by chain segments to model a hydrocarbon droplet. 5.1 Mean Radial Positions of Groups The measurement of mean radial positions of the groups on the surfactant chains is one of the simplest, yet most telling, evaluations of micelle structure that can be made on simulation output. While this measurement cannot be obtained reliably for all groups by experiment, the position of groups is a fundamental attribute of structural models of micelles (Gruen, 1981). The mean radial position for group i (i=l for head, i=9 for tail) is given by R,= J The quantity Ni(r) is the number of occurrences of the center of a group i at a radial position (relative to the center of mass of the micelle) r at time t and the angle brackets denote the average over the time period of the simulation. The unsubscripted N is the number of molecules in the system, equal to the sum of Ni(r) over all values of r and i. The standard deviation corresponding to this mean is given by i 4 [ The mean radial positions of the centers of the nine groups bracketed by one standard deviation in each radial direction are plotted for the four simulations in Figures 51 through 54. Assuming a normal distribution of a group's position over time, the range shown for each group in the figures includes 68 percent of its positions, while twice this range includes 95 percent.  i IIIr INI I 1 I I I 1 1 1 I 1 1 I i I I i I 1 I CN   'V 'sn!poDj 0 0a Lo CD C re 4 0 0 O) Cro 1fd ,i 0 0 0 4 01 0 4) 4 "a 0 I X 4i rH 0nJ. ,t i FokI (Cl ( a E (U ^ ~i~ I K i N I IN., I,, 4 I I I I I I I I I I I I I I I *V 'sn!pDoI 0 0j 03 I I 0 oC 41 4a) 0 rJ 4Q) , ( 1) I00 0 ot ao H 0) o ! k P 4) 4 0 0 C I En U) ic cd In ^ 10 ( 3 41  CT>^ > *rl 3 cr M T! N I I  I I I I I I I I I I I I I I "  p  .V 'snipDtj 0 I I I I I I I I  W tO ." rnO N 0 co lI L. CD 0 i 0 t i o < 0 r. C)c 0a 0 0 CC 0 ) 0 0 1 H U ?i3 ca n) +J inn 4> (UUo &< (fl TS Nt N) I i I I I I I i i I I I I I I I i i S '  *V 'snlpDj CL LCO0 0L L CD 4) 0 & 0 0 m 0 EU)~ i 01 p p0 M ( 0) 4( o P E44 0) oa .,I H 41U z 0 0) (a a) H 0 to 01g P ".4 *il *C *r( .LO) pIl 0 0 p 9) ".4 >1 4j 44 4 M 10 0 0) 3 C C U 0 (C(U >r 0) ro + B 0 BS (0 M r) ( P 0 *O k ^ tn3 ( "< >ii r ^3 i Although no group has its mean position in the center of the micelle, this does not indicate a void within the micelle. Rather, it is an artifact of the coordinate system chosen. Each chain segment has a diameter of 3.56 Angstroms, and subsequently excludes a spherical volume of this diameter to any other group. Since the spherical coordinate system provides an available volume which decreases with r, the excluded volume effect results in low values of Ni(r) near the center of the micelle. Thus, no group has its average position near the center, but the center is always within the excluded volume of one of several different groups. A range of three standard deviations about the mean includes approximately all of a group's positions, and, in Figures 51 through 54, approaches within a segment diameter of the center for several chain segments in each simulation. The mean position results for Runs 1 and 2 are virtually identical. This is to be expected, since the two models differ only in the mass of the head group, a difference that should manifest itself in the dynamic behavior of the system, but not in a timeaveraged result. In these two simulations, the head groups are found at the surface of the micelle with a relatively small deviation from the mean. Progressing along the chain toward the tail, the mean radial positions decrease and the deviations remain roughly constant through group five. The mean positions of the last four groups are at approximately the same radius, but toward the tail the deviations increase, a consequence of the excluded volume effect near the center of the micelle. The results for Run 3 are markedly different. Figure 53 shows larger values of mean positions for the chain segments and greater deviations for the head group and its adjacent three groups. This simulation differs from Run 2 in the size of the head group and the nature of the headchain bond. With a diameter of 8.72 Angstroms, the larger head group creates a steric effect which brings about greater disorder in the micelle. The head groups are spread over a greater range of radial positions, altering the positions of the chain segments from those observed in Runs 1 and 2. Figure 54 shows the results for Run 4, the hydrocarbon droplet. The mean positions and deviations are symmetric about group 5, the center of the molecule. This is the proper result, since the two ends of a chain are indistinguishable from one another. The chain ends have a slightly larger mean position and deviation, a result of the greater mobility of the chain end relative to the interior groups of the chain. Overall, the means and deviations are quite uniform, indicating a random arrangement of the molecules within the droplet. Again, the excluded volume effect places the mean positions of all groups over two segment diameters away from the center of the droplet. In a simulation of a micelle of forty dodecyl surfactant chains, the same shell model was used to surround the micelle and the head groups were identical to Run 1 in size, mass, and potentials (Woods et al., 1986). Comparison with this simulation reveals the effect of chain length on the micelle model. In Figure 55, the average radial positions of Woods' simulation are shown next to those of Run 1, pairing the groups of Run 1 with the nine groups furthest from the head group of the longer chain's thirteen. The dodecyl micelle being larger than those of the present work, its average radial positions are greater. Qualitatively, however, the trend observed in the radial positions of the groups and their standard deviations is quite similar in the two simulations. In Figure 56, the average radial positions of Run 1 are compared to those resulting from a statistical model of 30 twelvemember chains with one end of each fixed at the surface (Gruen, 1981). In this model, all possible configurations of the surfactant chains were sampled. Though Gruen's micelle is again larger than that of Run 1, the results of the two models show notably similar radial position behavior on a qualitative basis. 0i ? n c / a (a I o u 0 0 0  0 , U 0 SCD 0 Ii5i cs _ _ 0 0 1o C, ^ 4 4 U 02 0 I "I 0 1 I t Dr  , ,i , ii0 V 'Sni p  0 I0r ( O t n CC 0 ( O r 0rC I p3 NZ I I I:B ' ' I I I I I I I I I I 0 M0 I' O O.4 ) N " 0 "V 'snipoD 0  C (I C  I c I c I  (NJ .r r  o C Tyi r 0 3 c c ri tt 0 0 0 . o o 10 4 i4 (C C) 4 *Hl 3 Q) h P U N I N I N I3 N Il  I I I I I I I I I I I I I I I I I  "  *V 'snipDc 0  cI . i  I  E tN  I C r;  1 c I I I I r'~c'4~O 0 (a 0 t u4 0 0r= o 4o m~i . 0 o11 r(C *ri 0 l1 (0 C4 0 o E Lu 0 r. 4 ) F ) C,4 1 F 3 3 rt a &'1r 1 A model micelle of 15 sodium octanoate molecules, along with surrounding molecular water, was the subject of another molecular dynamics simulation (Jonsson et al., 1986). The sodium carboxylate head group was explicitly modeled and surrounding water molecules were included in the simulation. Jonsson found it necessary to reduce the charge on the head groups to produce acceptable results. Run 3, with the model sulfate head group, is compared to Jonsson's reduced charge simulation in Figure 57. Accounting for the difference in the size of the two micelles, the radial positions of the chain segments compare favorably. 5.2 Probability Distributions of Group Positions The positions of groups within the micelle can be studied in more detail by evaluating, at each radial position, the probability of finding a given group at any given time. The probability density function used in this analysis is the timeaveraged number of occurrences of the center of a group i within a spherical volume element centered at r: p,(r) = 4nr2Ar Evaluating this function over all values of r generates a probability density distribution for each molecular group. The results of this analysis for each of the four simulations are given in Figures 58 through 511. The L m C9'S 00 0 0~ C)3 tc aa,  EC 0D 4) a) En 0 ~ ~ S Z 0]4 411U CC 0 r. 4 CD) 4 U3 l a^^Q a) p sa :J 4 tp F LO ~ ~ ~ ~ ~ ~ ^ d n C4 0 C 0 r CO O K K) y q i 84 q~/v^ )^ Z;; q q q qF40( 0 ~~"^ T 0 0 0 o> a^o a a ao a) a ^~~ I k ( v ,d4 ao pi > ~ o  Q ) _________________""^ s t  I 1 I I \ I I     I I I I 3  in ^ rO CM ^ O OT OO ^ co n '+ rO CM ' OS' r' ' . '  0 0 0 0 0 0 0S o r o o o o oo oo o o o omur d o d d d d d o o d o *v! O *d C a 0< C.1 0 o Q" C4  0 0 00 dd tn 0 0 C QC r. c44W 0m 0 i 4). S.l M 04 r 4 0 0 3 ( 4) 0 ptt rl >fa 3) p 44 0 ( Fr el r) c3 cN  M c r0 ED 4O  f S0 00 0 0 0 0 66 o 6o 6 d o d o o o "v ",d '1 CN 0 0r  o O Co I, cO Lo  r') N 0   0 0 0 0 0 0 0 0 0 0 .v '!d 0 C. 0 4J 44 0i w <0 m 0 04C o 0 03 4r m1 24 44 . cn il 3 Eo~ 00 d d 5 5 o d 