Prediction of Surface and Bulk Partition of Nonionic Surfactants Using

Nov 30, 2015 - Prediction of Surface and Bulk Partition of Nonionic Surfactants Using Free Energy Calculations ... E-mail address: [email protected]...
0 downloads 13 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Prediction of Surface and Bulk Partition of Nonionic Surfactants Using Free Energy Calculations Zhe Shen† 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 ¶ State Key Laboratory of Inorganic Synthesis & Preparative Chemistry, College of Chemistry, Jilin University, Changchun 130012, China S Supporting Information *

ABSTRACT: The surface-bulk partition of nonionic surfactants was predicted by calculating chemical potential differences using two simulation boxes representing the surface and the bulk phases separately. A published coarse grained force field was modified and validated for this application. Thermodynamic integration (TI) was applied to compute excess chemical potentials. The high concentration surface was stabilized by applying an external harmonic potential, and the bulk was treated as ideal solution, which was confirmed by simulation using a lattice model at conditions near the critical micelle concentration (∼10−5 mol/L). Based on the calculated chemical potential differences with precision of ca. 1 kJ/mol, the equilibria of surface-bulk concentration was predicted well in comparison with the experimental data. accurate. However, for SDS with a CMC on the order of 10−3 mol/L, simulations of the microsecond time scale were not sufficient, and the results had large uncertainties and deviations from the experimental data. For surfactants with even lower CMCs, such as nonionic surfactants with CMCs on the order of 10−4 to 10−6 mol/L, direct simulation of surface-bulk partition is not feasible. Therefore, more efficient methods need to be developed. One of the potential solutions to this problem is to use the Gibbs Ensemble Monte Carlo (GEMC) method,14,15 which can predict phase equilibria. However, the efficiency of GEMC simulations depends on the states of the two phases. For the surface-bulk equilibria of the surfactant, both phases are condensed states where the molecules are arranged very differently, making it highly challenging to design efficient MC moves that for GEMC simulations. In this work, we predict the surface-bulk phase equilibria by calculating the chemical potential difference of transferring a surfactant molecule between the two phases using two uncoupled simulation boxes representing each of the two phases. Because the simulations are carried out in smaller simulation boxes, this approach is more efficient than direct simulation and GEMC simulation methods. The challenge with this method is calculating the chemical potentials accurately.

I. INTRODUCTION The surface-bulk partition1−4 of surfactant is one of the most important properties determining the performance of surfactant molecules. The surface concentration, e.g., the adsorbed amount of surfactants per unit area at the interface, is usually deduced from a surface tension profile measured as a function of the total (bulk) concentration based on Gibbs adsorption theory.1 Due to the approximate and complex nature of experimental measurements, the surface concentration is often reported with large variations.5 With computational chemical techniques, the surface concentration is commonly a prerequisite6−12 for predicting surfactant properties such as surface tension and surface elasticity.1 For novel surfactant molecules in which the surface concentration is unknown, there is a practical motivation for being able to predict the surfacebulk partition independently. In a previous work, 13 we used molecular dynamics simulations of up to several microseconds (μs) in duration to predict the surface-bulk partition of three ionic surfactants, sodium hexyl sulfate (SHS), sodium nonyl sulfate (SNS), sodium dodecyl sulfate (SDS) using a bilayer model representing a thin film. The surface adsorption isotherms, spatial arrangement, and critical micelle concentrations (CMC) of the surfactant solutions were predicted with reasonable agreement to experimental data. However, the accuracy of these results were correlated with the solubility or the CMC of the surfactant. For SHS and SNS, with CMCs on the order of 10−1 to 10−2 mol/L, the simulations generated a sufficient amount of data to predict results, which were statistically reliable and © XXXX American Chemical Society

Received: October 19, 2015 Revised: November 22, 2015

A

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B The chemical potential of a component in phase Γ can be generally written as a function of the concentration X of the component: μ Γ = μ Γo + RT ln XΓ

II. MODELS AND METHODS The calculations in this study were conducted for a group of nonionic surfactants, denoted as C12(EO)n, where n = 3, 5, 6, 8 and corresponds to the number of ethylene oxide (EO) units and with a fixed alkyl chain length of 12. Experimentally, these molecules exhibit low CMCs in the range of 5 × 10−5 to 1 × 10−4 mol/L. The SDK coarse grained force field19−21 was used in this work. In this model, one CG bead typically represents three heavy atoms with corresponding H atoms. As shown in Figure 1, a 10-bead model (1 OA, 5 EO, 3 CM and 1 CT2)

(1)

where the first term on the right corresponds to potential of a reference state. Using the letter S to denote the surface and B for bulk solution, the equilibrium conditions specify that μSo + RT ln XS = μBo + RT ln XB

(2)

which gives the relationship between the surface concentration XS and the bulk concentration XB. However, eq 2 is not very useful for computation, as the reference state potentials (μoS and μoB) are not readily available. The chemical potential can also be divided into ideal and excess contributions:

μ = μid + μex

(3)

The ideal contribution can be calculated analytically from statistical mechanics theory, and the excess chemical potential can be determined by using various simulation techniques such as thermodynamic integration (TI),16 Widom insertion,17 or perturbation18 methods. In this work, we applied the TI method to calculate excess chemical potentials. By calculating the chemical potentials of components in different phases using eq 3, we obtained the free energy difference for transferring a surfactant between the two phases. A general problem of the TI method is that the thermodynamic state to be modeled must be stable. For a system in which the surface and bulk phases are coexistent, both phases are unstable or metastable unless the partition between them is set exactly to that at equilibrium. This problem is less pronounced for the bulk phase, as the concentration of interest is lower or slightly above the CMC, but it is serious for micellar and lamellar phases. For example, a surface represented by a bilayer model would be unstable because excess molecules on the surface would tend to enter the bulk phase. In order to calculate chemical potentials at various surface concentrations, the unstable or metastable states must be stabilized. We accomplish this by applying an external force to stabilize the surface concentrations at the partition, which was removed from the calculated chemical potentials by using an exponential reweighting technique. Since the concentration difference is exponentially proportional to the free energy difference, it was imperative to achieve high accuracy in our calculations. Given that the simulation models are consistent throughout, the accuracy of the simulations is primarily determined by the underlying force field quality and the efficiency of sampling. The SDK force field was validated and optimized by adding hydration free energy (HFE) data to those used in the original development. Generally speaking, sampling efficiency is more problematic for the bulk phase due to the low concentration of surfactant. However, the solution at low concentration (∼10−5 mol/L) may well be treated as an ideal solution. A systematical way to evaluate the nonideal contribution for the bulk phase was developed in this work. By calculating the chemical potentials of both surfactant and water molecules carefully, we obtained chemical potential differences and predicted the surface-bulk equilibria in close agreement with experimental phase equilibria data.

Figure 1. CG models for C12EO5 and water. Different colors and symbols represent different CG beads.

represents C12EO5, and one W bead represents three water molecules. The valence terms for bond and angle interactions are described by harmonic potentials: Angle

Bond

Uvalence =

∑ k b(ri − r0)2 + ∑ i

ka(θj − θ0)2 (4)

j

where kb and ka are force constants, and r0 and θ0 are the distance and angle at energy minimum. Nonbonded forces are computed using a Lennard-Jones-9−6 (LJ-9−6) potential for all pair interactions except water: ULJ9 − 6(r ) =

9 ⎛ σ ⎞6 ⎫ 27 ⎧⎛⎜ σ ⎞⎟ ε⎨ −⎜ ⎟⎬ ⎝r⎠ ⎭ 4 ⎩⎝ r ⎠

(5)

A LJ-12−4 potential is used to describe pair interactions with water: ULJ12 − 4(r ) =

3 3 2

⎧⎛ σ ⎞12 ⎛ σ ⎞4 ⎫ ε⎨⎜ ⎟ − ⎜ ⎟ ⎬ ⎝r⎠ ⎭ ⎩⎝ r ⎠

(6)

Electrostatic interactions are included in the VDW parameters, therefore charge interactions are not considered explicitly. Chemical potential differences were calculated between surfactant molecules on the surface of a bilayer and those in bulk solution. The surfactant bilayer was modeled with a rectangular box of dimensions 5 × 5 × 40 nm3. There are two surfaces in the simulation box, separated by a 16 nm water layer, with the bilayer surface normal aligned with the long edge of the box. The surfactant molecules occupied about 2 nm of space on each surface. There was a 10 nm layer of vacuum space on top of each surface. Various numbers of surfactant molecules (4−60) were placed on each side, representing different surface concentrations. For surfactants in bulk solution, an infinitely dilute solution was simulated by placing one molecule in the middle of a cubic box with size of 5 × 5 × 5 nm3. The box was large enough to avoid any interactions between the surfactant and its periodic images. The box was filled with water molecules and its density dynamically determined. B

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

where the subscript Γ can be either the surface phase S or the bulk B. The subscripts “sur” and “wat” represent surfactant and water, and n indicates the number of water molecules corresponding to the same volume of one surfactant molecule. At equilibrium, the free energy change on both sides should be equal, therefore μS,sur − nμS,wat = μB,sur − nμB,wat (10)

NVT MD simulations were performed for the surface model described above and NPT MD simulations were performed for the bulk solution model. The LJ interactions were calculated with truncation at 1.5 nm without tail correction, which is the same condition used in the parametrization of the force field.19−21 The Nose−Hoover thermostat/barostat22,23 was used to control the temperature at 298 K, and the pressure at 1 atm. The time step used was 10 fs (fs). The initial configurations were generated by randomly placing the molecules in specified regions of the simulation boxes and equilibrated for at least 20 ns before any TI calculations were performed. The simulations were carried out using the LAMMPS24 software package. The chemical potentials were calculated using eq 3. The ideal contribution, μid on the surface and in the bulk can be written as25,26 μSid = RT ln ρS Λ3

(7)

μBid = RT ln ρB Λ3

(8)

Decomposing eq 10 into its excess and ideal parts of contribution derives the following equation: μSex + RT ln

μex S

XS XB = μBex + RT ln 1 − XS 1 − XB

(11)

μex B

Here and represent the chemical potential differences for the process of replacing the same volume of water by one surfactant molecule on the surface or in the bulk, and XS and XB are mole fractions of surfactant in the surface region and in the bulk region, respectively. XS and XB were calculated analytically depending on the concentration, and μSex and μBex were calculated by decoupling one surfactant molecule from the surface or the bulk using the TI method:

where ρS and ρB are number density of particles in the surface and bulk regions, Λ denotes the de Broglie wavelength. Both potentials were directly calculated for given number of molecules and surface area or volume. Since the surface and the bulk were modeled separately by two uncoupled simulation boxes, the free energy changes of moving one surfactant from the surface to the bulk or vice versa must be evaluated by considering both surfactant and water contributions. This is illustrated in Figure 2, which shows a

μex =

⎛ ∂F ⎞ ⎜ ⎟ = ⎝ ∂N ⎠ VT

∫0

1

∂F(λ) dλ = ∂λ

∫0

1

∂H ∂λ

dλ λ

(12)

Here, the bracket represents the ensemble average. The coupling factor λ = 0 stands for the initial state and λ = 1 stands for the end state in which the target molecule is decoupled, as the Hamiltonian is H(λ) = H0 + λ(H1 − H0)

(13)

27

The soft core model was utilized in order to avoid large fluctuations when the molecule was close to completely vanishing because of the close proximity between particles. Therefore, in this case, the real distance r was modified as r0 = (ασ06λ p + r 6)1/6

(14)

r1 = (ασ16(1 − λ) p + r 6)1/6

(15)

where α, and p are set parameters, and σ0, σ1 are equal to the VDW radii of the potential functions H0, H1 or an input parameter σ when the VDW radius is smaller than σ. Then eq 13 transforms into H(λ) = H0(r0) + λ(H1(r1) − H0(r0))

(16)

The parameters were chosen as α = 0.5, p = 1, σ = 3.0 Å. The λ points were chosen with an interval of 0.05 from 0 to 1. For each λ point, the equilibration time was 20 ns to ensure the system was fully relaxed, and the data collection time was 10 ns. The standard errors σsim(λ) at each decoupling point λ were estimated by using block average.29 The error in the calculated free energy difference was estimated by the integral 28

Figure 2. Moving of one surfactant molecule from the “surface” simulation box to the “bulk” box requires moving an equivalent number of water molecules in the opposite direction.

σΔF = (

1

2 σsim (λ) dλ)1/2

(17)

In the surface model, a harmonic potential was applied to the surfactant heads to keep the molecules on the surface (shown in Figure 3). This constrain was necessary for preventing the surfactant from leaving the surface due to either high surface concentration or vanishing interactions when the coupling parameter λ approached a value of 1. The minimum of the harmonic potential was set at the average position of the

process of moving one surfactant from the surface (left) to the bulk (right). The process consists of moving a surfactant from the surface box to the bulk box, and moving the same volume of water molecules in the opposite direction. Therefore, the free energy change for phase Γ can be written as ΔG Γ = μ Γ,sur − nμ Γ,wat

∫0

(9) C

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

agreement but are much closer to the experimental values, which is not too surprising, as the alkane parameters were optimized by fitting to a series of thermodynamic properties in the original SDK development. The VDW interaction parameters between the EO and OA beads with the water bead were determined mainly based on structural information (such as density profile, bilayer thickness, etc.)19 and this is not sufficient to accurately determine thermodynamic properties. Therefore, we chose to optimize the VDW interaction parameters between EO, OA, and water beads using the experimental HFE data. We found that a scaling ε by a factor of 0.86 (in eq 6) yields good agreement with the experimental data as shown in Table 1. To validate this modification, we also calculated the surface tensions of a EO surfactant layer-water solution at an experimental surface concentration for validation. The surface tension was calculated from ensemble averaged pressure tensors using the following equation:

Figure 3. Scheme of the bilayer simulation box. The dashed lines represent the external harmonic potential acting on surfactants.

surfactant heads extracted from a low concentration simulation. The spring constant for the harmonic potential was tested and chosen to be 1.0 kcal/Å2. The calculated free energy values were reweighted to remove the effect of the imposed harmonic potentials by Boltzmann exponential reweighting.

γ=

III. RESULTS AND DISCUSSION Validation of the Simulation Models. The SDK force field parameters were validated first for the purpose of this work. Essentially, it is the hydrophilicity and hydrophobicity of the surfactant that determines its partition between the surface and bulk phases. Therefore, we carried out calculations to determine the hydration free energies (HFE) of EO surfactants in addition to the other properties (surface tension and surface structures) that were used in the original SDK development. Using the TI method, the HFEs of a series of EO surfactants were calculated by removing one surfactant molecule from a cubic water box of 5 × 5 × 5 nm3. Results are shown in Table 1

model

original 7.8 8.1 8.5 9.2 9.4 9.8 10.1 10.7 11.3 −36.7 −50.3 −65.5 −80.5 −98.8 −112.2 −128.1

modified

exp.

−25.9 −34.3 −45.7 −53.6 −65.5 −75.8 −86

9 9.6 10.6 10.8 12 12.2 12.8 10.8 14.2 −25.2 −35.5 −45.7 −56.1 −66.3 −76.4 −86.9

(18)

Since EO surfactants have quite large heads and exhibit high solubility in water, the monolayer structure is unstable, especially at relatively high surface concentrations; therefore, the surfactants were constrained on the surface using an external harmonic potential, the same used in the TI calculations later. The calculated data are given in Table 2, which are consistent with the experimental data and follow the same trend. Table 2. Calculated Surface Tensions (mN/m) Using Modified SDK Force Fields, with the Experimental Data for Comparisona

Table 1. Calculated HFE Data (kJ/mol) for Alkanes and EO Surfactants Using Original and Modified SDK Force Fields, with the Experimental Data for Comparison C4H10 C5H12 C6H14 C7H16 C8H18 C9H20 C10H22 C11H24 C12H26 C12EO C12EO2 C12EO3 C12EO4 C12EO5 C12EO6 C12EO7

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

model C12EO3 C12EO4 C12EO5 C12EO6 C12EO8 a

surface tension 26.8 29.4 31.9 34.0 38.8

(6.1) (6.3) (6.4) (6.8) (5.3)

exp. 27.6 28.6 30.0 33.5 36.1

The calculated uncertainties are given in the parentheses.

The modification to the SDK parameters does not violate other predictions originally made by the authors. The number density for each of the CG beads along the C12EO2/water bilayer normal is almost identical to that reported in the literature19 (shown in Figure S1). In addition, the equilibrated cross-sectional area of the C12EO2 lamellar system in the NPT ensemble is calculated to be 0.294 nm2, which is close to the experimental value of 0.300 nm2 and the value calculated by the original SDK parameter set.19 The surface structure of the lamellar surfactant system provides another aspect to evaluate the force field parameters. Figure 4 shows the two distributions of chain tilt angles calculated for a C12EO5/water system at the experimental surface concentration (50 Å2 per molecule). A comparison is shown between a coarse-grained force field (CGFF) and all atom force field (AAFF)9 representations of the system. The chain tilt angle is defined as the angle between the tail axis and the normal to the surface. Inspection of the curves reveals that the CGFF profile is smoother and slightly overestimated in the small angle region (the molecules are more up-right). However,

with comparison to experimental data30−32 derived from air/ water partition coefficients. The calculated values for the HFE of EO surfactants using the original SDK parameters are systematically about 40% lower than the experimental data, which means the surfactant model is more hydrophilic than reality. To isolate the problem in parameters, we calculated the HFE for a series of alkanes, the results are not in perfect D

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

concentrations, the external harmonic potentials has little effect on the surface structures. At high concentration, there is a shoulder peak in the unconstrained system, indicating some surfactants are not on the surface, which can be avoided by using the aforementioned constraint. Computation of Chemical Potentials. The chemical potentials at both sides of eq 11 were calculated separately. The surface concentration is expressly included in the ideal contribution and implicitly included in the excess contribution, e.g., the excess contributions at different surface concentrations were calculated for given surface concentrations. Table 4 lists Figure 4. Normalized tilt angle distributions (a) calculated using CGFF and (b) calculated using AAFF.

Table 4. Calculated Excess Chemical Potential (kJ/mol) and Its Energetic and Entropic Components (kJ/mol) of C12EO5 on Surface as a Function of Surface Concentrations Γ (×1010 mol/cm2)

the figure also shows the most probable tilt angles in both the CGFF and AAFF profiles are about 40°, with average values of about 52°, which agrees well with the experimental value2 of ca. 50°. The average thicknesses of surfactant tails and heads were calculated by fitting the equilibrium density profile to a Gaussian function:33 ⎛ −4(Z − Z )2 ⎞ 0 ⎟ ρ = ρ0 exp⎜ 2 σ ⎝ ⎠

(19)

where Z0 denotes the equilibrium position and σ represents the average width of the tail or head layer. Table 3 lists the Table 3. Average Thickness (nm) of the Tail and Head Groups of Surfactants tail head

CG

AA

Exp.

1.03 1.35

0.90 1.35

∼ 1.2 ∼ 1.3

Γ

μSS

μSS(U)

μSS(−TS)

0.27 1.33 2.39 2.66 2.92 3.19 3.45 3.72 3.98

−115.4 (1.5) −110.0 (2.4) −103.8 (2.9) −102.1 (3.6) −100.2 (3.9) −98.4 (4.0) −96.8 (3.0) −94.7 (4.5) −93.0 (5.0)

−147.2 −237 −316.7 −311.5 −320.6 −335.3 −351.2 −354.1 −373.8

31.8 127.0 212.9 209.4 220.4 236.9 254.4 259.4 280.8

the calculated values for C12EO5 (similar results were obtained for other surfactants as shown in Figure S2). In the evaluation of the ideal contribution, the surface thickness was estimated using the surfactant molecule length37 listed in Table 5. The Table 5. Surfactant Length (nm) and Calculated Excess Chemical Potential (kJ/mol) of Surfactant in Infinite Dilute Solution

calculated thicknesses compared with the results obtained using simulations using an AAFF34 and also experimental data.35,36 It can be seen that the CGFF results are consistent with both the experimental data and the AAFF results. The impact of the external harmonic potential was investigated by comparing density profiles of the EO molecules at different packing densities. Figure 5 shows the density profiles of the surfactant heads with different surface concentrations (with 4, 36, 68 surfactants on each side). The right side of this figure is for free systems without an external harmonic potential, and the left is for the constrained systems. Comparing the distribution profiles of the low and medium

surfactant

length

μex

C12EO3 C12EO5 C12EO6 C12EO8

2.4 2.9 3.1 3.5

−43.3 −65.5 −75.8 −97.1

calculated surface chemical potential (μS) increases monotonically with the surface concentration because the high surface concentrations were maintained by the imposed external harmonic potential. In experiment, the excess molecules would migrate to the bulk phase and the surface concentration is stabilized at the saturation point of 2.65 × 10−10 mol/cm2. The free energy change can be decomposed into the energetic and entropic contributions: ΔF = ΔU − T ΔS

(20)

Table 4 lists the estimated enthalpic contribution and entropic contribution of the surface chemical potential at various surface concentrations. The decomposition shows that, within the concentration range considered, adding surfactants to the surface causes both enthalpy and entropy to decrease. Therefore, increasing the surface concentration is an enthalpically favored and entropically disfavored process. The bulk phase chemical potentials were also calculated using eq 11. Again, the ideal contribution is a function of bulk concentration (ρB) explicitly, and the excess contribution is

Figure 5. Density distribution profiles of surfactant heads at low, middle, and high surface concentrations for (a) free system and (b) constrained system. E

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B implicitly dependent on concentration. However, the excess contribution can be treated as a constant because of the low concentration (∼10−5 mol/L). The simulation box has to be at least 38 × 38 × 38 nm3 (nearly 144 times greater than the model used in this work) in order to put two C12EO5 molecules in one simulation box at the CMC (6 × 10−5 mol/L). In this scenario, the average distance (about 60 nm) between two C12EO5 surfactants at the CMC is far beyond the simulation cutoff value of 1.5 nm. We estimated the nonideal contribution for C12EO5 based on a lattice model and potential of mean force (PMF) calculation.38 Assuming only monomers, dimers and trimers exist in the solution, we set the lattice grid to be 25 × 25 × 25 Å3 and each lattice site is either empty or occupied by one monomer or one of the small oligomers (of 2 or 3 monomers). This is based on the fact that the PMF curves show that the attractive forces between two or three molecules vanish beyond 25 Å (shown in Figure 6). Given that the grid size is fixed, the

Figure 7. (Left) surface chemical potential (μSS) as a function of surface concentration. (Right) bulk chemical potential (μSB) and the chemical potentials derived from the experimental equilibrium data as shown in red line.

experimental μexp B curve for comparison. It can be seen that the calculated data are very close to the experimental data, with a small deviation of about 1 kJ/mol. The causes of the small deviation can be the approximation of surface thickness estimation, the force field quality, or the simulation model could be responsible. From the comparison of chemical potential curves between μS and μB, the equilibrium surface-bulk concentrations can be obtained. Despite the fact that the partition between the two phases is exponentially related to the free energy difference, reasonable agreements are obtained for the four surfactants studied. The results shown in Figure 8 agree well with the experimental data.3

Figure 6. PMF curves of C12EO5 dimer and trimer in water solution.

numbers of molecules and lattice sites are variables of the bulk concentration. The binding energies calculated using energy minimization for the dimers and trimers are ∼8.1 and ∼15.6 kJ/mol, respectively. With these data, the Boltzmann weighted probabilities (ρ2,ρ3) of forming dimers and trimers at 10−5 mol/L concentration are ∼2.4 × 10−3 and ∼4.6 × 10−6, respectively. The results are similar using either a lattice Monte Carlo simulation or analytical calculation (permutationcombination with Boltzmann weighting factor) as shown in Figure S3. The additional free energy relative to that of the monomer is then 3

Δμadd =

∑ ρi Δμi i=2

Figure 8. Predicted equilibrium concentration between surface and bulk of EO surfactant systems compared with experimental data.

(21)

where Δμi is the free binding energy of an oligomer of size i, which is ∼7.4 kJ/mol for a dimer, and ∼13.1 kJ/mol for a trimer, estimated from the PMF curves shown in Figure 6. Using eq 21, the nonideal contribution to the chemical potential (Δμadd) for C12EO5 is only about 1.8 × 10−2 kJ/mol, negligible in comparison with the ideal contribution of −65.5 kJ/mol. The excess chemical potentials for other EO surfactants were calculated in the same way, and the results are listed in Table 5. Prediction of Equilibrium. At equilibrium, eq 11 must be satisfied to keep the chemical potentials of the two phases identical. Figure 7 shows the surface chemical potential μS on the left and bulk chemical potential μB on the right. Using the experimental surface-bulk concentration diagram,3 we draw the

IV. SUMMARY In this work, we demonstrated a method to predict surface-bulk equilibria by calculating the chemical potentials of participating components using uncoupled simulation boxes to represent the surface and bulk phases separately. In comparison with direct simulation methods, this approach is more efficient and reliable, and especially suitable for systems with very low solution concentration. The method can be applied to any phase equilibria. However, it is important to note that in order to link the uncoupled simulations together, evaluations of chemical F

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(3) Lu, J. R.; Su, T. J.; Li, Z. X.; Thomas, R. K.; Staples, E. J.; Tucker, I.; Penfold, J. Structure of Monolayers of Monododecyl Dodecaethylene Glycol at the Air−Water Interface Studied by Neutron Reflection. J. Phys. Chem. B 1997, 101 (49), 10332−10339. (4) Meyers, D. Surfactant Science and Technology, Emulsion. WileyVCH: Weinheim, Germany, 1946. (5) 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. (6) 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 (5), 1393−1402. (7) Wijmans, C. M.; Linse, P. Surfactant self-assembly at a hydrophilic surface. A Monte Carlo simulation study. J. Phys. Chem. 1996, 100 (30), 12583−12591. (8) 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 (40), 8493− 8501. (9) 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 (2), 114−118. (10) Bandyopadhyay, S.; Chanda, J. Monolayer of Monododecyl Diethylene Glycol Surfactants Adsorbed at the Air/Water Interface: A Molecular Dynamics Study. Langmuir 2003, 19 (24), 10443−10448. (11) Chanda, J.; Bandyopadhyay, S. Molecular dynamics study of a surfactant monolayer adsorbed at the air/water interface. J. Chem. Theory Comput. 2005, 1 (5), 963−971. (12) 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 (1), 471−479. (13) Yang, C.; Sun, H. Surface−Bulk Partition of Surfactants Predicted by Molecular Dynamics Simulations. J. Phys. Chem. B 2014, 118 (36), 10695−10703. (14) Panagiotopoulos, A. Z. Direct Determination of Fluid Phase Equilibria by Simulation in the Gibbs Ensemble: A Review. Mol. Simul. 1992, 9 (1), 1−23. (15) Panagiotopoulos, A. Z. Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble. Mol. Phys. 1987, 61 (4), 813−826. (16) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications; Academic Press: San Diego, CA, 2001. (17) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39 (11), 2808. (18) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 (8), 1420−1426. (19) Shinoda, W.; DeVane, R.; Klein, M. L. Multi-property fitting and parameterization of a coarse grained model for aqueous surfactants. Mol. Simul. 2007, 33 (1−2), 27−36. (20) Shinoda, W.; DeVane, R.; Klein, M. L. Coarse-grained molecular modeling of non-ionic surfactant self-assembly. Soft Matter 2008, 4 (12), 2454. (21) 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 (20), 6836−6849. (22) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31 (3), 1695− 1697. (23) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511. (24) Plimpton, S.; Crozier, P.; Thompson, A. LAMMPS: Large-Scale Atomic/Molecular Massively Parallel Simulator; Sandia National Laboratories: Albuquerque, MN, 2007.

differences for all species must be conducted. In this work, the calculations must include both surfactant and water molecules. For the TI calculations, an external force was applied to maintain otherwise unstable or metastable states. The external forces can be removed by using an exponential reweighting technique. The force field quality is doubtlessly one of the most critical factors that determines the quality of the prediction. We found that the SDK parameters required a minor correction to describe the hydrophilicity of surfactant by fitting to experimental HFE data. To calculate the excess chemical potential of nonideal solutions, we evaluated the probability of forming small oligomers (up to a trimer) using statistical analyses and Monte Carlo simulations based on a lattice model, and the free energies of formation of the small oligomers using PMF calculations. It is confirmed that the nonionic surfactant solution near the CMC (∼10−5 mol/L) can be treated as an ideal solution. However, this technique is general, and can be applied to treat nonideal solutions. The calculated chemical potential curves are very close to that derived from experimental data. The little discrepancy that exists may come from approximations made during the calculations, or the inaccuracy of the force field and simulation model. Although the exponential dependence on free energy difference makes prediction of concentration equilibrium extremely difficult, based on the calculated chemical potential curve, the predicted equilibria of the surface-bulk concentration is in reasonable agreement with experiments.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.5b10239. Validation of the modified force field, calculated chemical potentials for C12EO3, C12EO6, and C12EO8, probability of dimer formation using analytical calculation and MC simulation, and optimized SDK force field parameters (PDF)

■

AUTHOR INFORMATION

Corresponding Author

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

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work has been supported by the National Natural Science Foundation of China, Numbers 21403138, 21173146 and 21473112, the National Basic Research Program of China (973 Program), Number 2014CB239702, and computational resources from the Center for High Performance Computing at Shanghai Jiao Tong University.

■

REFERENCES

(1) Rosen, M. J.; Kunjappu, J. T. Surfactants and Interfacial Phenomena; Wiley: Hoboken, NJ, 2012. (2) Lu, J. R.; Li, Z. X.; Thomas, R. K.; Binks, B. P.; Crichton, D.; Fletcher, P. D. I.; McNab, J. R.; Penfold, J. The Structure of Monododecyl Pentaethylene Glycol Monolayers with and without Added Dodecane at the Air/Solution Interface: A Neutron Reflection Study. J. Phys. Chem. B 1998, 102 (30), 5785−5793. G

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (25) Reiss, H.; Frisch, H. L.; Lebowitz, J. L. Statistical Mechanics of Rigid Spheres. J. Chem. Phys. 1959, 31 (2), 369. (26) Helfand, E.; Frisch, H. L.; Lebowitz, J. L. Theory of the Twoand One-Dimensional Rigid Sphere Fluids. J. Chem. Phys. 1961, 34 (3), 1037. (27) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. GROMACS: fast, flexible, and free. J. Comput. Chem. 2005, 26 (16), 1701−1718. (28) Hess, B.; van der Vegt, N. F. Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models. J. Phys. Chem. B 2006, 110 (35), 17616−17626. (29) Lawrenz, M.; Baron, R.; McCammon, J. A. IndependentTrajectories Thermodynamic-Integration Free-Energy Changes for Biomolecular Systems: Determinants of H5N1 Avian Influenza Virus Neuraminidase Inhibition by Peramivir. J. Chem. Theory Comput. 2009, 5 (4), 1106−1116. (30) CSID:17269, http://www.chemspider.com/Chemical-Structure. 17269.html (accessed 15:24, Aug 13, 2015). (31) CSID:68929, http://www.chemspider.com/Chemical-Structure. 68929.html (accessed 15:22, Aug 13, 2015). (32) CSID:71267, http://www.chemspider.com/Chemical-Structure. 71267.html (accessed 15:19, Aug 13, 2015). (33) Domínguez, H. Computer Simulations of Surfactant Mixtures at the Liquid/Liquid Interface. J. Phys. Chem. B 2002, 106 (23), 5915− 5924. (34) Li, X.; Li, Y.; Wu, Z.; Xie, Z.; Fan, K. Molecular Dynamics Simulations of Nonionic Surfactant at the Air/water Interface. Acta Chim. Sin. 2011, 69 (19), 2235−2240. (35) Lu, J.; Li, Z.; Thomas, R.; Staples, E.; Thompson, L.; Tucker, I.; Penfold, J. Neutron Reflection from a Layer of Monododecyl Octaethylene Glycol Adsorbed at the Air-Liquid Interface: The Structure of the Layer and the Effects of Temperature. J. Phys. Chem. 1994, 98 (26), 6559−6567. (36) Lu, J.; Li, Z.; Su, T.; Thomas, R.; Penfold, J. Structure of adsorbed layers of ethylene glycol monododecyl ether surfactants with one, two, and four ethylene oxide groups, as determined by neutron reflection. Langmuir 1993, 9 (9), 2408−2416. (37) 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 (2), 172−185. (38) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23 (2), 187−199.

H

DOI: 10.1021/acs.jpcb.5b10239 J. Phys. Chem. B XXXX, XXX, XXX−XXX