Thermodynamic modeling and molecular dynamics simulation of surfactant micelles

MISSING IMAGE

Material Information

Title:
Thermodynamic modeling and molecular dynamics simulation of surfactant micelles
Physical Description:
vii, 276 leaves : ill. ; 28 cm.
Language:
English
Creator:
Farrell, Robert Anthony, 1956-
Publication Date:

Subjects

Subjects / Keywords:
Micelles   ( lcsh )
Surface active agents   ( lcsh )
Micelles -- Computer simulation   ( lcsh )
Chemical Engineering thesis Ph. D
Dissertations, Academic -- Chemical Engineering -- UF
Genre:
bibliography   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1988.
Bibliography:
Includes bibliographical references.
Statement of Responsibility:
by Robert Anthony Farrell.
General Note:
Typescript.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001108624
oclc - 19761685
notis - AFK5070
sobekcm - AA00004800_00001
System ID:
AA00004800:00001

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
Single-Component 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 SINGLE-COMPONENT
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

seven-step 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 single-component 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 two-fold 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 single-component 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 models--three micelles and

one hydrocarbon droplet--is 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., polar-nonpolar) while maximiz-

ing the favored contact (i.e., polar-polar). 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 water--a favored interac-

tion--while the tails reduce their contact with the

water--an 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

size-dependent.

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 small-systems 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; Ben-Naim 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 semi-empir-

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 density--a 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 seven-step process shown in Figure 2-1.

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 all-trans chain. With the addition of the next




















C )=
C1) !-4
01 Q)O


4-r-i r-i 4
00 4J-H

000 au)
-r-i 0 -ri r3 Z 0
4J -r-.4 0) 4-) 0
-r-1 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) r-4 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 r-1 0

S-4+J o uM 0 tn
U2 0 4-4 U] V0 FO

4 L o p0 m
.0 0 4 z w
I-- 4-C I 0
1 3 ) 4- Z
4 4-) 4-1M (1)
0) -H U a) 0) U U
p-l 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 material-free 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, rod-shaped, and even disc-shaped, 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 hydrocarbon-sol-

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 weak-to-moderate dipole

moments. The dipole-dipole 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 micelle-free solution and

the final Poisson-Boltzmann distribution of counterions

around the surface of the micelle with a fraction bound in a

Guoy-Chapman 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 n-alkanes, 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 hydrocarbon-water 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 I-dimensional

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(01-CMC 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 Ben-Shaul

and coworkers (Ben-Shaul, Szleifer and Gelbart, 1985;

Szleifer, Ben-Shaul 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-- N-* c)N











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 3-1 and 3-2 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 3-3. 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 ) 4-4
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 4-4 1 1i X






04-4 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- 4-3 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
-- -e-l











.- I-r4 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 ( ^
*^ kl-la,





















0







C )
0 >4 l
-H E-I ,J

o 0 C0



C0 C
0' -H (o







( oo

r. 0
0 4-4 9
Q H 03
Cr z a)



gIn -H
C: ) -H 0)



N4J
4-)



-H


0 0 4










0 o oo o0 ,r
S4-I




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 Single-Component
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 3-4. 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.0E-04 a 73 a

308.15 6.6E-04 260

318.15 6.4E-04 640

C12H25(OC2H4)60H 298.15 8.7E-05 b 400 c

C14H29(OC2H4)60H 298.15 1.1E-05 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
4-C


0 a) 4-



n 0 -0 C
0) r--!
e w u z


SC~







0-4 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 E-i 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= .0549T-20.20

n=-.0137T+ 1.421

d =.1 125T-45.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)T-45.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,)nc-27.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

3-5 through 3-8. 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 3-5 and 3-8 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




U-4 in
UC>








C0
faQ

o >


0) .4^ U
C -3

0 )
4-4; i

r.?

S-i~
U

1 C)(


4c) 4
0 0c U





0 *C 03
O J S-




4-Jo
0 i -o-l







0- (a 0-S.u
0 CL a)

4- 0
0) 0 -1 0






4 05-
CC- Q) 44S


>(0 04-C









:- 41 1 (Z
1 Q) 4C 4
-P 0r














-H 4-3 U) (0
4C (a -
rcu



nU a la)a
0 *P


CP ) ^r
-*- 4- i
fE (0 3


















c a0
S-o 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



4-1
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-) -
1-rr1

CE0



0


C m
U









r4 0
ca, U

4c-) 04







cra)
0 *P (








0 0



4--t 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 3-9

and 3-10. Below the CMC, the model calculates a mean

aggregate size approximately equal to Ntrans, present as the

initial plateau in figure 3-10. 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 rr-4


_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
0--4

I I n IIII I -





4-4
) 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 N-C 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 3-11. 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 3-12 shows the

concentration derivative of equation (2.21) for the data of

Figure 3-11 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 r-o
I (U 0
Li D '7
o L. U C
o --i U




0 -0-3

OZ 0


o
0 U

iC 4-0 *
O L

J ir-



0 r- 0 m
i < ,-i0
O -n

o o-i
01 0
-0 4 PH 4-0

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 chain-length

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 3-13 and 3-14. 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 >








r-4
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


4-4
1 2


.) u U I I "- 0

I I "


II --


4.4 -r -

) 0






Q)2)
c
0 a) -0-u
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 3-15 and 3-16. 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'N-X)- m-


c0







S-4



0
03







4 30
C




4 U








0 J-






.c 3

E-4
cO












a) m


-,. '0.
Q)r c1
g 3I
=) 0c
Fp
i-c
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











small-angle 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).

Quasi-elastic 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 second-order 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 fifth-order

predictor-corrector 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 predictor-corrector 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 10-20 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 dr r-r


U,(r) = U(r) U(r,) -(r r L dr Jr-r I

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 4-1.

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 4-1. The model surfactant molecule, consisting of an
eight-member alkane chain attached to a polar head group,
shown in the all-trans 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 n-butane. 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. b-b -, (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 cosO6-cosB (-- -::
b 1 (bb I- b

~ +l [b b S(4. 1 .)












ib2(Oi)=Yb(COS oCOS ) b,+ 1 bi l b
,-ybtco-cos 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.116-1.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 4-2. 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) 4-3



/ '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 4-J >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 rmr-lP r.
(3 -0 014)
p O (a 4 0
3 I C *"-4

n .t .e) 14 ia ) (04-4 ) 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 group-shell

2. Alkane group-shell

3. Head group-head group

4. Alkane group-alkane group

5. Head group-alkane 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-

tial-field 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
chain-solvent 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(ri-ro)2 (4.22)

F(ri)=-2yh(ri-r)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 rw-ri, rw being the radial
position of the spherical repulsive shell, and rm and E are
the radius and energy of minimum potential for the
alkane-alkane intermolecular potential.
The interaction between two head groups is described by
a potential comprised of both dipole (r-3) and soft-sphere
(r-12) 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-
ed-force 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:

head-alkane 14.
m 2(rm +rhh) (4.34)

Figure 4-3 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 all-trans

conformation. With the bond rotation force constant reduced

to one-tenth 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 4-3. The intermolecular potentials for segment-seg-
ment (6,9), segment-shell (12), and head-head (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 +

Head-shell



Segment-shell



Head-head



Segment-segment



Head-segment


YV

bo

Yb

0o

Yr

Yh

lo

E

rm

E

rfhh



rm



h-a
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 nine-segment molecules by setting

the head group to the size and mass of a methylene and

replacing head-head and head-shell 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 all-trans configuration. The groups

were given initial velocities, the rotational force constant

was reduced to one-tenth 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. Ten-thousand 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- r-4 H c4-)
5 I 0-ay



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 1-4Cl to cp
r-i (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 <- r-l 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, bond-angle 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,

neighbor-listing 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 Vector-optimizing 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

eight-byte 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 variable-type 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

time-averaged 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 rdr (5.1)


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 -[- (r-R,)2dr (5.2)

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

5-1 through 5-4. 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
4-i -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


4-1
4a)


0
rJ


















4-Q) ,


( 1)

I00

0
ot
ao


















-H 0)
o !
k -P




















4) 4


0-
0 C












I En
U)
















ic cd
In
^ 10 (
3 4-1 -
CT>^ >
*r-l 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






































N--t--







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


4--1U
z 0 0)

(a a) -H

0 to






01g
P ".4
*il *C *-r(














.LO)




pI-l 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 5-1 through 5-4,

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 time-averaged 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

5-3 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

head-chain 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 5-4 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 5-5, 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 5-6, the average radial positions of Run 1 are

compared to those resulting from a statistical model of 30

twelve-member 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.





















0-i
?-





n c
/-- a (a






I--- o u

0 0






0

-- 0
,- --U-- -0









SCD 0
I---i--5------i 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 p--3----






N-Z


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 I--l-- -


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



4-o





m~-i









. 0
o11



r(C
*ri 0




l-1 (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 5-7. 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 time-averaged number of occurrences of the

center of a group i within a spherical volume element

centered at r:


p,(r) = (5.3)
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 5-8 through 5-11. 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 '- O-S' 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.
c4-4W

0m 0





i-

4-).
S-.l M





0-4
-r 4












0 0
3 (

4-) 0

ptt
-r-l









>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





4-4

0i w
<0




















m 0

04-C o





0 03
4r-









m1




























-24 4-4 .
cn il 3
Eo~


00
d d


5 5
o d