Free energy of metal-organic framework self-assembly

A variety of functional materials can now be prepared by relying on the self-assembly of distinct building blocks. In principle, the building blocks can be engineered and programmed to create a material with a pre-specified function.¹ Examples include complex structures prepared by the self-assembly from DNA functionalized nanoparticles,² patchy nanoparticles and colloids,³ perfectly ordered structures prepared through the directed self-assembly of block copolymers,⁴ and many others.^5–7 Metal-organic frameworks (MOFs) also take advantage of this building block approach.⁸ 

MOFs are nanoporous crystalline structures formed through a coordination-driven process of self-assembly that includes inorganic nodes and organic linkers. Multifarious building blocks enable control of pore size, geometry, and chemistry. Concretely, through careful selection of specific building blocks, it is possible to design structures that can then be subjected to iso-reticular expansion.^9,10 As such, MOFs can be tailored toward various applications. Examples include carbon capture,^11,12 hydrogen storage,¹³ light-harvesting,^14,15 catalysis,¹⁶ sensing,¹⁷ drug delivery and imaging,^18,19 and separations.^20,21

Molecular modeling has played an important role in the development of MOFs. Past efforts have sought to develop strategies to screen existing and hypothetical materials in order to identify optimal performers by taking advantage of the “instructions” encoded into MOF building blocks and considering the vast number of possible combinations of such blocks. Algorithms have been developed to enumerate large numbers of crystalline structures,^22,23 characterize them, and assess their performance in a high-throughput manner.²⁴ Some bottom-up approaches²⁵ have proposed sequentially connecting building blocks until a crystal is formed, while others have adopted a crystal topology as an input to generate a final structure.^26,27 Databases of existing and hypothetical structures have been further investigated through screening and machine learning techniques to gain valuable insights and identify promising materials for subsequent synthesis and implementation.^12,28–31 In some cases, these techniques have successfully predicted novel structures and their performance.^12,30 However, as influential as past modeling work on MOFs has been, past computational studies have done relatively little to elucidate the process or mechanisms of self-assembly. As such, it has been difficult to address crucial questions regarding the experimental feasibility of synthesizing computationally generated crystals, as well as the relative stability of different structures and the issue of polymorphism in MOFs.

Recent experiments to elucidate the process of self-assembly in MOFs have revealed a rich and complex range of possibilities in which a series of steps eventually result in the formation of a crystal. In situ X-ray scattering and powder X-ray diffraction experiments indicate that some MOFs are formed by nucleation and growth, while others involve intermediate metastable structures that depend on synthetic conditions (solvent, temperature, etc.).^32–38 Similarly, varying the concentration of the building blocks of interest can lead to different structures.³⁹ As alluded to above, computational efforts to model MOF self-assembly processes have been scarce and have met with limited success. The earliest attempts to model the MOF self-assembly were by Mellot-Draznieks and co-workers,^40,41 who assigned interaction sites to building blocks that were placed in a simulation box. They then relied on simulated annealing to equilibrate their system into a stable porous structure. While successful at predicting several crystal structures, little information was gleaned into the self-assembly mechanism itself. Similarly, other work by Metere et al.⁴² found a new MOF topology by relying on a spherically symmetric potential but did not provide information about the mechanism of self-assembly.

To understand some of the mechanistic steps involved in MOF formation, Cantu and co-workers⁴³ modeled the self-assembly of an individual MOF node using ab initio molecular dynamics (MD) calculations and a polarizable continuum model to account for solvent effects. Important mechanistic information was extracted from that effort, but ab initio calculations are currently too computationally expensive to access the self-assembly of a complete MOF structure in an explicit solvent. Yoneya et al.⁴⁴ simulated MOF building blocks, also in an implicit solvent, and varied the interaction between the building blocks until a MOF spontaneously formed over the course of a long molecular dynamics (MD) simulation. The MOF that was formed is not currently synthesizable, and it is unclear how the tuning of building block interactions influences the self-assembly mechanism. Most recently, Biswal and Kusalik^45,46 performed long molecular dynamics simulations of Zn²⁺ ions and benzene dicarboxylate (BDC) linkers in explicit dimethylformamide (DMF) and an implicit solvent and observed amorphous structures containing local motifs that are believed to be involved in the MOF self-assembly. Some of their simulations were carried out to time scales beyond 1 µs, but were unable to reveal the formation of a MOF structure, serving to highlight the difficulty in modeling the MOF self-assembly. Such simulations might have been trapped in metastable states, unable to overcome free energy barriers to form a stable MOF structure. How the MOF forms from those amorphous states remains an open question.

All of the computational work to date provides information on snippets of the full picture of the MOF self-assembly but does not tell a complete story. Some studies have focused on the final structure, others on the formation of building units, and others on transient metastable states. In this work, we build upon past studies by relying on enhanced sampling techniques to study the self-assembly process of MOFs and the free energy barriers that arise along the pathway to stable MOF structure formation. We model the self-assembly starting from dissociated building blocks and amorphous chain-like states, to the final MOF structure, and identify plausible mechanisms for both scenarios and their corresponding free energies.

Our study focuses on the self-assembly process of the quintessential metal-organic framework, MOF-5.^47,48 Replica exchange with solute scaling (REST2) simulations are used to validate our force field, and a finite-temperature string (FTS) method is used to delineate the self-assembly process and the underlying free energy profile. We find that the process of self-assembly involves multiple metastable intermediates and that crystalline nuclei must involve four or more unit cells in order to be thermodynamically stable. We also find that the formation of an ordered crystal is in fact favored entropically, through the contributions of the solvent to the assembly process.

Born-Oppenheimer molecular dynamics simulations of MOF-5 clusters in an explicit DMF solvent show that the coordination to the inorganic building block, Zn₄O, is dynamic.⁴⁹ Solvent molecules can coordinate to the building block, and BDC linkers can detach and reattach to Zn₄O. To capture this type of behavior without the need of a reactive force field, we chose a model that does not include explicit bonding terms between the inorganic and organic building blocks. MOF-5 was modeled using parameters from Dubbeldam et al.,⁵⁰ who parameterized bonding and nonbonding potentials to reproduce adsorption isotherms. The model is also known to capture the negative thermal expansion of the structure. We made a slight modification in order to hold the distance between the oxygen and zinc atoms in the inorganic node constant at 1.93 Å, based on the crystallographic structure of MOF-5. Comparisons between classical simulations and quantum mechanical calculations of the structure of the metal cluster and linker binding energies can be found in the supplementary material. Figure S1 shows the building blocks that comprise the MOF-5 structure and the configuration of MOF-5. The synthesis of MOF-5 takes place in dimethylformamide (DMF) at room temperature (298 K).⁵¹ Therefore, DMF was included and modeled explicitly in our simulations. Parameters were taken from the OPLS-AA forcefield, and cross-terms were calculated using Lorentz-Berthelot mixing rules. Configurations and parameters were taken from the literature.^52,53 MOF configurations were generated by placing individual inorganic nodes and organic linkers to form the cubic structure of MOF-5. Once formed, enough BDC linkers were randomly placed near nodes so as to satisfy charge neutrality. Two different systems were studied, a single unit cell in an explicit DMF solvent and a four-unit cell in an explicit DMF solvent. A cubic simulation box of a pre-equilibrated solvent was used for both systems with a length of 5.08 nm. The MOF systems were then solvated resulting in 889 DMF molecules and 766 DMF molecules used for the single and four-unit cell systems, respectively. The single unit cell system contains 8 Zn₄O nodes and 24 BDC linkers; the four-unit cell system has 18 Zn₄O nodes and 54 BDC linkers. NPT simulations at 1 bar and 298 K were then performed to equilibrate the systems, and the resulting configurations were used in the rest of the simulations in the NVT ensemble.

In general, MOFs are described by their chemistry and their unit cell dimensions. In our case, we are not modeling periodic structures of MOF-5 but rather discrete configurations of the material in a cubic simulation box with explicit DMF solvent molecules. Therefore, to compare simulations with experimental structures, we describe our system in terms of collective variables: average angle between adjacent nodes in the unit cell (three angles per node), average distance between adjacent nodes in the unit cell, and coordination. We chose these CVs because they allow us to monitor changes in geometry and coordination between nodes and linkers during the structural transformations while also allowing for simple comparisons to the unit cell of the experimental crystalline structure.

Coordination was calculated according to

where r_ij is the distance between the center of mass of the node and the center of mass of the linker and d₀ is a parameter which was set to 0.85 nm,⁵⁴ roughly the distance between the center of mass of a linker and the center of mass of an inorganic node in the fully assembled MOF. Values of the average distance between nodes and average angle between nodes for the assembled MOFs, single and four-unit cell systems, are 1.3 nm and 90°, respectively. The single unit cell system has 12 linkers that are shared between two nodes (12 edges of a cube). Since each linker is coordinated to two nodes, the coordination number CV is 24. There are also additional 12 linkers that are placed randomly coordinated to single nodes of the MOF, bringing the total coordination number CV of the single unit cell system to 36. Similarly, the four-unit cell system has 33 linkers that are shared between nodes and 21 linkers coordinated to single nodes, totaling 87 for the coordination number CV.

We performed replica exchange with solute scaling simulations to validate the MOF-5 force field of choice and verify that the MOF structure represents a metastable state in an explicit DMF solvent. In replica exchange or parallel tempering,^55,56 multiple copies of the system are generated and held at different temperatures. At a set number of steps, exchanges between neighboring configurations are attempted and accepted according to a Metropolis criterion. Conformations at high temperature are able to overcome energetic barriers, which would be difficult at lower temperatures. Through the exchange of the configurations, the low temperature systems can then explore different conformations difficult to observe in standard simulations.

This method scales poorly with system size, however, because it requires a large number of replicas to facilitate exchanges, which is problematic for our study since it includes the MOF structure as well as explicit DMF solvent molecules. Replica Exchange with Solute Tempering (REST)^57–60 seeks to solve this problem by eliminating the dependence of the acceptance probability on the explicit solvent molecules, therefore improving the sampling efficiency of the system. Replica Exchange with Solute Scaling (REST2),⁶¹ the technique used in this investigation, builds on this idea by scaling the intramolecular potential of the solute: MOF-5 in our case. So, instead of having replicas at different temperatures and different potential energy surfaces (REST), replicas in REST2 are all at the same temperature but on different potential energy surfaces. This is achieved using Hamiltonian Replica Exchange (HREX). In our particular case, the scaling of the intramolecular potential is performed in the dihedral and nonbonded potential of the BDC linkers.

REST2 was used to determine the free energy minimum configuration for our single unit cell and four-unit cell MOF configurations. Sixteen replicas geometrically separated in the equivalent temperature of the MOF from 300 to 500 K were used in our simulation. The initial configurations were taken from an energy minimization of the MOF solvated in DMF. Exchange attempts took place every 1 ps. Simulations ran for 60 ns, and statistics were gathered every picosecond for the last 50 ns. We track the average distance between adjacent nodes in a unit cell and the average angle between adjacent nodes in a unit cell. Probability distributions are estimated from histograms of the collective variables. Free energies are then calculated according to

where F(s) is the free energy, k_B is Boltzmann’s constant, T is the temperature, and P(s) is the probability distribution.

String methods are used here to investigate transitions between stable basins. The Finite Temperature String (FTS) method, in particular, allows for thermal fluctuations to influence the final path that is identified. As such, the final string depends on the global features of the free energy landscape.^62,63

Briefly, an initial string is generated between two chosen states. It is then discretized into nodes which are equally spaced along the string. Each node will explore a Voronoi cell associated with it. In string method implementations, the edges of the Voronoi cells have a hard or soft potential so that each node will only explore its respective cell. Then, after a prescribed number of steps, the running averages of the variables of interest are updated. The string is then updated toward the running averages while also guaranteeing that the nodes on the curve are an equal arc length apart and that the curve is smooth (usually using a cubic spline). This process is repeated until the string converges.

FTS calculations for the MOF self-assembly used as collective variables (CVs) the average distance between adjacent nodes, the average angle between adjacent nodes, and the coordination to the nodes. The time step for evolution and the smoothing coefficient were set to 0.1, and the M tensor was assumed to be the identity tensor.

FTS calculations were performed for two different scenarios: (1) an amorphous state (formed from uncoordinated building blocks) assembling into a MOF and (2) uncoordinated building blocks assembling into a MOF. In the first scenario, unbiased molecular dynamics simulations of uncoordinated building blocks are allowed to run for 50 ns. The final configuration is then used as one of the ends of the string with the other being the MOF structure. For this scenario, all nodes of the string, including the bookend/endmost nodes, were allowed to evolve until convergence. In the second scenario, we are interested in probing the self-assembly process starting from a disassembled state, which is inherently unstable. For this reason, the configurations of the bookend nodes are “pinned.” That is, the position of the nodes at the ends of the strings is not updated; only the nodes between the ends change position in CV-space.

Convergence was determined by monitoring the root-mean-square distance (RMSD) of the images along the string compared to the initial string. The equation is as follows:

where θ_α is the CV coordinates of node α, t is the iteration, N is the total number of nodes in the string, and K is the total number of CVs. Convergence was assumed once D(t) plateaued or the RMSD difference between two subsequent iterations was 10⁻⁶.⁶³ 

Once a string has converged, the free energy of each node may be calculated by performing a restrained molecular dynamics simulation. The restraining potential is

where i is the CV, k_i is the force constant for CV i, θ_i(x) is the value of CV i during the MD trajectory, and θ_α,i is the value of CV i in the node α of the converged string. The gradient of the free energy can be approximated as

where T is the total time of restrained sampling, which should be large enough so the average values of the CVs have converged. The change in free energy with respect to the end node of the string, F_α − F₀, can then be calculated using the trapezoidal rule.⁶³ 

REST2 calculations were performed using GROMACS 4.6.5⁶⁴ and PLUMED 2.2.⁶⁵ FTS method calculations were performed using GROMACS 5.1.3⁶⁶ and SSAGES.⁶⁷ MD simulations were carried out in NPT ensemble at 1 bar and 298 K using the Berendsen barostat and the velocity rescale thermostat and the NVT ensemble at room temperature (298 K) using the Nose-Hoover⁶⁸ thermostat. Particle-Mesh Ewald (PME)⁶⁹ with a cut-off of 0.9 nm and a Fourier spacing of 0.33 nm was used for the electrostatics. Bonds to hydrogen atoms and Zn–O bonds in the node of the MOF were held rigid using the LINCS algorithm.⁷⁰ Equations of motion were integrated using velocity Verlet and a time step of 2 fs. Nonbonded interactions were modeled by a Lennard-Jones plus Coulomb potential, and cross-terms were calculated using Lorentz-Berthelot mixing rules.

We first confirmed, using REST2 simulations, that the MOF structures are at least metastable for the chosen force field. REST2 simulations also show that the CV values for angle and distance are in agreement with the experimental values of the MOF-5 structure. These results are notable since the original force field was parameterized as a periodic structure in vacuum, whereas we are simulating a discrete nanoparticle in an explicit DMF solvent. Additional details are provided in the supplementary material. Then, we randomly placed Zn₄O nodes and BDC linkers, enough to form a single unit cell or four-unit cell and maintain charge neutrality, in a box of pre-equilibrated DMF solvent, and performed an unrestrained MD simulation, in the NPT ensemble, at 1 bar and 298 K for 50 ns. The simulations did not yield MOFs but, rather, amorphous chain-like structures resembling those observed in the literature.^45,46 We then used finite temperature string (FTS) method calculations to model the self-assembly process from the amorphous structures to the final MOF structure.

The string calculations take place along three collective variables: the average angle between adjacent nodes, the average distance between adjacent nodes, and the total coordination between nodes and linkers. The initial images of the strings for the single-unit cell and four-unit cell systems are generated using restrained molecular dynamics. Images are generated so that they are equally spaced in each of the collective variables between the amorphous state and the assembled MOF state. Five images were used for the single unit cell system and ten for the four-unit cell system. Figure 1 shows the nodes of the converged strings for both systems and images of the disordered and assembled states. The single unit cell system (left panel, Fig. 1) shows that the disordered state has the highest values of distance between nodes (1.8 nm) and coordination between nodes and linkers (41.7), while the assembled state has the lowest values in those CVs (1.28 nm and 36) and an average angle between nodes close to 90°. Similarly, the four-unit cell system shows that the disordered state has the highest values of distance between nodes (1.9 nm) and coordination (97.8), while the assembled state has the lowest values of those same CVs (1.28 nm and 91.7) and an average angle between nodes near 90°. The nodes on the strings for both systems show that the transition from the disordered state to the MOF state involves a decrease in the distance between the nodes and the coordination between nodes and linkers.

The free energies of the final strings can be calculated using Eq. (5); the corresponding values are reported in Fig. 2. To calculate the free energies for both the single-unit cell and four-unit cell systems, the converged string was further discretized by linear interpolation between the nodes on the string. The single unit cell system shows a sharp rise in free energy going from the amorphous state (α = 0) to the final MOF structure (α = 1). For this system size, the amorphous structure is more stable than the assembled MOF structure. The four-unit cell system, on the other hand, shows that the assembled MOF structure (α = 1) is more stable than the amorphous state (α = 0). As the four-unit cell assembles, there are multiple free energy barriers—ranging from 8 to 50 kJ/mol—that must be overcome in order for the MOF structure to form. Figure 2 also shows remarkably different behavior depending on the system size. Comparing a single-unit cell to a four-unit cell system, we see that the transition from an amorphous state to the ordered MOF structure goes from being unfavorable to being favorable. This result suggests that concentration must be sufficiently high for the formation of MOF structures and indicates that the formation pathway for a stable MOF material involves multiple unit cells. The results also suggest, given the free energy differences, that amorphous states are expected during the MOF self-assembly.

From Fig. 1, we know that the structure needs to decrease the distance between the nodes and also decrease its coordination to form the MOF. So, to better understand the rearrangements that take place and how they relate to the free energy transitions observed, it is helpful to look at images of selected nodes along the string. The images on the converged string (Fig. 3) for the single unit cell system show the evolution from the amorphous structure to the final MOF unit cell. The connections between nodes and linkers of the amorphous structure (α = 0) decrease (α = 0.25, 0.4375) to form a chain-like structure (α = 0.5625). This transition (α = 0.4375 to α = 0.5625) corresponds to the rise in free energy in Fig. 2; the increase corresponds to ∼50 kJ/mol. Then, the structure closes on itself (α = 0.75) and finally forms the MOF structure (α = 1). It is important to note that although the uncertainties in Fig. 2 for the MOF free energy of self-assembly are large, the features that are identified correspond to important transitions and structural transformations necessary for the MOF self-assembly to take place starting from the amorphous configuration.

Figure 4 shows images of the converged string for the self-assembly process of the four-unit cell system. The Zn₄O nodes and linkers are interconnected in the first amorphous structure (α = 0) and then separate (α = 0.111, 0.389, 0.444) to allow for the rearrangement that leads to the formation of individual unit cells (α = 0.5, 0.67, and 0.78), culminating in the formation of the four MOF-5 unit cells (α = 1). The free energy rise in Fig. 2 (bottom) corresponds to the separation of nodes and linkers (α = 0.389, 0.444, and 0.5) that allows the formation of the unit cells, which is consistent with what is observed in the single-unit cell case: free energy barriers start to appear as nodes and linkers separate to allow for a rearrangement of the structure. The highest barrier in the four-unit cell case corresponds to the beginning of the formation of the unit cells (α = 0.5). For the single unit cell, the full structure transitions into a chain-like structure to rearrange and form the MOF-5 framework. Meanwhile, the four-unit cell system retains part of the interconnections, while another portion separates the connections to allow for rearrangement. Following this procedure, the amorphous system is able to rearrange into the final MOF structure by overcoming smaller free energy barriers. This explains the size dependence observed in the calculated free energies.

We now turn our attention to the role of the solvent in the MOF self-assembly. Figure S7 reports the radial distribution functions of the centers of mass of the MOF building blocks Zn₄O nodes, and BDC linkers, with the center of mass of the DMF solvent in the amorphous disordered state and final MOF structures for the single and four-unit cell cases. Both scenarios exhibit similar behavior, with the fully formed MOF structure showing higher coordination with solvent molecules than the amorphous disordered state. The single unit cell system shows a sharp peak close to 0.5 nm for the Zn₄O nodes that is higher in the formed MOF structure than in the amorphous state. Similarly, the BDC linkers show a broad peak close to 0.5 nm that is higher in the formed MOF than in the amorphous structure. The four-unit cell system shows similar behavior to the single unit cell system, but the difference in the peak heights for nodes and linkers in the formed MOF and amorphous structures is less pronounced. In both systems, the solvent has higher coordination with the formed MOF structure than with the amorphous state. This higher coordination of the solvent with the formed MOF is also reflected in the corresponding interaction energies.

Figure 5 [panels (a) and (b)] shows the nonbonded MOF-solvent energy of interaction (solid lines) and nonbonded MOF-MOF energy of interaction (dashed lines) corresponding to the final string configurations for the single and four-unit cell systems. The energies are shown as a function of reaction coordinate (alpha). As suggested by the radial distribution functions, the MOF-solvent interactions are stronger in the fully formed MOF structures (α = 1) compared to the amorphous disordered states (α = 0). The MOF-solvent nonbonded interactions decrease as the MOF structure forms and the MOF-MOF interactions increase. Interestingly, the interaction energies at the formed MOF (α = 1) in the single unit cell case increase considerably; that is not the case for the four-unit cell case, where the energy first increases and then decreases as the MOF formation is completed.

Figure 5 also shows [panels (c) and (d)] the free energy, internal energy, and entropy as a function of the reaction coordinate. The free energy (F) is the same as in Fig. 2, the internal energy (U) is the total potential energy, and the entropy is calculated from their difference: TS = F − U. Interestingly, the formed MOF structure shows favorable entropic contributions to the free energy compared to the amorphous structure in the single and four-unit cell systems. The single unit cell system [panel (c)] shows a sharp transition in the internal energy and the entropy at α = 0.4375, which corresponds to the disruption of the interconnected amorphous structure (Fig. 3). The four-unit cell case [panel (d)] shows a sharp decrease in the internal energy and increase in the entropic contribution at α = 0.6 corresponding to the formation of the MOF-5 unit cells. The range in the reaction coordinate (α = 0.4–0.6) that shows the free energy barrier also corresponds to favorable MOF-solvent interactions [Fig. 5(b)]. This indicates that the observed free energy barriers are in part associated with the MOF-solvent interactions that need to be disrupted for the MOF self-assembly to take place.

We now proceed to study the MOF self-assembly starting from uncoordinated building blocks, also using the string method. Initial images are generated using restrained MD from the fully formed MOF into a completely dissociated system, so the coordination number is as close to zero as possible. The initial string contains images that are linearly spaced in collective variable space, from the MOF structure to the dissociated state. The string for the single-unit cell system contains 10 images, while the four-unit cell system contains 20 images.

To model the self-assembly process of MOF-5 from uncoordinated building blocks to the final structure, the collective variables of the end states are held fixed. The string is then allowed to evolve between these two states. Figure S8 shows the final strings for both systems. Both strings start at a coordination near zero and an average distance between nodes close to 2.5 nm. As the MOF formation progresses from this state, the coordination collective variable steadily increases, while the average distance between nodes decreases until the final structure is formed. As expected for the cubic unit cells, the average angle between nodes reaches 90°.

The free energies of the final strings can be calculated using Eq. (5), and the corresponding values are reported in Fig. 6. We clearly observe that the disassembled state (α = 0) is the most unstable for both scenarios, and the assembled MOF (α = 1) states are the most stable. The self-assembly processes of the MOF structures proceed in a downhill manner in both the single-unit cell and four-unit cell cases, with no free energy barriers. The lack of free energy barriers would imply that an unrestrained molecular dynamics simulation that starts from the uncoordinated state should end at the assembled MOF state. As has been previously noted, this is not what happens. Rather, chain-like amorphous structures are obtained after an unrestrained MD simulation of the uncoordinated building blocks. This discrepancy is due to “pinning” the ends of the string for the FTS method calculation, which is necessary to observe the self-assembly process from an inherently unstable state as is the disassembled state.

To visualize the nodes along the final string for the single unit cell system, six snapshots are shown in Fig. S9 which illustrate the MOF formation mechanism starting from a fully disassembled state. From the uncoordinated building blocks in solution (α = 0), dimers and trimers of Zn₄O nodes with BDC linkers are formed (α = 0.11). Subsequently, clusters containing multiple nodes appear (α = 0.33) which coalesce into a single chain-like cluster (α = 0.78). This chain then forms the final MOF-5 structure (α = 1). Similarly, Fig. S10 shows various configurations along the converged string for the four-unit-cell system. Like in the single-unit-cell case, the ends of the string were held fixed at the dissociated and MOF configurations. The uncoordinated nodes and linkers (α = 0) aggregate to form dimers and trimers (α = 0.37). Then, clusters start to form (α = 0.68) which amass into chain-like structures (α = 0.79). The structures start to form the unit cells (α = 0.89) of the final MOF structure (α = 1). Both systems show the same self-assembly pattern: dimers and trimers to clusters, leading to a chain-like structure. In the single unit cell case, the chain-like structure rearranges to form the final MOF structure. Meanwhile, for the four-unit cell, the chain-like structure rearranges and gradually forms the individual unit cells, ending in the four-unit cell MOF system. Note that the chain-like structures observed here resemble those identified in the literature.^45,46

To further analyze the role of the solvent in this self-assembly scenario, we calculate the radial distribution function of the centers of mass of the MOF building blocks (nodes and linkers) with respect to the center of mass of the DMF solvent (Fig. S11). The Zn₄O building blocks in the dissociated state (dashed lines) show higher association with the solvent molecules than in the formed MOF structure (solid lines). The BDC linkers show similar association with the solvent molecules in the single unit cell case, while they show stronger association with the formed MOF structure than the dissociated building blocks in the four-unit cell system. As with the previous scenario, the MOF-solvent association is reflected in the interaction energies.

Figure 7 [panels (a) and (b)] shows the MOF-solvent (solid lines) and MOF-MOF (dashed lines) nonbonded interactions for the single unit cell and four-unit cell cases along the reaction coordinate of the final string. As suggested by the radial distribution function, the MOF-solvent interactions are lowest at the dissociated state (α = 0) and increase as the MOF forms (α = 1), opposite to the previous scenario. Figure 7 [panels (c) and (d)] also shows the free energy, internal energy, and entropy along the reaction coordinate. Notably, as MOF formation proceeds in the single unit cell case, the entropic contribution to the free energy dominates, while the internal energy contribution is unfavorable. In the four-unit cell case, both the entropic and internal energy contributions to the free energy are favorable.

Before concluding, a brief comment on the choice of CVs is appropriate. This work and others in the literature have shown that the MOF self-assembly is an activated process; long, unrestrained MD simulations of MOF building blocks do not form the MOF unit cells but instead remain in a disordered configuration. Enhanced sampling techniques facilitate overcoming free energy barriers that prevent the self-assembly process from being observed. These techniques use collective variables or order parameters of interest. In our case, we have chosen CVs that we can relate back to the MOF structure. That being said, determining which CVs to use is not a trivial matter. In our case, non-MOF configurations may be found that share the same CV values as the ordered MOF, indicating that the chosen CVs may be improved. Determining the ideal CVs for the self-assembly process will be the focus of future work. Furthermore, our study considers pathways that are dependent on the basins; in our case, these are the disordered and ordered MOF states. Potentially, there are other basins of the amorphous disordered state in particular that have not been considered. This could lead to different pathways being observed as well as different free energy profiles.

Another important consideration for the MOF self-assembly is that the building blocks used for MOF synthesis are salts. The Zn₄O arises from zinc acetate dihydrate while BDC arises from terephthalic acid. So, there are components that we have not included in our self-assembly simulations. Future efforts will aim to address the role of these counterions in the assembly process.

The self-assembly of MOF-5 was successfully modeled using a combination of enhanced sampling techniques. MOF-5 systems of a single discrete unit cell as well as four-unit cells were simulated in an explicit DMF solvent. REST2 simulations were used to validate the model; the calculated free energies exhibit a minimum for structures that are consistent with those observed in experiments. FTS calculations were performed for the two MOF-5 systems under two different scenarios. The first considered the self-assembly process from an amorphous state to the final MOF-5 structure, while the second examined the self-assembly process from an uncoordinated state to the final MOF-5 structure. In the first scenario, the single unit cell system is less stable than the amorphous system. By contrast, the four-unit cell system is shown to be more stable than the amorphous system and the free energy of the self-assembly process shows various free energy barriers in the way to the final MOF structure, associated with structural rearrangements of the MOFs and MOF-solvent interactions. In the second scenario, dimers and trimers are formed, followed by clusters which form a chain-like structure, culminating in the MOF-5 unit cells. Free energies of the process show the MOF-5 structures to be more stable than the fully uncoordinated system, and the process is downhill. Notably, in all cases considered, the fully formed MOF structures showed favorable entropic contributions to the free energy.

Simulations of the MOF assembly pathway in a solvent have not been conducted before in the context of methods capable of revealing the underlying free energy surface. It is therefore difficult to put our findings in the context of previous reports. We do note, however, that past results from long molecular dynamics simulations did reveal some of the transient structures that have also been identified in our own simulations. Those past calculations, however, were unable to overcome some of the large free energy barriers that are encountered along the way to full stable MOF structure formation. As such, our results at this point should be viewed as a proof of principle that one can indeed follow the formation from beginning to end (full MOF formation) by relying on advanced sampling techniques. Three new and important insights to have emerged from this study are that (1) the MOF nucleus must be sufficiently large (at least four unit cells) to be stable, (2) the solvent plays an important role in the formation process, and (3) the formation of a stable crystal nucleus is entropically favorable. The results presented here also serve as a starting point for consideration of larger systems and different materials. We also hope that the transition states reported in this work will stimulate new experimental studies aimed at establishing their validity.

Future work will focus on generating trajectories that form the MOF and using techniques to investigate new collective variables for the MOF self-assembly. Rates of the self-assembly can also be investigated using the generated trajectories. Future work will also consider multiple amorphous configurations and self-assembly pathways. Reactive force fields will also be used to account for the bond formation energy between the metal cluster and organic linker. Finally, future work will also consider the role of counterions which are present in the experimental assembly of MOFs.

See supplementary material for details regarding the MOF-5 model, REST2 simulations, and additional details regarding string method calculations.

This work was supported by the Department of Energy, Basic Energy Sciences, Materials Science and Engineering Division. The SSAGES software employed here was developed with support from the Midwest Integrated Center for Computational Materials (MICCoM).