Micelle Formation in Alkyl Sulfate Surfactants ... - ACS Publications

Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Condensed Matter, Interfaces, and Materials

Micelle Formation in Alkyl Sulphate Surfactants Using Dissipative Particle Dynamics Richard L. Anderson, David J. Bray, Annalaura Del Regno, Michael A. Seaton, Andrea S. Ferrante, and Patrick B Warren J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00075 • Publication Date (Web): 23 Mar 2018 Downloaded from http://pubs.acs.org on March 25, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Micelle Formation in Alkyl Sulphate 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 E-mail: [email protected]; [email protected] Abstract We use dissipative particle dynamics (DPD) to study micelle formation in alkyl sulphate surfactants, with alkyl chain lengths ranging six to twelve carbon atoms. We extend our recent DPD force-field [J. Chem. Phys., 2017, 147, 094503] to include a charged sulphate chemical group, and aqueous sodium ions. With this model we achieve good agreement with the experimental 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 pseudo-phase 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 the surfactant is classed as micellar and half as monomers and sub-micellar 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 SDS micelles, and determine that it corresponds to SDS concentrations in the region 300–500 mM.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction Used as detergents, foaming and wetting agents, 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 behaviour. 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. 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 structures, 14,15 and free energy of micelle formation. 16 Recent simulations on the micellisation of sodium alkyl sulphates (hexyl, heptyl, octyl, and nonyl sulphates) have calculated equilibrium properties such as CMC and mean aggregation number (Nagg ) although the method 2


Page 2 of 42

Page 3 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 which is becoming increasingly popular for the simulation of micelle behaviour 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 coarse-grained into particles (beads hereafter) which obey Newton’s 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 literature, 29–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 an 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 behaviour of non-ionic 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 behaviour of sodium dodecyl sulphate (SDS) and cetyl trimethyl ammonium bromide (CTAB) and were able to reproduce the experimental CMC to within a factor 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

3



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 parameterisation 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 parameterisation strategy, and show that these parameters are transferable and sufficient to match experimental values and trends for CMC and Nagg for sodium alkyl sulphate 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 sulphates ( denoted S2nS where n = 3, 4, 5, 6) and compare our results favourably to the CMCs and Nagg values reported in the literature. For S12S (i. e. sodium dodecyl sulphate - SDS) we extend the study to observe the sphere-to-rod morphology transition, and again compare favourably to experiment. These results further validate and support our overall DPD parameterisation strategy, and we anticipate studying more complex behaviours 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.

Developing the DPD model In this work we adopt the Dissipative Particle Dynamics (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

4


Page 4 of 42

Page 5 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Repulsion amplitudes (Aij ) and cut-off distances (Rij ) between all bead pairs in the model. bead i

bead j

Aij Rij

H2 O H2 O H2 O H2 O H2 O CH2 CH2 CH2 CH2 CH2 CH2 CH2 CH2 CH3 CH3 CH3 Na+ CH2 OSO− 3 CH2 OSO− 3

H2 O CH2 CH2 CH3 CH2 OSO− 3 Na+ CH2 CH2 CH3 CH2 OSO− 3 Na+ CH3 CH2 OSO− 3 Na+ Na+ Na+ CH2 OSO− 3

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

(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 Espa˜ nol and Warren. 31 In our DPD simulations we use a standard choice for reduced units in which the beads have unit mass, the system is governed by temperature kB T = 1, and the baseline cut-off distance for the short-range soft pairwise repulsions between solvent beads is set as rc = 1 (note that we allow deviations from this for non-solvent beads, as described below).

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 (H2 O) which is treated super-molecularly. 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 5



water beads is set to ρrc3 = 3 in reduced DPD units, where rc is the above-mentioned cut-off distance. The mapping ρNm vm ≡ 1, where vm ≈ 30 ˚ A3 is the molecular volume of liquid water, then gives rc ≈ 5.65 ˚ A in physical units. Alkane chains are constructed from connected (i. e. bonded) beads comprising (i) CH2 CH2 groups of atoms and (ii) CH3 , a terminal methyl group. To model the sulphate moiety in this work we use a bead comprised of 6 heavy atoms, i. e. CH2 OSO− 3 , which carries a negative charge of −1. In an attempt to improve transferability, we incorporate a carbon atom in the sulphate bead as this atom will be subject to polarisation by the sulphate chemical group. The sodium counterion is modelled in a partially-hydrated form, consisting of an Na atom and two water molecules. We denote this bead by [Na+ ]. It carries a charge of +1. To summarise therefore, the CG mapping for alkyl sulphate surfactants in our model is [CH3 ][CH2 CH2 ]m [CH2 OSO3 − ][Na+ ], where the contents of the square brackets denote the different bead types, and m = 2, 3, 4, 5 determines the alkyl chain length. As mentioned, for dynamical purposes all beads are assumed to have unit mass in reduced units.

Non-bonded interactions The non-bonded interactions are defined by the usual DPD pairwise soft repulsions, UijC = 1 A R (1 2 ij ij

− rij /Rij )2 for rij ≤ Rij and UijC = 0 for rij > Rij . Here Aij is the repulsion

amplitude of the conservative force, Rij the cut-off 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 since DPD is an example of a ‘mean-field fluid’, 43 and the thermodynamic properties are expected to be 3 largely determined by Aij Rij . 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 6


Page 6 of 42

Page 7 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


under consideration. Hence we choose to separately specify Aij and Rij , and use the selfrepulsion cut-offs 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 volume, 44 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 cut-off between dissimilar bead types we adopt a simple arithmetic ‘mixing rule’ Rij = 12 (Rii + Rjj ). The cut-off for the Na+ bead was assigned to be the same as for water beads (Rij = 1). To complete the description, the cut-off for the dissipative and random forces was assigned equal to the maximum cut-off 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 were 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 parameterisation of charged beads in DPD to date. 11 As such, we follow their pragmatic 3 = 25 for all interactions involving two ionic approach so that as Rij varies we set Aij Rij

beads or one ionic bead and water. This choice of scaling reflects the above-mentioned meanfield behaviour 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 [CH2 OSO− 3 ] − [alkyl] interactions we adopt the same head-group to tail-group interactions as 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 sulphate [CH2 OSO− 3] bead and the water beads, although we do not attempt a systematic parameterisation for reasons discussed later. Table 1 lists the baseline values of Rij and Aij for all beads used in the present study. 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, ex-

7



Page 8 of 42

plained 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 González-Melchor et al. 46 in which the Coulomb potential between point charges is modified to eliminate the divergence at overlap. The electrostatic pair potential between ions is

UijE (rij ) =

Γqi qj ∗ [1 − (1 + β ∗ rij )e−2β rij ] , 4πrij

(1)

where rij is the ion separation, qi and qj are the ion charges (valencies), Γ = e2 /(kB T 0 r rc ) is a dimensionless electrostatic coupling parameter which includes the relative dielectric permittivity r of the background (assumed uniform), and β ∗ is the tuneable Slater smearing parameter which we set as β ∗ = 0.929 rc−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 cut-off as 3.0 rc . 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.

Bonded interactions A simple harmonic potential UijB = 12 kb (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 [CH2 CH2 ]-[CH2 CH2 ] bond, we set r0 = 0.39. We use a similar approach to set the length of the alkyl-sulphate bead interactions. Here there are a total of 8 heavy atoms involved in the bonded pair, so that r0 would be 0.79 following the above convention. We break with this though, and

8


Page 9 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 A potential used by Smit and collaborators, 48,49 viz. Uijk = 21 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).

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 sub-micellar region corresponding to monomers and sub-micellar aggregates (Pfree ), and a peak at much larger mean aggregation number corresponding to micelles (Pmicelle ). An example is shown in Fig. 1a. In principle this allows us to distinguish between sub-micellar ‘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, whilst the mean aggregation number Nagg is determined as the number weighted average of the micellar aggregate size distribution. The distinction between monomers and sub-micellar aggregates on the one hand, and micelles on the other, can be defined by a cut-off point in the distribution Ptot (N ) such that the two regions are cleanly separated. This point is usually one of lowest (or at least low) probability between the two regions. For some surfactants and concentrations, assigning a cut-off is straightforward as the two regions are separated by a clear minimum. However,

9



0.05

(a)

P(N)

0.04

Ptotal Pfree Pmicelle

0.03 0.02 0.01 0.00 0

20

40

60

80

100

120

140

15

20

25

30

35

N

0.20

(b)

P(N)

0.15 0.10 0.05 0.00 0

0.60

10

5

N

(c)

0.50

P(N)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.40 0.30 0.20 0.10 0.00 0

2

4

6

N

8

10

12

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 sub-micellar (Pfree , red dashed line) and micellar contributions (Pmicelle , blue dashed line). As the concentration increases above the CMC the distinction between sub-micellar and micellar surfactant becomes clearer (a), but around the CMC a significant overlap is observed (b).

10


Page 10 of 42

Page 11 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


frequently the choice of a cut-off becomes less obvious and may even depend on the total concentration CT . 36 This can be particularly problematic for some surfactants around the CMC (see Fig. 1b). For systems where there is no clear separation between free surfactant and micelle distributions, small changes in the choice of cut-off potentially alter the calculated CMC values. This method of determining 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 non-ionic surfactants, there is a weak dependence due to a 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 Appendix). In order to calculate the CMC for an ionic surfactant based on simulated CF , one must account both for the surfactant packing fraction correction, and the effect of the micellar equilibrium. Starting from a charged psuedo-phase model (see Appendix), Quina and co-workers, 55 proposed a method to account for the micellar equilibrium that was further modified to account for crowding effect. 56,57 Calculated CF values can be used to determine

11



Page 12 of 42

the system CMC using the following expression

log CF + β log

βC + (1 − β)C F T = (1 + β) log CMC. (1 − φ)

(2)

Here the factor 1/(1 − φ) accounts for the crowding effect on the unassociated counterions. In this φ = Vm CT 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. Eq. (2) has been applied by a number of authors in the determination of surfactant CMC from simulation. 17,39,58 We note that Eq. (2) implies 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 sub-micellar 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 ) (Fig. 1c). This gives the free surfactant distribution, Pf ree 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 P∞ P∞ 1 N =1 Pfree (N ) = N =1 Pmicelle (N ). This is equivalent to asserting CF = CM = 2 CT . 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.

General simulation & analysis conditions For all simulations, cubic boxes with periodic boundary conditions were used, in the NVT ensemble. The box size, number of alkyl sulphate molecules, and number of water beads, were chosen to obtain the desired concentrations whilst maintaining the pressure to within

12


Page 13 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 timeunits (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 simulation box size (to the nearest integer) to ensure 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 dl meso simulation package. 60 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 neighbour 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 were performed using the auxiliary ummap analysis package. 61

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 non-micellar surfactant which 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 5 000 DPD time-units) were smaller than the error in the reported value. On the whole CF 13



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 over 5 blocks of 5 000 DPD time units following equilibration. The same procedure was followed for both methods A & 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 pre-selected sampling concentration. We therefore determine the point where CF = CM based on a linear interpolation between the nearest sampled concentrations.

Mean aggregation number determination The mean aggregation number Nagg was calculated at a fixed concentration of 7.5 wt % for all molecules, which 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 which initially contained 500 surfactant monomers. The box size and the number of water beads was 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 5 000 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 14


Page 14 of 42

Page 15 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1000

1000

(a) S6S

(b) S8S

100

100 100

10 10

1000

100

1000

100 100

(c) S10S

(d) S12S

10

10 10

1 1

100

10

100

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). Nagg beyond the five blocks (each system was run for a minimum of 3 ×104 DPD time-units). P z The Nagg we report is Nagg = hN 2 i/hN i, where hN z i = ∞ 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 method B for the CMC calculation to decouple the distribution into monomer (and sub-micellar) and micellar contributions.

15



Identifying the sphere-to-rod transition of S12S (SDS) Simulations were carried for S12S (i. e. SDS) for selected concentrations in the range 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 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 each DPD bead has the same mass. For a spherical micelle, Rg is expected to be directly proportional to its radius. Since the volume of the hydrophobic core is proportional to Nagg , in general we expect spherical aggregates to be signaled by Rg ∝ N 1/3 , but a larger power law to be obtained for rodand 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 ∝ N 1/3 .

Results and Discussion 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 whilst 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 CF value, corresponds

16


Page 16 of 42

Page 17 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 Fig. 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 0.4–0.8. 17 We observe the same insensitivity in our results. We note that to achieve this level of quantitative agreement with the theory it is crucial to ensure that the pressure remains within 2% of 23.7 (DPD units) as set in the original parameterisation 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 & B are compared to experimental CMC values in Fig. 3 and a detailed list of values obtained are 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 for 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 the surfactant is in micelles (i. e. CF = 21 CT ). As we show in the Appendix, this implies ∆ log10 CMC = 0.254, which is in accord with the results in Fig. 3. We note that the experimental CMC values are taken from Aniansson et al. and are from conductometry, so the definition of the CMC in method B may be the 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 de-

17



3.0

expt method A method B

2.5

log10(CMC/mM)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.0 1.5

1.0 0.5

6

8

n

10

12

Figure 3: Experimental versus calculated critical micelle concentrations (CMCs) for n-alkyl sulphate surfactants for n = 6, 8, 10, 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. scribing the experimental data (−0.285). The experimental CMCs show a slightly stronger dependence on alkyl chain length than the calculated CMCs which underpredicts the values for the shorter surfactants. The fit we achieve with our model to experimental CMC values is pleasing as we made no effort to tailor the sulphate bead interactions. Altering the [CH2 OSO− 3 ] – [H2 O] interaction parameter to be more hydrophilic (Aij =12 versus Aij =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 improved the correlation of the gradients of the best fit lines for the calculated CMCs compared to the experimental CMCs. A more systematic parameterisation 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 the Eq. (2) provides an extremely good fit to CF so that a reliable CMC can be calculated using data from simulations at concentrations

18


Page 18 of 42

Page 19 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 have 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 cut-off of five monomers to define the free surfactant region; Mai et al. used ten. Using the models of Mai et al. and Mao et al. we re-computed the CMC values to be 8.5 ± 0.5 mM and 27 ± 2 mM respectively. The CMC we calculate for the model of Mao et al. is higher than 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 Fig. 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 only seen for the models which include electrostatic interactions. The uncharged model of Mai et al. only presents a weak decrease which we attribute to the crowding effect mentioned in the main text, which is typically observed for non-ionic surfactants.

19



100 Mai et al. (2014) Mao et al. (2015) this work (method A)

10

1 1

10

100

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 arrowed. The dashed line is CF = CT and the solid line is Eq. (2) with CMC = 6.7 mM and β = 0.7.

Mean aggregation numbers Experimental and calculated Nagg 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 sulphates, Nagg is underpredicted. The magnitude of the underprediction in the present study for S8S–S12S is similar to previous work. 17,26,39 Underprediction has also been observed for non-ionic surfactants. 38 This behaviour has frequently been attributed incomplete equilibration since 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 occurred. As such we believe our underprediction

20


Page 20 of 42

Page 21 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Calculated CMCs and Nagg values compared to experimental data. Experimental CMCs and S6S experimental Nagg are from Aniansson et al., 68 and S8S–S12S experimental Nagg 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 to S6S, respectively. These concentrations are approximately (26, 8.5, 2.8 and 1 × calculated CMC resulting from method B) For calculated results, the estimated error is shown as the figure in brackets.

Surfactant

S6S S8S S10S S12S (SDS)

expt.

CMC / mM method A method B

420 130 33 8.2

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

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

Nagg (7 wt %) expt. calc. 17 36 56 87

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

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 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 coworkers identified a dielectric effect entering in the Stern layer (i. e. head group region) which may be a contributory factor in the micelle size underestimation. 72,73 It is also conceivable that our parameterisation strategy and DPD model could be further optimised to reproduce better the experimental Nagg , since 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 so it is hard to draw conclusions.

Sphere-to-rod transition of SDS In Fig. 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

21



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. Fig. 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 Fig. 5b (bottom plot) confirms the aggregates are all roughly spherical at CT = 250 mM, since Rg3 /N is approximately constant as discussed previously. Above CT ≈ 500 mM we observe (Fig. 5a) that Nmax grows above the linear regime observed at lower CT . This growth occurs as micelles lengthen significantly and start to 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 Fig. 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 Fig. 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 CT ≈ 350–500 mM is where we start to observe the formation of elongated, non-spherical micelles. In the analysis of Rg shown in Fig. 5b (middle plot), one can identify a region conforming to spherical micelles (Rg3 /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 semi-flexible 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

22


Page 22 of 42

Page 23 of 42

1000 800

N

(a)

(b)

Nmax Nagg

600

Rg3/N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


400 200 0 0

200

400

600

800

CT /mM

N

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. theory, 75 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-to-rod transition (see structures in Fig. 5a).

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 sulphate surfactants, with alkyl chain lengths in the range 6–12. We have extended our recently published DPD parameterisation strategy, 38 to include both charged sulphate beads and aqueous sodium ions.

23



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 psuedophase 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 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 pseudo-phase 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 parameterisation 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 this poor equilibration may be the case for atomistic simulations, 71 we do not believe 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 sub-optimal for reproducing aggregation numbers as Mai et al. were able to obtain good calculated values when aggregation number was used explicitly as a parameterisation target. 66 It is difficult

24


Page 24 of 42

Page 25 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to compare our study to 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 sodium dodecyl sulphate (S12S) occurs in the region of 350–500 mM (which compares favourably to the experimentally observed 250 mM. 63,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

Acknowledgement This work was in part supported by the STFC Hartree Centre’s Innovation: Return on Research programme, funded by the Department for Business, Energy & Industrial Strategy. RLA and DJB thank Michael Johnston and William Swope for many stimulating discussions that helped shape the work contained in this article. PBW appreciatively recalls learning solution thermodynamics from Denver Hall many years ago.

Appendix: charged pseudo-phase model We summarise the essential elements of the charged pseudo-phase model which 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, CF CI β ≤ 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 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 which 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 25



Page 26 of 42

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,

log CF + β log(βCF + (1 − β)CT ) = log K .

(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 the CMC is defined as the point where CF = 21 CT , thus in terms of Eq. (3) one has

log( 12 CT ) + β log((1 − 12 β)CT ) = (1 + β) log CMC .

(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 Fig. 3).

References (1) Kwon, G. S.; Okano, T. Polymeric micelles as new drug carriers. Adv. Drug Deliv. Rev. 1996, 21, 107 – 116. (2) Kedar, U.; Phutane, P.; Shidhaye, S.; Kadam, V. Advances in polymeric micelles for drug delivery and tumor targeting. Nanomed. Nanotech. Biol. Med. 2010, 6, 714 – 729.

26


Page 27 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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. Control. 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.; chi Jao, T. 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. Trans. 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.

27



(12) Thompson, G. R. Absorption of fat-soluble vitamins and sterols. J. Clin. Pathol. 1971, 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. Arab. 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. (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 SelfAssembly 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 Journal 2008, 70, 1–5. 28


Page 28 of 42

Page 29 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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. Sphere-to-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) Drew Bennett, W. F.; 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. 29



(30) Groot, P. B., Robert D.Warren Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 1997, 107, 4423–4435. (31) Espa˜ nol, 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, USA, 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.

30


Page 30 of 42

Page 31 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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) Espa˜ nol, 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 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 Gracià, R. Coarse-Grained Models for Automated Fragmentation and Parametrization of Molecular Databases. J. Chem. Inf. Model. 2016, 56, 2361–2377. (46) González-Melchor, M.; Mayoral, E.; Velá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. Phys. Chem. Comm. 1999, 2, 45–49. 31



(49) Venturoli, M.; Sperotto, 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. 32


Page 32 of 42

Page 33 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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. DLMESO: highly scalable mesoscale simulations. Mol. Simul. 2013, 39, 796–821. (61) Bray, D. J. https://www.scd.stfc.ac.uk/Pages/UMMAP.aspx 2017, (62) Ranganathan, R.; Tran, L.; Bales, B. L. Surfactant- and Salt-Induced 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. 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 33



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.

34


Page 34 of 42

Page 35 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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.

35


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Probability


Aggregation number

36


Page 36 of 42

Page Journal 37 of 42of Chemical Theory and Computation 0.05

(a)

1 0.04

P(N)

2 3 0.03 4 5 0.02 6 7 8 0.01 9 10 0.00 0 11 12 13 0.20 14 15 16 0.15 17 18 19 0.10 20 21 22 0.05 23 24 25 0.000 26 27 28 29 0.60 30 31 0.50 32 0.40 33 34 0.30 35 36 0.20 37 38 0.10 39 0.00 0 40 41

Ptotal Pfree Pmicelle

20

40

60

80

100

120

140

15

20

25

30

35

N

P(N)

(b)

10

5

N

P(N)

(c)


2

4

6

N

8

10

12

1000

1 2 3 4 5 6 100 7 8 9 10 11 12 13 100 14 15 16 17 18 19 20 10 21 22 23 24

(a) S6S

1000 Journal of Chemical Theory and Computation (b) S8S

Page 38 of 42

100

100

1000

10 10

100

1000

100

(c) S10S

(d) S12S

10

ACS Paragon Plus Environment 10

100

1 1

10

100

3.0

Page 39 Journal of 42 of Chemical Theory and Computation

expt method A method B

2.5 1 2 3 2.0 4 5 1.5 6 7 8 1.0 9 10 11 0.5


6

8

10

12


1000

N

800

(b)

Nmax Nagg

600

Rg3/N

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

(a)

Page 40 of 42

400 200 0 0

200

400


600

CT /mM

800

N

Page 41 of 42 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


100 Mai et al. (2014) Mao et al. (2015) this work (method A)

10

1 1

10


100


ToC Image 79x39mm (300 x 300 DPI)


Page 42 of 42

Micelle Formation in Alkyl Sulfate Surfactants ... - ACS Publications

Recommend Documents