Prediction of thermal conductivity in CALF-20 with first-principles accuracy via machine learning interatomic potentials

CALF-20 (Calgary Framework 20) belongs to the Metal-Organic Framework (MOF) family of materials, known for their high-capacity CO₂ capture and diverse mechanical and thermal properties¹. Similar to other MOFs, CALF-20 features a porous crystalline structure in which Zn atoms and organic ligands create an ordered, periodic network. This arrangement gives rise to intrinsic and periodic voids or pores, leading to a high surface area that is essential for adsorption and catalytic processes. The soft porous crystal structure of CALF-20 makes it highly effective for the physisorption of CO₂ from post-combustion gases, demonstrating excellent CO₂/N₂ selectivity², and maintaining CO₂ capture performance even at high relative humidity (RH)³. The recently achieved scalable synthesis of CALF-20 has captured the interest of researchers, prompting studies into its physical properties and potential successful integration into industrial applications for gas adsorption and separation processes. Alongside the highly selective gas adsorption properties⁴, hydrophobic nature, and outstanding CO₂ capture abilities of CALF-20^5,6 and its derivatives⁷, recent theoretical investigations⁸ have uncovered unconventional thermal and mechanical characteristics of such MOF, exemplified by phenomena like negative area compression (NAC) and negative thermal expansion (NTC). The structure of CALF-20 [Zn₂(1,2,4-triazolate)₂(oxalate)] consists of layers of zinc(II) ions bridged by 1,2,4-triazolate and pillared by oxalate ions, forming a three-dimensional (3D) crystalline lattice (Fig. 1a).

a Primitive unit cell of CALF-20 3D crystal viewed from different angles, highlighting the triazole linker and oxalate pillar of Zn atoms. b HFAC of CALF-20 at 300 K calculated from the raw (gray line) and renormalized (red line) z-component of heat flux using the PW-aiMD trained DeePMD potential. Note the significant noise and oscillatory behavior in the HFAC derived from raw heat flux. The inset shows that the HFAC with renormalized heat flux decays within a few ps. c z-component of the thermal conductivity tensor of CALF-20 calculated from the GK integral of raw (gray line) and renormalized (red line) heat flux at 300 K. Removing non-contributing terms from the heat flux through the renormalization process results in a much smoother κ curve over integration time. Inset: The converged κ of CALF-20, obtained by averaging the cumulative κ-curves over the GK integration time from multiple independent MD runs.

Full size image

According to recent reviews^{9,10,11,12,13}, the study of the thermal transport properties of MOFs has surged exponentially in recent years, driven by their promising applications in greenhouse gas capture^14,15,16 and various technological fields^17,18. Despite being composed of light atoms, MOFs typically exhibit poor thermal conductivity, generally below 1 Wm⁻¹K⁻¹ at normal room temperature, due to their highly porous crystal structure¹⁹. Some of the most studied MOF materials for thermal transport properties include MOF-5 (also known as IRMOF-1)²⁰, HKUST-1²¹, UiO-66²², and various types of zeolitic imidazolate frameworks (ZIFs)^23,24,25,26. The thermal conductivity of these conventional MOFs has been measured experimentally (in both single crystals and powder form) and using molecular dynamics (MD) simulations. The effect of temperature, pressure, various gas adsorbates^{21,27,28,29,30,31,32}, tensile strain^33,34, defects³⁵, linker defects^35,36, interpenetration³⁷, and inorganic material composition³⁸ on the thermal transport properties of MOFs has been extensively studied and reported.

The theoretical approach commonly employed for studying thermal transport in MOFs relies on the Green-Kubo (GK) formalism^39,40 within linear response theory, applied through equilibrium molecular dynamics simulations. Alternative approaches for numerical study of thermal conductivity in crystalline materials include the Peierls-Boltzmann transport equation (PBTE) formalism⁴¹, the direct heating-cooling method in nonequilibrium MD (NEMD) simulations^42,43, and the homogeneous nonequilibrium MD (HNEMD) method^44,45. While the PBTE-based formalism, which treats higher-order anharmonicity perturbatively to calculate phonon–phonon scattering rates, provides deeper physical insight into the phonon thermal transport mechanism, this perturbative approach fails when there is significant anharmonicity in the interatomic potential. The numerical cost of calculating higher-order phonon scattering is prohibitively high for complex materials like MOFs, which have a large number of basis atoms per unit cell. GK formalism, employing a non-perturbative approach, encompasses all the anharmonicity present in the underlying potential of MD simulation and avoids the limitations of the PBTE formalism, particularly when studying materials with very low thermal conductivity. In practice, GK formalism presents certain challenges, including its noisy nature in predicting thermal transport coefficient due to the thermodynamic fluctuations and finite system-size effects, which we will systematically address in the current study. An accurate interatomic potential for MD simulations is essential for the precise and reliable prediction of intrinsic properties such as thermal conductivity in materials. In this context, calculating ab initio heat flux in GK formalism using ab initio molecular dynamics (aiMD) trajectories is paramount. This first-principles approach to calculating thermal conductivity using the GK formalism, known as aiGK, has been successfully implemented and applied to various materials^{46,47,48,49,50}. However, the computational cost of aiMD is significantly higher and can be prohibitive for complex materials like MOFs. Thus, one typically relies on a classical force field (FF) for MD simulations of MOF, such as MOF-FF⁵¹ or the more general universal force field (UFF)⁵². Although MOF-FF has been parameterized and utilized for studying conventional MOFs, it has yet to be adapted for newer materials like CALF-20. On the other hand, the generic UFF presents challenges in terms of accuracy and reliability when calculating physical properties. Machine learning potential (MLP), a very recent development in the field of molecular simulation, offers a promising solution to overcome this challenge. The per-step computational cost of MLP-based MD is several orders of magnitude lower than that of aiMD, while still preserving accuracy and incorporating all the many-body and long-range interactions present in aiMD. The most common ML-based potential that has been successfully implemented and widely adopted for a variety of materials is the several flavors of artificial neural network potential (ANNP)^{53,54,55,56,57,58}. The effectiveness of the ANNP has been demonstrated for various MOFs in the accurate prediction of several physical properties^{59,60,61,62,63,64,65}, including thermal conductivity^66,67. Another class of machine-learned potentials includes Gaussian approximation potentials (GAPs)⁶⁸ and moment tensor potentials (MTPs)⁶⁹, which have also been successfully implemented for MOFs⁷⁰.

In this study, we investigated the thermal transport properties of CALF-20 through the GK formalism based on MD simulations utilizing MLP. We employed the smooth edition of Deep Potential (DeepPot-SE)⁷¹ from the DeePMD-kit^72,73 to construct an ANNP for CALF-20. In this approach, a neural network potential energy surface (NN-PES) for the system is generated from a training dataset that includes atomic positions, forces, potential energy, and virial obtained through aiMD. The training dataset must encompass a sufficient aiMD trajectory, ensuring a substantial number of points in the configuration space are used to train the neural network. This thorough training ensures the correct physical description of the systems and prevents the presence of any spurious extrema on the NN-PES. From this point forward, we refer to the MLP generated using this approach as the DeePMD potential.

We developed two types of DeePMD potentials, each corresponding to aiMD training data generated using different levels of ab initio theory. The first is derived from aiMD with a plane wave (PW) basis set, referred to as the PW-aiMD trained DeePMD potential. The second is based on aiMD with a linear combination of atomic orbitals (LCAO) basis set, which we denote as the LCAO-aiMD trained DeePMD potential. The LAMMPS⁷⁴ package is used for all the MD simulations. Note that LAMMPS is compatible with the DeePMD potential, with the interface provided by the DeePMD-kit, which also includes a GPU-accelerated version⁷⁵.

The rest of the paper is organized as follows: In the next section, we discuss the method of heat flux renormalization and its significance for the precise and reliable prediction of the thermal conductivity (κ) of CALF-20 with DeePMD potential. Next, we present the calculated κ values of CALF-20 using both the PW-aiMD trained and LCAO-aiMD trained DeePMD potential. We also present a comparative analysis of the κ of CALF-20 calculated using the UFF-based MD and the DeePMD approach. Finally, we examine the effects of temperature and pressure on the κ of CALF-20, summarize the key findings, discuss the limitations of the current study, and offer perspectives for future research.

Heat flux renormalization

The microscopic definition of the system’s heat flux can be expressed as a combination of convective and virial contributions, as follows⁷⁶:

Here, Ω is the crystal volume and ({e}_{i}=frac{1}{2}{m}_{i}{leftvert {{{{boldsymbol{v}}}}}_{i}rightvert }^{2}+{u}_{i}) represents the total energy of the atom i, with m_i being atomic mass, v_i being velocity, and u_i is the potential energy of the atom i. The term of the above expression for J(t) that includes the atomic stress tensor ({bar{S}}_{i}={sum}_{j}({{{{boldsymbol{r}}}}}_{i}-{{{{boldsymbol{r}}}}}_{j})otimes frac{partial {u}_{j}}{partial {{{{boldsymbol{r}}}}}_{i}}), where the index j run over all the atoms in the system, constitute the virial contribution to the heat flux. Here, ⊗ denotes the tensor product of two vectors, and ({bar{S}}_{i}) is a 3 × 3 matrix representing the stress tensor associated with the ith atom. Note that the centroid method^77,78 is employed to account for the effects of many-body interactions in the per-atom stress tensor. The thermal conductivity values presented in this paper include this multi-body correction (refer to Supplementary Note 11 for more details on the centroid method). It is important to note that the convective component of the heat flux J(t), which includes the energy term e_i, does not contribute to thermal transport in crystalline solids⁴⁹, such as CALF-20 (see Supplementary Note 5). The expression of thermal conductivity tensor (({kappa }_{alpha beta }^{Xi (0)})) for an initial phase-space configuration Ξ(0) of the system is given by the GK integral of heat flux:

where α and β represent the cartesian direction (x, y, z). Here, 〈. . . , . . . 〉 denotes the autocorrelation of heat flux (HFAC) over the simulation time. Note that here (Xi equiv {{{{{{boldsymbol{R}}}}}_{{{{boldsymbol{i}}}}};{{{{boldsymbol{P}}}}}_{{{{boldsymbol{i}}}}}}}_{i = 1,ldots ,N}) represents phase-space configuration of the system, where R_i is the position vector and P_i denotes the momentum vector of atom i. J^Ξ(t) is the heat flux corresponding to the phase-space point Ξ(t), which is obtained from the evolution of the system in time through the Hamiltonian ({{{mathcal{H}}}}(Xi )) on initial configuration Ξ(0). The thermal conductivity (κ_αβ) is defined as the phase-space integral of ({kappa }_{alpha beta }^{Xi }) as: ({kappa }_{alpha beta }=frac{1}{{{{mathcal{Z}}}}}{int}_{Xi }{kappa }_{alpha beta }^{Xi }{e}^{-frac{{{{mathcal{H}}}}(Xi )}{{k}_{B}T}}dXi), where ({{{mathcal{Z}}}}) represents the canonical partition function. This expression of κ_αβ can be approximated to apply through the MD simulation as the average of ({{{mathcal{N}}}}) independent MD runs with different starting phase-space configurations (Ξ^j):

It can be demonstrated that the above expression of the heat flux J(t) (Eq. (1)) encompasses non-contributing terms in the GK integral (Eq. (2)) of the thermal conductivity tensor. In practice, these terms in the heat flux from MD runs on a finite-size system generate noise into the κ value. Figure 1b illustrates the HFAC for the cartesian z-component of heat flux of CALF-20 at a temperature of 300 K. Notice the pronounced oscillations and noise in the HFAC (gray curve) with correlation time, which arises from the raw heat flux obtained during the MD simulation using the DeePMD potential. HFAC derived from the raw heat flux does not exhibit any decay trend even after a high correlation time of 100 ps. This presents a problem for the convergence of the thermal conductivity (κ) value, which depends on the time integration of the HFAC. Consequently, the κ value computed from the raw heat flux data using the DeePMD potential exhibits significant oscillations over the integration time, making accurate and reliable prediction of the κ value unfeasible using this approach. This is evident in Fig. 1c, which shows the z-component of κ from the raw HFAC for CALF-20 at a temperature of 300 K (gray curve). The noisy nature of the κ calculated from the GK integral (Eq. (2)) of raw heat flux from ab initio MD has also been previously observed^49,79. As described in refs. ^80,81, the time-dependent per-atom potential energy (u_i) and per-atom stress tensor (({bar{S}}_{i})) can be expanded in a Taylor series around the initial ({u}_{i}^{0}) and ({bar{S}}_{i}^{0}) of the relaxed structure. This initial constant per-atom potential energy (({u}_{i}^{0})) and stress tensor (({bar{S}}_{i}^{0})) contribute as a linear term of the atoms’ velocity (v_i) in the heat flux of the system. The time integral of this linear velocity term, which appears in the GK integral of the thermal conductivity tensor, becomes zero. The heat flux calculated from the MD trajectories implicitly includes this zeroth-order term. In certain crystalline solids with basis atoms consisting of only a single element and possessing certain crystalline symmetry, this linear velocity term effectively vanishes in the GK integral⁸². However, this is not the case for the material studied in this work. The non-contributing part of heat flux in the thermal transport is expressed as: ({{{{boldsymbol{J}}}}}_{0}=frac{1}{Omega }{sum}_{i}[{u}_{i}^{0}{{{{boldsymbol{v}}}}}_{i}-{bar{S}}_{i}^{0}cdot {{{{boldsymbol{v}}}}}_{i}]). Here, we define the heat flux excluding this inert component toward κ as the renormalized heat flux: J_ren.(t) = J(t) − J₀. Subtracting these non-contributing terms from the heat flux at each timestep would significantly reduce numerical noise in the HFAC. This is evident from Fig. 1b, which compares the HFAC calculated from the renormalized heat flux (red curve) of CALF-20 at 300 K to the raw HFAC (gray curve). The inset plot of Fig. 1b highlights the significant reduction in the amplitude of the HFAC oscillations after renormalizing the raw heat flux. The renormalized HFAC decays to zero after a few picoseconds, resulting in the thermal conductivity becoming a smooth function of the integration time in the GK integral. Figure 1c shows the κ from the GK integral of the renormalized heat flux for CALF-20 (red curve). Compared to the κ calculated from the raw heat flux (gray curve), which is exceedingly oscillatory and even has negative values for the present material with low thermal conductivity, the renormalized heat flux produces a much smoother κ curve over the integration time. Thus, the renormalization process, which eliminates non-contributing terms in the heat flux, results in a more reliable prediction of the thermal transport coefficient using the GK formalism. For completeness, here we discuss another technique to discard the non-contributing part of the heat flux towards the thermal transport coefficient, known as gauge invariance principles, as described in refs. ^79,83. According to the gauge invariance principle, incorporating a linear combination of nondiffusive fluxes—those whose GK integral vanishes—into the actual heat flux does not affect the final thermal conductivity value. It turns out that parameter minimization on the HFAC via a variation principle yields the nondiffusive flux, as described by the expression: ({{{{boldsymbol{J}}}}}_{gauge}=frac{1}{Omega }{sum}_{i}[{langle {e}_{i}rangle }_{t}{{{{boldsymbol{v}}}}}_{i}-{langle {bar{S}}_{i}rangle }_{t}cdot {{{{boldsymbol{v}}}}}_{i}]), and the normalized heat flux: J_ren.(t) = J(t) − J_gauge. Here, 〈 〉_t denotes the time-averaged quantity, and the convective part can be simplified as ({langle {e}_{i}rangle }_{t}=frac{5}{2}{k}_{B}T+{langle {u}_{i}rangle }_{t}). A similar level of noise reduction in the value of the thermal transport coefficient is achieved by discarding the non-contributing part of the raw heat flux using the gauge theorem. This “gauge filtering” on heat flux has been adopted and successfully applied to the thermal transport study of crystalline solids, as described in the refs. ^49,80. However, we observed that the numerical value of κ differs negligibly between the renormalization of raw heat flux, as detailed in the first case, and using the gauge theorem. This marginal deviation cancels out when averaging several independent MD runs, as done in this study (see Supplementary Note 4). Therefore, we will adhere to the first kind of renormalization of the raw heat flux from the MD run using the DeePMD potential on CALF-20. Excluding the passive terms from the raw heat flux results in a less noisy κ curve, but it still does not guarantee the convergence of κ towards a definitive value. This is apparent from Fig. 1c, which plots κ against integration time for CALF-20 at 300 K (red curve). The non-converging nature of thermal conductivity within the GK framework poses challenges in determining the κ value of the material with certainty. Even averaging the κ over a specific range of integration times from a single MD run is unreliable, as it results in a significant error bar, and the values may vary substantially when another independent MD run is conducted. We note that there are approaches to numerically evaluating the thermal transport coefficient from the cumulative κ curve of the HFAC integral while relying on a minimal number of tunable parameters. The notable one is the “first dip (FD) method”, described in ref. ⁸⁴. In the FD criterion, the κ(t_c) is chosen as the thermal conductivity of the material, where t_c corresponds to the correlation time at which HFAC falls to zero. Beyond this time point, the accumulated numerical noise leads to a drop in the κ curve. While this approach is effective for simple crystalline solids⁴⁹, it is not suitable for the current case of CALF-20, a complex polyatomic crystal, as evidenced by the κ curve in Figure 1c. In this study, we average the κ values from several independent MD runs at a specific temperature to achieve a convergent κ curve over the integration time. This averaging technique is typically used to predict the GK-theory-oriented thermal transport coefficient from MD runs based on classical empirical force fields^85,86, where the lower computational cost allows for multiple independent MD runs. The computational cost of the MD based on the MLP falls between that of aiMD and classical MD, enabling us to use the averaging technique while maintaining ab initio level accuracy. Furthermore, thermal conductivity in GK formalism, by definition, is an averaged quantity over the initial phase-space configurations, as previously mentioned (Eq. (3)). In the inset plot in Fig. 1c, we demonstrate that the κ curve obtained through the averaging of multiple independent MD runs (green curve) of CALF-20 at a temperature of 300 K converges toward a specific value with the integration time. The final thermal conductivity value is determined from the mean of κ over the last 30 ps of the cumulative κ in the GK integral, with the error bar representing the standard deviation within this time window.

Comparison of DeePMD and UFF

With the robust framework for the accurate prediction of the thermal conductivity for CALF-20 using the GK formalism applied through equilibrium MD based on DeePMD-potential, as outlined in the previous section, we will now discuss the results and compare them with various levels of theory in both ab initio and classical FF. We first discuss the κ value obtained from the DeePMD potential, which was trained using the trajectories of ab initio MD based on the plane-wave basis set used for the expansion of electronic wave functions (PW-aiMD trained). Refer to the “Methods” section for more detailed information on the ab initio calculations. The values of the thermal conductivity tensor of CALF-20 at 300 K temperature and atmospheric pressure, obtained using the PW-aiMD trained DeePMD potential, are κ_xx = 0.85 ± 0.01 Wm⁻¹K⁻¹, κ_yy = 0.69 ± 0.01 Wm⁻¹K⁻¹, and κ_zz = 0.57 ± 0.02 Wm⁻¹K⁻¹. The off-diagonal components of the κ tensor (κ_xy, κ_xz, and κ_yz) are zero, as shown in the Supplementary Note 6. Notice the anisotropic behavior of the thermal transport coefficient in such crystalline material. Direction-dependent lattice thermal conductivity is expected when cubic symmetry is absent in the crystalline solid, as is the case with CALF-20. Similar anisotropic thermal transport characteristics have also been observed in Zn-MOF-74⁸⁷. In addition, note the ultra-low thermal conductivity value from our calculations. This extremely low thermal conductivity is advantageous for thermoelectric energy conversion devices, as the figure of merit (ZT) of a material for such application is inversely related to κ^{88,89,90,91,92,93}. To our knowledge, the experimentally measured thermal conductivity of CALF-20 has not been reported, and no theoretical predictions exist either. As such, we compare our calculated results with the existing κ value of well-studied MOFs. The experimentally reported κ values for some common MOFs at room temperature are 0.69 Wm⁻¹K⁻¹ for HKUST-1 (single crystal)²¹, 0.32 Wm⁻¹K⁻¹ for IRMOF-1 (single crystal)²⁰, and 0.32 Wm⁻¹K⁻¹ for ZIF-8 (thin films)²³. Our prediction of the thermal conductivity of CALF-20, which is a Zn-based MOF, using the GK formalism through MLP-based MD, aligns with the trend of low κ values observed in other similar MOFs. The very low κ value of CALF-20, like other MOFs, is attributed to their highly porous structure and low density. For CALF-20, the pore volume is 0.32 cm³ gm⁻¹, the void fraction is 0.51, the largest cavity diameter (LCD) is 4.2 Å, and the pore limiting diameter (PLD) is 2.8 Å⁷. The correlation between the thermal transport coefficient calculated in this study and the quantitative values of the pore structure of CALF-20 aligns with the findings of the ref. ⁹⁴ that extensively studied the relationship between κ and the pore structure of various MOFs. To further validate our calculated κ value for CALF-20 using PW-aiMD trained DeePMD potential and to study the impact of the choice of a different basis set in first-principles calculations for generating DeePMD training data, we also performed aiMD with the linear combination of atomic orbital (LCAO) basis set. We acronym the DeePMD potential of CALF-20 generated from this level of ab initio theory as LCAO-aiMD trained. See the “Methods” section for further details on the ab initio calculations in this case. The GK theory predicted κ values from the LCAO-aiMD trained DeePMD potential of CALF-20 at 300 K and atmospheric pressure are κ_xx = 1.23 ± 0.03 Wm⁻¹K⁻¹, κ_yy = 0.68 ± 0.01 Wm⁻¹K⁻¹, and κ_zz = 0.57 ± 0.01 Wm⁻¹K⁻¹. Figure 2a presents a comparison of the three diagonal components of the κ tensor calculated using the PW-aiMD train and LCAO-aiMD trained DeePMD potentials for CALF-20. The HFAC and the convergence of κ with integration time with the LCAO-aiMD trained DeePMD potential are shown in Supplementary Note 7. It is worth noting that the x-direction κ of CALF-20 is ~45% higher with the LCAO-aiMD trained DeePMD potential compared to PW-aiMD trained potential, while the other two components are roughly the same within the error bars for both cases. The 45% higher prediction of κ_xx is consistently reproduced using an independently generated MLP trained on a separate dataset comprising LCAO-aiMD trajectories (refer to Supplementary Note 3). To ensure that the DeePMD potential accurately reflects the characteristics of the underlying ab initio theory of the aiMD training data, we plotted the forces on each atom of the CALF-20 supercell as calculated from the ab initio level and the corresponding DeePMD potential. These values align perfectly along the agreement line (y = x) with a root mean square error (RMSE) of 0.025 eVÅ⁻¹ for the LCAO-aiMD trained DeePMD potential, as shown in Supplementary Note 1. In addition, the RMSE of energy and force for both the training and validation datasets, which are provided in Supplementary Note 1, confirms the excellent training quality. Therefore, the underlying theory and approximations of the ab initio method used have an impact on the prediction of the thermal transport coefficient of materials. In the case of CALF-20, our investigation revealed a slightly higher κ value prediction with the LCAO-aiMD trained DeePMD potential, while the computational cost of aiMD based on the LCAO basis set (~O(N)) is significantly lower than that of the PW basis set. We additionally investigate the thermal transport properties of CALF-20 based on the universal force field (UFF)⁵², which is a classical empirical force field. The parameters of bonded and non-bonded terms in the UFF of CALF-20 include some modifications, which we adapted from the ref. ⁹⁵. The thermal conductivity of CALF-20 at 300 K and atmospheric pressure, as predicted by UFF, is κ_xx = 10.39 ± 0.3 Wm⁻¹K⁻¹, κ_yy = 6.20 ± 0.2 Wm⁻¹K⁻¹, and κ_zz = 4.14 ± 0.1 Wm⁻¹K⁻¹, as shown in Fig. 2b. Figure 2c demonstrates the convergence of κ_zz from classical UFF-based MD with integration time (see Supplementary Note 9 for HFAC and convergence of κ_xx and κ_yy in case of UFF). It is worth noting that our calculation includes the modification in the per-atom stress tensor due to the many-body interaction of classical FF, as detailed in refs. ^77,96 (see Supplementary Note 12 for κ values without this many-body correction). The built-in heat flux calculator of the LAMMPS code was utilized for the UFF, and interestingly, the renormalization of heat flux from UFF-based calculations led to only a slight reduction in noise, as the HFAC and thermal conductivity curves were already sufficiently smooth (see Supplementary Note 10). Note that the κ value predicted by UFF for CALF-20 is roughly an order of magnitude higher than that obtained from DeePMD potential-based MD. Such a high thermal transport coefficient is uncommon in most of the stable and experimentally feasible MOFs. The overestimation of CALF-20’s elastic moduli with UFF compared to the values predicted by DeePMD has also been reported in a previous study on mechanical properties⁸. Although the GK theory, which leverages the fluctuation-dissipation principle, implicitly accounts for all the anharmonicity in the materials, the UFF is inherently limited to four-body (dihedral and improper) interactions. In contrast, the DeePMD potential, based on artificial neural networks, encompasses all the many-body interactions present in the aiMD training data. Consequently, higher-order phonon scattering, mediated by higher-order anharmonicity, is present with the DeePMD potential, which tends to reduce the thermal transport coefficient. In addition, generic empirical force fields like UFF are not always suitable for predicting intrinsic properties of materials, such as thermal conductivity, unless the parameters are optimized specifically for that purpose. This leads to an overprediction of the κ value of CALF-20 in the case of UFF compared to DeePMD. Although the absolute κ values differ between UFF and DeePMD, both show a similar direction-dependent nature of κ-tensor (κ_xx > κ_yy > κ_zz) for CALF-20. This anisotropic thermal conductivity is beneficial for designing devices that require an easy axis for heat dissipation while providing heat insulation in other directions.

a Bar plot comparing the three diagonal components of the κ-tensor of CALF-20 at 300 K, calculated using the PW-aiMD trained and LCAO-aiMD trained DeePMD potentials. Note the anisotropic nature of κ. The κ_xx predicted by the LCAO-aiMD trained DeePMD potential is ~45% higher than that of the PW-aiMD trained DeePMD, while the other two components (κ_yy and κ_zz) remain similar between the two types of DeePMD potentials. b Three diagonal components of κ of CALF-20 at 300 K calculated using the UFF from classical MD. UFF predicted κ values are approximately an order of magnitude higher compared to DeePMD, yet both exhibit a similar directional dependence in the κ-tensor (κ_xx > κ_yy > κ_zz). c Convergence of the z-direction thermal conductivity (κ_zz) calculated using the UFF within classical MD.

Full size image

Temperature-dependent thermal conductivity

Figure 3a shows the temperature-dependent thermal conductivity of CALF-20 over a range of 200 K to 500 K, in 25 K increments. Since MOF materials are primarily used in industrial-scale post-combustion gas adsorption processes, they often operate at temperatures other than room temperature. Thus, it is crucial to study how κ behaves with temperature variation. Our results indicate an overall decrease in the κ value of CALF-20 as the temperature (T) rises, though this trend is not uniformly monotonic across the entire temperature range. It is important to note that, from this point forward, we have consistently used the PW-aiMD trained DeePMD potential for CALF-20. The anisotropic, direction-dependent characteristics of the κ-tensor (κ_xx > κ_yy > κ_zz) are consistently maintained throughout the temperature range studied. As the temperature increases from 200 K to 500 K, the z-directional κ of CALF-20 decreases from 0.66 Wm⁻¹K⁻¹ to 0.35 Wm⁻¹K⁻¹, representing a reduction of ~47%. The average thermal conductivity (κ_ave.) exhibits a much smoother decreasing trend compared to the individual components, as illustrated by the gray curve in Fig. 3a. The κ_ave. is defined as: ({kappa }_{ave.}=frac{1}{3}{sum}_{alpha =x,y,z}{kappa }_{alpha alpha }). The κ_ave. of CALF-20 decreases from 0.87 Wm⁻¹K⁻¹ to 0.53 Wm⁻¹K⁻¹, with a value of 0.71 Wm⁻¹K⁻¹ at 300 K. To quantify the temperature-dependent effect on the thermal transport of CALF-20, we fitted κ_ave. to a power law of temperature (κ ∝ 1/T ^α). We obtained the exponent α as 0.56 for CALF-20, as shown in the inset plot of Fig. 3a. For typical crystalline semiconductors and insulators, α is close to 1, as Umklapp phonon-scattering dominates at higher temperatures and limits κ⁴¹. However, this is not the case for CALF-20, and there is a substantial deviation from the ideal case of  ~1/T. Our analysis of the electronic band structure for CALF-20 reveals a band gap of 3.4 eV (see Supplementary Note 14), identifying it as an electrical insulator. Consequently, the contribution of free electrons to thermal conduction remains negligible, even at higher temperatures. The weak temperature dependence of the thermal transport coefficient in several MOFs has been previously studied and reported both experimentally and through theoretical calculations. For instance, HKUST-1⁹⁷, ZIF-4^25,26, ZIF-62²⁵, ZIF-8^98,99, and MOF-5^66,100 exhibit minimal to no effect of temperature on the κ value over a similar temperature range as studied here. The exponent α for MOF-5 was reported to be ~0.6, as determined by MLP-based MD simulations detailed in ref. ⁶⁶. The temperature-independent glass-like thermal conductivity of ZIF-4 and ZIF-62, two Zn-based MOFs, was analyzed in ref. ²⁵. Their findings indicate that the large mass difference between Zn atoms and the organic sites results in enhanced phonon scattering, limiting the mean free path (MFP) to below 1 nm and leading to a low κ value. In addition, the contributions of propagating and non-propagating modes to κ exhibit anomalous temperature behavior, resulting in a glass-like temperature-independent κ. However, in the case of CALF-20, which is also a Zn-based MOF, we observe that κ decreases with temperature despite the dependency being relatively weak. To gain further insight into the effect of temperature on the thermal transport properties of CALF-20, we calculated the vibrational density of states (VDOS) from the MD simulation trajectory. Figure 3b presents the VDOS of CALF-20 at four different temperatures. It is important to note that the sharp peaks in the VDOS indicate the crystalline solid characteristics of CALF-20, in contrast to glass-like materials, where the VDOS is much smoother. In addition, as temperature increases, the peaks in the VDOS shift toward lower frequencies (redshift), as highlighted in the inset plot of Fig. 3b. To examine the effect of T on the phonon modes of CALF-20, we calculated the phonon dispersion at all temperatures considered in this study, as shown in Fig. 3c. Phonon dispersion of CALF-20 was computed using the DeePMD potential, incorporating the non-analytic term in the second-order force constant arising due to dipole-dipole interactions. The absence of imaginary phonon frequencies highlights the MLP’s effectiveness in accurately modeling the stable crystal structure of CALF-20. The substantial mass difference among the constituent atoms of CALF-20 results in pronounced frequency gaps between the phonon branches. In Fig. 3c, we focus on the low-lying branches, encompassing the acoustic modes and select higher-frequency optical modes. We observed a softening of phonon modes with increasing temperature, with the effect being more pronounced in the higher-frequency optical branches compared to the low-lying branches, as shown in the figure. This softening aligns with the redshift of VDOS peaks calculated from the velocity autocorrelation function in MD simulations. This softening of phonon modes can explain the κ reduction of CALF-20 with temperature. Thermal conductivity is directly proportional to the phonon lifetime (τ) and phonon group velocity (v ^g) in relaxation time approximation (RTA) of PBTE formalism⁴¹:

Here, α denotes the cartesian direction, ℏ represents the reduced Planck constant, N is the number of unit cells, V is the unit cell volume, λ represents the phonon modes, n_λ corresponds to the Bose-Einstein distribution at temperature T, and ω(λ) is the phonon frequency. Calculating τ numerically requires the evaluation of scattering matrix elements from third- and fourth-order anharmonic force constants. For materials like CALF-20, with a large number of basis atoms per unit cell and low symmetry, these calculations become computationally infeasible. Therefore, we interpret the behavior of κ using the quantitative values of phonon frequency and group velocity. The phonon group velocity of CALF-20 exhibits minimal variation with temperature, as illustrated in Figure 3d. However, the redshift of vibrational modes results in more strengthened phonon-scattering events, where the energy conservation condition among participating phonon modes is satisfied⁴¹. This leads to a shorter phonon lifetime (τ). As such, the decrease in κ with temperature is driven by the softening of phonon frequencies (ω(λ)) and the reduction in phonon lifetimes (τ(λ)), even though κ for CALF-20 shows only a weak dependence on T, approximately scaling as 1/T^0.56.

a Temperature (T) variation of thermal conductivity of CALF-20 from 200 K to 500 K in 25 K increments. The gray diamond symbol represents the average (κ_ave.) of the three non-zero directional κ-tensor components (blue square for the κ_xx, red pentagon for the κ_yy, and green star represents κ_zz). The error bars are defined in the “Heat flux renormalization” subsection of the Results and discussion section. Inset: Fitting of κ_ave. with a power law of T (κ ~ 1/T^α), showing the exponent α as 0.56, indicating weak temperature dependence. b Vibrational density of states (VDOS) of CALF-20 at four different temperatures. Notice the redshift of VDOS peaks as temperature rises. c Temperature-dependent phonon dispersion of CALF-20 along the high-symmetry path of the Brillouin zone, calculated using the PW-aiMD trained DeePMD potential. Frequency gaps between phonon branches, caused by significant mass differences among constituent atoms, are omitted, with a focus on low-frequency branches, including acoustic modes, and select higher-frequency optical branches. Pronounced softening of higher-frequency phonon modes with increasing temperature is evident. d Phonon group velocity of CALF-20 at different temperatures, showing minimal variation with temperature changes.

Full size image

Pressure-dependent thermal conductivity

In addition to temperature, pressure can affect the gas loading capacity of CALF-20, with CO₂/N₂ uptake generally increasing with higher pressure³. Therefore, examining how κ responds to varying pressure is crucial for its efficient application in diverse external settings. Figure 4a depicts the pressure-dependent thermal conductivity of CALF-20 at 300 K under hydrostatic pressure ranging from 0 to 1 GPa. Note that the isotropic pressure is applied during the equilibration phase of MD simulation through the Nosé-Hoover barostat (see the “Methods” section for further details). Our findings indicate that the κ of CALF-20, based on the DeePMD potential, shows minimal change with increasing pressure. The x-directional thermal conductivity (κ_xx) decreases from 0.85 ± 0.01 Wm⁻¹K⁻¹ at 0 GPa to 0.70 ± 0.01 Wm⁻¹K⁻¹ at 1 GPa, representing an approximate 18% decrease. In contrast, κ_zz decreases significantly from 0.57 ± 0.02 Wm⁻¹K⁻¹ to 0.34 ± 0.01 Wm⁻¹K⁻¹, a reduction of 40%. The κ_ave. (gray curve) exhibits an overall decline of about 24% as pressure increases, although this trend is not monotonic throughout the pressure range. For conventional semiconductors and insulators, κ can either decrease or increase with pressure, influenced by the mass ratio of the constituent atoms^41,101, and may exhibit a nonmonotonic trend when higher-order phonon scattering is considered¹⁰². For MOF materials, pressure-dependent thermal transport was studied for ZIF-8 for the same pressure range examined here, using MD simulations based on the MOF-FF⁹⁹. The reported study on ZIF-8 combined with different functional groups indicates a minimal decreasing trend in thermal conductivity up to a certain critical pressure, after which κ abruptly increases. This behavior is attributed to the phase transition of ZIF-8 from a porous to a dense phase as pressure rises. However, for CALF-20, we did not observe any such unconventional pressure dependence of κ. A metastable dense phase of CALF-20 is known to exist beyond 1 GPa⁸, which is outside the pressure range investigated in this study. The decreasing trend of thermal conductivity in crystalline solids with pressure is typically observed in materials with a large mass ratio of constituent atoms. These compounds exhibit a significant gap between the phonon frequencies of acoustic and optic branches. As pressure increases, this gap widens further, leading to enhanced phonon scattering due to more satisfied energy conservation conditions among the participating phonon modes¹⁰¹. This increased scattering among the phonon modes leads to the lowering of lattice thermal conductivity. CALF-20 also exhibits a significant frequency gap between its phonon branches, attributed to the considerable mass disparity between the Zn atoms and the organic sites. Figure 4b presents the phonon dispersion of CALF-20 under varying pressures. As pressure increases, phonon modes exhibit stiffening, with the effect being more prominent in the higher-frequency optical branches than in the low-lying branches. We also examined the phonon group velocity of CALF-20 under varying pressures, as depicted in Fig. 4c. While the absolute values show minimal variation, the group velocity shifts toward higher frequencies with increasing pressure, consistent with the observed stiffening of phonon modes. Consequently, the reduction in phonon lifetime, driven by increased scattering among phonon modes, leads to the decrease in the κ value with increasing pressure in CALF-20, following Eq. (4). Nonetheless, the reduction of κ in response to the pressure is very minimal and highly dependent on the heat transport directions in CALF-20.

a Thermal conductivity of CALF-20 as a function of pressure from 0 to 1 GPa, demonstrating a slight overall reduction in κ as pressure increases. (blue square for the κ_xx, red pentagon for the κ_yy, green star represents κ_zz, and the gray diamond represents κ_ave.). The error bars are defined in the “Heat flux renormalization” subsection of the Results and discussion section. b Phonon dispersion of CALF-20 along the high-symmetry path of the Brillouin zone, calculated at three different pressure conditions at 300 K. The stiffening of phonon modes with increasing pressure is evident, particularly in the higher-frequency optical branches. c Phonon group velocity of CALF-20 at three different pressures (Green crosses for 0 GPa, blue triangles for 0.3 GPa, and red plus signs for 0.7 GPa pressure). Notice that the variation in the absolute value of group velocity is minimal while they shift towards higher frequency due to the stiffening of phonon modes with increasing pressure.

Full size image

To the best of our knowledge, this is the first investigation of thermal transport phenomena in CALF-20. In summary, we systematically applied the Green-Kubo (GK) formalism via equilibrium MD, utilizing artificial neural network-based MLP to calculate thermal conductivity. The renormalization of raw heat flux in the GK framework to exclude non-contributing terms is implemented for the MLP to ensure a reliable and noise-free calculation of thermal conductivity. We developed and compared the calculated κ values with two types of MLP, based on different approximations and basis sets in the ab initio methods. The PW-aiMD trained DeePMD potential yields a κ of CALF-20 below 1 Wm⁻¹K⁻¹ at 300 K. This very low κ value is typical for most MOFs. The MLP for CALF-20 through aiMD based on the LCAO basis set, which is computationally more efficient, yields a slightly higher κ value compared to that predicted from the plane wave basis. In our study, the UFF-based classical MD simulation produces a κ value approximately an order of magnitude higher than that of the DeePMD-based MD. Such a higher κ value is uncommon for most stable MOFs. We examined the impact of temperature on the thermal conductivity of CALF-20 and discovered a weak temperature dependence of κ. Finally, we analyzed the effect of pressure on the thermal transport of CALF-20 and observed a nearly constant trend in κ with a marginal decrease as pressure rises. The primary application of CALF-20 is CO₂ capture at an industrial scale. Therefore, it would be advantageous to examine the impact of various gas adsorbates and structural defects on the thermal transport coefficient of CALF-20 for practical applications. However, this is beyond the scope of the current study. The synthesis of CALF-20, both in single-crystalline form³ and under standard room temperature and pressure conditions¹⁰³, has been recently achieved. Consequently, our study will guide the experimental thermal transport study of CALF-20 and more efficient practical heat management applications based on such material.

Ab initio MD

The aiMD runs for preparing the training dataset for CALF-20 were conducted using the SIESTA code¹⁰⁴, which utilizes the linear combination of atomic orbital (LCAO) basis for ab initio calculations. We expanded the wave functions using a polarized double-zeta basis set. Norm-conserving pseudopotentials¹⁰⁵ were used, with exchange-correlation handled via the GGA-PBEsol¹⁰⁶ approach. The vdW interactions were also incorporated through Grimme’s DFT-D3 dispersion model¹⁰⁷. During the self-consistent field (SCF) cycle, the tolerances for the density matrix and Hamiltonian were set to 10⁻⁹ and 10⁻⁸ eV, respectively. The temperature in the aiMD simulation was regulated using the Nosé-Hoover thermostat, while the pressure was maintained with the Parrinello-Rahman barostat. The integration timestep in the aiMD run was set to 0.5 fs. All aiMD runs were conducted on a geometry-optimized 2 × 2 × 2 supercell of CALF-20. The DeePMD training utilized aiMD trajectories from multiple temperature-pressure configurations of CALF-20, encompassing a total of 46,327 MD snapshots (Table 1).

The DeePMD-potential trained through the plane-wave (PW) basis set-oriented ab initio calculation was adapted from refs. ^8,108. Follow the corresponding reference for further details. Here, we mention that the aiMD runs for PW-trained DeePMD potential were carried out through the projected augmented wave (PAW) method implemented in the VASP code¹⁰⁹. The exchange-correlation functional was treated using the Perdew-Burke-Ernzerhof (PBE)¹¹⁰ generalized gradient approximation, and the van der Waals (vdW) interactions were also included via Grimme’s DFT-D3 dispersion model¹⁰⁷. The training dataset comprises aiMD runs for CALF-20 at various temperatures and pressures.

DeePMD training

The training and generation of the artificial neural network potential (ANNP) were conducted using the DeePMD-kit software package (version: 2.2.9). For more detailed information on the functionality of the DeePMD-kit, readers are advised to follow the refs. ^71,72,73. Here, we outline some key training parameters used for CALF-20. The descriptor used in this work, referred to as “se_e2_a” (also known as DeepPot-SE)⁷¹, is constructed from atomic distance data while encoding both angular and radial information and thus complete multi-body interaction in the given atomic configurations. The cutoff radius and smoothing radius were set to 8 Å and 0.6 Å, respectively. The embedding network size was {25, 50, 100}, and the fitting network size was {240, 240, 240}. 20% of the randomly selected aiMD snapshots were used for the validation dataset, while the remaining 80% were used for the training dataset. A total of one million training steps were performed to achieve convergence of the RMSE in energy and force (see Supplementary Note 1).

MD simulations

MD simulations using the ANNP for CALF-20 were conducted with the LAMMPS⁷⁴ package. It is important to note that the interface for the ANNP in LAMMPS was provided by DeePMD-kit, and thus, the version of LAMMPS (version: 2 Aug 2023) that accompanies DeePMD-kit was used throughout this study. All ANNP-MD simulations were conducted on a 3 × 3 × 3 supercell of CALF-20 unless otherwise specified. The integration timestep of the MD simulation was set to 0.5 fs. The Nosé-Hoover thermostat and barostat were used for the canonical (NVT) and isothermal-isobaric (NPT) ensembles, respectively. The damping parameter of the thermostat was 50 fs, and for the barostat, 500 fs. The system underwent an initial energy minimization using a conjugate-gradient algorithm. This was followed by a 1 ns NPT simulation to equilibrate at the specified pressure and temperature. Subsequently, an NVT simulation was conducted for 500 ps. Finally, the system was simulated under a microcanonical ensemble (NVE) for 2 ns, during which the heat flux was calculated every 1 fs for thermal conductivity determination using the GK theory. The same simulation protocol was followed for the UFF-based calculation with LAMMPS.