Article pubs.acs.org/Langmuir
Discrete Fractional Component Monte Carlo Simulation Study of Dilute Nonionic Surfactants at the Air−Water Interface Brian Yoo,†,§ Eliseo Marin-Rimoldi,† Ryan Gotchy Mullen,† Arben Jusufi,‡ and Edward J. Maginn*,† †
Department of Chemical and Biomolecular Engineering, University of Notre Dame, 182 Fitzpatrick Hall, Notre Dame, Indiana 46556-5637, United States ‡ Corporate Strategic Research, ExxonMobil Research and Engineering Company, 1545 U.S. 22, Annandale, New Jersey 08801, United States S Supporting Information *
ABSTRACT: We present a newly developed Monte Carlo scheme to predict bulk surfactant concentrations and surface tensions at the air− water interface for various surfactant interfacial coverages. Since the concentration regimes of these systems of interest are typically very dilute (≪10−5 mol. frac.), Monte Carlo simulations with the use of insertion/deletion moves can provide the ability to overcome finite system size limitations that often prohibit the use of modern molecular simulation techniques. In performing these simulations, we use the discrete fractional component Monte Carlo (DFCMC) method in the Gibbs ensemble framework, which allows us to separate the bulk and air−water interface into two separate boxes and efficiently swap tetraethylene glycol surfactants C10E4 between boxes. Combining this move with preferential translations, volume biased insertions, and Wang−Landau biasing vastly enhances sampling and helps overcome the classical “insertion problem”, often encountered in non-lattice Monte Carlo simulations. We demonstrate that this methodology is both consistent with the original molecular thermodynamic theory (MTT) of Blankschtein and co-workers, as well as their recently modified theory (MD/MTT), which incorporates the results of surfactant infinite dilution transfer free energies and surface tension calculations obtained from molecular dynamics simulations. simulations is typically fixed. Alternatively, Monte Carlo (MC) simulations in open ensembles could be performed to allow the system to reach lower concentrations than that predefined by the system size. In these techniques, molecules are inserted and deleted such that the equilibrium chemical potential is attained for each species in a given phase. In practice, however, efficient sampling for the insertion and deletion of large molecules in condensed phases is challenging and requires state-of-the-art biasing moves that help overcome the high energy penalties associated with deleting a large molecule and also with inserting the molecule into a crowded and strongly associated condensed phase. This is a classical problem often encountered in the Monte Carlo simulation community and is also commonly referred to as the “insertion problem”.17 Howes and Radke18,19 used MC simulations to estimate the bulk and interfacial concentrations for surfactants at the air− water interface. Their simulations were performed in the canonical ensemble (NVT) using a soft Lennard−Jones model for both pure component and mixed surfactant systems. To overcome the insertion problem, the authors implemented an
1. INTRODUCTION The study of surfactants at the air−water interface has received considerable attention over the past several decades. Also commonly referred to as the Langmuir monolayer, such systems have been widely used as a model to provide a fundamental understanding of surfactant behavior at fluid− vapor and fluid−fluid interfaces. This has been crucial to characterizing the phase behavior of lipid bilayers,1−3 a vital component in cell membranes, guiding the design of thin films,4−6 understanding the mechanisms of foam formation,7 and predicting the interfacial properties of water and oil systems8−11 which are ubiquitous in many industrially relevant applications. A major topic of recent studies has been on the development and improvement of existing theoretical models that can accurately predict the adsorption behavior and interfacial properties of these systems.12−16 In the past few years, there have been some attempts at using classical molecular simulations to predict the adsorption isotherms of surfactants at the air−water and water−oil interface. A major limitation in using molecular simulations has been the need to simulate the dilute concentration ranges of surfactants in the bulk phase, which would translate to very large system sizes using simulation methods such as molecular dynamics (MD), where the number of molecules in the © 2017 American Chemical Society
Received: June 16, 2017 Revised: August 25, 2017 Published: August 28, 2017 9793
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir identity swap move,20 in which a surfactant molecule was exchanged or swapped with a number of neighboring water beads. They mentioned that this move was instrumental in the efficient sampling of configurational space. Likewise, Tomassone et al.21,22 studied a similar soft Lennard−Jones based system of surfactants at the air−water interface using MD simulations in the canonical ensemble. Similarly, the authors were also able to obtain the equilibrium surface tensions along with the corresponding bulk concentrations and surfactant coverages. In both the works of Howes and Radke and Tomassone et al., bulk concentrations and coverages were estimated based on an accurate location of the Gibbs dividing surface (which defines the partitioning of the interface and the bulk). This surface was estimated by determining the position along the axis normal to the interface such that the integration of the solvent number density profiles resulted in zero excess. Subsequently, the bulk and interfacial concentrations were determined by integrating the number density profiles, starting from their respective sides, up to the position of the Gibbs dividing surface. For these concentrations to be accurate, the solvent slab size had to be large enough such that a constant equilibrium density profile was obtained in both the bulk and the interface. In another approach, Rekvig et al.23,24 used dissipative particle dynamics (DPD) combined with a MC method to predict the adsorption isotherm of surfactants at the water−oil interface. In this method, DPD was independently performed for a bulk aqueous box and a water−oil interface box. Note that by doing so, a Gibbs dividing surface did not need to be defined to estimate a bulk concentration or coverage, provided that surfactants in the interface box exclusively resided at the interface. The surfactants were present in both boxes and represented by soft DPD beads.25,26 MC moves were also performed in conjunction, in which surfactants were inserted or deleted into either box from a reservoir (similar to the grand canonical ensemble) such that the chemical potential for the surfactants in the two phases were equal. For this particular model, the configurational bias Monte Carlo (CBMC)27 move was sufficient for efficiently inserting and deleting the surfactant molecules. The authors were able to apply this method to estimate surface tensions and coverages as a function of bulk concentrations in the dilute limit. In the aforementioned studies, much of the success in obtaining reasonable sampling was attributed to the softness of the models. More accurate and realistic models have been developed since, such as those that have been parametrized on the basis of quantum mechanical calculations of the intermolecular and intramolecular interactions (e.g., potential of mean force, bond vibrations, angle bends, dihedral torsions, etc.). For these “hard” models, such as atomistic models or coarse-grained models that have been mapped from the results obtained from atomistic models, poor sampling becomes a major bottleneck with MC methodologies. Even special biasing moves such as the identity swap and CBMC moves are no longer as easily implementable or are insufficient. Recently, there have been increasing efforts to use such models that can accurately predict interfacial properties that are in good agreement with experimental data. However, since the concentration regimes of interest for these systems are often too low to simulate using unbiased MD simulations, several groups have resorted to combining the results from the atomistic or coarse-grained MD simulations with traditional equations of state, including one based on molecular
thermodynamic theory (MTT) developed by Blankschtein and co-workers.14−16,28−30 Shen and Sun28 used thermodynamic integration methods with MD to obtain the adsorption isotherms for surfactants at the air−water interface by simulating two uncoupled simulation boxes (with one box corresponding to the bulk phase and the other to the interface phase) using the coarse-grained SDK model.31 The differences in excess chemical potentials were then used to determine the bulk concentration of surfactants at a particular coverage. In a similar approach, Huston and Larson used a potential of mean force (PMF) method in conjunction with MD simulations to obtain the same adsorption free energy required in the MTT model of Blankschtein and co-workers for polyethoxylates at the water−oil interface.29 They further extended their work to investigate the reversible or irreversible adsorption and desorption of these oligomers at various coverages. Likewise in a very recent study, Sresht et al. demonstrated another method which combines PMF calculations from MD simulations with MTT (also referred to as MD/MTT) by rederiving the original theoretical equations to account for the difference in reference states between molecular simulations and experiments.30 In this work, we present a general MC scheme in the Gibbs ensemble32 to directly predict the adsorption isotherm of a dilute nonionic surfactant, tetraethylene glycol monodecyl ether (C10E4), at the air−water interface using the SDK model. The advantage of using a Gibbs ensemble MC framework is twofold. First, we can separate the interface from the bulk, removing any ambiguity on the location of the Gibbs dividing surface. Second, this method provides a framework which can allow the system to equilibrate to both the dilute and high (above the CMC) concentration regime without the need to simulate a large system size, provided that the MC moves are properly proposed. As previously discussed, developing moves that can overcome the insertion problem for hard models such as the SDK model is nontrivial. Shi and Maginn33,34 have developed the continuous fractional component (CFC) Monte Carlo method that can be applied to the Gibbs Ensemble. This method allows a pair of molecules to be swapped between two boxes through gradual progression by incorporating a coupling parameter (λ) which scales the interactions of a molecule in each of the boxes. One of these molecules is assigned a λ value, while the other is assigned 1 − λ, where a value of λ = 1.0 indicates a molecule is fully interacting and a value of λ = 0.0 indicates the molecule is noninteracting. All intermediate states or subensembles (0.0 < λ < 1.0) are considered “fractional”. When λ goes beyond a value of 1.0 or 0.0, the molecules are swapped and a new pair of molecules is selected and coupled with λ. This method also utilizes an additional biasing potential obtained iteratively using a Wang−Landau approach to help overcome the insertion problem.35,36 Recently, Poursaeidesfahani and co-workers37,38 modified this method so that only a single molecule is coupled with λ rather than a pair of molecules, allowing for the direct calculation of chemical potentials in each of the boxes or phases. For our system, we have implemented a modified version of the CFC method known as the discrete fractional component (DFC) method. This method is the discrete version of the CFC method, allowing for the rigorous property determination at the beginning and end states of the coupling parameter (λ = 1.0 or λ = 0.0), while also significantly enhancing the likelihood of inserting and deleting a large surfactant molecule in the 9794
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir aqueous phase.39 Like the work of Poursaeidesfahani and coworkers, a single molecule is coupled with λ at all times. This move allows us to efficiently sample the dilute concentration regime (or the Henry’s Law limit), in which case surfactant− surfactant interactions in the bulk are neglected. To our best knowledge, the DFC move and preceding variants of the move, have not yet been successfully applied to the case when a molecule with a moderate to large hydrophobic component is inserted into an aqueous phase.33,34,37,38 Similar to the work of Shen and Sun, we derive a simple thermodynamic relationship that incorporates a free energy of transfer to estimate the surfactant adsorption isotherm at various coverages. The determination of the free energy of transfer involves analysis of the Wang−Landau bias function iteratively generated during our simulations. The validation of this method has been discussed in our previous work for a simple hexane-water biphasic system.39 For the higher concentration regime, moves that can efficiently sample the aggregation behavior of surfactants in the bulk (i.e., estimation of the CMC) are required and should be the subject of future studies.
= 9 and m = 6. A harmonic bond angle potential was used for surfactant intramolecular angle interactions. Because the original framework of Cassandra was built on the basis of growth/regrowth algorithms of molecules via fragments libraries, we could only run our MC simulations with fixed bond lengths. Therefore, the model used in the present work does not account for bond flexibility in the surfactant. This results in small differences in properties, however, only at our highest coverage, compared to what one obtains from simulations of the formal SDK model. The effects of these constraints are discussed in more detail in the Supporting Information. 2.1. NVT Ensemble. Each system was initialized with a preequilibrated cubic water box and an orthorhombic water slab-vacuum (representative of air) interface box. The bulk water and slab-vacuum interface box dimensions were 4.16 × 4.16 × 4.16 nm3 and 4 × 4 × 14 nm3, respectively. Both boxes contained Nw = 800 water beads, and the cutoff radius was set to 1.5 nm. Varying numbers of C10E4 molecules with Ns = 10, 20, 30, 40, 50, 60, and 70 (corresponding to surfactant coverages of Γs = 0.31, 0.63, 0.94, 1.3, 1.6, 1.9, and 2.2 molecules/nm2) were placed symmetrically on both interfaces for each interface box, generating a total of seven independent systems. Note that Ns is the total number of surfactants in the system, and should be divided by two in order to obtain the number of surfactants per interface. Each interface box was then equilibrated in the canonical ensemble using 200 M MC steps, and another 300 M MC steps were used for the production run. The conventional translation, rotation, and also the configurational bias regrowth moves were used with equal attempt probabilities. We found that the slab consisting of 800 water beads was large enough such that surfactants on each side of the opposing monolayer could not interact with one another throughout the simulation. Furthermore, the system size was also large enough to provide accurate surface tension data for the majority of our systems. A more detailed investigation of system size effects on the surface tension is provided in the Supporting Information. Surface tensions were calculated using the equation
2. METHOD The simulations performed in this study were conducted using a development version of the open source Monte Carlo package, Cassandra, developed by our group.40,41 The algorithmic enhancements to Cassandra that were used to carry out our calculations will be released in a future version. The surfactant studied herein belongs to the class of single-tailed nonionic surfactants CxEy where x denotes the number of carbon atoms in the alkyl tail and y denotes the number of ethoxylated groups in the semipolar region of the surfactant headgroup. We used the SDK force field to model the surfactant, as this model not only reduces computational cost but can also provide accurate surface tensions and surface pressures as a function of surfactant coverage (see Figure 1B). 42 Figure 1A shows a
γ=
(2)
where γ is the surface tension, Lz is the box length in the normal direction, Pzz is the normal pressure tensor component, Pxx and Pyy are the components of the lateral pressure tensor, and the brackets denote the ensemble average. This interface box along with a pure water box were subsequently used as the initial configuration for the Gibbs ensemble simulations. 2.2. Gibbs Ensemble Method. The Gibbs ensemble was used to simulate a bulk and interface box simultaneously (Figure 1). To obtain efficient sampling in this ensemble, the DFC Monte Carlo move with Wang−Landau biasing was utilized (Figure 2). As mentioned earlier, this move significantly improves the efficiency of inserting and deleting long chain molecules of surfactants in both boxes. The DFC move is implemented by initially tagging or selecting a single surfactant molecule as “fractional” with an acceptance probability of
Figure 1. (a) Representation of the model used in the simulation. The SDK model for C10E4 and water implements a 3:1 mapping for each bead. (b) Snapshot of the simulation box of the Gibbs ensemble framework used in this work. Each box initial configuration is obtained by running a single NPT simulation for the bulk box and independent NVT simulations of an air−water interface box at various surfactant coverages. Water beads are represented by the blue isosurfaces. The bulk and interface box dimensions are roughly 4.16 × 4.16 × 4.16 nm3 and 4 × 4 × 14 nm3, respectively.
psel ∝
Ns, A Ns, B
exp(Δη) (3)
where Ns is the number of surfactants in one of two boxes, Δη is the difference in the Wang−Landau bias between neighboring λi states (e.g., Δη = ηi − ηi−1), and A and B denote either the bulk or interface box. When no surfactant molecules were present in the bulk box, a new surfactant molecule in the interface box was tagged as the new fractional molecule. For a fractional molecule, intermolecular interactions were scaled by the following relation
representation of the beads. The coarse-grained SDK model utilizes a 3:1 mapping of atoms to each bead, where intermolecular interactions are governed by the Mie potential defined by
⎛ n ⎞⎛ n ⎞m /(n − m) ⎡⎛ σ ⎞n ⎛ σ ⎞m⎤ ⎟⎜ ⎟ UMie(r ) = ⎜ ϵ⎢⎜ ⎟ − ⎜ ⎟ ⎥ ⎝r⎠ ⎦ ⎝ n − m ⎠⎝ m ⎠ ⎣⎝ r ⎠
⟨Pxx⟩+⟨Pyy⟩ ⎤ Lz ⎡ ⎥ ⎢⟨Pzz⟩− 2⎣ 2 ⎦
Uifrac(r ) = λiUMie(r ), Uifrac(r )
(1)
for
= λiUMie(r = 0.2 nm),
r ≥ 0.2 nm for
r < 0.2 nm
(4)
where λi ranges from 0.0 to 1.0 and i is the index corresponding to one of the values in λ = {1, 0.99, 0.98, ... 5 × 10−6, 5 × 10−6, ..., 0.98, 0.99, 1}. Our simulations included a total number of 283 λi values or
Here, n and m are the exponents for the repulsive and attractive terms of the interaction potential, respectively. For bead pairs involving water, n = 12 and m = 4. For surfactant−surfactant bead interactions, n 9795
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir
Figure 2. Schematic representation of the implemented Monte Carlo moves. The interfacial regions during the DFC swap move in step 3 (highlighted in red) are set by the Z-bounds obtained from the previous NVT simulations. Note that during this move, the surfactant headgroup fragment is first inserted into this region with the remaining molecular fragments of the surfactant grown using CBMC. subensembles, half of which correspond to the bulk box and the other half the interface box. Values of λi in the bulk and interface box are provided in the Supporting Information. Note that, for a pairwise distance less than 0.2 nm, the interparticle interaction energy was set to the value at r = 0.2 nm for all distances down to r = 0 nm. This is particularly important when λi approaches 0, as the overlapping of particles can cause large fluctuations in energy, which in turn can lead to very low acceptance probabilities in the Metropolis criteria and difficulty in converging the Wang−Landau bias function. Also note that intramolecular interactions for the fractional molecule are not scaled by λi. Once a molecule was tagged, the intermolecular interactions of the fractional molecule were scaled to a neighboring λi with an acceptance probability of ⎛ −ΔU ⎞ pscale ∝ exp⎜ + Δη⎟ ⎝ kBT ⎠
surfactant was swapped into the interface box, a surfactant headgroup was first inserted into one of the two specified regions, with the rest of the molecule subsequently grown using CBMC. Trial insertions were also implemented to increase the likelihood of insertion acceptance. The modification of the swap move ensured that surfactants were always swapped into the air−water interface while maintaining a conformation close to equilibrium. The Z-bounds were provided as input to the Gibbs ensemble simulations. Moreover, we rejected moves which caused the headgroup of the fractional molecule to move outside these bounds. Likewise, only surfactant molecules found within these bounds could be tagged as a fractional molecule in order to satisfy the condition of detailed balance. The surfactant molecules were not under these restraints as long as they were not tagged as a fractional molecule. An additional preferential bias translation move was implemented such that water beads neighboring a tagged molecule were preferentially translated.43,44 For each DFC move, a number of preferential translation moves were attempted for waters neighboring the fractional molecule within a spherical radius of 1 nm taken from the center of mass of the fractional surfactant. This ensured local thermal relaxation within the solvation shell of the fractional molecule, for each time the interparticle interactions of the fractional molecule was scaled to a different value. To satisfy detailed balance, moves that translated water outside the spherical radius were rejected. The probability of attempting each move was distributed by the following: DFC moves at 1% probability, biased translation moves at 35%, configurational bias regrowth moves at 9% for a fractional molecule and 9% for all other surfactant molecules, unbiased translation moves at 20%, rotation moves at 5%, dihedral moves at 5%, bond angle moves at 15%, and finally volume moves at 1%. The volume moves were only performed for the bulk box to maintain the pressure at 1 bar. Temperatures for all of the systems were set to 298 K. 2.3. Analysis of the Wang−Landau Bias. For each system, an independent equilibration run of 1500 M MC steps was performed to generate the Wang−Landau bias function. Each run was separated into three stages, with each stage consisting of 500 M MC steps. At the beginning of each stage, the configuration was reset to the original equilibrated conformation from the prior NVT simulation. This was to ensure that the process of generating or refining the bias function with a large bias factor did not perturb the system far away from its “equilibrium” state. In some instances, particularly for the low coverage systems, the starting bias factor was too large, causing an asymmetric
(5)
where ΔU is the difference in potential energy between two λ states. When λi reached its lowest value, either a scaling move to the greater neighboring λi value or a swap move was attempted such that the fractional surfactant molecule was swapped from one box to the other with an acceptance probability of
pswap ∝
⎛ − (ΔUA + ΔUB) ⎞ VA exp⎜ + Δη⎟ VB kBT ⎝ ⎠
(6)
where VA and VB are the bulk or interface volumes. These volumes are slightly modified as discussed below. Note that only a single fractional molecule was present throughout the entire simulation. For each fractional state or subensemble (λi), the Wang−Landau bias factor ηi was accumulated and used to improve the likelihood of walking across neighboring fractional states. We found that without the Wang− Landau biasing, the acceptance probability of walking across the entire span of λ was extremely low. 2.2.1. Additional Biasing Moves. The interbox swap moves of the DFC move were modified such that the volume of the interface box (VA or VB) was set to that of a rectangular region specified by the box widths in the X−Y direction (lateral to the interface) and a set of Zbounds for either interface. These Z-bounds were obtained by locating the position along the Z-direction corresponding to half the density maximum of the surfactant headgroup obtained from the previous NVT simulations (red region in Figure 2). For each time a fractional 9796
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir
Figure 3. Surface pressure (a) and surface tension (b) with surfactant coverage obtained from NVT simulations. The experimental value for the surface tension of water (red square) is shown in (b).
Figure 4. Surfactant and water density profile along the Z-axis for (a) low (Γs = 0.63 molecules/nm2), (b) intermediate (Γs = 1.3 molecules/nm2), and (c) high (Γs = 1.9 molecules/nm2) coverages. The corresponding snapshots for each system are shown above. Density profiles of the polar headgroup (yellow), surfactant ethoxy group (cyan), and alkyl tail (gray) have been plotted separately. Density profile for the water slab highlighted in light blue. Z-Distances at half the peak distance of the polar surfactant headgroup profiles indicated by the red arrows are used to obtain the bounds for the interfacial swap moves. distribution of surfactants at the two air−water interfaces in the interface box. Clearly, an asymmetric distribution would correspond to a nonequilibrium state in the slab box and likely cause a delay in the convergence of the algorithm. Hence, by resetting the simulation back to its original equilibrium conformation, we ensured that the Wang− Landau sampling was always performed near the equilibrium conformation of the slab. The starting Wang−Landau bias factor in all runs was set to 0.005, with a flatness criteria set to 0.75. The Wang−Landau bias factor was added to the bias function during each MC step. The simulation was ended when flatness was achieved after 9 to 10 Wang−Landau iterations corresponding to a bias tolerance ranging from 4.88 × 10−6 to 9.77 × 10−6. For the system at highest coverage (Γs = 2.2 molecules/nm2), we could only obtain 8 Wang− Landau iterations (corresponding to a bias tolerance of 3.91 × 10−6) within the number of MC steps. The criterion for convergence in the bias function is discussed in more detail in the Results and Discussion section. We have shown in our previous work that implementation of the Wang−Landau bias function can be easily analyzed to obtain both the
solvation free energy (ΔGsol) and the transfer free energy (ΔGtr) for an aqueous biphasic system (e.g., hexane and water box).39 In this system, the hexane rich phase and water rich phase were separated into two boxes where only a single solute molecule was transferred back and forth between the two phases. In this case, the Wang−Landau bias function equated to a free energy profile. If the bias function was converged, ΔGsol could be estimated by taking the difference in ηi values that correspond to λi = 1.0 (fully interacting) and λi = 0.0 (ideal gas) for the aqueous box. ΔGtr could also be obtained by taking the difference in values of ηi in the Wang−Landau bias function which correspond to λi = 1.0 in one box and λi = 1.0 in the other. In this work, we used the same methodology to obtain ΔGsol and ΔGtr for our surfactant systems, since near equilibrium, no more than a single surfactant swaps into the simulated bulk box.
3. RESULTS AND DISCUSSION 3.1. Effect of Surfactant Coverage on the Surface Tension. Simulations in the NVT ensemble were performed to 9797
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir
Figure 5. Subensembles with MC steps from the last 500 M MC steps of the production run. Subensembles are provided for systems at low (Γs = 0.63 molecules/nm2), intermediate (Γs = 1.3 molecules/nm2), and high (Γs = 1.9 molecules/nm2) coverages. A surfactant is swapped from the bulk to the interface when N goes from 1 to 283 contiguously, and vice versa for the opposite path.
swaps over the last 500 M steps is 15, 11, and 8 for the low, intermediate, and high coverages, respectively. (Note that these frequencies were obtained near the end of the production when the bias function was essentially converged, whereas much greater swap frequencies were obtained during the first 500 M steps of generating the bias function and the subsequent 500 M steps of the production run.) As expected, the frequency of surfactant swaps from the bulk to the interface tends to decrease with increasing coverage, as the likelihood of overlap between neighboring surfactants increases. Although the period for transfer can be quite large (∼30−60 M steps per complete swap), it is still significantly improved than when no Wang− Landau biasing is used or when only conventional MC moves including CBMC are used. Given the large number of λi values and low attempt probability of the DFC moves, these periods are quite reasonable. 3.3. Free Energy of Transfer and Estimation of Bulk Surfactant Concentration. The free energy profile along the reaction coordinate λ is shown in Figure 6. Swapping of the fractional surfactant between boxes occurs at λ141 and λ142, where the values for λi are equal and closest to zero. In the bulk box, the free energy barrier for turning on the intermolecular interactions for the fractional surfactant from the “ideal gas” (λ141) state is ∼50kBT along the reaction coordinate at λi ≈ 0.1. As λi is scaled down toward the ideal gas state (λ141), the free energy profile decreases sharply. This is because of the softening of interparticle interactions, which allows for a higher frequency of overlap between beads within the fractional surfactant molecule and neighboring water beads. Similarly, for the interface box, the free energy profile increases sharply as the interactions are turned on. Depending on the coverage, the free energy barrier for turning on the interactions ranges from ∼5− 20kBT at λi ≈ 0.02. This is consistent with the fact that the energy barrier is lower for inserting a molecule into a less crowded area. Qualitatively, the free energy barrier for the interface box along the reaction coordinate tends to increase with increasing surfactant coverage. This is expected since higher coverages should correspond to a greater likelihood of overlap or increased steric repulsion between surfactants. To estimate the surfactant concentrations in the bulk, we used a simple thermodynamic relation based on the free energy of transfer of surfactant from bulk into the interface. Here, we applied the same assumption as in the Gibbs adsorption isotherm, in that surfactant−surfactant interactions were neglected in the bulk (Henry’s law region) and the
obtain the surface pressure and surface tension as a function of surfactant coverage (Γs) (see Figure 3). The surface tension for the air−water interface when no surfactants are present (or surface tension of water) matches that of the experimental value and that of previous simulation work,31,45 and as expected, surface tension decreases monotonically with surfactant coverage. Figure 4 plots the density profiles along the Z-axis for systems including “low” (Γs = 0.63 molecules/nm2), “intermediate” (Γs = 1.3 molecules/nm2), and “high” (Γs = 1.9 molecules/nm2) surfactant coverage. The density profiles for polar head groups (yellow shaded regions), semipolar ethoxy groups (cyan shaded regions), and alkyl tails (gray shaded regions) have been plotted separately and averaged over the last 100 M MC steps with a total of 1000 frames. (Note that the initial water slab size of ∼4.5 nm in the Z-direction is sufficiently large such that surfactants from one side of the slab do not interact with those in the other side, given the cutoff radius of rcut = 1.5 nm.) The equilibrium distribution of the surfactant headgroup resides most closely at the hydrophilic region of the interface, while the alkyl tail resides at the hydrophobic (air) region. With increasing surfactant coverage, the density profiles for the surfactants broaden as does the density profile of water near the interface. For each simulation box, half the distance corresponding to the surfactant headgroup density profile peak was noted and used in the subsequent Gibbs ensemble simulations as input for the Zcoordinates in the swap moves. 3.2. Sampling in the Gibbs Ensemble. As mentioned previously, swapping a large molecule such as C10E4 can be quite challenging in the Gibbs ensemble without special biasing moves. For our systems, we find that, without the DFC move, no surfactants swap either from the bulk box into the interface box, or vice versa, even after a large number (>100 M) of MC steps. The difficulty in performing swaps is because of the large energy penalty for deleting a large molecule and creating a void space in a condensed phase such as liquid water, and also because of the difficulty in swapping the surfactant into a favorable region at the interface. Here, we show that we obtain reasonable swapping with implementation of the moves aforementioned. Figure 5 plots the visited subensemble i of λi during the last 500 M MC steps of our production run for the systems corresponding to the low, intermediate, and high coverages. A complete swap occurs when i moves contiguously from 1 (bulk) to 283 (interface). The total number of complete 9798
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir where μσid = μσ° + kBT ln(Γsas)
(9)
Above, μb° and μσ° are taken as the infinite dilution standard state chemical potentials with an ideal gas reference state. μidσ and μnon‑id are the ideal and nonideal chemical potential for an σ adsorbed surfactant molecule respectively, and kBT ln(Γsa) is the ideal entropic contribution to the chemical potential with as equal to the headgroup area of a surfactant molecule.30 Substituting eq 9 back into eq 8 yields μσ = μσ° + kBT ln(Γsas) + μσnon‐id (Γs)
(10)
μnon‑id σ
Note that is a term in the chemical potential that accounts for the nonideal interactions at the interface, which are dependent on both the surfactant coverage and the surface pressure. At equilibrium, the chemical potentials of the surfactants in the two boxes are equal such that μ b° + kBT ln(Xb) = μσ° + kBT ln(Γsas) + μσnon‐id (Γs) (11)
Figure 6. Free energy profile (ηi) vs subensembles λi. The vertical line indicates the value λi at which point swapping of surfactants occur between the bulk box (left side) and the interface box (right side). The free energy profile for the bulk box is approximately equal for all coverages, Γs in [nm−2], since a maximum of one surfactant is swapped into the box. The differences in ηi used to obtain ΔGsol/kBT and ΔGtr/ kBT are provided.
The surfactant concentration in the bulk was obtained from Xb = Γsas exp[− (μ b° − μσ° − μσnon‐id (Γs))/kBT ] = Γsas exp(ΔGtr /kBT )
where ΔGtr is the free energy of transfer of a single surfactant from the bulk to the interface. As can be seen from eq 12, ΔGtr = μσ° − μb° + μnon‑id (Γs), and is therefore coverage dependent. σ In our analysis, ΔGtr is the difference in value between the λ1 = 1 (fully interacting fractional molecule in the bulk box) and λ283 = 1 (fully interacting fractional molecule in the interface box). This has been shown in our previous work to be an accurate estimate for the free energy of transfer for hexane in water in an aqueous biphasic system.39 We set as = 0.22 nm2 as determined by Sresht et al. using MD simulations with the SDK model.14,30 A combined molecular simulation and MTT analysis on our surface tension results from NVT simulations also yielded a similar value (see the Supporting Information). Table 1 lists values of ΔGtr/kBT for each system at various coverages. The solvation free energy, ΔGsol of C10E4 in water can also be estimated as the difference in free energy at λi = 1 and λi = 0 in the bulk box, which corresponds to a fully interacting and a noninteracting (with intramolecular interactions fully turned on) C10E4 in water, respectively. We obtain ΔGsol = −14.7 ± 0.4 kBT, which is in good agreement with the observed trend in simulation solvation free energy calculations for nonionic surfactants obtained by Shen and Sun. We found good agreement in the values of ΔGsol obtained for each of the independently studied systems, which we used as an indicator for sufficient convergence in the Wang−Landau bias function.
concentration regime was assumed to be below the critical micelle concentration. This assumption was based on an experimental neutron reflectivity study of polyethylene glycol surfactant-air−water systems, which showed that the CMC occurs near the onset of saturation of surfactants at the interface.46 As can be seen from our own isotherm calculations (to be further discussed below), most of our MC simulations were performed at coverages below this saturation point. Independent MD simulations for systems at coverages near or above our highest coverage system showed no instances of surfactant desorption from the interface (see the Supporting Information). Still, it is expected that when the concentration reaches or goes beyond the CMC, the Gibbs adsorption isotherm assumption will break down. Hence, our analysis is only applicable below the CMC. We estimated the chemical potential for the surfactant in the bulk as μ b ≈ μ b° + kBT ln(Xb)
(7)
with Xb equal to the bulk mole fraction. The chemical potential for the surfactant in the interface is μσ = μσid + μσnon‐id
(12)
(8)
Table 1. Total Number of Surfactants, Surfactant Coverage, Surface Tension, Transfer Free Energy, and Bulk Surfactant Mole Fraction Ns 10 20 30 40 50 60 70
Γs [nm−2] 0.31 0.63 0.94 1.3 1.6 1.9 2.2
γ [mN/m] 70.5 66.9 62.9 61.7 56.6 51.2 40.7
± ± ± ± ± ± ±
ΔGtr/kBT −14.1 −14.8 −14.2 −13.2 −13.0 −12.4 −9.4
2.5 2.5 2.8 6.7 3.2 2.5 0.3 9799
± ± ± ± ± ± ±
0.5 0.5 0.5 0.5 0.5 0.5 0.5
Xb 5.3 5.1 1.5 5.2 7.7 1.6 3.9
× × × × × × ×
−8
10 10−8 10−7 10−7 10−7 10−6 10−5
± ± ± ± ± ± ±
3.4 3.3 9.3 3.3 4.9 1.0 2.5
× × × × × × ×
10−8 10−8 10−8 10−8 10−7 10−6 10−5
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir
also affected sampling for this particular coverage (see the Supporting Information for more detail on system size effects). Figure 8 shows the surface tension corresponding to each coverage with the bulk surfactant mole fractions. As mentioned
In our simulations, the bias function corresponding to the bulk box for all of the runs should have been nearly identical (as shown in Figure 6), since no more than a single surfactant molecule was present in the bulk during the production run. Furthermore, we found that for our lowest coverage system, running an additional 500 M MC steps did not change the ΔGsol as well as ΔGtr more than 0.1kBT, providing us additional confidence that the bias functions were close to or at convergence. (Note that the value for solvation free energy obtained from linear interpolation, as our lowest value for λi was not scaled down completely to 0 but rather to 5 × 10−6. A nonzero value λi was used to improve the swap acceptance rate, as it tends to generate more favorable trial positions in the CBMC method. This is not necessarily required if enough preferential biasing moves are performed, as was the case for our previously studied hexane-water system.39) Figure 7 shows surfactant coverage as a function of the bulk surfactant mole fraction and a comparison with the results from
Figure 8. Surface tension vs bulk surfactant mole fraction. Results from this work (blue markers) are plotted against values obtained using molecular thermodynamic theory (MTT) (solid black line), the combined MD/MTT approach (dashed black line), and experimental measurements (red markers) of Blankschtein and co-workers.14,30 The experimental CMC is indicated by the dotted vertical line. The open marker indicates that the system was subject to poor sampling.
previously, we assumed that at a particular coverage, the removal and addition of a single surfactant from bulk and interface did not perturb the surface tension significantly. While we can calculate the surface tension directly in the Gibbs Ensemble, we can only take statistics from states at which point the fractional molecule is at an integer λi state.37 Within the span of our simulation, we do not obtain enough integer states to obtain reasonable statistics. As seen in Figure 8, the predicted surfactant adsorption isotherm from our simulations is in good agreement with those predicted using the original MTT model of Blankschtein and co-workers (which was adjusted to an experimentally measured surface tension data point) as well as the recently modified MD/MTT approach.14,30 The only exception is for our datum at highest coverage, which is presumeably caused by the reasons previously mentioned. Note that, in both the MC approach and the MD/MTT approach, no experimental measurements are necessary to obtain the predicted isotherm. This agreement with both MTT and MD/MTT is quite remarkable, especially considering that no additional reparametrization which corrects for an overestimation of solvation free energies in the SDK model (as noted by Shen and Sun) was performed in this work, and also since the estimations of bulk concentrations are highly sensitive to ΔGtr.28 Still, the consistency in both our predicted bulk concentrations as a function of coverage and as a function of surface tension with the results of MD/MTT demonstrates that this methodology can be quite robust.
Figure 7. Surfactant coverage (Γs) vs bulk surfactant mole fraction at the air−water interface. The solid line represents the results obtained from the combined MD/MTT approach.30 The experimental CMC obtained is represented by the dotted vertical line.14 The open marker indicates that the system was possibly subject to poor sampling.
the combined MD/MTT approach recently developed by Blankschtein and co-workers.30 The error bars in the bulk mole fractions were determined by propagating the errors obtained in the solvation free energy calculations ΔGsol for each of the seven coverages examined. Coverage increases monotonically with bulk mole fraction, except at the lowest coverage of Γs = 0.31 molecules/nm2. Still, the curve predicted by the MD/ MTT approach lies within the errors of the majority our predicted bulk mole fractions. For our system, at highest coverage (Γs = 2.2 molecules/nm2), the bulk concentration begins to shift away from the semilog−linear trend as well as from the results of MD/MTT. Furthermore, at this coverage, the estimation of bulk concentration is above the experimental CMC. We suspect that this discrepancy was caused by poor sampling, as evident from the fewer number of Wang−Landau iterations attained (8 iterations) than compared to those for the rest of the systems (9 or 10 iterations). Moreover, finite system size effects and our imposed fixed bond constraints could have 9800
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir Notes
4. CONCLUSIONS We have presented a new Monte Carlo scheme that can be used to simultaneously predict the bulk surfactant concentrations and surfactant coverages at a fluid−air interface in the dilute limit. The algorithm incorporates the recently developed DFC move, which for the first time is successfully implemented for a relatively large amphiphilic molecule in a condensed aqueous medium. In the limit where no more than one surfactant is present in the bulk box, we have derived a simple thermodynamic relationship based on the free energy of transfer for a single surfactant system, which can be used to obtain bulk surfactant mole fractions or concentrations based on a converged Wang−Landau bias function. This approach enables the direct computation of bulk and interfacial concentrations without making any assumptions to the equation of state, which is particularly important when the equation of state is not known. The results presented also demonstrate that this methodology can be used to estimate bulk concentrations that are in excellent agreement with both MTT as well as the combined MD/MTT approach. In principle, this methodology can be further extended to systems containing dilute surfactant mixtures, charged surfactants (provided swap moves can maintain charge neutrality in both boxes), as well as fluid−fluid interfaces (e.g., water−oil interfaces). For higher concentrations (i.e., when it is more favorable to have more than one surfactant in the bulk box), the Gibbs− DFCMC framework should provide the bulk concentrations directly, without any additional analysis in the Wang−Landau bias function. However, sampling of surfactant−surfactant behavior accurately in the bulk is also quite challenging, as these interactions are often associated with self-assembly phenomena. Moves that can efficiently sample this phenomena are further necessary to accurately obtain bulk concentrations near or beyond the CMC. Some possible examples of such moves that can be used in conjunction with our scheme might include the aggregation bias volume Monte Carlo (AVBMC) and its variants,47−50 the multiparticle-move,51,52 and hybrid MD moves.
■
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS Funding for this work was provided by ExxonMobil Research and Engineering. The authors thank Dr. Daniel Cherney at ExxonMobil for helpful discussions. Computational resources were provided by Notre Dame’s Center for Research Computing. No competing financial interest was provided in this work.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.7b02058. MTT theory, MC/MTT fitting procedure, surface tension versus surfactant coverage from the MC/MTT model, effects of system size and bond stiffness on surface tension, and intermolecular energy scaling values for each DFC subensemble (PDF)
■
REFERENCES
(1) Siepmann, J.; Karaborni, S.; Klein, M. Monte Carlo simulation of the liquid-vapor coexistence in a Langmuir monolayer of pentadecanoic acid. J. Phys. Chem. 1994, 98, 6675−6678. (2) Barton, S.; Thomas, B.; Flom, E.; Rice, S. A.; Lin, B.; Peng, J.; Ketterson, J.; Dutta, P. X-ray diffraction study of a Langmuir monolayer of C21H43OH. J. Chem. Phys. 1988, 89, 2257−2270. (3) Kaganer, V. M.; Möhwald, H.; Dutta, P. Structure and phase transitions in Langmuir monolayers. Rev. Mod. Phys. 1999, 71, 779. (4) Maoz, R.; Sagiv, J. On the formation and structure of selfassembling monolayers. I. A comparative atr-wettability study of Langmuir-Blodgett and adsorbed films on flat substrates and glass microbeads. J. Colloid Interface Sci. 1984, 100, 465−496. (5) Ulman, A. An Introduction to Ultrathin Organic Films: From Langmuir-Blodgett to Self-Assembly; Academic Press, 2013. (6) Shih, C.-J.; Lin, S.; Sharma, R.; Strano, M. S.; Blankschtein, D. Understanding the pH-dependent behavior of graphene oxide aqueous solutions: a comparative experimental and molecular dynamics simulation study. Langmuir 2012, 28, 235−241. (7) Bergeron, V.; Langevin, D.; Asnacios, A. Thin-film forces in foam films containing anionic polyelectrolyte and charged surfactants. Langmuir 1996, 12, 1550−1556. (8) Murray, B. S.; Nelson, P. V. A novel Langmuir trough for equilibrium and dynamic measurements on air-water and oil-water monolayers. Langmuir 1996, 12, 5973−5976. (9) Aveyard, R.; Clint, J. H.; Nees, D.; Paunov, V. N. Compression and structure of monolayers of charged latex particles at air/water and octane/water interfaces. Langmuir 2000, 16, 1969−1979. (10) Aveyard, R.; Binks, B.; Fletcher, P. Interfacial tensions and aggregate structure in pentaethylene glycol monododecyl ether/oil/ water microemulsion systems. Langmuir 1989, 5, 1210−1217. (11) Fraaije, J. G.; Tandon, K.; Jain, S.; Handgraaf, J.-W.; Buijse, M. Method of moments for computational microemulsion analysis and prediction in tertiary oil recovery. Langmuir 2013, 29, 2136−2151. (12) Eastoe, J.; Dalton, J. Dynamic surface tension and adsorption mechanisms of surfactants at the air-water interface. Adv. Colloid Interface Sci. 2000, 85, 103−144. (13) Chang, C.-H.; Franses, E. I. Adsorption dynamics of surfactants at the air/water interface: a critical review of mathematical models, data, and mechanisms. Colloids Surf., A 1995, 100, 1−45. (14) Nikas, Y.; Puvvada, S.; Blankschtein, D. Surface tensions of aqueous nonionic surfactant mixtures. Langmuir 1992, 8, 2680−2689. (15) Mulqueen, M.; Blankschtein, D. Prediction of equilibrium surface tension and surface adsorption of aqueous surfactant mixtures containing ionic surfactants. Langmuir 1999, 15, 8832−8848. (16) Mulqueen, M.; Blankschtein, D. Prediction of equilibrium surface tension and surface adsorption of aqueous surfactant mixtures containing zwitterionic surfactants. Langmuir 2000, 16, 7640−7654. (17) Theodorou, D. N. A reversible minimum-to-minimum mapping method for the calculation of free-energy differences. J. Chem. Phys. 2006, 124, 034109. (18) Howes, A.; Radke, C. Monte Carlo simulations of LennardJones nonionic surfactant adsorption at the liquid/vapor interface. Langmuir 2007, 23, 1835−1844. (19) Howes, A.; Radke, C. Monte Carlo simulation of mixed Lennard-Jones nonionic surfactant adsorption at the liquid/vapor interface. Langmuir 2007, 23, 11580−11586.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Brian Yoo: 0000-0002-0326-0831 Ryan Gotchy Mullen: 0000-0002-2099-6955 Edward J. Maginn: 0000-0002-6309-1347 Present Address §
B.Y.: BASF Corporation, 540 White Plains Road, Tarrytown, NY 10591-9005. 9801
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
Article
Langmuir (20) Frenkel, D.; Mooij, G.; Smit, B. Novel scheme to study structural and thermal properties of continuously deformable molecules. J. Phys.: Condens. Matter 1992, 4, 3053. (21) Tomassone, M.; Couzis, A.; Maldarelli, C.; Banavar, J.; Koplik, J. Phase Transitions of Soluble Surfactants at a Liquid- Vapor Interface. Langmuir 2001, 17, 6037−6040. (22) Tomassone, M.; Couzis, A.; Maldarelli, C.; Banavar, J.; Koplik, J. Molecular dynamics simulation of gaseous-liquid phase transitions of soluble and insoluble surfactants at a fluid interface. J. Chem. Phys. 2001, 115, 8634−8642. (23) Rekvig, L.; Kranenburg, M.; Hafskjold, B.; Smit, B. Effect of surfactant structure on interfacial properties. Europhys. Lett. 2003, 63, 902. (24) Rekvig, L.; Kranenburg, M.; Vreede, J.; Hafskjold, B.; Smit, B. Investigation of surfactant efficiency using dissipative particle dynamics. Langmuir 2003, 19, 8195−8205. (25) Groot, R. D. Mesoscopic simulation of polymer-surfactant aggregation. Langmuir 2000, 16, 7493−7502. (26) Groot, R. D.; Rabone, K. Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophys. J. 2001, 81, 725−736. (27) Siepmann, J. I.; Frenkel, D. Configurational bias Monte Carlo: a new sampling scheme for flexible chains. Mol. Phys. 1992, 75, 59−70. (28) 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. (29) 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. (30) Sresht, V.; Lewandowski, E. P.; Blankschtein, D.; Jusufi, A. Combined Molecular Dynamics Simulation-Molecular Thermodynamic Theory Framework for Predicting Surface Tensions. Langmuir 2017, 33, 8319−8329. (31) 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. (32) Panagiotopoulos, A.; Quirke, N.; Stapleton, M.; Tildesley, D. Phase equilibria by simulation in the Gibbs ensemble: alternative derivation, generalization and application to mixture and membrane equilibria. Mol. Phys. 1988, 63, 527−545. (33) Shi, W.; Maginn, E. J. Continuous fractional component Monte Carlo: an adaptive biasing method for open system atomistic simulations. J. Chem. Theory Comput. 2007, 3, 1451−1463. (34) Shi, W.; Maginn, E. J. Improvement in molecule exchange efficiency in Gibbs ensemble Monte Carlo: development and implementation of the continuous fractional component move. J. Comput. Chem. 2008, 29, 2520−2530. (35) Landau, D.; Wang, F. Determining the density of states for classical statistical models by a flat-histogram random walk. Comput. Phys. Commun. 2002, 147, 674−677. (36) Wang, F.; Landau, D. Efficient, multiple-range random walk algorithm to calculate the density of states. Phys. Rev. Lett. 2001, 86, 2050. (37) Poursaeidesfahani, A.; Torres-Knoop, A.; Dubbeldam, D.; Vlugt, T. J. Direct free energy calculation in the continuous fractional component Gibbs ensemble. J. Chem. Theory Comput. 2016, 12, 1481− 1490. (38) Poursaeidesfahani, A.; Rahbari, A.; Torres-Knoop, A.; Dubbeldam, D.; Vlugt, T. J. Computation of thermodynamic properties in the continuous fractional component Monte Carlo Gibbs ensemble. Mol. Simul. 2017, 43, 189−195. (39) Marin-Rimoldi, E. Monte Carlo simulations for phase equilibria and software development. Ph.D. Dissertation, University of Notre Dame, Notre Dame, IN, 2017. (40) Shah, J. K.; Maginn, E. J. A general and efficient Monte Carlo method for sampling intramolecular degrees of freedom of branched and cyclic molecules. J. Chem. Phys. 2011, 135, 134121.
(41) Shah, J. K.; Marin-Rimoldi, E.; Mullen, R. G.; Keene, B. P.; Khan, S.; Paluch, A. S.; Rai, N.; Romanielo, L. L.; Rosch, T. W.; Yoo, B.; Maginn, E. J. Cassandra: An open source Monte Carlo package for molecular simulation. J. Comput. Chem. 2017, 38, 1727−1739. (42) Shinoda, W.; Discher, D. E.; Klein, M. L.; Loverde, S. M. Probing the structure of PEGylated-lipid assemblies by coarse-grained molecular dynamics. Soft Matter 2013, 9, 11549−11556. (43) Rosch, T. W.; Maginn, E. J. Reaction ensemble Monte Carlo simulation of complex molecular systems. J. Chem. Theory Comput. 2011, 7, 269−279. (44) Owicki, J.; Scheraga, H. Preferential sampling near solutes in Monte Carlo calculations on dilute solutions. Chem. Phys. Lett. 1977, 47, 600−602. (45) Pallas, N.; Harrison, Y. An automated drop shape apparatus and the surface tension of pure water. Colloids Surf. 1990, 43, 169−194. (46) Li, P. X.; Li, Z. X.; Shen, H.-H.; Thomas, R. K.; Penfold, J.; Lu, J. R. Application of the Gibbs equation to the adsorption of nonionic surfactants and polymers at the air-water interface: comparison with surface excesses determined directly using neutron reflectivity. Langmuir 2013, 29, 9324−9334. (47) Chen, B.; Siepmann, J. I. A novel Monte Carlo algorithm for simulating strongly associating fluids: applications to water, hydrogen fluoride, and acetic acid. J. Phys. Chem. B 2000, 104, 8725−8734. (48) Chen, B.; Siepmann, J. I. Improving the Efficiency of the Aggregation- Volume- Bias Monte Carlo Algorithm. J. Phys. Chem. B 2001, 105, 11275−11282. (49) Chen, B.; Siepmann, J. I.; Oh, K. J.; Klein, M. L. Aggregationvolume-bias Monte Carlo simulations of vapor-liquid nucleation barriers for Lennard-Jonesium. J. Chem. Phys. 2001, 115, 10903− 10913. (50) Loeffler, T. D.; Sepehri, A.; Chen, B. Improved Monte Carlo scheme for efficient particle transfer in heterogeneous systems in the grand canonical ensemble: Application to vapor-liquid nucleation. J. Chem. Theory Comput. 2015, 11, 4023−4032. (51) Moučka, F.; Rouha, M.; Nezbeda, I. Efficient multiparticle sampling in Monte Carlo simulations on fluids: Application to polarizable models. J. Chem. Phys. 2007, 126, 224106. (52) Moučka, F.; Nezbeda, I. Multi-particle sampling in Monte Carlo simulations on fluids: Efficiency and extended implementations. Mol. Simul. 2009, 35, 660−672.
9802
DOI: 10.1021/acs.langmuir.7b02058 Langmuir 2017, 33, 9793−9802
