Micelle Formation in Alkyl Sulfate Surfactants Using Dissipative

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2018, 14, 2633−2643

Micelle Formation in Alkyl Sulfate Surfactants Using Dissipative Particle Dynamics Richard L. Anderson,*,† David J. Bray,† Annalaura Del Regno,† Michael A. Seaton,† Andrea S. Ferrante,*,‡ and Patrick B. Warren¶ †

STFC Hartree Centre, Scitech Daresbury, Warrington WA4 4AD, United Kingdom Ferrante Scientific Ltd., 5 Croft Lane, Bromborough CH62 2BX, United Kingdom ¶ Unilever R&D Port Sunlight, Quarry Road East, Bebington CH63 3JW, United Kingdom ‡

J. Chem. Theory Comput. 2018.14:2633-2643. Downloaded from pubs.acs.org by UNIV OF WINNIPEG on 11/01/18. For personal use only.

S Supporting Information *

ABSTRACT: We use dissipative particle dynamics (DPD) to study micelle formation in alkyl sulfate surfactants, with alkyl chain lengths ranging from 6 to 12 carbon atoms. We extend our recent DPD force field [J. Chem. Phys. 2017, 147, 094503] to include a charged sulfate chemical group and aqueous sodium ions. With this model, we achieve good agreement with the experimentally reported critical micelle concentrations (CMCs) and can match the trend in mean aggregation numbers versus alkyl chain length. We determine the CMC by fitting a charged pseudophase model to the dependence of the free surfactant on the total surfactant concentration above the CMC and compare it with a direct operational definition of the CMC as the point at which half of the surfactant is classed as micellar and half as monomers and submicellar aggregates. We find that the latter provides the best agreement with experimental results. Finally, with the same model, we are able to observe the sphere-to-rod morphological transition for sodium dodecyl sulfate (SDS) micelles and determine that it corresponds to SDS concentrations in the region of 300−500 mM.

1. INTRODUCTION Used as detergents, foaming and wetting agents, in addition to dispersants and emulsifiers, surfactants can be found in a variety of everyday products. Surfactant molecules are amphiphilic, consisting of hydrophobic tails and hydrophilic head groups, and in aqueous environments self-assemble into clusters known as micelles to minimize interactions between water molecules and their hydrophobic tails. Micelle formation takes place only when the surfactant concentration exceeds the critical micelle concentration (CMC), and depending on the surfactant chemistry, composition, and concentration, a wide variety of micelle morphologies may be encountered such as spheres, rods, or meandering worm-like structures. These structures dictate the properties of micellar structured liquids, so that, for example, rod- and worm-like micelles are associated with a significant increase in overall viscosity. The formation of micelles is an essential process in many settings. For example, micelles are used to control drug release,1−7 can act as friction modifiers,8,9 and are responsible for the properties of many household cleaning products.10,11 Micelles even play a role in the human body where they transport fat-soluble vitamins.12 Experimental investigation of micelles is commonplace, and a variety of techniques are available to study their behavior. Nuclear magnetic and electron paramagnetic resonance, fluorescence, light scattering, and small-angle neutron diffraction can be used to monitor micelle formation, structure, and composition.13 Computer simulation provides a complementary approach to experimental methods (or indeed theory or empirical-based approaches) in the study of micellar systems. © 2018 American Chemical Society

These simulation approaches provide molecular-level resolution, determination of thermodynamic quantities, and details of the dynamics of micelle formation that can guide development of new systems and products. Atomistic simulations (e.g., molecular dynamics) have been used, for example, to determine the micelle structures14,15 and the free energy of micelle formation.16 Recent simulations on the micellization of sodium alkyl sulfates (hexyl, heptyl, octyl, and nonyl sulfates) have calculated equilibrium properties such as CMC and mean aggregation number (Nagg), although the method was found to underpredict these quantities compared with experimental observation.17 Surfactant self-assembly is difficult to simulate with all-atom methods due to the high computational cost associated with accessing appropriate length and time scales. Consequently the self-assembly of at most only a few micelles can currently be studied with these techniques. As such, coarse-graining (CG) approaches are often applied in this field, such as the MARTINI force field.11,18−28 Dissipative particle dynamics (DPD) is another CG approach that is becoming increasingly popular for the simulation of micelle behavior as it allows significantly longer length and time scales to be accessed at reduced computational cost (albeit with a potential loss in accuracy).29−31 Like other CG approaches, in the DPD method, groups of atoms or molecules are coarsegrained into particles (beads hereafter) that obey Newton’s Received: January 25, 2018 Published: March 23, 2018 2633

DOI: 10.1021/acs.jctc.8b00075 J. Chem. Theory Comput. 2018, 14, 2633−2643

Article

Journal of Chemical Theory and Computation

standard choice for reduced units in which the beads have unit mass, the system is governed by temperature kBT = 1, and the baseline cutoff distance for the short-range soft pairwise repulsions between solvent beads is set as rc = 1 (note that we allow deviations from this for nonsolvent beads, as described below). 2.1. Atom to Bead Mapping. In this study, we extend our recent work in which DPD parameters were determined by fitting to experimental partition coefficients (log P values), mutual solubilities, and liquid density data, for a range of small molecules.38 In this approach, DPD beads represent molecular fragments comprising n = 1−3 “heavy atoms” (i.e., carbon, oxygen, and nitrogen in this work), with the exception of water (H2O), which is treated supermolecularly. Water beads in the model are defined by a mapping number of Nm = 2, so that each water bead corresponds on average to two water molecules.42 Following well-established protocol,30 the density of water beads is set to ρr3c = 3 in reduced DPD units, where rc is the above-mentioned cutoff distance. The mapping ρNmvm ≡ 1, where vm ≈ 30 Å3 is the molecular volume of liquid water, then gives rc ≈ 5.65 Å in physical units. Alkane chains are constructed from connected (i.e., bonded) beads comprising (i) CH2CH2 groups of atoms and (ii) CH3, a terminal methyl group. To model the sulfate moiety in this work, we use a bead comprised of six heavy atoms, i.e., CH2OSO3−, which carries a negative charge of −1. In an attempt to improve transferability, we incorporate a carbon atom in the sulfate bead as this atom will be subject to polarization by the sulfate chemical group. The sodium counterion is modeled in a partially hydrated form, consisting of a Na atom and two water molecules. We denote this bead by [Na+]. It carries a charge of +1. To summarize, the CG mapping for alkyl sulfate surfactants in our model is [CH3][CH2CH2]m[CH2OSO3−][Na+], where the contents of the square brackets denote the different bead types and m = 2− 5 determines the alkyl chain length. 2.2. Nonbonded Interactions. The nonbonded interactions are defined by the usual DPD pairwise soft repulsions, 1 UijC = 2 Aij R ij(1 − rij/R ij)2 for rij ≤ Rij and UCij = 0 for rij > Rij. Here Aij is the repulsion amplitude of the conservative force, Rij the cutoff distance, and rij = |rj⃗ − ri⃗ | the separation between beads i and j located at ri⃗ and rj⃗ , respectively. In this study, we use the values of Aij and Rij for the uncharged species as set out in our earlier work.38 We note that there is an element of trade-off between Aij and Rij because DPD is an example of a “mean-field fluid”,43 and the thermodynamic properties are expected to be largely determined by AijRij3. On this basis, Rij ≡ rc = 1 is frequently assumed in the literature, and attempts to accommodate differences in bead “volumes” are confined to the Aij matrix. However, for our chosen coarse-grained mapping, beads contain different numbers of atoms and therefore contribute unequally to the total molar volume of the molecules under consideration. Hence, we choose to separately specify Aij and Rij and use the self-repulsion cutoffs Rii to capture the contribution of the molecular fragments to the overall molar volumes. To do this, we utilize the rules of Durchschlag and Zipper for individual atom contributions to molar volume44 and assign the value of Rii3 for different beads in proportion to the fragment molar volume, taking the molar volume of water as a reference (for which Rii3 = rc3 = 1). To deal with the cutoff between dissimilar bead types, we adopt a simple arithmetic “mixing

laws of motion. Contrary to the many CG approaches though, such as the MARTINI force field, DPD employs very soft conservative potentials, which facilitates the use of large time steps. Further details of the DPD method can be found in the literature29−32 and in the Supporting Information. DPD simulations have been applied to the study of micelle systems by a number of authors. For example, for worm-like micelles, Tang et al. calculated a number of micelle properties for model body-wash systems.11 Chunyang et al. studied the pH-responsive self-assembly of amphiphilic dendritic polymers into micelles of different morphologies depending on the simulated pH.33 However, only a small number of DPD studies have reported the CMC and Nagg properties of micelles. The majority of these concern the behavior of nonionic surfactants.34−38 To the best of our knowledge, there is currently only one study where explicit electrostatic interactions have been used to probe the CMC and Nagg of micelles of ionic surfactants: Mao et al. investigated the behavior of sodium dodecyl sulfate (SDS) and cetyltrimethylammonium bromide (CTAB) and were able to reproduce the experimental CMC to within a factor of 2−3 for both surfactants.39 In this pioneering work, Mao and co-workers went on to examine the morphology of micelles comprising both SDS and CTAB, observing the appearance of rod-like micelles with the addition of NaCl. The sphere-to-rod transition has also been explored through CG and DPD simulations by a handful of other authors.24,39,40 Recently we developed a new DPD parametrization scheme based on fitting parameters to thermodynamic quantities characteristic of liquids and liquid mixtures. Specifically, the fitting procedure adopted mutual solubilities, partition coefficients, and molar volumes (densities).38 In the current article, we utilize and extend this parametrization strategy and show that these parameters are transferable and sufficient to match experimental values and trends for the CMC and Nagg for sodium alkyl sulfate surfactants. These surfactants are representative of a broad class of anionic surfactants widely used in industry. Specifically we consider sodium hexyl, octyl, decyl, and dodecyl sulfates (denoted S2nS, where n = 3−6) and compare our results favorably to the CMCs and Nagg values reported in the literature. For S12S (i.e., SDS), we extend the study to observe the sphere-to-rod morphology transition and again compare favorably to experiment. These results further validate and support our overall DPD parametrization strategy, and we anticipate studying more complex behaviors of surfactants in the future. The remainder of the article is arranged as follows. In section 2, the adopted DPD model is presented. Section 3 gives an overview of the simulation and analysis methodology used to determine the CMC, Nagg, and sphere-to-rod transition. Results and corresponding discussion appear in section 4. Finally, in section 5, we give some thoughts on the accuracy of the models and directions for further work.

2. DEVELOPING THE DPD MODEL In this work, we adopt the DPD approach in which particles (DPD beads) interact via soft repulsions and with local pairwise dissipative and random forces, which together work as a thermostat.30 Rather than repeat here the details of what is now a fairly standard simulation method, we point the reader to the Supporting Information (SI), chapter 17 of the textbook by Frenkel and Smit,32 and the original DPD literature.29,30,41 For an up to date perspective on the DPD methodology, see Español and Warren.31 In our DPD simulations, we use a 2634


Article

Journal of Chemical Theory and Computation 1

modified to eliminate the divergence at overlap. The electrostatic pair potential between ions is

rule” R ij = 2 (R ii + R jj). The cutoff for the Na+ bead was assigned to be the same as that for water beads (Rij = 1). To complete the description, the cutoff for the dissipative and random forces was assigned equal to the maximum cutoff distance in the system (i.e., set as the maximum value of Rii), and the dissipative friction amplitude was set at γ = 4.5. In our previous communication, a combination of mutual solubilities, partition coefficients (log P), and liquid densities was used to determine the values for Aij.38 However, we concur with Tang et al., who point out that there has been no systematic approach described for the parametrization of charged beads in DPD to date.11 As such, we follow their pragmatic approach so that as Rij varies we set AijRij3 = 25 for all interactions involving two ionic beads or one ionic bead and water. This choice of scaling reflects the above-mentioned mean-field behavior of DPD when operated in the usual parameter regime of interest.31,43 For [Na+]−[alkyl] interactions, we adopt the same parameters between water and alkanes beads. Finally, for [CH2OSO3−]−[alkyl] interactions, we adopt the same headgroup to tail-group interactions as those used in our previous work.38 We shall comment though on how the calculated CMC values change in response to altering the interaction between the sulfate [CH2OSO3−] bead and the water beads, although we do not attempt a systematic parametrization for reasons discussed later. Table 1 lists the baseline values of Rij and Aij for all beads used in the present study.

UijE(rij) =

bead i

bead j

Aij

Rij

H2O CH2CH2 CH3 CH2OSO3− Na+ CH2CH2 CH3 CH2OSO3− Na+ CH3 CH2OSO3− Na+ Na+ Na+ CH2OSO3−

25.0 45.0 45.0 17.9 25.0 22.0 23.0 28.5 45.0 24.0 28.5 45.0 25.0 17.9 13.3

1.0000 1.0370 0.9775 1.1170 1.0000 1.0740 1.0145 1.1540 1.0370 0.9550 1.0945 0.9775 1.0000 1.1170 1.2340

4πrij

* [1 − (1 + β*rij)e−2β rij]

(1)

where rij is the ion separation, qi and qj are the ion charges (valencies), Γ = e2/(kBTϵ0ϵrrc) is a dimensionless electrostatic coupling parameter that includes the relative dielectric permittivity ϵr of the background (assumed uniform), and β* is the tunable Slater smearing parameter, which we set as β* = 0.929rc−1. Our implementation of Slater model electrostatics uses standard Ewald methods and is identical to that recently described by Vaiwala et al.47 We set the real-space Ewald cutoff as 3.0rc. Using ϵr = 78.3 for water at T = 298 K and our value of rc, we find Γ = 15.94. We caution that this is specific to our choice of mapping number (Nm = 2, which fixes rc), and Γ should be recalculated if a different choice is made. 2.3. Bonded Interactions. A simple harmonic potential 1 UijB = 2 k b(rij − r0)2 is used to represent bonds between connected DPD beads. Bond lengths are set according to the number of heavy atoms ni and nj in the bonded beads,38 as r0 = 0.1(ni + nj) − 0.01. For example, for the [CH2CH2]− [CH2CH2] bond, we set r0 = 0.39. We use a similar approach to set the length of the alkyl−sulfate bead interactions. Here, there are a total of eight heavy atoms involved in the bonded pair, so that r0 would be 0.79 following the above convention. We break with this though and specify r0 = 0.59 to reflect the fact that only six heavy atoms are actually bonded linearly (the other two being “branched” from the main chain). A single bond constant kb = 150 (DPD units) was adopted throughout. We explicitly introduce rigidity by including a harmonic angular potential between pairs of bonds. We here adopt the same three-body angular potential used by Smit and collabora1 tors,48,49 viz., UijkA = 2 ka(θijk − θ0)2 , where θijk is the angle between adjoining bonds. In the present work, we set θ0 = 180° and ka = 5 (DPD units).

Table 1. Repulsion Amplitudes (Aij) and Cutoff Distances (Rij) between All Bead Pairs in the Model H2O H2O H2O H2O H2O CH2CH2 CH2CH2 CH2CH2 CH2CH2 CH3 CH3 CH3 Na+ CH2OSO3− CH2OSO3−

Γqiqj

3. CALCULATING MICELLE PROPERTIES Where computer simulations have been used to determine CMC and Nagg, a common approach has been to simulate at concentrations well above the CMC of the surfactant in question. Here, the aggregate size distribution Ptot(N) typically shows two distinct regions: a monotonically decreasing submicellar region corresponding to monomers and submicellar aggregates (Pfree) and a peak at a much larger mean aggregation number corresponding to micelles (Pmicelle). An example is shown in Figure 1a. In principle, this allows us to distinguish between submicellar “free” surfactant at concentration CF and micellar surfactant at concentration CM, such that the total surfactant concentration CT = CF + CM. From this, we can also define the mole fraction of free and micellar surfactant as CF/ CT and CM/CT, respectively. The CMC can be computed from CF and CT, as described below, while the mean aggregation number Nagg is determined as the number weighted average of the micellar aggregate size distribution. The distinction between monomers and submicellar aggregates on the one hand and micelles on the other can be defined by a cutoff point in the distribution Ptot(N) such that the two regions are cleanly separated. This point is usually one of lowest probability between the two regions. For some surfactants and concentrations, assigning a cutoff is straightfor-

We do not attempt to match the compressibility of water in our model, and we note that we are not alone in regarding an exact match as unnecessarily restrictive.45 Our view, explained in more detail in Appendix B in Anderson et al.,38 is that for the choice of baseline density ρrc3 = 3 and for Aij ≳ O(10), the DPD fluid can be regarded as being thermodynamically incompressible on length scales ≳ rc. Fraaije et al. recently also expressed a similar opinion: “for practical mixing thermodynamics calculations the exact value of the compressibility is in fact completely irrelevant”.45 For the electrostatic interactions, we adopt the Slater-type charge smearing proposed by González-Melchor et al.46 in which the Coulomb potential between point charges is 2635


Article


reduction in accessible volume for the free surfactant. This “crowding” effect can be corrected by factoring in the volume (packing fraction) occupied by surfactant in the simulation box. For ionic surfactants, there is additionally a much stronger dependence of CF on CM above the CMC. This is largely a consequence of the incomplete association of counterions with micellar aggregates, such that adding more (neutral) surfactant drives the micellar equilibrium in the direction of incorporating more surfactant monomers into micelles (see the Appendix). In order to calculate the CMC for an ionic surfactant based on simulated CF, one must account for both the surfactant packing fraction correction and the effect of the micellar equilibrium. Starting from a charged pseudophase model (see the Appendix), Quina and co-workers,55 proposed a method to account for the micellar equilibrium that was further modified to account for the crowding effect.56,57 Calculated CF values can be used to determine the system CMC using the following expression ⎛ βC + (1 − β)C T ⎞ log C F + β log⎜ F ⎟ = (1 + β)log CMC (1 − ϕ) ⎝ ⎠ (2)

Here the factor 1/(1 − ϕ) accounts for the crowding effect on the unassociated counterions. In this, ϕ = VmCT is the surfactant packing fraction and Vm is the surfactant molar volume, which we calculate for each surfactant assuming a density of 1 g cm−3. The parameter β in eq 2 is the degree of counterion binding to the micelles. Equation 2 has been applied by a number of authors in the determination of the surfactant CMC from simulation.17,39,58 We note that eq 2 implies that the CMC is the total surfactant concentration at the point where the first micelles start appearing, where CF = CT. In the present study, we also introduce a second approach for determining the CMC of surfactants. In the second approach (method B hereafter), we define the CMC as the point at which half of the surfactant material is free (i.e., as monomers and submicellar aggregates) and half exists as micelles.36 In this method, the amount of free surfactant in the system is obtained by fitting an exponential decay to the probability of monomers, dimers, and trimers in the system, in the overall aggregate size distribution Ptot(N) (Figure 1c). This gives the free surfactant distribution, Pfree, for the system. We subtract Pfree from the overall distribution to obtain the micelle distribution, Pmicelle. The CMC is then defined as the point at which ∑∞ N=1 Pfree(N) = ∑ N∞= 1 P m i c e l l e (N). This is equivalent to asserting 1 C F = CM = 2 C T . Figure 1 shows how the total distribution can be separated into free surfactant and micellar aggregates for three concentrations (below, at, and above the CMC) for the S8S molecule. 3.1. General Simulation and Analysis Conditions. For all simulations, cubic boxes with periodic boundary conditions were used in the NVT ensemble. The box size, number of alkyl sulfate molecules, and number of water beads were chosen to obtain the desired concentrations while maintaining the pressure to within 2% of 23.7 (DPD units). This pressure, for which our parameters were originally developed, corresponds to that of pure water beads at a reduced density ρrc3 = 3.38 A DPD time step of Δt = 0.01 was adopted throughout, and trajectory data was collected every 20 DPD time units (equivalent to 2000 time steps). Electrostatics were solved using a smoothed-particle mesh Ewald (SPME) algorithm,59 and in all cases, the number of k⃗ vectors was set equal to the

Figure 1. Aggregation number distributions for S8S calculated at (a) 4.0, (b) 1.0, and (c) 0.25× the experimental CMC. The total distribution (Ptot, black line) can be separated into submicellar (Pfree, red dashed line) and micellar contributions (Pmicelle, blue dashed line). As the concentration increases above the CMC, the distinction between submicellar and micellar surfactant becomes clearer (a), but around the CMC, a significant overlap is observed (b).

ward as the two regions are separated by a clear minimum. However, frequently the choice of a cutoff becomes less obvious and may even depend on the total concentration CT.36 This can be particularly problematic for some surfactants when sampling at concentrations around the CMC (see Figure 1b). For systems where there is no clear separation between free surfactant and micelle distributions, small changes in the choice of cutoff potentially alter the calculated CMC values. This method of determining the CMC and Nagg (method A hereafter) is convenient as it makes it possible to calculate the CMC by simulating at relatively high surfactant concentrations, thereby avoiding the need for large simulations of dilute systems. However, this method must be used with caution. At its simplest, the CMC can be read off as the plateau value of CF as CT is increased above the CMC. This supposes that the added surfactant goes into building more micelles and the free surfactant concentration remains approximately constant, reflecting the approximate constancy of the surfactant chemical potential. However, it is known that CF does not necessarily remain constant as CT is increased above the CMC.17,50−54 For nonionic surfactants, there is a weak dependence due to a 2636


Article


Figure 2. Free monomer concentration CF versus the total surfactant concentration CT for (a) S6S, (b) S8S, (c) S10S, and (d) S12S. Black circles are CF calculated using method A. The red squares are the CMCs then obtained from eq 2. These are fit to a horizontal straight line (red dashed) to obtain a consensus CMC. Finally, the consensus CMC is used to produce the black solid lines from eq 2.

over 5 blocks of 5000 DPD time units following equilibration. The same procedure was followed for both methods A and B. When determining the CMC using method B, it is obviously not always possible to hit the targeted equality between free and micellar surfactant at a preselected sampling concentration. We therefore determine the point where CF = CM based on a linear interpolation between the nearest sampled concentrations. 3.3. Mean Aggregation Number Determination. The mean aggregation number Nagg was calculated at a fixed concentration of 7.5 wt % for all molecules. This corresponds to ∼26, 8.5, 2.8, and 1× the calculated CMC for S12S, S10S, S8S, and S6S, respectively. It is known from experiment that Nagg increases with total surfactant concentration CT. We therefore compare our calculated Nagg values to aggregation numbers calculated using the model for the concentration dependence of Nagg on CT reported by Ranganathan et al., which was determined from experimental observations.62 The value of Nagg was determined from simulation boxes that initially contained 500 surfactant monomers. The box size and the number of water beads were adjusted to achieve the desired concentration. The aggregation number distribution is slower to equilibrate than the CMC. Therefore, longer simulation times were adopted to ensure complete equilibration. Equilibrium in Nagg was defined as the point at which block averages (where a block is 5000 DPD time units) of the number of micelles fluctuate about a constant value within the error in the measurement. The time take to equilibrate Nagg was 500, 4500, 6700, and 12 400 DPD time units for the S6S, S8S, S10S, and S12S molecules, respectively. Data was collected over five blocks following equilibration, and results were monitored to check for slow growth in Nagg beyond the five blocks (each system was run for a minimum of 3 × 104 DPD time units). The Nagg that we report is Nagg = ⟨N2⟩/⟨N⟩, where ⟨Nz⟩ = z ∑∞ N=1 N Pmicelle(N), in other words, the number weighted average of the micellar component of the aggregate size distribution. In order to do this, we adopt the protocol of

simulation box size (to the nearest integer) to ensure that the relative error in calculating electrostatic energies was kept below 1%. All simulations were started from initial random configurations. DPD simulations were performed using the 60 DL_MESO simulation package. To identify surfactant aggregates, each simulation snapshot in the saved simulation trajectory was analyzed to identify connected clusters of surfactants. A surfactant was said to be in contact with a neighbor if the separation distance between any of the hydrophobic alkane beads was less than 1 DPD unit. All analysis of aggregate distributions, CMC, and Nagg was performed using the auxiliary UMMAP analysis package.61 3.2. CMC Determination. To determine the calculated CMC, we sampled each system at multiple concentrations, both below and above the experimentally reported CMC value. The sampled concentrations correspond approximately to 0.25, 0.5, 0.75, 1, 1.5, 2, 3, and 4× the experimental CMC. A number of simulation box sizes were used to ensure the presence of sufficient surfactant material. For the molecule with the highest experimental CMC (S6S), a box of 203 was used (24 000 beads). For each additional alkyl bead in the surfactant tail, the box was expanded by an additional 10 DPD units in each dimension. Consequently, the S12S systems were simulated in a box of size 503 (375 000 beads). By using such sizes, we ensure that there is always some nonmicellar surfactant that contributes to CF. Varying simulation times were also required to achieve satisfactory equilibration for each molecule studied. Equilibrium in CF was defined as the point at which variations in the block averaged CF (where a block is 5000 DPD time units) were smaller than the error in the reported value. On the whole, CF equilibrates relatively rapidly (compared to Nagg), although when sampling surfactants with a lower CMC at low concentrations, longer simulations are required. Thus CF equilibrates after 25 000 DPD time units (or 2.5 × 106 time steps) for the slowest equilibrating system (S12S) at concentrations 1.0−1.5× the experimental CMC. The calculated CMC values were then determined using data collected 2637


Article

Journal of Chemical Theory and Computation method B for the CMC calculation to decouple the distribution into monomer (and submicellar) and micellar contributions. 3.4. Identifying the Sphere-to-Rod Transition of S12S (SDS). Simulations were carried for S12S (i.e., SDS) for selected concentrations in the range of 30−100× the experimental CMC to identify the sphere-to-rod morphology transition point; the literature quoted value for this is ∼250 mM.63,64 Simulation box sizes were determined using the method described in the previous section for Nagg. Simulations were run for 6 × 104 DPD time units, and the final 4 × 104 DPD time units were used for data collection (we found that shorter equilibration times were sufficient to equilibrate at these higher concentrations). From these, we compute Nagg as defined above, the maximum micelle aggregation number Nmax, and the radius of gyration, Rg (defined below) for all micelle aggregates. Three repeats of each simulation were carried out to reduce noise in reported values. The radius of gyration Rg is calculated for each aggregate from the positions of hydrophobic beads that make up the core of the micelle, assuming that each DPD bead has the same mass. For a spherical micelle, Rg is expected to be directly proportional to its radius. Because the volume of the hydrophobic core is proportional to Nagg, in general, we expect spherical aggregates to be signaled by Rg ∝ N1/3 but a larger power law to be obtained for rod- and worm-like structures. We therefore identify the calculated sphere-to-rod transition as the concentration at which the observed relationship between N and Rg begins to deviate significantly from Rg ∝ N1/3.

Figure 3. Experimental versus calculated CMCs for n-alkyl sulfate surfactants for n = 6, 8, 10, and 12. Experimental values show an exponential dependence on alkyl length (black squares and solid line). Both methods adopted to calculate the CMC (method A: red line and circles; method B: blue line and crosses) also show an exponential dependence. Calculated values are generally in good agreement with experiment. Error bars are omitted in this plot to prevent clutter. Error estimates for the calculated CMCs are given in Table 2.

Table 2. Calculated CMCs and Nagg Values Compared to Experimental Dataa Nagg (7 wt %)

CMC/mM surfactant

expt.

method A

method B

expt.

calc.

S6S S8S S10S S12S (SDS)

420 130 33 8.2

210(5) 78(2) 22.5(5) 6.7(1)

360(9) 116(3) 34(2) 10(1)

17 36 56 87

16(1) 31(2) 39(4) 47(3)

a

Experimental CMCs and S6S experimental Nagg values are from Aniansson et al.,68 and S8S−S12S experimental Nagg values were determined according to Ranganathan et al.,62 to correspond to 7.5 wt %. This is equivalent to ∼260, 290, 320, and 360 mM for S12S−S6S, respectively. These concentrations are approximately 26, 8.5, 2.8, and 1× the calculated CMC resulting from method B. For calculated results, the estimated error is shown as the figure in brackets.

4. RESULTS AND DISCUSSION 4.1. Critical Micelle Concentrations. Figure 2 presents CF as a function of CT, calculated using method A. For all systems, we find that CF increases linearly with CT while the system is below the CMC and no micelles form. At higher CT values, CF declines again. The location of the peak in CF with respect to CT, i.e., separating the regions of increasing and decreasing trends in the CF value, corresponds well with the CMC calculated from eq 2. Note that we adopt a value of β = 0.7 in eq 2 for all systems. Using this value of β minimizes the dependence of the calculated CMC values on CT, for all alkyl chain lengths (see Figure 2 red points), to be consistent with the theory.55−57 We find that deviating from β = 0.7 creates minor changes (within error) in the predicted value of CMC. Sanders et al. have previously demonstrated that the CMC is insensitive to the precise value of β in the range of 0.4−0.8.17 We observe the same insensitivity in our results. We note that to achieve this level of quantitative agreement with theory it is crucial to ensure that the pressure remains within 2% of 23.7 (DPD units) as set in the original parametrization strategy.38 If this level of control is not maintained, then the resulting calculated CMCs do not remain constant with increasing CT. Higher pressures lead to an underestimation and lower pressures to an overestimation of the CMC. The values for the CMC determined using methods A and B are compared to experimental CMC values in Figure 3, and a detailed list of values obtained is given in Table 2. For both methods, we are able to reproduce quantitatively the experimental trend of decreasing CMC values with increasing alkyl chain length. The choice of method has an effect on the computed values. Method A results in lower CMC values than method B, and the latter has a better overall fit to experiment. As discussed in the Appendix, this difference is because the

CMC computed by method A is the point at which the first micelles appear (i.e. CF = CT), whereas the CMC in method B is the point at which half of the surfactant is in micelles (i.e., 1 C F = 2 C T ). As we show in the Appendix, this implies Δlog10 CMC = 0.254, which is in accord with the results in Figure 3. We note that the experimental CMC values are taken from Aniansson et al. and are from conductometry; therefore, the definition of the CMC in method B may be more appropriate. Both experimental and calculated CMCs decrease exponentially with alkyl chain length.65 The slope of the log10 CMC versus n best-fit lines for the calculated results (−0.252 and −0.257 for methods A and B, respectively) are in reasonable agreement with the slope describing the experimental data (−0.285). The experimental CMCs show a slightly stronger dependence on alkyl chain length than the calculated CMCs, which underpredict the values for the shorter surfactants. The fit that we achieve with our model to experimental CMC values is pleasing as we made no effort to tailor the sulfate bead interactions. Altering the [CH2OSO3−]−[H2O] interaction parameter to be more hydrophilic (Aij = 12 versus 17.9) results in a positive shift in calculated CMC values for both methods A and B. The shift is in the order of 25% (Δlog10 CMC ≈ 0.1) for each surfactant. However, changing this parameter in isolation does not improve the correlation of the gradients of the best-fit lines for the calculated CMCs compared to the experimental CMCs. A more systematic 2638


Article


only presents a weak decrease, which we attribute to the crowding effect mentioned in the main text, which is typically observed for nonionic surfactants. 4.2. Mean Aggregation Numbers. Experimental and calculated Nagg values are compared in Table 2. Experimental values for S8S−S12S were determined using eq 3 of Ranganathan et al.62 The experimental Nagg for S6S was taken from Aniansson.68 Our calculated Nagg values increase with increasing alkyl chain length, but, with the exception of S6S, they fall short of the experimental Nagg values by varying amounts. For S8S, the calculated Nagg is 14% too small. For S10S and S12S, the shortfall is 30 and 46%, respectively. In almost all simulation studies of alkyl sulfates, Nagg is underpredicted. The magnitude of the underprediction in the present study for S8S−S12S is similar to that in previous work.17,26,39 Underprediction has also been observed for nonionic surfactants.38 This behavior has frequently been attributed to incomplete equilibration because it takes much longer for the micelle size distribution to equilibrate than for the free surfactant to reach a steady-state value.17,26,68 The work of Fujimoto and others adds credence to this theory.69−71 However, great care was taken in the present study to avoid problems with poor equilibration by ensuring that the simulation times are sufficiently long (see previous sections) for the number of micelles present to be equilibrated well before data collection occurs. As such, we believe that our underpredictions of Nagg are not a consequence of poor equilibration. One possible explanation is the assumption of a uniform background dielectric.39 However, only the oily cores of the micelles present regions of low dielectric permittivity, and it is a basic principle of electrostatics that there is no electric field inside of a uniformly charged shell. Therefore, if there is an effect of the lowered dielectric permittivity in the micelle cores, it must be subtle. For example, Alargova and co-workers identified a dielectric effect entering in the Stern layer (i.e., headgroup region), which may be a contributory factor in the micelle size underestimation.72,73 It is also conceivable that our parametrization strategy and DPD model could be further optimized to reproduce better the experimental Nagg values because DPD simulations have been shown to be able to reproduce this when used as a fitting target.66 However, this cited study did not incorporate explicit electrostatic interactions; therefore, it is hard to draw conclusions. 4.3. Sphere-to-Rod Transition of SDS. In Figure 5a, we show the maximum (Nmax) and mean (Nagg) aggregation numbers for S12S (SDS) as a function of CT up to around 800 mM, in addition to representative snapshots of the micelle structures. At concentrations up to CT ≈ 250 mM, the increase in both Nagg and Nmax is linear with CT. In this concentration range, micelles are broadly spherical in shape. At each concentration, there is a distribution in individual Rg values corresponding to the ever-present shape fluctuations, on top of the aggregation number distribution. Figure 5b represents this distribution at three values of CT as box-whisker plots of Rg3/N, computed from the raw Rg data, binned by N. For example, Figure 5b (bottom plot) confirms that the aggregates are all roughly spherical at CT = 250 mM because Rg3/N is approximately constant, as discussed previously. Above CT ≈ 500 mM, we observe (Figure 5a) that Nmax grows above the linear regime observed at lower CT. This growth occurs as micelles lengthen significantly and start to

parametrization effort is needed to improve this correlation or to better understand the observed discrepancy. This is the first simulation study to our knowledge that samples CF both above and below the CMC. In doing this, we have been able to verify that eq 2 provides an extremely good fit to CF so that a reliable CMC can be calculated using data from simulations at concentrations well above the actual CMC (provided the pressure is maintained at the correct value). In particular, the CMCs computed from eq 2 do not change appreciably (i.e., withing error) when sampling CF at significantly higher values of CT. For example, from simulations at 30, 40, and 50× the experimental CMC (of S12S), we calculate the CMC from eq 2 to be 6.1 ± 1.0 mM, which is within error of the CMC calculated using lower CT values. This reduces the computational cost considerably in two ways: first, one can use smaller simulation boxes for large concentrations, and second, the equilibration time is reduced if the micelle density is increased. To the best of our knowledge, prior to our work, only two groups attempted to calculate the CMC of SDS using DPD. Mai et al. introduced an uncharged model in which explicit electrostatic interactions were not included (in fact, the sodium counterion was neglected entirely).66 With this approach, the authors calculated a CMC of 9 mM, which is, however, in excellent agreement with experiment. Later, Mao et al. were the first to propose a DPD model of SDS, which incorporated electrostatic interactions.39 They reported a CMC of 19 mM. Both groups of authors determined the CMC using approaches analogous to method A in our present work: Mao et al. used a cutoff of 5 monomers to define the free surfactant region; Mai et al. used 10. Using the models of Mai et al. and Mao et al., we recomputed the CMC values to be 8.5 ± 0.5 and 27 ± 2 mM, respectively. The CMC that we calculate for the model of Mao et al. is higher than that reported in the original work, which we attribute to a correction in the implementation of the Slater model electrostatics in the DL_MESO simulation engine.67 Figure 4 illustrates how CF varies as a function of CT for the two models in addition to our own results obtained using method A. Strikingly, but expectedly, the significant decrease in CF with increasing CT is seen only for the models that include electrostatic interactions. The uncharged model of Mai et al.

Figure 4. SDS: CF versus CT for DPD models of Mai et al.66 (red squares), Mao et al.39 (blue triangles), and the present work using method A (black circles). The CMCs are indicated with arrows. The dashed line is CF = CT, and the solid line is eq 2 with CMC = 6.7 mM and β = 0.7. 2639


Article


Figure 5. (a) Maximum (Nmax) and mean aggregation numbers (Nagg) from simulation as a function of total concentration CT for SDS. The concentration range spans ∼4−100× the experimental CMC. Snapshots of typical structures for the largest aggregates present are shown as insets. These structures are preserved over at least 5 × 104 DPD time units. (b) Radius of gyration of individual aggregates plotted as box-whisker plots of Rg3/N as a function of binned aggregation number N (bin size 10). Results are shown for three values of CT, as indicated. The color scheme is arbitrary.

become rod-like. We note that Nagg itself only begins to appreciably deviate from the linear growth regime at concentrations above ∼600 mM. This is because Nagg is dominated by the larger number of smaller micelles in the simulation box and is therefore slower to increase. Similarly, Figure 5b (top plot) demonstrates how at ∼580 mM the simulation box contains two different types of micelle structures: broadly spherical aggregates for N ≲ 50−100, similar to Figure 5b (bottom), and elongated structures for N ≳ 100 where Rg3/N increases roughly quadratically with N, indicating the formation of rod-like micelles. The range of CT ≈ 350−500 mM is where we start to observe the formation of elongated, nonspherical micelles. In the analysis of Rg shown in Figure 5b (middle plot), one can identify a region conforming to spherical micelles (R3g/N constant) and a small population of larger nonspherical micelles. We therefore estimate that the sphere-to-rod transition in our SDS model occurs at CT ≈ 350−500 mM. This is ∼1.4−2.0× the experimental value of 250 mM.63,64 This overestimation is consistent with our underestimation of Nagg for SDS. Cryo-TEM studies of a number of surfactants have indicated that above the sphere-to-rod transition micelles fall into two distinct categories; globular and cylindrical.74 Cylindrical micelles present as semiflexible rods with spherical end-caps that are larger in diameter than the main cylinder body. This dumbbell shape for rod micelles has been predicted by theory75 and observed in molecular dynamics simulations.76 In this work, we also observe the formation of these dumbbell-type structures at concentrations around and above the sphere-torod transition (see the structures in Figure 5a).

We explored two methods to calculate the CMC from simulation. In method A, we compute the CMC as the point at which the first micelles appear using the dependence of the free surfactant concentration on the total concentration, using the charged pseudophase model for ionic surfactants. In method B, we decouple the free and micellar surfactant aggregation number distributions and define the CMC as the point at which 50% of surfactant is micellar. With both methods, we get good agreement with the experimentally reported CMCs and can match the trend in Nagg versus alkyl chain length. We find that method B provides the best agreement with experiment. In the context of method A, we have also shown that above the CMC the free surfactant (i.e., not assigned to micellar aggregates) declines with increasing total surfactant concentration, in good agreement with the theoretical expression in eq 2 (charged pseudophase model). As well as providing a pleasing confirmation of our model, this allows us to calculate the CMC of ionic surfactants using simulations at total concentrations well above the actual CMC. In this approach though, we caution that the pressure should be maintained at a value to correspond to the parametrization conditions, either by adjusting the contents of the simulation box (as we do here) or by simulating under constant-pressure (NPT) conditions. If the pressure deviates systematically from the set point, this can lead to spurious results. As with all authors who have studied the self-assembly of charged surfactants using explicit electrostatics, we are unable to reproduce Nagg to a level that is comparable to experimental reported values (with the exception of S6S). Although there is some recent evidence to suggest that this poor equilibration may be the case for atomistic simulations,71 we do not believe that this is the problem with our simulations. The exact reason for disparity in aggregation number is unknown. Potentially, it is a consequence of the assumed uniform dielectric permittivity that is often assumed in electrostatic CG simulation methods.39 There is also the possibility that the DPD parameters employed in this work are suboptimal for reproducing aggregation numbers as Mai et al. were able to obtain good calculated values when the aggregation number was used explicitly as a parametrization target.66 It is difficult to compare our study to

5. CONCLUSIONS In this work, we used the DPD simulation method to calculate the critical micelle concentrations (CMCs) and mean aggregation numbers (Nagg) of a family of alkyl sulfate surfactants, with alkyl chain lengths in the range of 6−12. We have extended our recently published DPD parametrization strategy38 to include both charged sulfate beads and aqueous sodium ions. 2640



■

the work of Mai as in the latter case explicit electrostatics were not included in their model. Finally, we have estimated that the calculated sphere-to-rod transition for SDS (S12S) occurs in the region of 350−500 mM (which compares favorably to the experimentally observed 250 mM63,64). Above this concentration, we observe rod-like micelles that display a dumbbell-type structure in which the end-caps have higher curvature than the body of the micelle. This observation is also in agreement with experimental observations and theory.74,75

■

*E-mail: [email protected] (R.L.A.). *E-mail: andrea@ferrantescientific.com (A.S.F.). ORCID

Richard L. Anderson: 0000-0002-4320-8472 Funding

This work was in part supported by the STFC Hartree Centre’s Innovation: Return on Research programme, funded by the UK Department for Business, Energy & Industrial Strategy. Notes

The authors declare the following competing financial interest(s): Patrick B. Warren declares a substantive stock holding in Unilever PLC.

Charged Pseudophase Model

■

We summarize the essential elements of the charged pseudophase model that leads to eq 2.57,77,78 In this, the free surfactant CF and free counterion CI concentrations are constrained by what is, in effect, a solubility product, CFCβI ≤ K, where K is a constant and β plays the role of the stoichiometry. This can be derived as the N ≫ 1 limit of a quasi-chemical association model for micelle formation, supposing that N surfactant monomers and βN counterions combine to form one micelle. As mentioned, β (not to be confused with the Slater smearing parameter) is an important parameter that empirically accounts for electrostatic interactions in these models. From this perspective, the decline in CF for surfactant concentrations above the CMC is the analogue of the common ion effect for sparingly soluble salts. By simple mass balance and using the fact that the total counterion concentration should be equal to CT by charge neutrality, one has CI = βCF + (1 − β)CT. We insert this into the logarithm of the solubility product to find that above the CMC

ACKNOWLEDGMENTS R.L.A. and D.J.B. thank Michael Johnston and William Swope for many stimulating discussions that helped shape the work contained in this article. P.B.W. appreciatively recalls learning solution thermodynamics from Denver Hall many years ago.

■

REFERENCES

(1) Kwon, G. S.; Okano, T. Polymeric micelles as new drug carriers. Adv. Drug Delivery Rev. 1996, 21, 107−116. (2) Kedar, U.; Phutane, P.; Shidhaye, S.; Kadam, V. Advances in polymeric micelles for drug delivery and tumor targeting. Nanomedicine 2010, 6, 714−729. (3) Ahmad, Z.; Shah, A.; Siddiq, M.; Kraatz, H.-B. Polymeric micelles as drug delivery vehicles. RSC Adv. 2014, 4, 17028−17038. (4) Yokoyama, M. Polymeric micelles as drug carriers: their lights and shadows. J. Drug Target. 2014, 22, 576−583. (5) Biswas, S.; Kumari, P.; Lakhani, P. M.; Ghosh, B. Recent advances in polymeric micelles for anti-cancer drug delivery. Eur. J. Pharm. Sci. 2016, 83, 184−202. (6) Mandal, A.; Bisht, R.; Rupenthal, I. D.; Mitra, A. K. Polymeric micelles for ocular drug delivery: From structural frameworks to recent preclinical studies. J. Controlled Release 2017, 248, 96−116. (7) Cagel, M.; Tesan, F. C.; Bernabeu, E.; Salgueiro, M. J.; Zubillaga, M. B.; Moretton, M. A.; Chiappetta, D. A. Polymeric mixed micelles as nanomedicines: Achievements and perspectives. Eur. J. Pharm. Biopharm. 2017, 113, 211−228. (8) Zheng, R.; Liu, G.; Devlin, M.; Hux, K.; Jao, T.-c. Friction Reduction of Lubricant Base Oil by Micelles and Crosslinked Micelles of Block Copolymers. Tribol. Trans. 2009, 53, 97−107. (9) Spikes, H. Friction Modifier Additives. Tribol. Lett. 2015, 60, 5. (10) Skoglund, S.; Lowe, T. A.; Hedberg, J.; Blomberg, E.; Wallinder, I. O.; Wold, S.; Lundin, M. Effect of Laundry Surfactants on Surface Charge and Colloidal Stability of Silver Nanoparticles. Langmuir 2013, 29, 8882−8891. (11) Tang, X.; Zou, W.; Koenig, P. H.; McConaughy, S. D.; Weaver, M. R.; Eike, D. M.; Schmidt, M. J.; Larson, R. G. Multiscale Modeling of the Effects of Salt and Perfume Raw Materials on the Rheological Properties of Commercial Threadlike Micellar Solutions. J. Phys. Chem. B 2017, 121, 2468−2485. (12) Thompson, G. R. Absorption of fat-soluble vitamins and sterols. J. Clin. Pathol. 1971, s3-5, 85−89. (13) Chakraborty, T.; Chakraborty, I.; Ghosh, S. The methods of determination of critical micellar concentrations of the amphiphilic systems in aqueous medium. Arabian J. Chem. 2011, 4, 265−270. (14) Shelley, J.; Watanabe, K.; Klein, M. L. Simulation of a sodium dodecylsulfate micelle in aqueous solution. Int. J. Quantum Chem. 1990, 38, 103−117. (15) Bruce, C. D.; Berkowitz, M. L.; Perera, L.; Forbes, M. D. E. Molecular Dynamics Simulation of Sodium Dodecyl Sulfate Micelle in Water: Micellar Structural Characteristics and Counterion Distribution. J. Phys. Chem. B 2002, 106, 3788−3793.

(3)

We now note that the point where the first micelles appear has CF = CT (=CMC), and thus, log K = (1 + β) log CMC. This essentially recovers eq 2. The factor 1/(1 − ϕ) that additionally appears (in the free counterion concentration) was empirically introduced by Bales.57 In principle, this means that a similar factor should appear in the definition of the CMC, but if the CMC is small, this correction can be neglected. We can now also see more clearly the connection between method A and method B used in the main text. In method B, 1 the CMC is defined as the point where C F = 2 C T ; thus, in terms of eq 3, one has ⎛⎛ ⎛1 ⎞ 1 ⎞ ⎞ log⎜ C T⎟ + β log⎜⎜1 − β ⎟C T⎟ = (1 + β) log CMC ⎝2 ⎠ ⎝ ⎝ 2 ⎠ ⎠ (4)

This solves to CT = f(β) × CMC. For β = 0.7, for instance, f(β) ≈ 1.8. This gives rise to the prediction Δlog10 CMC ≈ log10(1.8) ≈ 0.254 used in the main text (see Figure 3).

■

AUTHOR INFORMATION

Corresponding Authors

APPENDIX

log C F + β log(βC F + (1 − β)C T) = log K

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00075. Further details of the DPD method (PDF) 2641


Article

Journal of Chemical Theory and Computation (16) Yoshii, N.; Okazaki, S. Free energy of water permeation into hydrophobic core of sodium dodecyl sulfate micelle by molecular dynamics calculation. J. Chem. Phys. 2007, 126, 096101. (17) Sanders, S. A.; Sammalkorpi, M.; Panagiotopoulos, A. Z. Atomistic Simulations of Micellization of Sodium Hexyl, Heptyl, Octyl, and Nonyl Sulfates. J. Phys. Chem. B 2012, 116, 2430−2437. (18) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (19) Wang, S.; Larson, R. G. Coarse-Grained Molecular Dynamics Simulation of Self-Assembly and Surface Adsorption of Ionic Surfactants Using an Implicit Water Model. Langmuir 2015, 31, 1262−1271. (20) Burov, S. V.; Obrezkov, N. P.; Vanin, A. A.; Piotrovskaya, E. M. Molecular dynamic simulation of micellar solutions: A coarse-grain model. Colloid J. 2008, 70, 1−5. (21) Wu, R.; Deng, M.; Kong, B.; Yang, X. Coarse-Grained Molecular Dynamics Simulation of Ammonium Surfactant Self-Assemblies: Micelles and Vesicles. J. Phys. Chem. B 2009, 113, 15010−15016. (22) Jalili, S.; Akhavan, M. A coarse-grained molecular dynamics simulation of a sodium dodecyl sulfate micelle in aqueous solution. Colloids Surf., A 2009, 352, 99−102. (23) Sanders, S. A.; Panagiotopoulos, A. Z. Micellization behavior of coarse grained surfactant models. J. Chem. Phys. 2010, 132, 114902. (24) Velinova, M.; Sengupta, D.; Tadjer, A. V.; Marrink, S.-J. Sphereto-Rod Transitions of Nonionic Surfactant Micelles in Aqueous Solution Modeled by Molecular Dynamics Simulations. Langmuir 2011, 27, 14071−14077. (25) Kraft, J. F.; Vestergaard, M.; Schiøtt, B.; Thøgersen, L. Modeling the Self-Assembly and Stability of DHPC Micelles Using Atomic Resolution and Coarse Grained MD Simulations. J. Chem. Theory Comput. 2012, 8, 1556−1569. (26) LeBard, D. N.; Levine, B. G.; Mertmann, P.; Barr, S. A.; Jusufi, A.; Sanders, S.; Klein, M. L.; Panagiotopoulos, A. Z. Self-assembly of coarse-grained ionic surfactants accelerated by graphics processing units. Soft Matter 2012, 8, 2385−2397. (27) Bennett, W. F. D.; Chen, A. W.; Donnini, S.; Groenhof, G.; Tieleman, D. P. Constant pH simulations with the coarse-grained MARTINI model - Application to oleic acid aggregates. Can. J. Chem. 2013, 91, 839−846. (28) Tang, X.; Koenig, P. H.; Larson, R. G. Molecular Dynamics Simulations of Sodium Dodecyl Sulfate Micelles in Water - The Effect of the Force Field. J. Phys. Chem. B 2014, 118, 3864−3880. (29) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155. (30) Groot, R. D.; Warren, P. B. Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 1997, 107, 4423−4435. (31) Español, P.; Warren, P. B. Perspective: Dissipative particle dynamics. J. Chem. Phys. 2017, 146, 150901. (32) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press, Inc.: Orlando, FL, 2001. (33) Yu, C.; Ma, L.; Li, K.; Li, S.; Liu, Y.; Liu, L.; Zhou, Y.; Yan, D. Computer Simulation Studies on the pH-Responsive Self-Assembly of Amphiphilic Carboxy-Terminated Polyester Dendrimers in Aqueous Solution. Langmuir 2017, 33, 388−399. (34) Vishnyakov, A.; Lee, M.-T.; Neimark, A. V. Prediction of the Critical Micelle Concentration of Nonionic Surfactants by Dissipative Particle Dynamics Simulations. J. Phys. Chem. Lett. 2013, 4, 797−802. (35) Lee, M.-T.; Vishnyakov, A.; Neimark, A. V. Calculations of Critical Micelle Concentration by Dissipative Particle Dynamics Simulations: The Role of Chain Rigidity. J. Phys. Chem. B 2013, 117, 10304−10310. (36) Johnston, M. A.; Swope, W. C.; Jordan, K. E.; Warren, P. B.; Noro, M. G.; Bray, D. J.; Anderson, R. L. Toward a Standard Protocol for Micelle Simulation. J. Phys. Chem. B 2016, 120, 6337−6351.

(37) Lee, M.-T.; Mao, R.; Vishnyakov, A.; Neimark, A. V. Parametrization of Chain Molecules in Dissipative Particle Dynamics. J. Phys. Chem. B 2016, 120, 4980−4991. (38) Anderson, R. L.; Bray, D. J.; Ferrante, A. S.; Noro, M. G.; Stott, I. P.; Warren, P. B. Dissipative particle dynamics: Systematic parametrization using water-octanol partition coefficients. J. Chem. Phys. 2017, 147, 094503. (39) Mao, R.; Lee, M.-T.; Vishnyakov, A.; Neimark, A. V. Modeling Aggregation of Ionic Surfactants Using a Smeared Charge Approximation in Dissipative Particle Dynamics Simulations. J. Phys. Chem. B 2015, 119, 11673−11683. (40) Sangwai, A. V.; Sureshkumar, R. Coarse-Grained Molecular Dynamics Simulations of the Sphere to Rod Transition in Surfactant Micelles. Langmuir 2011, 27, 6628−6638. (41) Español, P.; Warren, P. B. Statistical mechanics of dissipative particle dynamics. Europhys. Lett. 1995, 30, 191. (42) Groot, R.; Rabone, K. Mesoscopic Simulation of Cell Membrane Damage, Morphology Change and Rupture by Nonionic Surfactants. Biophys. J. 2001, 81, 725−736. (43) Louis, A. A.; Bolhuis, P. G.; Hansen, J. P. Mean-field fluid behavior of the Gaussian core model. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2000, 62, 7961−7972. (44) Durchschlag, H.; Zipper, P. Calculation of the partial volume of organic compounds and polymers. Prog. Colloid Polym. Sci. 1994, 94, 20−39. (45) Fraaije, J. G. E. M.; van Male, J.; Becherer, P.; Serral Gracià, R. Coarse-Grained Models for Automated Fragmentation and Parametrization of Molecular Databases. J. Chem. Inf. Model. 2016, 56, 2361−2377. (46) González-Melchor, M.; Mayoral, E.; Velázquez, M. E.; Alejandre, J. Electrostatic interactions in dissipative particle dynamics using the Ewald sums. J. Chem. Phys. 2006, 125, 224107. (47) Vaiwala, R.; Jadhav, S.; Thaokar, R. Electrostatic interactions in dissipative particle dynamicsEwald-like formalism, error analysis, and pressure computation. J. Chem. Phys. 2017, 146, 124904. (48) Venturoli, M.; Smit, B. Simulating the self-assembly of model membranes. PhysChemComm 1999, 2, 45−49. (49) Venturoli, M.; Maddalenasperotto, M. M.; Kranenburg, M.; Smit, B. Mesoscopic models of biological membranes. Phys. Rep. 2006, 437, 1−54. (50) Floriano, M. A.; Caponetti, E.; Panagiotopoulos, A. Z. Micellization in Model Surfactant Systems. Langmuir 1999, 15, 3143−3151. (51) Cheong, D. W.; Panagiotopoulos, A. Z. Monte Carlo Simulations of Micellization in Model Ionic Surfactants: Application to Sodium Dodecyl Sulfate. Langmuir 2006, 22, 4076−4083. (52) Kahlweit, M.; Teubner, M. On the kinetics of micellization in aqueous solutions. Adv. Colloid Interface Sci. 1980, 13, 1−64. (53) Johnson, I.; Olofsson, G.; Jonsson, B. Micelle formation of ionic amphiphiles. Thermochemical test of a thermodyanamic model. J. Chem. Soc., Faraday Trans. 1 1987, 83, 3331−3344. (54) Polacek, R.; Kaatze, U. Monomer Exchange Kinetics, Radial Diffusion, and Hydrocarbon Chain Isomerization of Sodium Dodecylsulfate Micelles in Water. J. Phys. Chem. B 2007, 111, 1625−1631. (55) Quina, F. H.; Nassar, P. M.; Bonilha, J. B. S.; Bales, B. L. Growth of Sodium Dodecyl Sulfate Micelles with Detergent Concentration. J. Phys. Chem. 1995, 99, 17028−17031. (56) Soldi, V.; Keiper, J.; Romsted, L. S.; Cuccovia, I. M.; Chaimovich, H. Arenediazonium Salts: New Probes of the Interfacial Compositions of Association Colloids. 6. Relationships between Interfacial Counterion and Water Concentrations and Surfactant Headgroup Size, Sphere-to-Rod Transitions, and Chemical Reactivity in Cationic Micelles. Langmuir 2000, 16, 59−71. (57) Bales, B. L. A. Definition of the Degree of Ionization of a Micelle Based on Its Aggregation Number. J. Phys. Chem. B 2001, 105, 6798−6804. 2642


Article

Journal of Chemical Theory and Computation (58) Jusufi, A.; Panagiotopoulos, A. Z. Explicit- and Implicit-Solvent Simulations of Micellization in Surfactant Solutions. Langmuir 2015, 31, 3283−3292. (59) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577−8593. (60) Seaton, M. A.; Anderson, R. L.; Metz, S.; Smith, W. DL_MESO: highly scalable mesoscale simulations. Mol. Simul. 2013, 39, 796−821. (61) Bray, D. J. UMMAP. https://www.scd.stfc.ac.uk/Pages/ UMMAP.aspx (2017). (62) Ranganathan, R.; Tran, L.; Bales, B. L. Surfactant- and SaltInduced Growth of Normal Sodium Alkyl Sulfate Micelles Well above Their Critical Micelle Concentrations. J. Phys. Chem. B 2000, 104, 2260−2264. (63) Reiss-Husson, F.; Luzzati, V. The Structure of the Micellar Solutions of Some Amphiphilic Compounds in Pure Water as Determined by Absolute Small-Angle X-Ray Scattering Techniques. J. Phys. Chem. 1964, 68, 3504−3511. (64) Christov, N. C.; Denkov, N. D.; Kralchevsky, P. A.; Ananthapadmanabhan, K. P.; Lips, A. Synergistic Sphere-to-Rod Micelle Transition in Mixed Solutions of Sodium Dodecyl Sulfate and Cocoamidopropyl Betaine. Langmuir 2004, 20, 565−571. (65) Klevens, H. B. Critical micelle concentrations as determined by refraction. J. Phys. Colloid Chem. 1948, 52, 130−148. (66) Mai, Z.; Couallier, E.; Rakib, M.; Rousseau, B. Parameterization of a mesoscopic model for the self-assembly of linear sodium alkyl sulfates. J. Chem. Phys. 2014, 140, 204902. (67) Seaton, M. A. DL_MESO INFOMAIL mail shots. https://www. scd.stfc.ac.uk/Pages/DL_MESO-infomail.aspx (2016). (68) Aniansson, E. A. G.; Wall, S. N.; Almgren, M.; Hoffmann, H.; Kielmann, I.; Ulbricht, W.; Zana, R.; Lang, J.; Tondre, C. Theory of the kinetics of micellar equilibria and quantitative interpretation of chemical relaxation studies of micellar solutions of ionic surfactants. J. Phys. Chem. 1976, 80, 905−922. (69) Kawada, S.; Komori, M.; Fujimoto, K.; Yoshii, N.; Okazaki, S. Molecular dynamics study of the formation mechanisms of ionic SDS and nonionic C12E8 micelles and n-dodecane droplets. Chem. Phys. Lett. 2016, 646, 36−40. (70) Kawada, S.; Fujimoto, K.; Yoshii, N.; Okazaki, S. Molecular dynamics study of the potential of mean force of SDS aggregates. J. Chem. Phys. 2017, 147, 084903. (71) Fujimoto, K.; Kubo, Y.; Kawada, S.; Yoshii, N.; Okazaki, S. Molecular dynamics study of the aggregation rate for zwitterionic dodecyldimethylamine oxide and cationic dodecyltrimethylammonium chloride micelles. Mol. Simul. 2017, 43, 1331−1337. (72) Alargova, R. G.; Danov, K. D.; Kralchevsky, P. A.; Broze, G.; Mehreteab, A. Growth of Giant Rodlike Micelles of Ionic Surfactant in the Presence of Al3+ Counterions. Langmuir 1998, 14, 4036−4049. (73) Alargova, R.; Ivanova, V.; Kralchevsky, P.; Mehreteab, A.; Broze, G. Growth of rod-like micelles in anionic surfactant solutions in the presence of Ca2+ counterions. Colloids Surf., A 1998, 142, 201−218. (74) Zana, R. Giant micelles; Taylor & Francis, 2007. (75) Majhi, P. R.; Dubin, P. L.; Feng, X.; Guo, X.; Leermakers, F. A. M.; Tribet, C. Coexistence of Spheres and Rods in Micellar Solution of Dodecyldimethylamine Oxide. J. Phys. Chem. B 2004, 108, 5980− 5988. (76) Poghosyan, A. H.; Arsenyan, L. H.; Shahinyan, A. A. Shape of Long Chain Alkyl Sulfonate Micelle upon Salt Addition: a Molecular Dynamics Study. J. Surfactants Deterg. 2015, 18, 755−760. (77) Shinoda, K.; Hutchinson, E. Pseudo-phase separation model for thermodynamic calculations on micellar solutions. J. Phys. Chem. 1962, 66, 577−582. (78) Hall, D. G. Thermodynamics of solutions of polyelectrolytes, ionic surfactants and other charged colloidal systems. J. Chem. Soc., Faraday Trans. 1 1981, 77, 1121−1156.

2643


Micelle Formation in Alkyl Sulfate Surfactants Using Dissipative

Recommend Documents