SurfaceâBulk Partition of Surfactants Predicted by Molecular

Article pubs.acs.org/JPCB

Surface−Bulk Partition of Surfactants Predicted by Molecular Dynamics Simulations Chunwei Yang and Huai Sun* School of Chemistry and Chemical Engineering and Ministry of Education Key Laboratory of Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai, 200240, China ABSTRACT: The surface−bulk partition of surfactants is a piece of important information for understanding the performance of surfactants. In this work, we predicted surface−bulk partitions for three ionic surfactants (sodium hexyl sulfate, sodium nonyl sulfate, sodium dodecyl sulfate) using molecular dynamics simulations and a coarse-grained force field. The simulations were performed with the bilayer model that represents thin films. By setting initial configurations using chemical potentials calculated for infinitely dilute solutions, we obtained a sufficient number of independent samples for each simulation which ran at least 2 μs. Predicted surface concentrations and surface tensions are in good agreement with the experimental data. The surface adsorption isotherms, surface structures, and critical micelle concentrations (CMCs) derived from the simulated data are in reasonable agreement with the experimental data and previously calculated data. For surfactants with a short alkyl chain (SHS), stable oligomers are formed before micelles in the bulk. measurements.6 For computational studies of surfactant-related phenomena and properties such as surface tension and surface elasticity,7 the surface concentration is usually a required input.8−15 Using experimentally determined surface concentration is not feasible for computations with newly designed surfactants. Therefore, in addition to fundamental interest, there is a practical motivation for being able to predict the surface−bulk partition. Computer simulations have been carried out to study the micellizations of surfactants in solutions by several research groups.16−24 The latest advances are large-scale molecular dynamics (MD),19 Monte Carlo (MC),18 and dissipative particle dynamics (DPD)24 simulations based on either a united atom force field (UAFF)17 or a coarse-grained force field (CGFF).25,26 Simulations using a CGFF are about 3 orders of magnitude faster than those based on an atom-specific force field. With good parametrization, quantitative or semiquantitative predictions of surfactant induced properties have been obtained using CGFF.20−22 However, predicting the surface− bulk partition of surfactants raises additional challenges, as the equilibrium between the surface and the bulk involves two topologically very different domains. To our best knowledge, direct simulation of surfactant distributions between the surface and the bulk has not been reported yet. In this work, we carried out direct MD simulations using the bilayer model that represents thin films to study the surface− bulk equilibria of sodium alkyl sulfates (SASs) with variable

1. INTRODUCTION Surfactants are a group of molecules each containing hydrophobic groups (tails) and hydrophilic groups (heads). Due to their amphiphilic feature, surfactants are found in many academic and industrial applications such as porous material syntheses, personal care chemicals, drug delivery agents, manufacturing solutions, and oil recovering additives. In aqueous solutions, surfactant molecules diffuse in water and aggregate at surfaces. The aggregation of surfactants on the surface decreases the surface tension of the solution. As the total concentration of surfactants increases, the number of surfactants adsorbed at the surface increases and the surface tension decreases. This trend continues until the critical micelle concentration (CMC) of the solution is reached. At the CMC, micelles are formed in the solution and the surface concentration reaches the maximum and is stabilized.1 The partition of surfactants between the surface and the bulk is a key property that characterizes the performance of surfactants. The surface concentration, namely, the adsorbed amount of surfactants in unit area, can be deduced by measuring the surface tensions against the total concentrations of the solution. By plotting the measured surface tensions against the total concentrations, the surface and bulk concentrations can be calculated from the slopes of the curve on the basis of Gibbs adsorption theory.2 Another, less popular but more sophisticated method to measure the surface concentrations is to use the neutron reflection technique.3−5 In this approach, the surfactant molecules need to be deuterated and placed in neutron null-reflecting solvent. Due to the approximations and complexity in experiments, surface tensions are often reported with large variations from different © 2014 American Chemical Society

Received: July 7, 2014 Revised: August 2, 2014 Published: August 8, 2014 10695

dx.doi.org/10.1021/jp506768b | J. Phys. Chem. B 2014, 118, 10695−10703

The Journal of Physical Chemistry B

Article

Figure 1. Snapshot of the bilayer model representing both surface and bulk phases of surfactant solution using SNS as an example. The purple, cyan, pink, and lime balls denote beads of SO4−, −CH2CH2CH2−, CH3CH2−, and Na+(H2O)3, respectively. The mauve dots denote water beads. The surfaces are covered by surfactants, and a micelle is formed in the bulk.

Table 1. Simulation Models Used in This Worka

alkyl chain lengths. The simulations were based on the Shinoda−DeVane−Klein (SDK) CGFF.26,27 We aimed to address the following questions: (1) Can we accurately predict surface−bulk equilibria for surfactants using direct MD simulations? (2) What are the structures and properties of SASs on the surface and in the bulk? (3) Is it possible to predict surface concentrations and CMCs from these simulations? In the following sections, we first present the simulation models and methods, then discuss the calculated data to answer the above questions, and finally summarize this work.

2. MODELS AND METHODS The bilayer model representing two vapor−liquid interfaces (see Figure 1) was used in this work. This model represents both surface and bulk phases of the surfactant solution. In each simulation box, the bulk was separated from the surfaces by two buffers in order to exclude interactions between surfactants on the surface and surfactants in the bulk. The length of each buffer was 1.5 nm, equal to the cutoff value for evaluation of nonbond interactions. The simulation boxes were fixed at 30 nm in length along the direction normal to the surfaces. The divisions between the surfaces and the buffers were set at where the water density was half of the pure water density (ca. 1.0 g/ cm3). The bulk phases were in the range of 8−14 nm thickness. The lengths in the x- and y-directions of the simulation boxes varied for different molecules. The three-dimensional periodic boundary condition was applied in the simulations, but the vapor phase was large enough (>10 nm) so that the model represents separated thin films. We carried out simulations with sodium hexyl sulfate (SHS), sodium nonyl sulfate (SNS), and sodium dodecyl sulfate (SDS). These surfactants were studied for micellization using the isotropic bulk model22 so that direct comparisons can be made for the bulk properties with our calculations. Table 1 lists the simulation models with different surface area (Lx × Ly) of the simulation boxes, the thickness of the solution which includes the bulk and buffers, the number of water beads (each represents three water molecules), the total number of surfactants (loads) in the simulation boxes, and the simulation time per configuration. The simulations were grouped into two sets. Set I was used for scanning different surfactant loads (the total numbers). For the low loads, the scan was performed with large increments of 20, 40, and 100 for SHS, SNS, and SDS, respectively. For loads close to surface saturation, a finer increment of 4 was used to scan the loads. Set II was designed for evaluating CMCs with larger bulk volumes. NVT MD simulations were carried out using the same conditions as reported in the literature for parametrization of the SDK model:20 10 fs time step, 1.5 nm cutoff without tail

set

surfactant

SA (nm2)

Lsol (nm)

NWB

I

SHS

16

8

1420

SNS

49

12

6520

SDS

196

14

30000

SHS SNS SDS

36 49 196

14 16 20

6420 8740 43700

II

Ntot 20−80 (20) 80−120 (4), 140 40−200 (40) 220−240 (4), 260, 280 100−900 (100) 900−920 (4), 940, 960, 1000 440, 480, 520 200, 220, 240 900, 920, 940

trun (μs) 2 3 2 3 2 4 4 3 2

a

The surface area (SA), thickness of the solution phase (Lsol), number of water beads (NWB), number of surfactants (Ntot),b and simulation time per configuration (trun) are listed. bLarge increments of 20, 40, and 100 for SHS, SNS, and SDS were set for scanning low surfactant loads; fine increments of 4 were used for high surfactant loads for simulation set I.

corrections for LJ interactions, PPPM method and 80 as the relative permittivity for evaluation of electrostatic interactions. The temperature was controlled using the Nosé−Hoover method28 (T = 298 K and τ = 10 ps). The simulations were conducted by loading up surfactants in a step-by-step manner. The initial configurations were constructed by placing surfactants randomly on the surface and in the bulk. The initial partitions between the surface and the bulk were determined by chemical potentials of surfactants at infinitely dilute conditions, estimated using the steered MD method.29 Starting from the lowest load, the simulations were carried out for the time specified in Table 1 to reach equilibration and then moved to the next step by adding the increment number of surfactants. At each load, the convergence was verified by measuring surfactant distributions in the bulk phases, and physical properties were calculated on the basis of the data collected after equilibration time. The statistical uncertainties were estimated using the block averaging method30 with blocks determined on the basis of the correlation time. The simulations were conducted using LAMMPS software31 running in parallel mode. The initial configurations were built by using the PACKMOL package,32 and visualizations were done by using VMD.33 The computational jobs were performed on the “π” system (the high performance computing center) of Shanghai Jiao Tong University. Over 2000 cores and over 450000 core hours were used, and more than 100 μs trajectory data were generated for analysis. 10696



Article

NsSNS:N bSNS = 40:1

3. RESULTS AND DISCUSSION Equilibration. Since the initial configuration strongly impacts the simulation time for reaching equilibrium, a reasonable partition of surfactants between the surface and the bulk in the initial configuration is important for reducing the computational expenses. Before the CMC is reached, the equilibrium in the bilayer simulation box is mainly between the surface and the solution. Therefore, we used the chemical potential of the surfactant between the surface and the bulk to determine the initial partition. At equilibrium, the chemical potential of an adsorbed solvent must be equal to its chemical potential in the bulk, which can be expressed as34 μs0 + RT ln γsCs = μ b0 + RT ln γbC b

NsSDS:N bSDS = 400:1

Figure 3 shows the numbers of surfactants in the bulk phases as functions of simulation time. Three sets of data with loads

(1)

Under dilute condition γ ≈ 1.0, at constant temperature and pressure, we have ⎛C ⎞ Δμ0 = μs0 − μ b0 = −RT ln⎜ s ⎟ ⎝ Cb ⎠

(2)

where Δμ expresses the standard work for transferring a surfactant molecule from the bulk to the surface,35 which can be estimated using the steered MD method on the basis of the Jarzynski relation: 0

Δμ0 = ΔGAB = −kBT ln⟨e−βWAB⟩A

(3)

where WAB is the work performed to force the system to move along one path from state A to B and the angular brackets denote averaging over a canonical ensemble of the reference state A. Figure 2 shows the potential of mean force (PMF) curves calculated for SASs as functions of the distances from the

Figure 3. Numbers of SHS (top), SNS (middle), and SDS (bottom) in the bulk as functions of simulation time. For each surfactant, three loads around the surface saturation points are analyzed. The equilibration times are denoted by triangles. The average values, denoted by straight lines, are calculated from the data after the convergence times.

Figure 2. Potentials of mean force (PMF) calculated for three SASs as the distance from the surface to the bulk.

surface to the bulk. As expected, the free energy differences, e.g., the chemical potentials (eq 3), are higher in the bulk than on the surface. The difference reflects the hydrophobicity of the surfactant, which depends on the structures of the surfactant: the longer the alkyl chain, the larger the chemical potential difference. Taking the average values of the chemical potentials (3.7, 6.0, and 8.1 kcal/mol for SHS, SNS, and SDS, respectively) and the average lengths of surfactants as the surface thicknesses, we obtained the estimated partitions of surfactants as follows:

around the surface saturation points are drawn for each of the surfactants. For each simulation, an equilibration time is identified as the time when the calculated numbers are stabilized. The equilibration time depends on the alkyl chain length. They are less than 50 ns for SHS, less than 600 ns for SNS, and less than 800 ns for SDS. The equilibration time also depends on the load of surfactants. The smaller the load, the shorter the equilibration time that is required. Since the loads of surfactants presented in Figure 3 are around the surface saturation points, the equilibration times represent the high end

NsSHS:N bSHS = 4:1 10697



Article

in our simulations. For most simulations with lower loads (not shown), the equilibration times are shorter than those presented in Figure 3. Table 2 lists the statistical data of Figure 3. The running time minus the equilibration time is the data collection time. The Table 2. Number of Surfactants (Ntot), Running Time (trun), Equilibration Time (teq), Correlation Time (tcor), Number of Independent Samples (S), Average Number of Surfactants in the Bulk (Nb), and Uncertainty (σ) Calculated for Loads near the Surface Saturation Points Ntot

trun (μs)

teq (μs)

tcor (μs)

S

Nb

σ

SHS-80 SHS-100 SHS-120 SNS-200 SNS-240 SNS-280 SDS-800 SDS-920 SDS-1000

3 3 3 3 3 3 4 4 4

0.03 0.04 0.05 0.1 0.6 0.7 0.2 0.3 0.8

0.002 0.003 0.002 0.029 0.015 0.012 0.010 0.088 0.166

743 494 738 50 80 96 190 21 10

22.3 36.1 52.7 2.8 36.8 71.4 0.4 9.5 83.6

0.08 0.13 0.01 0.23 0.24 0.18 0.03 0.10 0.11

Figure 4. Surface concentration (top) and surface tension (bottom) as a function of the total load of surfactant. The vertical dashed lines indicate the transition points at Ntot = 100, Ntot = 220, and Ntot = 920 for SHS, SNS, and SDS, respectively.

correlation time is calculated by integrating the autocorrelation function of the number of surfactants in the bulk after the equilibration time. ∞

tcor =

∫τ=0 ⟨A(0)A(τ)⟩ dτ

the surface tension decreases monotonically with the load increases. In either the surface concentration or surface tension curve, a transition point near the end of the curve can be identified. The transition points correspond to the same load for both surface concentration and surface tension for each surfactant. It is not surprising, as the most sensitive parameter that impacts surface tension is the surface concentration, given the same molecule.38 The transition point corresponds to the surface saturation. Past the transition point, the surface concentration and surface tension become constant in experimental measurement. We did not get the plateaus from the simulations because of the concentration of molecules in the buffers. According to the experimental observation, the concentration of surfactants in the buffer regions should be a metastable state past the transition point, which implies that another transition for molecules in the buffers migrating to the bulk may exist. However, our simulation time was not long enough to observe this phenomenon. Nevertheless, saturated surface concentrations and surface tensions can be calculated using the data collected at the transition points. The calculated data are summarized in Table 3. The calculated surface concentrations are 3.3, 3.5, and 3.8 × 106 mol/m2 and the surface tensions are 42.5, 42.3, and 38.2 mN/m for SHS, SNS, and SDS, respectively. The calculated values agree reasonably well with the experimental data, which unfortunately exhibit large variations.39−42 The order parameter defined as the average of the second Legendre polynomial43 1 S = ⟨3 cos2 θ − 1⟩ (7) 2 was used to analyze surface structures. The brackets in eq 7 denote ensemble average, and θ is the angle between a vector connecting the sulfate bead (SU) and the terminal carbon bead (CT) of the surfactant molecule and the normal direction of surface. The order parameter S varies from −1/2 to 1. It is 0 for

(4)

The number of independent samples is calculated by dividing the data collection time by twice the correlation time. As shown in Table 2, for most simulations, we have obtained sufficient numbers of independent samples. However, for SDS-1000, the number of independent samples is only 10. The average numbers of surfactants in the bulk are calculated after the equilibration time, and the uncertainties of statistically independent samples are calculated by30 σEB =

2

tcor σraw (trun − teq)

(5)

where σraw is the standard deviation of the raw data calculated from the data after the equilibrium time. The uncertainties are lower than 0.24 across the table. Surface Concentrations and Properties. Figure 4 shows the dependence of calculated surface concentration and surface tension on the load (in logarithm of Ntot) of surfactant molecules. A close examination of the simulated data shows that molecules are concentrated in the buffer due to interactions with molecules in the surface, as shown in Figure 1. Therefore, the surface concentration was calculated by dividing the number of molecules in the surface plus the number of molecules in the buffer by the surface area, and averaged over two sides of the simulation box. The surface tensions were computed using the Kirkwood−Buff formula which is related to diagonal components of the pressure tensor as36,37 γ=

Pxx + Pyy ⎫ Lz ⎧ ⎨Pzz − ⎬ 2⎩ 2 ⎭

(6)

where Lz denotes the dimension of the simulation box and Pii denotes the pressure tensor of three directions x, y, and z. For each of the surfactants, the surface concentration increases and 10698



Article

Table 3. Comparison of Calculated and Experimental Surface Tensions and Surface Concentrations Γmax (106 mol/m2)

a

surfactant

calc.

SHS SNS SDS

3.3 ± 0.3 3.5 ± 0.2 3.8 ± 0.2

γ (mN/m) expt. 39

a

2.7, 3.8 4.140 3.1,39 4.2,40 4.041

calc.

expt.

42.5 ± 0.4 42.3 ± 0.5 38.2 ± 0.2

33.5,39 42.0a 41.540 38.4,42 40.041

The data are extrapolated from series data measured for sodium alkyl sulfates (C8−C13) as given in ref 40.

a completely random and isotropic sample, −1/2 for a perfectly parallel orientation, and 1 for a perfectly perpendicular orientation. A sharp drop in the calculated order parameters is usually associated with phase transition. The calculated results for the surfactants are listed in Figure 5. For three

Figure 5. Order parameters calculated for three surfactants on surfaces at different total loads of surfactants.

surfactants, the order parameters approach roughly the same value when the surfaces are saturated. The average value of roughly 0.45 corresponds to an averaged angle of 37°, close to the all-atom simulation results.38,44 The surfactants with different hydrophobic chain lengths show very different structures on the surfaces at low concentrations. The shorter the chain length is, the more “upright” the molecule on the surface. The order parameter for SDS is close to 0 at the low concentration end, indicating a random orientation of SDS on the surface, but the value for SHS at the lowest concentration is 0.4, which is close to the order parameter obtained at the saturated concentration. Bulk Structures and Properties. Figure 6 shows the average numbers of surfactants in the bulk as the total loads increase for the three surfactants. We define the solvation ratio as the slope of the curve: α = ΔNb/ΔNtot

Figure 6. Number of surfactants in the bulk phase as a function of the total load for the three surfactants. Two transition points TP1 and TP2 are identified for SHS and SNS, but one transition point is identified for SDS. The surface saturation points are denoted with vertical dashed lines.

(8)

coincides with the surface saturation point. The solvation ratios are about 0.001 before the transition and 0.92 after the transition. The fact that solvation ratios are close to 1 after the surface saturation points is not surprising, as additional molecules are mostly in the bulk phases after the surfaces are packed with surfactants. In order to understand how the solvation ratios change in the entire range, we selected three loads for each surfactant to analyze corresponding liquid structures. The three loads were 40, 80, and 120 for SHS and 160, 212, and 280 for SNS, representing three states before TP1, between TP1 and TP2, and after TP2 for SHS. The three loads were selected as before the transition point at 800, right after the transition at

For SHS, the solvation ratio can be divided into three linear segments separated by two transition points, TP1 and TP2. TP1 is at an early stage of the total load at 50 and TP2 is at the total load of 100 which coincides with the surface saturation point denoted by a vertical line in Figure 6. The solvation ratio is 0.16 before TP1, 0.58 between TP1 and TP2, and 0.83 after TP2. For SNS, the first segment is linear with a low solvation ratio of 0.01. The second segment is nonlinear, the solvation ratio increases rapidly as the loads change from 200 to 220, and the third segment is linear again with a solvation ratio of 0.92. In the case of SDS, only two segments can be identified with the transition point located at the total load of 920, which 10699



Article

Figure 7. Cluster size distribution of three surfactants before TP1, between TP1 and TP2, and after TP2. An enlarged chart is shown for clarity for SHS-120.

1% of nothing before TP1. At this stage, molecules are dissolved in the bulk phase mostly as monomers; the solvation ratio is determined by the solubility of the monomer in water. After TP1 and before TP2, the aggregation is obviously visible, the cluster sizes span from 1 to 10, and the probability of the cluster sizes decreases rapidly from 70% for monomers to 0.1% for 10-mers. The formation of oligomers with their hydrophilic heads toward the water increases the solvation ratio. After TP2, as the cluster sizes increase further, the percent of monomer is reduced to 58% and more molecules are aggregated to form clusters in sizes larger than 10. The maximum aggregation number is 19, which is in agreement with previous atomic simulation data of 2017 and with experimental data of 17.45 In addition, a valley at 14 is close to the value of 11 as reported from the bulk simulations.22 The cluster size distributions of SNS in three loads are shown in the second row. Before TP1, the probability of finding nothing in the bulk is 45% and the probability of finding monomers is 55%. Again, the solvation ratio is determined by the solubility of monomers in water. After TP1 and before TP2, the percentage of monomers increases to 90%, small oligomers are formed (the largest is 6-mers with 0.05% probability), and there is a 1% chance of finding nothing in the bulk. Because the solubility of a monomer is low, the change of solvation ratio takes place in a very short range (ca. 20 in total load), resulting in a nonlinear curve in the calculated solvation ratio. After TP2, the probability of finding a monomer is reduced to 37% and

920, and far after the transition at 1000 for SDS. For each load, we calculated the probability of finding a cluster in size k in the bulk phase by A

P(k) =

∑i Pi(k) A

(9)

where A is the number of samples (frames in the trajectory) collected from the trajectory and Pi(k) is the probability of a cluster in size k in sample i. If there is no surfactant in the sample: P(0) =1 i

(10)

Otherwise: Pi(k) =

Ni(k) ∑k Ni(k)

(11)

where Ni(k) is the number of clusters with size k in sample i. The clusters were identified by calculating the radial distribution function (RDF) for the terminal carbon beads (CT). The first valley at 7.5 Å was used as the minimum distance when two surfactants are considered as in one cluster.27 Figure 7 shows the probability of observed cluster sizes for the three states selected as above. The first row in Figure 7 shows the probability of finding different cluster sizes at the three loads for SHS. The trajectory contains 92% of monomers, 6% of dimers, 1% of trimers, and 10700



Article

Figure 8. Bead density distributions of surfactants in three load states: before TP1, between TP1 and TP2, and after TP2 along the direction normal to the surfaces. CT, CM, SU, SOD, and W denote SO4−, −CH2CH2CH2−, CH3CH2−, Na+(H2O)3, and (H2O)3 beads, respectively.

larger clusters with sizes ranging from 21 to 46 are found. The maximum aggregation number (36) is consistent with the experimental result of 3345 and the simulated result of 37.22 The last row in Figure 7 shows the probability of finding clusters in different sizes for SDS in the three loads. Before the TPs, the probabilities are 92% for finding nothing and 8% for finding monomers in the bulk. On average, the solvation ratio (0.001) is very low. After the TPs, the surface is saturated and the solvation ratio becomes 1.0. We see 21% monomers, and the rest (79%) are oligomers with sizes ranging from 7 to 10. As the total load increases to 1000, 3% are monomers and the rest (97%) are distributed in four clusters with sizes of 17, 18, 32, and 33. The clusters are fairly stable; we did not observe any fission or fusion of the clusters in the simulation time up to 4 μs. As shown in Table 2, the number of independent samples in this load is only 19. In order to observe a continuous distribution like that found for the SHS and SNS systems, the simulation time needs to be extended much further. The density distributions along the normal direction of the simulation boxes are shown in Figure 8. For each surfactant, distributions are analyzed for the three different loads discussed above. The charts show that most surfactants are located at the two surfaces, with a few in the bulk. The increased population of surfactants in the bulk is more obviously seen by the dent in the water population. For SHS, a small number of surfactants are found in the bulk at a load of 40, but the number increases obviously at loads 80 and 120. For SNS and SDS, there is

almost nothing in the bulk at loads of 200 and 800. The populations of surfactant are visible in the bulk at loads of 212 and 920, and obviously visible at loads of 280 and 1000. Finally, we estimated CMCs for the surfactants using the concentrations of monomers in the bulk phases. It has been reported 46,47 that the counterion association must be considered for determining the CMCs from the raw simulation data. We calculated the CMCs using the same theoretical correction method47 that has been applied to predict CMCs for the same surfactants, the same force field but a different simulation model.22 The calculated results are compared with the experimental data and previous calculation data22 in Figure 9. For comparison, results obtained from the raw data are also displayed in Figure 9. Figure 9 shows that our data, obtained from the bulk phases of the bilayer model, are quite close to those obtained from the isotropic bulk model.22 Despite that the surface tensions and surface concentrations were predicted very well in comparison with the experimental data, all calculated CMC results are lower than the experimental data and the discrepancies increase as the alkyl chain increases; the problem is more pronounced for SDS which has the lowest CMC among the three surfactants studied. The problem is unlikely due to the force field quality, as the SDK model has been validated by calculating the surface and bulk properties in this and previous works.26,27 The discrepancies are likely due to the fact that the simulation times are not enough to equilibrate the micelle size distributions 10701



■

Article

AUTHOR INFORMATION

Corresponding Author

*Phone: +86-21-5474-8987. Fax: +86-21-5474-1297. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was supported in part by the Nature Science Foundation of China (No. 21173146) and National Basic Research Program (No. 2014CB239702). Figure 9. Comparison of calculated and experimental CMCs. The black squares are calculated by using the raw data of free monomer concentrations. The red circles are calculated using the counterion corrections.47 The blue up-triangles are calculated using the isotropic bulk model,22 and the pink down-triangles are experimental data.

(1) Prosser, A. J.; Franses, E. I. Adsorption and Surface Tension of Ionic Surfactants at the Air−water Interface: Review and Evaluation of Equilibrium Models. Colloids Surf., A 2001, 178, 1−40. (2) Chattoraj, D. K.; Birdi, K. Adsorption and the Gibbs Surface Excess; Plenum Press: New York, 1984. (3) Taylor, D.; Thomas, R.; Penfold, J. The Adsorption of Oppositely Charged Polyelectrolyte/Surfactant Mixtures: Neutron Reflection From Dodecyl trimethylammonium Bromide and Sodium Poly (Styrene Sulfonate) at the Air/water Interface. Langmuir 2002, 18, 4748−4757. (4) Simister, E.; Thomas, R.; Penfold, J.; Aveyard, R.; Binks, B.; Cooper, P.; Fletcher, P.; Lu, J.; Sokolowski, A. Comparison of Neutron Reflection and Surface Tension Measurements of the Surface Excess of Tetradecyltrimethylammonium Bromide Layers at the Air/Water Interface. J. Phys. Chem. 1992, 96, 1383−1388. (5) Li, Z.; Dong, C.; Thomas, R. Neutron Reflectivity Studies of the Surface Excess of Gemini Surfactants at the Air-water Interface. Langmuir 1999, 15, 4392−4396. (6) Burlatsky, S. F.; Atrazhev, V. V.; Dmitriev, D. V.; Sultanov, V. I.; Timokhina, E. N.; Ugolkova, E. A.; Tulyani, S.; Vincitore, A. Surface Tension Model for Surfactant Solutions at the Critical Micelle Concentration. J. Colloid Interface Sci. 2013, 393, 151−160. (7) Rosen, M. J.; Kunjappu, J. T. Surfactants and Interfacial Phenomena; John Wiley & Sons: New York, 2012. (8) Tarek, M.; Tobias, D. J.; Klein, M. L. Molecular Dynamics Simulation of Tetradecyltrimethylammonium Bromide Monolayers at the Air/Water Interface. J. Phys. Chem. 1995, 99, 1393−1402. (9) Wijmans, C. M.; Linse, P. Surfactant Self-assembly at a Hydrophilic Surface. A Monte Carlo Simulation Study. J. Phys. Chem. 1996, 100, 12583−12591. (10) Kuhn, H.; Rehage, H. Molecular Dynamics Computer Simulations of Surfactant Monolayers: Monododecyl Pentaethylene Glycol at the Surface between Air and Water. J. Phys. Chem. B 1999, 103, 8493−8501. (11) Kuhn, H.; Rehage, H. Molecular Orientation of Monododecyl Pentaethylene Glycol (C12E5) Surfactants at Infinite Dilution at the Air/Water Interface. A Molecular Dynamics Computer Simulation Study. Phys. Chem. Chem. Phys. 2000, 2, 1023−1028. (12) Kuhn, H.; Rehage, H. Molecular Orientation of Monododecyl Pentaethylene Glycol at Water/Air and Water/Oil Interfaces. A Molecular Dynamics Computer Simulation Study. Colloid Polym. Sci. 2000, 278, 114−118. (13) Bandyopadhyay, S.; Chanda, J. Monolayer of Monododecyl Diethylene Glycol Surfactants Adsorbed at the Air/Water Interface: A Molecular Dynamics Study. Langmuir 2003, 19, 10443−10448. (14) Chanda, J.; Bandyopadhyay, S. Molecular Dynamics Study of a Surfactant Monolayer Adsorbed at the Air/Water Interface. J. Chem. Theory Comput. 2005, 1, 963−971. (15) Chanda, J.; Chakraborty, S.; Bandyopadhyay, S. Monolayer of Aerosol-OT Surfactants Adsorbed at the Air/Water Interface: An Atomistic Computer Simulation Study. J. Phys. Chem. B 2005, 109, 471−479.

which involve micelle fissions and fusions. Generally speaking, the ergodicity of simulation of low concentration solution is poor; advanced sampling techniques such as hyperdynamics48 and replica exchange methods49,50 may be helpful to improve the predictions.

4. CONCLUSIONS Large scale MD simulations using the SDK coarse-grained force field were carried out to predict the surface−bulk partitions of surfactants using the bilayer model successfully. By setting reasonable initial configurations according to computed chemical potentials in infinite dilute solution, the simulation times for reaching surface−bulk equilibria were less than 0.8 μs for the three surfactants and different loads considered in this work. In most cases, sufficient numbers of samples were obtained, which ensured a reasonably good statistics in the data analysis. The surface saturation points were identified from the surface−bulk partition curves. The surface concentrations and surface tensions calculated at the saturated surfaces agree well with the experimental data. The surface structures, the density profiles, and the cluster size distributions of surfactant agree well with previous calculations using different simulation models. From the analysis of the solvation ratios, we found that the surface−bulk partitions depend on the solubility or hydrophilicity of the surfactants. For the most hydrophilic surfactant (SHS), many oligomers are formed in the bulk far before the CMC. For SNS, oligomers are not very stable but can be observed just before the CMC. For SDS, the oligomers are not stable and only the formation of micelles stabilizes the molecules in the bulk. By counting free monomers in the bulk, CMCs can be calculated using the bilayer model. The predicted CMC values agree well with those predicted using the isotropic bulk model. Overall, the predicted CMCs are lower than the experimental value, and the discrepancies increase with the increase of alkyl chain length. Simulation efficiency needs to be further enhanced in order to predict low CMCs accurately. 10702



Article

(16) Mackie, A. D.; Panagiotopoulos, A. Z.; Szleifer, I. Aggregation Behavior of a Lattice Model for Amphiphiles. Langmuir 1997, 13, 5022−5031. (17) Sammalkorpi, M.; Sanders, S.; Panagiotopoulos, A. Z.; Karttunen, M.; Haataja, M. Simulations of Micellization of Sodium Hexyl Sulfate. J. Phys. Chem. B 2011, 115, 1403−1410. (18) Jusufi, A.; Sanders, S.; Klein, M. L.; Panagiotopoulos, A. Z. Implicit-solvent Models for Micellization: Nonionic Surfactants and Temperature-dependent Properties. J. Phys. Chem. B 2011, 115, 990− 1001. (19) Sanders, S. A.; Panagiotopoulos, A. Z. Micellization Behavior of Coarse Grained Surfactant Models. J. Chem. Phys. 2010, 132, 1149021−114902-9. (20) Shinoda, W.; DeVane, R.; Klein, M. L. Zwitterionic Lipid Assemblies: Molecular Dynamics Studies of Monolayers, Bilayers, and Vesicles Using a New Coarse Grain Force Field. J. Phys. Chem. B 2010, 114, 6836−6849. (21) Levine, B. G.; LeBard, D. N.; DeVane, R.; Shinoda, W.; Kohlmeyer, A.; Klein, M. L. Micellization Studied by Gpu-accelerated Coarse-grained Molecular Dynamics. J. Chem. Theory Comput. 2011, 7, 4135−4145. (22) 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. (23) Shinoda, W.; DeVane, R.; Klein, M. L. Computer Simulation Studies of Self-assembling Macromolecules. Curr. Opin. Struct. Biol. 2012, 22, 175−186. (24) 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. (25) 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. (26) Shinoda, W.; Devane, R.; Klein, M. Multi-property Fitting and Parameterization of a Coarse Grained Model for Aqueous Surfactants. Mol. Simul. 2007, 33, 27−36. (27) Shinoda, W.; DeVane, R.; Klein, M. L. Coarse-grained Force Field for Ionic Surfactants. Soft Matter 2011, 7, 6178−6186. (28) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (29) Park, S.; Schulten, K. Calculating Potentials of Mean Force From Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 5946−5961. (30) Flyvbjerg, H.; Petersen, H. G. Error Estimates on Averages of Correlated Data. J. Chem. Phys. 1989, 91, 461−466. (31) Plimpton, S. Fast Parallel Algorithms for Short-range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (32) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (33) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (34) Wei, Y.; Latour, R. A. Determination of the Adsorption Free Energy for Peptide−Surface Interactions by SPR Spectroscopy. Langmuir 2008, 24, 6721−6729. (35) Danov, K. D.; Kralchevsky, P. A. The Standard Free Energy of Surfactant Adsorption at Air/Water and Oil/Water Interfaces: Theoretical vs. Empirical Approaches. Colloid J. 2012, 74, 172−185. (36) Rao, M.; Berne, B. On the Location of Surface of Tension in the Planar Interface between Liquid and Vapour. Mol. Phys. 1979, 37, 455−461. (37) Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Surface Tension. J. Chem. Phys. 1949, 17, 338−343. (38) Shen, Z.; Sun, H.; Liu, X.; Liu, W.; Tang, M. Stability of Newton Black Films Under Mechanical Stretch − A Molecular Dynamics Study. Langmuir 2013, 29, 11300−11309.

(39) Suárez, M. J.; López-Fontán, J. L.; Sarmiento, F.; Mosquera, V. Thermodynamic Study of the Aggregation Behavior of Sodium nHexyl Sulfate in Aqueous Solution. Langmuir 1999, 15, 5265−5270. (40) Varga, I.; Mészáros, R.; Gilányi, T. Adsorption of Sodium Alkyl Sulfate Homologues at the Air/Solution Interface. J. Phys. Chem. B 2007, 111, 7160−7168. (41) Lu, J.; Thomas, R.; Penfold, J. Surfactant Layers at the Air/ Water Interface: Structure and Composition. Adv. Colloid Interface Sci. 2000, 84, 143−304. (42) Jang, S. S.; Goddard, W. A. Structures and Properties of Newton Black Films Characterized Using Molecular Dynamics Simulations. J. Phys. Chem. B 2006, 110, 7992−8001. (43) Egberts, E.; Berendsen, H. Molecular Dynamics Simulation of a Smectic Liquid Crystal with Atomic Detail. J. Chem. Phys. 1988, 89, 3718−3732. (44) Domínguez, H. Computer Simulations of Surfactant Mixtures at the Liquid/Liquid Interface. J. Phys. Chem. B 2002, 106, 5915−5924. (45) Rassing, J.; Sams, P.; Wyn-Jones, E. Kinetics of Micellization from Ultrasonic Relaxation Studies. J. Chem. Soc., Faraday Trans. 2 1974, 70, 1247−1258. (46) 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. (47) Jusufi, A.; LeBard, D. N.; Levine, B. G.; Klein, M. L. Surfactant Concentration Effects on Micellar Properties. J. Phys. Chem. B 2012, 116, 987−991. (48) Voter, A. F. Hyperdynamics: Accelerated Molecular Dynamics of Infrequent Events. Phys. Rev. Lett. 1997, 78, 3908−3911. (49) Earl, D. J.; Deem, M. W. Parallel Tempering: Theory, Applications, and New Perspectives. Phys. Chem. Chem. Phys. 2005, 7, 3910−3916. (50) Kührová, P.; Banás,̌ P.; Best, R. B.; Šponer, J.; Otyepka, M. Computer Folding of RNA Tetraloops? Are We There Yet? J. Chem. Theory Comput. 2013, 9, 2115−2125.

10703


SurfaceâBulk Partition of Surfactants Predicted by Molecular

Recommend Documents

SurfaceâBulk Partition of Surfactants Predicted by Molecular