Chapter 1: Multi-scale Molecular Modeling Approaches for Designing/Selecting Polymers used for Developing Novel Membranes

During the last decade computational chemistry and numerical simulations have had a favorable impact in almost all branches of materials research ranging from phase determination to structural characterization and property prediction.^1–7  An important effort has been focused on developing simulation tools to describe thermodynamic and transport properties of confined fluids.^7–14  The present contribution illustrates the benefit of coupling experiment to molecular modeling for selecting novel membrane materials with better separation properties for given gas mixtures as well as the limitations of the existent computational methodologies. New modeling and simulation tools based on multi-scale hierarchical modeling are needed to cope with the complexity of materials and associated phenomena at different length and time scales.^2–4 

Transport properties of small molecules in amorphous polymer matrices play an important role in many industrial applications such as gas separation of mixtures, packaging applications ranging from food conservation to controlled drug and cosmetics release, to special coatings for protecting specific substrates from gases.

Different aspects of technology and industrial application of polymer membranes from materials research to permeator design to their optimal configuration to enhance processes performance have been discussed in detail in a previous lecture and recently reviewed in the literature.^15–17  The potential application of a polymer as a separation membrane depends upon the selectivity towards the gas to be separated and the permeate flux. The selectivity determines the product purity and recovery whereas the permeability is related to the productivity of the membranes. This means that both the permeability and the selectivity should be as large as possible. The control of gas permeability and permselectivity of polymer membranes has become a subject of active research with worldwide participation in both industrial and academic laboratories.¹⁷  The design and optimization of polymer membranes used in gas separation applications would be possible if reliable predictions of transport properties could be made rapidly in advance of synthesis and experiment. The actual status of available commercial software for modeling transport phenomena in polymer membranes,^{2–4,18–20}  does not allow the development of de novo material design approach. This is due not only to formidable time and length scales involved, but also to lack of detailed information on time evolution of the free volume and its distribution as a function of processing history during the manufacturing process. Rapid progress in computational methodology and validation of new simulation tools is improving the understanding of different facets of gas transport in polymer membranes and building the necessary tools for their effective use in materials design.^{1–4,21–23}  The possibility to predict transport properties of small molecules through polymer matrices permits the rational selection of polymer materials used in these applications and their optimal design. Although there has been reported an increasing number of studies on this subject over the last years,^24–43  the prediction of transport properties of gas molecules through glassy polymer membranes remains a difficult target. In many of the recent studies reporting molecular simulation predictions of diffusion and solubility of small gas molecules in several membrane models of the same glassy polymer a great scatter of the predicted values is observed. These results clearly indicate that the quality of the packing of the polymer chain into an amorphous cell membrane model strongly impacts the predicted gas transport properties.

The potential application of a polymer as a separation membrane depends upon the possible throughput and the purity of product.^{15–17,44,45}  This means that both the permeability of the gas that is transported more rapidly and the selectivity should be as large as possible. The permeability coefficient, Pe, of a small molecule through a polymer membrane is defined as:

the product of the diffusion coefficient, D (kinetic parameter), and of the solubility coefficient, S (thermodynamic parameter). The estimation of these coefficients can be done, either by molecular dynamics (MD) and grand canonical Monte Carlo simulations, or by applying the transition state theory (TST) approach provided that the quality of the membrane amorphous cells used in the calculation represent the real distribution of torsion angles, of the free volume and its distribution, as well as the structural, conformational and volumetric properties of polymer membranes. The selectivity of a polymer membrane for a pair (i, j) of gas molecules is characterized by the ideal separation factor α_ij defined in eqn (1.2):

Following this definition the selectivity of a membrane is the product of diffusion selectivity (D_i/D_j) and of solubility selectivity (S_i/S_j). In the case of glassy polymer membranes the overall selectivity is mainly controlled by diffusion selectivity. Two types of membranes are used commercially in gas separation technology. Glassy polymer membranes are made from stiff chain polymers and operate below their glass transition temperature. These membranes have moderately high free volume and separate gases predominantly based on differences in the sizes of the gas molecules. The smaller molecules (H₂, He, O₂) permeate more easily through the membrane than the larger ones (CH₄, C₂H₆, N₂). The second class of membranes is made either from highly flexible rubbery polymers or ultra-high free volume, glassy substituted polyacetylenes.^16,42  These membranes separate gases principally by differences in the solubility of gas molecules in these polymers. Larger and more soluble penetrants permeate faster than smaller and less soluble ones.

In the past 20 years, the control of gas permeability and permselectivity of polymer membranes has become a subject of active research with worldwide participation in both industrial and academic laboratories.^{15–17,44–51}  However, it has been found that simple structural modifications, which usually lead to an increase in polymer permeability, cause loss in permselectivity and vice versa. This so-called ‘trade-off’ relationship has been well described in the literature.^{17,44,45,48–51}  Here the log of the separation factor α versus the log of the higher permeability gas Pe yielded a limit for achieving the desired results. The upper bound limit is not fixed in the α–Pe space but moves with time as new polymers with optimized structures become available.

Recent studies of glassy polymer membranes indicate that in addition to the free volume content, gas transport parameters depend upon the backbone chain rigidity, its segmental mobility, the inter-chain distance and the chain interactions. For example, the introduction of bulky alkyl substituents opens up the polymer matrix resulting in a greater permeability. Also the reported introduction of n-alkyl side groups on a polymer backbone increases the side-chain flexibility as well as the membrane free volume with an overall increase of the permeability.^17,44,50,51  The free volume of a polymer, which corresponds to the unoccupied regions accessible to segmental motions, is an important parameter for understanding and predicting many of its characteristic properties. The free volume and the free volume distribution influence the molecular mobility and the transport properties of low molecular substances and gases in polymers.^52–55  Atomistic simulations allow a detailed investigation of these geometric characteristics. A complete description of definition of the free volume of membrane atomistic model configurations as well as detailed approaches for analyzing it has been reviewed recently.¹⁸  Voronoi polyhedra and Delaunay simplices are used for analyzing the local arrangement of free volume in glassy polymer membranes cells. The calculated free volume and its distribution are used to describe the gas solubility and diffusion in amorphous polymer cells through a QSPR correlation as well as to control the quality of the constructed amorphous cells.

Many literature examples that use the QSPR approach associated with polymer permeability are based on the use of group contribution methods to establish a correlation between the structure of the repeat unit and some physical property of the polymer membrane such as the free volume, the mean segment distance or dielectric constant polarizability, which in its turn is used to predict permeation properties. As a result of such studies some practical criteria have emerged to guide synthetic researchers in improving the permselectivity of membranes which have evolved through extensive experimentation: (i) inhibition of inter-segmental packing meanwhile simultaneously inhibiting intra-segmental (backbone) mobility; and (ii) weakening of inter-chain interactions (reduction of charge transfer complexes). These design rules are based on phenomenological paradigms that provide guidelines for polymer selection. The QSPR approach that uses appropriate descriptors for representing at the same time the repeat unit, the unperturbed polymer chain and the packed amorphous cells, integrating data at different length scales including the some data from the processing history of the polymer membrane would lead to new criteria for materials improvement.

This contribution focuses on the need to shorten the research time of novel polymer materials used for membrane fabrication by combining several computational approaches and experimental techniques. It does not pretend to be a review or cover all aspects that might be considered traits of the computational materials design. Various examples are used to illustrate the use of existing numerical simulation and modeling tools for complementing experimental work. After a brief discussion of some computational approaches used to cover different aspects of polymer membrane simulation, methodologies used to characterize gas transport through polymer membranes as well as to identify factors that control thermodynamic and kinetic properties of confined fluids in membranes will be detailed. It is noticed that the successful application of modeling approaches to gas separation by membrane technology needs the development of models dealing with multi-component gas mixture transport through model membranes. Moreover, for a given polymer membrane, both gas diffusivity and gas solubility depend strongly on process parameters such as pressure difference, feed composition, and temperature. More information on the effects of process parameters on selectivity contribution should be thoroughly considered for identifying membrane materials suitable for each application.

Progress in numerical simulation techniques has followed the developments in chemical engineering science and recently has received tremendous attention in both scientific and industrial communities due to the possibility to integrate micro-scale information in investigations that can be carried out at various levels of resolution on different scales of time and length. Several multi-scale modeling methodologies that are based on the information transfer between different scales starting from the molecular level and ending up at the industrial scale have been recently reviewed in the literature³  indicate a clear trend towards coupling of the design of chemical engineering equipment/units with nanoscale modeling. It is generally expected that multi-scale modeling can lead to optimal unit design as well as to cost optimization of the final products. Fundamental research in materials and structure design now considers numerical simulation, as a complement to theory and experimentation. The integration of physical testing, advanced computations and system simulation would dramatically reduce the design and development time and costs.

In this contribution we do not aim to provide a detailed description of each of the numerical techniques contributing to the multi-scale modeling methodology but will only outline its basic principles, strengths and weaknesses, and potential applications. Interesting readers can refer to relevant books, reviews and research articles for details.^56–58 

The terminology used for characterizing multi-scale methods often varies with the application domain. Two general approaches have been developed for integrating different models at disparate spatial and temporal scales. The hierarchical modeling extracts information from lower-scale models and transfers it as parameters to the upper-scale (coarser) models of overlapping domains. The coarse-scale model is used over the entire computational domain but a higher-resolution modeling is applied to zoom into a particular sub-domain to obtain updated parameters within the corresponding grids of the coarse model.

Mesoscale modeling uses a basic unit just above the molecular scale, and is particularly useful for studying the behavior of polymers and soft materials. It can model even larger molecular systems, but with the commensurate trade-off in accuracy.^59–61 

Two methods dissipative particle dynamics (DPD) was initially devised by Hoogerbrugge and Koelman⁵⁹  as a particle-based off-lattice simulation method for the flow of complex fluids and to tackle hydrodynamic time and space scales beyond those available with MD. Since DPD is a coarse-grained model and individual atoms or molecules are not represented directly by the particles but they are grouped together into beads, these beads represent local ‘fluid packages’ able to move independently.

The fundamentals of the DPD method are now fairly well-established,^62–67  as are the technical subtleties^68,69  and the coarse-grained parametrization of the DPD particles, with improvements consistently being introduced.⁷⁰ 

Furthermore, it is expected that DPD will play an ever-increasing role in multi-scale modeling approaches through bridging of the atomistic and continuum scales. In such approaches, atomistic simulations are performed to build the DPD models, followed by DPD simulations which provide the necessary input to the continuum codes.

DDFT, which was developed by Fraaije et al. in 1997,⁶²  is a field-based theoretical method for studying complex fluids, their kinetics and their equilibrium structures at micrometer length and microsecond time scales.^71–73  DDFT has been applied to the study of the self-assembly of block copolymers in bulk, under shear and in confinement,^74–76  and to study polymer blend compatibility.⁷⁷  Compared to the DPD method, DDFT is computationally extremely fast since larger elements can be modelled. Moreover, since the fluid elements can freely penetrate, larger time steps can be used, and furthermore it is less likely to become trapped in a local minimum. Since DPD is a particle-based method, it can provide somewhat more detailed structural information. Nonetheless, they are both powerful tools in simulating phase separated phenomena that occurs at the mesoscale and the consistency of results from the two methods for the same coarse-grain model is evaluated in this work.

The alternative to the hierarchical modeling is known as concurrent modeling and consists of combining different numerical models that simultaneously describe different sub-domains. Each numerical model runs in its sub-domain and exchanges information with its parent/child sub-domain at defined boundaries. The critical issue here is to define the criteria and protocols in order to automate the application of more detailed numerical models in a simulation domain in time and realizing the so-called temporal multi-scale calculations or/and space referred as spatial multi-scale ones.

Spatial multi-scale methods are based on the paradigm that in many real situations the atomic description is only required within small parts of the simulation domain whereas for the majority the continuum model is still valid. This allows one to apply concurrent continuum and molecular simulations for the respective parts of the simulation domain using a coupling scheme that permits to connect between the two domains. The majority of the spatial domain is calculated by continuum solvers (computational fluid dynamics) which are very fast and only the ‘active part’ is calculated using molecular simulation methods. In some cases several other coarser-grained (mesoscale) methods than the atomic simulations ones are used as interfaces between the molecular simulation and the continuum domains. Such approaches are called hybrid molecular–continuum methods and allow the simulation of problems that are not accessible either by continuum or by pure molecular simulation methods.

Successful examples of application of hybrid molecular–continuum approaches are the fluid slippage past surfaces,⁷⁸  fracture propagation,⁷⁹  problems involving phase transitions,⁸⁰  fouling at surfaces, wetting, moving contact line between two immiscible liquids,⁸¹  mass transport through membranes, systems with locally high stresses or strong gradients, flows in micro- or nanochannels with specific surface features, flows along surfaces with no ideal non-slip or slip condition or where the liquid molecules interact with wall molecules.

Modeling of gas separation process using membranes implies, as schematically reported in Figure 1.1, using different numerical simulation techniques going from the atomistic representation of the amorphous polymer matrix representing the thin layer of the active membrane to the hydrodynamics of fluid moves upon the membrane module. The computational fluid dynamics (CFD) simulations can reveal details of flow patterns, velocity distributions and resistances in an industrial module made of hundred thousands hollow fibers and correlate the pressure distribution upon membrane surface with the fluid flow pattern on the hollow fibers. However, in more realistic cases such as existence of the polarization layer the CFD calculations need to include details from lower scales in order to realistically predict permeate flux through the membrane module by integrating of modified boundary condition in the simulation and predicting the concentration polarization profile that is developed at the membrane fluid interface. The above described schematic representation of the hybrid atomistic-continuum is only one of the aspects of the numerical simulations applied to the R&D activity in membrane technology. The numerical simulation approach for predicting the morphology of the membrane is mainly based on a combination of atomistic and mesoscale methods. The final morphology of the membrane depends on the properties of the polymer materials as well as on the processing conditions. Most commercial membranes are prepared by the immersion precipitation process. Recently, it has been showed experimentally that the process of membrane formation proceeds via spinodal decomposition to form the asymmetric membrane.⁸²  Lattice Boltzmann simulation of membrane formation in two dimensions⁸³  was able to capture the morphology of the interface between coagulation bath and polymer solution as well as the asymmetric morphology of the membrane. In a recent paper⁸⁴  the simulation has been revisited using a phase field model of the immersion precipitation process and membrane formation in three dimensions. The simulation approach looks very promising in terms of providing a deeper understanding of the immersion precipitation process and theoretical guidance on experiments design membrane formation.

Another class of spatial multi-scale methods concerns the quantum chemistry community where efforts have been focussed on the combination of quantum mechanics (QM) methods with continuum electrostatic theories in order to realistically represent the solvation free energy in a polar environment. These methods have been refined over the years and can now give a reasonable description of solvation properties of an isotropic and homogeneous medium.^85–87  However, these continuum models are not appropriate to represent the electrostatic and steric interactions of the structured environment with the active site. This is particularly true in the descriptions of complex systems like enzymes or catalysts. An appropriate description of such systems has been developed using a hybrid quantum mechanical/molecular mechanical (QM/MM) approaches^88–91  where the QM methods are used to describe the active site where chemical reactions or electronic excitations occur, and MM methods are employed to capture the effect of the environment on the active site.

In this contribution we will focus mainly on the methodology used for atomistic simulation of gas transport through models of amorphous together with the generation of atomistic models of polymer membranes and the prediction of gas transport properties of gas molecules through the membrane and briefly mention the recent contributions on coupling such calculations to coarser or continuum-based methods.

The modeling and simulation methods at molecular level usually employ atoms, molecules or their clusters as the basic units considered. Atoms or molecules interact with each other through a force field, or intermolecular potential energy, and the accuracy of the force field directly determines the accuracy of the resulting calculations. The common simulation methods dealing with many-body systems can be divided into stochastic and deterministic ones. The first class is represented by the Monte Carlo (MC) method whilst the second one by the molecular dynamics (MD) method. Modeling of polymer membranes at the atomistic scale is predominantly directed toward the realistic representation of atomistic models of amorphous cells representing the membrane layer. However, the applicability and effectiveness of these methods often depends on the ability to fulfil large-scale computations; and despite the availability of high-performance computers, these methods are still restricted to solving systems that are too small for even nanoscale problems. Computer-aided molecular design of polymer membrane models at a detailed atomistic level has been reported in the literature^{24–26,29–43,92–99}  for investigating the sorption and diffusion of small gas molecules. MD simulation is used for exploring the structure and properties of bulk amorphous polymers and the diffusion of the penetrant molecules to be followed by exploring the trajectories generated during the simulations. MC simulations probe the configuration space of the model membranes by trial moves of particles either during the phase of amorphous cell building or during the simulation of sorption properties of gas molecules inside the amorphous cell using the so-called Metropolis algorithm to monitor the energy change from in successive steps as a trigger for accepting or rejecting a new configuration. Configurations resulting with a lower energy are accepted whereas those resulting with a higher energy are accepted with a probability governed by Boltzmann statistics. The algorithm ensures the correct limiting distribution and properties of a given system can be calculated by averaging over all MC moves within a given statistical ensemble.

Molecular dynamics (MD) is an atomistic simulation method for studying a wide class of materials, such as polymers, metals, ceramics and biomolecules under ambient as well as extreme conditions. MD allows one to predict the time evolution of a system of interacting particles (e.g. atoms, molecules, granules, etc.) and estimate the relevant physical properties.^100–102  Specifically, it generates such information as atomic positions, velocities and forces from which the macroscopic properties (e.g. pressure, energy, heat capacities) can be derived by means of statistical mechanics. MD simulation usually consists of three constituents: (i) a set of initial conditions (e.g. initial positions and velocities of all particles in the system); (ii) the interaction potentials to represent the forces among all the particles; (iii) the evolution of the system in time by solving a set of classical equations of motion for all particles in the system. MD methods are governed by the system Hamiltonian and the Hamilton's equations of motion are integrated to move particles to new positions and to assign new velocities at these new positions.

Given a force field for the potential energy, the Hamiltonian of a system of N atoms can be written as:

where it is assumed that the kinetic energy, K, depends only on the momenta (p) and it is separable from the potential energy, ϕ, that depends only on atomic positions.

Particles in MD move naturally under their own intermolecular forces and follow Newton's second law:

where m_i, and r_i are the mass, acceleration and position of particle i, respectively. During the simulation the configuration space as well as the phase space is explored allowing extracting information on dynamics of the system. In order to simulate the gas diffusion in a polymer membrane a force field representing the interactions between all the atoms of the system composed of the polymer amorphous cell and penetrant molecules is required. The force field has to be tested against experimental results and theoretical constraints.

The gas diffusion coefficients can be estimated either from MD simulations by using the Einstein (eqn (1.5)) or by means of the Green–Kubo (eqn (1.6)) formulations:

The mechanism of diffusion in rubbery polymers is different from the glassy ones. The diffusion coefficients for small gas molecules in rubbery polymer membranes do not depend on concentration whilst in glassy polymer membranes they do depend and reach a constant value at relatively high concentrations. This is mainly due to the fact that glassy polymers are not in a thermodynamic equilibrium state. For these polymers the final ‘metastable’ chain configuration depends on the processing history of the membrane. This detail makes even more difficult the modeling of glassy polymer membranes due to lack of experimental structural data for validating computational approaches. The prediction of self-diffusion coefficients for non-polar small gas molecules in amorphous rubber polymer matrices is normally done through MD calculations.^27,103–105 

After construction of the amorphous cell using the method of Theodorou-Suter¹⁰⁶  and geometric free volume analysis of the cell several (four to six) penetrant molecules to improve sampling are inserted at the free volume positions. The cell is further relaxed by 100 ps of a NPT-MD (constant particle number, temperature and pressure) simulation at 1 bar and room temperature before starting a longer (nanoseconds) NVT dynamics. The recorded trajectories of each penetrant gas molecule are analyzed and the diffusion coefficient is determined by means of relation (eqn (1.5)). In Figures 1.2 and 1.3 the packed cell model of polydimethylsiloxane (PDMS) and the trajectory of N₂ molecules in the PDMS matrix are reported. The MD simulations show two types of motions of the N₂ molecules: jumps between cavities and local motion inside cavities.

The predicted self-diffusion coefficients depend principally on the quality of the force fields used to model the interactions not only between the penetrant and polymer matrix, but also intramolecular interactions between polymer chains. These last ones, strongly affect the quality of the amorphous polymer cell and in particular the total free volume its distribution and dynamics which in their turn affect the predicted values of diffusion coefficients. The role of chain relaxation and matrix fluctuations in explaining the diffusion mechanism of small gas penetrants as N₂ in rubber polymer membranes has been clearly demonstrated through MD calculations in which the polymer matrix has been kept fixed.^27,103–105  MD simulation of gas diffusion in polymer membranes generates a wealth of information on the mechanism of gas transport but its use is limited to high free-volume rubber matrices and small gas molecules due to prohibitive CPU times for diffusion coefficients smaller than 4 × 10⁻⁷ cm² s⁻¹.²⁷  For this reason the technique is not adequate for systematic screening of a large number of polymer candidates and generating data to be used in a materials design approach.

The slowness of MD is an obstacle to the study of polymer properties, but not a barrier. This is especially true when the properties of interest are localized, as in the case of glass transition or diffusivity of small molecules through a polymer matrix. At the far end of what can be achieved, MD can also be applied to study polymer dynamics on the time and length scales of polymer entanglement. For current applications, we focus on methods that have been illustrated with several polymers and are on the verge of being competitive with existing empirical methods.

Predictions of polymer properties underscore the difficulty of developing MD methods that are competitive with existing empirical methods. Simulations of polymers require greater computational effort. In the case of gas penetrants, for example, the equilibrium concentrations are very small. This necessitates simulating a large number of polymer atoms just to simulate a few gas penetrant molecules in order to get a reasonable average. Long simulations of large systems are extremely expensive. On the other hand, academic motivation is very weak for developing MD methods that would be conclusively superior to empirical methods, as demonstrated with similarly large databases. The computational expense for polymer systems exacerbates this problem.

The Monte Carlo technique (MC) is a stochastic simulation method designed to generate a long sequence, or ‘Markov chain’ of configurations that asymptotically sample the probability density of an equilibrium ensemble of statistical mechanics.^107,108  Since its development, MC has been used to test statistical mechanics theories. Today several advances have been reached in devising new statistical mechanical ensembles and designing new MC moves for the efficient sampling of complex configuration spaces. A comprehensive review on progress and outlook in MC simulations has been given by Theodorou.¹⁰⁹ 

MC implemented in the framework of the transition state theory (Section 1.2.1.3) can provide estimates of rate constants for infrequent events. For the construction of polymer membrane models, filling a basic cubic volume element under periodic boundary condition, a rotational isomeric state Monte Carlo technique incorporating long-range interactions¹¹⁰  can be used.

The solubility coefficient can be calculated via simulations in the canonical ensemble in which the chemical potential is calculated using the Widom particle insertion method. The interaction energy of a gas particle inserted within the accessible free volume of the polymer matrix is calculated and the excess thermodynamic potential μ_excess can be estimated from eqn (1.7):

The solubility, S, is then obtained from eqn (1.8):

The simulation of sorption properties of gas molecules in the amorphous cells of glassy polymer can also be estimated using grand canonical Monte Carlo calculations. This needs as input the structural model of the amorphous cell and the force field describing sorbate/sorbent and sorbate/sorbate interactions. For the prediction of gas sorption in the generated amorphous cells, the interaction potential is the most important ingredient. A simplified interaction potential including only a dispersive-repulsive short-range potential, represented by a Lennard-Jones 6–12 potential combined with electrostatic interactions between partial charges on the adsorbent and guest atoms is used. The multipole–multipole interactions are calculated according to:

where A_ij is the repulsion constant and B_ij is the dispersion constant and q_i the point partial charges located at the atomic positions of the adsorbent and sorbate molecules.

The prediction of diffusivity in polymer glasses and low-temperature rubbery polymers via direct MD simulations would require extremely long simulation times because the penetrant diffusion becomes too slow to be predictable by MD.

The TST is a well-established methodology for the calculation of the kinetics of infrequent events in numerous physical systems. According to the TST method the gas transport mechanism through a dense polymer system is described as a series of activated jumps. For each transition, a ‘reaction trajectory’ leading from a local energy minimum to another through a saddle point in the configuration space is tracked, and the transition rate constant is evaluated.

The general approach followed to obtain the diffusivity, based on atomistic TST-based determination of rate constants for individual jumps executed by the penetrant in the polymer matrix and subsequent use of these rate constants within a kinetic MC simulation to track displacement at long times, is another good example of hierarchical modeling.

Three different approaches have been applied for describing the coupling between the jumping of a penetrant and the motions of nearby polymer chains.

In the original Gusev–Suter TST method¹¹¹  a frozen polymer method has been used. All polymer chains are considered fixed in place, and TST-based rate constants are calculated from the energy barriers found for a penetrant to pass from one local potential energy minimum to another. This method is the most straightforward; however, in polymers it yields rate constants that lead to diffusion coefficients much lower (by factors of 10³–10⁶) than experimental values, because neglecting chain fluctuation contributions is physically unrealistic in a polymeric material. Gusev and Suter implemented the frozen polymer method taking into account the thermal vibrations of the polymer matrix^27,112  with the assumption that the polymer atoms, in a sorption site, execute uncorrelated harmonic vibrations around their equilibrium positions to accommodate the guest molecules. The magnitude of the fluctuations is controlled by a parameter similar to the Debye–Waller factor in X-ray scattering. In the most recent method,^20,113  referred as the explicit polymer method, the dimensionality of the jump path was increased to include explicitly both penetrant displacements and polymer chain motions.

Each resulting rate constant then captures the details of the particular chain motions that accompany the jump, rather than mean-field average motions, but at the expense of a much higher computational cost. This method has been applied to a few penetrants in polypropylene,^113,114  with reasonable agreement compared to experiment and correlation.

The TST approach permits the calculation of D by using Einstein's relation (eqn (1.4)) only and the solubility through eqn (1.10):

Computational techniques intended to automate generation and mining of virtual libraries of compounds have successfully been developed and used in recent years for molecular and materials discovery and optimization. This activity has been extensively used in the field of drug discovery where quantitative structure property relationships and quantitative structure activity relationships (QSPR/QSAR) have been used to build correlations between structural molecular features and the experimentally measured properties or activities of the molecules. QSPR/QSAR approaches have been reported in the literature to predict many physicochemical properties, such as vapor pressures,¹¹⁵  aqueous solubility^115–118  and boiling points.^119–121  They have recently been applied also to polymer materials for predicting properties in a number of applications.^122–128  The generation of practical design rules for polymers used in membrane technology would profit from the development of reliable QSAR/QSPR methodologies applicable to this field.

Most of the examples from the literature that use the QSAR approach for predicting gas permeability in polymer membranes use of group contribution methods to establish a correlation between the structure of the repeat unit and some physical property of the polymer membrane such as the free volume, the mean segment distance or dielectric constant polarizability, which in its turn is used to predict permeation properties.¹²⁹  Some practical criteria have emerged to guide synthetic researchers in improving the permselectivity of membranes going through extensive experimentation steps: (i) inhibition of inter-segmental packing whilst simultaneously inhibiting intra-segmental (backbone) mobility; (ii) weakening of inter-chain interactions (reduction of charge transfer complexes). These design rules are based on phenomenological paradigms that provide guidelines for polymer selection correlate the target property with other physical or chemical properties of the polymers, for example using group contribution properties.^130,131  Building a QSPR/QSAR approach that would integrate data at different length scales including processing conditions of the polymer membrane in order to capture information at a multi-scale level would lead to new criteria for polymer materials discovery. The first step in a QSPR/QSAR study is to collect data for a set of polymers that are used for making the membrane. Such data could be either measured experimentally or obtained by numerical simulations. The next step consists in calculating molecular descriptors that mirror fundamental physicochemical factors that in some way relate to gas solubility and diffusion are needed. It is desirable that the description is reversible, so that the model interpretation leads forward to an understanding of how the modification of chemical structure influences gas transport and equilibrium properties. In recent QSPR/QSAR studies^132–134  the following conditions are of fundamental importance for building reliable correlations: (i) the selection of the training set, (ii) the selection of adequate descriptors, (iii) the accuracy of the initial activity-data set used to generate the correlation, (iv) the algorithm used to reduce the number of descriptors to the best set and develop the property model by appropriate regression.

In Figure 1.4 a schematic diagram of the QSAR methodology is reported. The generation of the molecular models for a large number of polymer membranes as well as the collection of reliable simulated values of S and D for several gas molecules in these membranes is the first step towards building a useful QSAR. Generation of appropriate physically meaningful descriptors based on molecular information not only of the repeat unit but also of amorphous cells is key to developing and searching for associations between apparently disparate or disjointed datasets. Generation of different kind of descriptors (topological, geometrical and electronic) is followed by principal component analysis which is a multivariate statistical factor analysis technique. This leads to data reduction since it points to the possibility that only a limited number of properties have to be measured or calculated in order to explain the major part of the information concerning gas permeation through the membranes. Next genetic forms of statistical regression methods such as genetic function approximation or genetic partial least squares are used to generate multiple QSAR models for the collected experimental gas permeation data. The use of this methodology for building reliable correlations between the structure of the polymer, the characteristics of the amorphous cell and the properties of the penetrant molecule with gas transport properties through the polymer membrane is under progress. The development of rapid algorithms that would permit rapid construction and estimation of S and D for different gas molecules paralleled by accurate experimental characterization of well-defined polymer membranes is needed to accelerate this process.

Here we detail through a specific example the results of MD modeling, i.e. the numerical approach used for simulating gas transport properties and the use and limits of detailed atomistic description of polymer membrane models.⁴³ 

In particular we specify on a new protocol for packing glassy polymer with intermediate control of free volume distribution for four polyether-ether-ketones, an amorphous poly(oxa-p-phenylene-3,3-phthalido-p-phenylenxoxa-p-phenylenexoxi-p-phenylene) (PEEK-WC), the dimethyl PEEK-WC (DMPEEK), the tetramethyl PEEK-WC (TMPEEK) and the di-isopropyl dimethyl PEEK-WC (DIDMPEEK).

Modified polyether-ether-ketones^135–138  are of considerable interest due to their excellent mechanical toughness, thermo-oxidative stability, solvent resistance and high transition temperature. In the last decade considerable effort has been spent to introduce chemical modifications in this class of polymers in order to obtain better physical properties and to build up membranes for electro-dialysis, gas dehumidification and gas separation. Relatively few atomistic simulations have been performed on this class of polymers.^37,139  The monomer structure, experimental density and the glass transition temperature of the four poly(ether ether ketone)s is reported in Figure 1.5.

A force field provides a mathematical description of the potential energy as a function of the atomistic configuration of the polymer. The amorphous packing of the PEEK-WC polymer have been constructed using the force field COMPASS (condensed-phase optimized molecular potentials for atomistic simulation studies) force field¹⁴⁰  which is particularly adapted to the simulation of polymeric systems.

In order to minimize chain end effects, each simulation box contains only one single polymeric chain rather than several confined to the same volume, which would lead to increased density of chain ends. The use of single-chain polymers representing bulk amorphous systems is common and has been proven to be quite accurate in replicating the behavior of experimental polymeric systems.^27,139,142  Due to the limited lateral dimensions of the packing models of few nanometers makes it impossible to simulate complete membranes or other polymer-based samples.

For this reason the bulk models are typically cubic volume elements of some nanometers side length and they represent a part cut out of the interior of a polymer membrane. The cubic basic volume element first are filled with segments of a growing chain under periodic boundary conditions following a combination of the Theodorou–Suter¹⁰⁶  chain-generation approach and the Meirovitch's scanning method¹⁴³  reproducing the natural distribution of conformation angles.

Careful consideration of the construction and equilibration procedures is necessary to assure a physically realistic cell. In the case of polymers without aromatic moieties the size of the volume element can be chosen so as to reproduce the experimentally observed or theoretically predicted macroscopic density of the relevant polymer. For partly aromatic polymers, however, the chain packing stage has to be performed at very low densities to avoid ring-catenation and ‘spearing’ effects.

The incorporation of small ‘spacer’ molecules such as methane, carbon dioxide etc. inserted during the initial phase of amorphous cell and removed in subsequent equilibration steps after cell building⁴²  can minimize the effects of ring–ring catenation.

The spacer-free packing models, with a reduced density in comparison to the experimental value, have been subjected to extensive equilibration procedures of NVT-MD (constant particle number, temperature and volume) and annealing simulations combined with force field parameter (torsion, non-bonded and coulomb interactions) scaling steps.

The goodness of each model has been analyzed, before the complete equilibration procedure has been performed at an intermediate stage when the density of the packing boxes was at a value corresponding at the 90% of the experimental value. The variation of surface to accessible volume ratio as well as its gradient with the probe radius has been used to analyze the accessible volume of each cell.

The calculations have been repeated systematically by varying the probe radius in the interval between 1.2 and 2.1 Å, with a 0.1 Å as step.

Several cycles of NPT-MD (constant particle number, temperature and pressure) runs at pressures of thousands of bars the density of the systems has been increased. Besides, simulated annealing runs with temperatures up to 1000 K and NVT dynamics at 303 K were used to further relax the polymer structure. The successive step, besides the check on the quality of boxes at an intermediate stage, has been the reaching of the experimental density by increasing the pressure with several cycles of NPT-MD (constant particle number, temperature and pressure) runs at pressures of thousands of bars.

Moreover, simulated annealing runs with temperatures up to 1000 K and NVT dynamics at 303 K have been used to further relax the polymeric structures.

Equilibration and density adjustment of the polymer system have been achieved through a final 300 ps MD run (Figure 1.6).

The estimation of the diffusion coefficient, D, of the solubility coefficient, S (and consequently of permeability, Pe) of small permeant molecules in a polymeric membrane is strongly dependent on the quality of the amorphous cell used in the calculation. This is particularly true in the case of stiff chain polymers containing aromatic moieties.

The TST approach of Gusev–Suter methodology was used to calculate diffusion and solubility coefficients for N₂, O₂, and CO₂ gas molecules in the generated membrane models.

The comparison of modelled values of modified PEEK-WC with the experimental ones indicates that the new approach of packing produces amorphous cells show less scatter in predicted S and D values with respect to other previous published papers dealing with these membranes.^37,139  The simulated values are more self-consistent showing less scatter than in recently reported papers dealing with modeling of alkylated PEEK-WC membranes.^37,139  A relatively good agreement is obtained for O₂ and N₂ molecules for diffusion coefficient as well as the solubility.

A comparison of predicted to experimental diffusion and solubility selectivities reinforces the conclusions obtained from the analysis of Table 1.1. The predicted diffusion selectivity nicely reproduces the experimentally observed trend for O₂ and N₂ diffusion. The predicted trend for solubility selectivity for O₂ and N₂ is in perfect agreement with the experimental observation with a slightly higher estimation of relative O₂ solubility.

The optimal design of membrane units for gas separation processes is by its nature very complex due to different phenomena taking place at different scales and a real multidisciplinary approach is needed to deal with the different facets that are frontiers of active research in many interconnected fields like materials discovery, membrane fabrication process, optimal fluid flow/distribution inside the membrane module and membrane unit process simulation just to mention some of these. It is to notice that each of these fields of research is by itself complex and implies phenomena that comport a variety of time and length scales. In this contribution we have illustrated through the PEEK-WC polymer membrane simulations the use and limits of detailed atomistic description of gas-membrane system for predicting gas transport properties through polymer membrane models. The need to develop a multi-scale simulation methodology becomes evident even for this case where for estimating the diffusivity of gases through membrane amorphous cell models is difficult to capture using molecular dynamics simulations. The use of TST for estimating gas transport properties of small gas molecules through these glassy polymer membranes has been proven useful for the case studied but it is limited to relatively simple gas permeates that could be treated as single Lennard-Jones sphere without an explicit consideration of the effect of the electrostatic interactions. Recently, improved hierarchical modeling approaches taking into account not only the atomistic description of the gas molecules including the partial charges on each atom but also combine the analysis of accessible volume within model glassy polymer configurations, followed by the identification of transition paths for elementary diffusive jumps in the multidimensional configuration space allowing the polymer matrix to locally relax when the permeate molecule is placed in the transition state and estimate rate constants for these jumps based on TST.^20,113,114 

The field of atomistic modeling of gas transport properties in membranes has experienced a rapid progress due not only to the well known Moore's law on computational power growth but also to the development of clever and efficient algorithms dealing with the simulation of different time-length phenomena that are important for the realistic description of gas permeation through membranes. It is worth noting that even with the recent developments the methodology used for the generation of the realistic models of the models of the membrane and described in this contribution is computationally quite intensive and can be hardly applied for the systematic search of novel polymers with desired properties. Several attempts^144–146  on generalizing simulation protocols based on coarse-graining methodologies have been reported in the literature but to our knowledge the difficult part is to develop methods that automatically would define the coarse grained groups and define their interaction parameters from their composition.

The gas permeation through the polymer membrane as a function of the structural features of the polymer chain and its repeat unit is only one of the aspects of the research of better membranes for a given gas separation. The properties of the membranes are defined also by the process of its formation. This process controls the morphology of the membrane. Most commercial membranes are prepared by immersion precipitation in which the polymer solution is cast on a substrate that could be a flat or hollow fiber and then immersed into a coagulation bath containing the non-solvent (preferably water). The kinetics of diffusion of the water into the polymer solution and of the solvent into the coagulation bath defines the membrane formed by a phase separation process followed the solidification of the polymer rich phase. The performance of the membrane strongly depends on the morphology obtained during the phase separation and solidification. Mesoscale methods have been recently¹⁴⁷  used to simulate this process. The simulation results look very promising in terms of providing a deeper understanding of the immersion precipitation process and serve for the guidance of experiments for designing better membrane formation.

The design of optimal membrane modules for a given gas separation depends also on the fluid distribution through the numerical modeling of the hydrodynamics, mass transfer for single module configurations and especially the pressure distribution in the case of pressure-driven processes. The hydrodynamics allows the increase of the shear stress near the wall, thus allowing the enhancement of membrane separation processes. The commercial available numerical methods for treating fluid distribution and transport phenomena at the level of the module geometry are the finite element and the computational fluid dynamics (CFD) methods.^148–153  CFD methods have been used to optimise membrane separation processes.^148–151  General models of hollow-fiber and spiral-wound membrane modules reported^148,149  were developed from rigorous mass, momentum and energy balances. The CFD analysis of transport phenomena basically assumes the continuity of matter and deals with problems on the length scale significantly higher than molecular dimensions. For the cases where the hypothesis of continuity of the fluid matter cannot be made the adapted approach is based on the lattice-Boltzmann model.^152,153  The fundamental equations of the model stem from the kinetic theory of gases and constitute the basic law employed in the stochastic fluid mechanics. Many engineering problems are concerned with fluid flow and heat transfer in micro-scale systems dealing with length scales of sub-micron order. For these cases the MD simulation method helps to get microscopic insight to the region near the interface and transfer it to the upper level of the lattice-Boltzmann simulation. The influence of long-range intermolecular forces between the porous membrane surface and fluid in the immediate contact with the pores and the membrane on the can be evaluated using microscopic information and contribute towards understanding of the complex physics of fluid flow.

The rational research for better gas separation polymer membrane systems necessitates a comprehensive understanding of phenomena taking place at different time and length scales. Multi-scale simulation is emerging as the adequate technique for the description of such complex systems. During the last decades several computational methods have been successfully developed and used to complement experimental studies of different aspects of membrane science and engineering. In this connection, many traditional simulation techniques (Monte Carlo, molecular dynamics, transition state theory, lattice-Boltzmann, computational fluid dynamics, finite element) have been employed by several groups for obtaining responses for various time and length scales from atomistic, to mesoscale (coarse grains, particle beads) and to macroscale domains (continuum, computational fluid dynamics and finite element). During the last years hierarchical computational approaches for coupling methods used to describe different spatial or temporal domains have been developed. Their use needs long preparation protocols and the simulations are computationally very demanding, thus not so applicable to the prediction of properties of complex systems like membrane modules. The development of methodologies that would allow the automatic definition of each calculation domain and the coupling protocols between each adjacent domain. Despite the progress made over the past years, there are a number of challenges in computer modeling and simulation. In general, these challenges represent the work in two directions. First, there is a need to develop new and improved simulation techniques at individual time and length scales. Second, it is important to integrate the developed methods at wider range of time and length scales, spanning from quantum mechanical domain (a few atoms) to the molecular level. Developing such a multi-scale methodology is very challenging but it will open bright perspectives for numerical simulation, not only in polymer membrane systems but in many other fields concerned by a large span of spatial and temporal domains. New computational tools need be developed in the future to make multi-scale modeling a useful approach for industrial design and predicting properties of materials used in complex systems like the membranes for gas separation.

The authors are indebted to the European Commission for financial support through the PERMOD project (Project Number G5RD-CT-2000-00200) and MULTIMATDESIGN project (Project Number NMP3-CT-2005-013644) and to all partners of both projects for their collaboration. Interactions with Dr Dieter Hofmann and Dr Matthias Heuchel of GKSS, Teltow, Germany, are deeply appreciated.