Identifying pathways to metal–organic framework collapse during solvent activation with molecular simulations

Metal–organic framework (MOF) materials are a vast family of nanoporous solids with potential applications ranging from drug delivery to environmental remediation. Application of MOFs in these scenarios is hindered, however, by difficulties in MOF ‘activation’ after initial synthesis – removal of the synthesis solvent from the pores to make the pore space accessible – often leading to framework collapse if improperly performed. While experimental studies have correlated collapse to specific solvent properties and conditions, the mechanism of activation-collapse is currently unknown. Developing this understanding would enable researchers to create better activation protocols for MOFs, accelerating discovery and process intensification. To achieve this goal, we simulated solvent removal using grand-canonical Monte Carlo and free energy perturbation methods. By framing activation as a fluid desorption problem, we investigated activation processes in the isoreticular metal organic framework (IRMOF) family of MOFs for different solvents. We identified two pathways for solvent activation – the solvent either desorbs uniformly from each individual pore or forms coexisting phases during desorption. These mesophases in turn lead to large capillary stresses within the framework, corroborating experimental hypotheses for the cause of activation-collapse. Finally, we found that the activation energy of solvent removal increased with pore size and connectivity due to the increased stability of solvent mesophases, matching experimental findings. Using these simulations, it is possible to screen MOF activation procedures, enabling rapid identification of ideal solvents and conditions and thus enabling faster development of MOFs for practical applications.

Figure S1 shows the adsorption behaviour of nitrogen in IRMOF-1 at 77 K measured experimentally and simulated from both grand canonical Monte Carlo simulations (labelled unbiased GCMC) and transition matrix Monte Carlo (TMMC).The data from both simulation methods agrees within the error bars thus validating that TMMC can be used to describe adsorption.The shape of the two simulated isotherms is very similar to the experimental isotherm but they are shifted towards lower pressure.We ascribe this to the choice of forcefield used, which has previously been shown to overestimate solid-fluid interaction strengths causing the onset of adsorption at lower pressures than in experiments. 14Despite this, the experimental and computational data largely agree in other aspects (total quantity adsorbed, isotherm shape), demonstrating that the adsorbed fluid behaviour was correctly modelled by the simulations.Accordingly, only adsorption configurations and energies were analysed in the remainder of this study, and no analysis the specific pressure values of fluidphase coexistence was carried out.

Figure S1 -A comparison of nitrogen adosprtion isotherms in IRMOF-1 at 77 K taken from experimental data (black), unbiased GCMC (blue) and TMMC (orange). Error bars are one standard deviation around the mean for unbiased GCMC, and ± 2 kbT for TMMC. Experimental data traken
from reference 15.

Figure S2 -Comparison between experimental and TMMC-predicted bulk surface tension for solvent compounds used in this study. The black line represents parity between experimental and simulated results
, while the shaded are represents 15% uncertainty around the experimental value.

Figure S3 (in external file) -Animation of acetonitrile probability density in IRMOF-1 at densities between 160-200 kgm -3 (30-40 acetonitrile molecules per unit cell).
To calculate sorbate phase nonuniformity and clustering within the framework, we calculated instantaneous interaction energies broken down by framework atom.Specifically, we focused on the carbonyl carbon atoms in the linkers, due to their symmetrically equivalent location within all pore spaces, and proximity to the metal-linker bonds which are critical for holding the network together.The distribution of total interaction energies on carbonyl carbon atoms is shown in Figure S4, broken down by their axial orientation within the unit cell (illustrated in Figure S4a).In all cases, the distribution of interaction energies should be approximately identical for all framework atoms due to the symmetrical nature of the crystal unit cell.This is true even in the case of preferential adsorption into different pore spaces within the MOF, as the distribution of pore spaces also remains symmetrically equivalent along the three crystal axes.

Figure S4 -Breakdown of instantaneous MOF-solvent interaction energies on carbonyl carbons by their orientation within the unit cell. (a) Illustration of the a-and b-axis orientations, along the same projection as in Figure 3 of the main article, (b) histogram of instantaneous energies per carbonyl carbon atom in the case of dichloromethane desorption, (c) histogram of instantaneous energies per carbonyl carbon atom in the case of acetonitrile desorption. Interaction strengths were calculated at sorbate densities of 330-410 kgm -3 and 200-240 kgm -3 respectively (40-50 sorbate molecules per unit cell in each case).
The distributions of solvent-carbonyl atom interaction strengths are shown in Figure S4b for dichloromethane and Figure S4c for acetonitrile.Dichloromethane acted as expected within the framework, with approximately identical interactions for all of the carbonyl carbon atoms within the MOF (Figure S4b).Further, this can be visualised in Figures 3, S3, and S5 -sorbate density is periodically distributed throughout the pore space, leading to uniform adsorption stresses in all three crystal planes regardless of dichloromethane's preferential adsorption at sites close to the metal SBUs.
Conversely, histograms of acetonitrile interaction strength show a large difference depending on the linker's axial orientation, indicating heterogeneous phase behaviour within the MOF's pores.This behaviour is clearly asymmetrical, therefore cannot be explained by preferential filling of the larger pore spaces over the smaller (as that would lead to checkerboard-like solvent phases, which remain symmetrical).

Figure S6 -
Figure S6 -Comparison of fluid-fluid and fluid-MOF interactions strengths for (blue) acetonitrile and (orange) dichloromethane in IRMOF-1, normalised by number of adsorbed molecules.