Development and Application of a Single Neural Network Potential for IRMOF-n (n=1,4,6,7,10)
Metal-organic frameworks (MOFs) with their exceptional porous and organized structures have been subject of numerous applications. Predicting macroscopic properties from atomistic simulations require the most accurate force fields, which is still a major problem due to MOFs' hybrid structures governed by covalent, ionic and dispersion forces. Application of ab-initio molecular dynamics to such large periodic systems are thus beyond the current computational power. Therefore, alternative strategies must be developed to reduce computational cost without losing reliability. In this work, we describe the construction of a neural network potential (NNP) for IRMOF-n series (n=1,4,7,10) trained by PBE-D4/def2-TZVP reference data of MOF fragments. We validated the resulting NNP on both fragments and bulk MOF structures by prediction of properties such as equilibrium lattice constants, phonon density of states and linker orientation. The energy and force RMSE values for the fragments are only 0.0017 eV/atom and
INTRODUCTION<!>MATERIAL AND METHODS<!>RESULTS AND DISCUSSIONS<!>Neural<!>On the fly data sampling during the MD simulation<!>Validation and Test on Fragments<!>Performance on Bulk Properties<!>Lattice constants<!>Energy profiles on the organic linker's rotation in the bulk IRMOF-1<!>Phonon modes<!>Thermal behavior<!>Mechanical behavior
<p>Metal-organic frameworks (MOFs) with their highly ordered microscopic structure have attracted scientific interest due to their large surface areas and pore volumes, high flexibilities, and capabilities of reversible processes. For the past decade, they have been the subject of numerous applications ranging from energy storage materials to industrial gas adsorption [1][2][3][4][5], drug delivery [6][7][8][9][10], separation [11][12][13][14], and catalysis [15][16][17]. Quick and easy modification of metal sites and organic linkers allows rational design and versatile functionalization.</p><p>Most of the applications require well defined process conditions such as pressure and temperature. For instance, thermal expansion due to higher temperatures might lead to deformation of plasticity in materials [18][19][20]. In particular, adsorption and desorption processes of guest molecules without heat control system might cause drastic temperature changes in MOF structures.</p><p>On the other hand, adjusting temperature can be beneficial to controlled adsorption and desorption guest molecules for some applications such as temperature swing adsorption (TSA) [21,22]. Thus, thermal properties of MOFs are still interest of investigation. The most interesting thermal property is that most frameworks show negative thermal expansion (NTE), the decrease in the total volume upon heating. This property becomes critical and can be utilized when designing on-board energy related (CH4 and H2) storage materials [18]. It also creates thermal pressure that plays important role in the mechanical stability of the frameworks.</p><p>Adsorption and desorption properties have been extensively studied by both experimental and theoretical methods [23][24][25][26]. However, there are few theoretical studies on mechanical behaviors of MOFs such as applied temperature and applied pressure [27]. Computational studies on adsorption/desorption properties as well as thermal behaviors of MOFs rely on classical force fields. In order to understand the atomistic details of host-guest interaction, the reactivity and flexibility of MOFs, more accurate methods such as Density Functional Theory (DFT) are required. Therefore, alternative strategies must be introduced to caveat the problem.</p><p>The bottleneck with DFT methods is that MOF structures may include thousands of atoms in unit cells and first principles definition of such large systems are quite costly in terms of computation [28]. In addition, the existence of multiple transition metals with open shell structures (in inorganic subunits) further complicates the calculations. Moreover, dynamic definition of physical and chemical events requires extensive sampling of the system with MOFs and/or guestmolecules. Application of ab-initio molecular dynamics to such systems are thus beyond the current computational power. Therefore, alternative strategies must be developed to reduce computational cost without losing reliability.</p><p>Over the last decade, Machine Learning (ML) has been greatly applied in many different systems [28][29][30][31][32][33][34]. Highly accurate and reactive ML potentials can be built to describe MOF systems with both organic and inorganic subunits for simulations. Rather than building electronic structure from the first principles at every time step, their great capabilities of learning of multidimensional potential energy surfaces (PES) allows one to compute much faster simulations [32,[35][36][37][38][39][40].</p><p>Behler and coworkers were first to construct neural network potentials (NNPs) for MOF-5 (IRMOF-1) and successfully reproduced the electronic structure calculation [28]. They used the idea of fragmentation subsets of bulk MOF to train PES at the DFT level. The DFT reference energies belonging to non-equilibrium geometries of MOF-5 fragments have been used to train NNPs. Although they only used the non-equilibrium fragments rather than using bulk structures in the construction of NNPs, they could accurately predict bulk properties such as lattice parameters, elastic constants and negative thermal expansion. This study showed the applicability of NNPs on MOF like organic-inorganic hybrid systems. Furthermore, the NNPs were constructed by affordable DFT calculations of fragments rather than costly bulk structures [28].</p><p>In the construction of potentials, there are several artificial neural network (ANN) models that are mostly based on atomic contribution to the target property (e.g. potential energy surface) such as ANI-torch [28], HDNN [43], deep tensor neural network (DTNN) [28] and SchNet [32].</p><p>In this work, we describe the construction and validation of a generic neural network potential for IRMOF-n series (n=1,4,7,10) based on DFT reference data and predict bulk properties such as equilibrium lattice constants, phonon density of states and linker orientation. Our developed NNP enables to reveal unusual behaviors of MOFs such as the thermal expansion properties and the effect of mechanical strain on the adsorption of hydrogen and methane molecules for the selected MOFs.</p><!><p>Our NNP has been constructed via fragmentation of MOFs with an iterative process as illustrated in Figure 1. In the construction process, we used SchNet [32], which is an end-to-end machine learning model (open-source code is available in SchNetPack toolbox [43]). SchNet is a continuous-filter convolutional neural network model designed to predict molecular properties such as energies and atomic forces etc. modified version of deep tensor neural network (DTNN) framework [44]. It is a multilayer neural network consisting of atom embedding, interaction blocks and prediction blocks. At the embedding layer, each atom type is mapped to feature vectors 𝑥 𝑖 𝑙 = ℝ 𝐹 to create a randomly initialized representation. Here, 𝐹 denotes the feature dimension and 𝑙 is the number of layers for 𝑖'th atom. This representation is refined during training by the interaction block. The building block computes pairwise interactions between all 𝑥 𝑖 𝑙 within a certain cutoff and then refines atomic representations by a subsequent t steps up to T times in total. These atom-wise representations are processed to obtain the desired molecular property by the summation of all atomic contributions in the prediction block. In addition to calculating energies by reliable neural network (NN) models, SchNet also provides the force models obtained by differentiation of energy with respect to atomic positions. Rotationally and translationally invariant description of energy model allows conservation of both force field and energy. Details of the algorithm can be found in Schutt's works [32,43,44] Herein, we introduce a generic SchNet model for IRMOF-n (n=1, 4, 6, 7, 10) series using DFT energies and forces based on fragments of the MOFs. According to hyper-parametrization on metal organic hybrid molecular systems, 96 number of feature maps for atom embedding, a radial cutoff of 6.0 Å and T=3 interaction refinements have been used. The model was trained with a learning rate of 10 -4 using ADAM optimizer. Among the dataset, 80% was used for training while 10% was used for the validation and the remaining 10% was reserved for testing. All training processes have been performed using a Tesla P100 GPU.</p><p>The DFT energies and forces for all non-equilibrium geometries were calculated using the ORCA software package (version 4.2) [45] at the PBE/def2-TZVP level for all atoms. SCF calculations were accelerated with the RIJDX option available in the software to notably increase computational efficiency [46]. Tight SCF criteria were applied for all calculations and long range interactions were corrected by the use of density dependent atom-pairwise dispersion (D4) functions [47].</p><p>A custom-built python code was prepared using the Atomic Simulation Environment (ASE) [48] framework to perform calculation and store all geometries and their DFT results in an ASE database format. UCSF Chimera [49] and VESTA [50] software packages were used for visualization of 3D structures.</p><p>The simulations of fragments used to train NNP were carried out using SchNetPack software with the standard velocity Verlet integrator [43]. The time constant for the simulations was set to 0.5 fs, and Langevin thermostat was used. The final simulations of periodic MOF structures were performed to investigate bulk properties for phonon modes, thermal mismatch and mechanical stress using SchNetPack and ASE interchangeably. Phonons were produced by calculation of vibrational modes using finite differences as implemented in ASE. For thermal studies, MD simulations at different temperatures were performed ranging from 200 K to 400 K by 50 K increments. For each temperature value, the energy of the periodic bulk system was minimized.</p><p>After minimization, MD simulations of 500 ps with 0.5 fs time steps (total of 1 million steps) were carried out using NPT ensemble at 1 atm reference pressure. Langevin thermostat along with Nose -Hoover barostats were used. Damping constants for both thermostats were set to 0.1 ps. Last 300 ps of the simulations was used for the analysis.</p><p>Thermal expansion coefficient (𝛼 0 ), bulk modulus (𝛽 0 ) and thermal pressure coefficients (𝛾) are the main parameters that can be calculated for the thermal behavior [18].</p><p>Thermal expansion coefficient (𝛼 0 ) is defined as:</p><p>This can be re-written by integrating for volume and temperature;</p><p>where 𝑉 0 is the equilibrium volume at 0 K. By running several MD simulations with different temperatures, 𝛼 0 can be calculated from the linear correlation between temperature and natural logarithm of volume expansion. The term 〈ln (𝑉)〉 𝑃 can be obtained from the ensemble average of an MD simulation at the NPT ensemble.</p><p>Thermal pressure is given by:</p><p>where 𝛾 is the thermal pressure coefficient with:</p><p>where 𝛽 0 is the bulk modulus and defined as:</p><p>where 𝑘 𝐵 is the Boltzmann coefficient. The volumetric term in the equation corresponds to the relative variance, the fluctuation of the volume throughout the MD simulation at T=300 K.</p><p>For the mechanical behavior of MOFs, the uniaxial compression systems were introduced by assigning the [001] plane to a vector and its gradual decrease. In line with literature studies, the −3x10 9 s −1 strain rate was engineered as follows: During the MD simulations at the NVT ensemble, the box volume was reduced by decreasing the [001] plane by 0.03% after each ps simulation. A total 45% volume decrease in 150 ps total simulation time was achieved.</p><p>Adsorption properties of methane and hydrogen gases were studied using Grand Canonical Monte Carlo (GCMC) simulations with RASPA software code [51] using deformed/undeformed structures of IRMOF-1 and IRMOF-7 after extracting MD trajectories.</p><p>For both H2 and CH4 uptake, GCMC simulations for all deformed/undeformed structures were run for 1000 cycles with the first 500 being equilibration and the last 500 for the ensemble average. The methane adsorption calculations were carried out at 298 K and 35 bar while for hydrogen we used 77 K and 35 bar.</p><p>Using an N2 probe with 1.86 Å radius, probe occupied accessible volume (POAV) and accessible surface area (ASA) values of deformed/undeformed MOFs were computed using Zeo++ software [52,53].</p><p>The 2 × 2 × 2 supercell of MOFs was chosen for all studies except for the thermal expansion behaviors in which 1 × 1 × 1 was used to reduce computational cost.</p><!><p>This paper is organized in the following sections: First, we construct a single NNP for IRMOF-n series (n=1,4,6,7,10). Next, we test the NNP for both fragment and periodic bulk MOF structures.</p><p>Finally, we study the application of this NNP for thermal and mechanistic behavior of selected MOFs.  For training data set and construction of NNP, calculation of energies and forces of bulk structures at the DFT level is computationally too costly. Therefore, to reduce cost, we used the idea of fragmentation as introduced by Behler et al. for construction and development of the dataset [34].</p><!><p>Since the NN uses atomic contributions, the environment of atoms must be clearly defined within a certain cutoff. Thus, the MOF structure has been divided into pieces (fragments) by finding chemically inequivalent atoms that form the bulk structure. Similar to Behler's work, we generated the different possible fragments so that each fragment includes all different atom types and their neighbor atoms within a sphere of 6 Å cutoff radius. This cutoff distance was also used by the interaction block. In creating MOF fragments, cutting from cyclic structures was avoided. In addition, the missing bonds were saturated with hydrogen to achieve charge balance.</p><p>We started with the fragments of IRMOF-1 as it is the base for all other MOFs and thus the fragments of IRMOF-1 are common in all other MOFs studied. We defined IRMOF-1 with 5 fragments. Other MOFs differ from IRMOF-1 by only an additional linker. The organic linker in IRMOF-1 is 1,4-benzenedicarboxylate (BDC) while BDC links with n-propoxy, cyclobutyl and fused benzene functional groups, in IRMOF-4, IRMOF-6 and IRMOF-7, respectively. The linker is 4,4'-biphenyldicarboxylate (BPDC) in IRMOF-10. Therefore, an additional fragment for each of the other IRMOF-n (n=4, 6, 7, and10) was added to the five fragments of IRMOF-1 with a total of 9 fragments that cover all MOFs studied. The common five fragments for all MOFs and additional differing fragments belonging to each MOFs are present in Figure S1 (Supplementary Information).</p><p>An initial set of 6750 reference fragments were generated by the following way: First, all the bonds of each fragment (total 9 fragments) were scaled from 0.95 to 1.10 with 0.0002 steps size and 750 x 9 structures with non-equilibrium geometries were obtained. Then, the atomic positions were distorted randomly with a maximum displacement of 0.16 Å in the x, y, and z directions.</p><!><p>By performing single point energy DFT calculations for these fragments, the energies and forces were obtained in order to create the initial training dataset. The primitive NNP1 and NPP2 models were trained with this initial data set. The training is continued until ensuring that these two NNPs with different architectures produce nearly the same RMSE values on the same test dataset. To develop dataset, MD simulations of all fragments at different temperatures ranging from 50 to 600 K were carried out in the NVT ensemble using trial NNP1 as the energy and force calculator.</p><p>Since the configurational space is not adequately represented with initial dataset, missing structures from well-represented regions of configurational space must be added in order to develop the NNPs. Thus, wisely chosen new geometries were added to the training set as reference structures. We applied two different criteria to determine these new reference structures from the trajectories sampled from MD simulations. First, if a structure is not inside the training set, the prediction would be purely based on extrapolation and not reliable. Thus, the first criterion is that the structure must be included if the atom centered symmetry function (ACSF) values are outside the training set range. This criterion was carried out by adopting the idea of symmetry function values reported by Behler [41,54] although we used SchNet rather than HDNN. This is because the symmetry function vectors are translationally and rotationally invariant and it is much easier to compare different structures whether they are inside or outside the training set. Second criterion is that the energy produced by an NNP of different architecture (i.e. NNP2), which is trained by the same dataset, should yield similar result. Therefore, if the energy difference for two independent NNPs (Uncertainty=ΔENNP1-NNP2) is above a certain threshold value, that structure must be included in the training set. Starting with 0.5 eV, this threshold was gradually lowered to 0.05 eV, as the data set develops. The error will be larger in the NNP trained in the early data set, and this error rate decreases as the data set develops. Uncertainty improvement (i.e. RMSE/MAE) between two NNPs by addition of new configurations to the training set shows an exponential decay function (Figure S2). When the new structures based on these criteria have reached to a sufficient number (~5000 data points), these new structures were calculated by DFT and added to the reference set, which eventually created the next generation of NNPs. These processes were repeated several times until smooth NVT MD simulations for 1 ns could be run on the periodic bulk MOF structures.</p><p>Finally, ~170k reference data points were obtained when the NNPs were converged. The data were split into training (~126k, 80%) and validation (~17k, 10%) sets while the remainder were reserved for tests (~17k, 10%). The RMSE values on energies and forces were selected as metrics for the losses in training and validation. The loss consisted of weighted average RMSE of 99% forces and 1% energies (RMSELoss= RMSEenergy × 0.01 + RMSEforce × 0.99). Early stopping was carried out when the loss of the training and validation curves reached the plateau (Figure S3)</p><!><p>After NNP construction, energies and forces were tested on the fragments. The maximum errors for the energies are 0.036 eV/atom and 0.026 eV/atom for the training and test sets, respectively.</p><p>Figure 2 shows the plot of absolute errors between the NNP and DFT energies with respect to DFT energies for the training set and the test sets. The errors on most fragment structures have been estimated to lie within the narrow range of 0.005 eV/atom for both training and test sets. The vertical bar shows the counts of structures for the corresponding error range. Only 1.91% in the training set and 1.58% in the test set is estimated with error greater than 0.005 eV/atom. As the structures approach to the equilibrium geometries, the errors shrink whereas large errors are concentrated in geometries that are too far from the equilibrium, higher-lying energies on PES.</p><p>The distribution of the errors is random and fits to a normal Gaussian function (Figure S4).</p><!><p>Apart from testing on fragments, 1 ns NVT MD simulations at 300 K on each periodic bulk IRMOF-1,4,6,7,10 structures were also run following by energy minimization in order to test the stability of NNP. All the configurations visited by the MD simulations were in the well-represented regions of NNP. No structural changes apart from subtle thermal fluctuations were observed during the simulations.</p><!><p>The equilibrium geometries of fragments in bulk MOFs were not included in the development of NNP in order to avoid memorization of those values. Thus, reproducing the equilibrium lattice constants from the final NNP is one of the most critical test parameters for the success of NNP.</p><p>Table 1 shows the comparison of equilibrium lattice constants produced by NNP of bulk IRMOF structures with experimental single crystal x-ray diffraction (SCXRD) values at 100 K [55] along with literature calculations. Our developed NNP was utilized for carrying out NPT MD simulations at 100 K and 300 K of 500 ps to produce values of equilibrium lattice constants and bulk modulus.</p><p>The NNP calculated equilibrium lattice constants are in an excellent agreement with experiments and literature calculations, indicative of the reliability of the NNP. The difference between the values of NNP calculated and SCXRD ranges only from 0.04 Å to 0.64 Å for all MOFs trained. It should be stressed that NNP, which has been trained from non-equilibrium geometries of fragments, could successfully predicts the periodic equilibrium parameters. a: [28]; b: [56]; c: [57]; d: [58]; e: [59]; f: [60]; g:</p><!><p>Orientation of organic linkers is a key feature for the structural stability of IRMOF series and accurate calculation of torsional energies could be an important indicator for the success of the NNP [28]. The phenylene rings are common in all MOFs studied here. Therefore, as a case study, we compared NNP and DFT produced rotational energy barrier of a phenylene ring in IRMOF-1 to assess the reliability of the NNP. In order to reduce cost of the DFT calculation, we truncated the IRMOF-1 structure and saturated it with hydrogens (Figure S8). NNP and DFT energy profiles were produced by scanning the dihedral angle of one of the phenylene rings. Figure 4a shows the energy profile along a 180⁰ rotation. The DFT and NNP calculated values of 0.58 and 0.65 eV, respectively, diverges only by 0.07 eV. This difference might arise from a single NNP representation of all studied MOFs rather than developing an NNP for each of MOFs. In addition, both methods somewhat overestimate the experimental energy barrier obtained by 2 H-NMR studies ( 0.490±0.087 eV) [61]. The overestimation could be due to the use of truncated structure rather than the periodic one.</p><p>Similar to phenylene ring torsions, the same truncated structure was also used to reproduce the aromatic C-H bond distance energy profile of the phenylene ring. Figure 4b shows the relative energies corresponding the C-H distance scan calculated by DFT and NNP. Although NNP was developed using non-equilibrium geometries, it can perfectly predict the DFT calculated equilibrium bond distance of C-H to be 1.09 Å along the scan coordinate. DFT calculated distance was totally reproduced by NNP. Both methods agree reasonably with the experimental equilibrium aromatic C-H distance of 1.10846 Å by neutron powder diffraction (NPD) measurements [62].  [61] and NPD for equilibrium distance [62].</p><!><p>Reproducing phonon modes of a periodic bulk structure from fragment trained NNP is another critical test for the reliability of the potential. Calculating phonon modes directly from first principles can become quite costly for such large MOF systems. However, using the NNP, we utilized the calculation of vibrational modes available in ASE with much lower computational cost.</p><p>Experimental studies on IRMOF-1 vibrations showed that there are several spectral regions [63].</p><p>The dominant vibrations are C-H stretching (due to aromatic C-H on phenylene ring) at 2900-3100 cm -1 ; C=O stretching (due to carboxylate linker groups), C=C stretching (due to aromatic C-C on phenylene ring) at 1300-1600 cm -1 ; in-plane and out-plane deformation of aromatic ring (due to the phenylene ring) at 600-1200 cm -1 ; the Zn involved bands (due to inorganic subunit Zn-O stretches and bends) at 200-600 cm -1 ; and collective vibrations and lattice modes below 200 cm -1 [63]. The latter (also reported as below 23 meV) are typical for most MOFs and play important role in negative thermal expansion properties.</p><p>IRMOF-1 shows aromatic C-H symmetric and antisymmetric stretching bands at the 3126.4 cm -1 (23) and 3138.9 cm -1 (46), respectively (Figure 5 and Figure S9). These values are very consistent with previous DFT studies. Phonon modes of simplified model system of IRMOF-1 were reported in literature using two different functionals, B3LYP [64] and PBEsol-D3 [63] as 3230/3243 cm -1 and 3110/3122 cm -1 . Since our NNP was trained with PBE-D4, our values are in excellent agreement with PBEsol-D3 values. Not only the frequency values, but the intensities and splitting of 12 cm -1 between the symmetric and antisymmetric stretches were also totally reproduced by our NNP. The C-H stretches were also reported by HDNNP [28], which was also trained by PBE-D3, as 3104/3148 cm -1 . Our NNP calculated C-H stretching frequencies are also in agreement with previously reported experimental values. IRMOF-10 shows very similar C-H stretching frequencies to IRMOF-1 since all C-H bonds are aromatic (Figure S8). However, due to unequal chemical environments, C-H bonds at the interface of the rings slightly differ from those that are towards carboxylate groups. The C-H stretches in IRMOF-4 are more complicated due to the existence of aliphatic C-H bonds of n-propoxy groups linked in BDC. Similarly, additional aliphatic C-H stretches are observed in IRMOF-6 due to cyclobutyl group linked in BDC. In IRMOF-7, the C-H stretches of fused benzene ring are observed at slightly different frequencies than those on BDC. This difference can be rationalized by the fact that H atoms on the fused benzene are closer to the carboxylate oxygen atoms whereas the H atoms on BDC are more distant from carboxylate groups. Indeed, two isomeric forms as "twisted" and "planar" can conceivably exist in IRMOF-1 and IRMOF-7 and the stability of one of these forms were dominated by the steric repulsion of the fused benzene ring [65].</p><p>The second region of IRMOF-1 consists of the spectrum lying from 1514 cm -1 to 1625 cm -1 with several peaks including the one at 1612 cm -1 and 1595 cm -1 , corresponding to symmetric and antisymmetric stretches of aromatic carbon double bonds and carboxylate groups. This is in complete agreement with previously calculated DFT values of 1649 cm -1 and 1549 cm -1 and experimental IR peaks at 1612 cm -1 [55] and Raman peak at 1600 cm -1 . In IRMOF-7, these bands are calculated to be red shifted. For IRMOF-7, fused benzene rings' hydrogens that are proximal to the carboxylate oxygen atoms show more intense C-H bending vibrations.</p><p>Experimentally observed Zn-O (the carboxylate oxygen) stretch in IRMOF-1 at 578 cm -1 were also reproduced by the peaks in the range 582-612 cm -1 . The same mode was reported as 579-606 cm -1 in DFT calculations previously. In addition, the Zn-O (inorganic oxygen coordinated by four Zn atoms) is calculated to be 518-523 cm -1 . Experimentally, this peak is assigned to be 524 cm -1   while the DFT calculations reported the value of 512 cm -1 . Finally, Zn involved collective motions were reported to be below 287 cm -1 . Our calculations show 293 cm -1 and below are the corresponding motions including Zn-Zn vibrations. Since the additional fragments on the BDC linker increase the degree of freedom, more accessible vibrational states exist in lower lying spectral regions.</p><p>Overall, our calculated vibrational spectra illustrate a good agreement with those reported experimentally and those predicted by DFT density of state calculations.</p><!><p>In terms of thermal behaviors of pure MOF crystals, we have calculated thermal expansion coefficient, thermal pressure coefficient and bulk modulus. Using NNP for energy and force calculations, several MD simulations with temperature ranging from 200 K to 400 K were performed at the NPT ensemble. The mechanical pressure was set to 1 atm, and the data of MD simulations with five different temperature values from 200 to 400 K in steps of 50 K were fitted the equation ( 1 ) to find thermal expansion coefficient values. Similarly, thermal pressure coefficients and coefficient factors have been calculated from equations ( 2 ) and ( 3 ). Table -3 shows the comparison of these parameters calculated using NNP with literature. Using NNP, we find the thermal expansion coefficient values of IRMOF-1 and IRMOF-10 to be 47 MK -1 and 54 MK -1 , respectively. In literature, there is no consensus about these values. Experiments suggest two values to be -48 [20] and -39 MK -1 [66] while predicted values from simulations in literature range from -9 to -79 MK -1 [18,57,59,60,[67][68][69][70]. Our result for IRMOF-1 is very close to the one of the two experimentally reported data of -48 and -39 MK -1 . This might suggest the superiority of NNP over the classical simulations. Theoretical calculations on IRMOF-10 have also suggested different values ranging between -16 and -55 MK -1 [59,60,67,71], which is consistent with our prediction. To the best of our knowledge, the values of IRMOF-4, IRMOF-6 and IRMOF-7 have not been reported in literature. We predicted the value of IRMOF-6 as -35 MK -1 . However, we could not get a linear correlation for IRMOF-4 and IRMOF-7.</p><p>The origin of negative thermal expansion properties of MOFs is attributed to a few different scenarios. Shortening the bond lengths and phase changes, electronic effects, magnetorestriction and vibrational modes of bridging atoms. Bridging atoms can lead to high NTE with more than one contribution. First, the transverse vibrational modes on the bridging phenylene rings cause contraction of the Zn4O nodes. As temperature increases, the contraction gets larger, resulting with volume shrinking in MOF [28]. Second, the flexible carboxylate groups' vibrational modes lead shortening of the distance between the carboxylate groups and rigid O atoms of Zn4O [20].</p><p>The torsional motion of the phenylene ring may also contribute higher NTE. In particular, IRMOF-7 has fused benzene and IRMOF-4 has n-propoxy linked in BDE. Independent torsion of BDC linker proximate these functional groups in IRMOF-4 and IRMOF-7, resulting with greater energy barrier when compared to IRMOF-1, IRMOF-6, and IRMOF-10. This might explain unusual thermal behavior of the former two MOFs. At low temperatures, due to the greater torsional barrier, IRMOF-4 and IRMOF-6 first expand with temperature, i.e. positive thermal expansion (PTE). After reaching to a thermal energy which can compensate this barrier, the MOFs start contracting with increasing temperature, i.e., negative thermal expansion (NTE). In order to justify this hypothesis, we have scanned the rotation of two of the phenylene rings that are facing to each other in IRMOF-1 and IRMOF-7.</p><p>IRMOF-1 and IRMOF-7 structures are similar with some distinct differences. In IRMOF-1 structure, carboxylate groups form two planes that are coplanar to each other (0⁰). Although these carboxylate planes can be twisted (90⁰) in principle, they cannot be interconverted without breaking bonds of the network. In addition, the phenylene ring is also coplanar (0⁰) with these carboxylate groups due to the conjugation between the aromatic system and the carboxylate groups. In addition, the attraction between the proximal H and carboxylate O atoms is dominant in the coplanar carboxylate and phenylene groups since the minimum distance of 2.4 Å is reached. [65]. On the other hand, IRMOF-7 has fused benzene ring attached (ring B) on the phenylene ring (ring A) of the BDC linker (Figure 6). The two carboxylate groups were experimentally observed in the twisted (90⁰) form in IRMOF-7 [55]. For the naphthalene being coplanar with either one of the carboxylate planes would cause steric clashes between the proximal H atoms on the fused ring (ring B) with the O atoms of carboxylate groups. That is, the distance between proximal H atom on ring B and carboxylate O atom can be as low as 1.7Å, in which repulsions dominates. Thus, the equilibrium geometry of IRMOF-7 has the angle between naphthalene plane and carboxylate groups of 45⁰ in which the distance between proximal H of ring B and carboxylate O is ~2.3Å, very close to the distance between H atom of the phenylene ring (ring A) and carboxylate O in the IRMOF-1 with 2.4 Å. In addition, similar to IRMOF-1, H atom of ring A in IRMOF-7 has a distance of 2.4A with carboxylate O, which also contributes the stabilization. Steric hindrance is not limited only to the distance between the proximal H of ring B and carboxylate O but also exists between two neighboring rings. The distance between proximal H atoms between two neighboring B rings can be as low as 1.0 Å (=135⁰) in the case of a torsion of the naphthalene group (Figure 6). In order to monitor the effect of steric hindrance from these two scenarios, a rigid scan has been carried out by changing the angle between the phenylene ring and carboxylate planes on a single unit cell of the MOF using NNP in a two dimensional (2D) scan of two facing rings. Figure 7 shows the two dimensional potential energy surface (2D-PES) of IRMOF-1 and IRMOF-7.</p><p>In the case of IRMOF-1, an independent rotation around two facing phenylene rings could bring atoms to a distance that is no less than 2.4 Å (==90⁰) with an energy barrier up to ~1.2 eV. For a single dimensional scan, this barrier is almost halved to ~0.6 eV (=0⁰, =90⁰ or vice versa), which is in perfect agreement with experimental and DFT calculated values. On the other hand, due to the neighboring fused benzene rings, the atoms in IRMOF-7 can be as proximate as 1.0 Å, which gives rise to the energy barrier as much as ~6 eV (==135⁰). For a single dimensional scan, this barrier is almost halved to ~3 eV (=45⁰, =135⁰ or vice versa). Additionally, another energy barrier (with double shoulders) of ~2.5 eV (translation to 1D is ~0.6 eV) has been observed from the fused benzene H and carboxylate repulsions with 1.7 Å distance (=45⁰, =0⁰ and 90⁰). Switching from PTE to NTE were first reported in literature and this "switching temperature" were determined as 230 K for DUT-49 [19]. Although experimental observation of this phenomenon required solvated systems, underlying molecular basis can still be investigated for unsolvated systems by theoretical calculations. It should be noted that the thermal expansion coefficients of MOFs with dual thermal expansion characteristics cannot be simply produced by linear ln(<V>) vs. T plot in equation ( 2 ) since α0 is not constant and thus should be reported as a function of temperature (For the details, refer to the discussion in Supplementary Information). shown to be larger than that of IRMOF-1 in literature both experimentally [68] and from simulations [18]. Our calculations agree well with experiments and computations in terms of the magnitude and the decrease of 𝛼 0 from IRMOF-10 to IRMOF-1.</p><p>After reaching to the NTE region, IRMOF-4 and IRMOF-7 contract more than IRMOF-1 with increasing temperature. These higher NTE values can be explained by greater repulsion of the phenylene rings due to greater energy barrier discussed earlier.</p><p>We also calculated thermal pressure coefficients from MD simulations using NNP (Table S3).</p><p>Thermal pressure is an indicator for the rigidity and stability of MOFs. Thermal pressure coefficient, on the other hand, indicates how the pressure would change when a material in a closed system is subjected to temperature increase. Among the IRMOF series studied, IRMOF-10 shows the lower thermal pressure coefficient than IRMOF-1, which is qualitatively in agreement with simulations in literature [18].</p><!><p>In addition to thermal behaviors, after developing NNP and testing it with fragments and bulk materials, we have also investigated elastic response of IRMOFs to applied mechanical pressure.</p><p>Pressure can also be utilized to control gas storage and release from MOFs. Thus, we have explored methane and hydrogen gas storage properties of IRMOF-1 and IRMOF-7 under applied external pressure. To do that, we have performed MD simulations using our NNP of uniaxial compression by reducing one of the direction vector in several steps. At each step, MD frames were extracted for the calculation of density, accessible surface area (ASA), probe occupied pore volume (POAV) and POAV fraction using 1.86 Å probe (N2) by Zeo++ [52,53]. As expected, the stress applied increases the density of both MOFs smoothly while reducing POAV fraction. Reducing volume also causes the compression of the organic linkers (Figure 8). Due to the steric hindrance of fused benzene rings, the decrease of POAV fraction in IRMOF-7 is predicted to be less than IRMOF1.</p><p>The CH4 uptakes of IRMOF-1 under mechanical compression has been reported by Zheng et al [27]. Using similar protocol, we have also studied the effect of mechanical compression on gas storage ability of IRMOF-1 and IRMOF-7. We have computed the CH4 and H2 uptakes of both compressed and uncompressed structures that are retrieved from strain applied MD simulations using NNP. Our calculations for loading of CH4 on IRMOF-1 is in a good agreement with literature.</p><p>Figure 9 shows the CH4 and H2 uptakes of IRMOF-1 and IRMOF-7 against the strain applied on MOFs. Gravimetric methane uptake starts decreasing once the MOF starts deforming. Since the accessible surface area (ASA) and probe occupied accessible volume (POAV) decreases by the bending of organic linkers upon the stress applied, the total gravimetric uptake decreases (Figure S10). The fused benzene rings on IRMOF-7 results in less ASA and POAV.</p>