We have studied the clusters involved in the initial stages of nucleation of Zeolitic Imidazolate Frameworks, employing a wide range of computational techniques. In the pre-nucleating solution, the prevalent cluster is the ZnIm 4 cluster (formed by a zinc cation, Zn 2+ , and four imidazolate anions, Im -), although clusters such as ZnIm 3 , Zn 2 Im 5 , Zn 2 Im 7 , Zn 3 Im 9 , Zn 3 Im 10 , or Zn 4 Im 12 have energies that are not much higher, so they would also be present in solution at appreciable quantities. All these species, except ZnIm 3 , have a tetrahedrally coordinated Zn 2+ cation. Small Znx Imy clusters are less stable than the ZnIm 4 cluster. The first cluster that is found to be more stable than ZnIm 4 is the Zn 41 Im 88 cluster, which is a disordered cluster with glassy structure. Bulk-like clusters do not begin to be more stable than glassy clusters until much larger sizes, since the larger cluster we have studied (Zn 144 Im 288 ) is still less stable than the glassy Zn 41 Im 88 cluster, suggesting that Ostwald's rule (the less stable polymorph crystallizes first) could be fulfilled, not for kinetic, but for thermodynamic reasons. Our results suggest that the first clusters formed in the nucleation process would be glassy clusters, which then undergo transformation to any of the various crystal structures possible, depending on the kinetic routes provided by the synthesis conditions. Our study helps elucidate the way in which the various species present in solution interact, leading to nucleation and crystal growth.

Introduction

Nucleation is the first stage in the formation of solid-state matter and its control constitutes a fundamental cornerstone for the tailored design and synthesis of functional materials. Understanding the chemistry at play in nucleation processes is far from straightforward, even for simple systems, and benefits from microscopic insights and principles scalable from clusters to crystallites with seminal bulk properties. [START_REF] Schüth | Nucleation and crystallization of solids from solution[END_REF][START_REF] Jin | Recent experimental explorations of non-classical nucleation[END_REF] Advances in the observation of molecular nucleation processes are intensively pursued to bridge the length and time scales involved. Typically, a recent study has captured images of the nucleation processes of the benchmark NaCl system. [START_REF] Shoaee | Attfield. Crystal growth of the nanoporous metal-organic framework HKUST-1 revealed by in situ atomic force microscopy. Angew[END_REF] Similar experiments for the Metal-Organic Frameworks (MOFs) are extremely challenging due to the diversity of the pool of organic and metal-based species from which self-assembly occurs, and the many types of coordinative bonds and polymorphic structures. Different experimental studies can even lead to contradictory conclusions regarding the role of prenucleation and secondary building (SBUs) or growth units in the pathways of MOF assembly. [START_REF] Shoaee | Attfield. Crystal growth of the nanoporous metal-organic framework HKUST-1 revealed by in situ atomic force microscopy. Angew[END_REF][START_REF] Embrechts | Elucidation of the formation mechanism of metal-organic frameworks via in-situ raman and FTIR spectroscopy under solvothermal conditions[END_REF][START_REF] Cubillas | Materials discovery and crystal growth of zeolite A type zeolitic-imidazolate frameworks revealed by atomic force microscopy[END_REF][START_REF] Cubillas | Crystal growth mechanisms and morphological control of the prototypical metal-organic framework MOF-5 revealed by atomic force microscopy[END_REF] The vast majority of the theoretical and computational research in the field of MOFs has been devoted to adsorption, diffusion and separa-tion processes, [START_REF] Ongari | Too many materials and too many applications: An experimental problem waiting for a computational solution[END_REF][START_REF] Demuynck | Efficient construction of free energy profiles of breathing metalorganic frameworks using advanced molecular dynamics simulations[END_REF] in a direct response to the urgent demand for applications of nanoporous materials, leaving less than a dozen works on nucleation. [START_REF] Cubillas | Crystal growth mechanisms and morphological control of the prototypical metal-organic framework MOF-5 revealed by atomic force microscopy[END_REF][START_REF] Wells | Early stages of phase selection in MOF formation observed in molecular monte carlo simulations[END_REF][START_REF] Yoneya | Simulation of metal-organic framework self-assembly[END_REF][START_REF] Kollias | Molecular level understanding of the free energy landscape in early stages of metal-organic framework nucleation[END_REF][START_REF] Biswal | Probing molecular mechanisms of self-assembly in metal-organic frameworks[END_REF][START_REF] Biswal | Molecular simulations of self-assembly processes in metal-organic frameworks: Model dependence[END_REF][START_REF] Filez | Elucidation of the pre-nucleation phase directing metalorganic framework formation[END_REF][START_REF] Yamil | Free energy of metal-organic framework self-assembly[END_REF][START_REF] Cheetham | Thermodynamic and kinetic effects in the crystallization of metal-organic frameworks[END_REF] A direct consequence of the lack of knowledge in terms of MOF crystal growth is that most studies devoted to the development of synthesis routes of new MOFs are based on chemical intuition and trial and error procedures rather than on a rational analysis of the nucleation and crystal growth processes. It is in fact remarkable that while over 70,000 MOFs have been reported up to date, [START_REF] Peyman | Development of a Cambridge Structural Database Subset: A collection of Metal-Organic Frameworks for past, present, and future[END_REF] such little focus has been put on the stages of nucleation or crystallization. [START_REF] Van | In situ, time-resolved, and mechanistic studies of metal-organic framework nucleation and growth[END_REF] Therefore, it is often stated, not without despair, that each MOF family has to be considered as unique, despite the general commonalities in coordination chemistry that the members of this large family of materials share. [START_REF] Maarten | Small-angle x-ray scattering documents the growth of metal-organic frameworks[END_REF] The atomistic modeling of nucleation and crystal growth of MOFs faces several major challenges. The non-periodic structure models needed to model nucleation preclude the use of large scale, first-principle calculations at a reasonable cost, although advances in the field have been reported in recent years. [START_REF] Wells | Early stages of phase selection in MOF formation observed in molecular monte carlo simulations[END_REF][START_REF] Yoneya | Simulation of metal-organic framework self-assembly[END_REF][START_REF] Kollias | Molecular level understanding of the free energy landscape in early stages of metal-organic framework nucleation[END_REF][START_REF] Biswal | Probing molecular mechanisms of self-assembly in metal-organic frameworks[END_REF][START_REF] Yamil | Free energy of metal-organic framework self-assembly[END_REF][START_REF] David | Formation mechanism of the secondary building unit in a chromium terephthalate metal-organic framework[END_REF][START_REF] Boussouf | Structure, spectroscopy, and bonding within the Zn q+ -Imidazolen (q = 0, 1, 2; n = 1-4) clusters and implications for Zeolitic Imidazolate Frameworks and Zn-enzymes[END_REF] Cost-effective classical force field simulations typically face the limitation of not being suitable in general to model bond breaking and formation, and of 1 Nucleation of zeolitic imidazolate frameworks: from molecules to nanoparticles including a limited ensemble of structural motifs in the parameterization of the interatomic potentials. In consequence, force field-based calculations can lead to different conclusions depending on the choice of the interaction model, and thus suffer from lack of transferability. [START_REF] Biswal | Molecular simulations of self-assembly processes in metal-organic frameworks: Model dependence[END_REF] The selection of force field parameters, including atomic point charges, can be difficult in the continuously changing chemical environment in which nucleation takes place. [START_REF] Jawahery | Ab initio flexible force field for metal-organic frameworks using dummy model coordination bonds[END_REF][START_REF] Nguyen | Microscopic origins of poor crystallinity in the synthesis of covalent organic framework COF-5[END_REF] The use of reactive force fields is regarded as a viable approach, as shown recently in a study of melting and recrystallization, [START_REF] Yang | Enabling computational design of ZIFs using ReaxFF[END_REF][START_REF] Salvador | Computer simulation of the early stages of self-assembly and thermal decomposition of ZIF-8[END_REF] but the accuracy of these methods requires substantial validation for both crystalline and pre-nucleating species, in order to ensure that the simulations are realistic. [START_REF] Macià Escatllar | Modeling hydroxylated nanosilica: Testing the performance of ReaxFF and FFSiOH force fields[END_REF] New avenues are expected from large-scale Monte Carlo (MC) and Molecular Dynamics (MD) simulations of MOFs with enhanced sampling techniques, following successes in modeling more straightforward systems. [START_REF] Yoneya | Simulation of metal-organic framework self-assembly[END_REF][START_REF] Kollias | Molecular level understanding of the free energy landscape in early stages of metal-organic framework nucleation[END_REF][START_REF] Biswal | Probing molecular mechanisms of self-assembly in metal-organic frameworks[END_REF][START_REF] Biswal | Molecular simulations of self-assembly processes in metal-organic frameworks: Model dependence[END_REF][START_REF] Cheetham | Thermodynamic and kinetic effects in the crystallization of metal-organic frameworks[END_REF] DFT based, ab-initio MD (AIMD) calculations, combining explicit solvent molecules and enhancing sampling techniques have also been employed to explore the competing mechanisms of formation of SBUs leading to MIL-101, MIL-53 and MOF-235. [START_REF] Maarten | Small-angle x-ray scattering documents the growth of metal-organic frameworks[END_REF][START_REF] David | Formation mechanism of the secondary building unit in a chromium terephthalate metal-organic framework[END_REF][START_REF] Yang | Enabling computational design of ZIFs using ReaxFF[END_REF] Morris 27 drew analogies between the nucleation and crystallization of MOFs and those of zeolites, leading to extensive research on the so-called Zeolitic Imidazolate Frameworks (ZIFs). However, the field is still in its infancy in terms of identifying the SBUs that drive the process of crystal growth. ZIF pre-nucleation building units with up to four metal atoms have been identified for ZIF-8, [START_REF] Ivy | Insights into the molecular assembly of zeolitic imidazolate frameworks by ESI-MS[END_REF] and CdIF-4, 29 by mass spectrometry techniques. Atomic Force Microscopy (AFM) was also employed for CdIF-4, [START_REF] Wagia | Determination of the preassembled nucleating units that are critical for the crystal growth of the metal-organic framework CdIF-4[END_REF] and it was the first time that metal-ligand species relevant for the surface nucleation and growth in a MOF were identified. The evolution of metal-organic aggregates of one or only a few Zn atoms and intermediate species, up to bulkier cluster units of ∼ 2 nm has also been evidenced in the formation of ZIF-8 from mass spectrometry and X-ray scattering and diffraction. [START_REF] Terban | Early stage structural development of prototypical zeolitic imidazolate framework (ZIF) in solution[END_REF][START_REF] Yeung | Control of metal-organic framework crystallization by metastable intermediate pre-equilibrium species[END_REF][START_REF] Cravillon | Rapid roomtemperature synthesis and characterization of nanocrystals of a prototypical zeolitic imidazolate framework[END_REF][START_REF] Buzanich | Observation of early ZIF-8 crystallization stages with x-ray absorption spectroscopy[END_REF] However, the structure of the larger units could not be correlated with those appearing in ZIF-8 and thus they were assigned to be amorphous precursors. In a remarkable attempt to lay out a systematic framework of MOF growth, a recent investigation by Filez et al. [START_REF] Filez | Elucidation of the pre-nucleation phase directing metalorganic framework formation[END_REF] has combined a variety of experimental methods and computations to link the prenucleation, nucleation and crystallization regimes in Co ZIF-67, by discriminating those metal-organic complexes involved in initiating nucleation and the species required for oligomerization into extended MOF networks. A complex framework of non-classic multi-stage nucleation routes is laid out in that investigation, which is expected to be largely driven by the rich metal-linker reaction chemistries concurring in the amalgam of species in the synthesis of this MOF. Very recently, Balestra and Semino 25 have studied the self-assembly of ZIF-8 via a force field that includes the possibility of breaking and forming metal-ligand bonds to achieve molecular details of the nucleation of ZIF-8.

In this work, we provide an unprecedentedly detailed analysis of the possible pre-nucleation building units for a large series of ZIFs, considering the influence of the solvent and of counterions on the pre-nucleation mechanisms. We conducted DFT calculations on a large ensemble of cluster configurations, based on Imidazolate (Im -) and 2-methyl-Imidazolate (mIm -) linkers, considering the presence of solvent and counterions to closely account for experimental conditions. mIm -is the linker in ZIF-8, the prototypic ZIF with cubic sodalite (SOD) topology. [START_REF] Cravillon | Rapid roomtemperature synthesis and characterization of nanocrystals of a prototypical zeolitic imidazolate framework[END_REF][START_REF] Kyo | Exceptional chemical and thermal stability of zeolitic imidazolate frameworks[END_REF][START_REF] He | Facile synthesis of zeolitic imidazolate framework-8 from a concentrated aqueous solution[END_REF][START_REF] Bustamante | Influence of the solvent in the synthesis of zeolitic imidazolate framework-8 (ZIF-8) nanocrystals at room temperature[END_REF] Remarkably, mIm -does not lead to the crystallization of ZIFs with other topologies. In contrast, ZIFs based on Im -have been prepared with a wide range of topologies, such as BCT, DFT, GIS, MER, zec, zni, AFI, and CAN. [START_REF] Kyo | Exceptional chemical and thermal stability of zeolitic imidazolate frameworks[END_REF][START_REF] Phan | Synthesis, structure, and carbon dioxide capture properties of zeolitic imidazolate frameworks[END_REF][START_REF] Shi | Zeolite CAN and AFI-type zeolitic imidazolate frameworks with large 12-membered ring pore openings synthesized using bulky amides as structure-directing agents[END_REF] Porous Zn-Im-ZIFs are typically obtained either with very low yields or with the aid of added templating molecules. We therefore study a large set of clusters using Zn and Im -/mIm -ligands, devoting efforts both to examine the pre-nucleation stages of formation of clusters of a few Zn atoms and to explore efficient methods to scale the computation to large clusters entering the nucleation regime, with tens to hundreds of Zn atoms.

Computational details

In this study three length scales are treated for studying ZIF nucleation: smaller clusters are studied using non-periodic Density Functional Theory (DFT); for extended clusters we turn into the less expensive semiempirical tight-binding (TB) calculations; while periodic DFT and TB approaches are used to model infinite crystalline ZIFs. The calculations to study the stability of ZIF clusters with up to six Zn 2+ cations and a reduced number of explicit H 2 O or CH 3 OH solvent molecules, or NO - 3 counterions were performed with the Gaussian 09 (Revision A.02 code). [START_REF] Frisch | Gaussian-09 Revision A.02[END_REF] The computations used the long-range corrected hybrid density functional ωB97X, [START_REF] Chai | Systematic optimization of long-range corrected hybrid density functionals[END_REF] to account for van der Waals interactions and the large 6-311++G(d,p) triple-ζ basis set with polarization and diffuse functions, in order to reduce basis set superposition errors. Since energy minimization of zeolitic materials containing water or other strongly interacting solvent molecules requires very long computational times, 41 implicit solvation has also been considered (more details below) by means of the Self-Consistent Reaction Field, using the Polarizable Continuum Model (PCM), [START_REF] Tomasi | Quantum mechanical continuum solvation models[END_REF] which takes into account the solvent implicitly, -i.e. it does not include the explicit presence of solvent molecules, but it does include the electrostatic effects that they would have over the reactants. We have included the PCM parameters to mimic the effect of ethanol using the default parameters in Gaussian (ϵ = 24.852, and UFF atomic Van der Waals radius).

The stability of the clusters was evaluated using the free energies, calculated by including the zero-point (ZP) energies and the vibrational entropies, obtained from analytical calculations of the vibration modes at the energy minima configurations. The solvation free energies of the clusters were calculated using the thermodynamic cycle proposed by Dudev and Lim, [START_REF] Dudev | Tetrahedral vs octahedral zinc complexes with ligands of biological interest: a DFT/CDM study[END_REF] and da Silva et al., [START_REF] Eirik | Explicitly representing the solvation shell in continuum solvent calculations[END_REF] in which the same standard states are used for each species involved in the cluster formation reactions. For the calculation of the energy barriers and the finding of the transition states we have considered the NEB-TS method (from Nudged Elastic Band with TS optimization) as implemented in the ORCA program, [START_REF] Neese | The ORCA quantum chemistry program package[END_REF] (version 5.0.3). The same level of theory as in the previous Gaussian calculations was used, including the UFF parameters for the implicit solvent.

Tight binding calculations on large clusters, with up to 24 Zn 2+ cations, whose sizes make unfeasible the use of ωB97X

Table 1: Free energies of formation (in eV) for the complexes formed by a single Zn 2+ cation with up to four ethanol molecules (E ), water molecules (w ), nitrate anions (N) or imidazolate anions (Im -). As derived from computations in vacuo (v) or in implicit ethanol solvent (s). The thermodynamic cycle method, described in Sec. 2, is employed to calculate the free energies of formation. Values for complexes with net coordination numbers 1, 2, 3 and 4 are shaded in light yellow, blue, green and red color, respectively, to ease the interpretation of the energy trends. In italic, we show the label, with the format zlewn, where z, l, e, w , and n represent the number of Zn atoms and Im, ethanol, water, and nitrate molecules, respectively. These labels are employed in Figure 1. calculations, were performed with the extended semiempirical tight-binding method GFN2-xTB, [START_REF] Bannwarth | GFN2-xTB-an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions[END_REF] as implemented in the xTB package (version 6.4.0). [START_REF] Bannwarth | Extended tight-binding quantum chemistry methods[END_REF] In the GFN2-xTB calculations the implicit ethanol solvation was performed using the ALPB (from Analytical Linearized Poisson-Boltzmann) method. [START_REF] Ehlert | Robust and efficient implicit solvation model for fast semiempirical methods[END_REF] The periodic calculations with the GFN2-xTB were performed with the DFTB+ code (version 21.2). [START_REF] Hourahine | DFTB+, a software package for efficient approximate density functional theory based atomistic simulations[END_REF] Accurate DFT calculations of periodic systems are performed using the Meta-Generalized Gradient Approximation, non-local SCAN+rVV10 exchange-correlation functional, [START_REF] Peng | Versatile van der waals density functional based on a metageneralized gradient approximation[END_REF] as implemented in the Vienna ab initio Package (VASP, version 6.0). [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF] The plane wave kinetic energy cutoff was 500 eV. The projector augmented wave (PAW) method was used to describe the interactions between atomic cores and valence electrons. [START_REF] Blöchl | Projector augmented-wave method[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] Finally, the AIMD simulations were carried out with CP2K code, [START_REF] Kühne | CP2k: An electronic structure and molecular dynamics software package -quickstep: Efficient and accurate electronic structure calculations[END_REF] with a 1 fs timestep. The PBE exchange-correlation functional was employed, 55 using double-ζ basis sets with Goedecker-Teter-Hutter (GTH) pseudopotentials, 56 with 1, 4, 5, 6, and 12 valence electrons for H, C, N, O, and Zn atoms respectively. Dispersion was included via Grimme's D3 method. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-d) for the 94 elements hpu[END_REF] The NVT ensemble was employed, with the density of the system being that of liquid ethanol at the temperature of the AIMD run, namely 300 K.

Bare E 1 E 2 E 3 E 4 w 1 w 2 w 3 w 4 N 1 N 2 N 3 N 4 Zn 2+ (v) - - 5
The work carried out in this study involves a considerable amount of high-performance computing (HPC) resources, close to two million computational hours. Most of the computational effort is due to the large number of steps required to avoid imaginary frequencies in the optimization process for the high theory level ωB97X/6-311++G(d,p). Previous work showed that zeolitic structures having imaginary modes not only exhibit higher energies than those with all real modes, but also do not compare well with the results from high resolution experimental techniques. [START_REF] Ruiz-Salvador | Space group symmetry and Al-O-P bond angles in alpo 4 -5[END_REF] We have observed that clusters with Im -have a flatter energy surface than clusters with mIm -, resulting in slow convergence: in some cases the computational time for geometry relaxations and frequency calculations has been as long as months using 24 CPUs, -i.e. the Zn 4 Im 12 cluster or Zn 6 Im 18 , and more than 4 months with 40 CPUs were spent trying to optimize the Zn 8 Im 20 (D4R) cluster without success.

Results and Discussion

Effect of solvation and counterions

We devoted the initial stage of the calculations to an assessment of the effect of solvent and counterions on the structure and relative stability of the metal-organic clusters. On the one hand, we were interested in determining the differences between direct computations of the complexes in vacuo and computations with implicit solvent (ethanol) through the thermodynamic cycle mentioned in Sec. 2. On the other hand, in order to explore sensible aggregation and eventual nucleation routes, it is particularly relevant to determine to what extent the incorporation of explicit solvent molecules (ethanol, water) or counterions (nitrate, NO - 3 ) to the first coordination sphere of the Zn 2+ cation contributes to the stabilization of the complexes. Note that ethanol is a common solvent medium and the NO - 3 anion is a common counterion in the Zn salts used in ZIF synthesis. Water is included as zinc nitrate is usually included as the hydrated salt in synthesis mixtures.

The results obtained for the clusters formed by a single Zn 2+ cation serve to illustrate some of the main findings derived from these computations. Figure 1 depicts an illustrative ensemble of ZnIm n X m complexes (X = ethanol, water or nitrate, n = 1-4, m = 0-3) and Table 1 provides the corresponding free energies of formation in vacuo and under implicit solvation in ethanol, while highlighting the coordination number (2, 3 or 4) in each of the complexes.

Under implicit solvation (computed using the thermodynamic cycle), we found that the binding of the imidazolate (here denoted as Im, as well as Im -, when mentioning of the anionic nature of the ligand is relevant; for complexes we also 2) are indicated in parentheses next to each complex (in eV). Note that while the tetrahedral [ZnIm 4 ] 2-complex is stable, the incorporation of ethanol, water or nitrate is only stable up to a net three-fold coordination (see also Table 1). Atom color code: Zn (smoke), C (gray), N (blue), O (red), H (white). The labels of the clusters are explained in Table 1.

omit to mention the charge when it is not relevant) and nitrate (N) anions may proceed up to the tetrahedral complexes [ZnIm 4 ] 2-and [ZnN 4 ] 2-. Table 1 shows that the free energy of formation becomes increasingly negative for coordination numbers up to four in these complexes (five-and six-fold coordinations, with free energies of -4.11 eV and -3.39 eV respectively, are found to be unstable when compared with four-fold coordination (-4.85 eV), but of similar stabilities than two-and three-fold coordination (-3.74 eV and -4.41 eV respectively). This trend is, however, not reproduced by the calculations in the gas phase, in which the favored coordination is limited to 3 Im -ligands (the formation energies of [ZnIm 3 ] -solvated with ethanol or water is -27.9 eV, lower than that of the four-coordinated [ZnIm 4 ] 2-). Hence, implicit solvation turns out to be essential to account for the stable tetrahedral anionic coordination of the Zn 2+ cation. It seems timely to point out that the implicit solvent does not alter the geometry of the ZnIm n X m complexes appreciably with respect to that in vacuo: the Root-Mean-Square Deviations (RMSDs) between the superimposed atomic coordinates of the optimized clusters determined with both methods are negligible -e.g. < 0.01 Å in bond distances-. Changes in the energetics of these small clusters should therefore not be attributed to the inclusion of the implicit solvent, but rather to the unbalance of negative charge in the clusters and to the shielding of the long-range repulsion between the anionic ligands, which changes the energetics but does not affect the geometries, for small clusters. Differences in the structural features of the complexes predicted by the vacuo and solvated computations do become increasingly relevant with growing cluster size. For instance, the unshielded repulsions in vacuo between the Im -moieties restrict the "folding" (distortion of the structures to form more dense structures -see Figure S1 for an example) of the extended metal-organic networks with several Zn atoms, discussed below in Sec. 3.2 Differences between in vacuo and solvated computations also emerge for the coordination of Zn 2+ with explicit solvent molecules. Whereas both approaches agree qualitatively in the stability of the tetrahedral [Znw 4 ] 2+ water (w) complex, for ethanol (E), the three-fold complex [ZnE 3 ] 2+ is favored over the [ZnE 4 ] 2+ complex under implicit solvation, but it is not so pronounced in vacuo. Destabilization of the tetrahedral ethanol complex results from steric repulsion between the bulky ethyl groups, which is apparently compensated for by the unshielded, hence stronger, Zn 2+ . . . O interactions in the computations performed in the gas phase.

Imidazolate anions exhibit stronger interactions with the metal cation than nitrate or the solvent molecules and have a greater influence on the coordination number. Noticeably, it is found that, once Im -anions coordinate with the Zn 2+ , the incorporation of ethanol, water, or nitrate to the complex is only favored up to a total coordination number of 3. In other words, under implicit solvation, ZnIm n X m complexes are stable up to n + m = 3 (for n =1-3). For instance, the complexes ZnImX 2 and ZnIm 2 X are appreciably more stable than their ZnImX 3 , ZnIm 2 X 2 four-fold coordinated counterparts. Note that the nitrate anion exothermically displaces ethanol or water in these complexes (see Figure 1 and Table 1). These energetic trends also hold systematically for all the clusters with more than one metal cation explored in this 1 work. Under implicit solvation, the Zn 2+ sites in Zn z Im n X m clusters of any size create a four-fold tetrahedral coordination with Im -anions, but they bind explicit methanol and water molecules or nitrate anions only up to a three-fold coordination. This general behavior guides the initial stages of Zn-Im aggregation, as described in more detail in Sec. 3.2 below. In contrast, in vacuo, these trends only hold for nitrate, and not for ethanol or water, for which the four-fold coordinated complexes (n + m = 4) are predicted to be the most stable. We believe that this might be due to the fact that nitrate ions are negatively charged, and thus the Coulombic repulsions between ligands preclude the full coordination.

In order to test whether explicit solvation (including solvating molecules) is also required for the realistic modeling of the nucleation process, we sought insights into the reaction where an imidazolate reacts with a ZnIm 3 cluster, which is already interacting with an ethanol solvent molecule. The final state would be the tetrahedral ZnIm 4 cluster, in which the imidazolate has displaced the solvent molecule: [ZnIm 3 E] - + Im -→ [ZnIm 4 ] 2-+ E . First, we investigated the interactions of ZnIm 3 clusters in an AIMD simulation of ZnIm 3 mixed with 29 ethanol molecules in a cubic simulation cell of length 14.6 Å. A snapshot of the simulation is shown in Figure 2.

During the simulation we observe that the planar ZnIm 3 cluster binds to an ethanol molecule, but the interaction is labile, and the ethanol molecule can easily move. The structure of the ZnIm 3 cluster is kept in a nearly planar fashion, in agreement with the energy minimization study mentioned above, which suggests that only the presence of 4 imidazolates can induce the adoption of a tetrahedral coordination by the Zn atoms. These findings are in agreement with the geometries that would be predicted by the VSEPR theory, in which three charges (the ligands) would organize as a triangle in whose center lies the metal cation, while four charges would be localized in the vertices of a tetrahedron to minimize repulsion.

We further carried out NEB-TS calculations to assess the energy profile between the ZnIm 3 E and ZnIm 4 states. The results of these computations are depicted in Figure 3.

Energy / eV

Reaction coordinate, q / -1 Explicit Etha. and NO + Implicit Etha.

(a)

(b) (c) (d) (b) (c) (d) (a)
N] 2-+ Im -→ [ZnIm 4 ] 2-+ N -.
In the starting system (a), the imidazolate ion is negatively charged (the molecule is deprotonated), and the ethanol molecule is neutral, as these are the expected protonation states in the nucleating solution. Atoms color code as in Figure 1.

The main conclusion we can draw from this analysis is that the presence of the explicit solvating molecule is key to obtaining an accurate profile along the reaction path. For instance, the profile calculated with the explicit ethanol molecule shows a local minimum, corresponding to the system in which the H atom from the OH group of the ethanol molecule has been transferred to one of the imidazolates of the ZnIm 3 E cluster, while the profile of the ZnIm 3 cluster, with no explicit ethanol, does not show any local minima. The energy barrier for the latter process is very small at 0.16 eV, compared to the 0.98 eV for the former), -i.e. an MD simulation in which no ethanol molecules (nor artificially added energy barriers to mimic them) were present, [START_REF] Yang | Enabling computational design of ZIFs using ReaxFF[END_REF] would not be able to realistically describe the transition states involved. On the other hand, we see that the presence of an explicit solvent molecule is not needed for achieving a good description of the equilibrium energies, since the energies of the initial and final states are correctly described using only implicit solvation. We have also studied the barrier when the solvating molecule is a NO - 3 counter anion, instead of an ethanol molecule. Due to the repulsive interaction between the Im -and NO - 3 moieties, the energy of the system easily reaches a value close to that of the transition state, ∼ 1.3 eV above the energy of the initial state when an ethanol molecule is included explicitly, and ∼ 1.0 eV when no explicit ethanol is present. However, the equilibrium energies for reactants and products are quite similar.

From the above discussion, we find that implicit solvation is essential to accurately describe these clusters and their assembly and will now consider only such computations for describing the energetics.

Initial stages of cluster growth

The energetic and structural features outlined in Sec. 3.1 provide a computational framework to now consider a pool of Zn z Im n N m building blocks relevant to the nucleation of ZIF materials. Figure 4 depicts the most stable conformations of clusters with up to four Zn 2+ cations, which constitutes the largest cluster size explored in previous ab initio computational studies of works. [START_REF] Boussouf | Structure, spectroscopy, and bonding within the Zn q+ -Imidazolen (q = 0, 1, 2; n = 1-4) clusters and implications for Zeolitic Imidazolate Frameworks and Zn-enzymes[END_REF][START_REF] Terban | Early stage structural development of prototypical zeolitic imidazolate framework (ZIF) in solution[END_REF] The characterization of these clusters allows us to assess the potential routes leading to the simplest Zn-Im cyclic networks. We observe that nitrate anions incorporated into the clusters build two-or three-fold coordination units with Zn 2+ sites, thus stabilizing the imidazolate-poor clusters. We summarize now the fundamental trends derived from the cluster configurations compiled in Figure 4. From here onwards, energies are given in eV per Zn cation.

Zn1 clusters: In clusters with only one Zn 2+ cation, there is a clear stabilization as imidazolate coordination increases reaching a maximum for the four-fold coordinated [ZnIm 4 ] -2 complex. The solvent is in excess with respect to the nitrate, so initially every Zn 2+ cation will be solvated with ethanol molecules. But, if enough nitrate ions are available, these will gradually displace the coordinating ethanol molecules (see also energetics in Table 1). The following clustering route can be envisaged, considering stabilization due to nitrate (formation energies in parentheses): 4.85 eV). Note that, in agreement with our conclusions outlined in Sec. 3.1, the [ZnIm 3 N] -2 (-4.02 eV) complex is energetically unstable with respect to the [ZnIm 3 ] -complex. Furthermore, a subsequent increase in Zn-imidazolate coordination to achieve the (square pyramidal) five-and (octahedral) six-fold coordinated clusters [ZnIm 5 ] -3 and [ZnIm 6 ] -4 (-4.11 eV and -3.39 eV respectively) induces a destabilization of the cluster, supporting the argument for pre-nucleation being dominated by tetrahedral Zn species.

[ZnImN 2 ] - (-3.02 eV) → [ZnIm 2 N] -(-3.79 eV) → [ZnIm 3 ] -(-4.41 eV) → [ZnIm 4 ] -2 (-
Zn2 clusters: Imidazolate coordination with two Zn 2+ cations saturates in the doubly tetracoordinated [Zn 2 Im 7 ] -3 cluster. A minimum energy route is then predicted in which imidazolate sequentially displaces the nitrate counterions, as [Zn

2 ImN 4 ] -(-2.93 eV) → [Zn 2 Im 2 N 3 ] - (-3.32 eV) → [Zn 2 Im 3 N 2 ] -(-3.82 eV) → [. . . ] → [Zn 2 Im 5 ] -(-4.32 eV) → [. . . ] → [Zn 2 Im 7 ] -3 (-4.77 eV).
To adopt four-fold coordination of the Zn cations, large stress appears in the Zn2-cyclic clusters that causes instability with respect to the more stable planar three-fold coordinated structures (see complex [Zn 2 Im 5 ] -in Figure 4), in line with what is observed in ZIFs, which do not show Zn2-cyclic motifs.

It is pertinent to note that, under typical synthesis conditions, there is an excess of imidazolate anions in solution with respect to the metal cations. Hence, cluster growth by addition of Im -anions is kinetically more likely than by addition of cations. One further reason supporting that Im -binding to the growing clusters may be faster than Zn 2+ binding is related to the two-fold charge difference, which should lead to lower electrostatic energy barriers for imidazolate addition. Moreover, while Im -anions stabilize with two-fold coordination, Zn 2+ cations require a higher number of coordinating ligands to be stabilized. For these reasons, our calculations suggest that there will a faster growth by addition of imidazolate or nitrate anions (in the downward direction of the chart in Figure 1 than by addition of zinc cations (moving left to right in Figure 1).

Another interesting remark is about the potential coalescence of two ZnIm 4 clusters (the most stable Zn1 form) to yield Zn 2 Im 7 (the most stable Zn2 form). The calculated free energy of the reaction ZnIm 4 +ZnIm 4 → Zn 2 Im 7 +Im is actually positive (+0.17 eV), so there is no thermodynamic driving force for this reaction to occur. A plausible route would be the reaction ZnIm 4 + ZnIm 3 → Zn 2 Im 7 , with a favorable free-energy balance (-0.28 eV). This suggests an active role of the [ZnIm 3 ] -cluster in the pool of aggregation precursors of ZIFs, which would be in consonance with postulations from previous studies for analogous Co 2+ -imidazolate ZIF frameworks, [START_REF] Filez | Elucidation of the pre-nucleation phase directing metalorganic framework formation[END_REF] and from recent studies, [START_REF] Salvador | Computer simulation of the early stages of self-assembly and thermal decomposition of ZIF-8[END_REF][START_REF] Macià Escatllar | Modeling hydroxylated nanosilica: Testing the performance of ReaxFF and FFSiOH force fields[END_REF] which found that tri-coordinated species are still present over long simulation times. Another possible reaction would be ZnIm 4 +ZnIm 2 → Zn 2 Im 6 , followed by the exchange of coordinating ethanol or nitrate molecules by an imidazolate molecule to form Zn 2 Im 7 (free energy balance of -0.95 eV). This reaction has a higher free energy, probably due to the lower stability of the ZnIm 2 cluster with respect to the ZnIm 3 cluster. Zn3 clusters: While Zn3 ring structures are the most stable in the pool of precursor Zn 3 Im n X m clusters (n = 3-8, m = 0-3), the full four-fold coordination of three Zn 2+ cations with imidazolate anions is marginally more favorable in a cyclic [Zn 3 Im 9 ] -3 cluster (-4.74 eV) over the linear [Zn 3 Im 10 ] -4 cluster (-4.71 eV) ones. This again points to the effect of deformation of the local structure, since the stress induced by the constraint of the formation of cyclic structures reduces the stabilization gained by the higher coordination achieved when cycles are formed. The opening of the ring could be induced in the last stage upon the incorporation of an additional imidazolate anion to the Zn 3 Im 9 cluster, to form the linear Zn 3 Im 10 cluster, but this does not seem to be the case (∆g = + 0.08 eV). Rather, this latter cluster may be mainly produced through the aggregation reaction [Zn 2 From the structures reviewed so far, it can be inferred that, despite the tendency of Zn cations to eventually bind to four imidazolate molecules forming a tetrahedron, in a similar fashion as in ZIFs crystal structures, most of the reactions that drive cluster growth involve the presence of three-fold coordinated moieties. For instance, the low free energies of the [ZnIm 3 ] -and [Zn 2 Im 5 ] -clusters (-4.41 eV and -4.32 eV respectively) suggest that a large number of these threecoordinated, planar species are present in solution during the initial stages of aggregation.

Zn4 clusters: Figure 4 shows that, of the clusters incorporating four Zn 2+ cations, the lowest energy structures are dominated by cyclic structures, culminating in the fully coordinated [Zn 4 Im 12 ] -4 cycle, which constitutes the most stable structure of this class (-4.74 eV). Remarkably, the analogous tetrahedrally coordinated linear arrangement [Zn 4 Im 13 ] -5 , shown in Figure S1, right, is significantly less stable (-4.27 eV). Already cyclic clusters containing four Zn cations exhibit Zn-Im-Zn unstressed angles that allow the intermolecular forces stabilizing these more compact oligomers that are topologically similar to the extended ZIF structures. It can be noted that the NO - 3 counterions transiently stabilize the Zn 2+ sites, while imidazolate units progressively incorporate to the cluster, leading to three-fold, then four-fold coordination arrangements. An interesting situation arises in the [Zn 4 Im 8 ] 0 complex, where a cycle with two tetrahedral (ZnIm 4 ) sites, stabilized by two nitrate anions (-4.21 eV), is very similar in energy to the planar configuration in which all four Zn cations are three-fold coordinated to Im -(-4.20 eV). However, the difference is not significant enough to draw any particular conclusion here. Also noticeable is the observation of a [Zn 4 Im 9 ] -configuration in which an imidazolate cation bridges two Zn sites across the Zn4 ring.

The results described above provide a consistent rationalization of reported experimental results. A study carried by Lim et al. [START_REF] Ivy | Insights into the molecular assembly of zeolitic imidazolate frameworks by ESI-MS[END_REF] found that monomeric species are the most common in EIMS spectra recorded from the analysis of an in situ gel synthesis of ZIF-8 at the very early stages of the nucleation, followed by the appearance of other species, containing up to four Zn atoms, whose concentrations decrease over time after 20 min of reaction. In our calculations, the most stable monomeric species, Zn 1 L 4 , have similar stability to the more stable dimeric and trimeric ones (Zn 2 L 7 and Zn 3 L 10 ), which suggests that in solution the three species would have similar populations. From a kinetics perspective, it can be thus expected that crystal growth will more likely take place from monomeric species than from dimeric species, since it will be kinetically more favorable to adapt the structure of a monomer to be inserted into a growing surface than to rearrange a larger, more rigid dimeric structure. This view is in agreement with the conclusions Moh et al. [START_REF] Pak | Revelation of the molecular assembly of the nanoporous metal organic framework ZIF-8[END_REF] defined from the AFM analysis of steps heights on ZIF-8 surfaces, which suggest that crystal growth proceeds through the incorporation of monomeric species at the surface, as well as with observations by Balestra and Semino. [START_REF] Salvador | Computer simulation of the early stages of self-assembly and thermal decomposition of ZIF-8[END_REF]

Structure of ring clusters

In order to span the conformational landscape involved in the ZIF nucleation process, we have computed the formation free energies of a broad ensemble of large Zn clusters with either imidazolate or 2-methylimidazolate. Specifically, 42 Zn z Im l clusters and 20 Zn z mIm l clusters, with z ranging from 1 to 6 and l ranging from 1 to 18 have been characterized. The corresponding energies obtained from thermodynamic cycle computations are listed in Table S3 and Table S4, respectively. These calculations are expected to shed some light on the observation that mIm -leads readily to (only) the formation of ZIF-8, while Im -leads to other topologies.

The cyclic Zn4 clusters tetrahedrally coordinated to Im -or mIm -already resemble seminal core structural features present in most ZIFs. The configurational space of the Zn4 clusters becomes increasingly complex, with various possible configurations for the 4-membered rings (4MRs). We studied two types of 4MRs, namely a 4MR extracted from the SOD crystal structure, [START_REF] Kyo | Exceptional chemical and thermal stability of zeolitic imidazolate frameworks[END_REF] of ZIF-8 (labeled 4MR-SOD), and a 4MR extracted from the zni crystal structure, [START_REF] Lehnert | Darstellung und kristallstruktur des mangan(II)-und zink(II)-derivates des imidazols[END_REF] of ZIF-zni (labeled 4MR-zni), as shown in Figure 5a. It is worth noting that ZIF-zni has very low porosity, and is the most stable of the Im-bearing ZIFs, as observed both experimentally, [START_REF] Hughes | Thermochemistry of zeolitic imidazolate frameworks of varying porosity[END_REF][START_REF] Schröder | Metastable metal imidazolates: development of targeted syntheses by combining experimental and theoretical investigations of the formation mechanisms[END_REF] and in periodic DFT calculations. [START_REF] Baburin | Enumeration of not-yetsynthesized zeolitic zinc imidazolate MOF networks: A topological and DFT approach[END_REF][START_REF] Lewis | Zeolitic imidazole frameworks: structural and energetics trends compared with their zeolite analogues[END_REF] The synthesis of ZIFs with Im -ligands often yields low porosity ZIFs, mainly ZIFzni, and in order to obtain a ZIF with the sodalite topology a post-synthetic ligand exchange process must be undertaken. [START_REF] Karagiaridi | Opening ZIF-8: A catalytically active zeolitic imidazolate framework of sodalite topology with unsubstituted linkers[END_REF] Indeed, our results show that, for Im-bearing clusters, the 4MR-zni is the most stable of the ring clusters. It is slightly more stable than 4MR-SOD (formation energies of -4.74 eV and -4.70 eV respectively, see in Table S3), and, as expected, much more stable than the undercoordinated 4MR (with formation energies of less than 4 eV).

Our calculations suggest there are two factors suggestive of why the zni topology is more stable than the SOD topology for Im-based ZIFs. Firstly, the 4MR-zni cluster is slightly more stable than the 4MR-SOD cluster (-4.74 eV vs. -4.70 eV). But more importantly, the dihedral angle between the four Zn atoms in the 4MR-zni Im cluster (14.5 • , see Figure S2) is similar to that of the 4MR in the zni crystal structure (∼ 11 • ), in contrast to that of the mIm-based cluster 4MR-zni (3.3 • ). Thus, the formation of an Im-based SOD topology would lead to a strained, and thus less stable, crystal structure.

Interestingly, for clusters with mIm -ligands, the situation is reversed, with the 4MR-SOD cluster being more stable than 4MR-zni (see Table S4). This reversal in stability might be related to the repulsions between the methyl groups of the mIm -ligands, which in mIm-4MR-zni clusters are pointing in opposite directions, thus increasing the repulsions between the groups (see Figure 5c), with respect to the mIm-4MR-SOD cluster (Figure 5d). Note that in the Im-4MR-zni cluster, one imidazole group is rotated, to increase the distance between the H atoms (Figure 5a). This increase in H-H distance to reduce repulsions induces an increase in Zn-Zn distances, in both the mIm-based clusters and crystal structures, since the distance between first neighbor Zn atoms is 6.1 Å in both 4MR and 6MR mIm -clusters, which is similar to that of the SOD crystal structure, but larger than that of the Im-based zni crystal structure (5.85 Å). The dihedral angles in the 6MR clusters are ∼ 0 • (as in the SOD structure) for both Im -and mIm -clusters, explaining why both ligands can form ZIFs with 6MRs in their topologies.

Ligand rotations can occur freely in the clusters we are studying, but as the size of the clusters increases, approaching bulk-like systems, they would be hindered by long-range interactions. This observation provides another reason why the SOD structure of ZIF-8 is the preferred crystal formed with mIm -ligands, since the rotated ligands can point their methyl groups into the larger cavities that, in periodic structures, would be placed at the left and right of the clusters shown, thus avoiding close contacts. The influence of ligand rotation on the stabilities of the different topologies is a complex subject. A controlled rotation of ligands permits the synthesis of new topologies, [START_REF] Yang | Principles of designing extra-large pore openings and cages in zeolitic imidazolate frameworks[END_REF] but more research is needed to obtain a detailed understanding of the relation between ligand rotation and crystal structure stability, something that is out of the scope of this study.

Towards clusters with bulk-like structures: Sizes up to 3 nm

Prior to a study of the (relative) stability of larger clusters, we validated the calculation of the free energies using TB methods, based on the DFT values already calculated as benchmarks. Figure 6 represents the formation energies of all the clusters studied in the previous sections (with Im -and mIm -ligands) calculated with the DFT ωB97X method and the semiempirical GFN2-xTB method. We find a good correlation between the two levels of theory, which will allow us to perform GFN2-xTB/ALPB(ethanol) calculations to study larger clusters, with an accuracy that is very similar to that of the much more time consuming ωB97X/6- Figure 7: Geometries and energies (in eV, in parenthesis) of some relevant large clusters, with a number of Zn 2+ cations ranging from 5 to 144. The number of Zn 2+ cations and Im -ligands are shown using the (z, l) notation. The energy is calculated from the linear regression reported in Figure 6. Color code as in Figure 1.

- , calculated with three types of calculations, namely the periodic meta-GGA exchange-correlation functional SCAN+rVV10, with the GFN2-xTB method, and with the non-periodic, hybrid exchange-correlation functional ωB97X/6-311++G(d,p). Each line represents a fitting of the cluster energies corresponding to a particular value of z (with different number of ligands, l). The black circles correspond to the energies of the clusters with z > 8. Red empty circles correspond to the energies of the crystalline materials calculated with the GFN2-xTB method (zni, cri, coi, ZIF-1, ZIF-4, ZIF-6, ZIF-8, SALEM-2, and ZIF-10). Empty triangles correspond to fitted ∆g of the corresponding cluster calculated with the periodic meta-GGA exchange-correlation functional SCAN+rVV10. These values are obtained following the methods described in the ESI. The inset shows the region of larger clusters and of periodic systems.

311++G(d,p)/PCM(ethanol) calculations.

In Figure 7 we show the geometries of a selection of the largest clusters included in this study, with a number of cations ranging from 5 to 144. The energies of all the clusters are reported in Table S5. We find that fully coordinated Zn clusters have energies ranging between -4.00 and -5.54 eV (per Zn cation) while clusters with undercoordinated Zn atoms are less stable. There is an unexpected high energy for the cluster sod-Zn 24 (with all 4-fold Zn cations and forming 6MR), which can be ascribed to the large cavity inside the cluster, making it less stable, as a large number of solvent molecules must be excluded to make space for it, so that this high value might be an artifact of the method employed to calculate the energy. The formation of bulk-like clusters can occur after the clusters have achieved a critical size, while clusters smaller than that size feature non bulk-like geometries. The transition from small clusters to bulk-like clusters via intermediate polymorphs (structures which in bulk phase are metastable) is commonly observed in Zn-based materials, such as ZnS and ZnO, and other inorganic materials (TiO 2 , etc.), 67-69 and it is an example of Ostwald's Rule, which states that the first crystal structure formed would be the least stable polymorph, and it also applies to MOFs. [START_REF] Cheetham | Thermodynamic and kinetic effects in the crystallization of metal-organic frameworks[END_REF][START_REF] Van | In situ, time-resolved, and mechanistic studies of metal-organic framework nucleation and growth[END_REF][START_REF] Yeung | Control of metal-organic framework crystallization by metastable intermediate pre-equilibrium species[END_REF] It is only after a critical nanoparticle size is achieved that the stable crystalline phase appears. Indeed, for the larger cluster sod-Zn 144 , the formation energy is much lower (-5.49 eV) than that of the sod-Zn 2 4 cluster (-3.83 eV), and that of the periodic crystal structure SALEM-2 (sodalite) is -5.86 eV, while the most stable crystal structure is zni (-6.18 eV, see below for a discussion). It is interesting that the relatively stable large Zn 144 cluster (with both 4 and 6MRs) contains six relatively unstable Zn 24 clusters. But the assembly of the Zn 24 clusters results in further 4MRs being generated, which appears to stabilize the larger cluster. The role of the assembly of small units in forming large pore zeolites was postulated three decades ago, [START_REF] Bnmner | Framework density distribution of zeolite-type tetrahedral nets[END_REF] but not seen before in MOFs. Another salient point from Figure 7 is the large stability of the (41, 88) cluster (formation energy -5.54 eV). This structure was generated by MD simulations at high temperature; the resulting structure is a continuous random network containing most tetracoordinated and some undercoordinated Zn atoms. The core of this cluster is dense, having no internal porosity, which greatly increases stability, overcoming the penalty associated with the (small) number of undercoordinated Zn cations present in the structure. This is the first cluster we find to be relatively more stable (per Zn cation) than ZnIm 4 . It is important to note that this cluster contains a z = 41 value considerably lower than the others with similar energies (60 and 144). This is in agreement with the experimental observation that glassy clusters form before crystalline clusters become more stable. Likewise, it also explains the formation of glassy ZIFs under specific reaction conditions, 71 as was found by Balestra and Semino 25 by computer simulations.

In order to get a better understanding of the pre-nucleation phase of ZIFs growth we aimed to investigate how clusters increase their stabilities as they grow. Experimentally, the most stable system that it is possible to create with Zn 2+ cations is the infinite, periodic crystal structure of ZIF-zni, but the ωB97X/6-311++G(d,p) calculations we have carried out so far (based on cluster models) cannot be employed to study periodic systems, so we cannot see how the energies of the clusters approach the value of the ZIF-zni crystal structure. To overcome this problem, we carried out two types of periodic calculations: the (very time-consuming) nonlocal SCAN+rVV10 exchange-correlation functional (as implemented in VASP), and the (less costly) GFN2-xTB/ALPB method. Both of these calculations show (see Figure 8) that, in agreement with experiments, the most stable crystal structure is that of ZIF-zni, with a formation energy of -6.18 eV (calculated with the GFN2-xTB/ALPB method). This correct prediction of crystal structure stabilities provides further support for the validity of our calculations in clusters presented above.

Given that the plane wave methods used to calculate the energy of the ZIF-zni crystal structure are very different to the molecular methods used to calculate the GFN2-xTB/ωB97X energies of the clusters, it could be reasonable to think that these two types of energies cannot be compared. For that reason, we also plotted the energies of the neutral clusters, now obtained with the same meta-GGA plane wave calculations. We see that the energies of the clusters calculated with both types of calculations are very similar, so that we can indeed use the bulk energy as a reference energy towards which the cluster energies must converge, with some confidence. The results, shown in Figure 8, indicate that the clusters we are studying are still very far away from achieving the stability of the most stable bulk structure (that of ZIF-zni, with a free energy of -6.18 eV). In Table S6 we have listed the energies of the other calculated periodic structures. For that reason, our results suggest that clusters with bulk-like geometries are not the most stable at the cluster size considered here, so that further growth must take place before surface effects have a lower impact on the cluster structures and allow classical crystal growth to proceed. Interestingly, in Figure 8, we observe that the lines calculated by fitting the series of energies obtained for each value of z and for varying values of l, converge to a formation energy of ∼ -7.5 eV, for both l = Im and l = mIm, although there is no clear physical interpretation of this fact.

Conclusions

We can draw the following conclusions from this computational study: (a) The inclusion of implicit solvation (via dielectric embedding methods) is required to achieve chemically sound structures as opposed to in vacuo calculations. Explicit solvation is only needed if accurate energy barriers are to be calculated, but the energies of equilibrium structures do not seem to depend much on the inclusion of explicit solvent or counterions, (b) The most stable cluster (free energy of formation per Zn atom of -4.85 eV) in the initial process of cluster growth is ZnIm 4 , although clusters such as ZnIm 3 , Zn 2 Im 7 , Zn 2 Im 7 , Zn 3 Im 9 , Zn 3 Im 10 , or Zn 4 Im 12 do have energies close to -4.50 eV, and would therefore be present in solution at appreciable quantities. All these species, except ZnIm 3 , have a tetrahedrally coordinated Zn 2+ cation, (c) Zn x Im y clusters take longer than Zn x mIm y to be optimized, due to the shallower nature of their potential energy surfaces, which is related to the higher constraints induced by the larger sizes of the mIm groups, (d) The low porosity zni structure, which is the most stable crystal formed with Im -ligands, is less stable than the SOD structure of ZIF-8 when the Zn 2+ cations are connected with mIm -ligands. This might be due the presence of larger pores in SOD, which allow for rotation of the mIm -ligands to point their methyl groups into different places, thus avoiding close contacts. (e) The ZnIm 4 cluster is the most stable cluster determined for a wide range of sizes, and clusters do not start becoming more stable than this cluster until they have more than 40 Zn atoms, which suggests that the smaller clusters would not have long lifespans, since they are likely to undergo formation and dissolution processes to yield stable ZnIm 4 clusters, (f) Ostwald's rule of stages seems to be observed, as glassy clusters are formed before crystalline clusters become more stable. The crystal structures formed in the nucleation process are then determined by the kinetic factors that the synthesis conditions create, permitting the transformation of the glassy particles to crystalline particles.

To conclude, our survey of oligomeric and larger clusters that could be present during ZIF formation, is consistent with the available experimental data and provides additional insight into the key processes that lead to the formation of crystal nucleation species. We have shown how the variety of DFT, and tight binding methods are comparable in accuracy, and allow the consideration of a wide range of species. This study provides important insight into the pre-nucleation of ZIFs, and paves the way to achieving synthesis control by a fundamental understanding of its underlying physicalchemistry.

Figure 1 :

 1 Figure 1: Representation of the complexes formed by a single Zn 2+ cation with Im -, with added ethanol or water solvent molecules or nitrate counterions. The free energies of formation obtained under implicit solvation (thermodynamic cycle computations, see Sec.2) are indicated in parentheses next to each complex (in eV). Note that while the tetrahedral [ZnIm 4 ] 2-complex is stable, the incorporation of ethanol, water or nitrate is only stable up to a net three-fold coordination (see also Table1). Atom color code: Zn (smoke), C (gray), N (blue), O (red), H (white). The labels of the clusters are explained in Table1.

Figure 2 :

 2 Figure 2: Snapshot of an AIMD simulation of a ZnIm 3 cluster in liquid ethanol. Color code as in Figure 1

Figure 3 :

 3 Figure 3: NEB-TS calculations: In the bottom panel we show the energy vs. the reaction coordinate, q, for four systems whose geometries are shown in the top panel. The pink solid curve shows the NEB calculation between (a) ZnIm 3 E + Im and (d) ZnIm 4 + E (E for ethanol). The (c) cluster is the transition state, and (b) the cluster is an intermediate, local minimum state, which exhibits a proton transfer between the imidazolate and the explicit solvent molecule. The red curve corresponds to the calculation without the explicit ethanol molecule. The cyan curve corresponds to the calculation using as reactants [ZnIm 3 N] 2-(N for Nitrate) + E + Im -→ [ZnIm 4 ] 2-+ E + N -. The blue line corresponds to the calculation for [ZnIm 3N] 2-+ Im -→ [ZnIm 4 ] 2-+ N -.In the starting system (a), the imidazolate ion is negatively charged (the molecule is deprotonated), and the ethanol molecule is neutral, as these are the expected protonation states in the nucleating solution. Atoms color code as in Figure1.

Figure 4 :

 4 Figure 4: Optimized structures of Znz ImnNm, clusters, with z = 1-4. The values (z, n, m) are shown in parenthesis to label the clusters. The formation energies (in eV per Zn atom) are shown in each configuration, with values in red color indicating the most stable complexes within each (z, n) class. Atom color code as in Figure 1.

  Im 7 ] 3-+ [ZnIm 3 ] -→ [Zn 3 Im 10 ] 4-(∆g = -0.2 eV), again suggesting the potential role of [ZnIm 3 ] -and [ZnIm 2 ] -as building blocks in ZIF growth.

Figure 5 :

 5 Figure 5: Structures of various 4MR (with 4 cations and 12 ligands) clusters and 6MR clusters (with 6 cations and 18 ligands), either with Im -or mIm -, with total charges equal to -4 and -6, respectively. a) Zn 4 Im 12 with zni topology, b) Zn 4 mIm 12 with SOD topology, c) Zn 4 Im 12 with zni topology, d) Zn 4 mIm 12 SOD topology, e) Zn 6 Im 18 f) Zn 6 mIm 18 .

Figure 6 :

 6 Figure 6: Linear Regression using as target variable the formation free enthalpies per Zn 2+ cation, ∆g(ωB97X/6-311++G(d,p)), and as features: z, l, and ∆g(GFN2-xTB).

nat ( 6

 6

Figure 8 :

 8 Figure8: Values of formation free energy per Zn atom, ∆g, of Znz L l clusters (l =Im -and mIm -) as a function of the number of Zn atoms and the number of ligands, (z + l /2) -1 , calculated with three types of calculations, namely the periodic meta-GGA exchange-correlation functional SCAN+rVV10, with the GFN2-xTB method, and with the non-periodic, hybrid exchange-correlation functional ωB97X/6-311++G(d,p). Each line represents a fitting of the cluster energies corresponding to a particular value of z (with different number of ligands, l). The black circles correspond to the energies of the clusters with z > 8. Red empty circles correspond to the energies of the crystalline materials calculated with the GFN2-xTB method (zni, cri, coi, ZIF-1, ZIF-4, ZIF-6, ZIF-8, SALEM-2, and ZIF-10). Empty triangles correspond to fitted ∆g of the corresponding cluster calculated with the periodic meta-GGA exchange-correlation functional SCAN+rVV10. These values are obtained following the methods described in the ESI. The inset shows the region of larger clusters and of periodic systems.

-----------(g) -4.85 --------

-1 Explicit NO + Implicit Etha.

Acknowledgments

S. H. acknowledges funding from the Agencia Estatal de Investigación and the Ministerio de Ciencia, Innovación y Universidades, of Spain (PID2019-110430 G B-C22), and from the EU FEDER Framework 2014-2020 and Consejería de Conocimiento, Investigación y Universidad of the Andalusian Government (FEDER-UPO-1265695). S. R. G. B. was supported by grants FJC2018-035697-I funded by MCIN/AEI/10.13039/501100011033 (Ministerio de Ciencia e Innovación; Agencia Estatal de Investigación) and POSTDOC_21_00069 funded by Consejería de Universidades, Investigación e Innovación, Junta de Andalucía. We are thankful to C3UPO for the HPC facilities provided.

Data Availability Statement

The data that support the findings of this study are available within the article, and its available online Supplementary Information (Figures S1-S2, and Tables S1-S6) in https://github.com/salrodgom/ZIFNucleation/raw/ main/ESI/ESI.pdf.

Conflicts of interest

The authors have no conflicts to disclose.