{'fraction_non_alphanumeric': 0.029730854161999554, 'fraction_numerical': 0.016097545364354978, 'mean_word_length': 5.244604316546763, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accelerating first-principles estimation of thermal conductivity by machine-learning interatomic potentials: A MTP/ShengBTE solution

Accurate evaluation of the thermal conductivity of a material can be a challenging task from both experimental and theoretical points of view. In particular for the nanostructured materials, the experimental measurement of thermal conductivity is associated with diverse sources of uncertainty. As a viable alternative to experiment, the combination of density functional theory (DFT) simulations and the solution of Boltzmann transport equation is currently considered as the most trusted approach to examine thermal conductivity. The main bottleneck of the aforementioned method is to acquire the anharmonic interatomic force constants using the computationally demanding DFT calculations. In this work we propose a substantially accelerated approach for the evaluation of anharmonic interatomic force constants via employing machine-learning interatomic potentials (MLIPs) trained over short ab-initio molecular dynamics trajectories. The remarkable accuracy of the proposed accelerated method is confirmed by comparing the estimated thermal conductivities of several bulk and two-dimensional materials with those computed by the full-DFT approach. The MLIP-based method proposed in this study can be employed as a standard tool, which would substantially accelerate and facilitate the estimation of lattice thermal conductivity in comparison with the commonly used full-DFT solution.


Introduction
Thermal conductivity is an essential property of a material, which can play a pivotal role in the engineering design of a product. Nowadays, appropriate thermal management is a common challenge for rapidly-growing fields, like nanoelectronics and electric vehicles. In principle, for the majority of applications, materials with higher thermal conductivities are more desirable in order to facilitate the excessive heat dissipation and avoid overheating issues. For some specific applications like thermal insulating and thermoelectric materials, materials with lower thermal conductivities are nonetheless more favorable to reduce the thermal energy losses and improve thermoelectric figure of merit, respectively. In the engineering design of a product, to accurately examine the temperature evolution during the operation, the thermal conductivity of each building block ought to be known. Therefore, accurate estimation of the thermal conductivity is required to achieve the efficient design and avoid dangerous scenarios stemming from overheating. During the last decade, two-dimensional (2D) materials have gained remarkable attention as the novel class of materials to improve the design and efficiency of diverse advanced systems. Graphene 1,2 , a prototypical member of the 2D materials class, is known to show ultrahigh thermal conductivity 3,4 outperforming all known materials. The exceptionally high thermal conductivity offers the graphene as an excellent candidate to enhance the thermal management in diverse systems, like nanoelectronics and Li-ion batteries [5][6][7][8] . Despite highquality and symmetrical atomic lattice of graphene and after extensive studies accomplished during the last decade, the exact value of its thermal conductivity still remains debated from both experimental and theoretical points of view 9,10 . For the single-layer graphene, the experimentally measured thermal conductivities mostly fall within 1500-5300 W/mK 3,[11][12][13] . In fact, there exist numerous sources of uncertainty in the experimental measurements of the thermal conductivity of 2D materials, usually resulting in the remarkable scattering of reported values. Taking into account the clean and high quality of graphene crystals, it is clear that for more complicated 2D lattices, the uncertainties in the experimentally measured values can be more pronounced. As an alternative to experimental characterizations, the development of an accurate modeling approach is critical for the evaluation of the thermal conductivity of 2D materials. To the best of our knowledge, the combination of density functional theory (DFT) simulations and the solution of Boltzmann transport equation (BTE) is currently the most trusted approach to evaluate the thermal conductivity of bulk and 2D materials. The BTE solution of thermal conductivity normally requires the evaluation of second-and third-order (anharmonic) interatomic force constants, which are commonly acquired by DFT calculations over supercell lattices. The most computationally demanding step in this approach is to obtain the anharmonic force constants, which depending on the cutoff distance and lattice symmetry may require a few hundred/thousand single-point DFT force calculations. In this study, we propose the employment of machine-learning interatomic potentials (MLIP) to substitute the computationally demanding DFT calculations in obtaining the anharmonic force constants. Due to negligible computational costs of force constant calculations using MLIP, the proposed approach therefore offers a substantial acceleration in the evaluation of anharmonic interatomic force constants in comparison with the DFT-based solution. We used moment tensor potentials (MTPs) 14 , as an accurate and computationally efficient model of MLIPs [15][16][17] in the evaluation of interatomic forces. To the best of our knowledge, the ShengBTE 18 package is currently the most employed package for the full iterative solution of the BTE. In this paper, we present the integrated MTP/ShengBTE method for the practical calculation of the lattice thermal conductivity. The proposed approach was tested on several bulk and 2D materials and shows remarkably close agreement with the results available in the literature on the basis of full-DFT calculations.

Computational methods
In this work, first-principles DFT calculations were conducted to obtain the phonon dispersion relations, second-order interatomic force constants, and also to create the training sets for the MTPs. To this aim, Vienna Ab-initio Simulation Package (VASP) [19][20][21] and generalized gradient approximation (GGA) with either Perdew−Burke−Ernzerhof (PBE) 22 , Perdew-Burke-Ernzerhof revised for solids (PBEsol) 23 or revised Perdew-Burke-Ernzerhof (revPBE) 24 were adopted in the calculations. The geometry optimized structures were acquired using the conjugate gradient method with the convergence criteria of 10 -5 eV and 0.001 eV/Å for the energy and forces, respectively. The PHONOPY code 25 was employed to create the optimal sets of atomic positions for DFT or MTP force calculations and subsequently obtain phonon dispersion relations and second-order interatomic force constants with DFT and MTP-based results as inputs. The interatomic force constants were evaluated by considering the supercell structures. The details of MTP and PHONOPY interface can be found in our latest study 26 . Abinitio molecular dynamics (AIMD) simulations were performed with a time step of 1 fs. We employed Monkhorst-Pack 27 k-point grids of 3×3×1 and 2×2×2, for 2D and bulk lattices, respectively. The supercell sizes for the phonon and AIMD calculations and corresponding plane-wave cutoff energies will be explicitly mentioned for every considered example. MTP is a local potential in the sense that the total energy E of an atomistic sample with N atoms is the summation of contributions V of neighborhoods of each -th atom: . The neighborhood of a central atom is defined as a collection, via: here the -th atom is referred as a neighbor of the -th (central) atom within the preset cutoff radius cut , and are the types of the central and neighboring atoms, respectively, is the corresponding interatomic vector and neigh is the number of atoms in the neighborhood.
The contribution of each central atom and its associated neighborhood to the system's total energy shows: here, the first factor (| |, , ) is the radial part of potential MLIP, which only depends on the distance between atoms and and their atomic types and "⊗" is the outer product. The radial part is expanded by a set of radial basis functions (| |) multiplied by the ( cut− | |) 2 smoothing factor. are the corresponding values predicted by the MTP, and we, wf and ws are the optimization problem importance weights for the energies, forces and stresses, respectively, which are set to 1, 0.1 and 0.001, respectively. In this work, interatomic force constants are acquired using the MTPs 14 , trained over AIMD runs at different temperatures of 50, 300, 500 and 700 K, each with less than 1000 time steps. The training sets were sampled from the full trajectories (one configuration was sampled every five time steps). Finally, the lattice thermal conductivity is estimated using the full iterative solutions of the Boltzmann transport equation, as implemented in the ShengBTE 18 package. Harmonic force constants are obtained using the PHONOPY code 25 with DFT results as inputs. Anharmonic interatomic force constants are calculated using the trained MTPs. Since force evaluations with MTP have a negligible cost, we take the same (large) supercell samples for the harmonic and anharmonic force constant calculations. Isotope scattering is considered in all examples, which is essential for direct comparison with experimental data for naturally occurring samples. Born effective charges and dielectric constants contributions in the dynamical matrix are taken into account only for the bulk InAs structure. The convergence of the lattice thermal conductivity with respect to the q-mesh grid is also examined. To facilitate the practical application, in the data availability section, the full computational details with developed python scripts are included. There we provide the full set of the input files for all considered examples along with the numerical procedure to extract the anharmonic force constants using the trained MTPs.

Results and discussions
The main objective of the present study is to examine the accuracy of estimated lattice thermal conductivities with the MTP-based solution in obtaining anharmonic interatomic force constants. In our latest study 26 , obtained results for diverse 2D lattices confirm that the MTPbased solution can accurately reproduce the phononic properties in close agreement with density functional perturbation theory (DFPT) simulations. For the validation of the proposed approach, we consider several 2D and bulk lattices and compare MTP-based results with available data in the literature. We first examine the lattice thermal conductivity of graphene, which has been extensively studied experimentally and theoretically during the last decade. By having a glance at the available data in the literature, it becomes conspicuous that despite a rather simple bonding nature, highly symmetrical lattice and high quality of graphene crystals, the exact value of its thermal conductivity still remains debated from both experimental and theoretical points of view. According to full-DFT BTE based studies 28,29 , it has been suggested that the type of exchange-correlation functional can yield substantial effects on the estimated thermal conductivity. In the work by Qin et al. 28 , single-layer graphene's thermal conductivity at room temperature was predicted from 1936 W/mK (with PBE) to 4376 W/mK (with VdW-DF2). In the recent work by Taheri and co-workers 29 , they predicted distinctly higher thermal conductivities ranging from 5442 W/mK (with LDA) to 8677 W/mK (with PBEsol). To check the effects of exchange-correlation functional for the case of graphene, we thus consider PBE, PBEsol, and revPBE functionals. In our calculations, we consider a plane-wave cutoff energy of 500 eV and the lattice constant of graphene is found to be 2.467, 2.461, and 2.479 Å, respectively with PBE, PBEsol, and revPBE functionals. In our interatomic force constants calculations, we consider a 10×10×1 supercell, which includes 200 carbon atoms. The AIMD calculations for the preparation of training sets are conducted over 6×6×1 supercells. As discussed in our recent study 26 , incorporation of larger supercells in the AIMD calculations can however secure higher accuracy. MTPs with 901 parameters are passively trained and subsequently used to reproduce the phonon dispersion relations of graphene. In Fig. 1 the phonon dispersion relations of single-layer graphene by MTP and DFPT results are compared, which reveal close agreement for the all three considered functional. We next examine the lattice thermal conductivity using the trained MTPs to acquire anharmonic interatomic force constants, in which the interactions with elevenths nearest neighbors are considered. The thermal conductivity of graphene at the room temperature is predicted to be 3730, 3640, and 3600 W/mK, respectively, on the basis of trained MTPs with PBE, PBEsol, and revPBE functionals. It is noticeable that unlike the earlier reports 28,29 , the effect of exchange-correlation functional on the estimated thermal conductivity value is negligible. In the work by Fugallo et al. 30 they calculated the graphene's thermal conductivity to be 3600 W/mK, which is very close to our estimated values. Using the PBE functional and ShengBTE package for the BTE solution, the room temperature thermal conductivity of naturally occurring graphene has been reported to be 1936 28 , 3100 31 , 3550 32 , 3845 33 , 3590 34 , 3720 35 and 3288 36 , which are generally consistent with our results. One of the main sources of scatterings in different reports is related to the supercell size effect. We found that for the case of graphene the thermal conductivity converges for the supercell size of 8×8×1. In Fig. 1d the MTP-based results for the temperature-dependent thermal conductivity of single-layer graphene are compared with experimental and full-DFT based estimations, which show good agreement. Furthermore, using the MTP-based solution, the temperature power factor for the thermal conductivity of graphene is found to be 1.35, compatible with full-DFT based reports of 1.32 by Lindsay et al. 36 and 1.34 by Fugallo et al. 30 . Fig. 1e compares the cumulative lattice thermal conductivity of single-layer graphene at the room temperature as a function of mean free path by MTP and full-DFT calculations, which also reveal close trends. For the single-layer graphene, acoustic phonons are known as the dominant heat carriers. In Fig. 1f the contribution of ZA, TA and LA acoustic modes on the total lattice thermal conductivity of graphene by MTP and full-DFT calculations are compared, which also reveal close agreement. It is noticeable once again that MTP-based solutions show negligible effects for different functionals. Despite the remarkable scattering between different experimental and theoretical results for the case of graphene, the conducted comparison clearly reveals the high accuracy of the proposed MTP-based solution. Because of the remarkable scattering in the experimentally measured and theoretically predicted thermal conductivity of graphene and other 2D materials, we next examine the accuracy of the proposed approach for the bulk materials. In this case, we consider silicon and InAs, the examples also considered in the original ShengBTE manuscript by Li et al. 18 . We also consider bulk diamond and BAs for further examination. For the silicon, InAs, diamond and BAs plane-wave cutoff energies of 330, 300, 500, and 400 eV are employed, respectively. The lattice constants of silicon and InAs with PBE functional are found to be 5.47 and 6.06 Å, respectively, match closely with those reported by Li et al. 18 . For the diamond and BAs, we use the PBEsol functional because the obtained lattice constant of 3.572 and 4.779 Å, respectively, are in close agreement with corresponding experimentally measured values of 3.567 and 4.777 43 Å. We remind that all considered bulk structures show Zinc Blende diamond-like lattices. In line with the work by Li et al. 18 , harmonic interatomic force constants are calculated using the 5×5×5 supercell samples. AIMD calculations are conducted over rectangular supercells with 144 atoms to create the training sets. MTPs with 901 and 1009 parameters are trained for monoelemental and binary lattices, respectively. In Fig. 2, the predicted phonon dispersion relations for considered bulk lattices by passively fitted MTPs are compared with DFPT-based results, which reveal remarkably close agreement. In Fig. 3, the estimated phononic thermal conductivity of the considered bulk structures by the accelerated method are compared with experimental and full-DFT counterparts. The negligible computational costs of the force evaluations with MTP enables using relatively large supercells. Our calculations are done for 5×5×5 supercell structures for the anharmonic force constants calculations, whereas in work by Li et al. 18 44 , they predicted a room temperature thermal conductivity of 2300 W/mK for naturally occurring BAs, which is only by around 9% difference with our predicted value. Currently, 2D materials are one of the fastest-growing class of materials. The complexity and low-symmetry of novel 2D materials lattices demand for excessive computational costs for the evaluation of anharmonic force constants. In fact, the introduced alternative by this study is highly promising for the assessment of 2D materials thermal conductivity. From the modeling points of view, to deal with computational difficulties in the valuation of anharmonic force constants using the conventional DFT calculations, researchers may consider smaller supercells, interactions with fewer neighbors, coarser K-point grids and/or lower plane-wave cutoff energies. These simplifications can justify the remarkable scattering in the predicted values for the 2D materials' thermal conductivity. Next we examine the results for the thermal conductivity of several 2D materials by the MTP-based accelerated approach, taking into account that the exact values of thermal conductivities are debated. Here we consider singlelayer penta-graphene, silicene, phosphorene, F-diamane, and MoS2, as illustrated in Fig. 4. Plane-wave cutoff energies of 500, 330, 400, 500, and 400 eV are considered for pentagraphene, silicene, phosphorene, F-diamane, and MoS2, respectively. All the considered 2D lattices are isotropic, except phosphorene (find Fig. 4c), which shows anisotropic lattice constants of 4.624 and 3.299 Å along the armchair and zigzag directions, respectively. The lattice constants of penta-graphene, silicene, F-diamane, and MoS2 are found to be 3.642, 3.868, 2.546, and 3.184 Å, respectively, all consistent with the reported values in the literature. Worthy to note that F-diamane was most recently experimentally realized by Bakharev and coworkers 54 . For the preparation of training sets, AIMD simulations are conducted using 4×4×1, 6×6×1, 3×5×1, 4×4×1, and 5×5×1 supercells for penta-graphene, silicene, phosphorene, F-diamane, and MoS2, respectively. Similarly to the bulk structures, MTPs with 901 and 1009 parameters are trained for monoelemental and binary 2D lattices, respectively. For the evaluation of interatomic force constants, 4×4×1, 10×10×1, 4×6×1, 5×5×1 and 6×6×1 supercells are employed for penta-graphene, silicene, phosphorene, F-diamane, and MoS2 monolayers, respectively. The predicted phonon dispersion relations for the considered 2D lattices by the passively fitted MTPs are compared with those by the DFPT method in Fig. 5, which reveal excellent agreement and highlight the high accuracy of MTPs in describing the interatomic forces. We next examine the lattice thermal conductivity of the considered 2D lattices. Anharmonic interatomic force constants are acquired using the trained MTPs, in which the interactions with the twelfth, eleventh, twelfth, ninth, and eights nearest neighbors are considered for pentagraphene, silicene, phosphorene, F-diamane, and MoS2 monolayers, respectively. In Fig. 6, predicted lattice thermal conductivities of the considered 2D lattices by the MTP-based BTE solution are compared with full-DFT counterparts. Similar to the case of graphene, a remarkable scattering between different DFT-based BTE solutions causes ambiguity in finding the exact values of lattice thermal conductivity. With decreasing of the lattice symmetry, the complexity of the DFT-based solution increases. Therefore to decrease the computational costs, various simplifications are applied in different studies, resulting in scatting of reported values. Let's consider the MoS2, for which the in-plane lattice thermal conductivity of bulk structure is experimentally measured to be 85-110 W/mK 55 . For 2D materials, because of absence of van der Waals phonon scattering from the interactions with adjacent layers 4 , the thermal conductivity of monolayers is usually expected to be higher than those of bulk and few-layer structures. This way, the thermal conductivity of single-layer MoS2 is expected to be higher reported experimental range of 85-110 W/mK 55 for the bulk structure, as also theoretically confirmed by Gandi et al. 56 and Gu et al. 57 . However, due to computational setups, DFT-based reports may not satisfy this expectation. For the single-layer MoS2, our proposed approach yields a thermal conductivity of 152 W/mK, very close to the predicted value of ~155 W/mK by Peng et al. 58 and 138 W/mK by Gu et al. 57 . In the study by Gu et al. 57 , they reported the value for the sample with a length of 10 μm, which can be probably lower than the converged diffusive thermal conductivity. Comparisons illustrated in Fig. 6 clearly show that the estimated thermal conductivities for the considered monolayers using the MTPbased accelerated method are within the values acquired using the full-DFT based solutions, which confirm the remarkable accuracy of the proposed approach. Conducted comparisons between the MTP-based BTE solutions and full-DFT counterparts for considered bulk and 2D lattices confirm the remarkable accuracy of the accelerated approach. This is a highly promising finding, as the computational costs are substantially decreased, and at the same time, the sacrificial in accuracy remains marginal. In fact, since the MTP-based approach enables the consideration of larger supercells and cutoff distances in the evaluation of anharmonic interatomic force constants, they can yield more accurate results than simplified DFT-based solutions. By decreasing the symmetry in the atomic lattice and increasing the number of atoms in the primitive unit cell, the number of structures for the force constant calculations increases substantially. These calculations can get exceedingly expensive using the DFT-based approach, whereas with the MTP-based solution, the computational costs always stay negligible. In comparison with the DFT-based solution, the proposed approach is straightforward and does not require the convergence tests with respect to the K-point grid and plane-wave cutoff energy. In the MTP-based approach, independent of lattice complexity, relatively short AIMD trajectories of less than 4000 time steps are required for the training of accurate interatomic potentials. In our earlier study 26 , our extensive results for diverse and complex 2D lattices confirm that for the majority of cases the MTPs trained over only 1000 time steps of AIMD simulations at 50 K can very accurately reproduce the phonon dispersion relations and other phononic properties in comparison with DFPT results. However, to improve the accuracy for the evaluation of anharmonic interatomic force constants the inclusion of AIMD trajectories at high temperatures is highly recommended. We remind that even for the case of a highly symmetrical lattice of graphene, the proposed approach is highly efficient from the computational point of view. Nonetheless, the computational efficiency of MTP-based solution will be even more obvious when moving from high-to low-symmetry structures and in particular for 2D lattices that require the consideration of interactions with distant neighbours. Our results in another word suggest that MTP can accurately substitute the DFT method in the evaluation of interatomic force constants. As it has been confirmed in recent studies for the case of graphene 73 and bulk BAs 74 , the inclusion of four-phonon scattering can be a critical issue to accurately estimate the thermal conductivity of particular structures. In order to take into count the four-phonon scattering in the solution of thermal conductivity, MTP can offer an outstanding application prospect owing to its inehrent negligible computational cost. We therefore hope that MTPbased approach can serve as a standard and versatile tool to conveniently and accurately examine the lattice thermal conductivity.

Conclusion
Combination of DFT simulations and BTE solution is currently the most trusted first-principles method to examine the lattice thermal conductivity. The main computational bottleneck of this method is to acquire the anharmonic interatomic force constants using the DFT calculations. In this work, we show that MTPs trained over short ab-initio molecular dynamics trajectories can substantially accelerate the evaluation of anharmonic interatomic force constants. Comparison between the MTP-based BTE solutions and full-DFT counterparts for the thermal conductivity of several bulk and 2D structures confirm the remarkable accuracy of the proposed approach. In order to facilitate the practical employment of the proposed approach in conjunction with the ShengBTE package for the BTE solution, elaborated computational details and all considered examples are included in the data availability section. MTP/ShengBTE approach is therefore expected to serve as a standard tool to conveniently, efficiently, and accurately examine the lattice thermal conductivity, which may otherwise require excessive computational resources with commonly employed full-DFT counterpart.
Here "n_cores" is the number of cores used for parallel training of MTP, "p.mtp" is the input/output (curr-pot-name/trained-pot-name) MTP file, "train.cfg" is the training set in the internal *.cfg MLIP format, the option "max-iter" determines the maximum number of iterations in the optimization algorithm. The options "energy-weight", "forceweight", and "stress-weight", respectively, define the we, wf and ws weights explained earlier.
In our work, we conducted the passive training, by parameterizing the MTPs using the subsampled AIMD trajectories. In this approach, from the complete sets of AIMD configurations, only subsamples are selected for the training of first MTPs. Nonetheless, some critical configurations that could result in the improved accuracy of trained MTPs may have been missed in the created subsamples. Therefore, the accuracy of the developed MTP "p.mtp" over current subsampled training set "train.cfg" should once again be checked over the full AIMD configurations "trainF.cfg", and the configurations with high extrapolations grades 4 will be selected, and will written to the file "trainN.cfg", via the following command: ./mlp select-add p.mtp train.cfg trainF.cfg trainN.cfg The selected configurations "trainN.cfg" should be added to the original training sets "train.cfg" and the final MTP will developed by retraining of new clean potentials over the updated training set. This way, the efficient use of conducted AIMD simulations will be guaranteed. 1.4 Structure of MTPs. MTP belongs to the family of machine-learning interatomic potentials by which potentials show flexible functional form that allows for systematically increasing the accuracy with an increase in the number of parameters and the size of the training. In the folder "Untrained-MTPs", we included three samples of clean MTPs. Depending on the number of parameters, the appropriate MTP should be chosen. Prior to training, there are some parameters to be adjusted, such as the "species_count", "min_dist" and "max_dist" which, respectively, define the number of elements in the system, minimum atomic distance and cutoff distance of the potential. Like the classical potentials, by increasing the cutoff distance more neighbors will be included in the calculations which accordingly increase the computational costs. The number of parameters in a MTP can be calculated via: species_count ⋅ radial_basis_size ⋅ radial_funcs_count + alpha_scalar_moments + Note that "radial_funcs_count" and "alpha_scalar_moments" are the fixed features of a particular MTP and only "radial_basis_size" can be manually changed to adjust the number of constants.
2. Evaluation of phononic properties using the MTPs. In our previous work 5 , we included the full details and numerous examples for the evaluation of phononic properties using the MTP and PHONOPY 6 package in a public Mendeley dataset, please refer to: http://dx.doi.org/10.17632/7ppcf7cs27.1 3. MTP/ShengBTE interface. ShengBTE 7 is a package for computing the lattice thermal conductivity on the basis of a full iterative solution to the Boltzmann transport equation. Its main inputs are sets of second-and third-order interatomic force constants and a CONTROL file for the adjustment of computational details. In this work, the calculation of anharmonic interatomic force constants is substantially accelerated by substituting DFT simulations with the MTP-based solution. For the calculation of anharmonic interatomic force constants, ShengBTE 7 provides a script, "thirdorder.py", implementing a real-space supercell approach to anharmonic IFC calculations. In this approach, according to the defined supercell size and cutoff distance, the input geometries for the force constant calculations will be generated. For compatibility with "cfg"file format, the "thirdorder_mtp.py" is developed using the original "thirdorder_vasp.py".
Moreover, we developed an additional script "fake_vasp_calcs.py", which uses the MTP-based calculated forces and artificially create the VASP output files of "vasprun.xml". This approach provides the possibility of direct comparison of forces by MTP and VASP. These developed two python scripts are included the folder "MTP_ShengBTE_py". In the folder "Examples", complete input files are included for every structure. In this case the subfolder "ShengBTE-inputs" includes the complete input files for the ShengBTE solution (namely: CONTROL, FORCE_CONSTANTS_2ND and FORCE_CONSTANTS_3RD). Using the data provided in the subfolder called "Anharmonic-MTP", the anharmonic interatomic force constants can be obtained using the trained MTPs "p.mtp". MTP/ShengBTE interface follows the same routine as that of the VASP/ShengBTE, explained in the ShengBTE documentation. To facilitate the practical usage, for every example we included a shell script, named "getFC.sh". In the aforementioned script, "supcell" and "Cutoff" are respectively, the supercell size and cutoff neighbour for the evaluation of anharmonic force constants on the basis of primitive unitcell "POSCAR". Please note that "POSCAR", "p.mtp", "mlip.ini" and related python scripts should be located in this folder for complete calculations.