Thermal transport mechanisms in ZIFs

Domain: physics
Published: 2025-11-27
Editorial outcome: minor_revision

Abstract
Zeolitic imidazolate frameworks (ZIFs), a subclass of metal-organic frameworks (MOFs), exhibit tunable thermal conductivity, which is crucial for applications such as gas adsorption, catalysis, and energy storage. Despite its significance, thermal conductivity in MOFs remains less explored compared to other key performance indicators. In this work, we investigate the thermal transport properties of 196 ZIF structures with diverse net topologies and organic linkers. We develop a thermal circuit model that quantitatively integrates network topology and consolidates various atomic contributions into heat conduction units for thermal analysis. Our results reveal a strong correlation between circuit-estimated and simulated thermal conductivity, demonstrating the model’s predictive power. Additionally, we find that functional groups influence thermal transport through a competing interplay between atomic mass and mechanical stability. These findings provide a systematic approach for tailoring MOF thermal properties, offering insights into the rational design of materials with optimized thermal performance.

Introduction
Metal-organic frameworks (MOFs) are a class of crystalline materials composed of metal clusters and organic ligands. Their structural tunability and high porosity make them promising for various applications1 such as gas adsorption and separation2–4, catalysis5,6, and drug delivery7–9. Researchers have developed various key performance indicators (KPIs) to assess MOFs for specific applications and have extensively studied their potential10,11. New materials are designed using advanced machine-learning techniques based on these KPIs12. However, thermal conductivity is often overlooked despite its potential to impact the MOFs’ ranking. For example, poor heat dissipation within the structure can reduce the amount of gas adsorbed13–15. Similarly, the rate of heat conduction in MOFs can affect reaction rates and selectivity in catalytic processes16.

Recent studies have increasingly focused on the thermal conductivity of MOFs. Islamov et al.17 conducted high-throughput molecular simulations on 10194 hypothetical MOFs generated by ToBaCCo, revealing that high thermal conductivity is associated with high density, small pores, and four-coordinated metal nodes. Similarly, Lin et al.18 found that close metal atom arrangements and a high proportion of metal and transition metal atoms contribute to higher thermal conductivity in their simulations of 1214 structures from the CoRE MOF database. Comparable correlations between thermal conductivity and geometric characteristics, such as pore size19–21 and density22, have also been extensively investigated in covalent organic frameworks (COFs). Due to the more uniform covalent bonding environments in COFs, these structure-property relationships generally exhibit stronger consistency. While these high-throughput simulations18,21 have identified some common trends, such as the influence of pore size and metal connections17,18, other factors like mass mismatch17,23, and the impact of functional groups24, the net topology25, and pore shape26 remain less understood.

Zeolitic imidazolate frameworks (ZIFs) are a distinct subclass of MOFs, consisting of tetrahedrally coordinated transition metal ions linked by imidazolate (IM) ligands. Ying et al.24 studied ZIF-8 with four different functional groups and found that thermal conductivity κ followed the trend: κ(ZIF-8(CH3)) > κ(ZIF-8(H)) > κ(ZIF-8(Br)) > κ(ZIF-8(Cl)). To gain broader insights into the structure-property relationship, Cheng et al.25 investigated 18 ZIFs with the same organic ligands (2-methyl imidazolate (mIM)) but varying net topologies. They introduced two parameters: alignment tensor and pathway factor, to quantitatively evaluate the orientation and distribution of heat transfer pathways within frameworks. However, these parameters are not easily generalizable to a broader range of MOFs18.

In MOFs, we typically have three parameters we can tune: the metal node, linker, and network topology. In this work, we focus on both the network topology and functionalization of the linker. Experimentally, for ZIFs only, we have been able to obtain a large number of ZIFs with the same metal node and linker but with different network topologies27. We expanded the understanding of thermal transport in ZIFs by simulating 196 structures with various net topologies and organic linkers. We developed a thermal circuit model to analyze the net topology of a structure in heat conduction quantitatively. Furthermore, we highlight how the mechanical stability changes and mass variations introduced by functional groups jointly influence thermal transport. In principle, similar studies can be done with MOFs, but there are far fewer experimental examples in which we can synthesize MOFs with the same linker and metal node but different pore topologies28.

Results
Thermal conductivity of ZIFs
The ZIFs studied in this work are the same set that was used by Moosavi et al.29 to investigate the mechanical properties. To ensure the chemical consistency of all 196 ZIFs, we have run mofchecker30, and no issues were detected. We performed equilibrium molecular dynamics (MD) simulations for 196 ZIFs with 49 different net topologies and four different organic ligands. The organic ligands include IM, dichloro imidazolate (dcIM), mIM, and 2-nitro imidazolate (nIM), whose structures are shown in Fig. 1a. Each zinc atom is tetrahedrally coordinated with four organic ligands. The results of these simulations are summarized in Fig. 1b, c. The numerical results can be found on Zenodo31.

Previous studies have identified a linear relationship between thermal conductivity κ and geometric properties such as density13,25 and pore size25,26. More recently, Islamov et al.17 expanded on this finding by demonstrating a general trend between density and thermal conductivity through simulations of a large, diverse MOF dataset. Figure 1b shows that for the ZIFs we studied, the thermal conductivity follows a linear relationship against density for each metal-ligand combination (distinguished by color). However, no overall correlation is found when considering all data points together (Pearson coefficient r = −0.37). An increase in density does not always correspond to enhanced thermal conductivity. A similar analysis of the correlation between void fraction and thermal conductivity can be found in Section 1 of Supplementary Information (SI).

Figure 1c shows thermal conductivity as a function of volumetric atom density, resulting in an improved overall Pearson coefficient of 0.78. Despite this improvement compared to density, volumetric atom density may not be an ideal indicator. Indeed, Fig. 1c suggests that atoms in the backbone and those in functional groups (e.g., Cl in dcIM, C and H in mIM, and N and O in nIM) likely contribute differently to thermal conductivity.

The unsatisfactory correlations between density, volumetric atom density, and thermal conductivity can be attributed to two main factors. First, from our studies, the contribution of each atomic species to thermal conductivity depends on its mass and chemical environment. Density and volumetric atom density fail to capture the varying effects. Second, these simple density-based correlations do not explicitly account for the network and connectivity within the structures. These observations motivated the search for alternative indicators that consider these factors, with the potential to explain thermal conductivity more effectively.

Thermal circuit model
Heat conduction is analogous to electrical conduction: the former is driven by temperature gradients, while potential gradients drive the latter. This analogy has been exploited to study thermal transport in continuum materials32,33. However, its application in simulating heat transfer within molecules remains unexplored. By transforming a crystal network into a thermal circuit network, we calculate thermal conductivity using an analogy to electrical conductivity. In the remainder, we refer to this quantity as “thermal circuit conductivity”.

An example of a two-dimensional periodic crystal network is shown in Fig. 2a. In this representation, nodes correspond to atoms (or groups of atoms) located at the intersections of heat conduction paths, while edges denote the connecting atoms and chemical bonds. A temperature difference, e.g., T1 on the left and T2 on the right side of a unit cell, drives the heat flux. By simulating the resulting heat flow, the thermal conductivity of the crystal network can be determined. To model this, we use an electricity analogy, mapping the crystal structure to a circuit. On the left boundary of the unit cell, the nodes with bonds crossing the left unit cell boundary are selected as starting points, and their equivalent nodes in the adjacent unit cell serve as endpoints. These starting and endpoint nodes are treated as equipotential. Applying a voltage U across them creates a uniform electric field Ei = U/Li, corresponding to a constant potential gradient along the unit cell direction Li. The current density is given by Ji = I/A⊥i, where A⊥i is the cross-sectional area perpendicular to the transport direction. Figure 2b shows the circuit corresponding to the two-dimensional (2D) crystal network. In this illustrative case, with U = 1 V and each resistor set to a resistance of R = 1 Ω, the total current is I = 15 A. In the general three-dimensional case, Ohm’s law yields the circuit conductivity as

σi = Ji / Ei = (I/A⊥i) / (U/Li) = (I/U) * (Li/A⊥i), i = x, y, z (1)

In the study of ZIFs, we construct the circuit model by defining the zinc atoms as the network nodes. A resistor is placed between two nodes when they are connected by an organic linker, representing a heat conduction unit (Zn-linker-Zn, where the linker is IM, dcIM, mIM, or nIM), as shown in the inset of Fig. 2a. The corresponding circuit is generated using an in-house Python code31.

Since the resistance of an individual heat conduction unit is not known a priori, we employ a reduced (dimensionless) formulation. Specifically, we set U* = 1 and assign a resistance of unity to each heat conduction unit. The reduced current I* is then obtained by circuit simulation in MATLAB Simulink34. The reduced circuit conductivity σ*i is calculated as

σ*i = (I* / U*) * (L*i / A*⊥i), i = x, y, z (2)

where A*⊥i is the reduced cross-sectional area of the structure, L*i is the reduced lattice length in the computed direction. The reduced thermal circuit conductivity τ* equals the value of σ*i. Equation (2) gives us the equivalent reduced circuit conductivity for each lattice direction by constructing the circuit along the respective direction. The overall reduced thermal circuit conductivity is calculated as the average conductivity across x, y, and z directions. The correlation between reduced thermal circuit conductivity τ* and MD thermal conductivity κ is shown in Fig. S7 in the SI.

Once we have a set of data for the same type of heat conduction units, namely, reduced thermal circuit conductivity (τ*) and MD thermal conductivity (κ), we can statistically estimate the thermal resistance of the heat conduction unit. Specifically, we fitted a linear model between the reduced thermal circuit conductivity and the MD thermal conductivity for each group of ZIFs with the same linker. Using the linker-specific scaling factors obtained from the fitting, we can convert the reduced thermal circuit conductivity (τ*) to a thermal circuit conductivity (τ) in real units. Details about the fitting process and the scaling factors can be found in Section 7 of SI.

With these scaling factors, the thermal circuit conductivities (τ) of the four types of ZIFs converge onto a single universal linear correlation with the MD thermal conductivity. The thermal circuit conductivity exhibits a strong and universal correlation with the MD thermal conductivity, with a Pearson correlation coefficient of 0.92, as shown in Fig. 3. By treating atoms in different chemical environments as collective heat conduction units, the circuit model accounts for varying atomic contributions to thermal transport. This approach outperforms density (r = −0.37) and volumetric atom density (r = 0.78).

Effect of functional groups
Figure 1b shows that ZIFs with the same net topology but heavier atoms can exhibit higher density yet lower thermal conductivity. For example, replacing imidazolate with 2-nitroimidazolate increases the framework density but lowers the heat conductivity. In contrast, volumetric atom density cannot distinguish ZIFs with different thermal conductivity. For instance, ZIFs(dcIM) exhibit lower thermal conductivity despite having a similar volumetric atom density to ZIFs(IM). By considering collective atoms as a heat conduction unit, the true effect of adding a functional group becomes evident in Fig. 3. Functionalization has minimal impact on ZIFs(nIM), enhances thermal conductivity in ZIFs(mIM), and reduces it in ZIFs(dcIM). This trend, κ(mIM) > κ(nIM) > κ(IM) > κ(dcIM), is also reflected in the fitted scaling factors (see Table S2 in the SI).

If we add a functional group, two things can happen. First, the mass of the ligand changes, and it is well-known that differences in mass can have a large effect on thermal conductivity. A second effect is that a functional group can impact the material’s stiffness29. Structural stiffness and thermal conductivity are inherently linked, as strong bonds suppress atomic vibrations and enhance phonon transport, leading to both high mechanical stability and increased thermal conductivity.

To investigate these effects, we first carry out a simulation in which we remove the mass difference by modifying the ligands’ atomic masses so that all four organic linkers have the same mass as the lightest linker, IM. This is not a simulation that mimics an experiment one can do in the real world. Eleven net topologies were randomly selected for these simulations. The results, shown as half-filled markers in Fig. 4, reveal a distinct trend: when the mass of the linker is the same, the introduction of functional groups increases thermal conductivity compared to ZIFs(IM). This effect is most pronounced in ZIFs(mIM) and ZIFs(nIM), where the thermal conductivity ratio exceeds 1, whereas ZIFs(dcIM) exhibit a more modest increase.

Now that we have (artificially) eliminated the effect of mass, the remainder of the effect can be explained in terms of changes in the rigidity of the network. Indeed, previous work on the mechanical stability of MOFs has highlighted a similar impact of functional groups on network stiffness. Moosavi et al.29 demonstrated that strategically placed functional groups can enhance mechanical stability by forming a secondary network of nonbonded interactions. Figure S3 in the SI, which is based on data from ref. 29, shows that ZIFs with mIM and nIM linkers have a higher bulk modulus, and thus greater mechanical stability, compared to those with IM and dcIM linkers.

In Section 5 in the SI, we further show that ZIFs with functional groups exhibit a higher contribution of van der Waals interactions to total thermal conductivity than their IM-linked counterparts. These results support the idea that enhanced mechanical stability, driven by non-bonded interactions, increases thermal conductivity.

Furthermore, our calculations reveal that mass plays a non-negligible role in heat conduction. The heaviest linker, dcIM, results in the greatest reduction in thermal conductivity, as evidenced by the large difference between the half-filled and solid markers for dcIM in Fig. 4. In contrast, for nIM, the competing effects of mass and mechanical stability approximately cancel each other out. In ZIFs(mIM), the dominant influence of enhanced mechanical stability leads to the highest observed thermal conductivity. This interplay between mass and mechanical stability leads to the trend: κ(mIM) > κ(nIM) > κ(IM) > κ(dcIM).

Two exceptions are observed in the analysis: BCT and SOD in ZIFs(mIM) and ZIFs(nIM), which can be attributed to their relatively high bulk modulus29 and similar zigzag heat conduction pathways (see Section 6 in the SI for more details). In addition, these two topologies are characterized by the presence of 4-membered and 6-membered rings, structural motifs that tend to enhance the vibrational freedom of the metal centers35,36. These additional vibrations may lead to variations in the thermal resistance of the corresponding heat conduction units.

Methods
We performed MD simulations using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software37. The heat flux was computed using the algorithm of Boone et al.38. The input files were generated using the lammps-interface package39. The frameworks were equilibrated with the isothermal-isobaric (NPT) ensemble at 300 K and 1 bar for 10 ns, followed by a canonical (NVT) ensemble at 300 K for another 10 ns. The Nosé-Hoover thermostat and barostat, respectively, control the temperature and pressure. The velocity Verlet algorithm is used to integrate the equation of motion with a time step of 1 fs. For structures for which it is difficult to reach equilibration, the time step was reduced to 0.5 fs in the equilibration stage. After setting the temperature, another 10 ns of equilibration was performed using a microcanonical ensemble (NVE). Then, we recorded the heat flux in each direction within the structure for around 50 ns with an NVE ensemble simulation. The thermal conductivity in each direction was computed from the corresponding Green-Kubo relation40. An automatic fitting procedure to compute the converged thermal conductivity is provided in Section 2 in the SI.

We investigated the effects of system size and force fields (FF) on thermal conductivity in Fig. 5. We find that smaller sizes yield lower thermal conductivity due to finite-size effects. Simulations using UFF4MOF with 2×2×2 supercells (small yellow circles in Fig. 5) return lower thermal conductivity values than those using UFF4MOF with 4×4×4 supercells (big yellow circles). Among the ZIFs studied, SOD(mIM) corresponds to the experimental structure ZIF-8, for which the thermal conductivity of a single crystal has been experimentally measured at 0.64 ± 0.09 W m−1 K−1 (the red area in Fig. 5)41. For ZIF-8, simulations with 4×4×4 supercells produce results closer to the experimental value. We conducted a size effect study for all the ZIFs in our studies and ensured convergence to the bulk thermal conductivity.

We tested two force fields, DREIDING42 and UFF4MOF43,44, and compared them with other FFs reported in the literature. Our results show that the DREIDING force field (cyan squares) predicts an unusually high thermal conductivity. Simulation using other FFs, such as MOF-FF24 and the flexible ZIF-8 force field45 (gray triangles), yield even lower values than our UFF4MOF simulations with 2×2×2 supercells. While including Coulombic interactions in UFF4MOF (orange rhombus) improves consistency with experimental data, these interactions were excluded from our study due to their high computational cost. In addition, one advantage of UFF4MOF is its transferability, as it can be used for various MOFs.

We chose to use UFF4MOF and a sufficiently large system because the resulting combinations strikes a good balance between accuracy in matching the experimental results and computational costs.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors developed a current circuit model to predict the thermal conductivity of 196 ZIF structures and a\nstrong correlation between circuit model predicted and MD simulated thermal conductivity was revealed. In addition,it is also\ndemonstrated that the impacts of functional groups on ZIFs are through the interplay between atomic mass and mechanical\nstability based on chosen ZIF structures.In general, this is an interesting work and the manuscript is well-written. However,\nsome critical issues should be addressed before publication.\n1. It is known that the thermal conductivity of MOFs/ZIFs are closely related with the crystal structure. What are the criteria of\nchoosing the 196 ZIFs chosen for this study? How did the authors to guarantee the correctness of their ZIF structures. Since\nit has been reported previously that many structures in CORE_MOF have some problematic structures, e.g, some atoms of\nMOFs have unreasonable oxidation state numbers.\n2. The current circuit model is the key of this study. However, the details on how to predict the thermal conductivity by such a\nmodel was not provided in the manuscript. Although Eq.1 is provided, but the detailed approach on how to obtain each\nparameter should be given in the main text.\n3. According to Fig. 3, it is clear that the circuit model predicted thermal conductivity is generally higher than MD thermal\nconductivity. Is there any reason for such a phenomenon? Besides, for some ZIFs with high thermal conductivity above 0.5\nW/mK, their thermal conductivity is significantly underestimated. Please explain on such a trend.\n4. Compared with experimental measured thermal conductivity, which one is more accurate? Such discussion should be\nalso supplemented. The universality of such current circuit model should be demonstrated by predicting some other MOFs\nwith known experimental measured and MD simulated thermal conductivities.\n5. Regarding the impacts of functional groups, the authors could provide the results from circuit model and MD simulation to\nvalidate the applicability of circuit model.\n6. In Fig.4, the reduced thermal conductivity of dcIM-based ZIF was attributed to the big mass difference from IM. However,\nnIM-based ZIF exhibited a bigger mass difference from IM, but their thermal conductivity is higher. In addition, the\nmechanical stability of such ZIFs should be also supplemented to demonstrate the interplay between atomic mass and\nmechanical stability.", "reviewer_2": "This manuscript presented by Zhang and co-author explores the thermal transport mechanisms in ZIFs. Generally,\nunderstanding the insight influence of MOF structure on thermal properties is quite important for MOF design and thermal\nengineering. Herein, they develop a circuit model to calculate the reduced thermal circuit conductivity and reveal its\ncorrelation with thermal conductivity. I recommend this manuscript be published after clarifying the following issues.\n1. As noted in the introduction, the developed parameter from limited structures illustrates the structure-property relationship\nhardly generalizable to a broader range of MOFs. How about the generalizable performance of such a circuit model applied\nto other MOFs except these 196 ZIFs?\n2. The heat resistance network is well-established in thermal engineering. The difference in such reduced thermal circuit\nconductivity is the introduced parameters of Li and Ai. What’s the insight meaning of such parameters' influence on\nthermal conductivity? Without such parameters, how about the updated correlation results?\n3. For the MOF crystal, the phonon is leading the thermal transfer, how about the functional group influence on phonon\ntransport?\n4. As provided in Figure S3, most ZIFs with the SOD and BCT net topologies exhibited relatively high bulk modulus (7-10\nGPa), which is fairly close to the range of ZIFs(mIM) and ZIFs(nIM) (6-10 Gpa). So, why the BCT and SOD topology were\nexcept to obtain the Pearson coefficient. How about the overall Pearson correlation with such two topologies?\n5. According to the definition, each chemical bond is modeled as a resistor with a 1 Ω resistance. However, the influence of\nvarious chemical bonds (i.e., Zn-N, C-C) on thermal transport should be quite different. The thermal resistance of various\nchemical bonds can be explored to better illustrate the correlation. For example, the values for such 4 functional groups\n(mIM, nIM, IM, and dcIM) can be determined."}, "author_response": "Response to the reviewers\nReviewer 1\nComments to the Author:\nInthiswork, theauthorsdevelopedacurrentcircuitmodeltopredictthethermalcon-\nductivity of 196 ZIF structures and a strong correlation between circuit model predicted\nand MD simulated thermal conductivity was revealed. In addition,it is also demonstrated\nthat the impacts of functional groups on ZIFs are through the interplay between atomic\nmass and mechanical stability based on chosen ZIF structures.In general, this is an in-\nteresting work and the manuscript is well-written. However, some critical issues should\nbe addressed before publication.\nReviewer Point P1.1 — It is known that the thermal conductivity of MOFs/ZIFs\nare closely related with the crystal structure. What are the criteria of choosing the 196\nZIFs chosen for this study? How did the authors to guarantee the correctness of their ZIF\nstructures. Since it has been reported previously that many structures in CoRE MOF\nhave some problematic structures, e.g, some atoms of MOFs have unreasonable oxidation\nstate numbers.\nReply: We used the same set as in a previous study. To clarify this, we added to the\ntext:\nThe ZIFs studied in this work are the same set that was used by Moosavi\net al. [1] to investigate the mechanical properties.\nAs an extra check of the correctness, we have run the new version of mofchecker and\nadded to the revised manuscript:\nToensurethechemicalconsistencyofall196ZIFs, wehaverunmofchecker[2]\nand no issues were detected.\nthe details on how to predict the thermal conductivity by such a model was not provided\nin the manuscript. Although Eq.1 is provided, but the detailed approach on how to\nobtain each parameter should be given in the main text.\nReply:\nWe revised the main text, which now reads:\nAn example of a two-dimensional periodic crystal network is shown in Fig. 2a.\nIn this representation, nodes correspond to atoms (or groups of atoms) lo-\ncated at the intersections of heat conduction paths, while edges denote the\nconnecting atoms and chemical bonds. A temperature difference, e.g., T on\n1\nthe left and T on the right side of a unit cell, drives the heat flux. By simu-\n2\nlating the resulting heat flow, the thermal conductivity of the crystal network\ncan be determined.\nTo model this, we use an electricity analogy, mapping the crystal structure\nto a circuit. On the left boundary of the unit cell, the nodes with bonds\ncrossing the left unit cell boundary are selected as starting points, and their\nequivalent nodes in the adjacent unit cell serve as endpoints.\nThese starting and endpoint nodes are treated as equipotential. Applying a\nvoltage U across them creates a uniform electric field E = U/L , correspond-\ni i\ning to a constant potential gradient along the unit cell direction L . The\ni\ncurrent density is given by J = I/A , where A is the cross-sectional area\ni ⊥i ⊥i\nperpendicular to the transport direction. Figure 2b shows the circuit corre-\nsponding to the 2D crystal network. In this illustrative case, with U = 1V\nand each resistor set to R = 1Ω, the total current is I = 15A. In the general\nthree-dimensional case, Ohm’s law yields the circuit conductivity as\nJ I/A\ni ⊥i\nσ = = ,i = x,y,z. (1)\ni\nE U/L\ni i\nIn the study of ZIFs, we construct the circuit model by defining the zinc\n2\natoms as the network nodes. A resistor is placed between two nodes when\nthey are connected by an organic linker, representing a heat conduction unit\n(Zn-linker-Zn, where the linker is IM, dcIM, mIM, or nIM), as shown in the\ninset of Fig. 2a. The corresponding circuit is generated using an in-house\nPython code[3].\nSince the resistance of an individual heat conduction unit is not known a\npriori, we employ a reduced (dimensionless) formulation. Specifically, we\nset U∗ = 1 and assign a resistance of unity to each heat conduction unit.\nThe reduced current I∗ is then obtained by circuit simulation in MATLAB\nSimulink[4]. The reduced circuit conductivity σ∗ is calculated as\ni\n...\nReviewer Point P1.3 — According to Fig. 3, it is clear that the circuit model\npredictedthermalconductivityisgenerallyhigherthanMDthermalconductivity. Isthere\nanyreasonforsuchaphenomenon? Besides,forsomeZIFswithhighthermalconductivity\nabove 0.5Wm−1K−1, their thermal conductivity is significantly underestimated. Please\nexplain on such a trend.\nReply:\nWe analyzed the deviation distribution and included it in the SI:\nFigure 1: The deviation distribution of the thermal circuit predictions. The deviation is\ndefined as the difference between the thermal conductivity of the thermal circuit and\nthe MD thermal conductivity.\nIn Fig. S8, we present the deviation of the thermal circuit predictions, de-\nfined as the difference between the thermal circuit conductivity τ and the\nMD thermal conductivity κ. Overall, the deviation distributions exhibit an\napproximately normal shape with only a slight rightward shift, indicating\nthat the thermal circuit model generally provides unbiased predictions. How-\never, a few outliers, such as the MD thermal conductivity values of BCT and\nSOD, exert a disproportionate influence on the fitted scaling factors τ .\nlinker\nThese outliers enlarge the fitted values of τ , thereby shifting the pre-\nlinker\ndicted conductivities upward. This effect is most pronounced in IM, where\nthe distribution shows a stronger positive bias, while in dcIM the influence of\nsuch exceptional values is more moderate, leading to only a minor shift.\nThe structures for which the model underestimates the thermal conductivity corre-\nspond to two specific topologies, as discussed at the end of Section 2.3:\nTwo exceptions are observed in the network analysis: BCT in the z-direction\nand isotropic SOD in x-, y-, and z-directions, which can be attributed to their\n4\nrelatively high bulk modulus[1] and similar zigzag heat conduction pathways\n(see Section 6 in the SI for more details).\nSubsequently, further discussion is added:\nIn addition, these two topologies are characterized by the presence of 4-\nmembered and 6-membered rings, structural motifs that tend to enhance the\nvibrational freedom of the metal centers.[5, 6] These additional vibrations\nmay lead to variations in the thermal resistance of the corresponding heat\nconduction units.\nReviewer Point P1.4 — Compared with experimental measured thermal conduc-\ntivity, which one is more accurate? Such discussion should be also supplemented. The\nuniversality of such current circuit model should be demonstrated by predicting some\nother MOFs with known experimental measured and MD simulated thermal conductivi-\nties.\nReply: As to the generalisation of the circuit model to other MOFs, we added the\nfollowing to the introduction to clarify why we focus on ZIFs:\nIn the introduction:\nIn MOFs, we typically have three parameters we can tune: the metal node,\nlinker, and network topology. In this work, we focus on both the network\ntopology and functionalization of the linker. Experimentally, for ZIFs only,\nwe have been able to obtain a large number of ZIFs with the same metal node\nand linker but with different network topologies.[7]\n...\nIn principle, similar studies can be done with MOFs, but there are far fewer\nexperimentalexamples in whichwecan synthesizeMOFswith the same linker\nand metal node but different pore topologies.[8]\nTo the best of our knowledge, there is only one experimental value available for the\nstudied ZIFs, and we compared our simulations against this value. The good agreement\nindicates that our simulation provides a reliable prediction.\n5\nand in Discussion:\nAlthough our work focuses on ZIFs, the proposed model provides a gener-\nalizable and efficient strategy that holds strong potential for extension to a\nbroader range of MOFs, offering guidance for tailoring thermal conductivity\nthrough topology design and functionalization.\nReviewer PointP1.5 — Regardingtheimpactsoffunctionalgroups,theauthorscould\nprovide the results from circuit model and MD simulation to validate the applicability of\ncircuit model.\nReply: The impacts of the functional groups on MD simulations and circuit model\npredictions follow the same trend, as shown in Fig. 3; we add the marginal distribution\nof both in Fig. 3 to clarify this:\nFigure 2: The MD thermal conductivity κ versus the thermal circuit conductivity τ.\nHistograms of τ and κ, distinguished by linker type, are shown along the x- and y-axes,\nrespectively. Thermal resistance coefficients are fitted for ZIFs with identical linkers,\nyielding good linear correlations when a single coefficient is applied. Pearson coefficients\nare calculated after excluding two exceptional net topologies, BCT and SOD, which are\nrepresented by squares.\nWith these scaling factors, the thermal circuit conductivities (τ) of the four\ntypes of ZIFs converge onto a single universal linear correlation with the MD\nthermal conductivity.\nReviewer Point P1.6 — In Fig.4, the reduced thermal conductivity of dcIM-based\nZIFwasattributedtothebigmassdifferencefromIM.However, nIM-basedZIFexhibited\na bigger mass difference from IM, but their thermal conductivity is higher. In addition,\nthe mechanical stability of such ZIFs should be also supplemented to demonstrate the\ninterplay between atomic mass and mechanical stability.\nReply: dcIM linkers have a larger mass compared to nIM linkers. Note that in Fig. 1a,\ndcIM has two hydrogen atoms substituted by chlorine atoms.\nPrevious work has studied the mechanical stability of the ZIFs. We wrote in the main\ntext:\nMoosavietal.[1]demonstratedthatstrategicallyplacedfunctionalgroupscan\nenhance mechanical stability by forming a secondary network of nonbonded\ninteractions. Figure S3 in the SI, which is based on data from Moosavi et al.\n[1], shows that ZIFs with mIM and nIM linkers have a higher bulk modulus,\nand thus greater mechanical stability, compared to those with IM and dcIM\nlinkers.\nComments to the Author: This manuscript presented by Zhang and co-author ex-\nplores the thermal transport mechanisms in ZIFs. Generally, understanding the insight\ninfluence of MOF structure on thermal properties is quite important for MOF design\nand thermal engineering. Herein, they develop a circuit model to calculate the reduced\nthermal circuit conductivity and reveal its correlation with thermal conductivity. I rec-\nommend this manuscript be published after clarifying the following issues.\nReviewer Point P2.1 — As noted in the introduction, the developed parameter from\nlimited structures illustrates the structure-property relationship hardly generalizable to\na broader range of MOFs. How about the generalizable performance of such a circuit\nmodel applied to other MOFs except these 196 ZIFs?\nReply: A similar comment was made by Reviewer 1, please see our reply P1.4.\nReviewer Point P2.2 — The heat resistance network is well-established in thermal\nengineering. The difference in such reduced thermal circuit conductivity is the introduced\nparameters of Li and Ai. What’s the insight meaning of such parameters’ influence on\nthermal conductivity? Without such parameters, how about the updated correlation\nresults?\nReply: First, asmentionedinthemaintext, thecircuitmodelhasbeenutilizedtoinves-\ntigate thermal transport in continuum materials. However, its application in simulating\nheat transfer within molecules remains unexplored, to the best of our knowledge.\nIn addition, L and A are not arbitrary fitting parameters but quantities that have\ni ⊥i\nphysical meanings. Specifically, L determines the potential gradient (U/L ) along the\ni i\ntransport direction, and A provides the effective cross-section through which thermal\n⊥i\ncurrentpasses. Withouttheseparameters, themodelwouldreducetoapurelytopological\nnetwork. We have revised our manuscript to clarify this. Please find the details in our\nreply to Reviewer 1, point P1.2.\ntransfer, how about the functional group influence on phonon transport?\nReply: PhonontransportinZIFshasalreadybeenextensivelyinvestigatedbyYingetal.\n[9], Zhou et al. [10]. We did not specifically study phonon transport, as our focus was on\nevaluating thermal conductivity using the Green-Kubo (GK) approach that inherently\naccounts for all contributions to heat transfer, including phonon transport. Therefore,\nwe believe that the conclusions based on the GK results sufficiently capture the net\nimpact of functional groups on thermal transport. We acknowledge that GK does not\nresolve phonon mechanisms; this would require mode-resolved analysis which maybe part\nof future work.\nReviewer Point P2.4 — As provided in Figure S3, most ZIFs with the SOD and\nBCT net topologies exhibited relatively high bulk modulus (7-10 GPa), which is fairly\nclose to the range of ZIFs(mIM) and ZIFs(nIM) (6-10 Gpa). So, why the BCT and SOD\ntopology were except to obtain the Pearson coefficient. How about the overall Pearson\ncorrelation with such two topologies?\nReply: Figure S3 shows that ZIFs with SOD and BCT net topologies exhibit relatively\nhigher bulk moduli compared to other topologies with the same linker. For example, the\nred symbols representing SOD and BCT net topologies have a higher bulk modulus than\nthe other red circles. As a result of their enhanced mechanical stability, these topologies\nalso display higher MD thermal conductivity, as shown in Fig. 3.\nWe supplied the Pearson correlation coefficients calculated both including and exclud-\ning these two topologies in the SI:\nIn Figs. 3, S7a and S7b, we present the Pearson correlation coefficients calculated\nafter excluding the two exceptional net topologies, BCT and SOD,\nwhich exhibit atypical thermal transport behavior. For comparison, Table S3\nsummarizes the Pearson correlation coefficients computed both including and\nexcluding these two net topologies. Our analysis indicates that, given the\nlarge number of topologies studied, the fitted scaling factors remain rela-\n9\ntively robust to these outliers, whereas the Pearson correlation coefficients\nare more sensitive.\nLinker Including Excluding\nIM 0.67 0.69\ndcIM 0.83 0.84\nmIM 0.79 0.76\nnIM 0.78 0.88\nFitted τ 0.86 0.92\nTable 1: Comparison of Pearson correlation coefficients between thermal circuit predic-\ntions and MD thermal conductivity for different linkers, computed including and exclud-\ning the two exceptional net topologies, BCT and SOD.\nReviewer Point P2.5 — According to the definition, each chemical bond is modeled\nas a resistor with a 1Ω resistance. However, the influence of various chemical bonds (i.e.,\nZn-N, C-C) on thermal transport should be quite different. The thermal resistance of\nvarious chemical bonds can be explored to better illustrate the correlation. For example,\nthe values for such 4 functional groups (mIM, nIM, IM, and dcIM) can be determined.\nReply: We consider various chemical bonds as a group. Please see the details in our\nreply to point P1.2 of Reviewer 1.\nComments to the Author: Reviewer #3 has not provided a full review but they have\ninformed us that they have concerns about the novelty of your work as there are several\nrelated works on covalent-organic frameworks (COFs).\nReviewer Point P3.1 — Inparticulartheauthorsraiseimportantconcernsabouthow\ngeneral the model is, and Reviewer #3 indicates that similar work has been reported for\nCOFs. Indeed, itseemsthatfurtherwork, andamorecomprehensiveliteraturerevieware\nnecessary to establish the significance of your results and the strength of conclusions you\ndraw. Moreover, it is far from clear whether the work will continue to appear impressive\nin the light of this additional work. And although we would not rule out consideration\nof a revised manuscript, we feel that the present work is yet at too preliminary a stage\nto support publication in Nature Communications.\nReply: Reviewer #3 is correct that there is a large body of literature reporting experi-\nmental and computational studies of covalent-organic frameworks.\nHowever, our study is a subclass of MOFs, specifically ZIFs, which are unique in that\nthey mimic zeolite pore topologies. However, most importantly, this is a class of materials\nin which, similar to zeolites, one can experimentally change the network topology using\nthesamemetalnodeandlinker. Inaddition, onecanstudytheeffectoffunctionalgroups.\nHence, this provides a chemical design space with 49 distinct pore topologies and four\ndifferent functional groups. Such a subset of materials does not exist in COFs.\nHence, it is simply not possible to perform a similar study with COFs; therefore, we did\nnot include references to the COF literature.\nMost likely, Reviewer #3 missed this point that ZIFs do not have an analogy in COFs,\nbut without a complete report, this is just speculation on our side.\nIn addition, the novelty of our work, in the context of what has been done in MOFs\nand COFs:\n• Circuit model In the MOF literature, one can find screening studies that conclude\n11\nthat thermal conductivity correlates with density. In this work, we demonstrate\nthat a circuit model is a more natural and accurate approach to correlating ther-\nmal conductivity. Such a circuit model would also apply to COFs; however, upon\nreviewing the literature, we could not find any references to this approach. In fact,\nscreening studies that are reported also use density as the main descriptor.\n• Effect of pore topology in MOFs, for a given metal node and linker, one can\nhave many different pore topologies. Currently, we only have insights into how to\nsynthesize these different topologies for ZIFs, which are a special class of MOFs. In\nCOFs, one could envision similar effects, but to the best of our knowledge, this is\neven less studied.\n• Effect of functional groups This is something that can be studied systematically\nin COFs. However, we are not aware of any systematic study directly relevant to\nour research.\nOur conclusion, after a careful analysis of the literature, is that the studies carried out\nwere very similar to those in MOFs (except that many COFs are two-dimensional). Still,\nto the best of our knowledge, none of these studies were directly relevant to our work.\nBecause of these reasons, we did not refer to COFs in our article. In any case, we\nhave added some general references in the introduction to the COF literature on thermal\nconductivity, to avoid the suggestion that we are not aware of the COF literature, which\nmay have triggered the reaction of Reviewer #3.\n12"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The reply and revision have addressed all of my concerns. The current manuscript is acceptable.", "reviewer_2": "We do think insight understanding the thermal conductivity of MOFs, COFs is quite important. Unlike the COFs, the various metal oxide nodes are quite important for determining the thermal conductivity and pose a challenge for the prediction of MOFs. After the revisions made by the authors, all my concerns have been addressed, and I can recommend this manuscript for publication in Nature Communications.", "reviewer_3": "The authors have provided satisfactory responses to the reviewer comments. However, the generality of the proposed modeling approach remains a significant concern. For instance, the study only examines zeolites with thermal conductivities below 1 W m⁻¹ K⁻¹. It is unclear whether the model can be extended to zeolites or MOFs that exhibit higher thermal conductivities. For example, Islamov et al. (npj Computational Materials (2023) 9:11) have reported thermal conductivities exceeding 1 W m⁻¹ K⁻¹ for several MOFs. Can the authors’ model accurately predict such values? If not, the applicability of the model appears to be quite limited."}, "author_response": "Reviewer 1\nRemarks to the Authors:\nThe reply and revision have addressed all of my concerns. The current manuscript is\nacceptable.\nReply: Thank you for the comments and suggestions, which helped us improve the\nmanuscript.\nReviewer 2\nRemarks to the Authors:\nWe do think insight understanding the thermal conductivity of MOFs, COFs is quite\nimportant. Unlike the COFs, the various metal oxide nodes are quite important for\ndetermining the thermal conductivity and pose a challenge for the prediction of MOFs.\nAfter the revisions made by the authors, all my concerns have been addressed, and I can\nrecommend this manuscript for publication in Nature Communications.\nReply: We sincerely thank the reviewer for their positive evaluation and supportive\nrecommendation. We are pleased that the revisions have addressed all concerns, that we\nmanged to clarify the role of metal oxide nodes and the challenges in predicting the ZIFs\nthermal conductivity.\nReviewer 3\nRemarks to the Authors: The authors have provided satisfactory responses to the re-\nviewer comments. However, the generality of the proposed modeling approach remains a\nsignificant concern. For instance, the study only examines zeolites with thermal conduc-\ntivities below 1 W m−1 K−1. It is unclear whether the model can be extended to zeolites\n1\nor MOFs that exhibit higher thermal conductivities. For example, Islamov et al. (npj\nComputational Materials (2023) 9:11) have reported thermal conductivities exceeding 1\nW m−1 K−1 for several MOFs. Can the authors’ model accurately predict such values?\nIf not, the applicability of the model appears to be quite limited.\nReply: We thank the reviewer for this comment. We would like to clarify that zeolites\nand ZIFs, although structurally related, are distinct classes of materials. Although zeo-\nlites and ZIFs share structural similarities, they differ significantly in bonding character\nand vibrational properties. ZIFs are hybrid organic–inorganic materials connected by\nmetal–ligand coordination bonds, whereas zeolites are purely inorganic and covalently\nbonded. Therefore, the adaptation of the theoretical framework should take into account\nthose important differences. The present study specifically focuses on ZIFs and zeolites\nwere not within the scope of this work. We agree that applying the proposed model to\nzeolites, other MOFs, or porous materials with higher thermal conductivities would be\nan interesting direction for future research.\n2"}]

Decision letter