Unravelling the pore network and gas dynamics in highly adaptive rubbery organic frameworks

Rubbery organic frameworks-ROFs have recently emerged as an intriguing class of dynamers by virtue of reversible connections between their building units. Their highly adaptative features at the origin of their spectacular self-healing properties made them also attractive candidates for the development of gas-selective membranes combining high selectivity and fast permeability. So far, little is known on the origin of this unique trait and this clearly hampers the exploitation of this class of dynamers in many areas where stimuli-responsive pore dynamics is of great importance. To address this lack of fundamental knowledge, herein we unravel the self-assembly process of ROFs via the development of an advanced computational methodology combining quantum and force field molecular simulations that enable the description of reversible connections of building units and the long-range organization of the cross-linked ROF network. We demonstrate that both accurate energy barriers associated with the covalent bond formation between the building units and presence of solvent are key parameters to ensure the in silico construction of reliable ROF structure models that are supported by a set of experimental data collected on synthesized ROFs including density, connectivity and porosity. Atomistic insights into the unusual guest-responsive pore dynamics of this intriguing class of dynamers are further gained with a special attention paid to the tunability of this pore flexibility by controlling the chemical composition of the building units. As a further stage, the dynamics of CO2 in these compliance frameworks is scrutinized to shed light on the mechanism at the origin of their promising performance as CO2-selective membranes. We highlight that guest-triggered pore dynamics enables the creation of a diffusion pathway to ensure effective gas transport throughout the whole ROF. This knowledge of the pore structure and its guest-responsive dynamics at the microscopic level is unprecedented in the field of dynamers and it is expected to pave the way towards the optimization of this class of adaptive porous frameworks for many potential applications. Interestingly, this computational approach can be transferable to the exploration of any complex disordered systems showing a high degree of flexibility and guest induced structure/pore reorganization.


Introduction
Membrane technology is extensively used at the industrial level for a myriad of separation applications in liquid/gas phases as it offers many advantages, e.g. energy saving, simple design and easy scale-up. 1,2 Membranes are mostly constructed from polymers 3-9 while more recently various highly selective microporous materials including zeolites and metal organic frameworks MOFs have been processed as pure or composite membranes in association with glassy polymers to potentially address the threshold between permeability and selectivity for a range of gas mixtures. [10][11][12][13][14][15] Alternatively, rubbery organic frameworks (ROFs) constructed by reversibly connecting components under molecular control via a reversible chemical bond, can address the poor gas selectivity exhibited by classical rubbery polymeric membranes with a non-reversible covalent bond network, while maintaining high gas permeabilities. [16][17][18] The high exibility and mechanical compliance of this intriguing family of dynamic polymers named dynamers, confers excellent mechanical stability and both intrinsic and guest-responsive pore dynamics to the corresponding ROF membranes. [16][17][18] These key features open new avenues towards the fabrication of adaptive membranes with high processability and optimal gas transport performances. Furthermore, the high chemical variability of the assembled building blocks of different shapes/sizes and chemical features enables custom-made design of almost innite number of architectures, making this family of ROFs highly tunable for target gas separation. [19][20][21][22] However, ROFs are poorly crystalline and oen metastable most likely due to continuous reversible reorganization of their bond networks that is also at the origin All publication charges for this article have been paid for by the Royal Society of Chemistry of their unique self-healing or self-repair properties. 18 The complexity and evolutionary nature of this family of dynamers explain that not only their formation mechanism but also their structure and guest-responsive structural changes are still unknown to date. Atomistic insight into the ROF structure organization is therefore required for a ne-tuning of this class of dynamers with the optimum pore network for target application. This challenging objective calls for the deployment of an innovative modelling strategy integrating (i) quantum calculations to accurately assess the energy barriers associated with the reactive assembly of building blocks and (ii) a hybrid Monte Carlo/molecular dynamics scheme that accounts for the reversible long-range order assembly of building blocks to deliver a reliable atomistic structure model of the dynamers. Herein, we devised such a computational strategy (workow schematized in Fig. 1) that was applied to a prototypical dynameric ROF composed of a isophthalaldehyde core building block and two amine-terminated connectors, stick-shaped polytetrahydrofuran (polyTHF) and star-shaped glyceroltris(polypropylene glycol)ether (polyMePEG), which we previously demonstrated to act as an effective membrane for the selective CO 2 capture. 16 We unravel for the rst time a ROF structure model showcased with the self-assembly of the building blocks mentioned above in the presence of solvent by explicitly considering the reversibility of the HC]N imine bond formation/dissociation, which connects the N atom of the amine group of polyTHF or of polyMePEG to the C atom of the aldehyde function belonging to isophthalaldehyde. The structure models built in silico were carefully analyzed in terms of local and long-range arrangements of the building units, connectivity degree and pore distribution/dimension. A set of experimental data collected on this ROF including the fractional free volume, density and assessment of the connectivity served as a preliminary validation of the structure models.
As a further stage, these structure models were used to predict the CO 2 adsorption and transport properties of the ROFs. These Monte Carlo-Molecular dynamics simulations evidenced that these dynamers are highly adaptive and capable of adsorbing gas in evolutive emerging porosity owing to their reversible covalent bond network. Remarkably, our calculations revealed how the diffusion of CO 2 can be modulated by changing the composition of the ROF in terms of polyTHF and polyMePEG concentrations, paving the way towards the selection of the best ROF composition to ensure optimum CO 2 transport throughout the dynamers.
2 In silico construction of the atomistic ROF structure model The rst explored ROF system was made of 25% polyMePEG (M n $ 3000 g mol À1 ) and 75% polyTHF (M n $ 1100 g mol À1 ) connected via isophthalaldehyde centers (dialdehyde). Therefore, this dynamer denoted hereaer as S25, was modelled with the consideration of 150 diamine-terminated polyTHF building units and 50 triamine-terminated polyMePEG building units, corresponding to a total of 450 reactive sites (2 and 3 terminal functions for each polyTHF and polyMePEG respectively) distributed randomly with 225 isophthalaldehyde connectors in a simulation box of 50 Â 50 Â 50Å 3 dimension. Chloroform was considered as the standard solvent previously used for the synthesis of this ROF. This solvent was further removed by drying the membranes prior to testing their CO 2 capture performance. The assembly of the dynamer was performed with a MC scheme to account for the imine bond formation/ dissociation resulting from the reaction between the carbonyl moiety of isophthalaldehyde and the amine group of either polyTHF or polyMePEG (eqn (1)).

R-HC]O + H
The MC algorithm rst lists the pair of reactive sites along the MD trajectories, identied for new HC]N imine bond Fig. 1 Overall workflow of the computational strategy to gain insight into the assembly of ROF dynamers and the resulting structure in the presence of solvent. The selected ROF/solvent system is composed of a connection center (isophthalaldehyde) and two amine-terminated connectors (stick-shaped polyTHF and star-shaped polyMePEG) with chloroform considered as a model solvent currently used for the synthesis of this ROF. (i) The free energy surface for the imine reaction is first determined using quantum calculations coupled with metadynamics. This enables the determination of the energy barrier for the assembly of the building units that is further implemented into a hybrid (ii) MC and (iii) MD scheme to deliver an atomistic ROF structure model further characterized in terms of porosity and connectivity prior to exploring (iv) its gas adsorption and transport properties.
formation using a C-N distance criterion of 2.8Å. The HC]N imine bond formation implies the release of one water molecule as shown in eqn (1). A cut-off of 2.5Å for the distance between the C atom of the HC]N imine bond and the O atom of H 2 O is subsequently considered for the bond dissociation. The Metropolis acceptance criterion is then applied for the bond formation/dissociation with the following probability P ¼ min(1, e ÀbDE ), where DE is the energy barrier for either bond formation (DE f ) or bond dissociation (DE d ). The consideration of these energy barriers enables sampling of the equilibrium distribution of the dynamers. This critically calls for an accurate assessment of these corresponding energy barriers. To this purpose, quantum calculations were preliminarily conducted to determine the free energy proles for the HC]N imine bond formation. All calculations were performed at the density functional theory (DFT) level using the PBE functional as implemented in the CPMD code 23 and a cut-off energy of 90 Ry. A metadynamics scheme with gaussian-shaped biasing forces (width of 0.1 u.s. and height of 0.01 eV, frequency every 10 steps 24,25 ) was selected to effectively scan the free energy prole of the reaction. In this case, two collective variables were considered: (i) the intermolecular C-N distance to enhance the imine bond formation and (ii) the intramolecular C-O distance in the dialdehyde to mimic the release of the O atom in the form of H 2 O during the reaction. Fig. 2 shows the calculated free energy prole of the HC]N imine bond formation for the 1 isophthalaldehyde/1 polyTHF pair introduced in a simulation box of 20 Â 20 Â 20Å 3 dimension. This prole exhibits three minima: one for the non-bonded polyTHF and dialdehyde building units, one for each assembled building unit implying a hemiacetal HOCH-NH bond or the HC]N imine bond just before and aer releasing water respectively (see eqn (1)). Along this reaction path, the maximum energy barrier is 0.7 eV and 0.9 eV for the HC]N imine bond formation (DE f ) and dissociation (DE d ) respectively. Note that the CPMD calculations were performed without the inclusion of dispersion corrections. Since the energy barrier is rather high, the long-range dispersion energy contribution is expected to be negligible and therefore does not affect the energy difference between the energy barrier for the bond formation and bond dissociation. The free energy prole was equally calculated for the 1 isophthalaldehyde/1 polyMePEG and 1 isophthalaldehyde/2 polyTHF pairs leading to similar energy barriers (see Fig. S1 in the ESI † showing the energy proles).
Once the bonds are equilibrated, the connectivity of the dynamer is modied accordingly and the structure is further relaxed using subsequent MD steps. All MD simulations are performed with the LAMMPS package 26 using the GAFF force eld 27 for all atoms of the system treated as charged Lennard-Jones (LJ) sites. The atomic charges of the initial building units and chloroform molecules were extracted by the restrained electrostatic potential (RESP) 28 charge tting approach applied to the electrostatic potential calculated in Gaussian 29 using the HF/6-31G* basis set. 27 Moreover, when the components are connected, the charges of C and N atoms in the HC]N imine bond are modied using values calculated with the same RESP method computed on a system containing pol-yTHF bonded with isophthalaldehyde and a free water molecule (see ESI † for the values of charges in the bonded and nonbonded cases). The non-bonded contribution was treated as the sum of the LJ term with a cutoff of 10Å and an electrostatic interaction calculated by means of the Ewald summation method as implemented in LAMMPS with a k-space value of 0.0001. Bonded-contribution includes harmonic potentials to model the stretching and bending modes, and cosine-based functions for dihedral and improper dihedral torsions (all the parameters are from the GAFF force eld 27 ). These MD simulations were performed in the NPT ensemble with a time step of 0.5 fs and using a Nose-Hoover barostat and thermostat with a relaxation time of 0.5 ps and 0.05 ps respectively. The building units are rst equilibrated during 1 ns prior to starting the MC scheme and the MD simulations are further run for at least 50 ns. For every subsequent 1 ps MD simulation, one MC step is considered corresponding to 10 attempts to realize bond formation/dissociation according to the Metropolis criterion mentioned above (see movie of the bond association between polyTHF and isophthalaldehyde in the ESI †). For each attempt, the energy barrier of the HC]N bond formation/dissociation is considered to determine the probability of the reaction. 3 Analysis of the structural and textural features of the ROF model Fig. 3 shows the MD simulated time-dependence of the connectivity for the ROF S25 calculated as the fraction of carbonyl groups of isophthalaldehyde involved in the formation of a bond with polyTHF or polyMePEG. The connectivity converges to about 62% aer 50 ns MD simulations (see the black line in Fig. 3) with isophthalaldehyde slightly more connected to polyMePEG (about 65%, green line) compared to polyTHF (about 57%, blue line). This evidences that the isophthalaldehyde units reacted at least once via their reactive carbonyl function, consistent with previous experimental ndings. 16 The density of the so-constructed model is 1.10 g m À3 , which is in agreement with the experimental data previously reported for the synthesized ROF S25. 16 We further assessed how the energy barrier for the bond formation/dissociation affects the connectivity of the structure model for ROF S25. The test case implies barrier-less bond formation and an innite barrier for the bond dissociation, a typical scenario considered in most of the in silico polymerisation soware (Fig. 3, red line). This leads to a very high connectivity (95%) reached aer only 10 ns MD simulations. This resulting phase is so dense that the connectivity cannot be converted back to that of the initial simulated structure by again applying the energy barrier calculated at the DFT level (Fig. 3, black line). This observation emphasizes that an accurate evaluation of the energy barrier for the HC]N bond formation/dissociation is a key to achieve a reliable atomistic model for such dynamic ROF systems.
The resulting atomistic model of ROF S25 is highly disordered as illustrated in Fig. 4a with intertwined chains that still contain non-bonded amine and aldehyde reactive sites consistent with incomplete connectivity as reported in Fig. 3. This is reected in the radial distribution function (RDF) calculated for the C-N pair (Fig. 4b) that exhibits the rst peak at 1.5Å assigned to the intramolecular HC]N bond length and the second one at 2.5Å accompanied by additional contributions above 3Å corresponding to non-bonded distances between unreacted amine and dialdehyde groups. Fig. 4a also shows that the connections of the dialdehyde building units can be achieved via one polyTHF and one polyMePEG molecules. The length of the polyTHF chain in the ROF S25 (about 18Å) is reduced by about 10Å in average compared to its length exhibited in the gas phase (about 28Å), due to the angular coiling of the molecule in the ROF (Fig. 4c). This strongly suggests that the polyTHF building units are able to re-extend if an external stimulus is applied to the dynamer leading to a high exibility of this family of so solids. On the opposite, the star-shaped polyMePEG building units (see Fig. 4d), are only slightly shrunk (d avg ¼ 8.2Å) compared to their conformation in the gas phase (d avg ¼ 8.5Å).
The resulting ROF S25 structure model shows a fractional free volume of 16% (Fig. 5a -conguration (i)) calculated with a probe molecule of 3.3Å, corresponding to the pore regions accessible to CO 2 considered as a model molecule diffusing in  the corresponding membrane. 16 This fractional free volume is in line with the corresponding experimental data reported previously for this ROF (about 11-12%). 16 Analysis of the pore size distribution (PSD) reported in Fig. 5b shows that the pores of ROF S25 exhibit a predominant contribution at 3.8Å.
When the chloroform is released from the ROF S25 structure model, the MD simulations converge towards a denser phase as revealed by a much lower fractional free volume (6%) (Fig. 5a -conguration (ii)). This reduction of porosity is associated with a substantial decrease of the population of pores with a dimension of 3.8Å that are shied to small voids below 3Å as illustrated in the corresponding PSD plot (Fig. 5b). Remarkably, the re-incorporation of the solvent in this guest-free model leads to re-opening of the porosity as observed in Fig. 5a -conguration (iii), with an associated fractional pore volume (15%) and PSD prole (Fig. 5b) that t with the corresponding data obtained for the model constructed initially in the presence of chloroform. These overall computational ndings deliver an unprecedented atomistic insight into the responsive pore dynamics of the ROFs upon guest adsorption/desorption. This prediction is also in line with the unique self-healing properties of this family of dynamers previously revealed. 18 Note that experimentally the self-healing mechanism occurs within minutes aer two cut pieces are held together.
We equally envisaged the construction of an atomistic ROF structure model without the initial presence of solvent. While the resulting solvent-free structure (Fig. 5a -conguration (iv)) shows connectivity very similar to that of the structure model built with chloroform (see Fig. S2 in ESI †), it is almost nonporous (fractional free volume lower than 1%) and much denser than the structure obtained aer solvent removal (Fig. 5a -conguration (ii)). This conclusion emphasizes that the solvent is a key parameter to account for the modelling of highly exible dynamers. It is to be noted that the consideration of solvent is only rarely considered in the modelling of structure models for standard polymeric systems and therefore can be a severe drawback for the most complex polymers implying a high degree of exibility. In our specic case, chloroform is expected to facilitate the mobility of the building units in the simulation box in order to increase the probability of the reactive groups to be close to each other for initiating the HC]N bond formation. This is evidenced by a faster convergence of the constructed structure model in the presence of solvent (see connectivity in Fig. S2 †). In addition, chloroform as a non-polar solvent interacts with the organic part of the chains and directs the structural organization of the dynamer that affect its pore distribution and hence its resulting fractional pore volume.

Prediction of the CO 2 transport in the ROF
As a further stage, we explored the CO 2 transport in this generated ROF S25 structure model to gain unprecedented microscopic insight into the attractiveness of this family of dynamer membranes for CO 2 separation. First, Monte Carlo simulations using the CADSS soware 30 were performed to load CO 2 molecules on the ROF S25 structure obtained aer removing chloroform and further MD simulations were performed for 50 ns at room temperature to explore the guest diffusion (100 CO 2 molecules were adsorbed at 3 bar considered as a typical pressure). The simulated time-dependent mean square displacement (MSD) prole of the CO 2 molecule, shown in Fig. 6 exhibits initially a ballistic regime and further adopts a linear trend characteristic of the Fickian-diffusion regime. The Fig. 5 (a) Evolution of the fractional free volume of the ROF S25 depending on the conditions used in the MD simulations for the generation of the structural models. The section in green corresponds to a cycle starting with the structure constructed with the inclusion of the solvent (configuration (i)) followed by the subsequent removal (configuration (ii)) and reincorporation (configuration (iii)) of the solvent. The pores accessible to a probe molecule of 3.3Å (kinetic diameter of CO 2 ) are represented in pink. The non-porous structure generated without considering chloroform initially is illustrated in gold. (b) Pore size distribution (solid lines) and contribution to the total fractional free volume pore density (dashed lines) for the structures (configuration (i), black), (configuration (ii), red) and (configuration (iii), blue).

Fig. 6
Mean square displacement (MSD) calculated for CO 2 in ROF S25 (black) and ROF S75 (red) over 100 ns MD simulation runs at 300 K. The snapshots on the right show the diffusion path followed by one tagged CO 2 molecule in each ROF. The self-diffusion coefficient (D s ) of CO 2 in ROF S25 is calculated using the Einstein relation applied in the 8-10 ns range.
self-diffusion coefficient (D s ) of CO 2 calculated using the Einstein relationship was found to be 5 Â 10 À11 m 2 s À1 at 300 K (see Fig. 6) in excellent agreement with the previous experimental work on the corresponding membrane that estimated a diffusion coefficient of about 10 À11 m 2 s À1 . 16 For comparison we investigated the diffusion of CO 2 in a ROF structure model denoted as ROF S75 containing 75% of polyMePEG corresponding to 50 polyTHF building units, 150 polyMePEG building units and 275 dialdehyde building units introduced in the simulation box of the same size as that of ROF S25 (50 Â 50 Â 50Å 3 ). This ROF S75 shows a connectivity of 70%, which is slightly higher than that of ROF S25 (65%) while the fractional free volume of 12% is a bit lower (vs. 16% for ROF S25), consistent with a slight increase of the density from 1.08 to 1.10 g m À3 . In order to have a direct comparison with ROF S25, the same number of CO 2 molecules was added to the ROF S75 model. Fig. 6 shows that the mobility of CO 2 in this ROF S75 is rather limited and did not enable us to extract a self-diffusion coefficient. This prediction again is in excellent agreement with the experimental observation previously reported showing that a higher concentration of polyMePEG incorporated in the ROF membrane leads to a drop of the CO 2 diffusion. Interestingly the ratio of the MSD values for CO 2 between ROF S25 and ROF S75 aer 100 ns is about 4 : 1, in qualitative agreement with the experimental measurements revealing a ratio of permeability of 5 : 1 for the two corresponding membranes. 16 We further explored the microscopic origin of the distinct gas transport behavior exhibited by ROF S25 and ROF S75. Aer adsorption, the pore size distributions of ROF S25 and ROF S75 (see Fig. S4 †) and their free volumes (10% and 9% respectively) are similar. However, ROF S25 shows a higher degree of exibility as reected by its lower calculated bulk modulus (2.4 GPa) as compared to ROF S75 (4.9 GPa) (see ESI † for details of the calculations). This trend is consistent with a larger concentration of the exible polyTHF molecule present in ROF S25. Therefore the higher CO 2 permeability observed for ROF S25 stems from its higher exibility that enables to generate CO 2triggered adaptive paths for optimum CO 2 transport throughout the whole system.

Conclusion
Herein, we unraveled the rst structural model for the family of dynameric rubbery organic frameworks via the development of a computational methodology integrating quantum-and force eld molecular simulations to account for the reactivity of the assembly of building units and the long-range organization of the cross-linked network. The in silico constructed atomistic models were validated with experimental data collected on the synthesized ROFs in terms of porosity, connectivity and density. We demonstrated that the consideration of an accurate energy barrier associated with the chemical assembly of the building units via covalent bonding and the presence of solvent, both features most oen neglected in the standard assembly soware, is a key to achieve reliable structure models for this family of highly adaptive frameworks. Remarkably, their pore networks were revealed to be reversibly guest-responsive featuring preferential pathways for the guest molecules to diffuse at the origin of the very attractive CO 2 permeability performance of the corresponding membranes. This exible behaviour was shown to depend on the chemical composition of the frameworks controlled by the concentration of the different building units considered initially for the assembly process. This knowledge is expected to pave the way towards the optimization of the pore structure/dimension of this class of adaptive porous frameworks to tune their gas transport properties and guide the fabrication of membranes with improved performance for the separation of molecules of industrial importance. Beyond this class of dynamers, the advanced computational methodology developed in this study is transferable to the exploration of any complex disordered systems showing a high degree of exibility and guest-induced structure/pore reorganization.

Data availability
All the computational details and simulated extra data are provided in the ESI. †

Conflicts of interest
There are no conicts to declare.