A molecular simulation strategy allowing the computation of sorption isotherms, up to high penetrant activities, in glassy polymer matrices will be presented. The new scheme permits an equilibrium repartition of the penetrant molecules in the polymer matrix. Glassy matrices of polystyrene (PS) were obtained from reverse mapping of equilibrated coarse-grained configurations (Spyriouni et al., Macromolecules 2007, 40, 3876). The sorption isotherms of CO₂ in atactic PS and the induced polymer swelling have been calculated in the temperature range from 308 K to 405 K (below and above the glass transition temperature of the pure polymer) and for pressures up to 300 bar. The results compare favorably with available experimental data. Derivative properties, such as the partial molar volume and the partial molar enthalpy of CO₂, have been calculated in good agreement with experimental values. The segmental dynamics of PS have been analyzed and the pressure induced glass transition has been calculated for various temperatures below the T_g of pure PS.

The glassy matrices are loaded with CO₂ at various concentrations by insertions of CO₂ molecules, after mapping the accessible free volume. The loaded matrices are subjected to molecular dynamics (MD) simulations at constant temperature, pressure and number of molecules (N₁N₂PT ensemble). The penetrant is redistributed among the accessible cavities in order to achieve an equilibrium repartitioning in the thermally fluctuating polymer matrix. This is accomplished by selecting periodically one configuration and attempting a certain number of CO₂ displacements between the accessible cavities. The latter are determined by Delaunay tessellation followed by volume and connectivity analysis, as implemented by Greenfield and Theodorou (Macromolecules 1993, 26, 5461). Following acceptance or rejection of the attempted repartitioning moves, we proceed with the MD simulation of the polymer-penetrant system.

The excess chemical potential of CO₂ is calculated using the Direct Particle Deletion (DPD) method, a generalization of Staged Particle Deletion (Boulougouris et al. Mol. Phys. 1999, 96, 905; J. Chem. Phys., 2001, 115, 8231). This method is very powerful for dense systems and/or large solute molecules in comparison to traditional Widom insertion. Moreover, the DPD version of the method permits the calculation of the excess chemical potential from a single simulation run.

The fugacity of CO₂ in the polymer matrix is calculated from the excess chemical potential via the relationship: f_CO2 = N₂k_BT exp(bm^ex)/<V>, where N₂ is the number of CO₂ molecules in the polymer matrix and <V> is the average volume under NPT conditions. A new pressure is calculated via an equation of state for pure gaseous CO₂, such as SAFT, by setting the gas-phase fugacity equal to f_CO2. The calculations are repeated at the new pressure, until the CO₂ fugacity in the polymer matrix and in the gas are equal. The scheme converges rapidly, since the CO₂ fugacity in the polymer is not very sensitive to the pressure.