Combined Molecular Dynamics Simulation ... - ACS Publications

Jul 27, 2017 - Corporate Strategic Research, ExxonMobil Research & Engineering Company, 1545 Route 22 East, Annandale, New Jersey 08801,...
1 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/Langmuir

Combined Molecular Dynamics Simulation−MolecularThermodynamic Theory Framework for Predicting Surface Tensions Vishnu Sresht,† Eric P. Lewandowski,‡ Daniel Blankschtein,*,† and Arben Jusufi*,‡ †

Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139, United States ‡ Corporate Strategic Research, ExxonMobil Research & Engineering Company, 1545 Route 22 East, Annandale, New Jersey 08801, United States S Supporting Information *

ABSTRACT: A molecular modeling approach is presented with a focus on quantitative predictions of the surface tension of aqueous surfactant solutions. The approach combines classical Molecular Dynamics (MD) simulations with a molecular-thermodynamic theory (MTT) [Y. J. Nikas, S. Puvvada, D. Blankschtein, Langmuir 1992, 8, 2680]. The MD component is used to calculate thermodynamic and molecular parameters that are needed in the MTT model to determine the surface tension isotherm. The MD/MTT approach provides the important link between the surfactant bulk concentration, the experimental control parameter, and the surfactant surface concentration, the MD control parameter. We demonstrate the capability of the MD/MTT modeling approach on nonionic alkyl polyethylene glycol surfactants at the air−water interface and observe reasonable agreement of the predicted surface tensions and the experimental surface tension data over a wide range of surfactant concentrations below the critical micelle concentration. Our modeling approach can be extended to ionic surfactants and their mixtures with both ionic and nonionic surfactants at liquid− liquid interfaces.

■

INTRODUCTION Although the fundamental mechanisms underlying the adsorption of surfactants at interfaces and the subsequent alteration of interfacial properties are well understood, making quantitative predictions of interfacial tensions (IFTs) in surfactant systems remains an active research area in interfacial science. Determining the relationship between the IFT and the quantities and the structures of surfactant additives in different temperature, pressure, salinity, and pH regimes is critical to formulation design in many industry sectors, ranging from consumer chemicals to the oil and gas industry. Advances in force-field parametrization have allowed Molecular Dynamics (MD) simulations to play the role of computational microscopes, enabling atomistic-level insights into the interactions between surfactant moieties and their surroundings at fluid interfaces. Several simulation studies1−7 have examined the effect of surfactants on the interface structure and the IFT. These simulations assume an a priori knowledge of the surface concentration Γ of the surfactant molecules, permitting interfacial properties such as the IFT (γ) to be studied as a function of the surfactant concentration, that is, γ = γ(Γ). Unfortunately, this implies that predicting γ given the bulk concentration c of the surfactant requires the determination of the adsorption isotherm Γ(c). Because conventional MD simulations can only access length- and time-scales of tens of nanometers and nanoseconds, respec© 2017 American Chemical Society

tively, it remains a challenge to utilize MD simulations to model both the dynamics of surfactant molecules at the interface and the adsorption/desorption equilibrium between the bulk solution and the interface that controls Γ(c). In principle, one could use Grand-Canonical Monte Carlo simulations, where the quantity known a priori is the chemical potential μ of the surfactants rather than Γ. However, this method is computationally unfeasible for liquid systems modeled on atomistic length-scales, particularly for moderate- to highmolecular weight surfactants. Recently, a link between the bulk concentration and the surface concentration within an MD approach was generated through free-energy calculations at different surface coverages.7,8 This approach is generally exact but requires numerous free-energy computations. Adsorption isotherms and the structural features of the interface may also be determined using experimental methods such as neutron reflectivity,9,10 high-energy X-ray reflectivity,11 or EPR techniques.12 However, these methods do not scale well enough to enable the high-throughput screening of large libraries of surfactants necessary for the rational design of surfactant formulations. Alternatively, Γ(c) may be calculated from thermodynamic adsorption isotherms, such as the Received: March 31, 2017 Revised: July 25, 2017 Published: July 27, 2017 8319

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

difference between γ0, the IFT of a pure air−water system, and γ, the IFT due to the adsorption of the surfactants at the air−water interface; i.e., Π = γ0 − γ. As in the original paper by Nikas et al., we relate this surface pressure to the surface coverage of the surfactant adsorbed at the interface and to the temperature T of this interface through a twodimensional equation of state Π = Π(Γ, T). We construct this equation of state based on the following assumptions: (i) at low surfactant concentrations, there is a single monolayer of surfactant molecules adsorbed at the air−water interface, as shown in Figure 1,

Langmuir,13 Frumkin,14 Fainerman,15 or other Gibbs-based adsorption models that also include intermolecular interactions. However, the value of Γ obtained is strongly sensitive to the isotherm chosen and the applicability of the underlying assumptions to the system being studied.16−18 Furthermore, these isotherms are no more favorable than experimental characterization for high-throughput formulation design because they often (a) contain empirical parameters that must be determined by fitting to experimental data, and (b) cannot be applied in a thermodynamically consistent manner to surfactant mixtures.19 Molecular-thermodynamic theory (MTT) provides yet another framework for the calculation of adsorption isotherms. MTT models are predictive in capturing nonideal free-energy contributions and have also been extended beyond pure adsorption predictions, e.g., for capturing micellization and microemulsion formation.20−27 However, as with the adsorption isotherms discussed above, the utilization of MTT models requires the determination of numerous parameters, if not by experiment, then via empirical models. For example, one of the fundamental parameters in MTT is the driving force for surfactant adsorption, a parameter that, at the air−water interface, is mainly governed by the transfer free energy of the hydrophobic surfactant segment. There are reliable estimates of this free-energy contribution for aliphatic chains based on various models,20,28,29 but these estimates may not be applicable to surfactants with nonaliphatic chain architectures. In this paper, we demonstrate that the challenges addressed above can be overcome by computing the parameters inherent in MTT through MD simulations. We choose an MTT approach originally developed by Nikas, Puvvada, and Blankschtein for mixtures of nonionic surfactant solutions at air−water interfaces.23 This MTT approach has been further refined to also model ionic and zwitterionic surfactants, as well as their mixtures, at air−water and oil−water interfaces.24,26,30 The original MTT model requires a single experimental input quantity, i.e., the free energy of adsorption of a surfactant that is determined via IFT measurements. Instead, here, we predict all the parameters required by MTT without any experimental inputs by utilizing MD simulations, and then use MTT to calculate IFT values for a range of bulk concentrations below the critical micelle concentration (CMC). We showcase this hybrid MD/MTT approach for two nonionic alkyl surfactants derived from polyethylene glycol (PEG) in an air−water system, and we validate our predictions against experimental data. By combining MD with MTT, we provide a fully predictive model for the IFT as a function of the bulk solution concentration that would allow the high-throughput in-silico screening of surfactant formulations. The remainder of the manuscript is organized as follows: We first recap the MTT framework and provide details of the MD model settings, including the methods that are used for the calculation of the MTT input parameters. We then present adsorption isotherms and surface tension predictions, including comparing these with experimental values taken from the literature for alkyl PEG surfactants. We finally summarize the results and provide an outlook for future work in this field.

■

Figure 1. Modeling the monolayer of surfactant molecules adsorbed at the air−water interface as a 2D gas. An atomistic snapshot of a simulation containing C10E4 (n-decyl tetra (ethylene oxide)) molecules in SPC/E water is shown on the left. Water molecules are shown as red lines, and carbon, hydrogen, and oxygen atoms are shown as cyan, white, and red spheres, respectively. The two figures on the right-hand side show the side and top views of the air−water interface. Water is shown in red, and each surfactant molecule is modeled as a hard disk (shown in green) with short-range intermolecular attractions.

and (ii) at low surfactant concentrations, individual surfactant molecules in this monolayer may be modeled as hard disks (with radius r) with short-ranged attractive interactions. The consideration of exclusively nonionic surfactants simplifies the equation of state by eliminating the need to consider electrostatic interactions between the charged surfactant moieties, as well as counterion binding. The underlying picture for this MTT is that the surfactant headgroups are positioned at the Gibbs dividing surface at the air−water interface, whereas the tails stick out into the air phase. The latter assumption may not be correct in the dilute regime, where tails may exhibit configurations with increased contact with the water phase, which is energetically favorable compared to the contact with the gas phase without any enthalpic benefit. However, in the dilute regime those interactions are weak and do not significantly alter the overall surface pressure, as our results show. Under these assumptions, it follows that23

⎡ ⎤ γ −γ 1 Π = 0 = Γ⎢ + BΓ 2 2⎥ kBT kBT ⎣ (1 − η) ⎦

(1)

where η = Γa is the packing fraction of the surfactants adsorbed at the interface, a = πr2 is the excluded area of the surfactants, and B is a pairwise coefficient that accounts for the short-ranged van der Waals attraction between the adsorbed surfactants. Note that B is a function of the temperature and the pressure of the air−water interface. However, in the following treatment, the temperature and the pressure of the interface are assumed to be constant. Furthermore, since B is a measure of attractive intermolecular forces, it is expected that B ≤ 0. Finally, in this MTT, the effects of the interactions between the surfactant molecules and water (such as the hydrophobic effect31) are accounted for implicitly through the parameters r and B. Although we

METHODS

Molecular-Thermodynamic Theory (MTT). The MTT is applied to a single species of nonionic surfactant in an air−water system, and it is a slightly modified version of the theory from the original paper by Nikas et al.23 We begin with the surface pressure Π defined as the 8320

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

(related to c) of the nonionic surfactant, and (2) eq 1 to relate the IFT γ to Γ. Therefore, by simultaneously solving eqs 1 and 7, one can calculate the surface tension isotherm, γ(c), for a single-surfactant system. However, before these two equations can be used, one must determine the three parameters in these equations: (i) the excluded area a (or equivalently the radius r) of the surfactant headgroups, (ii) the pairwise interaction coefficient B, and (iii) Δμ0, the difference between the reference-state chemical potential of a surfactant molecule at infinite dilution at the surface and that of a surfactant molecule at infinite dilution in the bulk solution phase. As demonstrated in the next two subsections, we will determine these three parameters using MD simulations. Determining the Parameters r and B. The two parameters r and B are associated with the equation of state for a single nonionic surfactant adsorbed at the air−water interface (see eq 1), and in the MD/MTT framework, these parameters are estimated using computer simulations. Two methods can be used to estimate r and B. The first method (“Method 1”) is based on the recognition that the 2D equation of state derived above (eq 1) is based on a model for pairwise intermolecular interactions. Therefore, the parameters r and B are explicitly related to the intermolecular interactions between two surfactant molecules. This suggests that one can obtain r and B by estimating the intermolecular interactions between two adsorbed surfactant molecules using MD simulations. The second method (“Method 2”) is based on the recognition that the 2D equation of state represents a relationship between the surface pressure Π and the adsorption Γ of the surfactant molecules at the air−water interface. Therefore, r and B can be estimated by fitting the equation of state in eq 1 to the predicted equation of state obtained using MD simulations at different surfactant adsorption values. Note that Method 2 mimics the broadly used practice of deducing the equation of state parameters by fitting it to the experimentally deduced equation of state. Here, however, we use instead the simulation-predicted equation of state, thereby circumventing the need for experimental data. Below, we present and compare the results of both methods when utilized to model the alkyl nonionic surfactants C10E4 (n-decyl tetra (ethylene oxide)) and C12E6 (n-dodecyl hexa (ethylene oxide)). Method 1: r and B Parametrization Based on Potential of Mean Force (PMF) Calculations. The use of Method 1 to estimate r and B requires computing the interaction energy between a single pair of nonionic surfactant molecules at the air−water interface. To this end, we set up a 5800 SPC/E molecules and 2 C10E4 molecules in a rectilinear simulation box of dimensions 7 nm × 7 nm × 10 nm. Both C10E4 molecules were positioned, at random, within the bulk of the SPC/E molecules. This simulation box was then equilibrated in the NVT ensemble at 300 K with a Berendsen thermostat33 with a 0.1 ps time constant for 1 ns with a time step of 2 fs. During this equilibration, the two surfactant molecules diffused to, and equilibrated at, the air−water interface. Subsequently, an externally imposed harmonic potential (with a force constant of 1000 kJ mol−1 nm−2) was used to vary the separation between the carbon atoms residing at the ends of the surfactant heads on both surfactant molecules (this carbon atom, shown in green in Figure 2 below, is referred to hereafter as the “anchor” atom) from 0.35 to 2.50 nm. Attempts to bring the two anchor atoms closer than 0.35 nm resulted in simulation errors caused by extremely strong interatomic repulsions.

restrict ourselves to the single-component, nonionic surfactant case, eq 1 was originally formulated for surfactant mixtures.23 With an equation of state in hand, it is possible to derive an explicit expression for the chemical potential μσ of a surfactant molecule adsorbed on the air−water interface. The Gibbs Adsorption Equation states that 1 dΠ Γ

dμ σ =

(2) σ

By rewriting eq 2 in terms of deviations of μ from a reference-state value μσ,id, and of Π from its reference-state value Πid, one can derive the following equation to determine μσ for any given equation of state:23 μσ − μσ ,id =

Γ

∫Γ=0 Γ1

d (Π − Π id) dΓ dΓ

(3)

Contrary to the original Nikas et al. paper, here, we define our reference state as an air−water interface at infinite dilution, and the surface pressure corresponding to this reference state is Πid = 0 mN/ m, a value different from the reference surface pressure in the original papers Πid = 1 mN/m.23,24 Since intermolecular interactions can be ignored at these state points, both choices are legitimate. However, the difference in the selection of the reference point in the context of the MD/MTT approach is quite important. Because of the challenges associated with carrying out MD simulations in an ensemble where the surface pressure is maintained at a constant nonzero value for an arbitrary surfacant molecule, it is difficult to calculate μσ, id using MD with proper precision at a reference state defined by a surface pressure of Π0 = 1 mN/m. Recall that surface tension and pressure are calculated through stress components. In the nanometer-size systems that are typically simulated using MD, the stress components fluctuate strongly and therefore limit the precision of quantities, such as IFT, that are determined via stress computations. By utilizing the equation of state in eq 1 in eq 3, we obtain the following expression for the chemical potential of a surfactant molecule at the surface:

μσ μσ ,id 3η − 2η2 = − ln(1 − η) + + 2B Γ kBT kBT (1 − η)2

(4)

id

At our chosen reference state at infinite dilution (Π = 0 mN/m), we can calculate the reference-state chemical potential of an adsorbed surfactant molecule with a purely entropic term:

μσ ,id = μσ ,0 + kBT ln(η)

(5)

Assuming that the surfactants in the bulk do not interact with one another in the dilute surfactant concentration range, i.e., below the CMC, the corresponding surfactant bulk chemical potential is simply expressed as follows:

μb = μb,0 + kBT ln x where x =

c c + cw

(6)

denotes the bulk mole fraction of surfactants, and cw

represents the molar concentration of bulk water. The surfactant molecules may also be modeled as polymers using the Flory−Huggins solution theory.32 The results of using this model for the chemical potential of surfactant molecules in the bulk are presented in the Supporting Information. At thermodynamic equilibrium, μσ = μb, yielding the key equation for the adsorption isotherm

⎛ η ⎞ Δμ0 3η − 2η2 + ln⎜ + 2B Γ ⎟+ kBT ⎝1 − η ⎠ (1 − η)2

ln x =

σ,0

(7)

Figure 2. Schematic of a C10E4 molecule used in the MD simulations. The ethoxylated head is shown in the cyan box. The hydrocarbon tail is shown in the orange box. The anchor atom is indicated in green. The anchor atom is the heavy atom located at the Gibbs Dividing Surface between the air and water phases. See Supporting Information for details about how this atom was identified. Carbon, oxygen, and hydrogen atoms are shown in black, red, and white, respectively.

where Δμ ≡ μ − μ is the difference between the reference-state chemical potential of a surfactant molecule at infinite dilution at the interface and the reference-state chemical potential of a surfactant molecule at infinite dilution in the bulk solution phase. The MTT derivation above yields (1) eq 7 to relate the adsorption Γ at the air− water interface (recall that η is related to Γ) to the bulk mole fraction x 0

b,0

8321

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

Figure 3. Schematic of the simulations utilized to compute the PMF between two C10E4 molecules at the air−water interface. In the atomistic simulation snapshots, water molecules are shown as red lines, and carbon, oxygen, and hydrogen atoms are shown as black, red, and white spheres, respectively. The anchor carbon atoms (shown in green) in each C10E4 molecule are pulled apart using a harmonic potential, from a separation of 0.35 nm to a separation of 2.50 nm. Simulation snapshots taken at these two separations, and at 37 separations in between, are used as the initial configurations to carry out biased-sampling simulations. The top and side views of two such simulations, at 1.06 and 1.78 nm, are shown in the cyan and black boxes, respectively. The histograms corresponding to the separations between the anchor atoms during these two simulations are plotted in cyan and black, respectively, below the corresponding top views. between a pair of surfactant molecules at the air−water interface. The value of the PMF at a given separation, x, represents the work needed to pull the two molecules together from infinity to a separation of x. The negative of the gradient of the PMF at a given separation x is the force between the two molecules at that separation. The distribution of the separations between the two anchor carbon atoms obtained from these simulations is shown in the histograms in the lower half of Figure 4. The significant overlap between adjacent histograms indicates that there is a significant overlap between the intermolecular separations sampled by each of the biased-sampling simulations. This overlap is essential to the success of the WHAM algorithm. The upper half of Figure 4 shows the PMF curve obtained by applying the WHAM algorithm to the histograms. At small separations, the gradient of the PMF is very negative, and it indicates strongly repulsive interactions. This is consistent with the repulsive excluded-volume interactions that one would expect between two hard disks. However, because the C10E4 molecules are only approximately hard disks, the repulsion at small separations has a finite gradient. The hard-disk diameter of a C10E4 molecule can be estimated as the location of the minimum in the PMF, and it is approximately 0.48 nm. Note that the corresponding value for the hard-disk radius (0.24 nm) is close to the estimate of 0.3 nm obtained by Nikas et al.23 At separations ≈1.50 nm, the gradient of the PMF, in Figure 4, goes from being weakly positive to zero, suggesting that as the separation between the two C10E4 molecules increases beyond 2.25 nm, the force between the two C10E4 molecules changes from weakly attractive to nonexistent. The deduced trend in the intermolecular force versus the intermolecular separation is consistent with our model of short-ranged attractions between the hard disks. The value of the parameter B can be obtained by integrating the simulated PMF as a function of the intermolecular distance from the distance of closest approach between the two C10E4 molecules (0.48 nm) to a distance of 2.25 nm at which the PMF is zero:45

At separations above 2.50 nm, as shown in Figure 4, the intermolecular interactions were observed to be negligible. Simulation snapshots (i.e., a record of the instantaneous positions of all the atoms in the system) corresponding to 39 equispaced separations between, and including, 0.35 and 2.50 nm were selected and used as the initial configurations for a set of 39 biased-sampling simulations, consistent with the MD simulation literature.34 The two surfactant molecules in these simulations were represented using the all-atomistic OPLS-AA force field in a manner that has been previously validated by Blankschtein and co-workers.35−37 The OPLS-AA force field accurately reproduces the densities and heats of vaporization of a large number of hydrocarbon and ethoxylated hydrocarbon molecules and, therefore, can be assumed to accurately represent the intermolecular interactions between the C10E4 and the C12 E6 molecules.38 The SPC/E water model, where charges are placed on the two hydrogen atoms and the oxygen atom comprising each water molecule, and Lennard-Jones sites are placed on the oxygen atoms,39 was used in this study due to its relatively simple structure and satisfactory reproduction of the aqueous environment around surfactant molecules.40 These simulations were carried out using the GROMACS 4.6 simulation package.41 In each biased-sampling simulation,34 the anchor carbon atoms in the two surfactant molecules were restrained at their initial configuration using a harmonic potential (with a force constant of 1000 kJ mol−1 nm−2). The separation between the two carbon anchor atoms was measured every 0.5 ps over a sampling period of 5 ns. In this manner, each of the 39 simulations generated a histogram for the separation between the anchor atoms, as shown in Figure 3. The weighted-histogram analysis method42−44 (WHAM) was then applied to the time-series data from the 39 simulations to calculate the potential of mean force (PMF) acting between the two C10E4 molecules as a function of the separation between the anchor atoms. This PMF is a direct measure of the intermolecular interaction 8322

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

Figure 5. SDK-CG model of a C10E4 surfactant showing the CG beads (large transparent spheres) overlaying the all-atom structure of the molecule. The SDK-CG bead types are OA (red), EO (green), CM (cyan), and CT (violet). Details of the SDK-CG bead types for PEG surfactants are found in ref 46. The image was generated with VMD using the Tachyon ray tracing library.51,52

Figure 4. Results of the biased-sampling simulations of two C10E4 molecules at the air−water interface. The bottom half of the figure shows the histograms of the separations between the two C10E4 molecules obtained from each of the simulations. The top half of the figure shows the PMF between the two C10E4 molecules reconstructed from the histograms using the WHAM algorithm. The minimum in the PMF curve occurs at an intermolecular separation of 0.48 nm. This suggests that the two C10E4 molecules have a hard-disk radius of 0.24 nm.

⎛ ϕ(x) ⎞⎤ ⎢1 − exp⎜− ⎟⎥x dx ⎢⎣ ⎝ kBT ⎠⎥⎦

Table 1. Ns: Number of Surfactants per Interface; Nw: Number of CG Water Beads; Lx/y = Lx = Ly: Lateral Box Dimensions; Lz: Box Length in Normal Direction of Interface

2.25 ⎡

B=π

∫0.48

(8)

For a pair of C10E4 molecules, the application of eq 8 to the PMF in Figure 4 results in a value of B = −0.386 nm2. The corresponding value for C12E6 is reported in the Results section below. Method 2: r and B Parametrization Based on the 2D Equation of State Fitting. The parameters r and B can also be obtained through MD simulations that are used to explicitly compute the 2D equation of state, i.e., the variation of Π with Γ for both C10E4 and C12E6. The simulations used in Method 2 utilized the coarse-grained (CG) Shinoda−Devane−Klein (SDK) force field,46 which has been specifically parametrized to simulate interfacial behavior and can accurately reproduce the interfacial tensions of both pure water and the surfactant solutions considered.1,47,48 Nonbonded interactions of solute molecules, such as alcohols or ether-containing molecules, have been parametrized by matching partitioning or solvation free energies.46 The force-field parameters for this model of C10E4 and C12E6 have been taken from the literature.46 In this model, typically, three united-atoms are mapped into a single CG bead.46 The model has been recently refined regarding the PEG chains.8 However, for the proof-of-concept demonstration, we will use here the original SDK model. The atomistic structure and the CG mapping of a C10E4 surfactant are shown in Figure 5. All PEG surfactants consist of SDKCG bead types OA, EO, CM, and CT/CT2. The choice of CT vs CT2 depends on the number of methylene groups in the tail; for example, a C10 tail of a PEG surfactant terminates with a CT type, whereas a C12tail terminates with a CT2 type.46 The absence of explicit charges in the CG model of both water and the surfactants allows for very efficient CPU parallelization schemes on computer clusters. All CGMD calculations were performed with the LAMMPS-MD simulation package.49,50 The simulations for Method 2 consisted of 6300 water beads (corresponding to 18,900 water molecules) assembled into an orthogonal slab with vacuum on either side, and up to 220 surfactant molecules per air−water interface, in a simulation box of dimensions 9.4 nm × 9.4 nm × 20 nm. Three different surfactant species: C12E12 (for force-field validation), C12E6, and C10E4 were simulated. For C12E12, a surface concentration of Γ = 1.39 nm−2 was assumed a priori in accordance with neutron reflectivity measurements taken at the CMC.10 For the other two surfactant types, we simulated a range of surface concentrations by varying the number of surfactants at constant surface area; see Table 1. These simulations were carried out at constant volume with a Nose-Hoover thermostat (with a time

Type

Ns

Nw

Lx/y [nm]

Lz [nm]

C12E12 C12E6 C12E6 C10E4

72 50−180 72 4−220

22,249 6,300 22,000 6,300

7.2 9.4 9.4 9.4

80.0 20.0 80.0 20.0

constant of 1 ps),53 with a time step of 10 fs for a total of 120 ns per simulation. We verified that there are no system size effects on the results by increasing the system size for the C12E6 case; see Table 1. The interfacial tension was determined using the principal stress tensor components Pi in the Cartesian directions and Lz, the length of the simulation box in the direction normal to the air−water interface, where γ=

Px + Py Lz Pz − 2 2

(9)

with ⟨...⟩ indicating time averages over the last 110 ns of the simulations. Note that the prefactor accounts for the two interfaces and the box boundary Lz normal to the interface (z). The standard deviation of the computed surface tension values ranges between 0.4 and 0.8 mN m−1. We first tested the use of Method 2 with the SDK force field for C12E12 and obtained γ = (38. ± 0.6) mN m−1 at Γ = 1.39 nm−2, a value reported by neutron reflectivity measurements at the CMC.10 This value agrees very well with the experimentally determined surface tension of 38.5 mN m−1 for C12E12 at the CMC. This agreement was previously reported1 and substantiates the SDK force field as a useful one for verifying our MD/MTT concept as a tool to quantitatively study interfacial properties of surfactant systems. The agreement between the γ values determined by MD simulations and experiment reduces uncertainties stemming from the force field quality with regard to IFT predictions. Determining the Parameter Δμ0. As discussed above, Δμ0 represents the difference between the chemical potential of a surfactant molecule at infinite dilution at the air−water interface and the chemical potential of a surfactant molecule at infinite dilution in the bulk solution phase. This difference in chemical potentials is equal to the difference in the free energy of a surfactant molecule when it is adsorbed at the interface and when it is in the bulk solution phase, and it represents the affinity of the surfactant molecule for the interface. The more negative the value of Δμ0, the greater the affinity of the surfactant molecule for the interface. 8323

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

Figure 6. Anchor carbon atom of C10E4 (shown in green) pulled from a depth of 1.5 nm below the air−water interface to a height of 1 nm above the interface. At 26 equispaced positions between these two end points (for example, 0.2 nm above the air−water interface as illustrated in the green inset on the bottom-left), the anchor carbon atom is restrained using a harmonic potential, and the variation of the position of the anchor carbon atom with respect to the air−water interface is recorded over a 2 ns simulation. These variations from each of the 26 simulations (represented as the series of histograms in the bottom-right) were combined using the WHAM algorithm to determine a PMF with respect to the air−water interface. In the atomistic snapshots, water molecules are represented by red lines, and carbon, oxygen, and hydrogen atoms are shown as black, red, and white spheres, respectively. The black dashed vertical lines in the two graphs represent the air−water interface as determined by the point at which the local density of water is 5% of the bulk density of water. The vertical dotted red lines represent the position of the Gibbs dividing surface.54 Simulation snapshots were generated with VMD using the Tachyon ray tracing library.51,52 To obtain Δμ0, we compute the free-energy change associated with transferring a single surfactant molecule from the bulk solution phase to the interface. To be consistent with the two force fields used in Method 1 and Method 2, the free-energy change was computed using the all-atomistic OPLS-AA force field as well as the coarse-grained SDK force field. The simulations with the OPLS-AA force field were initialized with a simulation box of dimensions 4 nm × 4 nm × 13 nm containing a 4 nm × 4 nm × 7.5 nm slab of 4000 SPC/E water molecules and a single surfactant molecule with its anchor atom placed at a height of 1.5 nm below the surface of the SPC/E slab, in a manner similar to the start of the simulations used for the calculation of r and B. The surface of the slab was taken to be the locus of points where the local density of water molecules is 5% of the bulk density of water molecules, to be consistent with the definition commonly utilized in MD simulations of the air−water interface.54 Using an externally applied harmonic potential (with a force constant of 1000 kJ mol−1 nm−2), the surfactant molecule was pulled toward the empty space above the slab of SPC/E water molecules. No additional restraints were applied on the surfactant molecule in any other direction. Simulation snapshots were extracted at 26 equispaced positions between a depth of 1.5 nm below the surface of the SPC/E water slab and a height of 1.0 nm above the surface. These simulation snapshots were subsequently utilized as the initial configurations for a series of biased-sampling simulations where the anchor atom was restrained (with a harmonic potential with a force constant of 2000 kJ mol−1 nm−2) at a fixed height above, or depth below, the air−water interface. Each of these 26 simulations was carried out in the NVT ensemble using a NoseHoover barostat (with a time constant of 2 ps) at a temperature of 300 K, for 2 ns with a time step of 2 fs. The position of the anchor carbon atom relative to the air−water interface was recorded every 500 fs. The time-series data from all the 26 simulations were then processed using the WHAM algorithm to generate the change in the free energy of the

surfactant molecule as a function of its depth below the air−water interface. Convergence of the free energy profile was verified by comparing the profile obtained from the 1.5 ns trajectories with that obtained from the 2.0 ns trajectories (see Supporting Information for details). Figure 6 illustrates this simulation workflow for a C10E4 molecule modeled using the OPLS-AA framework. The free-energy change calculated using this method for C10E4 and C12E6 as a function of the height of the surfactant molecule above the interface is shown in Figure 7. Both surfactants have a significantly lower free energy at the interface than in the bulk solution phase (at a height of ≈1.5 nm above the interface) or in the ”air”’ above the water

Figure 7. Free-energy change obtained using the OPLS-AA force field for C10E4 and C12E6 as a function of the height of their anchor atoms above the air−water interface. The reference-state chemical potential difference is estimated based on the depth of the minima in the PMF. The black dashed vertical line represents the air−water interface as determined by the point at which the local density of water is 5% of the bulk density of water. The red dashed vertical line represents the position of the Gibbs dividing surface. 8324

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir slab (at a height of ≈1.0 nm above the interface). This result is consistent with the intuitive picture that surfactants prefer to adsorb at interfaces. The ideal equilibrium position for C10E4 and C12E6 (i.e., the position of the minima in the free-energy profiles) coincides with the position of the Gibbs dividing surface (indicated by the red dashed vertical line in Figure 7). Accordingly, the desired Δμ0 values for C10E4 and C12E6 can be conveniently estimated by the depths of the corresponding minima, or equivalently, by the free-energy changes when C10E4 and C12E6 are located at the corresponding Gibbs dividing surfaces. Figure 7 also suggests that C12E6 has a higher affinity for the interface and, therefore, a more negative Δμ0 value than C10E4.23 In the umbrella sampling simulations carried out using the SDK force field, the stiffness of the umbrella potential was kept at 1255 kJ mol−1 nm−2, and neighboring distance windows were separated by 0.05 nm, corresponding to about 50 bins for each free energy profile calculation. Convergence of the free energy profile was verified by extension of the simulation time from 10 to 100 ns per bin (see Supporting Information for details).

standard deviation of the set of r and B values obtained from these 10 different PMF curves is reported in Table 2. The results of the simulations carried out in Method 2 to estimate γ(Γ) for C10E4 and C12E6 are shown in Figure 9.

■

RESULTS Estimates for the Parameters r and B. The PMF between a pair of C12E6 molecules adsorbed at the air−water interface is compared against the PMF between a pair of C10E4 molecules in Figure 8. The increasing size of C12E6 relative to

Figure 9. Interfacial tension at the air−water interface γ = γ0 − Π versus the surfactant surface coverage Γ for C10E4 and C12E6. The SDK force field used for these simulations has been parametrized to yield a value of γ0 = 72 mN m−1 for water at 298 K. The circles and squares denote the values of γ computed using MD simulations. The solid lines represent the fits of the MD simulation data to eq 1. For the fit procedure, boundaries were imposed for a > 0 and B < 0 (attraction). The dashed lines represent the fit to eq 1 with B set to zero. The error in the MD simulation estimates of the interfacial tension is typically 0.2 mN m−1 (i.e., smaller than the size of the symbols).

Figure 9 indicates that we obtain reasonable fits of the simulation data for both surfactants over a wide range of surface concentrations, Γ. Note that the surface tension in the absence of surfactants is γ0 = 72 mN m−1, a value that the SDK CG water model yields by construction.48 By fitting eq 1 to the IFT data obtained from the MD simulations (see Figure 9), one can obtain estimates of the hard-disk radius r and the pairwise attractive intermolecular interaction term B. The fit values for r and B for both surfactants are given in Table 2 under Method 2. The error bars in r are quite small, whereas the standard deviations in B can be quite large if one does not include sufficient data points in the fitting procedure. In order to extract precise values of B, one has therefore to ensure small error bars in the calculated surface tensions (≲0.2 mN m−1) and/or obtain many data points for γ(Γ). We note, however, that the tail attraction contribution to the surface pressure is quite weak, as can be seen from the single parameter fit in Figure 9 (dashed line), where B is set to zero. In this case, the fit parameters for r are 0.253 and 0.283 nm for C10E4 and C12E6, respectively, only slightly different from the values shown in Table 2 under Method 2. The values of r and B obtained using Method 1 are consistent with those obtained using Method 2, suggesting that

Figure 8. PMF between pairs of surfactant molecules computed at the air−water interface using MD simulations for the surfactants C10E4 and C12E6. The minima in the PMF curves of C10E4 and C12E6 occur at intermolecular separations of 0.48 and 0.56 nm, respectively. This suggests that the C10E4 and the C12E6 molecules have hard-disk radii of 0.24 and 0.28 nm, respectively.

C10E4 is reflected by the increase in the position of the minima in the PMF curves: C10E4 and C12E6 have hard-disk radii of 0.24 and 0.28 nm, respectively. Further, the application of eq 8 to these PMF curves yields values of −0.386 nm2 and −0.130 nm2 for the B parameters for C10E4 and C12E6, respectively. These values are reported in Table 2 under Method 1. Error estimates for r and B were obtained by extracting (with replacement) 10 subsets of 1000 values of the intermolecular separations from each of the 39 biased-sampling windows and then applying the WHAM algorithm to generate 10 different PMF curves. The

Table 2. Parameters for the 2D Equation of State for C10E4 and C12E6 Obtained Using (i) Method 1, by Calculating r and B Explicitly Using Eq 8 from the Simulation Data in Figure 8, and (ii) Method 2, by Fitting Eq 1 to the Simulation Data in Figure 9 Method 1

Method 2 2

Surfactant

r (nm)

B (nm )

r (nm)

B (nm2)

C10E4 C12E6

0.240 ± 0.015 0.280 ± 0.008

−0.386 ± 0.030 −0.130 ± 0.021

0.265 ± 0.001 0.287 ± 0.001

−0.332 ± 0.020 −0.141 ± 0.041

8325

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir both the all-atom force field used in Method 1 and the coarsegrained force field used in Method 2 are appropriate in their representation of the intermolecular interactions. This is particularly significant given that the force fields used for water and the two surfactants C10E4 and C12E6 in Method 1 were parametrized to reproduce bulk properties such as the densities and heats of vaporization of water and a variety of organic compounds, respectively, while the force field used in Method 2 was designed to reproduce the interfacial properties of water and surfactants.55,56 Comparing the MTT parameters listed in Table 2, we notice that the effective headgroup area of C12E6 is larger than that of C10E4, as can be expected due to the larger hydrophilic group in the former surfactant. Less obvious is the stronger attraction between the tails for C10E4 despite the lower molecular weight of the alkyl chain compared to that of C12E6. The main reason is that the average distance between the tails is restricted by the headgroup size, which in the C10E4 case is smaller than the size of a C12E6 surfactant. This can be readily seen from eq 8, where the lower integration limit is set by the headgroup diameter 2r. Evidently, large contributions to B stem from interaction potentials at small separations. Since the headgroup area of C10E4 is smaller than that of C12E6, |B| becomes larger (more attractive) despite the lower molecular weight of the tail. Therefore, one has to be cautious in ignoring tail attraction contributions to the surface pressure, i.e., setting B = 0. Although such an approximation seems to work for C12E6 and C10E4 (compare the dashed and solid lines in Figure 9), this may not always be the case, particularly for surfactants with small headgroup sizes and long alkyl tails. The values of r reported in Table 2 may be compared to estimates of r obtained from the molecular structures of the C10E4 and the C12E6 molecules. Nikas et al. estimated the harddisk radii of these molecules to be 0.307 and 0.340 nm, respectively,23 which are significantly larger than the values in Table 2. However, as noted by Nikas et al., these differences could be rationalized by considering the difficulties in (i) assigning a hard-disk area to a flexible chainlike macromolecule, (ii) accounting for changes in chain configurations when surfactant molecules are adsorbed at the air−water interface, and (iii) determining the dependence of the conformational entropy of chainlike surfactant molecules on surface coverage. Finally, we comment on the potential role of many-body effects that could influence the magnitude of B. Such effects are not captured in the surface pressure, as it is formulated in eq 1. The good quality of the fits shown in Figure 9 suggests that many-body effects are not significant in the surface concentration range that is considered here but may become important at higher surface concentrations. However, in practice, consideration of a high surface concentration range is not relevant because the onset of micellization above the CMC limits such high packing at the surface. Those limits are denoted in adsorption isotherms; see the next section. Estimates for the Parameter Δμ0. The Δμ0 values for C10E4 and C12E6, obtained using both the OPLS-AA and the SDK force fields, are shown in Table 3. Once again, it is noteworthy that, in spite of the difference in the physical properties that were used to parametrize each of the two force fields (bulk properties and interfacial properties, respectively), the two force fields produce estimates for Δμ0 that are within 1 kBT at room temperature. This difference corresponds to a factor of e (≈ 2.7) in the value of Γ for a given bulk concentration c and is within the factor of 5 difference between

Table 3. Values of the Reference-State Chemical Potential Differences, Δμ0, for C10E4 and C12E6 Obtained Using MD Simulations Based on the OPLS-AA Force Field and the Coarse-Grained Sdk Force Fielda

a

Surfactant

OPLS-AA Force Field

SDK Force Field

C10E4 C12E6

−34.2 kJ mol−1 −42.8 kJ mol−1

−36.9 kJ mol−1 −44.6 kJ mol−1

As can be seen, good agreement is obtained.

experimental values and theoretical predictions that previous studies have used as a benchmark for accuracy.24 The obtained parameter values for C10E4 and C12E6 are qualitatively in agreement with those reported in the original papers.23,24 There are, however, quantitative differences which, apart from influence of the force field on the MD simulations, stem from fundamental differences in the parametrization methods. For example, in the original paper by Nikas et al.,23 the headgroup area and the second virial coefficient were estimated on the basis of polymer theory or Monte Carlo simulation of two hydrocarbon chains end-grafted at a wall. Further, the free energy of adsorption was calculated from experimental measurements of the surface pressure at one state point. Our approach relies solely on a single MD model and, hence, exhibits a consistent parametrization method of the MTT model. In general, the downside of an MD/MTT approach is that potential flaws in the force field quality would translate into poor predictions of the surface tension. Here, the calculated transfer free energy seems to be accurate, because the predicted surface tension agrees well with the experimental data, as will be demonstrated in the next section. Predicting the Surface Tension Isotherms. Having estimated the parameters required for the implementation of the MD/MTT framework, one can readily predict adsorption isotherms, and air−water interfacial tensions as a function of bulk solution surfactant concentration. To predict the adsorption isotherm, eq 7 was solved for the adsorption at each given bulk solution mole fraction x using the MD/MTT parameters shown in Table 2 and Table 3. Figure 10a shows the adsorption Γ as a function of the bulk solution mole fraction x for C10E4 and C12E6 using Methods 1 and 2, and the OPLS-AA and the SDK force fields, respectively, and Figure 10b juxtaposes the predicted and the experimental IFTs as a function of x. Figure 10a indicates that, due to its higher surface affinity, as reflected in a more negative value of Δμ0 (see Table 3), C12E6 shows higher surface adsorption than C10E4 at the same x value, for values of x below the CMC of C12E6 (the region of applicability of the MD/MTT framework). Furthermore, Figure 10a suggests that above x ≈ 1 × 10−6 (where the black and blue Γ(x) curves intersect), it becomes possible for C10E4 to pack better than C12E6 due to the former’s smaller headgroup area. As seen in Figure 10b, the agreement between the predictions of the MD/MTT framework and Method 1 is very good over the region of applicability of the MD/MTT framework. Recall that the OPLS-AA force field was parametrized to match several bulk properties of organic moieties and, as such, was not designed to accurately reproduce interfacial properties. Nevertheless, using the MD/MTT framework, along with this force field and Method 1, yields a reasonable prediction of the adsorption isotherms for C10E4 and C12E6. 8326

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

C10E4 is larger than that of C12E6, the IFT of the former surfactant levels off at a lower γ value.

■

CONCLUSIONS AND FUTURE DIRECTIONS In this paper, we presented a novel theoretical framework to predict interfacial properties, such as the adsorption isotherms, of nonionic surfactant solutions based solely on the chemical structures of the surfactant molecules and in the absence of any empirical fitted parameters. To summarize this framework, in order to predict the interfacial tension isotherm for a given nonionic surfactant molecule, the following steps are involved: 1. The surfactant molecule is represented using either an all-atom or a coarse-grained MD force field (such as the OPLS/AA or the SDK force fields). 2. The molecular-thermodynamic parameters r and B are computed using MD simulations via either Method 1 (where pairwise interactions are calculated from the PMF between a pair of surfactant molecules adsorbed at the air−water interface) or Method 2 (where pairwise interactions are calculated from the relationship between the surface pressure Π and the adsorption Γ of the surfactant molecules at the air−water interface, as obtained from a series of MD simulations with different values of Γ). 3. The chemical potential difference Δμ0 is determined by using MD simulations to calculate the free-energy change associated with transferring a surfactant molecule from the bulk solution phase to the air−water interface. 4. Equation 7 is then solved for the adsorption Γ at each given bulk solution mole fraction x using the MD/MTT parameters obtained in the previous two steps. Equation 1 may then be used to obtain the surface tension isotherm γ(x). The utilization of conventional MD simulations together with MTT calculations provides a powerful modeling tool for quantitative predictions of adsorption isotherms and IFT predictions as a function of bulk surfactant concentration. In this hybrid MD/MTT approach, the molecular-thermodynamic theory acts as a bridging tool to relate the interfacial tension to the bulk solution surfactant concentration, while the MD simulations determine the molecular and thermodynamic (r, B, and Δμ0) parameters that appear in the molecular-thermodynamic theory. This enables the MD/MTT framework to accurately capture the partitioning of the surfactant molecules between the bulk solution and the interface. Therefore, the MD/MTT approach overcomes a major limitation of atomistic MD simulationsbecause they only cover limited length- and time-scales, atomistic simulations cannot accurately capture the equilibrium between the bulk and the interface. We also demonstrated the predictive power of the MD/ MTT framework for nonionic alkyl ethoxylated surfactants adsorbed at the air−water interface, and we showed that one can obtain very good agreement between the MD/MTT predictions of surfactant-induced interfacial tension lowering for C10E4 and C12E6, using both all-atomistic and coarse-grained MD force fields, and the experimental results (Figure 10a and Figure 10b). The MD/MTT framework does have some limitations. Although we have demonstrated that the framework yields reasonable predictions for both the OPLS-AA and the SDK force fields for the C10E4 and C12E6 alkyl ethoxylated surfactants, it is yet to be determined whether a similarly

Figure 10. (a) Surfactant adsorption Γ at the air−water interface as a function of the bulk solution surfactant mole fraction x predicted by the MD/MTT framework. (b) The surface tension isotherm γ(x) relating the interfacial tension at the air−water interface as a function of x for C10E4 and C12E6. The circles and squares denote experimental measurements23 for C10E4 and C12E6, respectively. The dashed black and blue lines represent predictions for C10E4 and C12E6, respectively, using Method 1 and the OPLS-AA force field. The solid black and blue lines represent predictions for C10E4 and C12E6, respectively, using Method 2 and the SDK force field. The green and red vertical lines represent experimental critical micelle concentrations and provide an upper limit for the range of applicability of the MTT framework, for C10E4 and C12E6, respectively.

Figure 10b also shows that there is excellent agreement between the surface tension isotherms predicted by the MD/ MTT framework with the SDK force field and Method 2 and the surface tension isotherms determined experimentally for surfactant concentrations below the CMC, which is the surfactant concentration range for which the MD/MTT theory was developed. This is quite exciting, considering that the MD/ MTT framework relies solely on the accurate parametrization of the SDK force fields using the interfacial properties of a family of organic moieties, and does not contain any other empirical parameters. We note that, regardless of the force field used in the MD/ MTT framework, while the experimental data level off above the CMC, the predicted MD/MTT curve continues to decline since we did not include micellization into the calculations. In reality, the CMC limits the packing of surfactants at the surface and, therefore, also limits the monotonic surface tension decline shown in Figure 10(b). Note that since the CMC of 8327

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir

(5) Tang, X.; Huston, K. J.; Larson, R. G. Molecular Dynamics Simulations of Structure-Property Relationships of Tween 80 Surfactants in Water and at Interfaces. J. Phys. Chem. B 2014, 118, 12907−12918. (6) Isele-Holder, R. E.; Ismail, A. E. Atomistic Potentials for Trisiloxane, Alkyl Ethoxylate, and Perfluoroalkane-based Surfactants with TIP4P/2005 and Appliction to Simulations at the Air-Water Interface. J. Phys. Chem. B 2014, 118, 9284−9297. (7) Huston, K. J.; Larson, R. G. Reversible and Irreversible Adsorption Energetics of Poly(ethylene glycol) and Sorbitan Poly(ethoxylate) at a Water/Alkane Interface. Langmuir 2015, 31, 7503− 7511. (8) Shen, Z.; Sun, H. Prediction of Surface and Bulk Partition of Nonionic Surfactants Using Free Energy Calculations. J. Phys. Chem. B 2015, 119, 15623−15630. (9) Lu, J. R.; Marrocco, A.; Su, T. J.; Li, Z. X.; Thomas, R. K.; Penfold, J. Adsorption of Dodecyl Sulfate Surfactants with Monovalent Metal Counterions at the Air-Water Interface Studied by Neutron Reflection and Surface Tension. J. Colloid Interface Sci. 1993, 158, 303−316. (10) 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, 10332−10339. (11) Tamam, L.; Pontoni, D.; Sapir, Z.; Yefet, S.; Sloutskin, E.; Ocko, B. M.; Reichert, H.; Deutsch, M. Modification of Deeply Buried Hydrophobic Interfaces by Ionic Surfactants. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 5522−5525. (12) Dzikovski, B. G.; Livshits, V. A. EPR Spin Probe Study of Molecular Ordering and Dynamics in Monolayers at Oil/Water Interfaces. Phys. Chem. Chem. Phys. 2003, 5, 5271−5278. (13) Israelachvili, J. N. Intermolecular and Surface Forces, 3rd ed.; Academic Press, 2011. (14) Frumkin, A. N. Electrocapillary Curve of Higher Aliphatic Acids and the State Equation of the Surface Layer. Int. J. Res. Phys. Chem. Chem. Phys. 1925, 116, 466−488. (15) Fainerman, V. B.; Miller, R.; Aksenenko, E. V. Simple Model for Prediction of Surface Tension of Mixed Surfactant Solutions. Adv. Colloid Interface Sci. 2002, 96, 339−59. (16) 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. (17) Kolev, V. L.; Danov, K. D.; Kralchevsky, P. A.; Broze, G.; Mehreteab, A. Comparison of the van der Waals and Frumkin Adsorption Isotherms for Sodium Dodecyl Sulfate at Various Salt Concentrations. Langmuir 2002, 18, 9106−9109. (18) Rusanov, A. I. Novel Approach to the Equation of State for Fluid Systems, Based on the Concept of the Exclusion Factor. Russ. Chem. Rev. 2005, 74, 117−127. (19) Franses, E. I.; Siddiqui, F. A.; Ahn, D. J.; Chang, C.-H.; Wang, N.-H. L. Thermodynamically Consistent Equilibrium Adsorption Isotherms for Mixtures of Different-Sized Molecules. Langmuir 1995, 11, 3177−3183. (20) Puvvada, S.; Blankschtein, D. Molecular-Thermodynamic Approach to Predict Micellization, Phase Behavior, and Phase Separation of Micellar Solutions. I. Application to Nonionic Surfactants. J. Chem. Phys. 1990, 92, 3710−3724. (21) Nagarajan, R.; Ruckenstein, E. Theory of Surfactant SelfAssembly: A Predictive Molecular Thermodynamic Approach. Langmuir 1991, 7, 2934−2969. (22) Peck, D. G.; Schechter, R. S.; Johnston, K. P. Unified Classical and Molecular Thermodynamic Theory of Spherical Water-in-Oil Microemulsions. J. Phys. Chem. 1991, 95, 9541−9549. (23) Nikas, Y. J.; Puvvada, S.; Blankschtein, D. Surface Tensions of Aqueous Nonionic Surfactant Mixtures. Langmuir 1992, 8, 2680− 2689. (24) Mulqueen, M.; Blankschtein, D. Prediction of Equilibrium Surface Tension and Surface Adsorption of Aqueous Surfactant

satisfactory agreement may be obtained for more complex nonionic surfactants or using other force fields. Furthermore, this molecular-thermodynamic model, where the adsorbed surfactant molecules are treated as hard disks on a planar surface, with pairwise tail contributions added separately, may also break down for more complex surfactant architectures, or at higher surfactant surface coverages, where third-order and higher-order molecular interactions become important. Finally, it should be noted that although the applicability of the MD/MTT framework was demonstrated for single surfactant systems of C10E4 and C12E6 nonionic surfactants, a molecular-thermodynamic formulation is already available for surfactant mixtures. 30 This implies that the molecular parameters evaluated here (see Table 2 and Table 3) could be used to determine the surface tension of mixtures of these nonionic ethoxylated surfactants without carrying out any additional MD simulations. Additionally, the effect of temperature on the robustness of the predictive capabilities of this MD/MTT framework needs to be evaluated. This aspect, along with the development of a new MD/MTT framework for ionic and zwitterionic surfactants at air−water and oil−water interfaces, is the subject of ongoing investigations.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.7b01073. Effect of sampling time on the convergence of freeenergy calculations for Δμ0, procedure for the identification of the anchor atom, and results obtained using the Flory−Huggins model for the chemical potential of surfactants in the bulk (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: arben.jusufi@exxonmobil.com. ORCID

Daniel Blankschtein: 0000-0002-7836-415X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Daniel P. Cherney, Aruna Mohan, Aditya Jaishankar, and Mohsen Yeganeh, all at ExxonMobil Research & Engineering Co., for helpful discussions.

■

REFERENCES

(1) Shinoda, W.; DeVane, R.; Klein, M. L. Coarse-grained Molecular Modeling of Non-ionic Surfactant Self-assembly. Soft Matter 2008, 4, 2454−2462. (2) Yang, W.; Wu, R.; Kong, B.; Zhang, X.; Yang, X. Molecular Dynamics Simulations of Film Rupture in Water/Surfactant Systems. J. Phys. Chem. B 2009, 113, 8332−8338. (3) Cheng, T.; Chen, Q.; Li, F.; Sun, H. Classical Force Field for Predicting Surface Tension and Interfacial Properties of Sodium Dodecyl Sulfate. J. Phys. Chem. B 2010, 114, 13736−13744. (4) Martínez, H.; Chacón, E.; Tarazona, P.; Bresme, F. The Intrinsic Interfacial Structure of Ionic Surfactant Monolayers at Water-Oil and Water-Vapour Interfaces. Proc. R. Soc. London, Ser. A 2011, 467, 1939− 1958. 8328

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329

Article

Langmuir Mixtures Containing Ionic Surfactants. Langmuir 1999, 15, 8832− 8848. (25) Nagarajan, R.; Ruckenstein, E. Molecular Theory of Microemulsions. Langmuir 2000, 16, 6400−6415. (26) Mulqueen, M.; Blankschtein, D. Theoretical and Experimental Investigation of the Equilibrium Oil-Water Interfacial Tensions of Solutions Containing Surfactant Mixtures. Langmuir 2002, 18, 365− 376. (27) Moreira, L. A.; Firoozabadi, A. Molecular Thermodynamic Modeling of Droplet-Type Microemulsions. Langmuir 2012, 28, 1738−1752. (28) Nagarajan, R.; Ruckenstein, E. Aggregation of Amphiphiles as Micelles or Vesicles in Aqueous Media. J. Colloid Interface Sci. 1979, 71, 580−604. (29) Maibaum, L.; Dinner, A. R.; Chandler, D. Micelle Formation and the Hydrophobic Effect. J. Phys. Chem. B 2004, 108, 6778−6781. (30) Mulqueen, M.; Blankschtein, D. Prediction of Equilibrium Surface Tension and Surface Adsorption of Aqueous Surfactant Mixtures Containing Zwitterionic Surfactants. Langmuir 2000, 16, 7640−7654. (31) Stephenson, B. C.; Goldsipe, A.; Beers, K. J.; Blankschtein, D. Quantifying the Hydrophobic Effect. 1. A Computer SimulationMolecular-Thermodynamic Model for the Self-Assembly of Hydrophobic and Amphiphilic solutes in Aqueous Solution. J. Phys. Chem. B 2007, 111, 1025−1044. (32) Holmberg, K.; Bo, J.; Kronberg, B. Surfactants and Polymers in Aqueous Solutions; John Wiley & Sons Ltd., 2002. (33) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (34) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations using Jarzynski’s Equality. J. Chem. Phys. 2003, 119, 3559−3566. (35) Stephenson, B. C.; Beers, K.; Blankschtein, D. Complementary Use of Simulations and Molecular-Thermodynamic Theory to Model Micellization. Langmuir 2006, 22, 1500−1513. (36) Stephenson, B. C.; Stafford, K. A.; Beers, K. J.; Blankschtein, D. Application of Computer Simulation Free-Energy Methods to Compute the Free Energy of Micellization as a Function of Micelle Composition. 1. Theory. J. Phys. Chem. B 2008, 112, 1634−40. (37) Iyer, J.; Mendenhall, J. D.; Blankschtein, D. Computer Simulation-Molecular-Thermodynamic Framework to Predict the Micellization Behavior of Mixtures of Surfactants: Application to Binary Surfactant Mixtures. J. Phys. Chem. B 2013, 117, 6430−6442. (38) Siu, S. W. I.; Pluhackova, K.; Böckmann, R. A. Optimization of the OPLS-AA Force Field for Long Hydrocarbons. J. Chem. Theory Comput. 2012, 8, 1459−1470. (39) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (40) Rideg, N. A.; Darvas, M.; Varga, I.; Jedlovszky, P. Lateral Dynamics of Surfactants at the Free Water Surface: A Computer Simulation Study. Langmuir 2012, 28, 14944−14953. (41) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GRGMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (42) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (43) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollmann, P. A. Multidimensional Free-Energy Calculations using the Weighted Histogram Analysis Method. J. Comput. Chem. 1995, 16, 1339−1350. (44) Grossfield, A. WHAM: The Weighted Histogram Analysis Method, v. 2.0.2; http://membrane.urmc.rochester.edu/content/ wham. 2008. (45) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313.

(46) Shinoda, W.; DeVane, R.; Klein, M. L. Multi-Property Fitting and Parameterization of a Coarse Grained Model for Aqueous Surfactants. Mol. Simul. 2007, 33, 27−36. (47) Klein, M. L.; Shinoda, W. Large-scale Molecular Dynamics Simulations of Self-Assembling Systems. Science 2008, 321, 798−800. (48) He, X.; Shinoda, W.; DeVane, R.; Klein, M. L. Exploring the Utility of Coarse-Grained Water Models for Computational Studies of Interfacial Systems. Mol. Phys. 2010, 108, 2007−2020. (49) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (50) LAMMPS; http://lammps.sandia.gov, 2015. (51) Humphrey, W.; Dalke, A.; Schulten, K. VMD − Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (52) Stone, J. An Efficient Library for Parallel Ray Tracing and Animation. M.Sc. thesis, Computer Science Department, University of Missouri-Rolla, 1998. (53) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (54) Horinek, D.; Herz, A.; Vrbka, L.; Sedlmeier, F.; Mamatkulov, S. I.; Netz, R. R. Specific Ion Adsorption at the Air/Water Interface: The Role of Hydrophobic Solvation. Chem. Phys. Lett. 2009, 479, 173−183. (55) Zielkiewicz, J. Structural Properties of Water: Comparison of the SPC, SPCE, TIP4P, and TIP5P Models of Water. J. Chem. Phys. 2005, 123, 104501. (56) Shinoda, W.; DeVane, R.; Klein, M. L. Coarse-grained Force Field for Ionic Surfactants. Soft Matter 2011, 7, 6178.

8329

DOI: 10.1021/acs.langmuir.7b01073 Langmuir 2017, 33, 8319−8329