Discrete Fractional Component Monte Carlo Simulation Study of Dilute

Corporate Strategic Research, ExxonMobil Research and Engineering Company, 1545 U.S. 22, Annandale, New Jersey 08801, United States. Langmuir , 2017, ...
0 downloads 12 Views 4MB Size
Subscriber access provided by University of Texas Libraries

Article

A Discrete Fractional Component Monte Carlo Simulation Study of Dilute Nonionic Surfactants at the Air-Water Interface Brian Yoo, Eliseo Marin-Rimoldi, Ryan G. Mullen, Arben Jusufi, and Edward J. Maginn Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.7b02058 • Publication Date (Web): 28 Aug 2017 Downloaded from http://pubs.acs.org on August 29, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

A 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, IN 46556-5637 ‡Corporate Strategic Research, ExxonMobil Research and Engineering Company, 1545 U.S. 22, Annandale, NJ 08801, USA ¶Current address: BASF Corporation, 540 White Plains Road, Tarrytown, NY 10591-9005 E-mail: [email protected]

1 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

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 regime of these systems of interest are typically very dilute ( + < Pyy > γ= < Pzz > − 2 2

(2)

where γ is the surface tension, Lz is the box length, 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.

9 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 33

2.2 Gibbs Ensemble Method The Gibbs ensemble was used to simulate a bulk and interface box simultaneously (Fig. 1). To obtain efficient sampling in this ensemble, the DFC Monte Carlo move with WangLandau biasing was utilized (Fig. 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

psel ∝

Ns,A exp (∆η) , Ns,B

(3)

where Ns is the number of surfactants in one of two boxes, ∆η is difference in the WangLandau 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 Uif rac (r) = λi UM ie (r), Uif rac (r) = λi UM ie (r = 2.0 Å),

f or r >= 2.0 Å

(4)

f or r < 2.0 Å.

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, ... 5E-6, 5E-6, ..., 0.98, 0.99, 1}. Our simulations included a total number of 283 λi values or 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 2 Å, the interparticle interaction energy was set to the value at r = 2 Å for all distances down to r = 0 Å. 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 10 ACS Paragon Plus Environment

Page 11 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 33

one box to the other with an acceptance probability of

pswap

VA ∝ exp VB



 −(∆UA + ∆UB ) + ∆η . kB T

(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 Z-bounds 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 head group obtained from the previous NVT simulations (red region in Fig. 2). For each time a fractional surfactant was swapped into the interface box, a surfactant head group 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 insertions. 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 head group 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. 12 ACS Paragon Plus Environment

Page 13 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

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 10 Å 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 10% for a fractional molecule and 10% for all other surfactant molecules, unbiased translation moves at 20%, rotation moves at 5%, dihedral moves at 15%, 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 1500M MC steps was performed to generate the Wang-Landau bias function. Each run was separated into three stages, with each stage consisting of 500M 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 distribution of surfactants at the two air-water interfaces in the interface box. Clearly, an asymmetric distribution would correspond to a non-equilibrium state in the slab box and likely cause a delay in the convergence of the algorithm. Hence, by resetting

13 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the simulation back to its original equilibrium conformation, we ensured that the WangLandau 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.88x10−6 to 9.77x10−6 . For the system at highest coverage (Γs = 2.2 molec/nm2 ), we could only obtain 8 Wang-Landau iterations (corresponding to a bias tolerance of 3.91x10−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.

14 ACS Paragon Plus Environment

Page 14 of 33

Page 15 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 16 of 33

Page 17 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

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 Fig. 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 ∼50 kB T along the reaction coordinate at λi ≈0.1. As λi is scaled down towards 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–20 kB T 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 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

18 ACS Paragon Plus Environment

Page 18 of 33

Page 19 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 33

We estimated the chemical potential for the surfactant in the bulk as

µb ≈ µob + kB T ln(Xb ),

(7)

with Xb equal to the bulk mole fraction. The chemical potential for surfactant in the interface is non−id µσ = µid σ + µσ

(8)

o µid σ = µσ + kB T ln(Γs as ).

(9)

where

Above, µob and µoσ are taken as the infinite dilution standard state chemical potentials with non−id an ideal gas reference state. µid are the ideal and non-ideal chemical potential σ and µσ

for an adsorbed surfactant molecule respectively, and kB T ln(Γs a) is the ideal entropic contribution to the chemical potential with as equal to the head-group area of a surfactant molecule. 30 Substituting Eq. (9) back into Eq. (8) yields

µσ = µoσ + kB T ln(Γs as ) + µσnon−id (Γs ).

(10)

Note that µσnon−id is a term in the chemical potential that accounts for the non-ideal 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

µob + kB T ln(Xb ) = µoσ + kB T ln(Γs as ) + µσnon−id (Γs ).

(11)

The surfactant concentration in the bulk was obtained from

  Xb = Γs as exp −(µob − µoσ − µσnon−id (Γs ))/kB T = Γs as exp (−∆Gtr /kB T ) 20 ACS Paragon Plus Environment

(12)

Page 21 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

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 = µob − µoσ − µσ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 Supporting Information). Table 1 lists values of -∆Gtr /kB T for each system at various coverages. The solvation free energy, ∆Gsol of C10 E4 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 non-interacting (with intramolecular interactions fully turned on) C10 E4 in water, respectively. We obtain ∆Gsol =-14.7±0.4 kB T, 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 a good agreement in values of ∆Gsol obtained for each of the independently studied systems which we used as an indicator that the Wang-Landau bias function for each system was sufficiently converged. In our simulations, the bias function corresponding to the bulk box for all of the runs should have been nearly identical (as shown in Fig. 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 500M MC steps did not change the ∆Gsol =-14.7±0.4 kB T as well as ∆Gtr more than 0.1 kB T, 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 5E-6. A non-zero value λi was used to improve the swap acceptance rate, as the 21 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 33

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−1 70.5 ± 2.5 66.9 ± 2.5 62.9 ± 2.8 61.7 ± 6.7 56.6 ± 3.2 51.2 ± 2.5 40.7 ± 0.3

-∆Gtr /kB T

Xb

-14.1 ± 0.5 -14.8 ± 0.5 -14.2 ± 0.5 -13.2 ± 0.5 -13.0 ± 0.5 -12.4 ± 0.5 -9.4 ± 0.5

5.3E-8 ± 3.4E-8 5.1E-8 ± 3.3E-8 1.5E-7 ± 9.3E-8 5.2E-7 ± 3.3E-8 7.7E-7 ± 4.9E-7 1.6E-6 ± 1.0E-6 3.9E-5 ± 2.5E-5

nonzero value for λi generates 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 ) Fig. 7 shows surfactant coverage as a function of the bulk surfactant mole fraction and a comparison with the results from the combined MD/MTT approach recently developed by Blankschtein and coworkers. 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 molec/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 molec/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 poorly 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 also affected sampling for this particular coverage (see Supporting Information for more detail on system size

22 ACS Paragon Plus Environment

Page 23 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 24 of 33

Page 25 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

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.

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 surfac-

25 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tant 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 multi-particle-move, 51,52 and hybrid MD moves.

Acknowledgement 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.

Supporting Information Available 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. This material is available free of charge via the Internet at http://pubs.acs.org/.

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. 26 ACS Paragon Plus Environment

Page 26 of 33

Page 27 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(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 self-assembling 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 2011, 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. 27 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(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 Lennard-Jones 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. (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. 28 ACS Paragon Plus Environment

Page 28 of 33

Page 29 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(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.

29 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31) Shinoda, W.; DeVane, R.; Klein, M. L. Multi-property fitting and parameterization of a coarse grained model for aqueous surfactants. Mol. Sim. 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. Sim. 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.

30 ACS Paragon Plus Environment

Page 30 of 33

Page 31 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(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.

31 ACS Paragon Plus Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(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. Aggregation-volume-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ˇcka, 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ˇcka, F.; Nezbeda, I. Multi-particle sampling in Monte Carlo simulations on fluids: Efficiency and extended implementations. Mol. Sim. 2009, 35, 660–672.

32 ACS Paragon Plus Environment

Page 32 of 33

Page 33 of 33

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment