Metrics-based computational method selection for the prediction of a physical property

Abstract

Examples are disclosed that relate to the selection of a method to compute a physical property based upon metrics obtained using a universal machine learning force field. One disclosed example provides a computing system comprising a logic subsystem, and a storage subsystem comprising instructions executable by the logic subsystem. The instructions are executable to obtain one or more metrics computed based upon an energy and a force determined for a material using a universal machine learning force field (MLFF), and based at least upon the one or more metrics, determine to use one of a first computational method or a second computational method to compute a predicted physical property value for the material.

Description

METRICS-BASED COMPUTATIONAL METHOD SELECTION FOR THE PREDICTION OF A PHYSICAL PROPERTYBACKGROUND

[0001] Information on the physical properties of materials can be researched using a wide variety of sources. For example, physical properties for specific materials can sometimes be found in the scientific literature. Further, information on the physical properties of materials has been compiled in various online databases.SUMMARY

[0002] Examples are disclosed that relate to the selection of a method to compute a physical property based upon metrics obtained using a universal machine learning force field. One disclosed example provides a computing system comprising a logic subsystem, and a storage subsystem comprising instructions executable by the logic subsystem. The instructions are executable to obtain one or more metrics computed based upon an energy and a force, and optionally a stress determined for a material using a universal machine learning force field (MLFF) , and based at least upon the one or more metrics, determine to use one of a first computational method or a second computational method to compute a predicted physical property value for the material.

[0003] This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject  matter. Furthermore, the claimed subject matter is not limited to implementations that solve any or all disadvantages noted in any part of this disclosure.BRIEF DESCRIPTION OF THE DRAWINGS

[0004] FIG. 1A schematically shows an example system for determining a physical property by using metrics obtained from an output of a universal machine learning force field (uMLFF) to compute a predicted physical property.

[0005] FIG. 1B schematically shows the system of FIG. 1 used to determine a first predicted physical property that is missing from a materials database.

[0006] FIG. 1C schematically shows the system of FIG. 1 used to determine a second predicted physical property that is missing from a materials database.

[0007] FIG. 2 schematically shows a flow diagram that illustrates an example method for the determination of a computational method to predict the thermal conductivity of a material based on one or more metrics.

[0008] FIG. 3 schematically shows a flow diagram that illustrates an example method for the determination of a computational method to predict the electronic structure of a material based on one or more metrics.

[0009] FIG. 4 schematically shows another example system for determining a physical property by using metrics obtained from an output of a universal machine learning force field to compute a predicted physical property.

[0010] FIG. 5 shows an example parity plot between a thermal conductivity predicted by using a uMLFF and equilibrium molecular dynamics (EMD) and reference thermal conductivity values for 24 example materials.

[0011] FIG. 6 shows an example parity plot between a thermal conductivity predicted by using a uMLFF and three-phonon Boltzmann transport equation (BTE) and reference thermal conductivity values for 120 example materials.

[0012] FIG. 7 shows example phonon spectra and density of states (DOS) for bulk silicon as determined using different methods.

[0013] FIGS. 8A-8B show example comparisons of experimentally determined and EMD or DFT-BTE predicted thermal transport properties for bulk silicon and silicon clathrates.

[0014] FIGS. 9A-9C show example thermal conductivity predictions for various silicon clathrates as a function of correlation time.

[0015] FIG. 10 shows an example summary table of predicted and experimentally determined thermal conductivity values of various complex materials.

[0016] FIG. 11 shows a block diagram of an example computing system.DETAILED DESCRIPTION

[0017] Information on the physical properties of materials can be used to identify suitable materials for the development of new technologies. For example, the thermal conductivity value of a material can inform whether it is a suitable material for use in an energy-efficient building, or for use in heat dissipation in electronic devices. As mentioned above, the physical properties of some materials are currently available in scientific literature or compiled into searchable databases. The databases and literature can be populated by experimentally determined data. However, experimentally determined data may not be available for some materials of possible interest. As such, searching for suitable materials for a desired application by filtering  using experimentally determined properties may fail to discover many potentially well-suited materials.

[0018] To address this issue, experiments can be performed to measure physical properties of materials for which experimental data is currently missing. However, performing such experiments for known materials with incomplete physical property data sets would be impracticable on a reasonable time scale. Further, some materials are difficult to synthesize in suitably pure form and / or at a suitable quantity for performing such experiments (e.g. diamond and boron arsenide, as examples) . Also, processing these materials into useable forms can introduce defects and impurities, which can alter (e.g. reduce) their real magnitudes in pure and single-crystalline format. thermal conductivities. Such problems may be even larger for the universe of unknown materials yet to be discovered and studied.

[0019] Thus, another approach to expanding the amount of physical property data available for known and to-be-discovered materials is to use computational methods to predict the physical properties for a material. For many physical properties and materials, suitably accurate computational methods exist for predicting the physical property values. As an example, numerous theoretical methods have been employed to predict thermal conductivity and bolster the understanding of heat transfer mechanisms, ranging from early semi-empirical models such as Callaway and Allen-Feldman (AF) models to solving the Boltzmann transport equation (BTE) , running approach-to-equilibrium molecular dynamics (AEMD) , nonequilibrium molecular dynamics (NEMD) , or applying the Green-Kubo formalism using classical molecular dynamic (MD) or first-principles MD (FPMD) . Theoretical investigation of thermal conductivity has advanced the discovery of materials significantly, with the re-discovery of Boron arsenide (BAs) as an ultra-high thermal conductivity material being a notable example.  However, the theoretical study of thermal conductivity involves a delicate balance between the accuracy of the models and the computational feasibility. While classical force fields have been developed for decades to predict thermal conductivity for simple crystals (e.g., diamond, silicon) , low-dimension and complex systems (e.g., silicon nanowires, carbon nanotubes, amorphous silicon and porous silicon) , as well as liquids, the lack of accurate force fields limited its applicability to only simple systems with well-tested force fields, such as Tersoff for silicon and TIP3P for water. Recently, advancements in computational techniques and the growing power of supercomputers offer promising avenues for further progress in this field by approaching the problem using BTE and AIMD with density functional theory (DFT) . However, due to the demanding computational costs, these studies have been limited to case-by-case. The largest high-throughput calculation with BTE at DFT level in a single study was only 92 perovskites. Further, the largest high-throughput calculation with FPMD was only performed for 24 materials.

[0020] In recent years, time has witnessed the advancement of artificial intelligence (AI) in the field of computational chemistry and materials science, such as machine learning force fields, machine-learned exchange-correlation functional, quantum Monte Carlo, excitation spectroscopies, and inverse design for new material generation. To accelerate the prediction of physical properties of materials using AI, current approaches involve training end-to-end models to predict physical properties. For example, the prediction of thermal conductivity, the training is based on the macroscopic observables of materials, and training machine-learned force fields from a limited set of DFT calculations for specific materials, followed by the application of BTE or MD-based methods. However, little success has been made for transferability  to new, unseen material discoveries due to the lack of massive labels across compositions and operation conditions such as temperature and pressure.

[0021] To address the above challenges, examples are disclosed that relate to the automated selection of a computational method to use for the prediction of a realistic physical property value of a material. Briefly, the disclosed examples utilize a universal machine learning force field (uMLFF) for the determination of an energy, a force, and optionally a stress for a material. The uMLFF may comprise, for example, a deep-learning model for materials simulation such as MatterSim, described in Yang, Han, et al."MatterSim: A Deep Learning Atomistic Model Across Elements, Temperatures and Pressures. " arXiv preprint arXiv: 2405.04967 (2024) . Such a uMLFF can be trained using supervised training methods based upon a materials database having a wide range of known structures of materials with elements from across the periodic table with known energies, forces, and optionally stresses at various temperatures and pressures. The energy, force, and optionally stress, determined by the uMLFF then can be used to determine one or more metrics for the material. Based upon one or more metrics, the disclosed examples determine whether to use one of at least a first computational method or a second computational method to compute a predicted physical property value for the material. For the prediction of thermal conductivity values, the determination may be generalized to incorporate non-equilibrium molecular dynamics (NEMD) , approach-to-equilibrium molecular dynamics (AEMD) , and Wigner transport equation (WTE) for amorphous and low-dimensional materials, as well as liquids.

[0022] As described in more detail below, the disclosed examples may allow the selection of an efficient computational method that will yield a suitably accurate and realistic value for a physical property for a material. This provides a solution to  simulating a physical property (e.g. thermal conductivity) with the same accuracy of quantum mechanical simulations, but 100 times faster than the traditional method. The selection of a computational method for a material can depend on the physical properties of each material. For example, a predicted thermal conductivity of a material can be calculated using EMD or BTE methods. The efficiency and accuracy of the thermal conductivity prediction from each of these computational methods will vary depending on the energy, force, and optionally stress, of the material as computed using the uMLFF. In some cases, a method can yield a result that may be somewhat less accurate, but still suitably accurate, and the method is more computationally efficient. In such instances, the more efficient method can be chosen. On the other hand, where a more efficient method is likely to yield an inaccurate predicted value, a less efficient but more accurate method can be used.

[0023] Returning again to the example of thermal conductivity, the metrics computed using the output of a uMLFF may comprise phonon modes, anharmonicity score, Debye temperature, and 3-phonon scattering phase space volume (V_3ph) . Analysis of these metrics can be used to determine whether to use equilibrium molecular dynamics (EMD) , 3-phonon Boltzmann transport equation with Bose-Einstein distribution (BTE) , or the Boltzmann transport equation with temperature-dependent force constants (BTE-FC) to compute the thermal conductivity. While the BTE method may be most computationally efficient of these methods, it may not produce suitably accurate, and thus realistic, results for some materials, such as where the metrics determined from the output of the uMLFF indicate that the phonon modes for the material comprise imaginary frequencies. Therefore, in such instances, a less efficient but more accurate method such as EMD may be used to compute the predicted  thermal conductivity based upon this metric. This and other examples are discussed in more detail below.

[0024] FIG. 1A schematically shows an example computing system 100 configured to compute a physical property for a material. In this example, the computing system 100 is in communication with a client computing device 102, for example, via an application programming interface (API) of the computing system 100. The client device is in communication with a materials database 104. The materials database 104 includes a compilation of physical property data for a collection of materials, shown here as material 1 through material n, where n is an arbitrary integer. A materials database may have physical property information for hundreds of thousands, or even millions, of materials.

[0025] More specifically, the materials database 104 comprises identifications of the materials (e.g. “material 1” ) , structural information comprising atomic coordinates specifying the 3-dimensional location and identification of each atom of the material, (e.g. “mat1 coord. ” ) , and a collection of available material properties for the listed materials (e.g. property 1, property n …) . In materials database 104, experimentally determined values of each property may or may not be available. In FIG. 1, the materials database 104 comprises values of a property 1 for a material 1 and a material 2, but the value of property 1 for a material n is not available. Further, the materials database 104 comprises values of a property n for material 1 and material n, but not for material 2. Thus, a search for materials with potentially suitable values for property n may miss material 2, even though material 2 may actually have a value of physical property n within the range of suitable values.

[0026] To help fill in the missing physical property values of the materials database 104, the client device 102 can provide the atomic coordinates for a material  with a missing physical property value to the computing system 100 for computation of a predicted value for the physical property. As mentioned above, the computing system 100 comprises a uMLFF 106 configured to output values for energy, force, and optionally stress, for a material based at least upon the atomic coordinates of the material.

[0027] The computing system further includes a metric module 108 configured to determine one or more metrics based upon the energy, force, and optionally stress, for the material determined by the uMLFF. Different metrics can be determined for different physical properties in some examples, as indicated by “Property 1 metrics” and “Property n metrics” in FIG. 1A. Examples of metrics that can be computed for various physical properties are described below. The computing system 100 further includes a computation selection module 110. The computation selection module is configured to analyze the metrics obtained using metric module 108 against one or more metric thresholds for the property ( “property thresholds” ) . Based upon the analysis, computation selection module 110 determines a selected computational method to use to compute the predicted physical property value for the physical property. For example, the computation selection module may compare each of one or more of the property 1 metrics to a corresponding threshold to determine which computational method to select. Different metrics can be used for different physical properties, as shown by “Property 1 thresholds” and “Property n thresholds” in FIG. 1A.

[0028] The physical property value is then computed using the selected computational method. For example, for the computation of physical property 1 at 114, one of computational method 1 at 116 through computational method N at 118 is used based on the selection results of the computation selection module 110, where N is an integer representing a number of available computational methods for the physical  property. Likewise, for the computation of physical property N at 120, one of computational method 1 122 through computational method N 124 can be used based upon the selection results of the computation selection module 110. The computed physical property value may then be communicated back to client device 102 from computing system 100. The result may then be added to materials database 104.

[0029] FIG. 1B schematically shows the computing system of FIG. 1A used to predict the missing property n value for material 2. Property n may comprise, for example, thermal conductivity. As shown in bold on FIG. 1B, the client device 102 provides the atomic coordinates of Material 2 (Mat2) from materials database 104 to computing system 100 for the selection of a computational method to compute a predicted physical property value (e.g. thermal conductivity) . The atomic coordinates can take various forms based upon type of material. For example, the atomic coordinates can take the form of the locations and identifications of atoms in a unit cell for crystalline or polycrystalline solids. As another example, the atomic coordinates can take the form of a larger, randomly arranged set of locations and identification of atoms in an amorphous solid or liquid.

[0030] Based on the Mat2 atomic coordinates, the uMLFF 106 outputs the energy, force, and optionally stress computed for Mat2. The energy, force, and optionally stress for Mat2 are then input into the metric module 108 of computing system 100. The metric module determines the property n metrics for Mat2 based on the energy, force, and optionally stress for Mat2 output from the uMLFF. The computation selection module 110 then analyzes the metrics from metric module 108 against Property n thresholds. For example, the computation selection module can comprise thresholds for determining a computation to use for a predicted thermal conductivity, as described in further detail in FIG. 2. In the example shown in FIG. 1B  for Mat2, computational method n is chosen for the computation of a predicted value for physical property n. The predicted value for physical property n for Mat2 is then communicated back to client device 102 from the computing system 100. The predicted property n value is entered into materials database 104 from client device 102.

[0031] FIG. 2 schematically shows an example workflow that can be followed by the computation selection module of FIG. 1B for the selection of a computational method for the prediction of a thermal conductivity value based upon a series of metrics thresholds. In this example, the computational methods comprise EMD, 3-phonon Boltzmann transport equation with Bose-Einstein distribution (3ph-BTE) , and 3-phonon Boltzmann transport equation with temperature-dependent force constants (3ph-BTE FC) . The metrics 200 in this example comprise phonon modes, anharmonicity score, Debye temperature, and 3-phonon scattering phase space volume (V_3ph) , and are computed based upon outputs of the uMLFF. Example methods for computing metrics are described in more detail in the Experimental section below.

[0032] One of the metric thresholds may comprise, at 202, whether the phonon modes comprise an imaginary frequency. If the phonon modes do not comprise an imaginary frequency, the workflow leads to another metric threshold 204, which determines whether the anharmonicity score meets a threshold anharmonicity score. If the anharmonicity score does not meet a threshold anharmonicity score, then the workflow leads to EMD as the computational method at 206 to compute the accurate predicted thermal conductivity value 212 of the material, shown by a dotted line arrow in FIG. 2.

[0033] On the other hand, if it is determined, at metric threshold 204, that the anharmonicity score meets a threshold anharmonicity score, the workflow leads to another metric threshold 208, which determines whether the 3-phonon scattering phase  space volume V_3ph meets a threshold space proportion. If the V_3ph meets the threshold space proportion for the material, then it is determined that the 3ph-BTE computational method at 210 is a suitably accurate and efficient method, and is used to compute the predicted thermal conductivity 212, as shown by the dotted line arrow in FIG. 2. Referring back to metric threshold 208, if the V_3ph does not meet the threshold space proportion for the material, then the workflow determines that the 3-ph BTE method would produce a less accurate result. Therefore, the EMD computational method 206 is used to compute a more accurate predicted thermal conductivity 212.

[0034] Referring back to metric threshold 202 in FIG. 2, if the phonon modes do comprise an imaginary frequency, then the workflow leads to metric threshold 214, which determines whether the quantum effect meets a threshold Debye temperature. If the quantum effect does not meet the threshold Debye temperature, then EMD is the chosen computational method to predict the realistic thermal conductivity value of the material. If the quantum effect does meet the threshold Debye temperature, then the 3-ph BTE-FC method 216 is determined to be the more efficient and suitably accurate method and is used to predict the thermal conductivity value of the material, as shown by the dotted arrow leading to the predicted thermal conductivity value 212.

[0035] In this way, the computation selection module in this example can compare the metrics computed by the metric module to metric threshold (s) to determine which computational method may be used to predict the thermal conductivity value for a material.

[0036] Continuing, FIG. 1C schematically shows the computing system of FIG. 1A used to predict the missing property 1 value for a material n (Matn) , and to subsequently enter the property value into a materials database 104. Property 1 may comprise, for example, an electronic structure value (s) of a material. As described  above and shown in FIG. 1C, the atomic coordinates of Matn are retrieved from the materials database 104 by client device 102 and communicated to the computing system 100. As shown in bold on FIG. 1C, the client device 102 provides the atomic coordinates of Matn from materials database 104 to computing system 100 for the selection of a computational method to compute a predicted physical property value (e.g. electronic structure) .

[0037] Based on the Matn atomic coordinates, the uMLFF 106 outputs the energy, force, and optionally stress computed for Matn. The energy, force, and optionally stress for Matn are then input into the metric module 108 of computing system 100. The metric module determines the property 1 metrics for Matn based on the energy, force, and optionally stress output from the uMLFF 106. The computation selection module 110 then analyzes the metrics from metric module 108 against property 1 thresholds.

[0038] The computation selection module can comprise thresholds for determining a computation to use for a predicted electronic structure value (s) , as described in further detail in FIG. 3. In the example shown in FIG. 1C for Matn, computational method 1 is chosen for the computation of a predicted value for physical property 1. The predicted value for physical property 1 for Matn is then communicated back to client device 102 from the computing system 100. In this example, the predicted property 1 value is entered into materials database 104 from client device 102.

[0039] FIG 3 schematically shows an example workflow that can be followed by the computation selection module of 1C for the selection of one of a first computational method or a second computational method for the prediction of the electronic structure value (s) of a material. In this example, the computational methods comprise EMD and a displacement-based computational method. The metrics 300 in  this example comprise phonon modes and anharmonicity score, and are computed based upon outputs of the uMLFF, as described above.

[0040] One of the metric thresholds may comprise, at 302, whether the phonon modes comprise an imaginary frequency. If the phonon modes comprise an imaginary frequency, then EMD is the chosen computational method at 306 to compute the more accurate predicted electronic structure value (s) 310 of a material, as shown by a dotted arrow in FIG. 3. If, at metric threshold 302, the phonon modes do not comprise an imaginary frequency, the workflow leads to another metric threshold 304, which determines whether the anharmonicity score meets a threshold anharmonicity score.

[0041] If the anharmonicity score does not meet a threshold anharmonicity score at metric threshold 304, then EMD is the chosen computational method at 306 to compute the predicted electronic structure value (s) of the material 310, as indicated by a dotted arrow in FIG. 3. Referring back to metric threshold 304, if the anharmonicity score for a material does meet a threshold anharmonicity score, then a displacement-based computational method 312 is used to compute the predicted electronic structure value (s) 310 of the material. In this way, the computation selection module in this example compares the metrics computed by the metric module to metric thresholds to determine which computational method may be used to predict the electronic structure value (s) for a material more accurately.

[0042] FIGs. 1A-1C show an example computing system in which the uMLFF, which outputs the energy, force, and optionally stress computed for a material, is incorporated in the computing system. However, the uMLFF may also be external to the computing system, as shown in FIG. 4. More particularly, FIG. 4 shows an example computing system 400 in which a uMLFF 402 operating on remote computing system  separate from the computing system 400 is used to provide the energy, force, and optionally stress.

[0043] In this example, the computing system 400 is in communication with a client device 404 via an application programming interface (API) of the computing system 400. The client device is in communication with a materials database 406, computing system 400, and the uMLFF 402. Like the materials database described for FIG. 1A, the materials database 406 includes a compilation of physical property data for a collection of materials, shown in FIG. 4 as material 1 through material n, where n is an arbitrary integer. Experimentally determined values of each property may or may not be available in materials database 406.

[0044] In the same manner as described for FIGS. 1A-1C, in order to help fill in the missing physical property values of the materials database 406, the client device 404 can search the materials database 406 for the atomic coordinates of a material. The client device 404 can then provide the atomic coordinates for the material to the uMLFF for an output of an energy, a force, and optionally a stress of the material, which may then be communicated back to client device 404. The client device 404 receives the output and enters the energy, force, and optionally stress into computing system 400 for the computation of a physical property value. In other examples, the computing system 400 can communicate with the uMLFF 402 on the remote computing system.

[0045] Similarly to the computing systems in FIGS. 1A-1C, the computing system 400 comprises metric module 408 configured to determine or more metrics based upon an energy, a force, and optionally a stress for the material determined by the external uMLFF 402 and communicated to computing system 400 via client device 404. Different metrics can be determined for different physical properties in some examples, as described above.

[0046] The computing system 400 shown in FIG. 4 further includes computation selection module 410. As described above, the computation selection module is configured to analyze the metrics obtained using metric module 408 against one or more metric thresholds for the property ( “property thresholds” ) . Based upon the analysis, computation selection module 410 determines a selected computational method to use to compute the predicted physical property value for the physical property, as described above for FIGs 1A-1C. The physical property value is then computed using the selected computational method.

[0047] For example, for the computation of physical property n at 412, one of computational method 1 at 414 through computational method n at 416 is used based on the selection results of the computation selection module 410, where n is an integer representing a number of available computational methods for the physical property. In a similar manner as described in FIGs. 1A-1C, the computed physical property value may then be communicated back to client device 404 from computing system 400. The result may then be added to materials database 406. Thus, a predicted physical property value can be computed using the computing system 400 with values from the external uMLFF 402.

[0048] The examples described above can help to balance computational efficiency and accuracy when predicting physical properties of materials. This can help to populate a database with suitably realistic values for a physical property more efficiently than by always using a more accurate, but less efficient method of computing the physical property, and more accurately than by always using a more efficient, but less accurate method. Results of experiments to compute predicted thermal conductivity values are summarized below.

[0049] EXPERIMENTAL RESULTS

[0050] I. Benchmark Performance for EMD and Three-Phonon BTE

[0051] A. Equilibrium Molecular Dynamics

[0052] The examples described above can help to balance computational efficiency and accuracy when predicting physical properties of materials. This can help to populate a database with suitably realistic values for a physical property more efficiently than by always using a more accurate, but less efficient method of computing the physical property, and more accurately than by always using a more efficient, but less accurate method. Results of experiments to compute predicted thermal conductivity values are summarized below. To illustrate the accuracy of applying the uMLFF MatterSim to predict thermal conductivity κ using equilibrium molecular dynamics (EMD) , the performance of the application of MatterSim is first benchmarked as disclosed herein with respect to 24 materials reported in literature. As shown in FIG. 5, the predictions by MatterSim record an R2 of 0.89 and mean absolute error (MAE) of below 3Wm-1K-1 when comparing against the experimental measurements compiled from the literature. Such an agreement is favorable when considering order of magnitude differences of κprocessed by these materials, as well as the underlying accuracy of force predictions for such materials. When κ is predicted with MatterSim and EMD, a uniformed protocol is employed, such that three simulations with hundreds of atoms in size, as well as nanosecond (s) in length, are performed for each material. With sufficiently converged simulations, κ predictions can be made readily available without imposing any post corrections. Fig 5 reveals that κ predictions for materials across 20 compositions (H, Mg, O, Li, Na, F, Cl, K, Rb, Zn, Ba, Cd, Ca, Br, Sn, Se, Ga, Ag, I, Sr) are suitably accurately predicted at the first principle level.

[0053] B. Three-Phonon Boltzmann Transport Equation

[0054] After gauging the accuracy of κ predictions with EMD, a thorough and systematic benchmark on the κ predictions is performed when coupling MatterSim with three-phonon Boltzmann transport equation (BTE) . It is worth noting that compared to EMD, BTE captures nuclear quantum effects and is much more computationally affordable. Thus, BTE can be the method of choice for materials with high thermal conductivity.

[0055] A set of 120 materials is compared, of which the thermal conductivities are either probed experimentally or validated with first-principles-based BTE calculations. To standardize κ values across different orders of magnitude, the error metric is reported as the average factor difference (AFD) . As shown in FIG. 6, the κpredictions for the entire dataset obtained from BTE yield an AFD score of 1.33. This performance is on par with, if not exceeding, end-to-end models in the literature. The majority of κ predictions fall within the ±50%margin, as shown in FIG. 6, and are evenly distributed along the parity line. This uniform distribution is consistent across materials with different space group symmetries, indicating the robustness and unbiasedness of κ predictions obtained by combining MatterSim and BTE. However, when seeking κ reference values to compare against the computationally predicted values, it is often a challenge to establish a fair assignment, as the forces used in constructing the 2nd and 3rd force constants are modeled approximately at the level of Perdew-Burke-Ernzerhof (PBE) functional, and the BTE is known to fall short for materials with non-negligible contributions to κ from higher-order scattering processes.

[0056] Following the benchmark of the κ predictions for MatterSim against commonly known crystals, its application in modeling κ for various complex materials is considered. In this disclosure, κ is examined in materials such as bulk silicon, type I  silicon clathrate (Si46) and its barium intercalated counterpart (Ba8Si46) , type II silicon (Si136) , and amorphous silicon (aSi) . These substances exhibit a wide range of κ values and are promising candidates for semiconductors, solar cells, and thin film transistors. By delivering accurate κ predictions for these materials, the transferability of the potential is demonstrated for various systems across different phases.

[0057] To better understand the harmonic properties of the materials under study, the phonon spectra and total density of states (DOS) for bulk silicon, is computed using the MatterSim, as depicted in FIG. 7. The MatterSim results are compared with results derived from the empirical bond-order Tersoff potential (see e.g. Tersoff, J.: New empirical approach for the structure and energy of covalent systems. Physical review B 37 (12) , 6991 (1988) ) , which has demonstrated its effectiveness in capturing the thermal transport properties across amorphous, bulk, and nanostructured silicon. For the case of bulk silicon, the phonon spectra and DOS are further supplemented with harmonic force constants, computed using Vienna Ab Initio Simulation Package (VASP) under the same density functionals and computational settings as in the training dataset. The comparison, as illustrated in FIG. 7 shows a close agreement argument between the MatterSim calculations, Tersoff potential predictions, and DFT results for bulk Si. More specifically, FIG. 7 shows example phonon spectra for PBE, MatterSim and Tersoff respectively at 702, 704, and 706, and densities of state respectively at 712, 714, and 716. Close agreement was also seen for Si46 and Si136. Such an agreement emphasizes the reliability of the MatterSim computed phonon spectra and DOS predictions when employing consistent initial structures, supercell dimensions, and reciprocal q-grid settings.

[0058] Referring to FIG. 7, in the case of bulk silicon, there is a general agreement in the optical branches of the phonon spectra and DOS between MatterSim  and DFT, with some minor deviations observed in the transverse acoustic branches where MatterSim slightly overestimates the frequencies. Additionally, a comparison between the Tersoff potential and MatterSim reveals subtle differences in the optical branches, with the Tersoff potential consistently predicting higher phonon frequencies for all the complex materials studied. The Tersoff potential also predicts higher transverse acoustic frequencies than both MatterSim and DFT for bulk silicon. These discrepancies can be attributed to the inherently short-range nature of the Tersoff potential, which may result in an overestimation of κ.

[0059] Upon validating the phonon spectra and DOS, a comparison is made between phonon group velocities, lifetime under RTA (relaxation time approximation) , and mean free path, as well as κ as a function of frequency and temperature for bulk silicon, type I, and type II silicon clathrate. FIGS. 8A –8B illustrate thermal transport properties for bulk Si, Si46 and Si136. As shown respectively at 802, 804, and 806 in FIG. 8A, for all three materials, the temperature dependence of κ exhibits an inversely proportional trend, which is in line with T-1 characteristic for crystals. Notably, the κ predictions for bulk silicon closely match experimental results up to approximately 700 K. Beyond this temperature, minor discrepancies emerge, for example, due to the ineligible impact of high-order anharmonicity. Such an outstanding agreement reassures the accuracy of the potential and justifies the computational settings. Not only does the MatterSim provide robust κ predictions for bulk silicon, it also accurately predicts κ for different silicon-based clathrate cages. Triangles 810 in FIG. 8A denote the κ obtained with EMD using MatterSim. Circle, square and diamond symbols 812 represent DFT-BTE calculations from the literature (see, e.g.,  V.J., Karttunen, A.J.: Ab initio studies on the lattice thermal conductivity of silicon clathrate frameworks ii and viii. Physical Review B 93 (2) , 024307 (2016) ; Norouzzadeh, P.,  Krasinski, J.S., Tadano, T.: Thermal conductivity of type-i, type-ii, and type-viii pristine silicon clathrates: A first-principles study. Physical Review B 96 (24) , 245201 (2017) ; Chen, C., Zhang, Z., Chen, J.: Revisit to the impacts of rattlers on thermal conductivity of clathrates. Frontiers in Energy Research 6, 34 (2018) ) .

[0060] As depicted in FIG. 8B, the κ from BTE converges at 145.1, 49.5 and 49.8 W m-1K-1 at 300 K for silicon 822, Si46 824 and Si136 826. These κ values, alongside with the one-third reduction in κ from Si to Si46 and Si136, are consistent with values reported by various thermal conductivity studies of silicon clathrate compounds with DFT-based BTE calculations (see e.g.  V.J., Karttunen, A.J.: Ab initio studies on the lattice thermal conductivity of silicon clathrate frameworks ii and viii. Physical Review B 93 (2) , 024307 (2016) ; see also Norouzzadeh, P., Krasinski, J.S., Tadano, T.: Thermal conductivity of type-i, type-ii, and type-viii pristine silicon clathrates: A first-principles study. Physical Review B 96 (24) , 245201 (2017) 39; see also Chen, C., Zhang, Z., Chen, J.: Revisit to the impacts of rattlers on thermal conductivity of clathrates. Frontiers in Energy Research 6, 34 (2018) ) . Despite the underrepresentation of cage-like structures in the training dataset, the excellent agreement between the disclosed results and literature results exemplifies the transferability of MatterSim from modeling bulk silicon to silicon-based nanostructures. This versatility is notable, given the complexity of obtaining accurate κ values from BTE, which depends on the accurate predictions of second and third derivatives of interatomic forces, a task that is inherently more challenging than predicting forces itself.

[0061] However, while the MatterSim potential reasonably captures the overall thermal conductivity values for Si46 and Si136, it may not suitably distinguish subtle differences in thermal transport between the two clathrates. As investigated by Norouzzadeh et al. (see, e.g. Norouzzadeh, P., Krasinski, J.S., Tadano, T.: Thermal  conductivity of type-i, type-ii, and type-viii pristine silicon clathrates: A first-principles study. Physical Review B 96 (24) , 245201 (2017) ) , different crystal structure and arrangement of silicon atoms give rise to different phonon lifetime, which results in different lattice thermal conductivity. However, it was found that the group velocity norm and phonon lifetime as well as mean free path are indistinguishable from Si46 to Si136, which accounts for the minor differences in κ predicted by MatterSim. Considering an under-representation of cage-like structures in the training dataset, such a lack of differentiability is unsurprising and can be improved by imposing pretrain-finetune techniques.

[0062] C. Thermal Conductivity for Complex Systems

[0063] After demonstrating the transferability of MatterSim across bulk and nanostructured silicon, the limit is further pushed by investigating κ for Ba8Si46 and amorphous silicon (aSi) . Unlike crystalline structures, amorphous materials exhibit shorter mean free paths, leading to the breakdown of quasi-particle assumptions inherent in BTE. Consequently, molecular dynamics (MD) may be a preferred method for calculating lattice thermal conductivity in disordered and amorphous materials. As a similar example, when barium is encapsulated as an alkaline metal host, the anharmonicity of Ba8Si46 increases noticeably even at room temperature (not shown in FIGS. 9A-9C) . Thus, MD is a more proper method over three-phonon BTE as it includes full order of anharmonicity. EMD is therefore employed as an approach to model κ for aSi. FIGS. 9A-9C illustrates the running thermal conductivity as a function of correlation time τ for each of aSi, , Si46, and Si136, respectively. Shaded regions imply the standard error derived from three separate nanosecond MD simulations. The predicted κ values are shown respectively by lines 902, 912, and 922 in FIGS. 9A-9C. The three separate simulations for each of aSi, Si46, and Si136 are depicted by lines 904,  914, and 924 in FIGS 9A-9C, respectively. The long, converged tails exhibited in the EMD calculations indicates the choice of correlation time to be sufficient and suitable.

[0064] In terms of the shear values, MatterSim records a κ of 1.6 ± 0.2 Wm-1 K-1 for aSi. As noted in FIG. 10, such a κ prediction is in reasonable agreement with those from other machine learning models developed for silicon in different phases (see e.g., He, Y., Donadio, D., Galli, G.: Heat transport in amorphous silicon: Interplay between morphology and disorder, Applied Physics Letters 98 (14) (2011) ( “He” in FIG. 10); see also, Qian, X., Peng, S., Li, X., Wei, Y., Yang, R.: Thermal conductivity modeling using machine learning potentials: application to crystalline and amorphous silicon. Materials Today Physics 10, 100140 (2019) ( “Qian” in FIG. 10) ; see also, Li, R., Lee, E., Luo, T.: A unified deep neural network potential capable of predicting thermal conductivity of silicon in different phases. Materials Today Physics 12, 100181 (2020) ( “Li” in FIG. 10) ; see also, Barbalinardo, G., Chen, Z., Lundgren, N.W., Donadio, D.: Efficient anharmonic lattice dynamics calculations of thermal transport in crystalline and disordered solids. Journal of Applied Physics 128 (13) (2020) ( “Barbalinardo” in FIG. 10) ; see also, Glassbrenner, C.J., Slack, G.A.: Thermal conductivity of silicon and germanium from 3 k to the melting point. Physical review 134 (4A) , 1058 (1964) , Stuckes, A.D.: The thermal conductivity of germanium, silicon and indium arsenide from 40 c to 425 c. Philosophical Magazine 5 (49) , 84–99 (1960) , Shanks, H., Maycock, P., Sidles, P., Danielson, G.: Thermal conductivity of silicon from 300 to 1400 k. Physical Review 130 (5) , 1743 (1963) , Capinski, W., Maris, H., Bauser, E., Silier, I., Asen-Palmer, M., Ruf, T., Cardona, M., Gmelin, E.: Thermal conductivity of isotopically enriched si. Applied physics letters 71 (15) , 2109–2111 (1997) , Kremer, R., Graf, K., Cardona, M., Devyatykh, G., Gusev, A., Gibin, A., Inyushkin, A., Taldenkov, A., Pohl, H.-J.: Thermal conductivity of isotopically  enriched 28si: revisited. Solid state communications 131 (8) , 499–503 (2004) , Inyushkin, A., Taldenkov, A., Gibin, A., Gusev, A., Pohl, H.-J.: On the isotope effect in thermal conductivity of silicon. physica status solidi (c) 1 (11) , 2995–2998 (2004) (“Glassbrenner, Stuckes, Shanks, Capinski, Kremer, Inyushkin” in FIG. 10) ; see also, Chen, C., Zhang, Z., Chen, J.: Revisit to the impacts of rattlers on thermal conductivity of clathrates. Frontiers in Energy Research 6, 34 (2018) ( “Chen” in FIG. 10) ; see also, Norouzzadeh, P., Krasinski, J.S., Tadano, T.: Thermal conductivity of type-i, type-ii, and type-viii pristine silicon clathrates: A first-principles study. Physical Review B 96(24) , 245201 (2017) ( “Norouzzadeh” in FIG. 10) ; see also,  V.J., Karttunen, A.J.: Ab initio studies on the lattice thermal conductivity of silicon clathrate frameworks ii and viii. Physical Review B 93 (2) , 024307 (2016) ( in FIG. 10) ; see also, Wieczorek, L., Goldsmid, H., Paul, G.: In Thermal Conductivity 20, edited by DPH Hasselman and JR Thomas. Plenum, New York (1989) , Kuo, B., Li, J., Schmid, A.: Thermal conductivity and interface thermal resistance of si film on si substrate determined by photothermal displacement interferometry. Applied Physics A 55, 289–296 (1992) , Cahill, D.G., Katiyar, M., Abelson, J.: Thermal conductivity of a-si: H thin 40 films. Physical review B 50 (9) , 6077 (1994) , Wada, H.W.H., Kamijoh, T.K.T.: Thermal conductivity of amorphous silicon. Japanese journal of applied physics 35 (5B) , 648 (1996) , Moon, S., Hatano, M., Lee, M., Grigoropoulos, C.P.: Thermal conductivity of amorphous silicon thin films. International Journal of Heat and Mass Transfer 45 (12) , 2439–2447 (2002) , Zink, B., Pietri, R., Hellman, F.: Thermal conductivity and specific heat of thin-film amorphous silicon. Physical review letters 96(5) , 055902 (2006) , Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Chiarotti, G.L., Cococcioni, M., Dabo, I., et al.: Quantum espresso: a modular and open-source software project for quantum simulations of  materials. Journal of physics: Condensed matter 21 (39) , 395502 (2009) , Yang, H.-S., Cahill, D.G., Liu, X., Feldman, J., Crandall, R., Sperling, B., Abelson, J.: Anomalously high thermal conductivity of amorphous si deposited by hot-wire chemical vapor deposition. Physical Review B 81 (10) , 104203 (2010) ( “Wieczorek, Kuo, Cahill, Wada, Moon, Zink, Giannozzi, Yang” in FIG. 10) ) . Different factors, such as the accuracy of the force-field, the size of the amorphous sample, and the quench rate used in preparing the sample, can have subtle impacts on the resultant κ prediction. As shown in FIG. 10, different experimental conditionals can result in a wide range of κ measurements. To achieve statistically significant results with MD, amorphous silicon samples in hundreds to thousands of atoms, as well as with simulation runs beyond pico-second in length are often required. Thus, machine learning or empirical force fields may be applied in such studies.

[0065] Although the κ predictions for aSi discussed above are consistent with values derived from non-propagating vibrational modes using the AF theory (see, e.g., Larkin, J.M., McGaughey, A.J.: Thermal conductivity accumulation in amorphous silica and amorphous silicon. Physical Review B 89 (14) , 144303 (2014) ; also see Isaeva, L., Barbalinardo, G., Donadio, D., Baroni, S.: Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach. Nature communications 10 (1) , 3853 (2019) ; also see and Fiorentino, A., Pegolo, P., Baroni, S.: Hydrodynamic finite-size scaling of the thermal conductivity in glasses. arXiv preprint arXiv: 2303.07010 (2023) ) , the predictions represent a conservative agreement due to finite-size effects. In comparison, the EMD simulations are also conducted for aSi samples with merely 512 atoms. When evaluating the κ predicted by MatterSim for aSi, it can be seen in FIG. 10 that κ is sensitive to sample size, as κ increases from 1.2 ± 0.2 to 1.6 ± 0.2 W m-1K-1 as the number of atoms triples in between the two amorphous samples. MatterSim’s robust  performance across various forms of silicon, from bulk to nanostructured and amorphous, suggests its potential in providing valuable insights into the lattice thermal conductivity of aSi with varying sizes, quenching rates, and hydrogen concentrations.

[0066] For the κ prediction in FIG. 10 for Ba8Si46, MatterSim yields a result of 1.2 ± 0.4 W m-1K-1, which aligns well with measurements done on Ba8Si46 single crystals using scanning thermal microscopy experiments (see, e.g., Pailhès, S., Euchner, H., Giordano, V.M., Debord, R., Assy, A., Gomès, S., Bosak, A., Machon, D., Paschen, S., De Boissieu, M.: Localization of propagative phonons in a perfectly crystalline solid. Physical review letters 113 (2) , 025506 (2014) ; also see Schopf, D., Euchner, H., Trebin, H.-R.: Effective potentials for simulations of the thermal conductivity of type-i semiconductor clathrate systems. Physical Review B 89 (21) , 214306 (2014) ( “Pailhès, Schopf” in FIG. 10) ) . However, the κ prediction remains around one-half to one-third of κ obtained with Green-Kubo formalism in conjugate with first-principle MD simulations (see, e.g., Tse, J.S., English, N.J., Yin, K., Iitaka, T.: Thermal conductivity of solids 26 from first-principles molecular dynamics calculations. The Journal of Physical Chemistry C 122 (20) , 10682–10690 (2018) ( “Tse” in FIG. 10) ) . Such a difference in factors can be attributed to both differences in system size and simulation length. It is important to highlight that sufficient convergence on both size and time scales remain critical to attain statistically sound κ values when deploying EMD. Owing to the MatterSim’s versatility and computational efficiency, it may be feasible to predict κ for compounds in the silicon clathrate family at the first principle level with simulations using thousands of atoms in size and nanosecond (s) in length.

[0067] Referring again to FIG. 10, EMD simulations are also conducted for the two pure silicon clathrates. While little difference is observed for κ of two different types of clathrate, a twenty-fold reduction is observed when comparing the κ from Si46  to Ba8Si46. Such a reduction can be explained by the fact that the inclusion of a metal host enhances the phonon-phonon scattering, which may in turn reduce the thermal conductivity. Overall, the remarkable performance of MatterSim in modeling κ for complex materials suggests its transferability in different classes of materials and sets its foundation in spotting next generation energy-harvesting materials. Besides predicting κ from EMD for aSi, when including quantum effects and three-phonon interactions for these materials, it is also insightful to investigate their κ at the classical limit and when incorporating high-order anharmonicities.

[0068] II. High-Throughput Screening

[0069] After the assessment of the performance of MatterSim when incorporating with three-phonon BTE and EMD, the high throughput screening of thermal conductivity of crystals preexisting in Material project database is performed utilizing BTE and EMD for specific type of systems of interest (see, e.g., Jain, A., Ong, S.P., Hautier, G., Chen, W., Richards, W.D., Dacek, S., Cholia, S., Gunter, D., Skinner, D., Ceder, G., et al.: Commentary: The materials project: A materials genome approach to accelerating materials innovation. APL materials 1 (1) (2013) ) . For this purpose, crystals are selected with the only constraint being on the number of atoms in the unit cell to 15. After predictions are performed on over 60,0000 materials, systems that exhibit imaginary frequencies at 0 K are filtered out at the first round. Imaginary frequencies at 0 K indicate the instability of systems modeled, and would require temperature re-normalization of both second and third order force constants, often by fitting a temperature dependent effective potentials or via self-consistent approach. Later, for materials exhibiting imaginary frequencies, a specific type of material pool can be selected to further conduct 3ph-BTE FC or EMD by going through the computational system as described above. After quality-checking, a set of over 10,000  predictions is established, which will be beneficial to the community of thermal transport and management.

[0070] III. Methods

[0071] A. Data Preparation and Model Training`

[0072] A hurdle of the success of AI-based emulator in predicting thermal conductivity is how to achieve extreme accuracy at near-equilibrium (close to local minimum) states. The term “emulator” represents the uMLFF. The extreme accuracy at near-equilibrium is a different requirement from previous universal materials emulators, which mostly require sparse samples from vast space. To help ensure the accuracy of the emulator to realize its predictive power, a specially designed data generation algorithm is devised for thermal conductivity prediction. The data generation scheme includes two parts: (a) the uncertainty based selection rules, as well as the molecular dynamics based data exploration technique, as disclosed in International Application No. PCT / CN2024 / 081751, and (b) to further enhance the accuracy of phonon prediction, the data set is augmented with perturbative materials structures with displaced atoms along phonon modes.

[0073] B. MD Simulations using MatterSim

[0074] Molecular Dynamic simulations were performed using the Large-scale Atomic / Molecular Massively Parallel Simulator (LAMMPS) software package (see Thompson, A.P., Aktulga, H.M., Berger, R., Bolintineanu, D.S., Brown, W.M., Crozier, P.S., Veld, P.J., Kohlmeyer, A., Moore, S.G., Nguyen, T.D., et al.: Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales. Computer Physics Communications 271, 108171 (2022) ) . In particular, the energy, force, and optionally stress of the target systems are inferred, at each step of the MD simulation, by a pre-trained MatterSim model via an in-house implementation of pair style. A time step of 1 fs (femtosecond) is used to propagate the  equation of motion. Before the production run, a 200 ps (picosecond) equilibration run is performed under the canonical (NVT) ensemble, during which the temperature is controlled by the Nosé-Hoover thermostat (see, e.g., Capinski, W., Maris, H., Bauser, E., Silier, I., Asen-Palmer, M., Ruf, T., Cardona, M., Gmelin, E.: Thermal conductivity of isotopically enriched si. Applied physics letters 71 (15) , 2109–2111 (1997) ; also see Parrinello, M., Rahman, A.: Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied physics 52 (12) , 7182–7190 (1981) ; also see Nosé, S.: A unified formulation of the constant temperature molecular dynamics methods. The Journal of chemical physics 81 (1) , 511–519 (1984) ; also see Martyna, G.J., Tobias, D.J., Klein, M.L.: Constant pressure molecular dynamics algorithms. The Journal of chemical physics 101 (5) , 4177–4189 (1994) ; also see Shinoda, W., Shiga, M., Mikami, M.: Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B 69 (13) , 134103 (2004) ) with a relaxation time of 0.1 ps. After the equilibration run, the NVT production is run for 1 ns with similar settings as the equilibration. During the production runs, configurations are sampled every 10 fs, accumulating 100000 samples for subsequent κ predictions with equilibrium molecular dynamics.

[0075] C. Thermal conductivity predictions from Equilibrium Molecular Dynamics.

[0076] To obtain thermal conductivity κ with equilibrium molecular dynamics (EMD) , the running thermal conductivity (RTC) tensor καβ is formulated at a given correlation time t:

[0077] In equation 1, Cαβ (t') represents the heat current autocorrelation function (HCACF) . The HCACF can be expressed in terms of heat flux J:

[0078] Cαβ(t′) =<J (0) J (t′) >.           (2)

[0079] Equation 2 implies an accumulative HCACF between the heat flux at its time origin and at a designated correlation time. Analogous to the fundamental Hardy formula (see, e.g., Hardy, R.J.: Energy-flux operator for a lattice. Physical Review 132 (1) , 168 (1963) ) , J can be decomposed into potential Jpot and convective Jconv components:

[0080] E and v shown in equation 3 stand for per-atom energies (the sum of per-atom potential and kinetic energy) and velocities of atom i. Both quantities are available at every step of the MD simulations, and thus Jconv can be evaluated. Although the contribution to κ from <Jconv (0) Jconv (t′) > is negligible for solids, the contribution to κfrom cross terms between Jconv and Jpot can be over 10%. It may thereby be critical to retain the Jconv component in the heat flux calculation so as to accurately predict κ. As for Jpot, its implementation is an active field of research and particularly challenging for message passing neural networks (MPNNs) . Given the semi-local nature of MPNNs, the long-range interactions can be partially and implicitly included. Recent work by Langer et al (see, e.g., Langer, M.F., Knoop, F., Carbogno, C., Scheffler, M., Rupp, M.: Heat flux for semi-local machine-learning potentials. arXiv preprint arXiv: 2303.14434 (2023) ; see also Langer, M.F., Frank, J.T., Knoop, F.: Stress and heat flux via automatic differentiation. arXiv preprint arXiv: 2305.01401 (2023) ) has proven that the Jpot implementation, without properly accounting for the minimum image convention (MIC) for direct neighbor nodes in MPNNs, leads to an underestimation of κ by 40%.

[0081] In this regard, the semi-local Jpot implementation from Langer et al is adopted, which was originally compatible with SchNet potential (see, e.g., Schütt, K., Kindermans, P.-J., Sauceda Felix, H.E., Chmiela, S., Tkatchenko, A., Müller, K.-R.: Schnet: A continuous-filter convolutional neural network for modeling quantum interactions. Advances in neural information processing systems 30 (2017) and Schütt, K.T., Sauceda, H.E., Kindermans, P.-J., Tkatchenko, A., Müller, K.-R.: Schnet–adeep learning architecture for molecules and materials. The Journal of Chemical Physics 148(24) , 241722 (2018) ) . Specifically, Jpot takes the following form:

[0082] In equation (4) , I, J, K are the atom indices, ν and ε symbolize the vertices and edges for the graph representation of the target systems, and indicates the pairwise distances between atom I and J under MIC. Despite yielding accurate Jpot when incorporating into MPNN, the above implementation requires explicit evaluations of the Jacobian and its counterpart which results in a quadratic scaling. Therein, to avoid the insurmountable computational cost, Jpot is procured in a post-process manner, so that the computation of heat flux for the sample configurations can be equally partitioned and spontaneously parallel among hundreds of GPU cards. With Jpot and Jconv being defined, κ can be obtained by performing the time integral over HCACF. However, as proven by Grasselli et al. (see, e.g., Grasselli, F., Baroni, S.: Invariance principles in the theory and computation of transport coefficients. The European Physical Journal B 94, 1–14 (2021) ) the transport coefficients derived from the Helfand-Einstein (HE) integral (see Einstein, A.: Annalen der physik. Nr 10, 891 (1905) ) share a similar asymptotic limit as the Greek-Kubo (GK) integral, yet preserve much  smoother convergence. Therefore, the GK integral shown in equation 1 is reformulated to a HE integral:

[0083] To obtain καβ following equations 1 ～ 5, three separate MD simulations are performed. The average value of κ and its corresponding standard error from these runs are applied to access the κ predictions for MatterSim via equilibrium molecular dynamics.

[0084] D. κ Predictions from Boltzmann Transport Equation

[0085] To account for nuclear quantum effects while predicting κ, MatterSim is integrated with the Boltzmann Transport Equation (BTE) . In this work, the BTE calculations are performed using phono3py (see, e.g., Togo, A., Chaput, L., Tadano, T., Tanaka, I.: Implementation strategies in phonopy and phono3py. J. Phys. Condens. Matter 35 (35) , 353001 (2023) https:  / / doi. org / 10.1088 / 1361-648X / acd831) . Interested readers can refer to (Togo, A.: First-principles phonon calculations with phonopy and phono3py. Journal of the Physical Society of Japan 92 (1) , 012001 (2023) and Togo, A., Chaput, L., Tadano, T., Tanaka, I.: Implementation strategies in phonopy and phono3py. Journal of Physics: Condensed Matter (2023) ) for full details and implementations of lattice dynamics and Boltzman Transport Equation. Here the fundamentals are summarized, and the integration between phono3py and MatterSim is highlighted. When computing κ with BTE, its core inputs, 2nd order (harmonic) and 3rd order (anharmonic) force constants can be obtained by solving linear system of equations between the displaced structures and their forces:

[0086] In Equation 6 and 7,  represents the force for atom k in lattice l along direction α, and φlkα, l′k′α′and φlkα, l′k′α′, l″k″α″stand for the 2nd and 3rd order force constant matrices.  describes the super-cell displacements relative to equilibrium structures. Substituting the expensive first-principle calculations, MatterSim is deployed to infer forces for displaced super-cell structures. Given the harmonic force constants, the phonon frequency can be computed by diagonalizing the dynamical matrices (mass re-scaled 2nd order force constants) , as shown in equation 8.

[0087] The species masses involved in the construction of the second order force constants, the eigenvectors, and the corresponding eigenvalues are represented by mk and mk′,  and ω2 accordingly. Solving ω along the high-symmetry Brillouin zone recovers dispersion relations compatible to inelastic neutron scattering measurements and reflects the vibrational energy propagation through crystal lattices. Translating indices from real to reciprocal space, the phonon group velocities, constant volume heat capacity and density of states are derived:

[0088] Where N means the total number of atoms, and δ showcases the Dirac delta function. From equation 9 and 11, DOS and the group velocities result from the integration and derivatives of phonon frequencies along reciprocal q-grids. Equation 10  indicates that C is composed of phonon energy  and the derivative of phonon population (nqμ) with respect to temperature. It is worth noting that obeys Bose-Einstein distributions and its equal partitioned counterpart,  can be coupled to equation 10 in modeling phonon populations at classical limit, which can be particularly insightful when comparing κ from MD simulations. Since phonons with different wavelengths and frequencies do not interact under harmonic approximation, anharmonic force constants may be necessary to resolve a finite thermal conductivity. Compositing indices at the reciprocal space qμ to ν, one reaches the 3rd order projections and expresses the imaginary component of phonon self energy in a manner analogous to the many body perturbation theory.

[0089] In equation 12, φ0κα, l′κ′α′, l″κ″α″illustrates the 3rd order force constants within the unit-cell and R0 implies the equilibrium positions. With C, υqv and Γν, the compacted lattice conductivity expression is then:

[0090] Equation 14 presents the solution of the linearized BTE using the relaxation time approximation (RTA) , where the phonon lifetime under RTA is obtained by substituting equilibrium phonon populations into equation 13. Despite recent developments in solving the BTE beyond RTA (see, e.g., Barbalinardo, G., Chen, Z., Lundgren, N.W., Donadio, D.: Efficient anharmonic lattice dynamics calculations of thermal transport in crystalline and disordered solids. Journal of Applied Physics 128(13) (2020) ; see also Omini, M., Sparavigna, A.: Beyond the isotropic-model approximation in the theory of thermal conductivity. Physical Review B 53 (14) , 9064 (1996) ; see also Chaput, L., Togo, A., Tanaka, I., Hug, G.: Phonon-phonon interactions in transition metals. Physical Review B 84 (9) , 094302 (2011) ) , these methods are more computationally demanding, making them impractical for investigating thermal conductivity across a wide range of compositions, unit cell sizes, and crystal structures. Therefore, most recent high-throughput screening studies of κ for various materials may employ RTA. In this study, κ is also predicted under RTA, along with force constants computed using a uniform displacement of

[0091] The ongoing search for high thermal conductivity materials is driven by the need to meet the ever-increasing demands of various industries, from electronics to energy. As the boundaries of material performance are continually pushed, a deeper understanding of thermal conductivity will undoubtedly play a pivotal role in shaping the future of technological advancement. The theoretical study of thermal conductivity involves a delicate balance between the accuracy of the models and the computational feasibility. Despite these challenges, advancements in computational techniques and the growing power of supercomputers offer promising avenues for further progress in this field. As the theoretical understanding of thermal conductivity continues to be  refined, the above work further approaches the rational design of materials with tailored thermal properties, opening up new horizons in fields ranging from energy to electronics. Further, the examples disclosed herein can be extended to other material properties than thermal conductivity, such as electronic structure.

[0092] In some embodiments, the methods and processes described herein may be tied to a computing system of one or more computing devices. In particular, such methods and processes may be implemented as a computer-application program or service, an application-programming interface (API) , a library, and / or other computer-program product.

[0093] FIG. 11 schematically shows a non-limiting embodiment of a computing system 1100 that can enact one or more of the methods and processes described above. Computing system 1100 is shown in simplified form. Computing system 1100 may embody the computing systems 100, 400, 402, client devices 102, 404, and materials databases 104, 406 described above. Components of computing system 1100 may be included in one or more personal computers, server computers, tablet computers, home-entertainment computers, network computing devices, video game devices, mobile computing devices, mobile communication devices (e.g., smartphone) , and / or other computing devices, and wearable computing devices such as smart wristwatches and head mounted augmented reality devices.

[0094] Computing system 1100 includes a logic processor 1102, volatile memory 1104, and a non-volatile storage device 1106. Computing system 1100 may optionally include a display subsystem 108, input subsystem 1110, communication subsystem 1112, and / or other components not shown in FIG. 11.

[0095] Logic processor 1102 includes one or more physical devices configured to execute instructions. For example, the logic processor may be configured to execute  instructions that are part of one or more applications, programs, routines, libraries, objects, components, data structures, or other logical constructs. Such instructions may be implemented to perform a task, implement a data type, transform the state of one or more components, achieve a technical effect, or otherwise arrive at a desired result.

[0096] The logic processor may include one or more physical processors configured to execute software instructions. Additionally or alternatively, the logic processor may include one or more hardware logic circuits or firmware devices configured to execute hardware-implemented logic or firmware instructions. Processors of the logic processor 1102 may be single-core or multi-core, and the instructions executed thereon may be configured for sequential, parallel, and / or distributed processing. Individual components of the logic processor optionally may be distributed among two or more separate devices, which may be remotely located and / or configured for coordinated processing. Aspects of the logic processor may be virtualized and executed by remotely accessible, networked computing devices configured in a cloud-computing configuration. In such a case, these virtualized aspects are run on different physical logic processors of various different machines, it will be understood.

[0097] Non-volatile storage device 1106 includes one or more physical devices configured to hold instructions executable by the logic processors to implement the methods and processes described herein. When such methods and processes are implemented, the state of non-volatile storage device 1106 may be transformed-e.g., to hold different data.

[0098] Non-volatile storage device 1106 may include physical devices that are removable and / or built in. Non-volatile storage device 1106 may include optical memory, semiconductor memory, and / or magnetic memory, or other mass storage  device technology. Non-volatile storage device 1106 may include nonvolatile, dynamic, static, read / write, read-only, sequential-access, location-addressable, file-addressable, and / or content-addressable devices. It will be appreciated that non-volatile storage device 1106 is configured to hold instructions even when power is cut to the non-volatile storage device 1106.

[0099] Volatile memory 1104 may include physical devices that include random access memory. Volatile memory 1104 is typically utilized by logic processor 1102 to temporarily store information during processing of software instructions. It will be appreciated that volatile memory 1104 typically does not continue to store instructions when power is cut to the volatile memory 1104.

[0100] Aspects of logic processor 1102, volatile memory 1104, and non-volatile storage device 1106 may be integrated together into one or more hardware-logic components. Such hardware-logic components may include field-programmable gate arrays (FPGAs) , program-and application-specific integrated circuits (PASIC  / ASICs) , program-and application-specific standard products (PSSP  / ASSPs) , system-on-a-chip (SOC) , and complex programmable logic devices (CPLDs) , for example.

[0101] The terms “module, ” “program, ” and “engine” may be used to describe an aspect of computing system 1100 typically implemented in software by a processor to perform a particular function using portions of volatile memory, which function involves transformative processing that specially configures the processor to perform the function. Thus, a module, program, or engine may be instantiated via logic processor 1102 executing instructions held by non-volatile storage device 1106, using portions of volatile memory 1104. It will be understood that different modules, programs, and / or engines may be instantiated from the same application, service, code block, object, library, routine, API, function, etc. Likewise, the same module, program,  and / or engine may be instantiated by different applications, services, code blocks, objects, routines, APIs, functions, etc. The terms “module, ” “program, ” and “engine” may encompass individual or groups of executable files, data files, libraries, drivers, scripts, database records, etc.

[0102] When included, display subsystem 1108 may be used to present a visual representation of data held by non-volatile storage device 1106. The visual representation may take the form of a graphical user interface (GUI) . As the herein described methods and processes change the data held by the non-volatile storage device, and thus transform the state of the non-volatile storage device, the state of display subsystem 1108 may likewise be transformed to visually represent changes in the underlying data. Display subsystem 1108 may include one or more display devices utilizing virtually any type of technology. Such display devices may be combined with logic processor 1102, volatile memory 1104, and / or non-volatile storage device 1106 in a shared enclosure, or such display devices may be peripheral display devices.

[0103] When included, input subsystem 1110 may comprise or interface with one or more user-input devices such as a keyboard, mouse, touch screen, camera, or microphone.

[0104] When included, communication subsystem 1112 may be configured to communicatively couple various computing devices described herein with each other, and with other devices. Communication subsystem 1112 may include wired and / or wireless communication devices compatible with one or more different communication protocols. As non-limiting examples, the communication subsystem may be configured for communication via a wired or wireless local-or wide-area network, broadband cellular network, etc. In some embodiments, the communication subsystem may allow  computing system 1100 to send and / or receive messages to and / or from other devices via a network such as the Internet.

[0105] PLACEHOLDER FOR FOREIGN CLAIM SUPPORT SECTION

[0106] “And / or” as used herein is defined as the inclusive or ∨, as specified by the following truth table:

[0107] It will be understood that the configurations and / or approaches described herein are exemplary in nature, and that these specific embodiments or examples are not to be considered in a limiting sense, because numerous variations are possible. The specific routines or methods described herein may represent one or more of any number of processing strategies. As such, various acts illustrated and / or described may be performed in the sequence illustrated and / or described, in other sequences, in parallel, or omitted. Likewise, the order of the above-described processes may be changed.

[0108] The subject matter of the present disclosure includes all novel and non-obvious combinations and sub-combinations of the various processes, systems and configurations, and other features, functions, acts, and / or properties disclosed herein, as well as any and all equivalents thereof.

Claims

1.A computing system, comprising:a logic subsystem; anda storage subsystem comprising instructions executable toobtain one or more metrics computed based upon an energy and a force determined for a material using a universal machine learning force field (uMLFF) ; andbased at least upon the one or more metrics, determine to use one of a first computational method or a second computational method to compute a predicted physical property value for the material.2.The computing system of claim 1, wherein the predicted physical property is a predicted thermal conductivity, wherein the first computational method is a three phonon (3ph) Boltzmann Transport Equation (BTE) with Bose-Einstein Distribution (3ph-BTE) , and wherein the second computational method is equilibrium molecular dynamics (EMD) .3.The computing system of claim 2, wherein the one or more metrics comprises a plurality of metrics, the plurality of metrics comprising phonon modes, anharmonicity score, Debye Temperature, and 3ph scattering phase space volume (V_3ph) .4.The computing system of claim 3 wherein, when the phonon modes comprise no imaginary frequencies, the anharmonicity score meets a threshold anharmonicity score, and when the V_3ph meets a threshold phase space proportion, then the instructions are executable to use the 3ph-BTE to compute the predicted thermal conductivity.5.The computing system of claim 3, wherein, when the phonon modes comprise no imaginary frequencies, the anharmonicity score meets a threshold anharmonicity score, and when the V_3ph does not meet a threshold phase space proportion, then the instructions are executable to use EMD to compute the predicted thermal conductivity.6.The computing system of claim 3, wherein, when the phonon modes comprise no imaginary frequencies and the anharmonicity score does not meet a threshold anharmonicity score, then the instructions are executable to use EMD to compute a predicted thermal conductivity.7.The computing system of claim 3, wherein the instructions executable to determine to use one of the first computational method or the second computational method further comprise instructions executable to determine whether to use a third computational method to compute the predicted physical property value, and wherein the third computational method is a three phonon (3ph) Boltzmann Transport Equation (BTE) with temperature-dependent force constants (3ph-BTE-FC) .8.The computing system of claim 7, wherein,when the phonon modes comprise imaginary frequencies and the quantum effect meets a threshold Debye temperature, then the instructions are executable to use 3-ph BTE-FC to compute a predicted thermal conductivity; andwhen the phonon modes comprise imaginary frequencies and the quantum effect does not meet a threshold Debye temperature, then the instructions are executable to use EMD to compute a predicted thermal conductivity.9.The computing system of claim 1, wherein the metrics are further computed based upon a stress of the material.10.The computing system of claim 1,wherein the physical property is an electronic structure,wherein the one or more metrics comprises a plurality of metrics, the plurality of metrics comprising phonon modes and anharmonicity score; andwherein the first computational method is a displacement-based method and the second computation method is EMD.11.A method of computing a predicted physical property value for a material, the method comprising:obtaining one or more metrics computed based upon an energy and a force determined for the material using a universal machine learning force field (uMLFF) ; andbased at least upon the one or more metrics, determining to use one of a first computational method or a second computational method to compute the predicted physical property value for the material.12.The method of claim 11, wherein the physical property is a thermal conductivity, wherein the first computational method is a three phonon (3ph) Boltzmann Transport Equation with Bose-Einstein Distribution (3ph-BTE) , and wherein the second computational method is equilibrium molecular dynamics (EMD) .13.The method of claim 12, wherein the one or more metrics comprises a plurality of metrics, the plurality of metrics comprising phonon modes, anharmonicity score, Debye Temperature, and 3ph scattering phase space volume (V_3ph) .14.The method of claim 13, wherein, when the phonon modes comprise no imaginary frequencies, the anharmonicity score meets a threshold anharmonicity score, and when the V_3ph meets a threshold phase space proportion, then using the 3ph-BTE to compute a predicted thermal conductivity.15.The method of claim 13, wherein, when the phonon modes comprise no imaginary frequencies, the anharmonicity score meets a threshold anharmonicity score, and when the V_3ph does not meet a threshold phase space proportion, then using EMD to compute a predicted thermal conductivity.16.The method of claim 13, wherein, when the phonon modes comprise no imaginary frequencies and the anharmonicity score does not meet a threshold anharmonicity score, then using EMD to compute a predicted thermal conductivity.17.The method of claim 13, wherein determining to use one of the first computational method or the second computational method to compute the predicted physical property value for the material further comprises determining whether to use a third computational method to compute the predicted physical property value, and wherein the third computational method is a three phonon (3ph) Boltzmann Transport Equation (BTE) with temperature-dependent force constants (3ph-BTE-FC) .18.The method of claim 17, wherein,when the phonon modes comprise imaginary frequencies and the quantum effect meets a threshold Debye temperature, then using 3-ph BTE-FC to compute a predicted thermal conductivity; andwhen the phonon modes comprise imaginary frequencies and the quantum effect does not meet a threshold Debye temperature, then using EMD to compute a predicted thermal conductivity.19.The method of claim 11, wherein the metrics are further computed based upon a stress of the material.20.A method of computing a predicted physical property value for a material, the method comprising:inputting structural information for a substance into a universal machine learning force field (uMLFF) ;receiving an output from the uMLFF comprising one or more of an energy, a stress, or a force;computing one or more metrics based upon the output from the uMLFF;based at least upon the one or more metrics, determining to use one of a first computational method or a second computational method to compute a predicted thermal conductivity for the material.

Images