Langevin Dynamics Simulations of Cationic Surfactants in Aqueous

N. Chennamsetty , H. Bock , M. Lísal , J. K. Brennan ... Alan R. Katritzky , Liliana M. Pacureanu , Svetoslav H. Slavov , Dimitar A. Dobchev , Dinesh...
1 downloads 0 Views 569KB Size
Download PDF
Langmuir 2004, 20, 2017-2025

2017

Langevin Dynamics Simulations of Cationic Surfactants in Aqueous Solutions Using Potentials of Mean Force Hiroyuki Shinto,* Shintaro Morisada, Minoru Miyahara, and Ko Higashitani Department of Chemical Engineering, Kyoto University, Nishikyo-ku, Kyoto 615-8510, Japan Received October 8, 2003. In Final Form: December 17, 2003 A Langevin dynamics (LD) technique is proposed for the simulation of surfactant molecules in aqueous solutions, where no water molecules of the solvent are explicitly treated, but the effects are incorporated using both the self-diffusion coefficients of the solutes and the potentials of mean force between them in water. The present LD technique is employed to simulate (i) the self-assembly of n-decyltrimethylammonium chloride surfactants in water and (ii) the micelles of different sizes in water and in a 0.5 M NaCl solution. In the first simulation starting from the disaggregate configuration, the surfactant aggregates grew by addition of monomers and small aggregates to them rather than by fusion of large aggregates into larger ones within 12 ns: this behavior corresponds to the early stage of surfactant self-assembly and indicates that the system did not completely reach the equilibrium states for the period of several nanoseconds. In the second simulation starting from the spherical micelles, the small micelle of 20 surfactants was unstable and less spherical in water compared with the larger micelles of 30 and 40 surfactants, but it exhibited a stable spherical shape in the NaCl solution where no significant difference among these micelles was found in the shape and the stability. The results are in fair agreement with those from the atomistic model molecular-dynamics simulations and the experimental measurements, demonstrating that our LD technique can capture the characteristics of ionic surfactants in aqueous solutions despite the implicit treatment of water molecules of the solvent.

1. Introduction Surfactant molecules have both hydrophilic and hydrophobic parts: the former has a strong affinity for water, while the latter has a poor affinity. In aqueous solutions, the surfactants associate into a variety of secondary structures such as spherical and cylindrical micelles, vesicles, and bilayers depending on their primary (or molecular) structures. If the surfactant concentration becomes high enough, the secondary structures often develop long-range order to become liquid crystalline.1,2 Various surfactants with the self-assembling nature have been utilized in the wide range of industrial and material processes; the choice and use of the surfactants are still more or less dependent on a rule of thumb, because the phenomena of surfactant self-assembly remain to be explored from a molecular point of view. Molecular dynamics (MD) and Monte Carlo (MC) simulations are very useful for investigating the microscopic features of matter and material.3,4 To investigate the molecular-level features of the surfactant solutions, several researchers have implemented the MD and MC simulations of surfactant systems.5-19 However, the atomistic model simulations of the self-assembly of * Author to whom correspondence should be addressed. Phone: +81-75-383-2672. Fax: +81-75-383-2652. E-mail: shinto@ cheme.kyoto-u.ac.jp. (1) Holmberg, K.; Jo¨nsson, B.; Kronberg, B.; Lindman, B. Surfactants and Polymers in Aqueous Solution, 2nd ed.; Wiley: Chichester, U.K., 2003. (2) Larson, R. G. The Structure and Rheology of Complex Fluids; Oxford University Press: New York, 1999. (3) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1987. (4) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2002. (5) Klein, M. L. J. Chem. Soc., Faraday Trans. 1992, 88, 1701. (6) Smit, B. In Computer Simulation in Chemical Physics; Allen, M. P., Tildesley, D. J., Eds.; NATO Advanced Study Institute Series C; Kluwer: Dordrecht, The Netherlands, 1993; Vol. 397, Chapter 12. (7) Pastor, R. W. Curr. Opin. Struct. Biol. 1994, 4, 486. (8) Esselink, K.; Hilbers, P. A. J.; Karaborni, S.; Siepmann, J. I.; Smit, B. Mol. Simul. 1995, 14, 259.

surfactants are still beyond the power of computers to realize because long-time simulations of the large-scale systems are required for studying the self-assembly of surfactants. A common alternative is the use of coarsegrained models to mimic the oil/water/surfactant systems. The coarse-grained model simulations, however, provide only the qualitative results, which are sufficiently suggestive but are not quantitatively comparable with the experimental results. Now pay attention to the fact that, in aqueous solutions of surfactants, the number of solvent molecules is extremely larger than that of solutes such as surfactants and salt ions. The second alternative is then to use Langevin dynamics (LD) simulations, where the effects of the solvent are considered via the frictional and random forces on the solutes and the effective forces between them.3 This implicit treatment of the solvent would realize the large-scale simulations, which are crucial for studying the self-assembly of surfactants. When the concentration of the solutes is not very high, the frictional and random forces are generated using the self-diffusion coefficient (SDC) of the solutes in the solvent at infinite dilution, and the effective interactions between the solutes are plausibly (9) Karaborni, S.; Smit, B. Curr. Opin. Colloid Interface Sci. 1996, 1, 411. (10) Larson, R. G. Curr. Opin. Colloid Interface Sci. 1997, 2, 361. (11) Jakobsson, E. Trends Biochem. Sci. 1997, 22, 339. (12) Tarek, M.; Bandyopadhyay, S.; Klein, M. L. J. Mol. Liq. 1998, 78, 1. (13) Bandyopadhyay, S.; Tarek, M.; Klein, M. L. Curr. Opin. Colloid Interface Sci. 1998, 3, 242. (14) Shinto, H.; Tsuji, S.; Miyahara, M.; Higashitani, K. Langmuir 1999, 15, 578. (15) Shelley, J. C.; Shelley, M. Y. Curr. Opin. Colloid Interface Sci. 2000, 5, 101. (16) Feller, S. E. Curr. Opin. Colloid Interface Sci. 2000, 5, 217. (17) Forrest, L. R.; Sansom, M. S. P. Curr. Opin. Struct. Biol. 2000, 10, 174. (18) Schmid, F. In Computational Methods in Surface and Colloid Science; Boro´wko, M., Ed.; Surfactant Science Series 89; Marcel Dekker: New York, 2000; Chapter 13. (19) Saiz, L.; Bandyopadhyay, S.; Klein, M. L. Biosci. Rep. 2002, 22, 151.

10.1021/la035874o CCC: $27.50 © 2004 American Chemical Society Published on Web 02/05/2004

2018

Langmuir, Vol. 20, No. 5, 2004

Shinto et al.

represented by the potential of mean force (PMF) at infinite dilution. The SDCs and the PMFs can be quantitatively computed using atomistic model MD or MC simulations. Such LD simulations with the SDCs and the PMFs will provide at least the semiquantitative results relevant to the experimental measurements of the corresponding systems. Similar LD simulations with the more simplified models of surfactants were already reported,20,21 but the resultant data are of course qualitative. Recently, the third alternative was proposed, which is a combination of the MC method applied only to the surfactant molecules and the reference-interaction-site-model theory for incorporating the solvent effects;22 however, this hybrid simulation technique has been applied only to the coarse-grained models. In this paper, we describe the LD simulation technique for surfactant solutions, where the SDCs of the solutes and the PMFs between them in water are used for incorporating the solvent effects. Using the present LD method, the aggregation of cationic surfactants at 0.1 M is investigated. Also, a micelle of the cationic surfactants in aqueous solutions is examined, where the micelles of different sizes and the solutions of 0 and 0.5 M NaCl are employed. Figure 1. PMFs for solute pairs in ambient water at infinite dilution. The solid lines are the guides for the eyes.

2. Simulation Methods 2.1. Langevin Equations of Motion. Let us first consider N solute particles in a solvent.3 The Langevin equations of motion for a solute particle i ()1, ..., N) are represented by

dri(t) ) vi(t) dt mi

∂W(r1, ..., rN ) dvi(t) )- miγivi(t) + Ri(t) dt ∂ri

(1a) (1b)

where mi is the mass of particle i and ri, vi, and γi are the position, velocity, and friction coefficient, respectively. In the right-hand side of eq 1b, the first term represents the interaction forces of the particle with other particles and the potential W(r1, ..., rN) is modeled by PMFs, which include the averaged effects of the solvent. The second term describes the frictional drag force on the particle due to the solvent, and this frictional force is proportional to the velocity vi(t). The third term is the force on the particle due to the random fluctuations caused by the interactions with the solvent molecules. The friction coefficient γi is independent of the time and the positions. The stochastic force Ri(t) is assumed to be a stationary Gaussian variable with zero mean and to be uncorrelated with prior velocities and the interparicle forces. These conditions require that Ri(t) be coupled with γi via a fluctuation-dissipation relation

〈RiR(t) Rjβ(t′)〉 ) 2miγikBTδijδRβδ(t - t′)

(2a)

〈RiR(t)〉 ) 0

(2b)

where 〈‚‚‚〉 denotes the ensemble average of the function enclosed. The subscripts R and β denote Cartesian coordinates (R, β ) x, y, z), i and j represent particle labels, (20) (a) von Gottberg, F. K.; Smith, K. A.; Hatton, T. A. J. Chem. Phys. 1997, 106, 9850. (b) von Gottberg, F. K.; Smith, K. A.; Hatton, T. A. J. Chem. Phys. 1998, 108, 2232. (21) Noguchi, H.; Takasu, M. Phys. Rev. E 2001, 64, 041913. (22) (a) Kinoshita, M.; Sugai, Y. Chem. Phys. Lett. 1999, 313, 685. (b) Kinoshita, M.; Sugai, Y. J. Comput. Chem. 2002, 23, 1445.

δ is the Kronecker delta or the Dirac delta function, kB is the Boltzmann constant, and T is the temperature. According to the Einstein relation, γi is given by

γi )

kBT miDi∞

(3)

where Di∞ is the SDC of particle i in the solvent at infinite dilution. 2.2. Implicit Solvent Models. In the present study, we consider n-decyltrimethylammonium chloride (CH3(CH2)9-N(CH3)3+Cl- or C10TAC), sodium chloride, and water as the cationic surfactant, salt, and solvent, respectively. We need the PMFs for every pair of three solutes, the surfactant, Na+, and Cl- in water at infinite dilution as well as the SDCs of the solutes. These PMFs and SDCs can be computed in principle using MD or MC simulations with the atomistic models. However, it is fatally difficult to compute the PMFs for the solute pairs including the surfactant(s) because the surfactant is a flexible molecule with many degrees of freedom. Therefore, we employed the following technique. The hydrophilic and the hydrophobic groups of the surfactant were approximated by a tetramethylammonium ion [(CH3)4N+ or TMA+] and a chain of methane molecules (CH4 or Me), respectively. The PMFs Wij(r) for 10 different binary combinations with TMA+, Me, Na+, and Cl- and the SDCs Di∞ of the four solutes in ambient water at infinite dilution were computed using MD simulations,23 and the results are summarized in Figure 1 and Table 1. The resultant PMFs and SDCs were used in eqs 1-3 for each interaction site, where the interactions were given by assuming the pair additivity of the PMFs. The surfactant of CH3(CH2)9-N(CH3)3+ was modeled as (Me)10-TMA+, that is, t10h+, where the hydrophobic and the hydrophilic sites of the surfactant are referred to as t and h+, (23) (a) Shinto, H.; Morisada, S.; Miyahara, M.; Higashitani, K. J. Chem. Eng. Jpn. 2003, 36, 57. (b) Shinto, H.; Morisada, S.; Miyahara, M.; Higashitani, K. Manuscript submitted for publication. (c) Shinto, H.; Morisada, S.; Miyahara, M.; Higashitani, K. Manuscript in preparation.

LD Simulations of Cationic Surfactants

Langmuir, Vol. 20, No. 5, 2004 2019

Table 1. SDCs Di∞ of Solutes in Ambient Water at Infinite Dilution Di∞ [10-9 m2 s-1]

Di∞ [10-9 m2 s-1]

solute

MDa

exptb

solute

MDa

exptb

CH4 (t) TMA+ (h+)

1.94 1.15

c 1.20

ClNa+

1.85 1.26

2.03 1.33

a Results from MD simulations (see ref 23c). b Experimental values (see ref 24). c Not available.

Table 2. Molecular Geometry of t10h+ Surfactant bond length [nm] t-t t-h+ t-t-t t-t-h+

bond angle [deg]

0.152 0.278 109.5 126.4 Table 3. Systems of LD Simulations

side concentration initial number of solutes length Csurf Csalt configuration [M] of surfactants system t10h+ Na+ Cl- L [nm] [M] A B.1 B.2 B.3 C.1 C.2 C.3

125 20 30 40 20 30 40

100 150 200

125 12.76 20 6.925 30 7.927 40 8.725 120 6.925 180 7.927 240 8.725

0.1 0.1 0.1 0.1 0.1 0.1 0.1

0.5 0.5 0.5

disaggregated spherical spherical spherical spherical spherical spherical

respectively. For intramolecular interactions, the bond lengths and the bond angles were fixed at constant values, as listed in Table 2, and only the PMFs between the sites separated by at least three bonds were included. It should be noted that the PMFs for separations beyond the range of the computed PMFs (i.e., r > 0.8 or 1.2 nm) were represented by

Wij(r) )

1 QiQj 4πr0 r

(4)

where Qi is the point charge of solute i, 0 is the permittivity of the vacuum, and r is the relative permittivity of water; the experimental value of r ) 77.6 at T ) 300 K was used.25 The details of the PMFs and the SDCs can be found elsewhere.23 2.3. Simulation Systems. The systems considered are summarized in Table 3: (A) a solution of the t10h+Clsurfactants in water at Csurf ) 0.1 M; (B) a micelle composed of different numbers of the surfactants in water at Csurf ) 0.1 M; and (C) the same as system B but in a NaCl solution of Csalt ) 0.5 M at Csurf ) 0.1 M. The initial configuration of the surfactants for system A is different from those for systems B and C, as will be described in sections 2.4.1 and 2.4.2, respectively. The concentrations of the surfactants for systems A-C are the same as Csurf ) 0.1 M, which is higher than the critical micelle concentration (cmc) of C10TAC in ambient water, 0.050.065 M.26,27 It is noted that the average size of the aggregates, Nav, for C10TAC is estimated as Nav ) 30, according to an extrapolation of the experimental results for CnTAC: Nav ) 44, 62, and 84 for n ) 12,28 14,29 and 16,30 respectively. (24) Lide, D. R. CRC Handbook of Chemistry and Physics, 80th ed.; CRC Press: Boca Raton, 1999. (25) The Chemical Society of Japan. Kagaku-binran Kiso-hen II, 4th ed.; Maruzen: Tokyo, Japan, 1993. (26) Bacaloglu, R.; Blasko´, A.; Bunton, C. A.; Cerichelli, G.; Ortega, F. J. Phys. Chem. 1990, 94, 5062. (27) Klevens, H. B. J. Am. Oil Chem. Soc. 1953, 30, 74. (28) Ozeki, S.; Ikeda, S. Bull. Chem. Soc. Jpn. 1981, 54, 552.

2.4. Simulation Details. We employed a cubic box with periodic boundary conditions, which includes the finite numbers of the solutes and has the side length L as listed in Table 3. The origin of the coordinates was taken at the center of the simulation boxes. The Langevin equations of motion (eq 1) for the sites of t10h+ surfactants, the Clcounterions, and the Na+ co-ions were integrated using the third-order algorithm proposed by van Gunsteren and Berendsen.31 The bond lengths and the bond angles of the surfactant molecules were fixed by the SHAKE method.32 The Coulomb contributions of eq 4 to the PMFs were handled by the Ewald sum technique,3 while the rest of the contributions were considered only within the range of the computed PMFs (i.e., r e 0.8 or 1.2 nm). The temperature was set at T ) 300 K, at which the PMFs and the SDCs used were evaluated.23 The time step of ∆t ) 1 fs was used unless specified otherwise. The configuration was stored every 100 fs for subsequent analysis. We prepared two types of the initial configuration for the LD simulations of systems A-C, as follows. 2.4.1. Surfactant Self-Assembling (System A). Selfassembly of surfactants in salt-free water was simulated using 125 t10h+Cl- surfactants at Csurf ) 0.1 M. Note that the number of the surfactants employed is about four times as large as Nav ) 30 for C10TAC. Their initial configuration was prepared as follows: (i) the center of mass of a fully extended t10h+ surfactant was placed at one of the 5 × 5 × 5 lattice points, and the surfactant was then randomly rotated around the lattice point; (ii) likewise, the rest of the surfactants were generated; and (iii) the Cl- counterions were randomly positioned in the simulation box such that they never overlapped with each other and the surfactant sites. This configuration was relaxed for 50 ps by scaling the velocities every 5 fs. Starting from the final configuration where the time was taken as t ) 0 (see part a of Figure 2), the system was allowed to evolve for 12 ns. 2.4.2. Surfactant Micelles of Different Sizes (Systems B and C). We simulated a micelle composed of 20, 30, and 40 t10h+Cl- surfactants in water at Csurf ) 0.1 M using systems B.1-B.3, respectively. Notice that these numbers of the surfactants employed are around Nav ) 30 for C10TAC. To investigate the influence of salt on the micelles, we also employed systems C.1-C.3, which are the same as systems B.1-B.3, respectively, but in a 0.5 M NaCl solution other than in water. For the initial configuration of a micelle, the fully extended t10h+ surfactants were arranged in a spherical shape following Tieleman et al.,33 where the surfactants radiated from the center of the simulation box with the tail ends (i.e., t10 sites) and the center being separated by 0.73 nm. The Cl- counterions and the Na+ co-ions were randomly located around the micelle such that they never overlapped with each other. This configuration was carefully relaxed by velocity scaling at every 2 fs with ∆t ) 0.2 fs for 30 ps; during the first 10 ps, the t10 site of each surfactant at rt10 was constrained toward the box center by an external potential

k Wext ) |rt10|2 2

(5)

with the force constant k, which was initially equal to 50 (29) Imae, T.; Ikeda, S. J. Phys. Chem. 1986, 90, 5216. (30) Reiss-Husson, F.; Luzzati, V. J. Phys. Chem. 1964, 68, 3504. (31) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Simul. 1988, 1, 173. (32) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (33) Tieleman, D. P.; van der Spoel, D.; Berendsen, H. J. C. J. Phys. Chem. B 2000, 104, 6380.

2020

Langmuir, Vol. 20, No. 5, 2004

Shinto et al.

Figure 2. Snapshots of 125 t10h+ surfactants and 125 Cl- counterions in water at Csurf ) 0.1 M (system A), at time t ) 0, 1, 3, 6, 9, and 12 ns. Green, red, and blue spheres denote t, h+, and Cl-, respectively. The simulation cell is depicted by the lines.

LD Simulations of Cationic Surfactants

Langmuir, Vol. 20, No. 5, 2004 2021

Figure 4. Normalized coordinate(s) x/L for the center(s) of mass of the instantaneously largest aggregate(s) of the surfactants in system A as a function of time t. The consecutive coordinates of four growing aggregates are denoted by AG1, AG2, AG3, and AG4. See also text.

Figure 3. Size distributions of surfactant aggregates in system A for different periods of 2 ns until time t ) 12 ns. The average concentration of the surfactants in aggregates of size N is denoted by XN. Each distribution was obtained by averaging over 2000 configurations taken every 1 ps for the period of 2 ns.

N/m and then decreased by 0.1 N/m every 20 fs to 0 after 10 ps. Thereafter, no velocity scaling was applied and ∆t ) 1 fs was used. Starting from the final configuration where the time was taken as t ) 0 (see parts a-c of Figure 6), the system was allowed to evolve over 10 ns. 3. Results and Discussion In the two following subsections, we investigate the selfassembling of surfactants in water (system A) and the surfactant micelles of different sizes in water (system B) and in a NaCl solution (system C), as summarized in Table 3. 3.1. Surfactant Self-Assembling. Figure 2 displays the snapshots of system A, where the surfactants and the counterions at t ) 0 are disaggregated, as expected from the procedure in section 2.4.1. Figure 3 shows the size distributions of the surfactant aggregates for different periods of 2 ns, where XN denotes the average concentration of the surfactants in aggregates of size N. Here, two surfactants were considered to belong to the same surfactant aggregate, when the minimum site-to-site separation between their tails was less than 0.48 nm at which the PMF for the CH4-CH4 pair equals 0, as shown in part b of Figure 1. The different aggregates were identified by the standard clustering procedure.34 It should be noted that the size distributions in Figure 3 are rather noisy as a result of the lack of sufficient sampling. As observed in part b of Figure 2 and part a of Figure 3, at the early stage of t ) 0-2 ns, most of the surfactants exist as monomers. At t ) 2-4 ns, surfactants spontaneously associate into small aggregates of N ) 2-7, as in part c of Figure 2 and part b of Figure 3. The aggregates then become larger with time, and by the end of the 12 ns simulation, they appear to be spherical micelles, which point the hydrophilic sites outward and the hydrophobic (34) Stoddard, S. D. J. Comput. Phys. 1978, 27, 291.

chains inward. Because of the continuous growth of the aggregates, as in Figure 3, system A did not reach the equilibrium states within the 12-ns simulation. On the other hand, Figure 3 shows that the concentration of monomers, X1, decreases with time until 10 ns and then reaches a constant value of about 0.01 M. This value is considered as a rough estimate of the cmc for the t10h+Clsurfactant, which is of the same order as the experimental values of the cmc for the corresponding real C10TAC surfactant, 0.05-0.065 M.26,27 Now we track the large aggregates of the surfactants to investigate their growth processes. Figure 4 shows the x coordinate(s) for the center(s) of mass of the largest aggregate(s) as a function of time t. At t < 0.7 ns, the coordinates are categorized into five groups that stay around the lattice points (i.e., x/L ) 0, (0.2, (0.4) for the initial configuration of the surfactants, signifying that all the surfactants exist as monomers. The coordinates are scattered during t ) 0.7-1.7 ns, indicating that some of the surfactants assemble into dimers and trimers (see part a of Figure 3) at several different places. Then, the consecutive coordinates of four growing aggregates appear, as denoted in Figure 4 by AG1, AG2, AG3, and AG4. Note that AG2 and AG4 exhibit similar coordinates in the x direction but have dissimilar coordinates in the other two directions (data not shown); hence, AG2 and AG4 are distinct aggregates. In contrast, two divided strings of the x coordinates for t ) 1.7-7.1 ns and t ) 11.5-12 ns in Figure 4 belong to the same aggregate, AG1. The aggregates of AG1, AG2, and AG3 grow to be the first three largest aggregates at 12 ns. The sizes of the three aggregates as a function of time, NAG1(t), NAG2(t), and NAG3(t), are displayed in Figure 5, where the maximum size of the aggregates, Nmax(t), is also plotted. The aggregate sizes often increase by 1 or 2, but they occasionally become larger by 3 or 4 at t ) 11.5, 10-10.7, and 7.1 ns for AG1, AG2, and AG3, respectively. By the occasional addition of a trimer or a tetramer, one of these aggregates exceeds the others in size to give Nmax. It is noted that the vertical bands in Figure 5 represent the fluctuation of the surfactant(s) approaching the aggregates. At 12 ns, NAG1, NAG2, and NAG3 equal 12, 10, and 9, respectively. These values are less than 1/2 of the average size Nav ) 30 of the C10TAC aggregates, indicating that system A was not completely equilibrated because of such a short duration of 12 ns. During the 12-ns simulation starting from the disaggregate configuration, the surfactant aggregates grew by addition of monomers and small aggregates of N ) 2-4 to them rather than by fusion of large aggregates into larger ones. This is because the

2022

Langmuir, Vol. 20, No. 5, 2004

Figure 5. Sizes of the aggregates of AG1, AG2, and AG3 (black lines) in system A as a function of time t: (a) NAG1(t); (b) NAG2(t); and (c) NAG3(t). The gray lines represent the variation of the maximum size, Nmax(t). See also Figure 4.

diffusion of a large aggregate is much slower than that of a monomer and a small aggregate. It should be noted that the time scale of several nanoseconds is too short to examine the important behavior of micellar solutions, such as the surfactant exchange and the micellar formation and breakdown that may occur on the time scale of milliseconds. Several research groups have reported MD simulations of 20-60 surfactants starting from random configurations in water. Maillet et al.35 employed 48 (>Nav ) 25)36 C9TAC surfactants in 2997 water molecules at Csurf ) 0.23 M (>cmc ) 0.12 M)37 and found two spherical micelles of 15 and 17 surfactants as well as a small number of monomers after 0.9 ns. This result is similar to our result just mentioned. On the other hand, Marrink et al.39 simulated 54 (≈Nav ) 50-60)40 zwitterionic surfactants of dodecylphosphocholine in 22 496 water molecules at Csurf ) 0.12 M (.cmc ) 0.001 M).40 They observed that all the surfactants spontaneously aggregated into a single spherical micelle after about 12 ns and then the micelle exhibited its equilibrium structure after another 5 ns. This successful observation of a large micelle would be achieved by the use of the extremely high concentration compared to the cmc: Csurf/cmc ≈ 120, whereas Csurf/cmc ≈ 2 for system A in the present study. 3.2. Surfactant Micelles of Different Sizes. The surfactant solution of Csurf ) 0.1 M starting from the disaggregate configuration (system A) did not completely reach the equilibrium states within the 12-ns simulation, as mentioned in section 3.1. Hence, we investigate the surfactant micelles of different sizes of N ) 20, 30, and 40 in water at Csurf ) 0.1 M using systems B.1-B.3, respectively. These numbers of the t10h+Cl- surfactants employed are around Nav ) 30 for C10TAC. Also, we examine the influence of salt on the micelles using systems (35) Maillet, J.-B.; Lachet, V.; Coveney, P. V. Phys. Chem. Chem. Phys. 1999, 1, 5277. (36) The average size of the aggregates for C9TAC was estimated in the same way as that for C10TAC. See text in section 2.3. (37) The cmc for C9TAC was estimated by an interpolation of the experimental values of CnTAC for the even numbers of n ) 8-18 (refs 26, 27, and 38). (38) Lindblom, G.; Lindman, B. J. Phys. Chem. 1973, 77, 2531. (39) Marrink, S. J.; Tieleman, D. P.; Mark, A. E. J. Phys. Chem. B 2000, 104, 12165. (40) Lauterwein, J.; Bo¨sch, C.; Brown, L. R.; Wu¨thrich, K. Biochim. Biophys. Acta 1979, 556, 244.

Shinto et al.

C.1-C.3, which are the same as systems B.1-B.3, respectively, but in a 0.5 M NaCl solution. 3.2.1. In Salt-Free Water. Figure 6 displays the snapshots of the micelles of different sizes for system B at t ) 0 and 10 ns. After 10 ns, the micelles shrink to some extent, and all the surfactants in each micelle point the hydrophilic heads outward and the hydrophobic tails inward, as in parts d-f of Figure 6. To examine the detailed structure of the micelle, we calculated the number densities of t sites, h+ sites, and Cl- counterions as a function of the distance from the center of mass of the micelle, r, using the last 2-ns configurations of t ) 8-10 ns, which are considered to be equilibrated. The results of N ) 30 for system B.2 are shown in Figure 7. The t sites form the hydrocarbon interior of the micelle at the distances of r < 1.4 nm, while the h+ sites exist at the distances of r ) 1.1-1.6 nm. The Cl- counterions are located outside of the micelle and favorably reside near the h+ sites. The number density profiles for systems B.1 and B.3 (data not shown) were similar to those for system B.2. The micelles have the radii of 1.23, 1.37, and 1.49 nm for systems B.1B.3, respectively, which were calculated as the average distance of the h+ sites from the center of mass of the micelle. The density profiles of N ) 30 for system B.2 quantitatively agree with those of the MD study,41 where a micelle of 30 C10TAC surfactants in 2166 water molecules at 0.67 M was simulated for 0.215 ns. Thus, although no water molecules of the solvent are explicitly considered in our LD technique, it can mimic the average structure of a micelle represented by the atomistic model MD simulation. We emphasize that our LD simulations computationally cost far less than the MD simulations, which require about 11 000-22 000 water molecules for the systems corresponding to systems B.1-B.3. Also, parts d-f of Figure 6 indicate that the micelle of N ) 20 is less spherical than those of N ) 30 and 40. The principal moments of inertia of the micelle, I1, I2, and I3 (I1 g I2 g I3), were calculated to obtain more detailed information of the micelle shape. Figure 8 shows the values of I1/I3 and I2/I3 as a function of time t, where each plot is the result averaged over 2 ps. In parts b and c of Figure 8, both I1/I3 and I2/I3 are close to unity (i.e., I1 ≈ I2 ≈ I3). Because the characteristic length lR of the micelle along the R axis (R ) 1, 2, 3) may be estimated as lR ) (IR/ Nmsurf)-1/2, where msurf is the mass of the surfactant employed,22b the micelles of N ) 30 and 40 are nearly spherical (i.e., l1 ≈ l2 ≈ l3). In part a of Figure 8, the result of I1/I3 ≈ I2/I3 > 1 suggests that the micelle of N ) 20 has a prolate ellipsoidal shape (i.e., l1 ≈ l2 < l3). Furthermore, because I1/I3 and I2/I3 for N ) 20 fluctuate more significantly than those for N ) 30 and 40, the micelle of N ) 20 is unstable and unfavorable compared with the larger micelles. Tieleman et al.33 carried out the MD simulations of dodecylphosphocholine micelles for 11 ns, where three different sizes were chosen as N ) 40, 54, and 65, which are around Nav ) 50-60.40 The smallest micelle of N < Nav was the least stable, and the larger micelles of N > Nav showed no significant difference in the shape and the stability.33 The result is similar to that of the MD study for nonionic surfactant micelles of octyl glucoside42 and the previously mentioned result of our LD study. In contrast, Maillet et al.35 successfully observed the fragmentation of a C9TAC micelle of N ) 48 (>Nav ) 25)36 into two small micelles of N ) 15 and 29 and four monomers (41) Bo¨cker, J.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1994, 98, 712. (42) Bogusz, S.; Venable, R. M.; Pastor, R. W. J. Phys. Chem. B 2000, 104, 5462.

LD Simulations of Cationic Surfactants

Langmuir, Vol. 20, No. 5, 2004 2023

Figure 6. Snapshots of a micelle of N surfactants of t10h+Cl- in water at Csurf ) 0.1 M: (a, d) N ) 20, system B.1; (b, e) N ) 30, system B.2; and (c, f) N ) 40, system B.3. Panels a-c and d-f show the configurations at t ) 0 and 10 ns, respectively. Green and red spheres denote t and h+, respectively. The Cl- counterions involved are not shown for the sake of visual clarity.

after 1.1 ns: unfortunately, the initial configuration of the surfactants was not explained in detail, although it would largely influence the stability of the micelles. Thus, the stability of surfactant micelles with different sizes remains to be scrutinized by molecular-level simulations. 3.2.2. In an Electrolyte Solution. According to numerous experiments, the addition of electrolytes to a surfactant solution causes the reduction of the cmc, where the electrolytes are considered to stabilize the surfactant aggregates. To investigate the effects of an electrolyte on the micelles, we employed systems C.1-C.3 in Table 3, which are respectively the same as systems B.1-B.3 but

in a 0.5 M NaCl solution. The number density profiles for system C.2 are shown in Figure 9. A comparison between Figures 7 and 9 indicates that the densities of the t and the h+ sites for system C.2 are almost the same as those for system B.2 and that the density of the Cl- counterions near the h+ sites is higher than that for system B.2. While the Cl- counterions reside near the positively charged surface of the micelle, the Na+ co-ions stay away from the surface, as in Figure 9. Accordingly, the electric double layer is formed around the micelle of the cationic surfactants, where the electrostatic repulsion between the neighboring h+ sites is screened by the Cl- counterions.

2024

Langmuir, Vol. 20, No. 5, 2004

Shinto et al.

Figure 7. Number densities of t, h+, and Cl- of a micelle of size N ) 30 in water (system B.2), as a function of the distance from the center of mass of the micelle. Each density is the result averaged over 20 000 configurations during the last 2 ns of simulation (i.e., t ) 8-10 ns). The inset focuses on the densities of the charged sites. Figure 10. Same as Figure 8 but in a 0.5 M NaCl solution (system C).

Figure 8. Values of I1/I3 (black lines) and I2/I3 (gray lines) of a micelle of N surfactants in water (system B) as a function of time t: (a) N ) 20; (b) N ) 30; and (c) N ) 40. The principal moments of inertia of the micelle are represented by I1, I2, and I3, where I1 g I2 g I3. Each plot is the result averaged over 2 ps.

Figure 11. Probabilities of the h+ sites with the coordination number nhfCl of the Cl- counterions in the shell around each h+ site for a micelle of N surfactants in solutions of Csalt ) 0 (open bars, system B) and Csalt ) 0.5 M (filled bars, systems C): (a) N ) 20, systems B.1 and C.1; (b) N ) 30, systems B.2 and C.2; and (c) N ) 40, systems B.3 and C.3. Each distribution is the result averaged over 20 000 configurations during the last 2 ns of simulation (i.e., t ) 8-10 ns).

Figure 9. Same as Figure 7 but in a 0.5 M NaCl solution (system C.2). The density profile of the Na+ co-ions is also shown.

Similar results were observed for systems C.1 and C.3 (data not shown). Figure 10 displays the variations of I1/I3 and I2/I3 for systems C.1-C.3. A comparison in part a between Figures 8 and 10 demonstrates that the values of I1/I3 and I2/I3 for system C.1 are close to unity and the fluctuations are significantly reduced, compared with system B.1. These results indicate that the micelle of N ) 20 in a 0.5 M NaCl solution becomes more spherical and stable than that in

water. For the larger micelles of N ) 30 and 40, there is no significant difference in the shape and the stability between the salt-free and the saline systems, as in parts b and c of Figures 8 and 10. To clarify the effects of salt addition, we calculated the number of the Cl- counterions in the coordination shell around each h+ site, nhfCl, and then evaluated the probability distribution of the h+ sites with nhfCl, P(nhfCl), and the average coordination number, n j hfCl, as shown in Figure 11 and Table 4, respectively. Here, the coordination shell is defined as the region of rhfCl e 0.9 nm, at which the PMF for the TMA+-Cl- pair exhibits the second barrier, as shown in part b of Figure 1; accordingly, the present shell consists of the first and the second coordination shells for the Cl- ions around the TMA+. In the salt-

LD Simulations of Cationic Surfactants

Langmuir, Vol. 20, No. 5, 2004 2025

Table 4. Average Numbers of Cl- Counterions in the Coordination Shell of the h+ Sitea system

n j hfCl

system

n j hfCl

difference

B.1 B.2 B.3

0.35 0.54 0.68

C.1 C.2 C.3

0.80 0.91 1.06

0.45 0.37 0.38

a

The definition of the coordination shell is given in text.

free systems (Csalt ) 0), the probability of nhfCl ) 0 for system B.1 is much higher than those for systems B.2 and B.3, while the probabilities of nhfCl > 0 for system B.1 are lower. These results indicate that the repulsion between the h+ sites is less screened by the Cl- counterions in system B.1 than in systems B.2 and B.3, which leads to the instability of the micelle in system B.1. In the saline systems (Csalt ) 0.5 M), a micelle of every size exhibits the reduced probability of nhfCl ) 0 and the increased probabilities of nhfCl > 0, compared with the salt-free systems; this tendency is remarkable for the smallest micelle of N ) 20. Table 4 also reveals that n j hfCl for N ) 20 is greatly increased by the salt addition, compared with those for N ) 30 and 40. Consequently, the addition of salt is more influential in stabilizing the smaller micelle of ionic surfactants. To our best knowledge, only Maillet et al.35 have reported the molecular simulation study on how the salt addition influences a surfactant micelle. They carried out MD simulations of a cylindrical micelle of 48 C9TAC surfactants in both water and a NaCl solution for a few nanoseconds. It was observed that the cylindrical micelle in water broke down into two small spherical micelles, while such fragmentation did not occur in the NaCl solution. These results also demonstrate the stabilization of micelles of ionic surfactants by salt addition. The next goal of the computer simulations of surfactant molecules will be to predict not only the cmc of a surfactant solution and the average size of the surfactant aggregates but also the salt effects on them, though it is rather ambitious within the current computational power. 4. Conclusion We have proposed a LD technique for the simulation of surfactant molecules in aqueous solutions, where no water molecules of the solvent are explicitly treated, but the

effects are incorporated using both the SDCs of the solutes and the PMFs between them in water. The present LD technique was employed to simulate (i) the self-assembly of C10TAC surfactants in water and (ii) the micelles of different sizes in water and in a 0.5 M NaCl solution. In the first simulation starting from the disaggregate configuration, the surfactant aggregates grew by addition of monomers and small aggregates to them rather than by fusion of large aggregates into larger ones within 12 ns: this behavior corresponds to the early stage of surfactant self-assembly and indicates that the system did not completely reach the equilibrium states for the period of several nanoseconds. In the second simulation starting from the spherical micelles, the small micelle of N ) 20 surfactants was unstable and less spherical in water compared with the larger micelles of N ) 30 and 40, but it exhibited a stable spherical shape in the NaCl solution where no significant difference among these micelles was found in the shape and the stability. It is worth noting that, as far as the average structure of a micelle of 30 C10TAC surfactants in water is concerned, the results of our LD simulation coincide with those of the MD simulation. These results demonstrate that our LD technique can capture the characteristics of ionic surfactants in aqueous solutions despite the implicit treatment of water molecules of the solvent. Our LD approach computationally costs far less than the MD and MC methods with atomistic models, and it allows us to change the concentration(s) of the surfactants and salt ions involved without changing the computational cost because of our implicit treatment of water molecules. While we implemented the present LD simulations of several nanoseconds using single central-processing-unit computers, the use of massive parallel computers realizes the simulations of microseconds. Thus, the present LD method will be useful for the long-time large-scale simulation of surfactant molecules in aqueous solutions. Acknowledgment. H.S. thanks the Association for the Progress of New Chemistry, Japan, for partly supporting the present study in 2000/2001. Visualization of simulation data was performed at the Academic Center for Computing and Media Studies, Kyoto University. LA035874O