The glassy matrices are loaded with CO2 at various concentrations by insertions of CO2 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 (N1N2PT 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 CO2 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 CO2 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 CO2 in the polymer matrix is calculated from the excess chemical potential via the relationship: fCO2 = N2kBT exp(bmex)/<V>, where N2 is the number of CO2 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 CO2, such as SAFT, by setting the gas-phase fugacity equal to fCO2. The calculations are repeated at the new pressure, until the CO2 fugacity in the polymer matrix and in the gas are equal. The scheme converges rapidly, since the CO2 fugacity in the polymer is not very sensitive to the pressure.