Large-Scale Free Energy Calculations on a Computational MOF Database: Toward Synthetic Likelihood Predictions
Metal-organic frameworks (MOFs) have captivated the research community due to a modular crystal structure that is tailorable for many applications. However, with millions of possible MOFs to be considered, it is challenging to identify the ideal MOF for the application of choice. Although computational screening of MOF databases has provided a fast way to evaluate MOF properties, validation experiments on predicted "exceptional" MOFs are not common due to uncertainties on the synthetic likelihood of computationally constructed MOFs, hence hindering material discovery. Aiming to leverage the perspective provided by large datasets, here we created and screened a topologically diverse database of 8,500 MOFs to interrogate whether thermodynamic stability metrics such as free energy could be used to generally predict the synthetic likelihood of computationally constructed MOFs. To this end, we first evaluated the suitability of two methods and three force fields to calculate free energies in MOFs at large scale, settling on the Frenkel-Ladd path thermodynamic integration method and the UFF4MOF force field. Upon defining a relative free energy, ∆LMFFL, that corrects for some force field artifacts specific to MOF nodes, we found that previously synthesized MOFs tended to cluster in a region below ∆LMFFL = 4.4 kJ/mol per atom, suggesting a general first filter to discriminate between synthetically likely and unlikely MOFs. However, a second filter is needed when several MOF isomorphs are below the ∆LMFFL threshold. In 84% of the cases, the synthetically accessible MOF within an isomorphic series presented the lowest predicted free energy. The present work suggests that crystal free energies could be key to understanding synthetic likelihood for MOFs in computational databases (and MOFs in general), and that the thermodynamics stability of the fully assembled MOF often determines synthetic accessibility.
INTRODUCTION<!>METHODS<!>RESULTS AND DISCUSSION<!>Comparison between free energies calculated with different force fields.<!>Can free energies be used to identify synthetically feasible MOFs?<!>How does MOF isomorphism affect synthesizability predictions?
<p>Nanoporous crystals, spearheaded by a class of materials known as metal-organic frameworks (MOFs), have the potential to provide technological breakthroughs in fields as diverse as catalysis, [1][2][3][4] chemical separations, [5][6][7][8] molecular sensing, [9][10][11] drug delivery, [12][13][14][15] and energy storage. [16][17][18][19] MOFs are promising in such diverse fields due to the tunability of their pore geometry and chemistry, a property which allows researchers to precisely control how molecules adsorb, diffuse, and react within the framework. This tunability is the direct result of the combinatoric nature of MOF construction, where each combination of MOF building blocks and topology defines a unique MOF, which can be thought of as a point in "MOF-space." Given the number of these combinations (i.e. the vast size of MOF space), synthesizing and experimentally testing every possible MOF is impossible. Therefore, discovering a singularly high-performing MOF for a given target application is a major challenge that requires new, sophisticated methods for sifting through MOF space.</p><p>In recent years, high throughput computational screening (HTCS) has emerged as an efficient method (relative to experimental screening) for the exploration of the MOF-space. A typical HTCS study consists of the computational synthesis of a database of MOF prototypes, followed by molecular simulation(s) to calculate MOF performance-defining properties for a target application. Ideally, some of the promising MOF prototypes in the database would then be selected for synthesis and experimental validation of their properties, thus accelerating material discovery. Indeed, a few highperforming MOFs have been synthesized as the result of an HTCS study. For example, the chemically robust NU-800 for methane storage 16 or of the small-pore SBMOF-1 for Xe/Kr separation. 20 However, many HTCS studies do not result in the synthesis of a high-performing MOF, which is perhaps due to the way databases of MOF prototypes are constructed.</p><p>MOFs in these databases are built by connecting building blocks primarily on the basis of geometrical rules informed by coordination chemistry. [21][22][23][24][25] Thus, although the majority of MOF prototypes in databases may look feasible, there is uncertainty as to whether a given prototype can actually be synthesized. Currently, synthesis decisions following HTCS are based on chemical intuition about MOF synthesizability. Given the cost of exploratory experiments, synthesis is not necessarily attempted for the best MOF prototypes. Rather, synthesis is preferably attempted for "good" prototypes that resemble previously synthesized MOFs (e.g. a prototype within a known isoreticular MOF series), at the expense of potentially missing truly exceptional MOFs. This hinders the effectiveness of simulation-guided MOF discovery, and effectively constrains the exploration of the MOF-space to regions in the vicinity of already existing MOFs. Incidentally, the tendency to constrain exploration around already existing materials also pervades purely experimental work. As noted in our earlier work, despite the appeal of MOF tunability, MOF-space has arguably been underexplored, with 50% of over 40,000 articles on MOFs (published at the time) focusing on only 30 MOFs. 24 Clearly, one way to modify typical HTCS protocols to truly accelerate MOF discovery (and encourage wider experimental exploration of MOF-space) is to predict the "synthesizability" of MOF prototypes in a database. The development of an effective criterion for synthesizability is tied to understanding how MOF synthesis occurs. However, although recent experiments have provided some insights into the fundamental processes controlling MOF nucleation and growth, 26 there remains a lack of consensus on exactly what factors control whether a MOF can be synthesized or not, and whether these factors are thermodynamic or kinetic in nature. Molecular simulations show potential for clarifying how and why certain MOFs are observed experimentally, and others are not. Yet, compared to the number of studies that focus on MOF property prediction, only a few simulation studies have investigated aspects relevant to MOF synthesis.  In 2004, Mellot-Draznieks et al. used Monte Carlo simulations to anneal building blocks (starting from a disordered arrangement) into MOFs. 27,28 Although the method was successful in predicting potential MOF structures for a given set of building blocks, few conclusions about MOF stability or formation mechanisms could be made due to the absence of dynamics and an unrealistic description of building block interactions. Since 2015, researchers have started to use molecular dynamics (MD) simulations to follow kinetic aspects of MOF self-assembly. Starting with disordered building blocks in an implicit solvent, Yoneya et al.  where able to tune metal-linker interactions during a long MD simulation until formation of crystalline framework was accomplished. 29 However, the effects of artificially tuning the building block interactions on the formation mechanism remained unclear. Biswal and Kusalik used equilibrium MD to observe the early formation stages of MOF-2/MOF-3 in explicit solvent. However, only a primarily amorphous phase was able to form during the simulation time. 30,31 Kollias et al. 32 used metadynamics simulations to probe the free-energy landscape during the formation of crystal-like secondary building units (SBUs). Notably, Colón et al. 33 used MD simulations and the finite string method to probe the free-energy landscape during the assembly of 27 SBUs into a ~2.5 nm MOF-5 crystallite. 33 Evidently, the use of biased MD to model MOF self-assembly is encouraging, but questions remain on how the choice of interaction parameters, reaction coordinates, and model simplifications could impact conclusions.</p><p>Another set of simulation studies have focused on interrogating the energetics of fully assembled MOFs. Many of these studies have focused on understanding the synthetic preference for certain isomorphs. MOF isomorphs, while all having the same composition, can have considerably different properties 16,24,34 due to differences in how the building blocks are arranged in space. That is, due to differences in crystal topology. [35][36][37] Accordingly, when a high-performing MOF prototype is identified during HTCS, there is a chance that a lowperforming isomorph is synthetically preferred. Using molecular mechanics, Amirjalayer et al. showed that the experimentally observed HKUST-1/Cu-BTC 38 (tbo topology) has a lower strain energies than its unobserved pto isomorph. Analogous findings were reported for the experimentally observed MOF-14 39 (of the catenated, pto-c topology) and its unobserved tbo-c isomorph. 40 Impeng et al. showed that three experimentally observed Cu-paddlewheel MOFs of the nbo-b topology had lower strain energy than four hypothesized isomorphs, respectively. 41 Using density functional theory (DFT) energies, Cai et al. 42 explained the preference between the tbo and fmj topology for four MOFs (one of them being Cu-BTC) depending on how the BTC linker was functionalized. Arhangelskis et al. used DFT energies to determine the most stable isomorph among eight zeolitic imidazole framework (ZIF) prototypes based on a CF3IM linker, which was subsequently synthesized. 43 Based on the similarity of DFT formation energies, Liu and Truhlar explained the often reported coexistence of the MOF isomorphs NU-901 and NU-1000 (scu and csq topology, respectively). 44 Curiously, an inspection of the literature [45][46][47] seem to indicate that synthesis more readily yields the reportedly less stable isomorph, NU-1000, which could be due to a number of synthetic factors not being considered in the calculations.</p><p>Although strain/lattice energies have shown some success in explaining synthetic preference for a given MOF isomorph, [40][41][42][43] these quantities do not account for thermal effects. On the other hand, there are instances where temperature (while keeping other synthesis conditions constant) has been reported to change the topology of the synthesized MOF. [48][49][50] Thus, in principle, a temperature-dependent quantity such as free energy should be considered in MOF synthetic accessibility predictions. Gee and Sholl 51 calculated the harmonic free energy and free energy of solvation for 25 ZIF isomorphs to assess their relative stability. However, large-scale, free energy calculations on MOF databases used for HTCS have not been attempted despite the potential insights that could arise from the large amount of simulation data (e.g. through the emergence of well-defined structure-property relationships relevant to MOF synthesizability).</p><p>Here we decided to assess the viability of classical MD simulations to estimate MOF free energies in a way that is (sufficiently) accurate, yet still amenable for largescale screening. Inspired by methodologies used in the pharmaceuticals field for crystal structure prediction of organic molecules, [52][53][54][55][56] we decided to apply and evaluate two methods to calculate free energies in a database of 8,500 MOF prototypes: the quasi-harmonic approximation (QHA) and Figure 1. The computational building blocks used to construct our 8,500-MOF database. a) nodular inorganic building blocks (NIBBs). b) 3-dimensional, nodular organic building blocks (NOIBs). c) 2-dimensional nodular organic building blocks. d) connecting building blocks (CBBs). NOIBs approximately correspond to so-called "secondary building units" (SBUs) in MOF synthesis literature. NOIBs correspond to the node of multitopic organic linkers. CBBs approximately correspond to ditopic linkers and/or to the arms of multitopic organic linkers. thermodynamic integration, specifically the Frenkel-Ladd (FL) path. [57][58][59] The MOF prototypes were computationally constructed using ToBaCCo-3.0, 24,25 using 34 building blocks combinatorically assembled according to 110 topological blueprints. As is typical of this procedure, 17,21,[60][61][62] one expects most generated MOFs to not have been experimentally reportedregardless of whether they are actually synthesizable or not-but one also expects to generate some experimentally reported MOFs. Additionally, with sufficient topological diversity in the database, one expects to find numerous, non-experimentally observed isomorphs of the previously synthesized MOFs. Here we aimed to exploit this composition of the database to address questions beyond the methodology to calculate MOF free energies. Specifically, using the free energy of previously synthesized MOFs as a point of reference, we investigated the efficacy of free energy as a criterion to predict i) the absolute synthesizability of a MOF of any composition, and/or ii) the relative stability between isomorphic MOFs.</p><!><p>2.1 MOF Database Construction. Our 8,500-MOF database was created using ToBaCCo-3.0, 24,25 with the building blocks shown in Fig. 1. These building blocks were selected to allow for topological and chemical diversity, allowing our database to span a wide range of textural and energetic properties. In contrast to databases we created for recent studies focused on the prediction of adsorption properties, 61,63 the current database does not include functionalized MOFs, which allowed us to explore more linker shapes and coordination schemes, and enhance the occurrence of isomorphism. The geometry of each built MOF prototype was optimized in LAMMPS 64 with the interatomic interactions being described by the relevant force field: UFF, 65 Dreiding, 66 or UFF4MOF. 67,68 Note, however, that since parameters for Cu in UFF are defined for tetrahedral coordination, we altered the equilibrium angle to be 90° to be consistent with the Cu-paddlewheel geometry. Similarly, in Dreiding the equilibrium bond lengths and angles are unavailable for some metals (Cu, Cr, Zr), so we adopted the equilibrium bonds/angles we used in our UFF implementation but otherwise retained the functional forms and force constants used for Dreiding. Note that these modifications to Dreiding and UFF have been previously validated for MOFs. 69 To optimize MOF geometries, we employed an iterative approach. First, only the atom positions were relaxed using the fast inertial relaxation engine (FIRE) algorithm 70 with a timestep of 10.0 fs. Then the atom positions and simulation cell parameters were relaxed simultaneously using a conjugate gradient algorithm. We repeated this process (relaxing only the atom positions and then the atom positions and cell parameters together) until the energy change from the previous iteration was less than 4.184 × 10 kJ/mol. Each individual optimization finished when the maximum force on any atom was less than 4.184 × 10 kJ/(mol Å) and the energy changed by less than 1 × 10 % from the previous step.</p><p>2.2 Free Energy Calculations. All free energy calculations were performed using molecular dynamics (MD) in LAMMPS with a timestep of 1 fs. We calculated quasi-harmonic approximation (QHA) free energies using:</p><p>max where 𝑈 is the potential energy, 𝜈 are vibrational frequencies, and 𝐷(𝜈) is the vibrational (phonon) density of states (PDOS). [71][72][73] 𝐷(𝜈) is normalized such that ∫ 𝐷(𝜈)d𝑣 = 3𝑁 where 𝑁 is the number of atoms in the sytem. 71 To apply Eq. 1 we needed to calculate the PDOS for each MOF, which we did by Fourier-transforming the velocity autocorrelation function (VACF). 72,73 Each VACF was computed from a short (10 ps) NPT ensemble MD simulation with a Nosé-Hoover thermostat (damping parameter of 100 timesteps) and barostat (damping parameter of 1000 timesteps). For these NPT simulations we used a temperature of 300 K and zero pressure to allow direct comparison with Frenkel-Ladd path free energies. First, the simulation supercell was equilibrated for 0.1 ns, then we collected the velocities of each MOF atom every 5 timesteps (5 fs) for 10,000 timesteps (0.01 ns). We used this velocity trajectory to compute the VACF for each MOF.</p><p>We calculated the Frenkel-Ladd (FL) path free energies using the method of Freitas et al. which is described in detail in ref. 59. Briefly, the FL path is a thermodynamic integration method in which the "real" Hamiltonian (where interatomic interactions are described by the relevant force field) is switched to the reference Einstein crystal Hamiltonian (where each atom is constrained to its defined equilibrium position by a harmonic potential) using a parametrical Hamiltonian:</p><p>where 𝐻 is the initial (or real) Hamiltonian, 𝐻 is the final (or Einstein crystal Hamiltonian), and 𝜆 is the switching parameter. For 𝜆 = 1, 𝐻(𝜆) is 𝐻 and for 𝜆 = 0, 𝐻(𝜆) is 𝐻 . 57,59,74 The difference in free energy (Δ𝐹) between the Einstein crystal representation and the force field representation of the MOF is equal to the work required to switch to the Einstein crystal, which is given by:</p><p>where angled brackets denote an ensemble average at each value 𝜆 . Since the free energy of an Einstein crystal (𝐹 ) can be computed analytically, an absolute free energy of the solid is given by:</p><p>) where:</p><p>with 𝜔 being the harmonic frequency of atom i. One way to calculate the Δ𝐹 defined by Eq. 3 is by performing a series of equilibrium MD simulations for each value of 𝜆, during which 〈𝜕𝐻/𝜕𝜆〉 is calculated as an ensemble average. These averages can then be numerically integrated from 𝜆 = 1 to 𝜆 = 0 to yield Δ𝐹. However, Freitas et al. showed (following work by de Koning 75 ) that Δ𝐹 can be calculated as the average work of two nonequilibrium switching processes. 59 First, the system is switched from the real Hamiltonian to the Einstein crystal Hamiltonian (i.e. 𝜆 is changed from 0 to 1) during the course of a single non-equilibrium MD trajectory. Then, the system is switched back to the real Hamiltonian after equilibration with the Einstein crystal Hamiltonian. The forward work is denoted 𝑊 → irr and the reverse work is denoted 𝑊 → irr , with the superscript indicating that the forward/reverse switching processes are irreversible. The key ideas are that, for sufficiently long switching trajectories, the heat dissipated during switching is small 76 and that the heat dissipated during the forward and reverse processes are the same and can be cancelled (see ref. 59 for a detailed justification), meaning:</p><p>Accordingly, from Eq. 4-6 one can write:</p><p>where the third term on the right hand side of Eq. 7 is included to correct for the required fixing of the system center of mass during the simulation. 59,74 In principle, 𝑊 → and 𝑊 → irr should be averages calculated from several replicate simulations. However, we found that a single simulation with long enough switching time was sufficiently accurate (see Section S1).</p><p>Practically, we implemented the FL path free energy calculations in a single MD simulation. First, the simulation supercell was equilibrated for 0.1 ns in the NPT ensemble (Nosé-Hoover thermostat with damping of 100 timesteps and barostat with damping of 1000 timesteps), with the average cell parameters being calculated over a subsequent 0.1 ns of simulation time. We then switched to the NVT ensemble (as the FL path is only valid for the canonical ensemble) using the average cell parameters calculated in the previous step. During this simulation stage we used the thermostat of Bussi et al. 77 (damping parameter of 100 timesteps) which can correctly treat harmonic oscillators. 58 We then calculated the mean-squared displacement of each forcefield atom type over 0.1 ns after an equilibration period of 0.03 ns (see below). We then ran the forward switching process for 0.4 ns (see details on switching time convergence in Section S1 and Figure S1). After forward switching, the system was allowed to equilibrate according to the Einstein crystal Hamiltonian for 0.03 ns, followed by the reverse switching process for 0.4 ns. The equilibrium positions used in the Einstein crystal representation where calculated by relaxing the atom positions according to the real Hamiltonian and using the cell parameters calculated from the NPT trajectory. The above procedure, including forward and reverse switching and Einstein crystal equilibration was implemented using the LAMMPS command fix ti/spring. 59 The "strength" of the springs for each atom in the Einstein crystal is defined by force constants. Here, we calculated a different force constant (𝛼 ) for each atom type i in the simulation (according to the relevant force field), using:</p><p>where 〈(∆𝒓) 〉 is the mean squared displacement (MSD) of atom type i (calculated as described above). This ensures that the MSDs for the real and Einstein crystals are similar, which makes the switching process more efficient. 59 We, however, enforced a lower bound on 𝛼 of 41.84 kJ/(mol Å 2 ) (see detailed justification in Section S1) in cases where 𝛼 was too small, as particularly small values of 𝛼 allow for dramatic displacements from the defined Einstein crystal equilibrium positions and lead to LAMMPS errors.</p><!><p>3.1 Comparison of MOF free energy calculations methods. The quasiharmonic approximation (QHA) assumes that each MOF vibrational mode can be reasonably modeled with a harmonic potential. By contrast, no such assumption is made in the Frenkel-Ladd (FL) path method. In principle, then, one can take the FL free energies to be the more accurate free energies. However, should QHA and FL free energies differ by a sufficiently small amount, one would prefer to calculate QHA free energies in MOF databases as they require shorter simulations and are technically easier to implement than the FL path (Fig. S3). Thus, we start our discussion by comparing the QHA and FL path free energies (both calculated using UFF4MOF) in our 8,500-MOF database. Fig. 2 is a parity plot comparing QHA (y-axis) and FL (x-axis) free energies. All free energies in the present work are calculated at 300 K and reported on a "per atom" basis unless specified otherwise. If the energies calculated using the two methods were equal, all points would lie on the line 𝑦 = 𝑥 (the parity line). Based on Fig. 2, there is systematic deviation between FL path and QHA free energies, which is consistent with observations by Pastorino and Gamba when comparing QHA and FL free energies in sulfur crystals. 72 Purely systematic deviation in QHA free energies would not matter when using free energy as a metric of relative thermodynamic stability across MOFs in a database, as free energy differences between individual structures would be exactly maintained. Such a scenario would be reflected by a perfectly linear correlation between QHA and FL free energies.</p><p>To quantify the correlation between QHA and FL free energies, we report the Pearson (R) and Spearman rank (S) correlation coefficients. An 𝑅 value of one indicates a perfectly linear correlation between QHA and FL free energies, whereas an S value of one indicates that ordering MOFs by relative stability based on QHA or FL path free energies would result on identical rankings. We found both R and S to be close to one (Fig. 2, 𝑅 = 0.996, S = 0.939), which indicates a strong correlation between QHA and FL free energies. These R and S values suggest that one would observe similar free energy vs. MOF property relationships with either the QHA or the FL method (i.e. broad trends would be preserved when switching between methods). However, the reproduction of broad trends has no bearing on whether QHA and FL path free energies are equally valid criteria for assessing the relative thermodynamic stability of MOFs. To further quantify differences between QHA and FL free energies, we calculated the errors between QHA free energies and the corresponding FL free energies, after removing the systematic deviation from the QHA values, which we accomplished by subtracting the y-intercept of a linear model fit to the data shown in Fig. 2. Once this adjustment was made, we found the absolute errors to span the 0-10 kJ/mol per atom range (Fig. S4), with the mean absolute error (MAE) being 1.3 kJ/mol per atom. We will see in Sections 3.3 and 3.4 that i) previously synthesized MOFs encountered in our database cluster below a 4.4 kJ/mol per atom free energy threshold, and ii) free energy differences between MOF isomorphs are typically small (with a mean difference of 3 kJ/mol and many differences of less than 1 kJ/mol). On these grounds, one can see that the use of QHA free energies to assess thermodynamic stability (and perhaps synthesizability) of MOFs may be problematic. Based on the above, in subsequent sections we will focus our discussion on FL free energies. However, at this point let us discuss the origin of the differences between QHA and FL free energies in MOFs. Considering the underlying QHA assumptions, one would expect that QHA free energies would be accurate for a MOF with only high-frequency (harmonic) vibrations, and inaccurate for a MOF with only low-frequency (anharmonic) vibrations. However, as illustrated by their phonon density of states (PDOS) (Fig. 3a and Fig. S5), MOFs present vibrations in a wide range of frequencies from low to high, with many vibrational modes in the 1000-2000 cm -1 range. Although all MOFs seem to present a single, large peak at ~3000 cm -1 in their PDOS, some MOFs do have a prominently larger number of vibrations at lower frequencies (< 1000 cm -1 ), including well-known MOFs such as UiO-66 (Fig. 3a) and IRMOF-1 (Fig. S4).</p><p>Thus, to elucidate which frequencies tended to collectively contribute the most to the QHA free energies of our MOFs, we split the PDOS for each MOF into five equal ranges and calculated the contribution of each region to the QHA vibrational free energy. This contribution depends on both the individual contribution of each vibrational mode and the number of vibrational modes in each frequency range. Knowing the percent contribution of each range to the total vibrational free energy for each MOF, we constructed a histogram (Fig. 3b) that reveals that, in most MOFs, high-frequency vibrations (> 2669 cm -1 ) contribute the most to the total QHA free energy, followed by 1334-2001 cm -1 and 667-1334 cm -1 vibrations, respectively. The lowest frequency vibrations (< 667 cm -1 ), on the other hand, contributed negligibly to the vibrational free energy of the MOFs in our database.</p><p>However, if we combine the 1334-2001 cm -1 and 667-1334 cm -1 distributions shown in Fig. 3b the combined distribution covers approximately the same range of contributions to QHA free energy as the high-frequency distribution. We refer to MOFs in the 90 th percentile of this combined distribution as "flexible" MOFs, in contrast to "stiff" MOFs, which are in the 90 th percentile of the high frequency distribution (see the full details of this categorization in Section S2). When differentiating between flexible and stiff MOFs on the parity plot of QHA versus FL path free energies (Fig. 3c), some interesting trends emerge. First, the previously discussed systematic error for QHA free energies is higher for stiff MOFs. Yet, for these MOFs, the correlation between QHA and FL free energies is higher and the shifted absolute error (i.e. the absolute error after removing the systematic error) lower (Table S1). On the other hand, while the systematic error is lower for flexible MOFs, the correlation between QHA and FL free energies is worse and the shifted absolute error higher (Table S1). From this we can conclude that the systematic error for QHA free energies is due to high-frequency vibrations, while the non-systematic error is due to mid-frequency vibrations.</p><!><p>A question that emerges in any classical simulation study is the sensitivity of calculations to force field selection. Up to this point, we have focused on free energies calculated using UFF4MOF, because this update of UFF by Addicoat et al. 67,68 was specifically made for MOFs. While additional MOF force fields (of varied generalizability) have emerged, [78][79][80] here we focus on a comparison of UFF4MOF with the generic force fields Dreiding 66 and UFF. 81 This was motivated by the tendency of HCTS studies to use these generic force fields (given their complete transferability), as well as work by Boyd et al. 69 who reported that Dreiding and UFF perform reasonably well to predict MOF properties such as bulk modulus and thermal expansion coefficients.</p><p>Parity plots comparing FL free energies among force fields are presented in Fig. 4. The closest agreement in free energies is between UFF and UFF4MOF, especially for MOFs with Cu (paddlewheel) nodes (Fig. S6). The main deviations are observed for MOFs with Cr, Zn, and especially Zr nodes, which can be rationalized by considering that UFF4MOF only modifies the UFF parameters that involve metals (e.g. M-O bonds, O-M-O, M-O-C angles, etc.), 67,68 and there is a large fraction of these bond/angles in our Cr, Zn, and Zr-MOFs. The largest disagreement in free energies is between Dreiding and UFF4MOF, which is expected, given that Dreiding uses different functional forms and force constants than UFF4MOF not only for interactions involving metals but also non-metals. Further inspection of outlier free energies in Dreiding reveal that they correspond to cases where Dreiding and UFF4MOF produce dramatically different MOF unit cell parameters.</p><p>Although the outliers in Fig. 4 visually stand out, a look to the density of points on the parity plots (Fig. S6) reveal a large number of points near the parity line. Thus, a more comprehensive comparison is given by the values of 𝑅, 𝑆 and MAE inset in Fig. 4, as well as histograms for the absolute errors (Fig. S7). In all cases R and S are close to one, indicating that similar free energy vs. MOF property relationships would be obtained with either force field. However, consistent with our discussion comparing QHA and FL free energies (MAE ~ 1 kJ/mol), we must conclude that the use of Dreiding (MAE ~ 3 kJ/mol) instead of UFF4MOF to scrutinize the thermodynamic stability of MOFs may be problematic. Although a similar conclusion regarding UFF is less clear (with a MAE of ~1 kJ/mol, but a median absolute error of 0.1 kJ/mol), it seems that the use of UFF4MOF to calculate MOF free energies is the best choice. Thus, from this point on our discussion are based only on calculation with this force field.</p><!><p>Inspired by prior studies that developed criteria for zeolite synthetic feasibility [82][83][84] we decided to explore similar criteria for MOFs. Such a criterion could be used as a preliminary step in HTCS studies to eliminate MOFs with low synthetic likelihood from a computational database before running simulations to predict relevant performance metrics. Thus, computer resources would be dedicated only to MOFs that have some chance to be synthesized following HTCS. Moreover, even when calculations on non-synthesizable MOFs are still valuable (e.g. for machine learning 63,85 ) a synthetic likelihood criterion could be used to downsize a database that is too large for simulations of a given computational cost.</p><p>Given that MOFs have diverse compositions, the development of a synthetic likelihood criterion is perhaps less straightforward than for zeolites. For example, while Li et al. identified a criterion for zeolite stability based on the relationship between T-site to T-site distances and Osite to O-site distances, 82 such a relationship was only possible because all zeolites are arrangements of TO4 tetrahedra. An analogous criterion certainly does not exist for MOFs. On the other hand, Bushuev et al. 84 and Dawson et al. 83 identified regions of synthetic feasibility in zeolites when plotting the energy of existing and hypothesized zeolites vs. zeolite geometric properties. They found that previously synthesized zeolites tended to cluster in particular regions of energy versus density plots. It is easier to imagine an extension of this approach to MOFs, so we plotted our MOF free energies against MOF properties to see whether previously synthesized MOFs tended to cluster in a specific region.</p><p>First, we needed to identify whether there were MOFs in our database that had been previously synthesized. To accomplish this, we calculated the MOFid (a systematic MOF identifier developed by Bucior et al. 86 ) for each MOF our database. We then compared the MOFids in our database to those of the CoRE MOFs (a database of previously synthesized MOFs). 87,88 Since the MOFid for a specific MOF is unique, 86 if a MOFid for a MOF in our database was also present in the CoRE MOFs, that means that the MOF in our database has been previously synthesized. Through this procedure, we found 54 previously synthesized MOFs in our database. For comparison, consider that Bucior et al. 86 found 50 and 18 previously synthesized MOFs from 10,159-and 1614-MOF subsets of earlier computationally constructed databases by Gómez-Gualdrón and coworkers 17,25 and Wilmer and coworkers. 21 After identifying the previously synthesized MOFs in our database, we proceeded to plot FL free energies against MOF properties (Fig. 5a-c). We found that previously synthesized MOFs do not cluster in a single region. Rather, we found that they cluster in different regions depending on their composition. For instance, Zr MOF free energies occupy higher free energy regions than the other MOFs. We found that this occupation of high free energy regions (Zr-MOFs > Cr-MOFs > Zn-MOFs ~ Cu-MOFs) correlates with the (strain) energy (per atom) of the corresponding isolated nodes (Table S2). The relative energies of the isolated nodes also explains the different behavior of free energy versus the metal/linker atom ratio for the MOF compositions-in calculating this ratio, any Cu, Cr, Zn, or Zr atom was counted as a metal, and any atom bonded to a carbon was counted as a linker atom. For instance, in Zr-MOFs, free energy increases pronouncedly with the metal/linker atom ratio as adding more Zr-nodes will always increase potential energy.</p><p>When examining the different regions separately, we can observe that the free energies of previously synthesized MOFs tend towards the lower (free energy) boundaries of the corresponding region. This supports the hypothesis that having a relatively low free energy-at least relative to same MOF types-is key to MOF synthetic feasibility. However, developing separate criteria for synthetic feasibility for different MOF types is inconvenient and less informative than a single criterion. The clearest pattern for previously synthesized MOFs emerges from Fig. 5c, where a linear relationship (for each MOF type) between free energy and metal /linker atom ratio is apparent for each MOF type, allowing us to fit four separate linear models to these free energies, as illustrated in Fig. 5d-g. Other synthetic factors aside, it is reasonable to assume that MOFs below the corresponding fit line have high synthetic likelihood, while those above have low synthetic likelihood.</p><p>The data shown in Fig. 5c can be transformed by subtracting the free energy output by the relevant linear model from each MOF original free energy. We refer to the free energy transformed in this manner as ∆ LM 𝐹 FL . The results of this transformation are presented in Fig. 6, where all the previously synthesized MOFs in our database now cluster within a small ∆ LM 𝐹 FL region. The maximum ∆ LM 𝐹 FL among experimentally observed MOFs in our database is 4.4 kJ/mol. If we classify all MOFs with ∆ LM 𝐹 FL higher than this value as synthetically unlikely, we could eliminate 1,615 MOFs from our database. Or a different number if some tolerance is used (e.g. 1,354 MOFs, if we use a 20% tolerance for the ∆ LM 𝐹 FL cutoff). Clearly, the exact ∆ LM 𝐹 FL criterion for synthetic feasibility presented here is specific to free energies calculated with UFF4MOF. However, an analogous criterion can be developed for another forcefield and extended to other MOF types using a similar strategy to the one outlined here. Note that we visually inspected many of the MOFs categorized as synthetically unlikely, and a large portion of these are not visibly strained (Figure S9). This indicates that visual inspection is likely insufficient to detect all synthetically unlikely MOFs, and that it is not obvious whether a geometry-based MOF synthesizability criterion could be developed.</p><!><p>Given the occurrence of isomorphism in MOFs, it is possible that a MOF that would be otherwise synthesizable (as per the criterion in Section 3.3), may not actually be experimentally observed simply because a more stable isomorph exists that is preferentially formed during synthesis. On the other hand, previous work 24 has shown that key MOF properties (e.g. mechanical stability, 24 adsorption selectivity, 24 and catalytic activity 34 ) can vary dramatically among isomorphic MOFs. Therefore, to inform synthetic attempts following HTCS, it is important to predict the most synthetically likely MOF within every isomorphic family in a database. Additionally, by only retaining the synthetically likely isomorphs one could further downsize a computational MOF database.</p><p>So, is it possible to predict the synthetically preferred isomorph? To address this question, we first determined which MOFs in our databases had isomorphs. To accomplish this, we identified all the unique linkers in our database by their canonical SMILES string, and then considered any MOFs with the same node and linker(s) to be isomorphs. We found that the free energy differences among isomorphs can frequently be well under 1 kJ/mol per atom (Fig. 7), which is less than the width of the "synthetic likelihood" region in Fig. 6, indicating that the previously discussed synthesizability criterion discussed may not be sufficient to eliminate all synthetically inaccessible isomorphs (although it could be useful as a first filter).</p><p>To directly assess whether free energies can accurately differentiate between synthetically accessible and non-accessible isomorphs, we identified all the isomorphs for each previously synthesized MOF in our database. This allowed us to identify 32 isomorphic series in our database that also included at least one previously synthesized MOF (20 previously synthesized MOFs in our database did not have any isomorph). Two isomorphic series contained two previously synthesized isomorphs each and are discussed in Section S5, along with additional case-studies. As for the 32 isomorphic series, we determined whether the experimentally observed isomorph in each series indeed had the lowest free energy. This was the case for 27 (84%) of the isomorphic series (Table S4), or 25 (83%) if we do not consider the series with two experimental MOFs. This accuracy is remarkable considering the small free energy differences between isomorphs. We show one of the 32 studied isomorphic series in Figure 8, where the free energy of the synthetically preferred NOTT-100 (fof topology) 89 is less than 0.5 kJ/mol per atom lower than three other possible isomorphs. Although for five isomorphic series the experimentally observed member did not have the lowest free energy, the free energy of the experimentally isomorph was within (at most) 1.7% of the lowest one. The infrequent discrepancies between free energy rankings and the experimentally reported isomorph may thus be due to i) synthetic factors (e.g. solvents, modulators, and kinetics in general) not accounted for by the free energy calculation, ii) experimental oversight where the lowest free energy isomorph has formed during synthesis but not reported-note that isomorphs that are close in free energy may coexist such as the NU-1000 and NU-901 cases discussed in Section S5-or iii) simply a prediction error inherent to either the forcefield or our free energy calculation procedure. Based on ii and iii, one may consider the use of a tolerance before deciding to discard a MOF isomorph as synthetically inaccessible. 3.5. Can relative MOF stability be assessed with alternatives to free energy? MOF free energies have shown considerable potential to be used as metric for synthetic likelihood of MOF prototypes. Thus, free energies could be used to eliminate synthetically unlikely MOFs from databases (reducing their size and increasing the computational efficiency of HTCS), and to inform synthesis efforts. However, since free energies would need to be calculated for all MOFs in a database, an important question to ask is whether more affordable metrics could be used as effectively. Although we have not envisioned a geometry-based metric at this point, we decided to briefly inspect whether metrics such as potential energy at 0 K (strain energy, U0K) and at the temperature of interest (here, U300K) could be as effective as the FL free energy as a metric of synthetic likelihood. This was partly motivated by the strong correlation between free, potential and strain energy (Table 1) with MAE values smaller than 0.8 kJ/mol (although these correlations may worsen for when considering higher temperatures). For the 300 K temperature considered here, it turns out that using potential or strain energy in an analogous manner to how we used free energy data in Fig. 5 to obtain the plot in Fig. 6, yields similar results. In other words, synthetic likelihood regions for MOFs can also be identified on the basis of strain or potential energy (Fig. S15), which could be used to discard a roughly similar number of MOFs deemed synthetically unlikely. However, as noted earlier, a more challenging task (since smaller free energy differences are involved) is to be able to discard otherwise synthetically likely MOFs that present a more thermodynamically stable isomorph.</p><p>Figure 9. a) Percentage of MOFs that can be discarded because they are synthetically inaccessible as a function of tolerance, the dashed line corresponds to the percent of MOFs discarded at the tolerance for free energy (1.7%) that retains all of the experimentally observed isomorphs in our database. b) Percentage of synthetically accessible isomorphs that are kept as a function of tolerance, the 100% dashed line corresponds to 34 isomorphs. The tolerance is defined as being within a certain percent of the lowest free, potential or strain energy predicted for an isomorphic series to be considered synthetic accessible.</p><p>From Section 3.4 we already know that on the basis of having the lowest free energy, we correctly predicted the experimentally observed isomorph for 27 (out of 32) isomorphic series to be synthetically accessible, and needed a 1.7% tolerance to classify all experimentally observed isomorphs as synthetically accessible (to the expense of a few false positives). If we were to use potential or strain energy as a proxy for free energy, we would correctly predict the experimentally observed isomorph for 25 and 22 isomorphic series, respectively. Thus, using potential (strain) energy results in two (four) more false negatives than free energy. A tolerance of 3.3% is needed to classify all the experimentally observed isomorphs in our database as synthetically accessible using potential energy, slightly higher than the tolerance required using free energy. However, a significantly higher tolerance of 93.0% is required when using strain energy. Thus, using potential or (especially) strain energy to identify synthetically accessible isomorphs would result in more false positives.</p><p>It would be possible to lower the tolerance when using strain or potential energy as a proxy for free energy, if one were to accept the occurrence of false negatives. As the balance between false negatives and false positives is a recurring subject in efforts to accelerate HTCS, we calculated the percentage of synthetically accessible isomorphs that would be kept in our database subset (composed by 32 isomorphic series) as a function of tolerance (Fig. 9b), as well as the percentage of MOFs that can be discarded due to synthetic inaccessibility (Fig. 9a). A loose tolerance would result in fewer MOFs being discarded, but more false positives are bound to occur. A tight tolerance allows for more MOFs to be discarded, but more false negatives are bound to occur. At the tolerance required to retain all the synthetically accessible isomorphs in our database (i.e. an estimate of the tolerance required to eliminate false negatives), the free energy metric enables the elimination of the largest number of MOFs (i.e. minimizes the number of false positives).</p><p>Based on the above observations, we conclude that free and potential energy definitively outperform strain energy for differentiation between synthetically accessible and inaccessible MOF isomorphs. Both free and potential energy require low tolerances (1.7 and 3.3%, respectively, to remove false negatives in our data), and would likely perform similarly in the vast majority of cases, at least for the 300 K temperature considered here. However, since there are instances where we found the experimentally observed isomorph to have the lowest free energy but not the lowest potential energy-and never the other way around-free energy may be safer metric to use at this point. Specially since we have not studied the correlation between free and potential energies at higher temperatures. However, in instances where isomorphism is not a factor, even the MOF location with respect to a synthetic likelihood region delimited by a strain energy threshold (Fig. S15a) may be acceptable to classify a MOF (at least preliminarily) as synthetically likely or unlikely.</p><p>As a final point discussing the effectiveness of free, potential and strain energies to characterize synthetic likelihood, let us note we identified the lowest free energy isomorph within each of the 960 isomorphic series we found in our 8,500-MOF database. Only 781 (81%) and 744 (77%) were also the lowest potential and strain energy isomorph, respectively. Consequently, if we consider free energy to be the most accurate identifier of a synthetically accessible isomorphs, using potential or strain energy as a proxy for free energy may lead to incorrectly missing the synthetic feasibility of an unacceptably large number of MOFs, at least if no tolerance is used.</p>