Molecular Dynamics Simulations of the “Breathing” Phase Transformation of MOF Nanocrystallites

The displacive phase transformation of metal‐organic frameworks (MOFs), referred to as “breathing,” is computationally investigated intensively within periodic boundary conditions (PBC). In contrast, the first‐principles parameterized force field MOF‐FF is used to investigate the thermal‐ and pressure‐induced transformations for non‐periodic nanocrystallites of DMOF‐1 (Zn2(bdc)2(dabco); bdc: 1,4‐benzenedicarboxylate; dabco: 1,4‐diazabicyclo[2.2.2]octane) as a model system to investigate the effect of the PBC approximation on the systems' kinetics and thermodynamics and to assess whether size effects can be captured by this kind of simulation. By the heating of differently sized closed pore nanocrystallites, a spontaneous opening is observed with an interface between the closed and open pore phase moving rapidly through the system. The nucleation temperature for the opening transition rises with size. By enforcing the phase transition with a distance restraint, the free energy can be quantified via umbrella sampling. The apparent barrier is substantially lower than for a concerted process under PBC. Interestingly, the barrier reduces with the size of the nanocrystallite, indicating a hindering surface effect. The results demonstrate that the actual free energy barriers and the importance of surface effects for the transformation under real conditions can only be studied beyond PBC.


Introduction
In their seminal review on functional porous coordination polymers, Kitagawa and coworkers referred to flexible systems that can undergo structural transformations (upon certain stimuli like, e.g., guest adsorption) as third generation soft porous crsytals. [1] The so-called "breathing" effect, where metal-organic frameworks (MOFs) show large volume changes, is one of the The copyright line for this article was changed on 25 October 2019 after original online publication.

DOI: 10.1002/adts.201900117
reasons for the ongoing interest in this new class of crystalline porous materials. [2][3][4] Such structural transformations, which can lead to hysteretic behavior and gate opening upon guest adsorption, can have enormous impact also for the applicability of MOFs in gas adsorption, [5] gas separation, and other potential applications such as sensing. [6,7] In order to understand, predict, and exploit this key feature of MOFs, a variety of experimental and theoretical investigations of such phenomena has been performed. [8][9][10][11] In general, they can be triggered not only by guest adsorption but also by pressure, temperature, and other stimuli like electric fields. [12,13] Albeit a number of network topologies allow for volumechanging transformations, [14] the so called "wine rack" deformation of pcu like topologies is most often observed and studied. In particular the MIL-53 series was most intensively investigated both experimentally and theoretically. [11,[15][16][17][18] The "wine rack" motion as observed in MIL-53 and other MOFs is a shear displacive phase transformation as schematically shown in Figure 1. The orthorhombic open pore (op) form converts to the closed pore (cp) form by a shearing of one plane, whereas the other is bent sideways. Most systems are symmetric with respect to two directions and the shear displacement can occur in either direction. Interestingly, for MIL-53 the system is in the op form when activated and closes for moderate loading of guest molecules like CO 2 which interact strongly with the inorganic rod like structures. [19] At higher loading the system opens up into the op form. By increasing the temperature, the entropy contributions gain weight and favor the op form, whereas hydrostatic pressure leads to a transformation to the cp form. [11,18] It is evident that the hysteretic behavior in an adsorption experiment indicating such a structural transformation is due to kinetic effects, where due to a barrier, the adsorption and desorption branch differ and the state of the system is historydependent. The actual mechanism for such transformations and its kinetics are still not well understood, and predictions of gating pressures are difficult to achieve. One of the limitations in the theoretical simulations is the use of periodic boundary conditions (PBC) together with limited system sizes in the vast majority of studies. For a real system it is unlikely that a complete MOF crystallite converts from one phase to the other in a synchronous concerted fashion. According to the schematic free Shear displacive wine rack motion of a MOF in pcu topology together with the schematic free energy of this transformation (left) and a sketch of a real-world example of a crystallite featuring surfaces and a phase interface, which will have an effect on the free energy curve as indicated by the respective colors. energy diagram shown in Figure 1a, this would result in a barrier from a state where the strain energy is already substantial, but no substantial stabilization from dispersive interactions is present for the entire system. In contrast, as shown in Figure 1b, a first-order phase transformation with an interface between the op and the cp phase that travels through the system could be expected. It is evident that a much smaller free energy would have to be overcome and the transformation kinetics would then be mainly determined by the speed in which a first transition seed is formed and the resulting interface travels through the crystal.
For various MOFs it has been observed that crystallite size plays an important role for these kinds of structural transformations. [20][21][22][23] For DUT-8, smaller crystallites below 500 nm do not close, whereas above 1000 nm reversible breathing is observed. [24] In the intermediate regime a non-reversible closing on activation can be observed. The same tendency holds for the pillared-layer MOF-508 Cu 2 (bdc) 2 bipy, where crystal downsizing also leads to a stabilization of the op form, which do not close to the cp form if the crystallites are below a critical size of about 60 nm. [25] Recently, a similar behavior was observed for thin films of functionalized pillared-layer MOFs. Only for surface mounted films of 80 or more layers of Cu 2 (DE-bdc) 2 dabco (DE-bdc: 2,5diethoxy-1,4-benzene-dicarboxylate) a transition to a narrow pore (np) could be observed upon removal of methanol. [26] In general, such a behavior can be explained in analogy to nucleation theory, where the energy of a particle is expressed in terms of a penalizing surface term and a stabilizing volumetric term. [27] In this context, the surface hampers closing of the system, resulting in only sufficiently large crystallites being able to overcome the barrier. Clearly, such surface effects can only be considered in simulations beyond PBC. Note that such surface effects hindering the transformation for smaller crystallites have previously been found in MD simulations of the phase transformations of molecular crystals. [28][29][30] In order to investigate this behavior on a fundamental level, we chose here Zn 2 (bdc) 2 dabco (DMOF-1) as a model system. [31] The pcu topology allows for a displacive shear transformation, [14] and even though this particular system has not been observed in a closed pore (cp) state, our force field predicts a stable cp phase, making DMOF-1 a suitable, yet simple system for our investigations, serving as a demonstrator system to develop and benchmark the methodology of nanocrystallite simulations in general.
Our aim is to investigate nanocrystallites (NCs) of different size to compare their phase transformation with the corresponding process simulated under PBC in order to prove the hypothesis of a first-order phase transformation, starting with a nucleation event and proceeding by a growth of the new phases domain. In such a scenario, the adsorption of guest molecules is not straightforward, and approaches similar to the non-equilibrium methods for adsorption in systems with an interface proposed recently by Parrinello and coworkers would have to be employed to maintain the reservoir partial pressure. [32] Thus, in this first step we have focused on temperature and pressure stimuli, or in case of the NCs the equivalent of a mechanical force to trigger the phase transformations.

Methods and Computational Details
The NCs of different size Figure 2 have been constructed by the reversed topological approach (RTA). [33][34][35][36] The pcu topology has been extended by 2-connected "vertices" on the x and y axes for the linker and a further 2-connected "vertex" on the z axis for the pillar. By slicing along specified hkl-planes of sufficiently large supercells, blueprints for the nanocrystal were generated as shown in Figure 3. We refer to the size of the NC by the number of paddle-wheel units in each spatial direction. This means a system of size 6 × 6 × 6, as shown in Figure 3, contains 5 × 5 × 5 orthorhombic pores. Connections between vertices cut  by the slicing process have been saturated by 1-connected stubs, representing surface groups. With our weaver code, the atomistic system was generated based on this blueprint. On the [001] and [00-1] planes, a dabco pillar was used to saturate the surface, whereas on the other surface planes, we used a simple phenyl unit to replace the stub. This is equivalent to a benzoate modulator saturating the surface on the xz and yz planes. Note that in our study of the surface termination of HKUST-1 we used instead acetate modulators, [37] which is motivated by the growth of HKUST-1 from copper-acetate. In the case of the pillared-layer (PL) type DMOF-1 investigated here, the effect of modulators on the crystal morphology has been studied. [38] Using benzoate is a straightforward choice, but investigating the influence of the chemical nature of the modulator on the phase transformation is out of the scope of this work.
All force field calculations have been performed using the LAMMPS molecular mechanics program package (http://www.lammps.sandia.gov, version 30 Mar 2018) where the MOF-FF potential terms have recently been added within the USER-MOFFF package. [39] The DMOF-1 Zn 2 (bdc) 2 (dabco) system was described by the first principles parameterized MOF-FF force field. [40,41] Force field parameters were assigned automatically using our hierarchical automated parameter assignment (HAPA) strategy, [42] where building blocks are detected by substructure search based on the connectivity graph and force field parameters are automatically assigned from an SQLdatabase hosted at the MOF+ website (http://www.mofplus.org). As a reference for the NC simulations, a periodic model was also constructed as 2 × 2 × 2 supercell from system rotated by 45 • around the z-axis to maintain a mainly orthorhombic unit cell during transformation. For the van-der-Waals interaction a distance cutoff of 12Å was used. Coulombic interactions between the Gaussian charge distributions used in MOF-FF were also computed by a distance-based cutoff scheme using the shifted force method proposed by Gezelter et al. for performance reasons with a cutoff of 12Å. [43] This approximation was validated for the periodic system comparing to results using Ewald summation (short range cutoff 12Å). Pressure volume curves p(V ) for the periodic systems were computed by the simulation protocol described by Rogge et al.: [15] Initial structures along an op-cp phase transition were obtained from a series of NPT simulations, where the pressure was ramped up to 10 kbar within 10 ns in order to trigger the phase transformation within the simulation time. The hydrostatic pressure needed for the p(V ) equation of state (EoS) was determined in a second step. For an equidistant set of volumes between 0.8 V cp and 1.2 V op initial structures were taken from the preceding NPT simulations, which are close to the respective target volume. For each structure, a simulation in the (N, V, σ a = 0, T ) ensemble was performed using the volume constrained version of the MTTK barostat [44,45] (the implementation was taken from the repository at Github https://github.com/stevenvdb/lammps.git (newbarostat branch)). This barostat samples cell fluctuations by performing volume-preserving steps in such a way that the off-diagonal elements of the pressure tensor σ a vanish on average. [15,46] In this way the hydrostatic pressure p(V ) for a given volume can be obtained as the average of the trace of the pressure tensor. From the p(V ) EoS, the Helmholtz free energy A(V ) can be obtained by integrating pressure over the volume starting from an arbitrary reference volume V 0 according to The inner energy U(V ) is obtained as the average potential energy at a given volume. The entropy change can thus be computed by All periodic (N, V, T ) and non-periodic constant temperature simulations were carried out using the Nosé-Hoover thermostat using a relaxation time of 0.1 ps. [47] Temperature ramp simulations were performed by incrementing the bath temperature of the thermostat at each step of the simulation.
The collective variables module colvars was used via the LAMMPS USER-COLVARS package. [48] Harmonic restraints of 2 were applied for each distance d between the centers of mass of metal dimers from the paddlewheel units located on two opposing z-directed NC edges. In other words, for each xy-layer of the NC one of the diagonals was restrained to a given distance.
The NCs constructed in the op form with the weaver code [33] were first optimized and then equilibrated for 100 ps without any restraint at T= 300 K. Then the restraint was activated and the reference distance d 0 was gradually reduced from 15Å ·(n − 1) to 6Å ·(n − 1) with a force constant of k = 20 kcal molÅ 2 within 1 ns of a steered MD simulation at a temperature of 300 K, where n is the size of the NC. The restraints are shown in orange in Figure 4 together with the other descriptors and variables of the NC with size 4 × 4 × 4.
The work required to transform the entire NC was then computed by umbrella sampling (US) simulations. [49] For this, snapshots at every d 0 = 0.5Å were taken from the steered MD simulation. For each of them, a 0.5 ns (N, V, T ) simulation with a fixed d 0 was performed using a force constant of k = 10 kcal molÅ 2 and the trajectory of the collective variable d(t|d 0 , k) was recorded every 100 steps. The work W, which in this ensemble corresponds to A(d), was computed using the Weighted Histogram Analysis Method (WHAM) using the wham code, [50] which was run on the colvar trajectories d(t|d 0 , k) using a convergence threshold of 10 −6 kcal mol . For WHAM, adjacent distributions ρ(d 0 , k) approximated by the histogram of the collective variable trajectory d(t|d 0 , k), require to overlap in order for WHAM to construct a valid free energy surface. This was visually confirmed by examining the distributions shown in Section S3.2, Supporting Information. In some cases, there was no sufficient overlap and additional runs were performed in between the adjacent 0.5Å windows using a larger force constant of 20 kcal molÅ 2 in order to restore the overlap in these critical regions.
The volumes for individual pores were computed by first determining the centers of mass of each paddle-wheel fragment along www.advancedsciencenews.com www.advtheorysimul.com the trajectory. Each set of eight mass centers adjacent to a pore was used to compute a convex hull using the qhull library within scipy, [51] which yields the volume of one pore as described in Section S3.1, Supporting Information.

Reference System under PBC
As a reference for the NC calculations we first investigated the thermodynamics of the volume change for the periodic system, implying a concerted transformation due to the rather small 2 × 2 × 2 supercell. The simulation cell volume in Figures 5 and  6 were normalized to one pore by dividing it by a factor of four in order to better compare it to the volumes of the individual NC pores. Since Ewald type calculations of the Coulombic interactions are prohibitive for the non-periodic NC simulations, we first tested the use of cutoff based schemes.
Therefore, we first investigated the effect of using the shifted force-damped, cutoff-based scheme proposed by by Gezelter et al. [43] on the thermodynamic features of the system in periodic boundary conditions.
The p(V ) EoSs shown in Figure 5 for the different methods show qualitatively the same behavior. However, in the transition region between 800 and 1100Å 3 oscillations are observed for the cutoff based scheme, which is reduced by increasing the cutoff and is absent for the Ewald summation results. This is mainly due to the fact that at some points of the phase transition additional paddle-wheel units enter the cutoff sphere. Interestingly, with respect to the reference computed by the Ewald method, the averaged pressure is over-or underestimated in an oscillatory fashion. This is why the integrated relative free energy curves do not differ too much apart from small oscillations at the barrier for the cutof-based cases. Also the difference between the free energy minima is only slightly affected by the cutoff scheme.
For comparison with the non-periodic simulations, the temperature-dependent analysis of the thermodynamic functions have been performed by the cutoff based scheme (r cut = 12Å) and the oscillations in the transition region in Figure 6 are an artifact due to this approximation. The thermodynamic functions change only very little with temperature. The main difference is the relative stability of the closed pore, which-as expected-becomes less favorable with increasing temperature since strain effects reduce the available phase space. [52] The main contribution to this difference is consequently the entropic term −T S(V ). Note that even at room temperature the op form is thermodynamically preferred. Experimentally DMOF-1 is received in the solvent-filled op form as synthesized and stays this way also during activation. Note that for the corresponding Cu 2 (bdc) 2 (dabco), amorphization is observed at elevated hydrostatic pressures. [53] However, for the larger ndc linker in DUT-8, a cp form is observed for various metals. [10] Our non-reactive MOF-FF is not able to describe and reproduce amorphization if its includes bond cleavage of any kind. The cp form of DMOF-1 is thus a virtually high pressure phase, serving as a model to study such phase transformations processes in general. In parallel to the destabilization of the cp form at higher temperature, also the barrier is reduced and makes the phase transformation more and more likely.

Thermal Opening
Due to the preference of the op form for all temperatures we started out to study the process of phase transformation of NCs by simulating the thermally induced opening. We prepared cp NCs taken from the steered MD runs and performed 1 ns simulations at constant T = 300 K without restraints. All the NCs in size between 3 × 3 × 3 and 16 × 16 × 16 remained in their cp form in this time. Only the smallest 2 × 2 × 2 system instantaneously transformed to the op form. This could either be due to larger surface contributions or simply to the likelihood of accumulating enough energy in the right phonon modes. Note also that we can not exclude a transformation of the smaller systems at T = 300 K at extended simulation times.
In order to initiate the transformation the temperature was ramped in the following MD simulation from 300 to 500 K with a ramp velocity of v T = 0.2 K ps . All systems from 3 × 3 × 3 to 9 × 9 × 9 in size transitioned within 1 ns toward the op form. Only the 16 × 16 × 16 system did not undergo a transformation even during further 0.5 ns of constant T = 500 K sampling. This system was further ramped up to 800 K with the same ramp velocity v T , where it finally transformed to the op form at T = 580 K. Table 1 lists the transition temperatures, which were obtained by identifying the first occurrence of a pore volume www.advancedsciencenews.com www.advtheorysimul.com   3 is also well past the free energy barrier computed in PBC shown in Figure 6. In addition to that, the transition onset was visually confirmed by inspecting the trajectories. The reported temperatures are the temperature of the heat bath of the thermostat at that particular time, not the instantaneous temperature of the system. Interestingly, the phase transformation always occurs via the same pathway for all NCs: The pores open up first at one edge of the crystallite (see snapshot in Figure 7 at t = 37.5 ps) where two sides are exposed to the surface, serving as a nucleation for the phase transformation. Note that the entire edge along the z-axis opens together, which implies an increasing number of pores for larger crystallite sizes. From this initial point the phase transition then spreads throughout the crystallite in both the x-and y-directions. This creates a strain in the next cp layers such that there is also a phase boundary traveling along the xy-diagonal as it can be seen in Figure 7 for the snapshots at t = 52.5 ps and t = 60.0 ps.
After the transition, the resulting momentum persists for some time in the respective phonon modes, resulting in a partial closure of the structure right after the cp-op transition. Such structural oscillations in the terahertz region have been observed by Raman spectroscopy for the related MOF DUT-8. [54] The larger the system, the longer it takes for the energy in this phonon mode to dissipate. In the small system of size 3 × 3 × 3 it takes only about 200 ps until the momentum generated by the phase transition is redistributed, whereas in the largest system, there is still a significant oscillatory motion visible even 0.5 ns after the transition has finished. Note in this context that the thermal opening is driven by a thermostat with rising target temperature. The equipartition is not fully maintained in this phase due to the rapid phase transition and the energy release into the breathing phonon mode.
In order to estimate the speed v PT at which the phase interface travels through the system we used the pore volumes computed for the trajectory of the 9 × 9 × 9 NC for a numerical approximation. For this the time evolution of the volume of the eight pores along the side of the NC (x-direction) and along the diagonal (xydirection) were averaged over the z-axis. As shown in Figure 8 www.advancedsciencenews.com www.advtheorysimul.com Figure 7. Snapshots from the thermal breathing transition of the 9 × 9 × 9 nanocrystal a) with timings and b) the pore size analysis. The time t = 0 is set to the simulation time 1 ps before the phase transition initiates (volume above 1050Å 3 , see Table 1).

Figure 8.
Pore volume trajectories for the 9 × 9 × 9 NC system close to the phase transition (averaged in z direction). The legend denotes the pore index along a) the NC side x-direction and b) the diagonal xy-direction of the NC. the volumes switch from below 700Å 3 to above 1100Å 3 one after the other over the course of the simulation. By fitting with a logistic function, the half-times and widths have been determined by which the transition velocity and thickness of the individual transformations have been determined (details are given in the SI Section S2.1). The velocity in both x-and xy-direction increase as the transition progresses. In the faster x-direction, the phase transition travels at an initial speed of v PT x,initial = 92 m s and speeds up to about v PT x = 700 m s , whereas in xy, the velocities are v PT xy,initial = 34 m s and v PT xy = 240 m s , respectively. Also the thickness of the boundary is very different for the two different directions. For the xy direction, the width on average is just about the distance of one neighboring pore in that direction (6-7Å), whereas the much faster x direction also spreads out farther toward 30-40Å, which is about three to four neighboring pores. (For details see Section S2.1, Supporting Information) In principle, from these velocities at different temperatures an activation barrier could be determined. However, these velocities are on one side modulated by the further increasing temperature during the simulation and on the other hand the different sizes of the NCs with the different surface-to-volume ratios will also affect the speed of propagation of the interface. Also the thermostat has to some degree affected the speed of propagation, which should in future studies be avoided by switching off the thermostat and continue in the microcanonic ensemble. For this reason we have not further analyzed the tendencies.

Closing by a Mechanical Force
In contrast to the thermal opening, the phase transition can also be initiated by a mechanical force like, for example, a hydrostatic pressure. However, for a non-periodic NC this is not as straightforward to mimic such an experiment, as for the system in PBC where simply the cell parameters need to be changed in order to change the volume and a force/pressure can be applied simply by adding it to the stress tensor. As it is clearly visible in Figure 7 the NCs do not necessarily adapt a well-defined shape during the transition and the overall volume is only defined via the sum of the polygonal volumes of the individual pores. As explained in Section 2 we have therefore used a simple collective variable to represent a uniform force for a volume reduction by the breathing mode: for each xy-layer a restraint keeps the centers of mass of the metal dimers in all paddle wheels along opposite edges at a certain distance. By reducing this distance to lower and lower values, a closing of the MOF by a volume reduction is enforced. In analogy to integrating the average pressure in the periodic simulations to get the free energy, in the aperiodic case the necessary work for the closing is determined by US along this distance. Due to the substantially larger numerical effort for US where many windows need to be sufficiently sampled, this was only feasible up to a NC with size 6 × 6 × 6.
The phase transformation process initiated by mechanical force is similar to the thermal opening with a distinct phase boundary traveling through the system, but with certain differences. For the thermal opening we analyzed the snapshots of www.advancedsciencenews.com www.advtheorysimul.com Figure 10. Free energies of the systems 2 × 2 × 2 − 6 × 6 × 6. a non-equilibrium process, whereas here we have access to an averaged equilibrium situation for a given restraint distance. In Figure 9 corresponding snapshots of the NC with size 6 × 6 × 6 are shown together with an average of the individual pore sizes for different values of the restraint length d 0 . The pore diagonal r = d 0 n−1 is used for a comparison between the differently sized NCs. Again, the first pores that open are located at the edge of the NC. For d 0 = 40.3Å, which is close to the global transition barrier where the first pore has to open, strong volume oscillations are observed for these pores. Note that in contrast to the thermal opening, where the front propagates along both xand y-direction, the system herein breaks symmetry and chooses spontaneously one direction. For d 0 = 44.5Å or d 0 = 51.0Å it is obvious that rows open pore by pore, with the largest volume oscillations for the pores that are in the process of transformation.
In Figure 10 the work profiles for the phase transformations are given as determined from the US with respect to the restraint distance. Because of the increasing size of the systems an increasing A for the overall NC transformation is found. In order to compare the results of the different NCs with each other and with the results from the simulations in PBC it needs to be normalized. In Figure 10a normalization with respect to both the number of pores n pores = (n − 1) 3 and with respect to the number of paddle-units n PWs = n 3 is shown. For larger NCs the difference between this normalization will become smaller and smaller, but for the smallest systems this makes a substantial difference.
The unnormalized free energies along the reaction coordinate for NCs of different sizes shown in Figure 10a reveal two fundamental tendencies: First of all, the total barrier of opening increases for larger particles, which can be attributed to the fact that more and more pores have to undergo the transformation for larger crystallites and is related to the tendency that larger crystals start to transform at higher temperatures. Each free energy curve furthermore features an elastic regime for r between about 15 and 12.5Å, where the entire framework is sheared but still all pores being close to the op form, as it is visible for example for r = 13.7Å in Figure 9. Note that in contrast to the PBC simulated free energy profiles, the ripples in the NC free energies are not only due to cutoff effects, but originate from entire rows transitioning, as for example in snapshots at r =8. 4, 9.7, and 11.3Å in Figure 9. These partially transformed intermediates are the reason for the flattening of the normalized free energy curves as the crystallites get larger.
As already mentioned, in order to analyze size effects it is necessary to properly normalize the per NC free energy curves, which is however, to some extent ambiguous. Either we consider the energy per pore or per deformed paddle-wheel unit. The www.advancedsciencenews.com www.advtheorysimul.com latter makes sense stoichiometrically, but paddle-wheels at the surface do not have to deform by the same degree as fragments in the bulk structure upon transformation. Despite the differences between these normalization methods we find in both cases similar tendencies. First of all, the differences in free energy between the op and cp form ( A(cp -op)) gets smaller and converges toward the energy difference computed for PBC, because for a surface-to-volume ratio close to zero the surface effects become negligible and the bulk free energy difference is restored. This agreement also indicates the validity of our reaction coordinate chosen for the US approach. In addition, the free energy barrier is (apart from the smallest 2 × 2 × 2 system and normalization per pore) always smaller than the barrier for the concerted process computed in PBC. As discussed above, this is due to the fact that the majority of the pores stay either in the cp or the op free energy minimum and only the formation of the interface implies a major free energy penalty, whereas moving the interface through the crystallite poses only a smaller barrier. Finally, in line with the free energy difference between op and cp form, also the barrier reduces with the size of the system, indicating a substantial surface effect hampering the transformation.
The size of the systems investigated here is still much smaller than real MOF NCs for which flexibility starts to be observed. In other words the surface effects prohibiting a pore closing must be substantially larger in real systems to prevent the phase transformation up to a larger size of the bulk volume. However, our NC models show an idealized surface termination allowing for an easier pore closing of the surface. Nevertheless, a visual inspection of the closed structures reveals a strain from the terminating benzoate groups in case of a surface of the cp form. Thus, also for the idealized small model NCs the corresponding surface-to-volume effect is qualitatively observed; however, in a smaller magnitude as this is expected for real systems, where defects and surface corrugation will lead to much wider layers of non-ideal structures at the surface.

Conclusions
We have performed the first simulations of a breathing phase transformation for nanocrystallite models of the pillared layer DMOF-1 beyond periodic boundary conditions, using the first principles parameterized force field MOF-FF. NCs with ideal surface termination and sizes up to 16 × 16 × 16 paddle-wheel units or 23 nm diagonal distance where thermally opened and closed by mechanical force. The NC simulations clearly revealed a mechanism of first order phase transition with an interface between the op and cp form, traveling through the crystallite, in contrast to the concerted process in PBC for small supercells (2 × 2 × 2). In both processes distinct size effects are observed. In the thermal opening, a nucleation at an edge along the pillar z-axis is found, propagating symmetrically in the xy-plane throughout the system. The critical temperature at which the process initiates was found to increase with size, which is attributed to a higher free energy barrier for the nucleation because of the longer edge for larger systems. Free energies profiles for the closing of the NCs by a mechanical force determined by US give substantially lower barriers as compared to the PBC results, again due to the non-concerted first order transformation. Despite the ambiguity in normalization, a convergence of the free energy difference for the cp to op phase to the PBC values is found with increasing size of the NCs. The larger energy differences and barriers of the smaller NCs are due to surface effects that penalize the surface of the cp form.
Our simulations are only a first step, but show that for an understanding of the breathing phase transformation mechanism in porous MOFs on a molecular level, simulations of large system sizes, and ideally beyond PBC, need to be performed. For computing the thermodynamic free energy differences between two phases, simulations in PBC with small cells are sufficient, since both methods-with and without PBC-converge to the same result. On the other hand, only with nanocrystallite simulations proper barriers of activation and thus kinetics of the phase transformation can be derived.
Despite the size of up to a quarter of a million atoms, the here investigated model NCs are still substantially smaller than the real systems where phase transformations begin to be observable. However, the salient features of the atomistic mechanism as well as the hampering surface effects-even though of smaller magnitude-are all present for the here investigated model NCs. This corroborates the value of large-scale atomistic MD simulations using accurate force field models on one side but also the limitations in the majority of previous theoretical investigations constrained by the PBC approximation. The ingredients for such aperiodic simulations are thus an accurate force field to allow for large system sizes and a reasonable model for the surface termination of the MOF NC. In this work a distance restraint was used to model the effect of mechanical pressure. In principle other choices to model the phase transition exist, as for example a different choice of the collective variable or even a bath of explicit soft spheres mediating the pressure to the NC, which would directly mimic a mercury intrusion experiment. We are currently extending our investigations to further PL MOFs with different linkers and to grand canonical investigations where the phase transformation is initiated by guest molecule adsorption.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.