{'fraction_non_alphanumeric': 0.02282396914390234, 'fraction_numerical': 0.010736013360372183, 'mean_word_length': 5.128167641325536, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, '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}}
Exploring Phononic Properties of Two-Dimensional Materials using Machine Learning Interatomic Potentials

Phononic properties are commonly studied by calculating force constants using the density functional theory (DFT) simulations. Although DFT simulations offer accurate estimations of phonon dispersion relations or thermal properties, but for low-symmetry and nanoporous structures the computational cost quickly becomes very demanding. Moreover, the computational setups may yield nonphysical imaginary frequencies in the phonon dispersion curves, impeding the assessment of phononic properties and the dynamical stability of the considered system. Here, we compute phonon dispersion relations and examine the dynamical stability of a large ensemble of novel materials and compositions. We propose a fast and convenient alternative to DFT simulations which derived from machine-learning interatomic potentials passively trained over computationally efficient ab-initio molecular dynamics trajectories. Our results for diverse two-dimensional (2D) nanomaterials confirm that the proposed computational strategy can reproduce fundamental thermal properties in close agreement with those obtained via the DFT approach. The presented method offers a stable, efficient, and convenient solution for the examination of dynamical stability and exploring the phononic properties of low-symmetry and porous 2D materials.


Introduction
Phonon dispersion relations (PDRs) are key components to study the lattice dynamics and atomic vibrations in a crystal 1,2 . They also provide useful information regarding the transport properties such as the thermal conductivity. In particular, PDRs are extensively employed to examine the dynamical stability of various compositions and structures. To acquire PDRs, the most popular theoretical approach is to calculate the force constants using the density functional theory (DFT) simulations. To assess the thermal properties, DFT calculations are commonly executed using supercell structures. DFT is a computationally efficient approach for the majority of highly symmetrical lattices. In recent years, two-dimensional (2D) materials 3,4 are gaining remarkable attentions because of their unique properties, suitable to address critical challenges in various advanced technologies like nanoelectronics, energy harvesting devices or biosensors. Since the isolation of graphene 1,2 , 2D materials family has been continuously and quickly extending. This large family of new materials includes highly symmetric lattices, such as the graphene-like structures 5,6 and 2H transition metal dichalcogenides 7,8 , as well as many structures with low-symmetry or nanoporous lattices like graphdiyne 9 and 1T' transition metal dichalcogenides 10 . To date, high-symmetry and isotropic lattices have been excessively studied. Recently, however, anisotropic and nanoporous 2D lattices are attracting considerable attention, because they offer new possibilities to achieve angle-dependent properties and design more efficient energy storage/conversion systems. In order to theoretically predict novel 2D systems, it is thus essential to carefully examine the dynamical stability and explore the vibrational properties on the basis of PDRs. For lowsymmetry and nanoporous 2D structures, currently employed DFT simulations exhibit limited flexibility due to severe computational issues. In these cases, to bypass the computational limitations, one usually treats smaller supercells and lower k-point girds and plane-wave cutoff energy within the DFT simulations. But such simplifications may actually result in poorer accuracy of the acquired PDRs. Moreover, in some cases nonphysical imaginary frequencies may appear in the PDRs thus impeding the examination of dynamical stability or thermal properties. Machine-learning based methods are currently among the most attractive solutions to address various challenges, also in materials science [11][12][13][14][15][16][17] . In this regard, machine-learning interatomic potentials (MLIP) 18 have shown outstanding efficiency in many applications of computational material science, such as predicting novel materials 19,20 lattice dynamics 21 and estimating the thermal conductivity 22,23 , etc. The main advantage of MLIPs is that they enable the efficient use of classical molecular dynamics simulations to evaluate the forces and energies with the DFT level of accuracy. The employment of MLIPs is particularly promising to study the large systems, for those DFT-based methods become infeasible due to the known computational limitations. In this work we propose a computationally efficient and accurate methodology to acquire PDRs and explore other critical phononic properties on the basis of passively fitted moment tensor potentials (MTPs) 24 . MTP belongs to the family of MLIPs 25-27 by which we understand potentials possessing flexible functional form that allows for systematically increasing of the accuracy with an increase in the number of parameters and the size of the training 28 . Similar to other MLIPs, MTP is able to approximate any model of interatomic interaction, offering a compromisable choice between accuracy and computational efficiency. The proposed approach only requires inexpensive and short ab-initio molecular dynamics trajectories, without the need in active learning or additional DFT calculations. We employed such an approach to evaluate the phononic properties of a wide variety of 2D materials, for which close agreements with density functional perturbation theory (DFPT) results are observed. We describe our method on the basis of passively-trained MTPs, which can be conveniently employed to explore the phononic properties of low-symmetry and nanoporous 2D materials, with enhanced stability and computational efficiency.

Computational methods
First-principles DFT calculations in this work were carried out using the Vienna Ab-initio Simulation Package (VASP) [29][30][31] . The generalized gradient approximation (GGA), namely the Perdew−Burke−Ernzerhof (PBE) 32 functional was employed in the calculations. We assumed the plane-wave cutoff energies of 600 eV and 400 eV for the carbon-based and carbon-free systems, respectively. For geometry optimization, the convergence tolerances for the energy and forces were set to 10 -5 eV and 0.001 eV/Å, respectively. The DFT simulations were performed over supercell samples using a 3×3×1 Monkhorst-Pack 33 k-point grid. The planewave cutoff energy of the DFPT simulations was set as the default value by VASP. The PHONOPY code 34 was utilized to create the optimal sets of atomic position for DFPT calculations and also to acquire phonon dispersions and group velocities with the DFPT results as inputs. Ab-initio molecular dynamics (AIMD) simulations were performed with a time step of 1 fs using a 3×3×1 k-point gird. A class of machine-learning interatomic potentials-moment tensor potentials 24 were used to describe interatomic interactions. Similar to classical potentials, MTPs include parameters that are optimized on a set of training configurations. In this work AIMD simulations were used to create the training sets. MTP was first proposed for single-component systems 24 and has been recently generalized to multiple component systems 19,35 . This potential is local, i.e., the total energy E of the system containing N atoms is partitioned into contributions V of neighborhoods of each -t atom: ∑ . We refer to the -th atom as a neighbor of the -th (central) atom if the distance between them is less than a predefined cut-off distance . The neighborhood is then expressed as a collection, where are the relative atomic positions (interatomic vector), and are the types of the central and neighboring atoms, respectively, and is the number of atoms in neighborhood. Each contribution to the total energy has the following form: ∑ where are the free parameters of the potential to be optimized, and are the basis functions. We construct the basis functions as all possible contractions of the moment tensor descriptors: yielding a scalar (see Ref. 24 for details). The first factor (| | ) in the aforementioned summation is the radial part depending only on the distance between atoms and and their types. We expand the radial part through a set of radial basis functions (| |) multiplied by a factor ( | |) for smoothing near the cut-off radius.
where are the radial coefficients. We denote by " " the outer product and refer to the second factor in Eq. 2 as the angular part. In order to optimize the and parameters of an MTP, one needs to solve the following minimization problem (training of MTP): where , and are the energy, atomic forces and stresses in the training set, respectively, and , and are the corresponding values calculated with the MTP, K is the number of the configurations in the training set, and w e , w f and w s are the nonnegative weights that express the importance of energies and forces and stresses in the optimization problem, respectively, which in our study were set to 1, 0.1 and 0.001, respectively. We note that the weights for the energy and force are the default values.
Here we used the PHONOPY code 34 to evaluate the phononic properties, in which MTP replaces VASP in the force calculation step 36 . To facilitate future studies, in the data availability section we present a guide to create the training set (VASP AIMD inputs), as well as the MTP training procedure and the integration code to calculate the MTP forces of the PHONOPY input structures. In particular, all the examples of this work are included. We hope that such information will facilitate the use of MTP for accurate evaluation of the phononic properties in structural systems of versatile complexity.

Results and discussions
The main objective of this study is to replace DFPT with MTP in the force-constant calculations for the evaluation of phononic properties. Although our work is focused on the 2D materials, the proposed approach is equally valid for 3D or 1D lattices. First, we study the monoelemental 2D structures. Since carbon shows the unique ability to form diverse stable 2D atomic lattices, we mostly consider carbon-based structures. In Fig. 1, the studied monoelemental 2D lattices are illustrated. Apart from the graphene (Gr), we consider pentagraphene 37 , various haeckelite 38 (Haeck) lattices, phagraphene 39 (Pha-Gr) and different graphyne 40 (GY) structures. We note that the haeckelite and phagraphene lattices include non-hexagonal carbon rings and can thus be considered as defective graphene lattices but with high degree of periodicities and close densities to that of the pristine graphene. On the other side, graphyne lattices are highly porous full carbon materials, which recently gained remarkable attentions because of experimental advances and their bright application prospect for the energy storage/conversion systems. Besides full-carbon structures, we also consider the single-layer black phosphorene (P) which has a rather complicated buckled lattice. When using the DFPT method, the phononic properties are acquired on the basis of applying small displacements in supercell lattices. Therefore, the preliminary training data sets were prepared by conducting the AIMD simulations at a temperature of 50 K, with an overall simulation time of less than 1 ps (1000 simulation time steps). In this case, the AIMD trajectories include configurations with small lattice deformations from the equilibrium condition, consistent with our employed DFPT procedure. Since the MTPs are trained with a cutoff distance of 4 or 5 Å, the AIMD simulations were conducted over supercells with periodic box sizes over 10 Å. Larger supercells can be useful to describe the local large deflections better. Nonetheless, supercells with the minimum sizes closer to 10 Å are more computationally efficient because the costs of AIMD simulations increase exponentially with the number of atoms. To facilitate future studies, all created training sets in this work are provided in the data availability section. Using the AIMD results, we trained MTPs with 901 parameters for the monoelemental systems. The cutoff distance for graphene, other carbon allotropes and phosphorene were assumed to be 3, 4 and 5 Å, respectively. In Fig. 2 we compare the PDRs predicted by MTPs with those by the DFPT. Remarkably close agreement between the MTP and DFT results are obtained. This is a highly promising outcome since the developed MTPs are first attempts without any long sequence of optimization. Noticeably for the penta-graphene, haeckelite lattices and phosphorene, acoustic and optical branches are reproduced with a high level of accuracy. The same observation is also valid for graphyne lattices. Notably, in some cases (Haeck-5, P, GY-2 and GY-3) slight imaginary frequencies appear in the PDRs calculated with DFPT, which are not observed in the MTP-based results. These imaginary branches in the DFPT results originate from computational artifacts, and are not representative of dynamical instabilities [41][42][43] and might be removed via more precise setups for the calculations, such as increasing the supercell size, increasing plane-wave cutoff energy or improving the resolution of k-point girds. These modifications can be affordably examined for highly symmetrical lattices like graphene, but for less symmetrical and nanoporous lattices may substantially increase the computational costs and lead to complexity of the problem. Despite close agreement between the predicted PDRs for the most of monoelemental structures, two exceptions exist in which the MTP-based results show some inaccuracies. Among all the studied samples, it is clear that the high-frequency optical modes are well reproduced by the MTP. Unexpectedly for the highest symmetrical lattice of graphene, optical modes with the highest frequencies are not well reproduced by the MTP (Fig. 2a). On the other side for the phagraphene (Fig. 2h), a slight soft mode appears in the MTP result which is absent in the DFPT result. These disagreements between the DFPT and MTP results reveal that some degree of extrapolation may have occurred and accordingly suggest that for these samples the created passive training sets may not have been as accurate as those for other samples. We remind that when conducting the AIMD simulations at the low temperature of 50 K for a very short simulation time, the occurring configurations remain highly correlated and some critical configurations may never get explored in the training set. To address this issue, a simple effective solution is to expand the training set by including additional uncorrelated AIMD trajectories at higher temperatures. To that end, additional 1 ps AIMD trajectories were included by conducting the simulations at 200, 300, 500, 700 and 900 K (each one for 0.2 ps). This way the total AIMD trajectories include 2 ps, and the new training sets were created accordingly. In Fig. 3, we include the predicted PDRs for singlelayer graphene and phagraphene by the improved training sets and compare the new results with the original results. It is clear that by incorporating higher temperature trajectories, the imaginary frequencies in the phagraphene's PDR vanishes, while the close agreement with DFPT results at higher frequencies is kept intact. Similarly, for the case of graphene, the accuracy of high-frequency optical modes is considerably improved without affecting the lower frequency optical and acoustic modes. The acquired results reveal that the incorporation of high-temperature trajectories can further enhance the accuracy and stability of MTPs. It is worth noting that in prior works, MLIPs have been successfully employed to accurately reproduce the phononic properties of graphene 44,45 . In comparison with these earlier works, the preparation of training sets in our approach is however more computationally efficient. In addition, when comparing with prior studies 44,45 , our method could more accurately reproduce the acoustic modes in graphene's PDR, which is an important finding because these modes are the dominant heat carriers in graphene. We next consider the binary 2D systems shown in Fig. 4. We note that various carbon-nitride 2D systems, like CN 41 , C 2 N 46 , C 3 N 4 47 and C 3 N 48 are among the most attractive 2D semiconductors that have been experimentally fabricated and exhibit promising performances for a wide range of applications. BC 3 49 has been also experimentally fabricated, C 4 N and C 7 N 6 , C 9 N 4 and C 10 N 3 have been theoretically predicted by Li et al. 50   The predicted PDRs by MTPs are compared with those by the DFPT method for the binary 2D lattices in Fig. 5. As it is clear, passively-trained MTPs can very accurately reproduce the PDRs for the considered binary lattices. It is noticeable that PDRs on the basis of MTPs are free of the slightly imaginary frequencies that have occurred around the Γ points in the DFPT results for some samples. The C 3 N 4 monolayer with a full flat structure shows conspicuous imaginary frequencies, accurately reproduced also by the MTP-based method. These results are highly promising taking into account the complexity of the considered structures, especially in view of rectangular SiP 2 and As 2 Se 3 exhibiting intricate buckled lattices. It should also be noted that while the acoustic modes are very closely reproduced for all considered samples, slight deviations are observed for some of the optical modes in a few samples, particularly for the C 2 N and BC 3 monolayers. As discussed for the case of graphene, to better reproduce the optical modes, the training set should be improved by including the AIMD trajectories at higher temperatures. We next explore the accuracy of the PDRs predictions by the MTP method for ternary 2D lattices. In Fig. 6 we compare the PDRs for six different ternary lattices by the DFPT and MTP methods. We note that BC 10 N 2 is recently predicted by Tromer et al. 54 , BC 6 N shows hexagonal and rectangular atomic lattices, BrCuTe 2 and ICuTe 2 are also predicted according to their bulk counterparts, and BC 6 N 6 is a novel ternary lattice whose stability we examine. It is noticeable that imaginary frequencies in the DFPT results are also well reproduced by the passively trained MTPs. These results, consistent with our previous observations, highlight that the accuracy of MTPs are insensitive to the number of elements' types and/or complexity of the structure. Although we have limited our AIMD simulations to less than 1 ps at 50 K, results presented up to this stage reveal that for the most of cases no additional AIMD trajectories are required to assess the phonon dispersions. Nonetheless, if required, the incorporation of additional high temperature AIMD trajectories in the training set is expected to improve the accuracy and stability as well. Besides the phonon dispersions, the phonons group velocity is another important phononic property, which can provide useful information concerning the lattice thermal conductivity. In Fig. 7 we compare the phonon group velocities predicted by the MTP method with those obtained within the DFPT using the PHONOPY code 34 . Remarkably, the passively trained MTPs can very closely reproduce the phonon group velocities for the studied samples. For some of the samples like graphene, penta-graphene and -SiP 2 , CN, BC 6 N 6 and transition metal dichalcogenides (MX 2 , M= Mo, W and X=S, Se, Te) excellent agreement between the MTP method and DFPT-based estimations are observed. It is clear that in the worst cases, like BrCuTe 2 , BC 3 , or C 2 N the group velocities are slightly shifted up or down and the results are free of substantial inaccuracies. Last but not least, we examine the agreement between the MTP-and DFPT-based approaches to calculate the free energy, heat capacity, and entropy from their statistical thermodynamic expressions using the PHONOPY code 34 . We compared the MTP-and DFPTbased estimations for the aforementioned thermal properties in Fig. 8. As it is clear, for these results the agreement between DFPT and MTP results is excellent. Presented results highlight the accuracy and stability of the MTP-based method to explore the phononic properties of complex 2D lattices. The required AIMD simulations to assess the phononic properties can be conducted relatively fast and are computationally less expensive than those commonly performed to examine the thermal stability (which usually requires more than 10 ps-long AIMD simulations). In comparison with DFPT method, AIMD simulations show lower sensitivity to the K-point grid and plane-wave cutoff energy. Moreover, the computational cost for the examination of supercell size effect on the PDRs is marginal for MTP-based approach, which is in contrast the bottleneck of DFPT method. One of the most salient advantages of the MTP-based method is that it typically yields smooth dispersions for the acoustic modes near the Γ-point. Differently when using the DFPT method for 2D materials, the convergence of flexural acoustic mode near the Γ-point may turn into a difficult task 55 . Nonphysical imaginary frequencies in the acquired PDRs may result in the instabilities for evaluating other thermal properties, particularly when calculating the lattice thermal conductivity. Worthy to point out that one of the reasons for appearance of imaginary frequencies in the PDRs might stem from slight inaccuracies in the atomic lattice. When using the DFPT method, small structural changes in the primitive unit cell, either in the atomic position or strain in the lattice size, demand entirely new calculations, whereas these effects can be checked with negligible computation costs using the MTP-based method. We should also observe that for highly symmetrical structures like graphene, transition metal dichalcogenides and penta-graphene, the standard DFPT method is by a large extent more computationally efficient. Actually the computational efficiency of MTP significantly enhances as the symmetry decreases, like that for haeckelites and C 7 N 6 structures or for materials with large primitive unit cells, such as: metal-or conductive-organic frameworks and graphyne/graphdiyne lattices. In the previous works 44,45 it was demonstrated that MLIPs can accurately reproduce the phononic properties of graphene as compared with the standard DFPT method, but with substantially higher computational costs, questioning the practical application of MLIPs for this purpose. However, our results for extensive 2D systems clearly confirm that for low-symmetry and highly porous atomic lattices, the MTP-based approach offers an accurate, simple and convenient alternative to commonly employed DFPT approach to access the phononic properties. This is therefore a novel application area of MLIPs in the materials science as they can accurately and efficiently substitute the standard DFPT method for the valuation of phononic properties of complex structures. Although, the majority of our results on the basis of short AIMD trajectories at the low temperature of 50 K reveal outstanding accuracy, it is nonetheless highly recommended to include AIMD trajectories at higher temperatures in the training set to ensure improved stability and accuracy as well. Moreover, to more accurately reproduce the high frequency optical modes, employment of larger supercells in the AIMD calculations is highly recommended. In this case, the increase in the computational cost can be partially compensated by shortening the AIMD simulation. In this work PBE method was employed in the AIMD calculations for the evaluation of phononic properties. However we note that the proposed MTP-based approach can be also effectively used for other exchange-correlation functionals.

Conclusion
Our extensive results for a wide variety of 2D materials highlight that machine-learning interatomic potentials trained over short ab-initio molecular dynamics trajectories are able to reproduce the phononic properties in a strikingly close agreement with those by DFPT-based results. The proposed methodology could therefore play a pivotal role to conveniently explore the phononic properties of a large variety of low-symmetry and porous nanomembranes, with a high level of accuracy and reproducibility. In recent years astonishing advances have been achieved on the synthesis of low-symmetry and highly porous atomic lattices, like metal-or conductive-organic frameworks and graphdiyne nanomembranes. It is clear that phononic properties of these novel systems can now be effectively explored using the approach proposed in this work, which would otherwise require excessive computational resources with DFT-based methods. Our study therefore proposes that MTP-based approach can accurately and efficiently substitute the standard DFPT method for the valuation of phononic properties of complex structures. To facilitate the practical application, the full computational details of the proposed approach are included in the data availability section.