Molecular Dynamics Simulation of Permeation in Polymers

Sorption and diffusion of small molecules through polymers is a topic of broad range of applications in many industrially important phenomena. There is a rapidly growing demand for polymers of specified permeabilities, such as selective memberanes for separation technologies, barrier membranes for packaging applications, foaming, and plasticization. Polymers with a high degree of permeability and permselectivity have been widely used for gas or liquid separation systems based on membranes (Kesting & Fritzsche, 1993), while those with a low degree of permeability have been used in barrier packaging films as containers (Koros, 1990). Therefore, understanding the underlying mechanism of solubility and/or diffusion of peneterants in macromolecular substances is very useful in obtaining a clear picture of the molecular level mechanism of polymer permeability and in design of new membranes.


Introduction
Sorption and diffusion of small molecules through polymers is a topic of broad range of applications in many industrially important phenomena.There is a rapidly growing demand for polymers of specified permeabilities, such as selective memberanes for separation technologies, barrier membranes for packaging applications, foaming, and plasticization.Polymers with a high degree of permeability and permselectivity have been widely used for gas or liquid separation systems based on membranes (Kesting & Fritzsche, 1993), while those with a low degree of permeability have been used in barrier packaging films as containers (Koros, 1990).Therefore, understanding the underlying mechanism of solubility and/or diffusion of peneterants in macromolecular substances is very useful in obtaining a clear picture of the molecular level mechanism of polymer permeability and in design of new membranes.
The first study of gas permeation through a polymer goes back to 1829 by Thomas Graham (Stannett, 1978).According to Graham's postulate, the permeation process involves the dissolution of penetrant followed by transmission of the dissolved species through the membrane.In the late 1870's, Stefan and Exner demonstrated that gas permeation through a membrane was proportional to the product of solubility and diffusion coefficients.Based on the findings of Stefan and Exner, von Wroblewski (1879) constructed a quantitative solution to the Graham's solution-diffusion model.The dissolution of gas was based on Henry's law of solubility, where the concentration of the gas in the membrane was directly proportional to the applied gas pressure.Later Daynes (1920) showed that it was impossible to evaluate both diffusion and solubility coefficients by steady-state permeation experiments.He presented a mathematical solution, using Fick's second law of diffusion, for calculating the diffusion coefficient, which was assumed to be independent of concentration.
Microscopic theories of the diffusion and permeation of penetrants in polymers, on the other hand, were developed during 1950-1970's.In these theories, to explain the mechanism behind the transport of gas molecules through the free volume present in the polymer matrix, the gas-polymer system was defined in terms of liquid molecules traveling through a liquid membrane.The Cohen and Turnbull (1959) theory considered transport through liquid molecules as hard spheres.In this model the gas diffusion in polymers occurs through the redistribution of free volume.Another theory, proposed by Dibenedetto (1963), used the same concept as that of Cohen and Turnbull (1959) regarding free volume distribution, but with a different chain packing theory at the molecular level.An elaborate theory was developed by Pace and Datyner (1979), assuming that the penetrant molecule is moving along the polymer chain bundles and being stopped only by chain entanglements or a crystallite.The penetrnat molecules were further assumed to jump into the adjacent bundles, similar to the Cohen and Turnbull's model.Such a jump event was considered to be the rate controlling step, with the diffusion along the bundle being three times faster than the perpendicular jump of the molecule.None of these theories, however, provides a molecular lever understanding of the permeability in polymers.
Significant advances in the understanding of gas permeation in polymers have been made only in the last few decades with the advent of powerful computers.Monte Carlo (MC) and molecular dynamics (MD) simulations on small-molecule permeation of amorphous polymers have become feasible in recent decades.The simulations now cover a range of different polymers with varieties of chemical complexities ranging from flexible polymers with simple chemical structure, like polyethylene (Boyd & Pant, 1991), to stiff polymers with detailed chemical structure, like polyamides (Eslami & Mehdipour, 2011).In the following we give a detailed discussion of the calculation methods as well as the polymeric samples employed to study the sorption/diffusion mechanism.The main difficulty in the calculation of gas permeabilities stems from the calculation of solubility coefficients.In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant (sorbate) in the gas phase and that dissolved in polymer.The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation.The chemical potential is, however, coupled to the number of particles and cannot be easily calculated using molecular simulation methods.Therefore, while there exist excellent reviews on the pearmation of gases in polymers, with the main emphasis on the diffusion process (Müller-Plathe, 1994), in this review we explain the polymer permeation putting more emphasis on the solubility process.

Sorption of gases in polymers
Historically, the dual-mode sorption theory (Stannett et al., 1978;Frederickson and Helfand, 1985) describes the concentration C (solubility) of the gas inside a polymer at equilibrium with the gas at a partial pressure p, i.e., 1 where k H is Henry's law solubility coefficient, C  is the saturation concentration of the gas, and b is the affinity coefficient.This model assumes that there are two distinct modes by which a glassy polymer can sorb gas molecules; Henry's law and a Langmuir mechanism which corresponds to the sorption of gases into specific sorption sites in the polymer.Henry's constant has the same physical meaning for glassy polymers as it does for rubbery polymers and liquids, whereas the Langmuir-type term is believed to account for gas sorption into interstitial spaces and microvoids, which are consequences of local heterogeneities and are intimately related to the slow relaxation processes associated with the glassy state of the polymer.Local equilibrium is assumed between the two modes.Figure 1 shows a schematic representation of Henry, Langmuir, and dualmode sorptions.In the low pressure region, Eq. 1 provides the following linear relationship against the pressure: where S 0 is called the apparent solubility coefficient in the zero-pressure limit in glassy polymers.
The determination of the free energy change in sorption of gases in polymer is a problem of primary importance, since the free energy is the main thermodynamic parameter governing the equilibrium between the gas in the pure state and the dissolved one in the polymer.For a system of N particles located at r 1 , r 2 , ...r N , the statistical mechanical expression for the Helmholtz free energy, A, reads as (McQuarrie, 1976): Here, Q is the canonical partition function, U N is the potential energy of the system, r N stands for the whole set of coordinates, r 1 , r 2 , …r N ,  is the de Broglie wavelength, h is the Planck's constant, m is the molecular mass, k B is the Boltzmann's constant, and T is the temperature.Assuming pair-wise additiveity of the potential energy between particles, we have: where u ij is the pair potential interacting between particles i and j and r ij is the interparticle distance.The expression for the chemical potential is obtained by taking the logarithm of the ratios of partition functions for a system composed of N particles, with the sorbate (penetrant) density of  s , and a system composed of N+1 particles, which is obtained by adding one sorbate particle to the previous N-particle system, i.e, Here  s and  s are the chemical potential and the density of sorbate, respectively, and the brackets indicate the ensemble average.Comparing  s with the chemical potential of the ideal gas at the same state, i.e.,

 
3 ln various contributions to the chemical potential can be interpreted according to Ben-Naim (1987).In Eq. ( 8)   to the chemical potential.In fact this term is the contribution of translational motion to the chemical potential.
Looking at the sorption process, we consider the pure sorbate gas s in equilibrium with the dissolved gas in the polymer.Again we connect the thermodynamic states of the system to their corresponding ideal gas states (at the same temperature, and the same density of sorbate s), as it is depicted in Scheme 1.
ideal gas s (T, ) ideal gas s (T,  ) gas s sorbed in polymer (T, ) gas s (T,  ) Scheme 1. Description of equilibrium between the gas s in the gaseous and polymeric phases.
Here,  and  indicate the densities of gas s in the polymeric and gaseous phases, respectively.The condition of equilibrium for gas s, between the gas phase and the polymer phase is the equality of the chemical potentials, i.e.,

   
,, ss TT where  and µ are the chemical potentials of the sorbed gas in polymer and in gas phase, respectively.According to the thermodynamic cycle shown in Scheme 1, the equilibrium condition between the pure gas and the gas sorbed in the polymer is formulated as: Equation 10 is the main equation governing the equilibrium between the sorbed gas molecules in the polymer and those in the gaseous phase.According to this equation, the main parameter to be determined for the calculation of solubilities is the excess chemical potential.In the following section we study the computational approach, with the main emphasis on the MD simulation methods, for calculation of chemical potentials.

Calculation of excess chemical potentials
From the historical point of view, traditional approaches such as the equation-of-state models (Lacombe and Sanchez, 1976 ;Sanchez and Rodgers, 1990 ;von Solms et al., 2005) and the activity coefficient models (Ozkan and Teja, 2005) have been applied to calculate the phase equilibria and sorption of penetrant molecules in polymers.Molecular simulations are the other attractive method for this type of calculation.These methods do not invoke any approximations and predictions are based on well-defined molecular characteristics.In the following sections we describe the application of molecular simulation methods in the case of sorption of gases in polymers.As stated above, the main problem is to calculate the excess chemical potentials, which is given address to in the following sections.

Widom's test particle insertion method
The Widom's test particle insertion method (Widom, 1963) is an elegant method for calculating the excess chemical potentials.In this method a test particle is momentarily inserted at random points in the simulation box and the interaction energy between the test particle and the host particles are averaged to yield the excess chemical potential.
Since the method relies on the statistical accuracy of sampling of configurations that permit the insertion of molecules with low values of the binding energies, its application is more feasible for small-sized solutes.In a dense fluid there is a small probability of finding a cavity to insert a particle in.Therefore, particle insertion becomes a rare event and long simulation times are needed to obtain a good statistics.Similarly, at high densities, the probability of particle deletion is also reduced.The deletion probability depends on the probability of generation of a high density configuration, from which to delete a particle.Due to several limitations, such as poor sampling at high densities or inadequacy of the method in the case of chemical potential calculation for polar systems, the utility of this method to many important applications is precluded.The problem of poor sampling has been solved to some extent by developing such techniques as the umbrella sampling (Shing and K. E. Gubbins, 1981) and map sampling (Deitrick et al., 1989), as the extensions of original Widom's method.The Widom method or its modifications can be used in both MC and MD simulations, but in this method the chemical potential will be calculated at the end of the simulation.This method has widely been applied in calculating the excess chemical potentials of small penetrant molecules in polymers including polypropylene (Müller-Plathe, 1991), polydimethyl siloxane (Sok et al., 1992), polyisobutylene (Müller-Plathe et al., 1993), polystyrene (Eslami & Müller-Plathe, 2007a), and polyimides (Neyertz & Brown, 2004;Pandiyan et al. 2010).However, the method is not applicable in the case of excess chemical potential calculation for polar penetrant in polar polymers.Moreover, the sampling probability gets poor in dense polymeric systems and in the case of big penetrants.

Thermodynamic integration and fast growth thermodynamic integration methods
Thermodynamic integration method is another useful method for the calculation of the excess chemical potentials.In this method, a parameter λ is used to couple the sorbate with the rest of the system and the excess chemical potential, or coupling work, is obtained by integrating along λ as: In this method, several simulations are required to obtain a series of points for the integration.To improve the sampling, recently Hess and van der Vegt ( 2008) have developed a fast-growth thermodynamic integration method.The approach is based on the non-equilibrium work theorem of Jarzynski (1997), which relates the work that is being performed on a system when going from an initial state to a final state along a coordinate λ with the free energy change.Through the applications of this method, polymer microstructure can be considered as a potential landscape as in Figure 2.
In this context, the free energy can be calculated by performing fast growth thermodynamic integration in order to calculate the work required to insert the sorbate in a solvent matrix from sampling the possible coordinates where interaction potentials determine the potential landscape.This method has been applied successfully to calculate the excess chemical potential of relatively large molecules in polymer matrices (Hess & van der Vegt 2008;Hess et al. 2008;Fritz et al., 2009).

Open system simulations
Since the chemical potential is conjugated with the number of particles, the statistical mechanical ensemble representative for its calculation is the grand canonical ensemble.
There are several open-system simulation techniques in the literature, allowing to set the target chemical potential as an independent thermodynamic variable.These methods are preferential to the methods explained above for the calculation of phase equilibria, in that the chemical potential in each phase is set from the beginning as an independent variable, while in closed system simulations the chemical potential is calculated at the end of the run.
In the following two subsections we describe the grand canonical ensemble simulations using MC and MD techniques.

Grand canonical ensemble Monte Carlo simulation
The first grand canonical simulation method (Norman & Filinov, 1969 ) was based on a MC simulation technique.This method was then developed by Adams (1974Adams ( , 1975)).In these methods full particles are inserted into or removed from the simulation box, in a Markov process with a probability, which depends on the target chemical potential.These kinds of movements are compatible with the stochastic nature of moves in MC technique.
The main difficulty with the MC simulation methods in the grand canonical ensemble is their inefficiency when applied to the condensed phases.Besides, the insertion or deletion of full particles highly perturbs the simulation box, and it takes some time for the fluid to relax.These methods envisage similar sampling problems, as addressed above in the case of Widom's method, at high densities.However, the sampling problem has been overcome somewhat by searching in the simulation box for suitable cavities and biasing the insertion event (Escobedo & de Pablo, 1995;Shelley & Patey, 1995;Boda et al., 1997).
Another problem with the MC simulation methods in the grand canonical ensemble is that these methods do not directly reflect the dynamic information about the system.Therefore, attention has been attracted toward development of the grand canonical ensemble simulations using MD technique.

Grand canonical ensemble molecular dynamics simulation
The discontinuity of the number of particles, introduced through particle insertion/deletion, makes it difficult to do MD simulation in the grand canonical ensemble.Therefore, in some of the existing MD-based simulations in the grand canonical ensemble, the sampling is performed through using the MC-based stochastic procedures (Boinepalli & Attard, 2003).Another alternative, is the introduction of the particle number as a continuous dynamic variable through the MD simulation (Eslami & Müller-Plathe, 2007b;Eslami et al., 2011).This approach is based on the so-called extended system method, in which a Hamiltonian including the kinetic and potential energy terms for the extension variable is constructed and the appropriate equations of motion are solved.
H e r e w e a s s u m e t h a t o u r p h y s i c a l s y s t e m i s c o m p o s e d o f N real particles plus one "fractional" particle and is coupled to a heat reservoir and a particle reservoir.The fractional particle is a particle whose potential energy is weighted by a variable  which varies from zero to one.The inclusion of the fractional particle in the system provides a variable number of particles in a dynamical way.The Hamiltonian of such a combined system; the real N particle system, the fractional particle, and the additional degrees of freedom, i.e., the material reservoir and the heat reservoir, is written as (Eslami & Müller-Plathe, 2007b): where H is the Hamiltonian, subscripts i and j refer to the real particles, subscript f refers to the fractional particle, p i is the momentum of the i th particle, s is the velocity scaling variable which couples the system to the thermostat, and  is the particle-number scaling variable which couples the fractional particle to the rest of the system and varies from zero to one.In this equation, the first four terms on the right hand side are the kinetic energy of the real particles, potential energy of interaction between real particles, kinetic energy of the fractional particle, and the potential energy of interaction between real particles and the fractional particle.The number extension variable, , and the temperature extension variable, s, have their own kinetic and potential energy terms, shown as the last four terms on the right hand side of Eq. ( 12).In fact U  is the work required to carry a fractional particle from a medium of zero potential to the medium of interest with the potential energy U; i.e., it is the so-called excess Helmholtz free energy.
Adopting one of the real particles (molecules) in the simulation box as the fractional particle, the fractional particle is grown to a full (real) particle or deleted from the system in a dynamical way, by solving the equations of motion.The relevant equation of motion, governing the dynamics of fractional particle (the magnitude of coupling to the ideal gas reservoir) is (Eslami & Müller-Plathe, 2007b): 1 where d λ   =d 2 /dt 2 and  id is the chemical potential of the ideal gas.In fact, combination of the last two terms on the right hand side of Eq. ( 13) is the target excess chemical potential.By adding/deleting the fractional particle into/from the simulation box, decision is made either to choose another real particle in the box as the new fractional particle, or to add a new fractional particle to the box.This is done according to the algorithm by Eslami and Müller-Plathe (2007b).On adding the new fractional particle to the box, the coupling parameter, , is set close to zero.Similarly, when a real particle is chosen as the new fractional particle,  is set close to 1. Repeating such a procedure, equilibrium is achieved, in which the density (number) of particles in the system oscillates about an average value, corresponding to the target chemical potential, temperature, and volume.

Phase equilibrium calculation
There are several methods for the calculation of phase coexistence points using molecular simulations, such as thermodynamic scaling method (Valleau, 1991;Valleau and Graham, 1990), histogram reweighting method (McDonald and Singer, 1967;Potoff and Panagiotopoulos, 1998), the Gibbs-Duhem integration method (Kofke, 1993a,b), NpT plus test particle method (Boda et al., 1995;1996), various extensions of it to other ensembles (Moller and Fischer, 1990;Boda et al., 2001), and the Gibbs eensemble MC method (Panagiotopoulos,1987).In the Gibbs eensemble MC method, which is one of the most popular methods for the calculation of phase equilibria, the simulation box is divided into two subsystems and simultaneous simulation of both subsystems are done.To establish the equilibrium (the equality of chemical potentials), particles are exchanged between the two subsystems.This technique has been applied to coexistence properties of simple systems, such as fluids of spherical Lennard-Jones or Yukawa particles (Panagiotopoulos et al., 1998;Rudisill and Cummings, 1989), as well as more complex systems, such as polyatomic hydrocarbons (de Pablo and Prausnitz, 1989;de Pablo et al., 1992) and chain molecules (de Pablo, 1995).There are also reports on the mixed methods in which the molecular simulation approaches have been utilized, to calculate the chemical potentials in the condensed phase and the results from equations-of-state predictions are used to calculate the phase coexistence point (Lim et al., 2003), or to calculate the interaction energy parameters of solvent and polymer, in combination with statistical-mechanical theories for the study of phase equilibria of polymer-solvent mixtures (Jang and Bae, 2002).
Many computational studies of the permeation of small gas molecules through polymers, using afore-cited techniques, have appeared, which were designed to analyze, on an atomic scale, diffusion mechanisms or to calculate the diffusion coefficient and the solubility parameters.Most of these studies have dealt with flexible polymer chains of relatively simple structure such as polyethylene, polypropylene, and poly-(isobutylene) (Müller-Plathe, 1991; Pant & Boyd, 1993;Tamai et al., 1995;Pricl and Fermeglia, 2003;Abu-Shargh, 2004).There are however some reports on polymers consisting of stiff chains.Of these we may address to the works by Mooney and MacElroy (1999) on the diffusion of small molecules in semicrystalline aromatic polymers, by Cuthbert et al. (1997) on the calculation of the Henry's law constant for a number of small molecules in polystyrene and the effect of box size on the calculated Henry's law constants, by Lim et al. (2003) on the sorption of methane and carbon dioxide in amorphous polyetherimide.In some of these reports the authors have used Gibbs ensemble Monte Carlo method (Vrabec and Fischer, 1995;Panagiotopoulos, 1987) to do the calculations, and in some other cases alternative methods, like the equation-of-state approaches are employed to describe the gas phase.In fact, as explained above, one needs to satisfy the equality of chemical potentials in both phases.
There are, however, some recent methods employing the grand canonical ensemble simulation formalism (to set the target chemical potentials in equilibrating phases), without the necessity to do simultaneous simulations on two boxes (unlike the Gibbs ensemble MC method) to calculate the gas solubilities in polymers.In the following section we explain these methods.

Grand equilibrium method
In the Gibbs ensemble simulation method one specifies the thermodynamic variables temperature, global composition, and global pressure for the simulation of both phases in separate volumes.Practically, this set of thermodynamic variables is in many cases not convenient and simultaneous simulation of both phases has the disadvantage that fluctuations occurring in one phase influence the other one.Recently a new method, grand equilibrium method, has been developed by Vrabec and Hasse (2002).This method circumvents the afore-cited problems for the study of phase equilibria.The specified thermodynamic variables are temperature and composition and two independent simulations are performed for the two phases without the need to exchange particles in the condensed phase.According to this method for a mixture composed of several components, it is possible to do a simulation in the isothermal-isobaric (NpT) ensemble at constant temperature, a constant composition of the condensed phase, and at an arbitrary constant pressure, preferably close to the pressure at the phase coexistence point, to obtain the density of the condensed phase.
In the grand equilibrium method, a simulation of the condensed phase is done to calculate the excess chemical potentials,  s ex , and the partial molar volumes, V s , of all components.
One may use the test-particle insertion method (Widom, 1963)  where p * is the target pressure at which the NpT ensemble simulation is done.Once the  s ex (p) is determined from Eq. ( 16) by one NpT ensemble simulation of the condensed phase, one vapour/gas simulation has to be performed in the pseudo-grand canonical ensemble (pseudo-VT).In a common grand canonical ensemble (Eslami and Müller-Plathe, 2007b) the parameters temperature, volume, and the chemical potential of all species are fixed, while in this pseudo-VT ensemble simulation, the parameters T and V are fixed in the common way, but instead of fixing the chemical potentials, they are set as a function of the pressure in the gas phase.This procedure ensures that equilibrium between the condensed phase and the gas phase is imposed.In a common VT ensemble simulation (Eslami and Müller-Plathe, 2007b), the chemical potentials are set through insertion and deletion of particles by a comparison between the resulting potential energy change and the desired fixed chemical potential.Here, starting from a low density state point, the gas phase simulation is forced to change its state to the corresponding phase equilibrium state point.

Applications
Recently, we have applied the grand equilibrium method for the calculation of solubilities of gases in polymers.The method has been well tested for the case of gas solubilities in polystyrene over a wide range of temperatures and pressures (Eslami and Müller-Plathe, 2007a).To calculate the phase equilibeium points, MD simulations were performed at a specified temperature and concentrations of sorbed gases in polymer.
After equilibration, several configurations are extracted at different times from the dynamic simulation and used to insert the test particles, to calculate the excess chemical potentials.Calculating the excess chemical potentials of the sorbed gases in polymer as a function of pressure, a simulation of the gas phase is done in the grand canonical ensemble, to find the phase coexistence point.We have shown in Figure 3 typical results for the zero-pressure limit solubility coefficients of methane in polystyrene (Eslami and Müller-Plathe, 2007a).Another more interesting application of the method is in the case of gas solubility calculation for polar penterants in polar polymers.This is a severe test, as most of the existing methods, like the Widom's test particle insertion method or its modifications, are not applicable in the case of polar fluids.We have recently applied this method to the case of water solubilities of poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009) and polyamide-6,6 (Eslami & Mehdipour, 2011).Preparing an initially relaxed configuration of the polymer, we have performed a long simulation, in the grand canonical ensemble, of the polymer phase to get the density of water, which corresponds to the target values of , V, and T. The number of molecules (density) during the simulation changes until achieving equilibrium, at which the density fluctuates about the average value.We have shown in Figure 4 the variation of number of water molecules in the simulation box at two different temperatures (one above and one below the glass transition temperature).
In Figure 4, the results at T= 280 K belong to  ex s = -18.5 kJmole -1 in a cubic box with a box length, L, of 3.082 nm.Calculating the partial molar volume of water from the results of this simulation, we expanded the chemical potential of water in the polymer phase as a function of pressure, according to Eq. ( 16).Then we have performed another simulation of the gas phase in the grand canonical ensemble, as explained above.Here the target chemical potential is varied during the simulation as a function of running-average pressure in the gas phase (see Eq. ( 16)).The result for time evolution of the density of water in the gas phase (L=50 nm) is also shown in Figure 4.The results in Figure 4 show that the gas phase quickly adjusts its state to the phase coexistence point.Also shown in the same figure are the results of similar calculations at 450 K ( ex s = -9.0kJmole -1 ).Here the values of L for simulations of polymer and gas phases are 3.128 nm and 7.0 nm, respectively.Similar calculations are done in the case of water solubility of polyamide-6,6 (Eslami & Mehdipour, 2011) over a wide range of temperatures and pressures.This is a particularly interesting example of the applicability of the method, since polyamide-6,6 is a highly polar polymer matrix with the ability of forming a strong hydrogen bond network and can dissolve water to about 10 wt %.We have shown in Figure 5 a typical calculated sorption isotherm for water in polyamide-6,6 at 300 K and its comparison with experiment.

Diffusion in polymers
Most theories describing the mechanism of diffusion in polymeric materials are based on the free-volume approximation (Vrentas & Duda, 1979).In the free-volume theories, there is a volume which is directly occupied by polymeric molecules, and there is the remainder of the volume, which is called the free volume.A part of the free volume is assumed to be uniformly distributed among the molecules and is identified as the interstitial free volume, which requires a large energy for redistribution and is not affected by random thermal fluctuations.The other part, which is called the hole free volume, is assumed to require negligible energy for its redistribution.Therefore, the hole free volume is being continuously redistributed due to random fluctuations, and is assumed to be occupied by penetrant molecules.This redistribution of hole free volume will move the penetrant molecule with it.According to this model, by movement of segments of the polymer chain, a void will be created adjacent to the penetrant molecule.If the size of this hole is sufficient to host a penetrant molecule and if the penetrant have sufficient energy to jump into the hole, a successful jump of the penetrant molecule is made into the hole.Although the free volume model has been used extensively to describe the mechanism of transport through molten or glassy polymers (Mauritz & Storey, 1990), this model does not show a microscopic view point of penetrant transport in polymers, since it just connects bulk transport properties, like diffusion coefficient, into bulk properties, like molecular volume or thermal expansion coefficient.
The transition-state theory (TST), introduced by Gusev et al. (1993), is another useful method for the calculation of diffusion coefficient of a low-molar-mass substance through the polymer matrix.In the TST it is assumed that the movement of the penetrant from an initial cavity to the saddle point and to a neighboring cavity is a unimolecular rearrangement.For such a transition the reaction trajectory in configuration space is tracked and the transition rate constant is evaluated.In the first studies on the application of TST to study the dynamics of light gases dissolved in rigid microstructures of glassy polycarbonate and rubbery polyisobutylene Gusev et al. (1993), the method was shown to be only capable to study just the dynamics of light gases like He.The method was developed by Gusev and Suter (1993), by Gusev et al. (1994), and by Karayiannis et al. (2004), to calculate the diffusion coefficient of bigger penetrants in glassy polymers.
Computer modeling of molecular systems at a detailed atomistic level has become a standard tool in investigation of sorption and diffusion of small molecules in polymeric media (Müller-Plathe, 1991, 1994, Milano et al., 2002, Mozaffari et al. 2010).MD simulation is a useful tool for exploring the structure and properties of bulk amorphous polymers.The length of the trajectories that can be generated in practice presently is on the order of many nanoseconds.Thus the range of properties that can be studied directly is limited to those that evolve over this time scale.One of the phenomena that appear to be suitable for investigation is the diffusion of small penetrant molecules in an amorphous polymeric matrix.That is, the diffusion coefficients of small penetrants in many rubbery or liquid polymers are such that, at temperatures close to room temperature and above, the average displacement of the diffusant is large enough in a nanosecond interval to be determined via molecular dynamics simulation.Performing such simulations is of practical importance in predicting diffusion coefficients and also in understanding the mechanism of diffusion.
Although there has been significant progress in the use of MD methods in the simulations of diffusion coefficients, early studies were focused on the simulation of gas diffusion in rubbery polymers which could be investigated using full atomistic or united-atom simulations in reasonable computational times (Boyd, 1991;Sok et al., 1992, Pant & Boyd, 1993;Gee & Boyd, 1995).Due to the recent development of improved force fields and the wider availability of sophisticated commercial softwares and high-speed computing facilities, attention is shifting to direct to the more challenging task of simulating the slower diffusional processes occurring in glassy polymers (Han & Boyd, 1996, Milano et al. , 2002, Lime at al., 2003, Mozaffari et al., 2010).In the following section we present the MD simulation results done recently by the authors on the diffusion of gases in polystyrene over a wide range of temperatures and pressures.

MD simulation of diffusion in polymers (case study)
Recently we have studied the diffusion mechanism of some gases in polystyrene over a wide range of temperatures, ranging from below the glass transition temperature to temperatures well above the glass transition temperature (Mozaffari et al., 2010).Center-ofmass mean-square displacements have been measured during the simulation to calculate the diffusion coefficients.In the limit of long times, which the penetrant molecules perform random walks in the polymer matrix, the mean-square-displacement becomes linear in time, and the diffusion coefficients can be calculated using the Einstein relation: The diffusion coefficients of penetrant gases in polystyrene have been calculated over a wide range of temperatures, 300 K-500 K.The motion pattern of penetrant gases in host polymer can be qualitatively studied by monitoring the penetrant's displacement |r(t)-r(0)|, from its initial position.We have shown in Figure 6 the displacements of argon and propane in polystyrene at 300 K.The curves are representative of a common hopping mechanism, showing that for a considerable time interval the penetrants dwell in existing voids in the polymer and occasionally do a jump into the neighboring voids.When dwelling in the voids, the penetrants just perform oscillatory motions around their equilibrium positions, therefore, no net motion of a penetrant molecule occurs with these positional fluctuations.The amplitude of oscillations varies according to the size of voids.From time to time, the penetrants can do a quick jump into their neighboring voids, see Figure 6.The jump frequency depends on penetrant's size, therefore the bigger penetrants (see Figure 6 in the case of propane) can rarely jump between the voids.
Particularly useful information can be obtained by analyzing the trajectories of penetrants in the polymer.The two-dimensional center-of-mass x-y trajectories of nitrogen, as a typical example, in polystyrene at temperatures below (300 K) and above (500 K) the glass transition temperature are indicated in Figure 7.
The results in Figure 7 are indicative of fas t e r m o v e m e n t o f p e n e t r a n t s a t h i g h e r temperatures, as represented by the broadened range of displacements at higher temperatures.This indicates that at higher temperatures the hole free volume redistributes faster and penetrant molecules have higher energy to overcome the activation energy required to jump into new voids.As a measure of diffusion coefficients, the center-of-mass mean-square displacements of carbon dioxide at temperatures below (340 K) and above (480 K) the glass transition temperature, are shown in Figure 8.The linearity of mean-square displacements versus time, as indicated in Figure 8, confirms Einstein diffusion.In the glassy polymer at times below 500 ps, the penetrant motion is highly anomalous and the diffusion regime begins at longer times.Intercavity jumps rarely occur at this time scale.However, as it can be seen from the results in Figure 8 at high temperatures, the diffusion regime sets in a shorter time.These findings clearly indicate the difference between the diffusion mechanisms in the glassy and melt polymers.

Permeation in polymers
For a permeant gas subject to a pressure gradient at two different sides of a membrane, the permeability, P, is defined as quantity of permeant film thickness P area time pressure drop across film The quantity of permeant gas is often expressed as the volume of gas at STP condition.In fact, the pressure difference between up-and down-stream sides of polymer membrane is the driving force of permeation phenomenon.As stated above, penetrant transport through membrane is described by solution-diffusion mechanism.According to this model the permeation occurs in three steps: 1-sorption of gas molecules at upstream side of polymer 2diffusion of the gas molecules into the polymer matrix 3-desorption of gas molecules at downstream side of polymer.This is indicated schematically in Figure 9.
Fig. 9.The schematic of gas permeation in a membrane The permeability coefficient of gas molecules across a polymer film of thickness L is described by the flux J: where Δc is the concentration gradient of penetrant molecules through the membrane.
Accordingly, the permeation coefficient can be expressed by the product of the diffusion coefficient D, and the solubility coefficient S, as: Since direct simulation of the permeation process is not possible due to the system size, Eq. ( 20) is a particularly useful expression connecting the permeability coefficient with the diffusion and solubility coefficients.Therefore, diffusion and solubility coefficients, calculated in MD simulation methods, are used to calculate the permeability coefficients according to Eq. ( 20).
We have shown in Figure 10 our recently calculated results on the permeability of water in poly(ethylene terephthalate) (Eslami & Müller-Plathe, 2009).The results, indicating the permeability coefficients over a wide range of temperatures, are shown in Figure (10).The Arrhenius plot of ln(P) vs. 1/T shows a break at the glass transition temperature.As a result of calculating a higher solubility coefficient, our calculated permeation coefficients are also higher than the experimental values.Solubility coefficients are higher than the experimental values, especially at lower temperatures.This is trivial for calculations of this type.Similar differences between experimentally and computed solubility coefficients have been observed in previous studies (Cuthbert et al., 1997;Lim et al., 2003, Eslami & Müller-Plathe, 2007a, 2009;Eslami & Mehdipour, 2011).According to Gusev and Suter (1993), errors of the order of (2-4)k B T are common for the calculation of Helmholtz energies by molecular simulations.It has been speculated that the main contribution to the solubility comes from single holes in the simulated polymer structure, which might not be present in similar proportion in real polymers (Müller-Plathe et al., 1993).Errors of this type become even more serious in cases where the solubility is very low.
We have also compared our calculated permeability coefficients of a number of gases in polystyrene (Mozaffari et al., 2010) with experimental values (Csernica et al., 1987) and with simulation results by Kucukpinar and Doruker (2003).The results in Table 1 indicate that the calculated values are in good agreement with experiment.The results indicate that polystyrene is much permeable to CO 2 , compared to other gases studied by Mozaffari et al. (2010), because of the higher solubility coefficient of CO 2 in polystyrene and its relatively higher diffusion coefficient.On the other hand, polystyrene is less permeable to the biggest penetrant molecule studied, propane, because of its very small diffusion coefficient.
An interesting test is to compare the calculated permeability coefficient ratios (selectivities) in the zero pressure limit with the corresponding experimental values (Csernica et al., 1987).This is done in Table 2, in which we have listed the ratios of permeability coefficients at 300 K to that of nitrogen's permeability and compared the results with experimental measurements (Csernica et al., 1987) and with the calculations by Kucukpinar and Doruker (2003).The results in Table 2 show that the calculated ratios are quite close to the experimental ratios.This shows that our calculated permeability coefficients are higher that the experimental values by nearly the same factor.The same conclusion was made in our work on the calculated solubility coefficients (Eslami & Müller-Plathe, 2007a).This indicates that there is an excessive free volume in the polymer sample compared to the experimental polymer samples.Better agreement of the results by Mozaffari et al. (2010) with experiment (Csernica et al., 1987) compared to those of Kucukpinar and Doruker (2003) has been attributed to using a united atom models for both polymer and diffusant gases, in the later case.

Summary
In this chapter the polymer permeability is discussed and reviewed from a computational point of view.Although it is impossible to directly simulate the permeability process, the simulation techniques are powerful tools to simulate and to give a molecular level insight to the solubility and diffusion mechanisms of gases in polymers.While there are lots of simulations in the literature, giving address to the diffusion mechanism of gases in glassy and/or rubbery polymers, the computations are less straightforward in the case of gas solubilities.Therefore, the main problem in calculation of gas perameability of polymers, using molecular simulation methods, stems from the calculation of solubility coefficients.In fact the calculation of gas solubilities necessitates the condition of equilibrium between the permeant in the gas/liquid phase and the permeant dissolved in polymer.The equality of the temperature, pressure, and chemical potentials of all species is the necessary condition to establish such an equilibrium situation.The chemical potential is, however, coupled to the number of particles and cannot easily be calculated employing molecular simulation methods.
The applicability and feasibility of different techniques for the calculation of chemical potentials is evaluated and discussed.It is explained that in some straightforward cases, like the thermodynamic integration method, several simulations are needed to calculate the chemical potential.Widom's test particle insertion method (Widom, 1963) is of practical importance in many cases, but the method encounters sampling problems at high densities and in the case of big solutes.The method is also not applicable in the case of chemical potential calculation for polar species.It is shown that recent grand canonical MD simulation methods (Eslami and Müller-Plathe, 2007b;Eslami et al., 2010).are trustable methods for the sake of phase equilibrium calculation.In these methods the chemical potential is set as an independent thermodynamic quantity, and the number of particles in the box is changed in a gradual and dynamical way.The applicability of the method to the challenging cases like gas solubilities in polymers with stiff chemical structures (Eslami and Müller-Plathe, 2007a) and solubility of water (as a polar penetrant) in polar polymers, like polyamide-6,6 (Eslami & Mehdipour, 2011) is reviewed and discussed.
The process of gas diffusion in polymers and different approaches for studying the gas diffusion are addressed to.It is explained that MD simulation well describes the jumping mechanism of diffusions in glassy polymers.According to this mechanism, in a glassy polymer the penetrants dwell in the polymer holes and occasionally do a jump into the neighboring voids (Mozaffari et al., 2010).Therefore, no net motion of a penetrant molecule occurs while the penterant just performs random oscillations inside a hole.From time to time, the penetrants can do a quick jump into their neighboring voids, with a jump frequency depending on penetrant's size, temperature, and chemical structure of the polymer.The permeation coefficients of gases in polymers are, therefore, reviewed in this chapter in terms of gas solubilities and diffusions.

Fig. 2 .
Fig. 2. A schematic one-dimensional representation of a free-energy landscape for an additive molecule (x is a component of its center of mass) in a polymer matrix as a function of the coupling parameter .During a fast-growth simulation, a molecule is forced along  from state A to B, while sampling in x; three example paths are given.With fast-growth thermodynamic integration, it is possible to explore the multiple minima in the free energy landscape by multiple trajectories starting from a flat free energy landscape (a solute not interacting with the matrix).Figure is taken from Hess and van der Vegt (2008) with permission.

Fig. 4 .
Fig. 4. Time evolution of the number of water molecules in the grand canonical ensemble simulation at 280 K and 450 K.The solid and dotted curves represent the results for water molecules in poly(ethylene terephthalate) and in gaseous phases, respectively.

Fig. 6 .
Fig.6.Displacement of argon (upper curve) and propane (lower curve) molecules from their initial positions at 300 K.In order to avoid overlapping between the curves, the displacements of argon are offset by 0.6 nm. Figure is taken fromMozaffari et al. (2010).

Fig. 8 .
Fig.8.Center-of-mass mean-square displacement for carbon dioxide at temperatures below and above the glass transition temperature.The curve at 480 K is scaled down by a factor of 15 for the sake of clarity.The dashed lines show the least-squares fits to the linear parts of the curves.Figure is taken fromMozaffari et al. (2010).