Multi-objective parametrization of interatomic potentials for large deformation pathways and fracture of two-dimensional materials

This investigation presents a generally applicable framework for parameterizing interatomic potentials to accurately capture large deformation pathways. It incorporates a multi-objective genetic algorithm, training and screening property sets, and correlation and principal component analyses. The framework enables iterative definition of properties in the training and screening sets, guided by correlation relationships between properties, aiming to achieve optimal parametrizations for properties of interest. Specifically, the performance of increasingly complex potentials, Buckingham, Stillinger-Weber, Tersoff, and modified reactive empirical bond-order potentials are compared. Using MoSe2 as a case study, we demonstrate good reproducibility of training/screening properties and superior transferability. For MoSe2, the best performance is achieved using the Tersoff potential, which is ascribed to its apparent higher flexibility embedded in its functional form. These results should facilitate the selection and parametrization of interatomic potentials for exploring mechanical and phononic properties of a large library of two-dimensional and bulk materials.


INTRODUCTION
Molecular dynamics (MD) simulation based on force fields is a powerful tool for studying the temporal behaviors of materials at submicron scales. With continual improvements in hardware and algorithms, MD simulations are becoming increasingly accurate and widely adopted in several frontier problems in materials science and biology 1 . While such advances have greatly expanded the capability of MD simulations in size, timescale, and complexity, their predictive powers rely heavily on the accuracy of empirical interatomic potentials in approximating the interactions between atoms. Given the rapid emergence of new two-dimensional (2D) materials 2,3 that have demonstrated promising electrical, chemical, optical, thermal, and mechanical properties, an increasing demand for accurate interatomic potentials needs to be fulfilled to facilitate mechanistic understandings of their behaviors at scales representative of those used in applications.
Typically, interatomic potentials are formulated for a specific class of materials and are parameterized for a selected list of properties. Consequently, their accuracies on materials beyond the target class or for properties not included in the parametrization need further validation by more accurate methods, i.e., ab initio calculations. When those validations are conducted, in general they fail to achieve accurate predictions and reparametrization or new formulations are necessary. For instance, a reactive many-body potential parameterized for molybdenum disulfide 4 was found to yield artificial stiffening at large strain, and an ad hoc parameter-tuning was conducted to correct such behavior 5 . A more systematic interatomic potential parametrization would be indispensable in this case, but becomes a complex and specialized task that requires strong domain expertise and in most cases deep chemical intuition. One of the complexity of the parametrizing procedure is to reasonably capture nonequilibrium properties, such as vacancy formation energies and uniaxial tension behaviors at the same time. Typical training data (structures, energies, and bond stiffness of atomic clusters or primitive cells [6][7][8] ) are found insufficient to accurately reflect such properties 5,9 . Furthermore, there is little guidance beyond chemical intuition for choosing more appropriate training data, thus posing potential limitations on the accuracy and efficiency of the parametrization. Another complication arises due to the fact that interatomic potentials are often parameterized for a finite number of target properties, and some multi-objective optimization schemes may inevitably rely on human interventions. Specifically, a common approach, the weighted-sum method, converts the multi-objective problems into single-objective problems with user-defined, objective-specific weights [10][11][12] . However, the choice of a priori weights may bias the optimization 13 , thus limiting a holistic evaluation of the performance of interatomic potentials on various properties. This motivated researchers to formulate other optimization approaches, e.g., the Pareto front approach 14,15 . The last problem of the parametrization is to obtain a set of parameters for a chosen potential form. The selection of a potential form for a new material requires a vast domain knowledge of not only the physics of the material at hand, but also the specific details of such a form. This limitation inevitably hinders the big-picture view of whether interatomic potentials can be parametrized to simulate the behaviors of some special class of materials, e.g., 2D materials, whose atomic structures and properties are distinct from bulk crystals. As a result, it prevents a direct comparison of performance between various potentials for the same material.
Despite several successful parametrizations over the past decade 4,10,[16][17][18][19][20][21][22] , an in-depth evaluation of the suitability of existing interatomic potentials for the prediction of phase transition and fracture of 2D materials is still lacking. Herein, we propose a robust parametrization method built upon density functional theory (DFT) data sets (considered as ground truth) and the evolutionary multi-objective optimization algorithm, NSGA-III 23 . Similar to other genetic algorithms, NSGA-III avoids the dependence on gradient computation, hence it can be applied to any functional form (potential). In addition, this algorithm enables a generation of more widely distributed points on the Pareto front in the criterion space, allowing a more thorough search for an optimum interval. As a result, this algorithm, along with adoption of a machine-learning-inspired protocol, shows good transferability and performance, and offer higher parametrization flexibility. The proposed method is applied to several interatomic potentials of increasing complexity, namely, Buckingham 24 , Stillinger-Weber (SW) 7 , and Tersoff 8 for the structural, mechanical, and thermal properties of monolayer 2D materials in both the equilibrium and nonequilibrium regimes (see "Parametrization method" for our classification of properties). The modified reactive empirical bond-order potential for transition metal dichalcogenides (referred to as REBO-TMDC) 4 is also considered for comparison. As a case study, we perform the parametrization of MoSe 2 and prioritize its mechanical behavior. We use the structure and stability of various surfaces, and thermal properties to examine the interatomic potentials' transferability. All of the parameterized potentials have better accuracy in nonequilibrium properties when compared with existing MoSe 2 potentials, highlighting the effectiveness of the proposed parametrization method. We further explore the parametrization flexibility of the selected interatomic potentials by conducting correlation and principal component analyses on their prediction errors, which reveals a positive correlation between the complexities of interatomic potentials, their flexibility, and their performances on MoSe 2 . Together, these results suggest a robust potential parametrization approach and a quantitative potential selection criterion, which may be generalized for a wide range of materials and materials properties beyond those explored in this study.

Selection of materials and interatomic potentials
Among existing 2D materials, TMDC are one group of materials described by MX 2 , where M is a transition metal (Mo, W, etc.) and X is from the oxygen family (S, Se, etc.). Most TMDCs in monolayer form are semiconducting with strong photoluminescence 25 , thus making them promising candidates for applications, such as transistors 26 , photodetectors 27 , supercapacitor electrodes 28 , and solar cells 29 . We set out to parameterize interatomic potentials for TMDCs with a focus on failure-related properties, which are critical to the stability and reliability of systems that require frequent mechanical deformation, e.g., flexible electronics. Specifically, we selected monolayer MoSe 2 and its stable phase 2H for which existing interatomic potentials, parameterized primarily for equilibrium properties (structures at equilibrium, phonon dispersion, etc.) 16,30 , show major deviations in comparison to DFT for nonequilibrium properties, including surface stability and uniaxial stress-strain response. This implies necessities to expand the parametrization to the nonequilibrium regime, defined as states with large perturbations from the equilibrium positions or systems possessing point defects.
We selected interatomic potentials according to their applicability for TMDCs and ability to describe atomic chemical environments. To the best of our knowledge, existing parametrizations for TMDCs include SW potentials for the mechanical and thermal properties of MoS 2 ( 10,16,18,20 ), MoSe 2 ( 16,31 ), and WSe 2 ( 31 ); a Tersoff potential for the thermal properties of WSe 2 21 ; a ReaxFF potential for the mechanical and transitional behaviors of MoS 2 ( 17,32 ); and a REBO-TMDC potential for the interfacial and mechanical properties of MoS 2 ( 4,5 ). Those interatomic potentials can be segmented into cluster potentials (SW), cluster functionals (Tersoff), and reactive cluster functionals (ReaxFF and REBO-TMDC) with increased levels of complexity and capabilities 33 . Specifically, from the simplest pair potentials (e.g., Lennard-Jones potentials), cluster potentials introduce many-body (>2) interactions, cluster functionals incorporate bond-order terms for coordination-dependent bond strength, and reactive cluster functionals enable simulation of chemical reactions. Herein, we chose SW, Tersoff, REBO-TMDC, and also the Buckingham potential 24 , a simple pair potential widely used for ionic crystals. The formulations of the selected interatomic potentials are detailed in Supplementary Note 4.

Parametrization method
Parametrization of the selected interatomic potentials was performed in an iterative manner. Each iteration consists of three steps, referred to as training, screening, and evaluation ( Fig. 1). In the training step, the parameters of the interatomic potentials are optimized with the multi-objective genetic algorithm to minimize the errors for a selected group of properties in comparison to ab initio data. Next, the optimized parameters are screened for the remaining properties with a set of user-specified maximum percentage errors to identify promising candidates. Such a protocol is inspired by machine-learning methods, in which the full data set is separated into training and validation sets to balance underfitting and overfitting 34 . The evaluation step, including correlation and principal component analyses, is deployed to identify the correlation relationships between properties and redundancy in them. The information is used to (a) guide the selection of training and screening properties for the next iteration, and (b) to quantify parametrization flexibilities of interatomic potentials.
To predict the failure of MoSe 2 , we hypothesized that an essential list of materials properties (Fig. 2), in both the equilibrium and nonequilibrium regime, needs to be captured by the interatomic potentials. In the equilibrium regime, we selected the lattice structure and cohesive energy of MoSe 2 at equilibrium, equation of state (near equilibrium), elastic constants (C 11 and C 12 ), surface energies (armchair (AC) and zigzag (ZZ) surfaces), and surface stability. For nonequilibrium properties, we selected the following: bond dissociation energy landscapes (along the AC and ZZ directions), vacancy formation energies (seven types), and uniaxial stress-strain curves under elastic instability and soft mode (along the AC and ZZ directions, see Supplementary Note 1 for discussion on elastic instability and soft mode). Moreover, the energy landscape of a 2H-1T phase transition was included to characterize phase changes observed in TMDCs under uniaxial and biaxial strain 35 , doping 36 , or vacancy reorganization 32 . For the Tersoff and REBO-TMDC potential, an additional set of properties for Se systems is needed. Following Chan et al. 21 , we selected the structures and cohesive energies of Se clusters (Se 2 , Se 3 , Se 6 , and Se 8 ), stability and dissociation of Se 6 and Se 8 , and expanded the nonequilibrium regime by including the dissociation energy landscapes for Se 3 , Se 6 , and Se 8 (see Supplementary Note 3). Ab initio calculations at the DFT level of theory were carried out on all the above properties, which were used as ground truth for parameterizing the interatomic potentials.
The properties are divided into optimization and validation sets. The former is further segmented into training and screening sets, and the optimized parameters are then applied into the validation set after the optimization is finalized. Some training properties imply screening properties and vice versa, although they contain some inherently different information, e.g., the near-equilibrium data points of uniaxial stress-strain curves versus elastic constants, uniaxial stress-strain curves under elastic instability versus curves under soft mode (see Supplementary Note 1). We started an optimizing process by selecting simple properties for training, such as bond dissociation energy landscapes, and more complex properties, e.g., vacancy formation energies, as screening properties 17,18 . Another factor that affected this choice is the size and type of calculation needed to be carried out for a given property. For example, a single-point energy calculation would be favored over one that requires MD equilibration, and a system consisting of a primitive unit cell will be prioritized unless a larger system with more atoms would provide more representative data. This criterion accelerates the optimization since the simpler properties often require single-point calculations that are faster  In the training step, the multi-objective genetic algorithm NSGA-III is used to optimize the parameters of the interatomic potentials for the training properties ( Fig. 2). MD simulations are integrated into the genetic algorithm workflow for the evaluation of errors. After reaching a predefined number of generations, the optimization terminates, and the optimized parameters are passed to the screening step for the evaluation of the screening properties (Fig. 2). The training and screening properties together are considered for the selection of promising candidates with percentage errors for all properties within a user-specified threshold. Meanwhile, correlation and principal component analyses are carried out in the evaluation step to reveal the correlation relationships between properties and the performance of interatomic potentials. This information guides the selection of training properties for the next iteration of parametrization.
to perform than those that require energy minimization or extra MD steps. These extra steps occasionally suffered from convergence issues due to the emergence of unphysical interatomic potential parameters during the first several steps of optimization. Moreover, properties of the same general attribute, i.e., uniaxial stress-strain curves, were used in both stages with different levels of perturbation. Specifically, stress-strain curves under elastic instability were used as training properties and the curves under soft mode were used as screening properties (see Supplementary Note 1). As such, we selected the following properties as training properties: lattice structure and cohesive energy at equilibrium, equation of state, bond dissociation energy landscapes along both the AC and ZZ directions, and uniaxial stress-strain curves under elastic instability along both the AC and ZZ directions. We note that such a choice is not fixed. Rather, its effect on the parametrization results can be determined from the correlation analysis, as discussed later.
Our training step involves solving a multi-objective optimization problem for which two major approaches prevail: scalarization and vector optimization methods 13 . Scalarization methods convert multi-objective optimization problems into single-objective optimization problems using methods such as weighted-sum, and are predominantly used for parameterizing interatomic potentials. On the other hand, vector optimization methods treat each objective independently and aim at exploring the Pareto optimal solutionssolutions that cannot be further improved without worsening at least one objective. These methods assign equal importance to each objective and allow users to emphasize or exclude certain objectives without biasing the optimization. This is helpful in preventing solutions from being swamped out, but may be a limitation if relative importance between objectives is known. We overcame this limitation in the screening step by assigning a lower objective-specific percentage error (discussed later), if one wants to focus on a specific property. Among various vector optimization methods, we adopted a multi-objective genetic algorithm named NSGA-III 23 . As a genetic algorithm, it conducts a global optimization through an iterative process that loosely mimics Darwin's theory of natural selection, i.e., via mutation, crossover, and selection operations on a population of individuals. Moreover, the algorithm incorporates a nondominated sorting procedure and a niche-preservation operator to identify nondominated individuals that are well-spread in the criterion space. As a result, this specific algorithm shows superiority by sampling more widely distributed points on the Pareto front. To the best of our knowledge, we report the first application of NSGA-III for parameterizing interatomic potentials. In most optimization problems with three to ten objectives, it outperforms its predecessor, NSGA-II 37 , which has been the predominant algorithm applied in similar problems 14,15 . Moreover, we found that it is more efficient in overcoming local minimum states during the optimization in comparison to a hierarchical genetic algorithm optimization framework 21 (see Supplementary Note 2).
As shown in Fig. 1, the optimization starts from an initial population of individuals (i.e., sets of parameters for a given interatomic potential) that are randomly generated within a predefined range according to potential-specific requirements (see Supplementary Note 4). In each generation, mutation and crossover operations are first conducted on the current population according to specified probabilities. The parameters after mutation and crossover are not allowed to exceed the predefined initialization range. The value to be minimized for each objective during the optimization, i.e., the fitness value, is defined as follows: where x represents an individual whose dimension equals the number of parameters to be optimized for a given interatomic potential, v i j (x) denotes predictions from interatomic potentials with parameter set x for objective i,v j i represents ab initio results for objective i, superscript j represents the jth point for objective i, and w i j is the corresponding weight for the jth point. The weight is used to emphasize certain regions of the training data, e.g., the elastic regime of the stress-strain curves. The values of the weights are provided in Supplementary Note 10, and are found to have an insignificant effect on parametrization results. The pointwise errors are squared, scaled, and summed over all points for objective i. Thus, the fitness value is the squared error for objectives with K = 1 (e.g., elastic constant C 11 etc.), and is the sum of (weighted) squared errors for objectives with K > 1 (e.g., stress-strain curves). For each parameter set x, there are M t fitness values corresponding to the M t training objectives (properties). Individuals of the current generation are ranked based on whether they are dominated by other individuals. An individual x is nondominated if and only if there does not exist another individual x Ã in the current population such that All nondominated individuals (Pareto optimal solutions) are ranked lowest and are selected first, followed by individuals that are only dominated by those Pareto optimal solutions and so on. The selected individuals are passed to the next generation, where the same protocol repeats until the predefined total number of generations is reached. For each interatomic potential, multiple optimizations are carried out concurrently with different random seeds to explore a large parameter space and various optimization paths.
Prior to the screening step, the optimized individuals from multiple optimizations are gathered, and the fitness values for the remaining properties (screening properties in Fig. 2) are calculated following Eq. (1). Those values are combined with the fitness values of the training properties to form a matrix of dimension N × M, where N is the total number of optimized individuals and M is the total number of properties. This matrix is then screened with the criterion defined as follows: where f i x ð Þ, w j i , andv j i follow the same definitions as in Eq. (1), and p i is an objective-specific percentage value. The criterion defines a maximum percentage error for objectives with one data point (K = 1), and resembles the sum of squared percentage error (SSPE) for objectives with multiple data points. Indeed, this criterion and the SSPE measurement agree well in magnitude and variation according to a sampling test on a harmonic function (see Supplementary Note 6). The parameters p i reflect the relative importance of each objective and provide a desirable level of flexibility for the parametrization without biasing the optimization, which is a major advantage over the prevalent weighted-sum method. The parameter sets that pass the screening step are deemed promising candidates.
Our approach is unique in several respects. In our implementation of screening sets, we used explicit criteria for selecting fittest interatomic potential parameters and evaluating the parametrization flexibility. This emerges when we evaluate the effect of one criterion on the validation properties or the entire set of properties. This approach is different to most, if not all, parametrization of empirical interatomic potentials in which the validation sets are used for examining the presumptive transferability instead of guiding the optimization. For instance, the validation test for a recently developed CHARMM force field 38 is conducted on chemically similar species with respect to the optimization sets. Such a validation should and indeed does reveal good transferability as it follows the underlying assumptions of the CHARMM General Force Field 39 . By contrast, several interatomic potentials for MoSe 2 , as reported in the literature, have poor transferability (see the next section). In the present approach, X. Zhang et al. the improvement on transferability can be done by redefining allowable errors for the screening sets (and training sets) that have more impact on a desired validation property (e.g., phonon dispersion and thermal conductivity). Thanks to this protocol, we were able to ascertain transferability and infer directions for balancing parametrization tradeoffs of specific properties. Lastly, we introduced correlation and principal component statistical analyses to ascertain correlation relationship between properties and infer parametrization flexibilities of interatomic potentials, both of which remain unexplored in the literature. We note that the statistical information closes the parametrization loop by guiding the selection of training and screening properties for the next iteration, as further elaborated in subsequent sections. Table 1 and Fig. 3 summarize the predictions of the parameterized interatomic potentials in comparison to ab initio calculations, as well as existing SW 16 and SNAP 30 potentials for monolayer MoSe 2 . In selecting the parameter sets for each interatomic potential with Eq. (2), we emphasize the accuracy of uniaxial stress-strain curves (under soft mode), while maintaining other properties to be at least within 90% error. The parameterized potentials are listed in Supplementary Note 5. Not surprisingly, since the SW potential parameterized by Kandemir et al. and the SNAP potential parameterized by Gu and Zhao are primarily trained for equilibrium properties, including structures and thermal transport, they show limited accuracy for the nonequilibrium properties studied herein: vacancy formation energies (Table 1) deviate from ab initio results (the surface and vacancies structures cannot be equilibrated with the SNAP potential); bond dissociation energy landscapes along the AC (Fig. 3a) and ZZ (Fig. 3b) directions deviate from ab initio curves at strain > 0.1; uniaxial stress-strain curves along both the AC (Fig. 3e) and ZZ (Fig. 3f) directions deviate significantly from ab initio results at strains larger than 0.05. In comparison, the parameterized interatomic potentials here reported yield predictions that are vastly more accurate for those properties due to the augmented ab initio training data in the nonequilibrium regime and an explicit screening step that further defines maximum allowable errors. Figure 3d shows the phase transition energy landscape identified by climbing image nudged elastic band simulations 40 (see Supplementary Video 1 for the movie of this simulation). Under the screening criteria, which prioritize stress-strain responses, all the parameterized interatomic potentials incorrectly predict the 1T phase to be the minimum energy state. For the Tersoff potential, an individual with the correct relative energy between the 2H and 1T phase have less accuracy on other properties (see Supplementary Note 7). The results suggest an important limitation of the Tersoff potential, which is the capturing of the phase transition energy landscape. This finding is confirmed by the correlation analysis discussed in the next section.

Parametrization results
Notably, the Tersoff potential has the overall best performance among the selected interatomic potentials. It provides a smooth uniaxial tension curve closely matching the ab initio results, whereas the Buckingham and SW potentials predict an artificial 2H-1T phase transition during uniaxial tension at~0.2 strain, manifested as kinks in Fig. 3e, f. This data seems to advocate a positive correlation between the complexity of the interatomic potential and its overall accuracy for monolayer MoSe 2 . Such observation is further corroborated by the correlation analysis, as described in the next section. We note the importance of direct force fitting (i.e., stress-strain curves as training data) for the Tersoff potential to achieve good accuracy for uniaxial tension response. Similarly, forces and other higher order derivatives of energies were found critical for accurate predictions of phonon dispersion and thermal transport in crystalline Si and Ge 11 . Thus, inclusion of force fitting during parametrization should result in better transferability of the interatomic potentials.  A dash denotes that the system with surface or vacancy undergoes significant structural changes, and cannot be equilibrated with the SNAP potential.

Correlation and principal component analyses
X. Zhang et al. in the M-dimensional design space enables a quantitative assessment of the relations between the M properties through statistical analysis. From the N × M matrix, we construct an M × M correlation matrix where each element R ij is the Pearson correlation coefficient between property i and property j. A graphical representation of those coefficients for the Tersoff potential is shown in Fig. 4. The correlation coefficients in Fig. 4 correspond to individuals with percentage errors <100% for all properties. Such setting aims at exploring regions in the criterion space where promising candidates are selected. R ij ranges from −1 to 1, denoting strongly negative and positive correlations, respectively. In the context of our problem, a positive R ij between property i and j indicates simultaneous increase or decrease of their prediction errors, and hence nonconflicting relations between the two properties. In comparison, a negative R ij suggests a conflicting relationship between property i and j: the errors for property i cannot be minimized without compromising the accuracy of property j.
As shown in Fig. 4, the Tersoff potential (the bottom triangle) possesses several pairs of properties that are strongly positively correlated. Conforming to chemical intuition, properties of similar nature have strong positive correlations, e.g., bond dissociation energies along the AC and ZZ directions (R = 0.7), AC and ZZ surface energies (R = 0.9), and vacancy formation energies of Mo and Mo 2 F (R = 1). Similar relations were identified for the Buckingham and SW potential (see Supplementary Note 8). The correlation matrix also reveals conflicting properties that cannot be directly deduced from chemical intuition. Specifically, the phase transition energy landscape is conflicting to almost all properties. This agrees with our observation of accuracy degradation of most properties when the phase transition energy landscape is prioritized (see Supplementary Note 7). Thus, the correlation matrix offers a direct gauge of the accuracy (and hence the parameters) of the parameterized interatomic potential on all properties should the relative importance between properties be different. Fig. 3 Prediction results of the parameterized interatomic potentials in comparison to ab initio results, as well as the existing SW 16 and SNAP potentials 30 . a, b Bond dissociation energy landscape along the armchair (a) and zigzag (b) directions. Two snapshots along the landscape are shown and correspond to strain of 0 and 0.6, respectively. c Equation of state. d 2H-1T phase transition energy landscape. The two snapshots from left to right correspond to reaction coordinates of 0 and 1, respectively (see Supplementary Video 1 for the movie of this simulation). e, f Uniaxial stress-strain curve at 1 K in comparison to ab initio curves under soft mode (see Supplementary Note 1) along the armchair (e) and zigzag (f) directions. The snapshots in e and f show the formation of the 1T phase during uniaxial tension for the Buckingham potential, which is also predicted by the SW potential during uniaxial tension along the zigzag direction. Legends of a, b, d-f are identical to that in c, and are omitted for clarity. In the atom snapshots, Mo atoms are colored cyan, and Se atoms are colored orange.
To assess the parametrization flexibility of any interatomic potential, we propose evaluating a quantity F defined as For an ideal interatomic potential approaching an ab initio level of accuracy, there should exist a region in the criterion space, where the prediction errors of all properties have strong positive correlations, and thus can be minimized simultaneously, which allows a large level of flexibility for parametrization. For the ideal interatomic potential, F ideal = (M − 1)/2. We sampled five regions in the criterion space near 100% percentage errors (p i = 80-120%) with enough individuals (>30) in each region. F Tersoff , F SW , and F Buckingham were found to be 1.44 ± 0.12, 0.71 ± 0.25, and 0.77 ± 0.12, respectively. Notably, the Tersoff potential has the highest flexibility, consistent with our observation that the parameterized Tersoff potential predicts the most accurate uniaxial stress-strain curves.
To further explore the intrinsic relationships between properties, we conducted principal component analysis on the correlation matrix. It was originally proposed to identify redundant objectives during multi-objective optimizations 41 , and is used herein to find redundant properties, i.e., properties that can be automatically captured if essential properties are captured with either the training or screening step. The protocol is detailed in Supplementary Note 9. The analysis reveals that all the training data are nonredundant for all the interatomic potentials, thus indicating optimized training properties. Furthermore, it shows that the formation energies of certain vacancies are redundant with respect to others (see Supplementary Note 9), in agreement with chemical intuition.

Validation and transferability of parameterized potentials
We evaluated the validity and transferability of the optimized Tersoff potential on the edge stability and thermal properties of monolayer MoSe 2 , which are relevant for applications in nanoelectronics 16,30 and catalysis 42 , but are not parameterized within the scope of this study. This can also be referred as the "test" data set in other machine learning frameworks 34 . Figure 5 shows ab initio molecular dynamics (AIMD) and Tersoff predictions on the stability of various edge configurations at 300 0 K and elevated temperatures. Those configurations correspond to the Mo-Klein, Mo-ZZ, Se-ZZ, and the AC edge, which were identified by scanning transmission electron microscopy (Fig. 5) in nanoporous MoS 2 films grown with molecular beam epitaxy under high Mo flux 42 . The Tersoff potential shows a decent level of transferability owing to its higher flexibility in its functional form: it Fig. 4 Correlation matrix of the Tersoff potential. The Pearson correlation coefficient for the corresponding pair of properties, ranging from −1 to 1, is shown in each cell. When calculating the correlation coefficients, individuals with percentage errors > 100% for any property were excluded to explore the region in the criterion space, where promising candidates are selected. Two correlation matrices are shown; the bottom corresponds to the as-optimized Tersoff potential and the top corresponds to the same population with additional screening on the three acoustic phonon modes, ZA, LA, and TA. Since each correlation matrix is symmetric, only half of the correlation matrix is shown and the diagonal components (always equal to 1) are removed for clarity. For visualization purposes, all correlation coefficients are rounded to a decimal. Thus, a value of "−0" means the true correlation coefficient is within [−0.05, 0]. The cells are colored according to their correlation coefficients. Uniaxial tension items herein refer to curves under soft mode (see Supplementary Note 1).  42 and are shown herein. Red dashed lines highlight the Mo atomic layers, and the scale bars represent 0.5 nm. The same edge configuration was equilibrated with ab initio MD (AIMD; for 1 ps) and MD (500 ps) at the specified temperature. An "unstable" configuration is defined as a configuration that undergoes breakage and reformation of chemical bonds during equilibration (see Supplementary Videos 2-5 for the movies of those simulations). We note a longer simulation time for MD to reveal the unstable edge configurations due to limitations of interatomic potentials to reflect electronic interactions embedded in the AIMD simulations. In the atomic snapshots, Mo atoms are colored cyan, and Se atoms are colored orange. The Tersoff-ZTL corresponds to the parameterized potential with the training data of "Tersoff" plus the three acoustic phonon modes, i.e., ZA, LA, and TA. Reprinted (adapted) with permission from Zhao et al. 42 Copyright (2018) American Chemical Society.
reproduces AIMD predictions for all the edges at 300 0 K, for the Mo-Klein and Se-ZZ edges at elevated temperatures (650 and 750 0 K, respectively), but tends to overstabilize the Mo-ZZ edges at 650 0 K (see Supplementary Videos 2-5 for the movies of those simulations). Figure 6a shows the phonon dispersions predicted by the parameterized Tersoff potential in comparison to the ab initio results. The Tersoff potential predicts no negative frequency and correct Γ point for acoustic bands albeit lower frequencies for the out-of-plane (ZA) mode and a smaller phonon band gap. Such inconsistency results in lower in-plane thermal conductivity in comparison to first-principles calculations 43 and experimental measurements on a suspended monolayer MoSe 2 membrane 44 , as shown in Fig. 6b. However, we note that the excellent accuracy of the longitudinal (LA) and in-plane transversal acoustic band (TA) was captured by the Tersoff potential, due to the inclusion of force fitting in training, i.e., uniaxial stress-strain curves along two directions.
We next discuss how the correlation and principal component analyses are employed to close the parametrization loop and improve the accuracy of phonon dispersion. We screened the optimized Tersoff population on the three acoustic phonon modes, and conducted a correlation analysis. The correlation matrix (Fig. 4, top triangle) shows that the ZA mode is more conflicting to other properties in comparison to the TA and LA modes, i.e., for C 11 , uniaxial stress-strain curves along the ZZ direction, and stability at 300 0 K. Such results suggest that adding the ZA mode into the training data will increase the accuracy of ZA mode at the expense of decreased accuracy on the aforementioned properties. However, the worsening effect may be mitigated by including the other two modes that possess relatively positive correlation relationships. We carried out two iterations, one with the addition of the ZA mode and the other with the addition of all three modes to the training data set. The results support the above statement (see Supplementary Note 11). Specifically, the iteration with all three modes, referred as Tersoff-ZTL, resulted in a more accurate ZA mode (Fig. 6a) and thermal conductivity (Fig. 6b) with minimum deterioration of other properties (Supplementary Fig. 9 and Fig. 5). We note that the bending rigidity of 2D materials is directly related to the ZA mode 45 . Thus, the Tersoff-ZTL should possess improved accuracy on bending rigidity.

DISCUSSION
We propose a robust approach of parameterizing interatomic potentials. It incorporates the multi-objective genetic algorithm NSGA-III, a machine-learning-inspired protocol, and a correlation and principal component analyses framework. Using monolayer MoSe 2 as a testbed, we demonstrate the effectiveness of the proposed approach in capturing properties of monolayer MoSe 2 in both the equilibrium and nonequilibrium regimes. Compared with existing parametrization methods, our approach incorporates a more efficient optimization algorithm, provides more flexibility for balancing the tradeoff and priority of specific properties without biasing the optimization, and shows good transferability for various interatomic potentials with different levels of complexity. In all cases, the method is straightforward to implement, given the appropriate computer codes. Moreover, this approach enables the exploration of the intrinsic relationships between properties through correlation and principal component analyses, which is absent in other parametrization frameworks, such as GARField 46 and Paramfit 47 . With the correlation matrix, one can (a) evaluate the feasibility of improving the parametrization of a given potential; (b) assess the parametrization flexibility of a given potential with the value of F. At this stage, the analyses are used to close the parametrization loop through selections made by the user. An automated workflow could be developed by adopting other machine learning algorithms, such as logistic regression and decision trees on a cross-validating set, which would further minimize human intervention 34 . Such approaches are left to future studies.
In the approach here presented, the choice of interatomic potentials and parametrization parameters constitute an iterative process. One should start by using computationally inexpensive properties (see the "Parametrization method" section) for training, to cover a wider range of configurations and achieve better optimization efficiency. The correlation and principal component analyses are then employed to guide the selection of training and screening properties for the next iteration, e.g., exclude redundant properties from the training set or include properties into the training set after evaluating the effect on other properties. We note that simple expansion of the data base, without considering their correlation with other properties of interest, is not a suitable approach.
We identified intrinsic conflicting relationships between certain properties for the parameterized interatomic potentials and attributed such behaviors to the limitations of the functional forms. To examine if more sophisticated functional forms can alleviate this issue, we parameterized a modified reactive emperical bond-order potential for TMDCs (referred to as REBO-TMDC) 4 , using the same training data as Tersoff-ZTL. In comparison to Tersoff-ZTL, the REBO-TMDC shows improved accuracy on phase transition, but decreased accuracy on several other properties (see Supplementary Note 12). This suggests that  43 , as well as experimental measurements on suspended monolayer MoSe 2 using optothermal Raman techniques 44 . The Tersoff-ZTL corresponds to the parameterized potential with the training data of "Tersoff" plus the three acoustic phonon modes, i.e., ZA, LA, and TA.
the intrinsic conflicting roles of properties are always present and that more sophisticated interatomic potentials do not necessarily translate into better accuracy across all properties. Rather, as can be seen from the example of Tersoff and Tersoff-ZTL, attention needs to be paid as to how properties correlate for each functional form (potential).
We leave the parametrization of more complex functional forms, with a larger number of parameters, e.g., ReaxFF 48 , for future work. Nevertheless, when equilibrium properties are the primary interest of parametrization, simpler interatomic potentials are worth exploring due to the ease of parametrization and better computational efficiency. Indeed, common interatomic potentials were found to be overdesigned for the purpose of exclusively modeling equilibrium properties, including phonon dispersion and thermal transport in crystalline Si and Ge 11 , thus suggesting sufficient flexibilities for parametrization. In either circumstance, our approach offers a framework for future studies aiming at expanding the capability of empirical interatomic potentials for quantifying unconventional chemical and physical phenomena in emerging new materials.
We highlight the better performance of the NSGA-III algorithm over several existing multi-objective global optimization algorithms. We also note that other multi-objective global optimization algorithms, e.g., MOES 49 and GARField 46 , have been developed for parameterizing ReaxFF for molecular crystals and SiC. Future work should perform comparison among existing approaches. Furthermore, for TMDCs, the transferability test can be expanded to other DFT calculations and in situ transmission electron microscopy observations, including vacancy induced phase transition 32 , formation of inversion domains 50 , atomic morphologies of the crack tip 5 , etc., which we leave for future exploration.

Ab initio calculations
The training data for the optimization were created by ab initio calculations. These simulations were carried out using the density functional approach via SIESTA 4.0.2 software 51 . We applied the nonspin-polarized generalized gradient approximation in the Perdew-Burke-Ernzerhof (PBE) form 52 together with the split polarized valence doublezeta (DZP) basis set 53 . For the molybdenum and selenium atoms, nonrelativistic norm-conserving Troullier-Martins pseudopotentials 54 were utilized. The energy shift and mesh cutoff were selected to be 250 eV and 300 Ry (~4081 eV), respectively, at which energy convergence was attained. Geometry optimization was conducted without any symmetry constraints, until the forces acting on the atoms became <0.01 eV Å −1 . The interaction between monolayers (or molecules) were prevented by a 40 Å vacuum layer. A monolayer thickness of 7.726 Å was used to calculate perarea quantities (e.g., monolayer stresses) 55 . To achieve accurate electronic structure calculations, we allowed a 15 Å cutoff for the set of k-points in the first Brillouin zone. The resultant k-grids were chosen in an optimal way, according to the method of Moreno and Soler (which utilized an effective super cell close to spherical shape, thus minimizing the number of k-points for a given precision) 56 . The self-consistent and the conjugate gradient minimization schemes were employed for the electronic structure calculation and for the geometry optimization, respectively. The cohesive energy of a compound was computed from where E pristine is the energy of the compound, E Mo and E Se are the energies of an isolated Mo and Se atom, and n Mo and n Se are the numbers of corresponding atoms in the compound. The elastic constants were extracted from uniaxial stress-strain curves in the small-deformation regime. A fitting procedure reported by Cooper et al. 57 was used to extract the polynomial of the finite-deformation Green tensor of different orders, and the second-order terms were used in the screening process. Vacancy formation energies were calculated with the following equation: where E defected is the energy of the defected system, E pristine is the energy of the pristine system, n Mo and n Se are the number of missing Mo and Se atoms in the vacancy, and μ Mo and μ Se are chemical potentials for Mo in its stable BCC structure and Se in Se 8 rings, respectively.

Phonon dispersion calculation
In order to obtain the phonon dispersion curves and densities of states, the theory of lattice dynamics on a Born-Oppenheimer surface was applied.
Assuming that atomic displacements are in the form of plane-wave functions with different wave numbers, lattice vibration turns into an eigenvalue problem. Siesta 51 offers utility functions (e.g., vibra and fcbuild) that displace each atom in the monolayer and measure the reaction of other atoms to form the stiffness matrix. The wave dispersion along the path connecting symmetry points of the hexagonal lattice Γ−M−K−Γ in the first Brillouin zone for MoSe 2 monolayers were investigated. These points are located at (0,0,0), (0.5,0,0), and (0.333,0.333,0) in reciprocal space. These symmetry directions stem from the hexagonal lattice 2H MoSe 2 , similar to graphene. To sample the dispersion within the first Brillouin zone, a super cell of 4 × 4 × 1 repetitive unit cells was utilized to include all possible attenuations of the real-space force constants within it. To avoid interactions between monolayers, a vacuum layer of 40 Å was set.

Molecular dynamics simulations
MD simulations were conducted with LAMMPS 58 (3Mar20, serial version for optimization and mpi version for simulating larger systems). To compare MD simulations with ab initio calculations, we used the same atomic systems for most objectives except for the lattice structures and cohesive energies, in which we enlarged the size of the system for better sampling. For energy landscapes (equation of states, bond dissociation, phase transition, and dissociation of Se clusters), single-point calculations were performed on the equilibrated structures from ab initio calculations without energy minimization. For the remaining objectives, an energy minimization step was carried out on the input structures with the conjugate gradient algorithm (energy tolerance 0 eV, force tolerance 10 −10 eV Å −1 ) before calculating the energies. For simulations with MD steps, a time step of 1 fs was used. Phonon dispersion calculations were performed with phonopy (2.4.1.post5) 59 . Thermal conductivity was calculated using the equilibrium Green-Kubo formalism 60,61 . A monolayer MoSe 2 flake of 2.3 by 2.3 nm was first equilibrated with an NVT ensemble for 0.1 ns, followed by an NVE step of 1 ns during which the ensemble average of the autocorrelation of the heat flux was measured for calculating the thermal conductivity. We computed in-plane thermal conductivities as the average over conductivities along the AC and ZZ direction. The thermal conductivity at a given temperature was further averaged over 6 replicas with different initial random velocities. Atomic visualizations were created with OVITO 62 .

Nudged elastic band simulations
Climbing image nudged elastic band simulations were carried out with LAMMPS 40 . Fourteen replicas including the initial and final configurations were used. The energy and force cutoff for energy minimization were selected to be 0.1 eV and 0.01 eV Å −1 . The spring constant for the nudging force was set as 1 eV Å −1 . The minimum energy path was found using the Tersoff potential and was adopted by ab initio simulations and other interatomic potentials.

Optimization parameters and system setup
The population size was set to be 156 for the genetic algorithm optimizations following Deb and Jain 23 . Each optimization was conducted for 500 generations with which the optimization converged. We submitted 20 runs concurrently with different random seeds. We used a simulated binary crossover operator with a crossover probability of 1 and a crowding degree, η of 30. For mutation operations, we used polynomial mutation with a mutation probability of 1 and a η value of 20. We stored and output statistics of the entire population for every certain number of generations to monitor the optimization progress. After all runs finished, we combined the optimized parameters from all runs for the subsequent screening process. The optimization environment was set up with Python (3.7.7) and the genetic algorithm optimization was based on the DEAP framework (1.3.0) 63 . Parameter initialization, genetic algorithm operations, and calculations of statistics were conducted with DEAP. When the evaluation for an individual was needed, the code initiated a LAMMPS calculation via a system-level call, and read the output from LAMMPS log files.
Parallel programing was enabled by SCOOP (0.7.1.1) 64 to offer accelerated performance on supercomputer clusters. The optimization times scaled with the complexity of the interatomic potentials and were in the range of several hours to several days.

DATA AVAILABILITY
The parameterized interatomic potentials, as well as the corresponding potential files, are documented and included in the Supplementary Information.