# Large scale molecular dynamics simulation of self-assembly processes in short and long chain cationic surfactants

## 1 Introduction

For the past two decades, amphiphilic systems have constituted a field of great interest, both from a fundamental and an industrial point of view . This is mainly due to the fact that they exhibit a rich phase diagram when mixed with water and other organic species. Their amphiphilic nature leads to the formation of self-assembled mesoscopic structures, for example spherical micelles at intermediate surfactant concentrations. At higher concentrations, or with the addition of electrolyte or organic compounds (co-surfactant), cylindrical wormlike micelles may be stabilized. These wormlike micelles confer visco-elastic properties to the fluid. Indeed, gels are frequently formed, reflecting the entanglement between worms, with viscosity dependent on temperature, salt concentration, etc. As the surfactant concentration increases, new structures appear, including vesicles and bilayers. The rheology of these amphiphilic systems is often analysed using a simple theory of self-assembly which combines both thermodynamics and geometrical arguments . According to this approach, the ability of surfactant molecules to form a particular type of structure is governed by the packing parameter $`P`$, which is the ratio of the volume of the hydrophobic tail $`v`$ to the length of this tail $`l`$ times the effective surface area per head group $`a`$:

$$P=\frac{v}{la}.$$
(1)   
Considering $`v`$ and $`l`$ to be constant for a particular molecule, the geometry of the self-assembled structure is controlled by the effective surface area $`a`$ of the head group (for example, the addition of electrolyte tends to screen the repulsive interaction between head groups, and so to decrease the effective head group area, leading to a transition from spherical to wormlike micelles).   
Molecular dynamics simulations have succeeded in resolving the detailed structure of spherical micelles, and their interactions with the solvent. These simulations have commonly involved 10 to 20 surfactant molecules with up to 1000 water molecules, and have been run for a few hundred picoseconds . Under these conditions, the dynamical properties of micelles are extremely difficult to study because the time scales of dynamical processes may vary from $`10^{11}`$ (presumed characteristic time of shape fluctuations) to $`10^4`$ seconds (relaxation time of some shear-induced structures). The length and timescales involved in self-assembly phenomena correspond typically to the domain of applicability for different mesoscopic models including lattice-gas , lattice-Boltzmann  and dissipative particle dynamics (DPD) . An outstanding challenge for these upscaling methods is the definition of a coherent link between the atomistic description of the sytem and the related meso and macroscopic behaviour. Some recent work  has shown the feasibility of such an approach based on a systematic coarse graining of the system in order to link DPD to the microscopic world. Smit *et al.*  earlier developed a molecular dynamics model using a simplified description of the component molecules. This allows a prediction of the size distribution of micelles which agrees qualitatively with both experiments and theory. Nevertheless, chemical specificity is not taken into account in this model. In this paper, we report a fully atomistic study of the structure and dynamics of micelles, based on trajectories calculated up to 3 nanoseconds. On this timescale, we can study some fast dynamical processes, such as monomer insertion or removal from a micelle as well as local growth and fragmentation for individual micelles, which occur much faster than the approach to global thermodynamic equilibrium. The kinetics of such processes may be modelled on the basis of a generalisation of Classical Nucleation Theory (CNT) . This theory considers clusters of different sizes, which exchange monomers with the surrounding medium (solvent + surfactant monomers), but neglects cluster-cluster interactions. It provides a description of the time evolution of the size distribution of clusters. A development of this model including inhibition phenomena has recently been published . The Becker-Döring equations of CNT are a special case of the more general discrete Smoluchowski coagulation-fragmentation equations which describe rate processes between two clusters of size $`r`$ and $`s`$ and a cluster of size $`r+s`$. Although the determination of the rate coefficients for each single process is not possible from a small number of molecular dynamics trajectories, we are mainly interested here in analysing the dominant mechanism (Becker-Döring or Smoluchowski) when systems are not at equilibrium. To the best of our knowledge, this work constitutes the first fully atomistic approach to such phenomena on this timescale, for systems containing up to 15000 atoms. This has been achieved with the conjugate use of large parallel computers and a highly scalable molecular dynamics program.

## 2 Description of the model

### 2.1 Surfactant and water molecules

Two surfactant molecules have been studied in aqueous solution, namely *n*-nonyltrimethylammonium chloride ($`\mathrm{C}_9\mathrm{TAC}`$), and erucyl *bis* \[2-hydroxyethyl\] methylammonium chloride (EMAC). The latter is known to form wormlike micellar viscoelastic fluids which find commercial use in hydraulic fracturing operations . The two surfactant molecules are shown in figure 1. In some simulations, salicylate co-surfactant molecules have been added to the EMAC solution. The force-field used to model these surfactant molecules is the “constant valence force-field” (cvff) , which includes explicitly all the atoms. The total energy is the sum of the intramolecular and intermolecular interactions. The intramolecular interactions are represented as a sum of four types of term:

$$\epsilon _{intra}=\epsilon _{stretch}+\epsilon _{bend}+\epsilon _{torsion}+\epsilon _{outofplane}.$$
(2)   
The explicit forms for each term in eq (2) are given in eqs (3-6). The bond stretching term is given by:

$$\epsilon _{stretch}=\underset{i}{}k_i^b(bb_0)^2,$$
(3)   
where $`k_i^b`$ is the bond stretching force constant, $`b_0`$ is the equilibrium bond length, and $`b`$ is the actual bond length. The bond bending term is given by eq (4) :

$$\epsilon _{bending}=\underset{i}{}k_i^\theta (\theta \theta _0)^2,$$
(4)   
where $`k_i^\theta `$ is the bond bending force constant, $`\theta _0`$ is the equilibrium bond angle, and $`\theta `$ is the actual bond angle. The torsional contribution to the intramolecular energy is represented by a cosine series:

$$\epsilon _{torsion}=\underset{i}{}k_i^\varphi (1+Scos(n\varphi )),$$
(5)   
where $`k_i^\varphi `$ is the torsional force constant, S is a phase factor (equal to $`1`$ or $`1`$ depending on the dihedral angle considered), and $`\varphi `$ is the torsional angle. The out-of-plane term describes the resistance to out-of-plane bending and is expressed by a quadratic distortion potential function:

$$\epsilon _{outofplane}=\underset{i}{}k_i^\chi \chi ^2,$$
(6)   
where $`k_i^\chi `$ is the bending constant and $`\chi `$ is the bending angle. The expression for the non-bonded interactions is:

$$\epsilon _{inter}=\epsilon _{vdW}+\epsilon _{coulombic}$$
(7)   
where the summations are performed over all the non-bonded pairs of atoms. A Lennard-Jones 12-6 pair interaction is used for the van der Waals energy, $`\epsilon _{vdW}`$, and the partial charges involved in the coulombic term were computed using the semi-empirical quantum mechanical program MOPAC , within the AM1 approximation. The charge on the $`N(CH_3)_3`$ head group of the $`\mathrm{C}_9\mathrm{TAC}`$ surfactant was found to be $`+0.88`$, and $`+0.90`$ for the EMAC head group $`N(C_2H_4OH)_2CH_3`$. Water molecules are represented using the Jorgensen  TIP3P model: interactions between water molecules are described by a Lennard-Jones potential between oxygen atoms and electrostatic contributions between all atoms (hydrogen and oxygen). All the parameters of the TIP3P model are shown in table 1. In order to check the performance of the TIP3P water force-field in predicting the bulk properties of water, we performed a molecular dynamics simulation of liquid water and compared the computed properties such as water density, diffusion coefficient and radial distribution function between oxygen atoms to the corresponding experimental data. More details about the simulation procedure can be found in ref . The calculated bulk water density and self-diffusion coefficient as averages over the stored time series of particle coordinates are listed in table 2. Both values compare well with experiments.

### 2.2 Molecular dynamics method

Several distinct initial configurations of surfactant molecules surrounded by water molecules were constructed. They consisted of an infinite wormlike micelle, a spherical micelle, or a random distribution of surfactant molecules. Forty-eight or fifty surfactant molecules were employed in the cases of $`\mathrm{C}_9\mathrm{TAC}`$ and EMAC respectively. The number of surrounding water molecules was approximately equal to 3000 in each simulation. In some cases, electrolyte (NaCl) and/or co-surfactant (salicylate) was also added to the solution.   
At the beginning of the dynamics simulation the total potential energy was minimised in order to generate a reasonable starting point. To carry out the minimisation, we used a truncated Newton-Raphson method requiring evaluation of the second derivative of the potential energy with respect to the atomic coordinates. After minimisation, random velocities selected according to a Maxwellian distribution at a temperature of 300 K were assigned to each atom. The pressure was set to 1 atm and the temperature was fixed at 300 K. Pressure and temperature were controlled using the Nosé-Hoover algorithm . In most of our studies, the total simulation time was larger than 1 ns. To integrate the Newtonian equations of motion for all atoms, we used the Verlet leapfrog algorithm with a timestep of 1 fs. Periodic boundary conditions were applied in all three spatial directions and Ewald summation was used to handle the long-range electrostatic interactions in conjunction with the particle-particle / particle mesh (PPPM) method, an $`𝒪(NlogN)`$ algorithm . A cut-off radius of $`10.0`$ $`\mathrm{\AA }`$ was used for non-coulombic interactions.

### 2.3 Parallel implementation of molecular dynamics

All MD simulations were carried out either on a Silicon Graphics Origin 2000 or on a Cray T3E using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code . LAMMPS is a highly scalable classical molecular dynamics code designed for simulating molecular and atomic systems on parallel supercomputers. To study large systems of molecules for a large number of time steps, an algorithm is required that has a very good speedup with the number of processors used. This speedup has been calculated up to 1024 processors on a 1500 node T3E, and the results display the desired linear scalability property. The same calculations have also been performed on a Silicon Graphics Origin 2000 using up to 32 processors. Results of the benchmarks are displayed in figure 2, where the speedup is given by:

$$\mathrm{speedup}(\mathrm{k})=\frac{\mathrm{time}(1\mathrm{processor})}{\mathrm{time}(\mathrm{k}\mathrm{processors})}.$$
(8)   
One can see that the parallel performance of the LAMMPS code is superior to that of the MD codes in the commercial package $`\mathrm{Cerius}^2`$ . The superlinear behaviour observed in the case of LAMMPS, related to the efficient calculation of the coulomb interactions and the use of spatial domain decomposition, can be understood in term of cache memory utilisation which is sub-optimal for one (or a very small number) of processors.

## 3 Structure and dynamics of $`𝐂_9\mathrm{TAC}`$ micelles

The instantaneous configurations of the system are displayed using the MSI $`\mathrm{Cerius}^2`$ package . Some images employ Connolly surfaces  calculated specifically for some molecules or fragments. A Connolly surface is the van der Waals surface of the molecule/fragment that is accessible to a solvent molecule. In all cases studied in this paper, the solvent molecule is a water molecule. The blue square appearing in each snapshot represents the simulation box.

### 3.1 Spherical micelle

A spherical micelle was built with 48 $`\mathrm{C}_9\mathrm{TAC}`$ surfactant molecules surrounded by 2997 water molecules in a cubic box with sides of length 70 $`\mathrm{\AA }`$. It took around $`50`$ ps of molecular dynamics at T=300 K for the volume of the simulation cell to reach a stable value $`V=1.1`$ $`10^{25}m^3`$, corresponding to a density of $`0.97`$ $`g/cm^3`$. The surfactant concentration in the simulated solution is thus $`C=0.72`$ $`mol/l`$. The values of the critical micelle concentration (cmc) for similar surfactant molecules like $`n`$-decyltrimethylammonium chloride ($`\mathrm{C}_{10}\mathrm{TAC}`$) or $`n`$-dodecyltrimethylammonium chloride ($`\mathrm{C}_{12}\mathrm{TAC}`$) in water are reported as $`0.05`$ $`mol/l`$  and $`0.0610.065`$ $`mol/l`$  respectively. Our simulations are performed at a much higher concentration than the cmc; micelles would thus be expected to form.   
Figure 3 shows snapshots of the system at different times during the simulation. Water molecules are not displayed in order to have a more detailed view of the structure of the micelle. The calculated values of the radius of gyration, the ratios of the lengths of the principal axes, and the radius of the micelle at different timesteps, are reported table 3. After $`50`$ ps, both the volume and the total energy of the system have reached stable values. Figure 3 gives evidence of the spherical shape of the micelle. The polar head groups are located on the micellar surface and are in direct contact with water molecules; the alkyl chains are directed into the hydrophobic core. Analogous views of the system at $`600`$ ps, $`1.1`$ ns and $`3`$ ns are displayed in figure 3. After $`600`$ ps, we observe that 3 surfactant molecules have left the micelle. At this stage, the micelle contains $`45`$ molecules. Its shape is ellipsoidal rather than spherical, as can be deduced from the ratios of the lengths of the principal axes (see table 3). After $`1.1`$ ns of simulation, the micelle has broken down into two smaller spherical micelles, one containing $`29`$ surfactant molecules and the other one $`15`$ molecules. The four remaining surfactant molecules are solubilized as monomers in the water. The observed mechanism through which the micelle breaks into two smaller entities is as follows: the initial micelle undergoes a structural change from spherical to ellipsoidal, one of the principal axis becoming twice as long as the other two. This ellipsoid looks like a small dumbbell (the density of surfactant molecules at its middle is small). Finally, the dumbbell separates into two spherical entities. A rapid reorganisation of the surfactant molecules occurs after the separation.   
Experimental results of Imae *et al.*  on several similar surfactants in aqueous solution indicate the presence of spherical micelles with average aggregation numbers of $`84`$ for $`\mathrm{C}_{16}\mathrm{TAC}`$ ($`n`$-hexadecyl trimethylammonium chloride), $`62`$ for $`\mathrm{C}_{14}\mathrm{TAC}`$ ($`n`$-tetradecyl trimethylammonium chloride), and $`44`$ for $`\mathrm{C}_{12}\mathrm{TAC}`$ ($`n`$-dodecyl trimethylammonium chloride). An extrapolation of their results would predict an average aggregation number of approximately $`25`$ for $`\mathrm{C}_9\mathrm{TAC}`$. At the end of the simulation (figure 3), several surfactant molecules have left the larger micelle while the smaller one has increased its size by adsorption of one further monomer. The two micelles finally contain $`24`$ and $`16`$ surfactant molecules, with $`8`$ isolated surfactant monomers remaining in the solution. Three nanoseconds of molecular dynamics simulation are clearly not long enough to guarantee that the system has reached thermodynamic equilibrium. At equilibrium, spherical micelles are known to exhibit a size distribution rather than one particular aggregation number. The size of a micelle thus evolves in time, over a range given by the polydispersity of the size distribution function. In this simulation, one micelle has adsorbed one surfactant monomer while the other micelle has desorbed a few monomers.

### 3.2 Infinite cylindrical micelle

A cylindrical micelle containing $`48`$ $`\mathrm{C}_9\mathrm{TAC}`$ surfactant molecules was constructed by putting together $`6`$ discs of $`8`$ molecules each. This cylindrical micelle was surrounded by $`4305`$ water molecules. The length of the initial cylinder was 30 $`\mathrm{\AA }`$ and 3D periodic boundary conditions were applied in order to simulate infinite cylinders. Figure 4 shows the initial configuration (viewed perpendicular to the axis of the cylinder), and instantaneous configurations taken after $`200`$ ps, $`700`$ ps and at the end of the simulation ($`1.1`$ ns); it takes approximately $`80`$ ps to reach stable values of the volume of the simulation cell and the total energy.   
The micelle evolves in time and some distortions from the initial configuration appear quickly. After $`200`$ ps of simulation, the micelle is still cylindrical, but the density of molecules along the axis of symmetry of the cylinder no longer appears to be uniform: the surfactant cations clump together, forming high density regions while other regions of the worm exhibit a small number of surfactant molecules. The latter correspond to “weaker” points, i.e. preferential zones for fragmentation of the worm, with weaker hydrophobic tail interactions. Such a phenomenon signals the incipient break-up of the infinite cylinder, as is clearly seen in figure 4 after $`700`$ ps of simulation: the initial cylinder has expelled a small spherical micelle comprising $`15`$ surfactant molecules. The characteristic value of this spherical micelle’s radius of gyration and the lengths of its principal axes are reported in table 3, providing evidence of its spherical shape. The rest of the cylindrical micelle contains $`30`$ cations and exhibits a non-spherical structure. At the end of the simulation ($`t=1.1`$ ns), the small spherical micelle contains $`14`$ surfactant cations, and the other one $`30`$ monomers. The remaining four monomers are located in the aqueous solvent. The different characteristic shapes of these micelles are reported table 3. This state is similar to the state reached in the previous simulation after fragmentation of the spherical micelle into two micelles of sizes 15 and 30. We thus conclude that the time evolution of this state would probably produce the same features as those described previously, tending to thermodynamic equilibrium.

### 3.3 Infinite cylindrical micelle with electrolyte

The addition of electrolyte to an aqueous solution of surfactant micelles is known empirically to preferentially stabilize the cylindrical shape . For a cationic surfactant, the negative ions of the added salt associate with the positively charged head groups of the surfactant. Such associations reduce the strong electrostatic repulsion between neighbouring head groups that exists in the cylindrical structure where the latter lie closer together than they do on the surface of a sphere (thus the packing parameter $`P`$ increases due to a decrease in the surface area per head group $`a`$). Our aim here was to investigate whether $`\mathrm{C}_9\mathrm{TAC}`$ cylindrical micelles could indeed be stabilized by adding an electrolyte.   
The model was constructed in the same way as the previous one for an infinite cylindrical micelle but now $`\mathrm{Na}^+,\mathrm{Cl}^{}`$ ion-pairs were added to the solution. Specifically, the system was composed of $`48`$ surfactant cations (and $`48`$ chloride counterions), $`2997`$ water molecules, and $`100`$ sodium chloride ion pairs, each ion being placed at random within the water molecules of the solvent. The concentration of chloride anions is thus $`C=1.64`$ $`mol/l`$. A view perpendicular to the principal axis of the cylinder is shown in figure 5. After $`100`$ ps of MD simulation (figure 5), the infinite micelle has broken down into a finite micelle. Its shape remains roughly cylindrical, as evidenced by the ratios of the principal axes reported in table 3. In order to compare the stability of this structure with that of the spherical micelles produced in the absence of electrolyte, the simulation was run up to $`2.5`$ ns. At the end of the simulation (figure 5), the micelle still exhibited a cylindrical shape, containing $`44`$ surfactant molecules (four monomers left the cylinder and were individually dissolved in the water): it is an example of a finite, rod-like micelle (self-assembled structures which have been observed experimentally). However, we cannot confirm that this system is near equilibrium, owing to the lack of either longer-time simulation data or relevant experimental results.

### 3.4 Random monomer distribution

A starting configuration of $`48`$ surfactant and $`2997`$ water molecules was prepared, the surfactant cations and chloride counterions being randomly distributed throughout the system (see figure 6). The simulation cell was a cubic box with sides of length $`70`$ $`\mathrm{\AA }`$, corresponding to an initial surfactant density of $`0.3`$ $`\mathrm{g}/\mathrm{cm}^3`$. The simulation was run up to $`900`$ ps in order to obtain information on the mechanism of surfactant micelle formation. At $`t=75`$ ps, the system has already become inhomogeneous with some low and high density regions (figure 6). In the high density regions, no particular arrangement of the surfactant molecules is discernible. A structural aggregation between the surfactant particles is seen at $`t=200`$ ps (figure 6), leading to the appearance of two micelles. Each micelle is composed of approximately $`15`$ molecules. By examining several instantaneous configurations between $`t=75`$ and $`t=200`$ ps, we can obtain insight into the dynamics of spherical micelle formation at the molecular level. It appears that in the first stage (between 0 and 100 ps), the surfactant molecules approach one another, forming aggregates without any well-defined organisation, both from a translational and rotational point of view (figure 6). The size of these disordered micelles is small ($`10`$ molecules). In the second step these random aggregates rearrange to form spherical micelles. The driving force for this rearrangement appears to be the minimisation of the repulsive interactions between head groups together with the hydrophobic attractions between the hydrocarbon tails. During this structural rearrangement, the aggregation number of the micellar clusters increases via the addition of small clusters of surfactant (typically 2 or 3 molecules).   
Finally, by the end of the simulation, both micelles exhibit a spherical shape, containing $`15`$ and $`17`$ surfactant molecules respectively. Their characteristic radii are reported in table 3. A small number of monomers remain solvated in the surrounding aqueous medium. The sizes of the micelles are consistent with the assumption of a mean aggregation number of around 15-20.

## 4 Structure and dynamics of EMAC micelles

### 4.1 Spherical micelle

A simulation cell was constructed, comprising $`40`$ EMAC molecules within a spherical micelle, together with $`8`$ solvated monomers, $`48`$ chloride counterions and 3497 water molecules. Molecular dynamics was performed on this system up to $`1.05`$ ns. It took $`120`$ picoseconds for the volume of the simulation cell and the total energy to reach stable values (the surfactant concentration is equal to $`0.54`$ $`mol/l`$). At this point in the simulation, the characteristic dimensions of the spherical micelle are as reported in table 4. During the entire simulation, the micelle maintained its spherical shape. Contrary to the case of $`\mathrm{C}_9\mathrm{TAC}`$ spherical micelles, no fragmentation occurred within this timescale. This is possibly due to the fact that the expected aggregation number for a spherical EMAC micelle may be greater than $`40`$ . Nevertheless, two surfactant molecules initially dissolved in water approached the micelle and adhered to its surface. This phenomenon is shown in figure 7, where the two adsorbing molecules are highlighted. We can see that the two molecules remain close to each other like a dimer, even after their adsorption at the micellar surface. Their hydrophobic tails are not directed toward the centre of the micelle; they rather remain in close proximity to water along their entire length. Figure 7 shows a view of the final configuration after $`1.05`$ ns. The shape of the micelle is spherical, and 6 surfactant molecules remain dispersed in water.

### 4.2 Infinite cylindrical micelle

A cylindrical micelle of EMAC molecules was constructed by combining five discs of 10 molecules each. Due to the periodic boundary conditions, the cylinder is effectively infinite in the direction of its principal axis. This micelle was surrounded by $`3119`$ water molecules and $`50`$ chloride ions. The system took $`100`$ ps to equilibrate (that is, to attain stable values of the simulation cell volume and the total energy). Over a total period of $`1.85`$ ns of molecular dynamics, the micelle retained its cylindrical shape. No surfactant monomers were present in the surrounding water at the beginning of the simulation and neither desorption of monomers, nor any other form of fragmentation was observed: the micelle remained an infinite cylinder throughout. Figure 8 shows a projection of the system on a plane perpendicular to the axis of the cylinder, at the end of the simulation. Water molecules are also shown. The cross section of the micelle is more ellipsoidal than circular, as can be deduced from the values of the ratios of the lengths of the principal axes (see table 4). Moreover, no penetration of water molecules in the micelle core was observed; water molecules remain at the micellar surface, in close proximity to and surrounding the head groups. The core of the micelle is indeed completely anhydrous. The dimensions of a cross section of the micelle are reported in table 4.   
A snapshot of the system is also shown in figure 8, perpendicular to the axis of the cylinder, where three images of the periodic simulation cell are displayed. The micelle exhibits regions of varying density along its principal axis. Some local regions of the worm are narrower than others, corresponding to “weaker” points. With our enhanced understanding of the behaviour of $`\mathrm{C}_9\mathrm{TAC}`$ micelles, this could be interpreted as the incipient site of rupture of the cylinder into smaller spherical micelles. Nevertheless, we believe that this system is still at some distance from equilibrium, because the final state in this simulation is very different from the one obtained in the previous simulation (spherical micelle) while the compositions of both systems are more or less identical.

### 4.3 Random EMAC monomer distribution

A simulation was set up starting from a random distribution of $`50`$ EMAC surfactant molecules in a simulation cell containing $`3119`$ water molecules and $`50`$ chloride ions. Figure 9 displays the initial configuration of the system from which molecular dynamics was perfomed for $`1.3`$ ns. At t = 500 ps, two surfactant clusters can already be seen, each containing $`20`$ molecules; these are shown in figure 9. The Connolly surface has been calculated for all atoms of the surfactant molecules and is displayed in yellow. The first micelle does not exhibit a dense spherical shape, as can be seen by the presence of several voids in the structure. The second micelle is spherical, however, with a uniform density of surfactant molecules in all directions, leading to a compact structure. The structural differences between the two micelles are summarized succinctly in terms of the various characteristic geometric parameters reported in table 4. By the end of the simulation, no significant change occurs, and the final configuration, consisting of two micelles of sizes $`22`$ and $`20`$, is shown in figure 9. One micelle has grown while the other has kept its size unchanged; eight monomers remain elsewhere in the solution. The final geometrical parameters for the two clusters are also reported in table 4.

### 4.4 Infinite cylindrical micelle with electrolyte and co-surfactant

An infinite cylindrical micelle was built as previously described. It contained $`50`$ molecules, and was surrounded by $`3322`$ water molecules, $`150`$ chloride ions, $`109`$ sodium ions, and $`9`$ salicylate co-surfactant molecules. The co-surfactant molecules were placed outside the micelle. Molecular dynamics was performed on this 3D periodic system for $`400`$ ps. During this period, the cylindrical shape remained stable (as in the previous study without additional electrolyte and co-surfactant, see section 4.2), and its geometrical parameters are reported in table 4. As in the case without additional electrolyte, the cross-section of the cylinder is not circular but rather ovoid, one principal axis being markedly greater than the other. Figure 10 shows two views of the micelle at the end of the simulation. In the first one, all atoms in the system are displayed. The Connolly surface has been calculated for the co-surfactant molecules and is displayed in yellow. We can see -as in the case without added salt-that no water molecules have penetrated the micelle. A large number (six out of nine) of co-surfactant molecules have entered the micelle structure, being adsorbed on its surface and surrounded by the surfactant head groups. The second snapshot shows the micelle in a view parallel to the axis of the cylinder; the Connolly surface has been calculated for the atoms of the head groups only. The surfactant molecules seem to be homogeneously dispersed along the structure, and there is no appearance of weak points, or an incipient site of fragmentation. Nevertheless, we can see that large voids remain in the structure, providing evidence of direct contacts between the hydrophobic tails of the surfactant and the water molecules at the micellar surface. These voids are not associated with inhomogeneous density regions but are more likely due to an insufficient initial density of surfactant and co-surfactant.

## 5 Discussion of results

### 5.1 $`\mathrm{C}_9\mathrm{TAC}`$ surfactant micelles

For micellar systems, it appears that the characteristic time for full thermodynamic equilibration is generally greater than the microsecond scale. Nevertheless, the study of the behaviour of such systems on nanosecond timescales can provide useful information on the dynamics and structure of micelles. Indeed, some processes like monomer absorption or desorption, or the fragmentation of cylindrical or spherical micelles, can be simulated on this timescale. Starting from different initial configurations (spherical, cylindrical or random distribution of surfactant molecules) of essentially similar composition, we observe that these systems tend to evolve to the same state, which is characterised by the existence of two quasi-spherical micellar clusters of small aggregation number (between $`15`$ and $`20`$) and a few monomers isolated in the water. The radius of gyration of each micelle is roughly equal to $`11`$ $`\mathrm{\AA }`$, and the radii of the micelles (calculated as the averaged distance between the centre of the micelle and the nitrogen atoms) is $`7.58`$ $`\mathrm{\AA }`$. These micelles are non spherical *on average*, the average ratio of the length of the principal axes being significantly different from unity.   
The way this final state is reached starting from distinct initial configurations is different: in the case of a spherical or cylindrical micelle, we observe fragmentation into two smaller clusters. In the case of a random distribution of surfactant cations, we directly observe the *formation* of micelles. This process starts with the appearance of a disordered aggregate of surfactant molecules, of small size. A reorganisation of the surfactant molecules occurs, the head groups becoming anchored at the micellar surface. The time needed for some kinds of intramicellar rearrangements (following fragmentation) can be very small. As an illustration, a simulation was performed to investigate the rearrangements of surfactant molecules at the end of a finite rodlike micelle. The initial configuration was a finite rodlike micelle lacking end-caps, thus exposing extensive amounts of its hydrophobic interior to direct contact with water. In the first $`10`$ ps of simulation, the surfactant molecules move to create hemispherical end-caps. The time needed for several molecules to rearrange inside a micelle is thus very short, at least when high energy configurations are involved.   
In general, growth of micelles is achieved by the absorption of groups (typically dimers or trimers) of surfactant cations rather than by stepwise adsorption of monomers. We believe that this mechanism of micelle growth, corresponding to a Smoluchowski kinetic model:

$$C_r+C_s\begin{array}{c}\stackrel{\stackrel{a_{r,s}}{}}{\stackrel{}{b_{r+s}}}\end{array}C_{r+s}$$
(9)   
where $`C_r`$ is a micellar cluster of aggregation number $`r`$ and $`a_{r,s}`$ and $`b_{r+s}`$ are the forward and backward rate coefficients for aggregation and fragmentation processes, is valid when the system is far from equilibrium and at surfactant concentrations well above the c.m.c., because it leads to a faster rate to reach equilibrium than with a Becker-Döring scheme:

$$C_r+C\begin{array}{c}\stackrel{\stackrel{a_r}{}}{\stackrel{}{b_{r+1}}}\end{array}C_{r+1}$$
(10)   
i.e. a purely stepwise addition or removal of surfactant monomer . On the other hand, starting from a spherical or cylindrical micelle, we can observe its break up. In both cases, the initial micelle breaks into two clusters of different size, one containing approximately $`15`$ molecules and the other one $`30`$ molecules. The fact that a fragmentation process is observed confirms the validity of the Smoluchowski mechanism when the system is far from equilibrium. The size of the smaller micelle corresponds approximately to the equilibrium size of a $`\mathrm{C}_9\mathrm{TAC}`$ micelle, which remains essentially unchanged for the rest of the simulation. The size of the larger micelle does not correspond to the equilibrium size, and we observe a decrease of its aggregation number with time. This decrease is a stepwise process, following the Becker-Döring scheme. It is observed that the initial micelle produces two micelles of different sizes (one of which appears to be closer to the mean equilibrium cluster size) rather than producing two micelles of equal size, but larger than the equilibrium size. For systems not far from equilibrium, and/or for surfactant concentrations not greatly exceeding the c.m.c., a Becker-Döring scheme is more likely to be observed , characterized by a stepwise change in the aggregation number.   
Figure 11 shows the evolution of the radius of gyration of the two micelles (initially one spherical micelle) during the last 2 ns of simulation. The lower curve is associated with the smaller micelle which grows from $`15`$ to $`16`$ molecules during this part of the simulation. The radius of gyration of this micelle fluctuates around its average value of $`10.63`$ $`\mathrm{\AA }`$, with a standard deviation equal to $`0.29`$. This may correspond to the behaviour of a micelle near equilibrium. The upper curve is associated with the bigger micelle and arrows indicate the moments at which individual surfactant monomers leave the micelle. The main feature exhibited by the radius of gyration is a tendency to decrease, and this decrease is associated with the loss of monomers. Moreover, we can see that the curve exhibits oscillations of large amplitude (greater than a typical fluctuation). This corresponds in fact to an expansion-contraction process. The expansion of the micelle corresponds to an elongation in one direction; the surfactant molecules leave the micelle during the contraction process. This contraction indeed induces an increase in the repulsive interactions between head groups, leading to the desorption of one molecule.   
In order to highlight these shape fluctuation phenomena in the two micelles, the autocorrelation functions of fluctuations in the ratio of the lengths of the principal axes have been computed as :

$$C(t)=\frac{\delta R(t)\delta R(0)}{\delta R(0)\delta R(0)}$$
(11)   
where $`\delta R(t)`$, defined in , is:

$$\delta R(t)=R(t)R(t),$$
(12)   
and $`R(t)`$ is the ratio of the lengths of the principal axes. These are displayed in figures 12 and 13. The curve pertaining to the smaller cluster exhibits quasi-periodic oscillations, suggesting a time scale for the shape fluctuations of about $`50`$ ps. This agrees well with what has been previously observed in a 100 ps molecular dynamics simulation of a sodium octanoate micelle containing $`15`$ surfactant anions in water , where the time scale for a shape fluctuation was found to be equal to 30 ps.   
The curve associated with the larger micelle (figure 13) shows a different behaviour: periodic oscillations of large amplitude are seen, with a periodicity of $`500`$ ps. This time interval is associated with the slow expansion-contraction process described above. In this system, the fast shape fluctuation process has disappeared. The quasi-periodicity observed in these oscillations is possibly due to the fact that the rate of monomer desorption varies slowly with the cluster size.   
These results indicate that a shape fluctuation can be coupled with a desorption process, leading to a considerable increase in its characteristic time scale. The time scale of shape fluctuations in spherical micelles thus depends on their sizes (and therefore on the proximity to equilibrium).   
When electrolyte is added to the system, the dynamic behaviour is changed: the initially infinite micelle breaks up, loosing its infinite length. It thus becomes a finite rod, but still retains cylindrical symmetry. This small rod-like micelle is stable over the $`2.5`$ ns of simulation. The addition of salt is known to reduce the effective surface area per head group , leading to the stabilization of the cylindrical shape. The head groups, in a rod-like micelle, are thus expected to be closer together than in a spherical micelle. Figure 14 shows the radial distribution function between nitrogen atoms, calculated over the last nanosecond of the simulations described in section 3, in the three cases of a spherical micelle, a cylindrical micelle, and an initially random configuration. The three curves exhibit their main peak at roughly the same distance ($`9\mathrm{\AA }`$), while some differences can be discerned at shorter distances. The curves pertaining to the spherical micelle and the initially random surfactant configuration are very similar. Calculations of the mean total energies and volumes of the simulation cells confirm the similarity between these two systems. By contrast, the curve associated with the cylindrical micelle exhibits a more pronounced secondary peak at short N-N distances, providing evidence of a large number of close head-group contacts for this micellar structure. Figure 15 shows the radial distribution functions calculated between nitrogen atoms and respectively chloride ions, sodium ions, and oxygen atoms of water molecules. The short distances of the N-O (water) and N-Cl radial distribution functions indicate a high degree of structuring around the polar head groups. The closest distances correspond to nitrogen-oxygen interactions. Moreover, the appearance of a second (outer) peak in the N-Cl and N-O (water) radial distribution functions is associated with the existence of a second solvation shell surrounding the polar head group. This second peak in these curves is at the same distance in both cases. The sodium cations are located between the two solvation shells, but remain closer to the inner one.

### 5.2 EMAC micelles

As in the case of $`\mathrm{C}_9\mathrm{TAC}`$ micelles, various starting configurations were used to investigate the dynamical behaviour of EMAC micelles. In the absence of electrolyte or co-surfactant, this molecule is known to form spherical micelles . Starting either from a spherical or a cylindrical micelle, the system was found to keep its initial shape over a few nanoseconds of molecular dynamics simulation. In the case of wormlike micelles, although some evidence of incipient fragmentation was detected, the simulations were not performed over a large enough time to confirm this. It was also found that cross-sections of the worm were not circular, but rather elliptical as in the case of ellipsoidal (initially spherical) sodium octanoate micelles . The dynamics of EMAC aggregation and fragmentation processes is found to be slower than those for $`\mathrm{C}_9\mathrm{TAC}`$ self-assembly. This might be interpreted in terms of the length of the hydrophobic tail of the surfactant (from steric and energetic considerations, a long hydrophobic tail in a liquid would be more difficult to displace than a short one, and its diffusion coefficient is smaller). Nevertheless, we can see some differences between the spherical and the cylindrical micelles of EMAC molecules in comparison with $`\mathrm{C}_9\mathrm{TAC}`$. Figure 16 displays the radial distribution function between the hydrogen atoms of the hydroxyl group of the EMAC surfactant cation and the chloride anions. The properties of this function pertaining to the spherical micelle and the initially random distribution of surfactant cations are similar. There is a close contact between these two atoms at about $`2.5`$ $`\mathrm{\AA }`$, followed by an exclusion domain. The close proximity is indicative of a strong association between these atoms; while the exclusion domain is due to the electrostatic repulsion between chloride anions. Beyond the exclusion domain, the probability of finding a chloride ion reaches a value equal to that for the bulk concentration of electrolyte. In the case of the wormlike micelle, the first peak is higher than those corresponding to other shapes, and a second peak is clearly present. The surfactant head groups within EMAC in a wormlike micelle are thus seen to exhibit a stronger interaction with the counterions than in the case of a spherical micelle.   
Figure 17 displays the radial distribution function between the terminal methyl group and the oxygen atoms of water. We can see that in the three cases, there is a first peak at about $`4`$ $`\mathrm{\AA }`$. In the case of the cylindrical micelle, it will be recalled (see figure 8) that there is no water penetration inside the micelle. The short C-O contact distance seen in the radial distribution function displayed in figure 17 is thus due to the fact that the terminal methyl group is not located at the centre of the micelle, but rather is in direct contact with surrounding water. We can see that the location of the first peak and the probability of finding a short contact between the terminal methyl group and the water molecule both increase when going from the cylindrical to the spherical micelle, and finally for the curve resulting from the initially random configuration of EMAC monomers. These trends can be interpreted in terms of the effective surface area per head group, and the compactness of the resulting structure. In micelles with a cylindrical geometry, the head groups are closer to each other than they are in a spherical micelle , so that the compactness of the former structure is greater, leading to reduced water penetration. Moreover, in the case of the initially random monomer distribution, the two resulting micellar clusters contain half the number of molecules as in the case of the spherical micelle studied here, while their radii (see table 4) are roughly identical to those of the spherical cluster. The structure of the two small clusters is thus less compact, leading to a larger probability of close contacts between the hydrophobic chain and the water molecules. Even if there are some differences concerning the water penetration process in $`\mathrm{C}_9\mathrm{TAC}`$ and EMAC micelles, in all cases we have observed that the core of the micelle is completely impervious to water.   
In the case of the cylindrical micelle with added electrolyte and salicylate co-surfactant, the radial distribution function between the hydrogen atoms of the hydroxyl group of the surfactant molecule and the oxygen atoms of the co-surfactant molecule is displayed in figure 18. We can see that there is a first peak at about $`2`$ $`\mathrm{\AA }`$, providing evidence of hydrogen bonding between these two species. As was also noticed earlier (see figure 10), a large number of co-surfactant anions are integrated into the surface of the micelle. These anions do not penetrate the core of the micelle but rather remain adsorbed on its surface, enjoying a strong interaction with the surfactant head groups. This interaction decreases the effective surface area per head group, and enhances the formation of rod-like micelles, precisely as proposed in simple geometrical theories . The existence of hydrogen bonds between the head groups of the surfactant and the solvent is evidently not a necessary condition for the formation of micelles . Nevertheless, in this particular systems, it seems to play a key role in the attachment of co-surfactant molecules to the micelle. We have also observed the presence of hydrogen bonding between the hydrogen atoms of the hydroxyl group of the EMAC surfactant and the oxygen atoms of water molecules. However, there is no evidence of hydrogen bonding involving the oxygen atom of the hydroxyl group of the EMAC surfactant and the hydrogen atom of water molecules nor between two hydroxyl groups of different surfactant head groups. No interaction between two salicylate molecules has been seen.

## 6 Conclusions

Large-scale molecular dynamics simulations have been performed to investigate the structural and dynamical properties of self-assembled cationic surfactants in aqueous solution. One of the surfactants was comparatively short-chained ($`\mathrm{C}_9\mathrm{TAC}`$), the other long-chained (EMAC). The nanosecond regime has been reached, allowing the study of the dynamics of various self-assembly processes (growth and fragmentation of micelles, surfactant monomer insertion or removal). The mechanism of micelle formation at the molecular level has been described. We have interpreted separately the kinetics of these dynamical processes when the system is variously far from or near equilibrium. In the former case, a Smoluchowski-type scheme is obeyed, according to which micelles coalesce or fragment. In the latter case, a Becker-Döring scheme is observed , where only step-by-step monomer exchanges take place. It was also found that, for an oversized micellar cluster, the step-by-step elimination of surfactant monomers is associated with a slow expansion-contraction process of the micelle, with a characteristic time period of $`500`$ ps. On the other hand, a characteristic time period of $`50`$ ps for shape fluctuations was found in the case of a spherical $`\mathrm{C}_9\mathrm{TAC}`$ micelle with size close to the mean cluster size. The dynamics of the EMAC molecule is slower, and the equilibrium state has not been reached starting from different configurations. This difference in the time scales of the dynamics between the $`\mathrm{C}_9\mathrm{TAC}`$ and the EMAC molecule is attributed mainly to the size difference of the two cationic surfactant molecules. The effect of co-surfactant has been investigated and hydrogen-bonding with the head groups of the surfactant molecule has been found to play an important and probably a key role in stabilizing the wormlike assembly of EMAC cations. Finally, we have found that the penetration of water molecules inside micelles is not important in any instance examined, at least over time scales up to a few nanoseconds.

## 7 Acknowledgements

This work has been done in collaboration with Silicon Graphics Inc. who have provided access to a number of large parallel machines. Daron Green and John Carpenter are gratefully acknowledged for their generous technical assistance. Fruitful discussions with Trevor Hughes and Edo Boek are gratefully acknowledged. Mike Stapleton, Andreas Bick and Richard Painter of Molecular Simulations Inc. are also thanked for their support of this work.