Grand Canonical Ensemble Simulations Research Articles

Understanding the origins of reversible and hysteretic pathways of adsorption phase transitions in metal-organic frameworks

Phase behavior of nanoconfined fluids adsorbed in metal–organic frameworks is of paramount importance for the design of advanced materials for energy and gas storage, separations, electrochemical devices, sensors, and drug delivery, as well as for the pore structure characterization. Phase transformations in adsorbed fluids often involve long-lasting metastable states and hysteresis that has been well-documented in gas adsorption–desorption and nonwetting fluid intrusion-extrusion experiments. However, theoretical prediction of the observed nanophase behavior remains a challenging problem. The mesoscopic canonical, or mesocanonical, ensemble (MCE) is devised to study the nanophase behavior under conditions of controlled fluctuations to stabilize metastable and labile states. Here, we implement and apply the MCE Monte Carlo (MCEMC) simulation scheme to predict the origins of reversible and hysteric adsorption phase transitions in a series of practical MOF materials, including IRMOF-1, ZIF-412, UiO-66, Cu-BTC, IRMOF-74-V, VII, and IX. The MCEMC method, called the gauge cell method, allows to produce Van der Waals type isotherms with distinctive swings around the phase transition regions. The constructed isotherms determine the positions of phase equilibrium and spinodals, as well as the nucleation barriers separating metastable states. We demonstrate the unique capabilities of the MCEMC method in quantitative predictions of experimental observations compared with the conventional grand canonical and canonical ensemble simulations. The MCEMC method is implemented in the open-source RASPA and LAMMPS software packages and recommended for studies of adsorption behavior and pore structure characterization of MOFs and other nanoporous materials.

Constant inner potential DFT for modelling electrochemical systems under constant potential and bias

Electrochemical systems play a decisive role in, e.g. clean energy conversion but understanding their complex chemistry remains an outstanding challenge. Constant potential and grand canonical ensemble (GCE) simulations are indispensable for unraveling the properties of electrochemical processes as a function of the electrode potential. Currently, GCE calculations performed at the density functional theory (DFT) level require fixing the Fermi level within the simulation cell. Here, we illustrate that this method is inadequate when modeling outer sphere reactions and a biased two-electrode cell. For these systems, the Fermi level obtained from DFT calculations does not accurately present the experimentally controlled electrode potential or describe the thermodynamic independent variable in GCE-DFT. To address this limitation, we developed and implemented a constant inner potential (CIP) method offering a more robust and general approach to conducting GCE-DFT simulations of electrochemical systems under constant potential or bias conditions. The primary advantage of CIP is that it uses the local electrode inner potential as the thermodynamic parameter for the electrode potential, as opposed to the global Fermi level. Through numerical and analytical studies, we demonstrate that the CIP and Fermi level GCE-DFT approaches are equivalent for metallic electrodes and inner-sphere reactions. However, CIP proves to be more versatile, as it can be applied to outer-sphere and two-electrode systems, addressing the limitations of the constant Fermi-level approach in these scenarios. Altogether, the CIP approach stands out as a general and efficient GCE-DFT method simulating electrochemical interfaces from first principles.

Water uptake by gecko β-keratin and the influence of relative humidity on its mechanical and volumetric properties.

Grand canonical ensemble molecular dynamics simulations are done to calculate the water content of gecko β-keratin as a function of relative humidity (RH). For comparison, we experimentally measured the water uptake of scales of the skin of cobra Naja nigricollis. The calculated sigmoidal sorption isotherm is in good agreement with experiment. To examine the softening effect of water on gecko keratin, we have calculated the mechanical properties of dry and wet keratin samples, and we have established relations between the mechanical properties and the RH. We found that a higher RH causes a decrease in the Young's modulus, the yield stress, the yield strain, the stress at failure and an increase in the strain at failure of the gecko keratin. At low RHs (less than 80%), the change in the mechanical properties is small, with most of the changes occurring at higher RHs. The changes in the macroscopic properties of the keratin are explained by the action of sorbed water on the molecular scale. It causes keratin to swell, thereby increasing the distances between amino acids. This has a weakening effect on amino acid interactions and softens the keratin material. The effect is more pronounced at higher RHs.

Thermodynamics of the uniform electron gas: Fermionic path integral Monte Carlo simulations in the restricted grand canonical ensemble

The uniform electron gas (UEG) is one of the key models for the understanding of warm dense matter—an exotic, highly compressed state of matter between solid and plasma phases. The difficulty in modelling the UEG arises from the need to simultaneously account for Coulomb correlations, quantum effects, and exchange effects, as well as finite temperature. The most accurate results so far were obtained from quantum Monte Carlo (QMC) simulations with a variety of representations. However, QMC for electrons is hampered by the fermion sign problem. Here, we present results from a novel fermionic propagator path integral Monte Carlo in the restricted grand canonical ensemble. The ab initio simulation results for the spin‐resolved pair distribution functions and static structure factor are reported for two isotherms (T in the units of the Fermi temperature). Furthermore, we combine the results from the linear response theory in the Singwi‐Tosi‐Land‐Sjölander scheme with the QMC data to remove finite‐size errors in the interaction energy. We present a new corrected parametrization for the interaction energy and the exchange–correlation free energy in the thermodynamic limit, and benchmark our results against the restricted path integral Monte Carlo by Brown et al. [Phys. Rev. Lett. 110, 146405 (2013)] and configuration path integral Monte Carlo/permutation‐blocking path integral Monte Carlo by Dornheim et al. [Phys. Rev. Lett. 117, 115701 (2016)].

Rimantadine Binds to and Inhibits the Influenza A M2 Proton Channel without Enantiomeric Specificity.

The influenza A M2 wild-type (WT) proton channel is the target of the anti-influenza drug rimantadine. Rimantadine has two enantiomers, though most investigations into drug binding and inhibition have used a racemic mixture. Solid-state NMR experiments using the full length-M2 WT have shown significant spectral differences that were interpreted to indicate tighter binding for (R)- vs (S)-rimantadine. However, it was unclear if this correlates with a functional difference in drug binding and inhibition. Using X-ray crystallography, we have determined that both (R)- and (S)-rimantadine bind to the M2 WT pore with slight differences in the hydration of each enantiomer. However, this does not result in a difference in potency or binding kinetics, as shown by similar values for kon, koff, and Kd in electrophysiological assays and for EC50 values in cellular assays. We concluded that the slight differences in hydration for the (R)- and (S)-rimantadine enantiomers are not relevant to drug binding or channel inhibition. To further explore the effect of the hydration of the M2 pore on binding affinity, the water structure was evaluated by grand canonical ensemble molecular dynamics simulations as a function of the chemical potential of the water. Initially, the two layers of ordered water molecules between the bound drug and the channel's gating His37 residues mask the drug's chirality. As the chemical potential becomes more unfavorable, the drug translocates down to the lower water layer, and the interaction becomes more sensitive to chirality. These studies suggest the feasibility of displacing the upper water layer and specifically recognizing the lower water layers in novel drugs.

Influence of electrical images and electrolyte concentration on capacitance of the electrode - molten salt interface

Results of Monte Carlo grand canonical ensemble simulations of electrical double layer (EDL), with and without electrical images on the electrode, are presented. Calculations were carried out for a wide range of electrode surface charge, σ from -0.7 to 0.7 Cm−2 and a wide range electrolyte concentrations, c = 1.00 – 8.67 M. The electrolyte concentration of 8.67 M corresponds to the density of molten salts. It is known from earlier studies that for low electrolyte concentrations and low electrode charges, the EDL capacitances obtained for electrical images included in the calculations are greater than those calculated excluding them. For high electrolyte concentrations, especially for those which correspond to the molten salt densities, the situation unexpectedly turns around, i.e. the EDL capacitances obtained for the case of excluded images are greater than those calculated including them. For the whole range of applied electrolyte concentrations and electrode charges greater than σ ≈ 0.4 Cm−2 and lower than σ ≈ -0.4 Cm−2 no visible influence of applying the electrical images is observed on the EDL capacitance values. The above observations are discussed, explained in terms of the relationship between EDL thickness and the capacitance value, and analysed in terms of the electrostatic interactions of ion charges with their images.

How Alcoholic Disinfectants Affect Coronavirus Model Membranes: Membrane Fluidity, Permeability, and Disintegration.

Atomistic molecular dynamics simulations have been carried out with a view to investigating the stability of the SARS-CoV-2 exterior membrane with respect to two common disinfectants, namely, aqueous solutions of ethanol and n-propanol. We used dipalmitoylphosphatidylcholine (DPPC) as a model membrane material and did simulations on both gel and liquid crystalline phases of membrane surrounded by aqueous solutions of varying alcohol concentrations (up to 17.5 mol %). While a moderate effect of alcohol on the gel phase of membrane is observed, its liquid crystalline phase is shown to be influenced dramatically by either alcohol. Our results show that aqueous solutions of only 5 and 10 mol % alcohol already have significant weakening effects on the membrane. The effects of n-propanol are always stronger than those of ethanol. The membrane changes its structure, when exposed to disinfectant solutions; uptake of alcohol causes it to swell laterally but to shrink vertically. At the same time, the orientational order of lipid tails decreases significantly. Metadynamics and grand-canonical ensemble simulations were done to calculate the free-energy profiles for permeation of alcohol and alcohol/water solubility in the DPPC. We found that the free-energy barrier to permeation of the DPPC liquid crystalline phase by all permeants is significantly lowered by alcohol uptake. At a disinfectant concentration of 10 mol %, it becomes insignificant enough to allow almost free passage of the disinfectant to the inside of the virus to cause damage there. It should be noted that the disinfectant also causes the barrier for water permeation to drop. Furthermore, the shrinking of the membrane thickness shortens the gap needed to be crossed by penetrants from outside the virus into its core. The lateral swelling also increases the average distance between head groups, which is a secondary barrier to membrane penetration, and hence further increases the penetration by disinfectants. At alcohol concentrations in the disinfectant solution above 15 mol %, we reliably observe disintegration of the DPPC membrane in its liquid crystalline phase.

H2S separation from biogas by adsorption on functionalized MIL-47-X (X = −OH and − OCH3): A simulation study

Monte Carlo molecular simulations in grand canonical ensemble were performed to study the H2S, CO2 and CH4 adsorption and separation on the MIL-47 and the functionalized MIL-47-X (X = -OH and -OCH3). The preferred situation for the adsorption of gas molecules on the studied MOFs was also investigated. The calculated adsorption isotherms of CO2 on the MIL-47 agreed well with the experimental uptakes. Furthermore, the simulated adsorption selectivities of three MOFs for the CO2/CH4, H2S/CH4 and H2S/CO2 mixtures were in the order H2S/CH4 > H2S/CO2 > CO2/CH4. These MOFs had higher selectivity for H2S over CH4 indicating that they can be good adsorbents for the H2S separation. The probability distribution of the adsorbed molecules showed that the guest molecules were more localized near the metal centers and functional groups.

Heat capacities of supercritical fluids via Grand Canonical ensemble simulations

In this work, suitable mathematical relationships to compute isobaric heat capacities from molecular simulations in the Grand Canonical (GC) ensemble are derived and tested via Monte Carlo methods. Using atomistic classical force fields, the residual isobaric heat capacities of pure carbon dioxide (CO2) and pure methanol (MeOH) were obtained at supercritical conditions (with critical properties estimated from a finite-size scaling analysis). The total isobaric heat capacity was determined by combining the residual isobaric heat capacity obtained from molecular simulations with the ideal gas contributions obtained from experimental correlations. Isobaric heat capacities generated from both GC and Isothermal–Isobaric ensemble simulations were compared to predictions from accurate equations of state (EOS)s for CO2 and MeOH at corresponding reduced temperatures and pressures. Isobaric heat capacities calculated from both ensembles were in good agreement with those obtained from the Span and Wagner EOS for CO2 and the IUPAC EOS for MeOH. For comparable computation times, simulations run in the GC ensemble generate results with significantly lower statistical uncertainty than those run in the Isothermal–Isobaric ensemble.

On the consistency of NVT, NPT, μVT and Gibbs ensembles in the framework of kinetic Monte Carlo – Fluid phase equilibria and adsorption of pure component systems

This paper aims to show the consistency between simulations of fluid phase properties, obtained with various ensembles, developed within the framework of kinetic Monte Carlo (kMC) simulation: NVT (canonical), NPT (isothermal-isobaric systems), μVT (grand canonical) and Gibbs ensembles, to ensure the reliability of the kMC methodology. The advantages of the kMC scheme, as compared to the conventional Metropolis Monte Carlo, are: (1) accurate determination of the chemical potentials compared to the Widom insertion method, and (2) a rejection-free algorithm, making the implementation of the kMC scheme simpler. For internal consistency in a grand canonical ensemble simulation, we have developed a means to calculate the intrinsic chemical potential of the system accurately, which must be the same (within statistical error of the simulation) as the specified chemical potential to ensure convergence to equilibrium. We test the consistency of canonical (NVT-kMC) and grand canonical (GC-kMC) ensembles for argon adsorption at 87K and 120K in a uniform open-ended slit pore, and hence derive governing factors affecting hysteresis in the isotherm and the microscopic mechanisms of condensation and evaporation.

An efficientmethod to determine chemical potential of mixtures in the isothermal and isobaric bulk phase with kineticMonte Carlo simulation

ABSTRACTA kinetic Monte Carlo (kMC) scheme in the NPT ensemble (constant number of molecules, pressure and temperature) has been developed to determine accurate chemical potentials for all components in a homogeneous mixture. The simulation requires two moves: (1) a displacement move and (2) a volume change move. In the former, the mobility rate of a selected molecule is determined by its interaction with all the other molecules in the system and is moved to a random position within the simulation box, according to the Rosenbluth algorithm, without any rejections (entropic sampling). The volume change move is decided by a comparison between either the instant pressure or the partial average pressure (with long-range correction) and the specified pressure and is carried out much less frequently than the displacement move. We applied this NPT scheme to a number of mixtures in both the gaseous and liquid phases, and show that the derived chemical potentials are accurate and reproducible. The method is recommended for obtaining chemical potentials in mixtures that are required as input in a grand canonical ensemble simulation.

Grand canonical Monte Carlo using solvent repacking: Application to phase behavior of hard disk mixtures.

A new "solvent repacking Monte Carlo" strategy for performing grand canonical ensemble simulations in condensed phases is introduced and applied to the study of hard-disk systems. The strategy is based on the configuration-bias approach, but uses an auxiliary biasing potential to improve the efficiency of packing multiple solvent particles in the cavity formed by removing one large solute. The method has been applied to study the coexistence of ordered and isotropic phases in three binary mixtures of hard disks with a small mole fraction (xL < 0.02) of the larger "solute" component. A chemical potential of 12.81 ± 0.01 kBT was found to correspond to the freezing transition of the pure hard disk "solvent." Simulations permitted the study of partitioning of large disks between ordered and isotropic phases, which showed a distinct non-monotonic dependence on size; the isotropic phase was enriched approximately 10-fold, 20-fold, and 5-fold over the coexisting ordered phases at diameter ratios d = 1.4, 2.5, and 3, respectively. Mixing of large and small disks within both phases near coexistence was strongly non-ideal in spite of the dilution. Structures of systems near coexistence were analyzed to determine correlations between large disks' positions within each phase, the orientational correlation length of small disks within the fluid phases, and the nature of translational order in the ordered phase. The analyses indicate that the ordered phase coexists with an isotropic phase resembling a nanoemulsion of ordered domains of small disks, with large disks enriched at the disordered domain interfaces.

Accurate and precise determination of critical properties from Gibbs ensemble Monte Carlo simulations.

Since the seminal paper by Panagiotopoulos [Mol. Phys. 61, 813 (1997)], the Gibbs ensemble Monte Carlo (GEMC) method has been the most popular particle-based simulation approach for the computation of vapor-liquid phase equilibria. However, the validity of GEMC simulations in the near-critical region has been questioned because rigorous finite-size scaling approaches cannot be applied to simulations with fluctuating volume. Valleau [Mol. Simul. 29, 627 (2003)] has argued that GEMC simulations would lead to a spurious overestimation of the critical temperature. More recently, Patel et al. [J. Chem. Phys. 134, 024101 (2011)] opined that the use of analytical tail corrections would be problematic in the near-critical region. To address these issues, we perform extensive GEMC simulations for Lennard-Jones particles in the near-critical region varying the system size, the overall system density, and the cutoff distance. For a system with N = 5500 particles, potential truncation at 8σ and analytical tail corrections, an extrapolation of GEMC simulation data at temperatures in the range from 1.27 to 1.305 yields T(c) = 1.3128 ± 0.0016, ρ(c) = 0.316 ± 0.004, and p(c) = 0.1274 ± 0.0013 in excellent agreement with the thermodynamic limit determined by Potoff and Panagiotopoulos [J. Chem. Phys. 109, 10914 (1998)] using grand canonical Monte Carlo simulations and finite-size scaling. Critical properties estimated using GEMC simulations with different overall system densities (0.296 ≤ ρ(t) ≤ 0.336) agree to within the statistical uncertainties. For simulations with tail corrections, data obtained using r(cut) = 3.5σ yield T(c) and p(c) that are higher by 0.2% and 1.4% than simulations with r(cut) = 5 and 8σ but still with overlapping 95% confidence intervals. In contrast, GEMC simulations with a truncated and shifted potential show that r(cut) = 8σ is insufficient to obtain accurate results. Additional GEMC simulations for hard-core square-well particles with various ranges of the attractive well and for n-decane molecules represented by the TraPPE force field yield data that support the trends observed for Lennard-Jones particles. The finite-size dependence of the critical properties obtained from GEMC simulations is significantly smaller than those from grand-canonical ensemble simulations. Thus, when resources are not available for a rigorous finite-size scaling study, GEMC simulations provide a straightforward route to determine fairly accurate critical properties using relatively small system sizes.

From determination of the fugacity coefficients to estimation of hydrogen storage capacity: A convenient theoretical method

The equation of state (EOS) from virial expansion (VE) is used in this work to pave the way for determining the fugacity coefficients of the hydrogen fluid at arbitrary temperature and pressure. The fugacity coefficients from our VE method have more physical meanings than the empirical values. In this way, the hydrogen storage capacity of a novel material model can be estimated by using few density functional theory (DFT) calculations with the aid of a continuum model. The efficient continuum model can provide a more accurate estimation of the hydrogen storage capacity than the pure DFT calculations. Furthermore, the expensive grand canonical ensemble (μNT) simulations combining with the quantum mechanics methods (i.e., QM/MD-μNT) are unnecessary within this method. The hydrogen fluid can be handled with our VE method at the temperature in the range of 160–773 K. The hydrogen storage capacity and the detailed thermodynamic information of a designed novel material can thereby be estimated by using this method with relatively high accuracy and low computing cost. As an example, the hydrogen storage capacities of the expanded bilayer graphene systems are presented. Our theoretical results agree with the experimental values very well.

Open boundary molecular dynamics

This contribution analyzes several strategies and combination of methodologies to perform molecular dynamic simulations in open systems. Here, the term open indicates that the total system has boundaries where transfer of mass, momentum and energy can take place. This formalism, which we call Open Boundary Molecular Dynamics (OBMD), can act as interface of different schemes, such as Adaptive Resolution Scheme (AdResS) and Hybrid continuum-particle dynamics to link atomistic, coarse-grained (CG) and continuum (Eulerian) fluid dynamics in the general framework of fluctuating Navier-Stokes equations. The core domain of the simulation box is solved using all-atom descriptions. The CG layer introduced using AdResS is located at the outer part of the open box to make feasible the insertion of large molecules into the system. Communications between the molecular system and the outer world are carried out in the outer layers, called buffers. These coupling preserve momentum and mass conservation laws and can thus be linked with Eulerian hydro- dynamic solvers. In its simpler form, OBMD allows, however, to impose a local pressure tensor and a heat flux across the system’s boundaries. For a one component molecular system, the external normal pressure and temperature determine the external chemical potential and thus the independent parameters of a grand-canonical ensemble simulation. Extended ensembles under non-equilibrium stationary states can also be simulated as well as time dependent forcings (e.g. oscillatory rheology). To illustrate the robustness of the combined OBMD-AdResS method, we present simulations of star-polymer melts at equilibrium and in sheared flow.

Adsorption of binary mixtures on two-dimensional surfaces: theory and Monte Carlo simulations

The adsorption of binary mixtures on square lattices is studied by combining theoretical modeling and Monte Carlo (MC) simulations in grand canonical ensemble. The adsorption thermodynamics is analyzed through the total and partial isotherms. Two theoretical models have been used in the present study: (i) the first, which we called cluster approximation (CA), is based on exact calculations of configurations on finite cells. An efficient algorithm allows us to calculate the detailed structure of the configuration space for $$m = l \times l $$ cells; and (ii) the second is a generalization of the classical quasi-chemical approximation (QCA) in which the adsorbate is a binary mixture of species $$a$$ and $$b$$ . Adsorbate–adsorbate lateral interactions are incorporated in the context of the two mentioned approximations. Results from CA and QCA are compared with MC simulations. Close agreement between simulated and theoretical data supports the validity of the theoretical models to describe the adsorption of mixed gases on two-dimensional surfaces.

Sampling of Organic Solutes in Aqueous and Heterogeneous Environments Using Oscillating Excess Chemical Potentials in Grand Canonical-like Monte Carlo-Molecular Dynamics Simulations.

Solute sampling of explicit bulk-phase aqueous environments in grand canonical (GC) ensemble simulations suffer from poor convergence due to low insertion probabilities of the solutes. To address this, we developed an iterative procedure involving Grand Canonical-like Monte Carlo (GCMC) and molecular dynamics (MD) simulations. Each iteration involves GCMC of both the solutes and water followed by MD, with the excess chemical potential (μex) of both the solute and the water oscillated to attain their target concentrations in the simulation system. By periodically varying the μex of the water and solutes over the GCMC-MD iterations, solute exchange probabilities and the spatial distributions of the solutes improved. The utility of the oscillating-μex GCMC-MD method is indicated by its ability to approximate the hydration free energy (HFE) of the individual solutes in aqueous solution as well as in dilute aqueous mixtures of multiple solutes. For seven organic solutes: benzene, propane, acetaldehyde, methanol, formamide, acetate, and methylammonium, the average μex of the solutes and the water converged close to their respective HFEs in both 1 M standard state and dilute aqueous mixture systems. The oscillating-μex GCMC methodology is also able to drive solute sampling in proteins in aqueous environments as shown using the occluded binding pocket of the T4 lysozyme L99A mutant as a model system. The approach was shown to satisfactorily reproduce the free energy of binding of benzene as well as sample the functional group requirements of the occluded pocket consistent with the crystal structures of known ligands bound to the L99A mutant as well as their relative binding affinities.

Sorption and diffusion of carbon dioxide and nitrogen in poly(methyl methacrylate)

Molecular dynamics simulations are performed to determine the solubility and diffusion coefficient of carbon dioxide and nitrogen in poly(methyl methacrylate) (PMMA). The solubilities of CO2 in the polymer are calculated employing our grand canonical ensemble simulation method, fixing the target excess chemical potential of CO2 in the polymer and varying the number of CO2 molecules in the polymer matrix till establishing equilibrium. It is shown that the calculated sorption isotherms of CO2 in PMMA, employing this method well agrees with experiment. Our results on the diffusion coefficients of CO2 and N2 in PMMA are shown to obey a common hopping mechanism. It is shown that the higher solubility of CO2 than that of N2 is a consequence of more attractive interactions between the carbonyl group of polymer and the sorbent. While the residence time of CO2 beside the carbonyl group of polymer is about three times higher than that of N2, the diffusion coefficient of CO2 in PMMA is higher than that of N2. The higher diffusion coefficient of CO2, compared to N2, in PMMA is shown to be due to the higher (≈3 times) swelling of polymer upon CO2 uptake.

Molecular Simulations of Water Adsorbed on Mesoporous Silica Thin Films

Grand canonical Monte Carlo and canonical ensemble molecular dynamics simulations were used to investigate water adsorption properties in mesoporous silica thin films. The effect of pore radius on the adsorption properties was assessed using two models of mesoporous silica thin films having different pore radius and film thickness (1.38 and 5.66 nm in Model 1, respectively, and 1.81 and 7.30 nm in Model 2, respectively). In the simulations, a water adsorption layer or water menisci were formed in a mesopore accompanying the growth or shrinkage of stable adsorption layers on the upper and lower surfaces. The thickness of a water adsorption layer immediately before the onset of capillary condensation and the curvature radius of a water meniscus prior to capillary evaporation were identified. In the initial stage of layer adsorption, the surface state and radius of the pore affected the coordination structure of water around silanol groups. In the pore-filling state, the hydrogen bond network of water in the pore was similar but not the same as that in the bulk liquid; the difference was because the pore wall affected the position of water molecules. This difference may confer unique dynamic properties upon the water in mesoporous silica.

Grand canonical ensemble molecular dynamics simulation of water solubility in polyamide-6,6

Grand canonical ensemble molecular dynamics simulation is employed to calculate the solubility of water in polyamide-6,6. It is shown that performing two separate simulations, one in the polymeric phase and one in the gaseous phase, is sufficient to find the phase coexistence point. In this method, the chemical potential of water in the polymer phase is expanded as a first-order Taylor series in terms of pressure. Knowing the chemical potential of water in the polymer phase in terms of pressure, another simulation for water in the gaseous phase, in the grand canonical ensemble, is done in which the target chemical potential is set in terms of pressure in the gas phase. The phase coexistence point can easily be calculated from the results of these two independent simulations. Our calculated sorption isotherms and solubility coefficients of water in polyamide-6,6, over a wide range of temperatures and pressures, agree with experimental data.