Over the last decades, metal-organic frameworks (MOFs) have emerged as an extraordinary class of materials as a result of their unique building block concept from metal nodes and organic linkers, giving rise to an almost limitless number of materials (Chui et al., 1999; Li et al., 1999; Férey, 2001; Yaghi et al., 2003; Kitagawa et al., 2004). MOFs may display anomalous responses to external triggers and exhibit for example negative linear compressibility, negative thermal expansion, negative gas adsorption, or large-amplitude structural transformations while retaining their structural integrity (Coudert, 2015; Balestra et al., 2016; Bennett et al., 2017; Liu et al., 2018; Burtch et al., 2019; Coudert and Evans, 2019; Evans et al., 2019). Such stimuli-responsive behavior has been explored in various research fields including adsorption, separation, catalysis, and drug delivery (Horike et al., 2009; Li et al., 2009; Horcajada et al., 2012; Schneemann et al., 2014). More in particular, the terminology “soft porous crystals” (SPCs) was coined by Kitagawa and colleagues referring to those materials that exhibit bistable or even multistable behavior, possessing the ability to undergo large structural transformations upon exposure to external stimuli while maintaining their long-range structural order (Horike et al., 2009). Enormous experimental and computational research efforts have been undertaken to understand and ultimately predict or tune this functional behavior. However, most of our insights so far originate from thermodynamic considerations, while the mechanistic details of such large-amplitude phase transformations are yet to be resolved. For example, it is still unclear how the transformation nucleates and how this nucleation is affected by the presence of defects or the size of the crystal. Experimentally, a series of in situ experiments have established that the dynamic response of these materials is strongly affected by the presence of defects and the crystal size, whereby crystal downsizing was found to suppress their ability to morph (Sakata et al., 2013; Miura et al., 2017; Krause et al., 2018; Ehrling et al., 2019; Kundu et al., 2019; Wannapaiboon et al., 2019; Bon et al., 2020; Krause et al., 2020; Ehrling et al., 2021).

Another major point of discussion regarding the transition mechanism concerns the degree to which the transition occurs in a collective, cooperative way. In many literature references, the hypothesis of collective behavior was assumed, where the entire framework transforms cooperatively (Schneemann et al., 2014; Wieme et al., 2018; Vanduyfhuys et al., 2018). Insights into the thermodynamics of the phase transformations have been obtained from atomistic simulations by constructing the Helmholtz free energy as function of the state variables governing the observed behavior (Rogge et al., 2015; Vanduyfhuys et al., 2018; Hobday and Kieslich, 2021). Figure 1B shows a series of conceptual free energy curves in terms of a collective variable that is able to uniquely distinguish between all (meta)stable or activated states. For many frameworks the unit cell volume proved to be an appropriate collective variable (Demuynck et al., 2018), however, this may not always be the case. For flexible linkers, more advanced procedures were necessary to find collective variables that allowed to describe the phase transformation. The study of such frameworks is beyond the scope of the current contribution. Herein, we will focus on MIL-53(Al), the prototypical flexible framework where transitions may be induced by temperature, pressure and gas adsorption (Demuynck et al., 2017, Demuynck et al., 2018). Experimentally determined structural parameters are presented in Supplementary Section S1. From the perspective of these Helmholtz free energy curves, collective behavior implies that any small barrier in the thermodynamic potential of a single unit cell would translate into a huge barrier for the entire system and the framework would only transition in case one of the minima in the thermodynamic potential would disappear (curves 1 and 3 in Figure 1B).

Nowadays, there are various indications that transitions in SPCs do not occur in a fully cooperative way. Early on, Triguero et al. proposed a theoretical layer-by-layer transition model for winerack frameworks like MIL-53(Al), where collective transformations between the large pore (lp) and closed pore (cp) phases only occur within a single layer at a time (Triguero et al., 2011). Their findings were based on a model that exploits the characteristic geometry of winerack frameworks, which are composed of rhombus-shaped cells in the directions perpendicular to the aluminum oxide chain. The shape of the rhombus cells can be characterized by the angle θ, which varies between 79° for the lp phase and 40° for the cp phase (see Figure 1C). The edge lengths do not change appreciably because the corresponding benzene-1,4-dicarboxylate (BDC) linkers do not allow for large variations in length. These specific geometrical features impose significant constraints on the possible deformation modes of the crystal if its crystallinity is to be preserved; as only those deformations are allowed which preserve the rhombus shape. With these assumptions, cells in a layer must deform coherently in order to preserve the lattice integrity. As such, it was concluded that the phase transformation occurs in an avalanche manner; once the lp-to-cp nucleation occurs at one point, the entire layer quickly transforms to the cp phase. Further lp-to-cp transformations were predicted to occur most likely in layers neighboring already transformed layers. More recent work has constructed energy versus volume curves for MIL-53(Al) using very accurate many-body dispersion calculations within the random-phase approximation (Wieme et al., 2018). Even at 0 K, a non-negligible barrier of about 8 kJ/mol was found between the lp and cp phase, indicating the existence of metastable lp crystals even at very low temperatures. This conclusion is corroborated by Mendt et al. who observed a fraction of lp material for MIL-53(Al) even at 9 K (Mendt et al., 2010).

Computational modeling can help to unravel such mechanistic details, provided that the employed models are realistic representations of experimentally observed MOFs, with dimensions comparable to real crystallites and with the explicit inclusion of the external surface and potential defects. Unfortunately, even with the massive amount of computing resources that is now available, such simulations are not yet feasible (Van Speybroeck et al., 2021). At present, the attainable length scales within the field of nanostructured materials are limited to a few nanometers and common molecular dynamics (MD) runs extend from the picoscale (for ab initio based methods) to the nanoscale range (for force field based methods). In the MOF field, the largest force field based simulations were performed on systems of about 10,000 atoms, with a length scale of tens of nanometers in one direction. Such spatial dimensions are still below experimental crystal sizes which may extend well into the micrometer range (see Figure 1A). As a consequence of these model limitations, some effects which are inherently related to (long-range) spatial disorder are not properly accounted for in current computational models (Van Speybroeck et al., 2021).

Recent work has attempted to simulate phase transformations in SPCs in a more realistic manner (see Figure 1A). Schmid et al. performed the first finite size simulations on nanocrystallites of DMOF-1(Zn) by setting up a crystallite model containing roughly 250,000 atoms, with explicit inclusion of the external surface (Keupp and Schmid, 2019). Such simulations abandon the periodic boundary conditions (PBCs) that are commonly employed in MOF computational modeling. A series of simulations were performed in a temperature ramp between 300 K and 500 K to observe the thermal opening, and constrained simulations were performed to observe the mechanical closing of the DMOF-1 system. They concluded that the transition nucleates near the external surface, where an interface between the cp and lp phase occurs which then travels dynamically through the lattice. The nucleation was found to be distinctively dependent on the surface-to-volume ratio. Such observations would not have been possible based on simulations with PBCs and small-sized unit cells, as spatial disorder is in that case not allowed and the external surface not taken into account. Simultaneously, Rogge et al. performed large scale force-field-based MD simulations on a series of MOFs, including MIL-53(Al), DMOF-1(Zn) and CoBDP using PBCs, with supercells that contained over 10,000 atoms (Rogge et al., 2019). The latter supercells were referred to as mesocells and their behavior was contrasted with nanocells having substantially smaller dimensions (see Figure 1A). Although the mesocells are still about an order of magnitude smaller than experimental crystal sizes, new physicochemical phenomena were observed. For the nanocell, the lp-to-cp transition occurs cooperatively and any form of spatial disorder is prevented by the enforced PBCs. On the other hand, when mesocells are exposed to pressures of 40 MPa at 300 K, part of the system undergoes an lp-to-cp transition whereas the rest temporarily remains in its original lp phase. Temporary interfacial defects were observed near the lp-cp phase boundary, which traverse the lattice until the transformation is complete. The phenomenon was also thermodynamically investigated by constructing the Helmholtz free energy profiles as a function of the volume, which strongly depend on the possibility of phase coexistence occurring in the larger mesocells. The latter phenomenon has not yet been confirmed experimentally, as this requires dedicated in situ cells where the temporal behavior of the material could be followed with very high spatiotemporal resolution.

From this discussion, it is clear that a full mechanistic understanding of the phase transformations in finite MOF crystals requires models that are substantially larger than the ones considered so far, as these would allow us to include spatial disorder at various length scales. Ultimately, simulations of finite-sized crystals containing millions of atoms and having sizes comparable to small experimental crystals would become feasible. In this work, we take an important step in this direction by exploiting the massively parallel architecture of state-of-the-art GPUs in order to simulate these phase transitions at unprecedented length scales. GPU acceleration is known to increase the attainable time and length scale of MD simulations by roughly two orders of magnitude as compared to a CPU, which is achieved by offloading the evaluation of interatomic forces and possibly the time integration onto the GPU. While most molecular mechanics engines nowadays provide some level of GPU acceleration, they are rarely compatible with the specialized force fields and thermodynamic ensembles that are required in common SPC simulation pipelines. In this work, we succeeded in porting the necessary force fields and sampling protocols to the highly extensible and GPU-oriented OpenMM software package (Eastman et al., 2017). In particular, this included the development and implementation of an extension of the existing pressure control algorithm towards fully anisotropic unit cell fluctuations. The proposed workflow enabled us to go far beyond currently accessible atomistic force field based MD simulations on MOFs. Herein, we perform simulations on models containing 1,040,440 atoms and a unit cell size of 54.9 × 6.6 × 45.1 nm³, thereby entering the range of small experimental crystals (Figure 1A). Visualization and analysis of the transition dynamics become increasingly challenging for such large systems. Herein, we present a dynamic two-dimensional lattice representation that provides a fundamentally new perspective of the framework dynamics. Based on the presented workflow, we are able to give more insight into operative transition mechanisms for different values of the external stimulus.

The remainder of the paper is organized as follows: Section 2 presents the methodological advances that were necessary to achieve this goal, with in particular the derivation and implementation of the proposed pressure control algorithm and the automated dynamic lattice representation. Section 3 then discusses the analysis of the observed transition mechanisms; the concluding discussion and perspectives for future work are provided in Section 4.

The phase transition mechanism is investigated using MD simulations on a massive scale, with framework models containing up to a million atoms and thousands of pores, and over timescales on the order of 10 ns. Such simulation sizes are unprecedented for SPCs, and require state-of-the-art GPUs in order to become computationally feasible. While GPU-accelerated MD has been around for over a decade, accurate computational models of the large-amplitude transition in winerack-type frameworks additionally require (i) an efficient implementation of a number of specialized force field interactions that are not commonly found elsewhere, and (ii) full flexibility of the simulation cell, as the box vector lengths and angles vary significantly during the transition (Rogge et al., 2018). In this work, we chose to employ the OpenMM software package because of its highly versatile and extensible API, and its flexibility in supporting custom interactions, thereby satisfying the first requirement. The second requirement—fully flexible unit cells—is more complicated because OpenMM does not natively support fully anisotropic unit cell fluctuations. Moreover, the implementation of existing anisotropic barostats in OpenMM would require an extensive rewrite of the entire codebase because they require the virial stress at every step in order to operate properly. To avoid this, we chose to extend an existing isotropic pressure control algorithm with fully anisotropic fluctuations, because it is relatively easy to implement on GPUs and is compatible with custom interaction potentials. The theoretical framework of the proposed method as well as a series of validation experiments are presented in Section 2.1.

The underlying physical mechanisms of the transition are exposed by visualizing trajectories using a dynamic two-dimensional lattice representation. This provides a fundamentally new perspective on the transition dynamics and the corresponding induced disorder. An overview of the automated and GPU-accelerated visualization procedure is given in Section 2.2.

We investigated the transition mechanism based on large-scale MD simulations at constant temperature and pressure. In order to eliminate PBC artefacts as much as possible and ensure that phase separation and/or coexistence is maximally allowed by the model, we performed the simulations on a 37 × 10 × 37 unit cell containing 1,040,440 atoms and just over 2,500 individual pores. Such unit cell sizes are unprecedented in computational research on MOFs and SPCs in particular, and represent a significant step forward with respect to the state-of-the-art (as visualized in Figure 1A). To identify the pertinent features in the framework dynamics, we simulated the lp-to-cp transition for this system at different thermodynamic conditions. An overview of all simulations performed with their specific control variables is given in Supplementary Table S2. Specifically, we considered three different temperatures (200, 300 and 500 K) and two different pressures for each temperature (100 and 300 MPa). While both pressures are much higher than the experimentally observed transition pressure for this material—which is estimated at about 13–18 MPa at room temperature (Yot et al., 2014)—these pressures ensured that the phase transition proceeds sufficiently fast as to make the simulations computationally feasible on single-GPU systems. Alternatively, enhanced sampling techniques may be used to speed up the transition, for example by biasing the dynamics of the framework along the θ angle indicated in Figure 3. As this is a first case study of the large-scale dynamics of winerack-type frameworks, we chose not to pursue this direction and consider regular unbiased dynamics at elevated transition pressures instead. Nevertheless, the transition times of the employed 37 × 10 × 37 system are observed to be roughly three orders of magnitude larger as compared to a fully cooperative transition in a 1×2×1 cell (∼10 ns versus ∼10 ps).

Figure 4 shows snapshots of the transition mechanisms at 300 K in the low- and high-pressure regimes, as obtained from simulations at 100 MPa (left) and 300 MPa (right). The final aim is to unravel more details about the transition mechanisms at various conditions. Before discussing the results in detail, we first provide some general remarks on the transition based on topological considerations. Considering the winerack-type structure of the framework and the very limited intrinsic flexibility of the BDC linkers in between aluminum chains, physically feasible transition mechanisms need to preserve the approximately rhombus shape for most of the pores in the system because the energy required to strongly deform the cells from their rhombus shape is expected to be very high. In addition, interatomic forces were modelled using a classical force field and as such we implicitly enforce strict covalent bond integrity throughout the entire transition. Deviations from the rhombus shape might eventually be possible provided that the mechanical energy supplied to the system is sufficiently high. Our simulations show that essentially two different transition mechanisms are active depending on the magnitude of the applied external pressure. At lower pressures, a layer-by-layer transition mechanism is observed whereas at higher pressures discrete nucleation points emerge within the lattice, which ultimately give rise to domain formation during the transformation. These results are schematically shown in Figure 4, where the left column shows the snapshots from the 100 MPa transition, and the right column shows the snapshots from the 300 MPa transition. While both transitions were recorded at 300 K, we observed entirely similar behavior at 200 K and at 500 K.

The simulations performed here constitute an important step forward towards a full mechanistic understanding of various phase transformations in realistic MOF particles. The implementation of the proposed pressure control algorithm in the GPU-accelerated OpenMM library as presented here allows to include fully anisotropic unit cell fluctuations in large-scale simulations containing millions of atoms. Its area of application is not limited to the breathing transition in winerack frameworks but is very general, and also includes e.g., crystalline-to-amorphous transitions. As such, we have taken a major leap forward in terms of the size of the systems that can be simulated in the field of MOFs. Based on the current simulations, we demonstrate that the transition mechanism may be critically dependent on various elements such as the size of the crystal but also the strength of the external stimulus, which is in this case the external pressure. For future simulations, it would be interesting to extend the protocol towards other stimuli such as temperature but eventually also guest molecules, and unravel whether similar transition mechanisms are observed. Such simulations might enable to construct a transition mechanism diagram that indicates the expected transition mechanism as a function of various control variables. Such a hypothetical diagram is illustrated in Figure 6A. The icons X,Y,Z indicate various mechanisms such as collective behavior, layer-by-layer behavior, or transitions starting from discrete nucleation points. Based on the simulations performed in this paper and literature data, we are already able to fix a number of points on the transition mechanism diagram in terms of the external pressure. However, more simulations are necessary to provide the complete transition mechanism diagram.

Despite the new insights obtained in this work, further methodological steps are necessary to enable the simulation of systems that are truly representative of realistic MOF particles, i.e., with crystal sizes similar to experimentally observed crystals, and with explicit inclusion of defects and the crystal surface. In what follows, we give some reflections on future perspectives in this direction.

First of all, the current approach will have to be extended towards finite nanocrystals with morphologies that are representative of experimentally observed crystallites and with length scales going beyond 50 nm. Subsequent embedding of this crystallite in a medium will then allow to apply an external pressure to the crystal in a very natural manner (Figure 6). Experimentally, both mercury as well as silicone oil have been used as effective pressure transducers in order to detect pressure-induced structural transitions in several winerack frameworks (Beurroies et al., 2010; Yot et al., 2012; Yot et al., 2014, 2016; Ramaswamy et al., 2017; Henke et al., 2018; Krause et al., 2018; Wahiduzzaman et al., 2018; Wieme et al., 2019; Yot et al., 2019) and other flexible frameworks such as ZIF-4 or DUT-8 (Kavoosi et al., 2017; Krylov et al., 2020); the simulation setup as described in Figure 6B would mirror such experiments.

Second, the current simulations were performed using an all-atom classical force field which does not allow any bond breakage to occur. It is not excluded that during a phase transition, bonds may temporarily break at specific locations within the framework as it is well known that linkers in MOFs may have a labile nature especially when being exposed to guest particles that allow to temporarily stabilize detached linkers (Hajek et al., 2018; Caratelli et al., 2019). To account for such effects, it is necessary to employ reactive force fields or more complex machine learning potentials (MLPs) that are trained based on quantum mechanical calculations in order to capture such effects. However, the systematic construction of MLPs for the complex systems under study here is highly nontrivial. To the best of our knowledge, only one MLP has so far been constructed for MOFs, by the group of Behler on MOF-5 (Eckhoff and Behler, 2019).

Finally, while current simulations allowed to deduce qualitatively new mechanistic details on the phase transition in MIL-53(Al), a next step would additionally aim to determine thermodynamic and kinetic properties associated with the transition, including its nucleation and growth. From a thermodynamic point of view, we could resort to the construction of Helmholtz free energy diagrams in terms of an appropriate collective variable. As noted in Section 2.3, the unit cell volume can no longer be regarded as an appropriate collective variable due to the absence of significant statistical fluctuations, and other variables will have to be considered such as the opening angle θ. The determination of kinetic properties, such as propagation rate constants for the phase boundaries, is highly challenging for these systems and may require specialized sampling protocols.

Previous considerations clearly illustrate the complexity associated with a full understanding of phase transformations in SPCs. The problem at hand is a prototypical example of a spatiotemporal process, where the dynamics of the MOF lattice is affected by spatial heterogeneities at various length and time scales (Van Speybroeck et al., 2021). A full understanding of the spatiotemporal response of MOFs will require a close partnership between the modeling and experimental community, whereby dedicated experimental in situ methods are necessary to track intermediate metastable states during their dynamic response towards external stimuli, and where theoreticians will have to explore new modeling avenues to tackle processes at various length and time scales.