Thermodynamics and partitioning of homopolymers into a slit-A grand canonical Monte Carlo simulation study

Grand canonical ensemble Monte Carlo simulation (GCMC) combined with the histogram reweighting technique was used to study the thermodynamic equilibrium of a homopolymer solution between a bulk and a slit pore. GCMC gives the partition coefficients that agree with those from canonical ensemble Monte Carlo simulations in a twin box, and it also gives results that are not accessible through the regular canonical ensemble simulation such as the osmotic pressure of the solution. In a bulk polymer solution, the calculated osmotic pressure agrees very well with the scaling theory predictions both for the athermal polymer solution and the theta solution. However, one cannot obtain the osmotic pressure of the confined solution in the same way since the osmotic pressure of the confined solution is anisotropic. The chemical potentials in GCMC simulations were found to differ by a translational term from the chemical potentials obtained from canonical ensemble Monte Carlo simulations with the chain insertion method. This confirms the equilibrium condition of a polymer solution partition between the bulk and a slit pore: the chemical potentials of the polymer chain including the translational term are equal at equilibrium. The histogram reweighting method enables us to obtain the partition coefficients in the whole range of concentrations based on a limited set of simulations. Those predicted bulk-pore partition coefficient data enable us to perform further theoretical analysis. Scaling predictions of the partition coefficient at different regimes were given and were confirmed by the simulation data.     

Thermodynamics of the liquid states of Langmuir monolayers

The liquid states and the liquid-liquid equilibrium of surfactant molecules forming an interphase between air and water have been considered using Monte Carlo computer simulations. Specifically, the expanded and compressed liquid phases observed for surfactant molecules were characterized as a function of pressure and temperature. Simple modified beadlike potentials were implemented in order to describe the interparticle forces between the hydrophobic and hydrophilic portions of surfactant molecules at the air/water interface. A simulation box was defined such that the monolayer was exposed to an externally applied lateral pressure in a modified isothermal-isobaric ensemble, whereas the water bath was modeled in a canonical ensemble. The simulation resembles the experimental setup used to measure lateral pressure (Pi) versus area isotherms obtained with Langmuir troughs. The applied lateral pressure-surface area phase diagram clearly showed the coexistence of the expanded and compressed liquid phases within certain temperature and pressure ranges. Distribution functions of distances and enthalpies for the monolayer were computed to clearly identify each liquid phase and the coexistence region.     

Simulation of chemical potentials and phase equilibria in two- and three-dimensional square-well fluids: finite size effects

We study the simulation cell size dependence of chemical potential isotherms in subcritical square-well fluids by means of series of canonical ensemble Monte Carlo simulations with increasing numbers of particles, for both three-dimensional bulk systems and two-dimensional planar layers, using Widom-like particle insertion methods. By estimating the corresponding vapor/liquid coexistence densities using a Maxwell-like equal area rule for the subcritical chemical potential isotherms, we are able to study the influence of system size not only on chemical potentials but also on the coexistence properties. The chemical potential versus density isotherms show van der Waals-like loops in the subcritical vapor/liquid coexistence range that exhibit distinct finite size effects for both two- and three-dimensional fluids. Generally, in agreement with recent findings for related studies of Lennard-Jones fluids, the loops shrink with increasing number of particles. In contrast to the subcritical isotherms themselves, the equilibrium vapor/liquid densities show only a weak system size dependence and agree quantitatively with the best-known literature values for three-dimensional fluids. This allows our approach to be used to accurately predict the phase coexistence properties. Our resulting phase equilibrium results for two-dimensional square-well fluids are new. Knowledge concerning finite size effects of square-well systems is important not only for the simulation of thermodynamic properties of simple fluids, but also for the simulation of models of more complex fluids (such as aqueous or polymer fluids) involving square-well interactions.     

Transferability of coarse-grained force fields: the polymer case

A key question for all coarse-graining methodologies is the degree of transferability of the resulting force field between various systems and thermodynamic conditions. Here we present a detailed study of the transferability over different thermodynamic states of a coarse-grained (CG) force field developed using the iterative Boltzmann inversion method. The force field is optimized against distribution functions obtained from atomistic simulations. We analyze the polymer case by investigating the bulk of polystyrene and polyamide-6,6 whose coarse-grained models differ in the chain length and in the number of atoms lumped in one bead. The effect of temperature and pressure on static, dynamic, and thermodynamic properties is tested by comparing systematically the coarse-grain results with the atomistic ones. We find that the CG model describing the polystyrene is transferable only in a narrow range of temperature and it fails in describing the change of the bulk density when temperature is 80 K lower than the optimization one. Moreover the calculation of the self-diffusion coefficient shows that the CG model is characterized by a faster dynamics than the atomistic one and that it overestimates the isothermal compressibility. On the contrary, the polyamide-6,6 CG model turns out to be fully transferable between different thermodynamic conditions. The transferability is checked by changing either the temperature or the pressure of the simulation. We find that, in this case, the CG model is able to follow all the intra- and interstructural rearrangements caused by the temperature changes. In addition, while at low temperature the difference between the CG and atomistic dynamics is remarkable due to the presence of hydrogen bonds in the atomistic systems, for high temperatures, the speedup of the CG dynamics is strongly reduced, leading to a CG diffusion coefficient only six times bigger than the atomistic one. Moreover, the isothermal compressibility calculated at different temperatures agrees very well with the experimental one. We find that the polymer chain length does not affect the transferability of the force field and we attribute such transferability mainly to the finer model used in describing the polyamide-6,6 than the polystyrene.     

Solubility of NaCl in water by molecular simulation revisited

In this paper, the solubility of NaCl in water is evaluated by using computer simulations for three different force fields. The condition of chemical equilibrium (i.e., equal chemical potential of the salt in the solid and in the solution) is obtained at room temperature and pressure to determine the solubility of the salt. We used the same methodology that was described in our previous work [E. Sanz and C. Vega, J. Chem. Phys. 126, 014507 (2007)] although several modifications were introduced to improve the accuracy of the calculations. It is found that the predictions of the solubility are quite sensitive to the details of the force field used. Certain force fields underestimate the experimental solubility of NaCl in water by a factor of four, whereas the predictions of other force fields are within 20% of the experimental value. Direct coexistence molecular dynamic simulations were also performed to determine the solubility of the salt. Reasonable agreement was found between the solubility obtained from free energy calculations and that obtained from direct coexistence simulations. This work shows that the evaluation of the solubility of salts in water can now be performed in computer simulations. The solubility depends on the ion-ion, ion-water, and water-water interactions. For this reason, the prediction of the solubility can be quite useful in future work to develop force fields for ions in water.     

Molecular dynamics simulation in the grand canonical ensemble

An extended system Hamiltonian is proposed to perform molecular dynamics (MD) simulation in the grand canonical ensemble. The Hamiltonian is similar to the one proposed by Lynch and Pettitt (Lynch and Pettitt, J Chem Phys 1997, 107, 8594), which consists of the kinetic and potential energies for real and fractional particles as well as the kinetic and potential energy terms for material and heat reservoirs interacting with the system. We perform a nonlinear scaling of the potential energy parameters of the fractional particle, as well as its mass to vary the number of particles dynamically. On the basis of the equations of motion derived from this Hamiltonian, an algorithm has been proposed for MD simulation at constant chemical potential. The algorithm has been tested for the ideal gas, for the Lennard-Jones fluid over a wide range of temperatures and densities, and for water. The results for the low-density Lennard-Jones fluid are compared with the predictions from a truncated virial equation of state. In the case of the dense Lennard-Jones fluid and water our predicted results are compared with the results reported using other available methods for the calculation of the chemical potential. The method is also applied to the case of vapor-liquid coexistence point predictions.     

Dissipative particle dynamics simulations in the grand canonical ensemble: applications to polymer brushes.

We have used the dissipative particle dynamics (DPD) method in the grand canonical ensemble to study the compression of grafted polymer brushes in good solvent conditions. The force-distance profiles calculated from DPD simulations in the grand canonical ensemble are in very good agreement with the self-consistent field (SCF) theoretical models and with experimental results for two polystyrene brush layers grafted onto mica surfaces in toluene.     

Isomolar semigrand ensemble molecular dynamics: development and application to liquid-liquid equilibria

An extended system molecular dynamics method for the isomolar semigrand ensemble (fixed number of particles, pressure, temperature, and fugacity fraction) is developed and applied to the calculation of liquid-liquid equilibria (LLE) for two Lennard-Jones mixtures. The method utilizes an extended system variable to dynamically control the fugacity fraction xi of the mixture by gradually transforming the identity of particles in the system. Two approaches are used to compute coexistence points. The first approach uses multiple-histogram reweighting techniques to determine the coexistence xi and compositions of each phase at temperatures near the upper critical solution temperature. The second approach, useful for cases in which there is no critical solution temperature, is based on principles of small system thermodynamics. In this case a coexistence point is found by running N-P-T-xi simulations at a common temperature and pressure and varying the fugacity fraction to map out the difference in chemical potential between the two species A and B (mu(A)-mu(B)) as a function of composition. Once this curve is known the equal-distance/equal-area criterion is used to determine the coexistence point. Both approaches give results that are comparable to those of previous Monte Carlo (MC) simulations. By formulating this approach in a molecular dynamics framework, it should be easier to compute the LLE of complex molecules whose intramolecular degrees of freedom are often difficult to properly sample with MC techniques.     

Reverse nonequilibrium molecular dynamics simulation of thermal conductivity in nanoconfined polyamide-6,6

A new molecular dynamics simulation method, with coupling to external baths, is used to perform equilibrium simulations on polyamide-6,6 trimers nanoconfined between graphene surfaces, in equilibrium with the bulk polymer. The method is coupled with the reverse nonequilibrium molecular dynamics simulation technique to exchange heat in the direction normal to the surfaces. To be able to study the effect of confinement on the heat conductance in nanoconfined pores, in this work a number of simulations on systems with different pore sizes are done. It is concluded that the coefficient of heat conductivity depends on the degree of polymer layering between the surfaces and on the pore width. Our results further indicate a considerable temperature drop at the interface between the surfaces and polymer. The calculated Kapitza lengths depend on the intersurface distance and on the layering of the polymer nanoconfined between the surfaces.     

A finite-size scaling study of a model of globular proteins

Grand canonical Monte Carlo simulations are used to explore the metastable fluid-fluid coexistence curve of the modified Lennard-Jones model of globular proteins of ten Wolde and Frenkel [Science, 277, 1975 (1997)]. Using both mixed-field finite-size scaling and histogram-reweighting methods, the joint distribution of density and energy fluctuations is analyzed at coexistence to accurately determine the critical-point parameters. The subcritical coexistence region is explored using the recently developed hyper parallel tempering Monte Carlo simulation method along with histogram reweighting to obtain the density distributions. The phase diagram for the metastable fluid-fluid coexistence curve is calculated in close proximity to the critical point, a region previously unattained by simulations.     

Monte Carlo molecular simulation of the hydration of K-montmorillonite at 353 K and 625 bar

Monte Carlo molecular simulations of the hydration of K-saturated Wyoming-type montmorillonite at constant stress in the NPzzT ensemble and at constant chemical potential in the grand canonical muVT ensemble, under basin-like conditions of 353 K and 625 bar, show a strong tendency of the K+ ions to adhere to the siloxane surface, forming predominant inner-sphere complexes with tetrahedral oxygen atoms and adsorbed water molecules. Simulations in the grand canonical ensemble predict that none of the K-montmorillonite hydrates, the one-, two-, and three-layer hydrates, are stable in this environment of high depth, temperature, and pressure. The most nearly stable configuration corresponds to the one-layer hydrate, characterized by a d001 spacing of 12.75 A, the adsorbed water being 60 molecules/layer or 180.83 mg of H2O/g of clay, an internal energy of -22.73 kcal/mol, an interlayer density of 0.365 g/mL, and a pressure tensor, Pzz, of 1999.9 bar. The interlayer structure consists of two close layers of water molecules 0.50 A from the midplane, with broad shoulders on the sides, the protons oriented toward the midplane and the siloxane surfaces, and the K+ ions close to the clay surfaces and on the interlayer midplane.     

Liquid-liquid phase transitions in supercooled water studied by computer simulations of various water models

Liquid-liquid and liquid-vapor coexistence regions of various water models were determined by Monte Carlo (MC) simulations of isotherms of density fluctuation-restricted systems and by Gibbs ensemble MC simulations. All studied water models show multiple liquid-liquid phase transitions in the supercooled region: we observe two transitions of the TIP4P, TIP5P, and SPCE models and three transitions of the ST2 model. The location of these phase transitions with respect to the liquid-vapor coexistence curve and the glass temperature is highly sensitive to the water model and its implementation. We suggest that the apparent thermodynamic singularity of real liquid water in the supercooled region at about 228 K is caused by an approach to the spinodal of the first (lowest density) liquid-liquid phase transition. The well-known density maximum of liquid water at 277 K is related to the second liquid-liquid phase transition, which is located at positive pressures with a critical point close to the maximum. A possible order parameter and the universality class of liquid-liquid phase transitions in one-component fluids are discussed.     

受限狭缝中CO_2-CH_4体系的汽液相平衡性质

受限条件下CO2-CH4体系的相平衡性质对化工工艺过程的设计具有非常重要的意义.采用Gibbs系综Monte Carlo模拟,对220K下CO2-CH4体系在主体相和受限狭缝中的相平衡性质进行了系统地研究.通过主体相模拟与实验结果比较,验证了流体分子势能参数的合理性;通过改变狭缝壁面原子的能量参数,研究了受限环境对CO2-CH4体系汽液相平衡性质的影响.与主体相相比,在硬壁狭缝中,CO2-CH4体系的露点压力增加,泡点压力降低,压力-组成相图变窄,且体系更容易达到超临界状态;在吸引狭缝中,随壁面原子能量参数的增大,CO2-CH4体系的压力-组成相图上移,临界点处CH4的摩尔分数减小,相图变窄.在体系汽液相总组成相同情况下,硬壁狭缝内体系的汽液相密度均比主体相中小;随壁面原子能量参数增大,气相密度变大、液相密度在CH4的摩尔分数较小时变大而当CH4的摩尔分数达到一定值后反而减小.在体系汽液相总组成相同时,受限环境下的汽化热比主体相的汽化热小且随壁面吸引势的增强越来越小;在主体相和硬壁狭缝中体系的汽化热随CH4含量的增加单调减小,而当壁面势能参数较大时汽化热随CH4含量增加先增大后减小.     

Monte Carlo simulation of polarizable systems: Early rejection scheme for improving the performance of adiabatic nuclear and electronic sampling Monte Carlo simulations

An early rejection scheme for trial moves in adiabatic nuclear and electronic sampling Monte Carlo simulation (ANES-MC) of polarizable intermolecular potential models is presented. The proposed algorithm is based on Swendsen–Wang filter functions for prediction of success or failure of trial moves in Monte Carlo simulations. The goal was to reduce the amount of calculations involved in ANES-MC electronic moves, by foreseeing the success of an attempt before making those moves. The new method was employed in Gibbs ensemble Monte Carlo (GEMC) simulations of the polarizable simple point charge-fluctuating charge (SPC-FQ) model of water. The overall improvement in GEMC depends on the number of swap attempts (transfer molecules between phases) in one Monte Carlo cycle. The proposed method allows this number to increase, enhancing the chemical potential equalization. For a system with 300 SPC-FQ water molecules, for example, the fractions of early rejected transfers were about 0.9998 and 0.9994 at 373 and 423 K, respectively. This means that the transfer moves consume only a very small part of the overall computing effort, making GEMC almost equivalent to a simulation in the canonical ensemble.     

Prediction of viscosities and vapor–liquid equilibria for five polyhydric alcohols by molecular simulation

Reverse nonequilibrium molecular dynamics in the canonical ensemble and coupled–decoupled configurational-bias Monte Carlo simulations in the Gibbs ensemble were used to predict the low-shear rate Newtonian viscosities and vapor–liquid coexistence curves for 1,2-butanediol, 1,3-butanediol, 1,4-butanediol, 2-methyl-1,3-propanediol, and 1,2,4-butanetriol modeled with the transferable potentials for phase equilibria-united atom (TraPPE-UA) force field. Comparison with available experimental data demonstrates that the TraPPE-UA force field yields very good predictions of the viscosities and vapor–liquid coexistence curves. A detailed analysis of liquid structure and hydrogen bonding is provided.     

Solubility in supercritical carbon dioxide: importance of the Poynting correction and entrainer effects

Gibbs ensemble Monte Carlo simulations were used to investigate the effect of pressure and of entrainers on the solubility of low-volatility species in CO 2. Two entrainers were examined, n-octane and methanol, as well as two solutes, hexamethylbenzene and benzoic acid. For the three pressures studied (12, 20, and 28 MPa), the simulations demonstrate that the increase in the solubility with increasing pressure is mostly due to an increase in the solute's chemical potential (as expressed by the Poynting correction) and not due to an increase in the solvent strength of supercritical CO 2. The presence of an entrainer enhances solubility, particularly when the solute and entrainer can form hydrogen bonds. The solubility of benzoic acid is enhanced by an order of magnitude upon addition of methanol entrainer, whereas the enhancements are less than 2 for the other systems.     

A simulation method for the calculation of chemical potentials in small, inhomogeneous, and dense systems

We present a modification of the gauge cell Monte Carlo simulation method [A. V. Neimark and A. Vishnyakov, Phys. Rev. E 62, 4611 (2000)] designed for chemical potential calculations in small confined inhomogeneous systems. To measure the chemical potential, the system under study is set in chemical equilibrium with the gauge cell, which represents a finite volume reservoir of ideal particles. The system and the gauge cell are immersed into the thermal bath of a given temperature. The size of the gauge cell controls the level of density fluctuations in the system. The chemical potential is rigorously calculated from the equilibrium distribution of particles between the system cell and the gauge cell and does not depend on the gauge cell size. This scheme, which we call a mesoscopic canonical ensemble, bridges the gap between the canonical and the grand canonical ensembles, which are known to be inconsistent for small systems. The ideal gas gauge cell method is illustrated with Monte Carlo simulations of Lennard-Jones fluid confined to spherical pores of different sizes. Special attention is paid to the case of extreme confinement of several molecular diameters in cross section where the inconsistency between the canonical ensemble and the grand canonical ensemble is most pronounced. For sufficiently large systems, the chemical potential can be reliably determined from the mean density in the gauge cell as it was implied in the original gauge cell method. The method is applied to study the transition from supercritical adsorption to subcritical capillary condensation, which is observed in nanoporous materials as the pore size increases.     

Monte Carlo simulation of the self-assembly and phase behavior of semiflexible equilibrium polymers

Grand canonical Monte Carlo simulations of a simple model semiflexible equilibrium polymer system, consisting of hard sphere monomers reversibly self-assembling into chains of arbitrary length, have been performed using a novel sampling method to add or remove multiple monomers during a single MC move. Systems with two different persistence lengths and a range of bond association constants have been studied. We find first-order lyotropic phase transitions between isotropic and nematic phases near the concentrations predicted by a statistical thermodynamic theory, but with significantly narrower coexistence regions. A possible contribution to the discrepancy between theory and simulation is that the length distribution of chains in the nematic phase is bi-exponential, differing from the simple exponential distribution found in the isotropic phase and predicted from a mean-field treatment of the nematic. The additional short length-scale characterizing the distribution appears to arise from the lower orientational order of short chains. The dependence of this length-scale on chemical potential, bond association constant, and total monomer concentration has been examined.     

Ordering in surfactant mixtures induced by polymers

We studied ternary mixtures of nonionic surfactant (C12E6, n-dodecyl hexaoxyethylene glycol monoether), polymer (PEG, polyethylene glycol), and water. A small amount of PEG induces demixing into the polymer-rich and surfactant-rich phases in the ternary PEG/C12E6/water mixture. Above a certain concentration and/or molecular weight of PEG, the surfactant-rich phase orders, even in a solution consisting of a few percent of surfactant. The phase boundary acts as a semipermeable membrane, and the equilibrium is determined by the chemical potential of water in two phases. The explicit expression for the amount of PEG needed to order C12E6 water solution is given and verified experimentally. The analysis of the coexistence conditions leads to the conjecture that only two oxygen atoms in the outward part of the hydrophilic surfactant head strongly affect the chemical potential of water. Our methodology is generic, i.e., on the same basis one can design a similar experiment for any surfactant/polymer/water system and find the right proportions of polymer that induce order in a surfactant-rich phase.     

Surface tension of the Widom-Rowlinson model

We consider the computation of the surface tension of the fluid-fluid interface for the Widom-Rowlinson [J. Chem. Phys. 52, 1670 (1970)] binary mixture from direct simulation of the inhomogeneous system. We make use of the standard mechanical route, in which the surface tension follows from the computation of the normal and tangential components of the pressure tensor of the system. In addition to the usual approach, which involves simulations of the inhomogeneous system in the canonical ensemble, we also consider the computation of the surface tension in an ensemble where the pressure perpendicular (normal) to the planar interface is kept fixed. Both approaches are seen to provide consistent values of the interfacial tension. The issue of the system-size dependence of the surface tension is addressed. In addition, simulations of the fluid-fluid coexistence properties of the mixture are performed in the semigrand canonical ensemble. Our results are compared with existing data of the Widom-Rowlinson mixture and are also examined in the light of the vapor-liquid equilibrium of the thermodynamically equivalent one-component penetrable sphere model.