Modeling and simulations for 2D materials: a ReaxFF perspective

Recent advancements in the field of two-dimensional (2D) materials have led to the discovery of a wide range of 2D materials with intriguing properties. Atomistic-scale simulation methods have played a key role in these discoveries. In this review, we provide an overview of the recent progress in ReaxFF force field developments and applications in modeling the following layered and nonlayered 2D materials: graphene, transition metal dichalcogenides, MXenes, hexagonal boron nitrides, groups III-, IV- and V-elemental materials, as well as the mixed dimensional van der Waals heterostructures. We further discuss knowledge gaps and challenges associated with synthesis and characterization of 2D materials. We close this review with an outlook addressing the challenges as well as plans regarding ReaxFF development and possible large-scale simulations, which should be helpful to guide experimental studies in a discovery of new materials and devices.


Introduction
The isolation of graphene in 2004 [1][2][3] has sparked a big wave of interest in two-dimensional (2D) materials, which have rapidly evolved into a broad range of emergent applications in the past decades. This discovery has theoretically and experimentally opened a wide design space to new synthetic '2D materials beyond graphene' , such as transition metal dichalcogenides (TMDs) [4], topological insulators [5][6][7][8], MXenes [9], hexagonal boron nitrides (h-BN) [10], 2D metals [11] and their nitrides (2D Ga, 2D GaN, etc) [12]. The atomically thin structure of these 2D materials offers a fertile ground for researchers to engineer and tailor the traits of 2D materials, such as bandgap, surface and edge reactivity, and electronic and optoelectronic properties, into unique physical properties. Especially, 2D materials offer an open canvas for the realization of heterostructures incorporating atomically and magnetically abrupt interfaces, giving rise to fascinating industrial applications.
To date, significant experimental efforts have been devoted to the synthesis and characterization of 2D materials, representing great progress. Further development of synthesis and characterization techniques is directly linked to our ability to manipulate and probe matter over a wide range of length scales. Therefore, the development of data-driven approaches that can shed light on atomic-scale events is of vital importance. Modeling and simulation can improve the atomic-level understanding of morphology and characteristics of 2D materials enhancing the ability to control the complex synthesis process and structural modulation. Quantum mechanics (QM) calculations are effective and accurate to predominantly explore the electronic, mechanical, and optical properties of materials; however, these methods are computationally expensive, limiting the system size to a few hundred atoms and the simulation time to a few nanometers. Empirical potentials, which can be classified into nonreactive and reactive potentials, have provided an alternative approach to simulating 2D materials at larger scales and with lower computational costs, not accessible by QM methods. Nonreactive potentials such as Lennard-Jones [13], AMBER [14], and CHARMM [15] provide computationally inexpensive methods for simulating molecular interactions in 2D-materials. However, the inability to simulate bond breaking/formation limits potential application of the nonreactive potentials. For example, the complex processes such as chemical vapor deposition (CVD) where gas-gas and gassolid reaction network is involved in the growth of a 2D material on a substrate cannot be investigated with use of the nonreactive potentials. Bond order based reactive force fields (i.e. Tersoff [16], Brenner [17] REBO [18] and ReaxFF [19,20]) can simulate the full dynamic evolution of chemical systems. Thus, the bond order based reactive methods bridge the gap between QM methods and nonreactive potentials. Among them, the ReaxFF has been the most extensively used to study the growth mechanism of 2D materials and their exotic mechanical, thermal, optical, and electrical properties that lead to various applications in nanoelectronics, insulators/semiconductor, catalysis, sensors, membrane separation, and other domains and will be further discussed in this review.
Recently, several theoretical reviews have been dedicated to provide a perspective on the modeling, simulations, and machine learning (ML) studies of 2D materials by various computational methods [21][22][23]. In this work, we particularly focus on a comprehensive review of recent ReaxFF force field developments and applications in the 2D materials' area, which is also expected to serve as a practical guide to users who are interested in studying the computational synthesis and processing of 2D materials at the ReaxFF level. As seen from table of contents above, in sections 2 and 3, we first briefly discuss the ReaxFF force field method, and compare traditional and ML techniques for optimizing ReaxFF parameters. We further review atomistic simulation techniques generally used with ReaxFF and other alternative approaches in modeling 2D materials, such as density functional theory (DFT), computational fluid dynamics (CFD), and phase field method. We dedicate sections 4 and 5 to layered (van der Waals) and non-layered materials, respectively, where each subsection focuses on a particular 2D material described at the ReaxFF level, as summarized in figure 1. This figure not only highlights the scope of this review but also summarizes recent progress in developing ReaxFF force fields for 2D materials and their applications. Then, in section 6 we focus the ReaxFF force field development and simulation efforts for emerging mixed dimensional architectures which are assembled by the combinations of 2D materials with other nD-materials where n is 0, 1, 2, 3. Sections 7 and 8 highlight recent progresses on the multiscale approaches that bridge spatial scales from nano to micro scales by combining ReaxFF, DFT, CFD, and phase field methods, and ReaxFF and ML combined hybrid approaches. Lastly, we provide our perspectives on the knowledge gaps and future directions in this field.

ReaxFF reactive force field
The ReaxFF reactive force field method developed by van Duin et al for hydrocarbons [19] is a bondorder dependent potential where charge on each atom is determined based on the polarizable charge method [47], potentially suitable to model any type of materials. Bond order is a measurement of the number of electrons involved in the bond between a pair of atoms. ReaxFF describes a functional relationship between the interatomic distance and bond order as well as bond order and potential energy. As the potential energy is a function of bond order, bond breaking and formation is possible, contrary to the case for any nonreactive potentials such as AMBER [14] and CHARMM [15] where bond connectivity is assumed as fixed. Therefore, ReaxFF simulations are highly effective in capturing detailed chemical events, reaction pathways, and product formation during the gas phase and gas/condensed phase simulations (e.g. CVD/MOCVD growth), for systems up to ∼1000 000 atoms and for timescales not accessible by QM-based techniques. The ReaxFF potential is an effective means of simulating the nucleation and growth of 2D materials, which involves gas phase and surface interactions affected by the local chemical environment, helping in predicting effective growth protocols [24,25,39,40,[48][49][50][51]. It features van der Waals interactions which particularly enable the simulation of multilayer van der Waals heteroand homostructures [24,25,39,48,50], besides nonlayered 2D materials [39,40]. Additionally, ReaxFF provides a thermodynamic and kinetic insight into fundamental solid-phase phenomena observed in 2D materials, such as phase transformation, grain boundaries (GBs), defect formation and diffusion, stress induced lattice distortions, and morphological evolution of 2D domains [48,[52][53][54][55]. The similar phenomena characteristic for any other non-2D systems could be investigated with use of ReaxFF, as well. Lastly, ReaxFF differs from the 'first generation' reactive force fields such as Tersoff [16] and Brenner [17], by applying a significantly longerranged bond-order relationship, which makes it possible to achieve accurate reaction kinetics. The functional form of potential energy equation, E system , deployed by ReaxFF can be defined as follows: (1) where the bond-order dependent E bond , E angle , and E tors denote the energy terms related to bond formation and dissociation, three-body valence angle strain, and four-body torsional angle strain, respectively. E Coulomb is electrostatic and E vdW is dispersive contributions. The energy term, E over, is used to prevent the over coordination of atoms. The calculation of each term contributing to the total energy uses an equation that considers the involvement of several force field parameters. The reader can find a more detailed form of equation (1) in [20].

Force field parameterization
ReaxFF force fields, like other empirical potentials, are system-specific and require tuning their parameters. Conceptually, the optimization of force field parameters can be divided into three main parts: (i) preparing a training set that consists of QM/experimentbased reference data. The quality and accuracy of a force field critically depends on the richness of a training set, meaning that physical and/or chemical phenomena that are expected to be reproduced by a force field should be adequately accounted for in the training data. This set is mainly composed of periodic data (e.g. the heat of formation, equations of state, structural impurities, diffusion kinetics, interface energies, and other properties) and non-periodic data (e.g. chemical reactions and energies, bond dissociation, angle distortion, and charges). (ii) Running an optimization algorithm that searches for optimal parameters minimizing the energy difference between reference and ReaxFF values. (iii) The validation of the developed parameters against test data by atomic simulations to ensure the reliability and accuracy of a force field. If a developed force field fails to produce target properties in atomistic simulations, more relevant QM/experiment data is added, and the force field training procedure is repeated. This cycle continues until every difference between the target QM/experimental and ReaxFF based values is not higher than assigned error and the expected properties can be obtained based on atomic simulations. The reliability and accuracy of a force field are critical to atomic simulations, and directly impacts the results generated over the course of simulations. Additionally, the ReaxFF force field includes a large number of bonded and non-bonded parameters for each atom. A presence of these many parameters significantly increases dimensions of the parameter space (PS) and degrees of freedom. Therefore, to avoid overfitting of parameters, a large training data set that is sensitive to be trained parameters is recommended. Due to these challenges, parametrizing the ReaxFF force field is a non-trivial task, and finding a local or global minimum in a high-dimensional PS requires a well-developed optimization algorithm. In the following section, we will discuss the common optimization algorithms utilized during the ReaxFF force field parametrization.

Conventional single parameter parabolic search
The single-parameter parabolic-search optimization scheme [56] is a brute-force optimization that is simple to implement but requires considerable computational resources. This method goes through all parameters one by one, iteratively through exhaustion, until the total error compared to reference data is below a set threshold. During the optimization, the error is computed based on the following equation where x i,ref is the reference and x i,calc is the ReaxFF predicted value. σ i (i > 0) is the weighting coefficient that reflects the significance of data of interest to optimization and target applications. Therefore, these weighting coefficients significantly impact optimization results. Optimal parameters are determined based on the best reproduction of training geometries obtained in successive iterations. This method first guesses a parabolic relationship between the total error and the value of a parameter of interest. Then, it determines the optimal value of a single parameter (i) by computing the total errors at three different values of a single parameter and (ii) by defining a new parabola that relates total errors to the updated value of a single parameter. This search process is repeated until the total error reaches a plateau. It is noteworthy that, during the optimization, one has to keep an eye on the optimization process because a parameter in search may reach an unrealistic and non-physical value, which needs the intervention of a developer.

ML based force field parameterization
A single-parameter parabolic search method inherits obstacles in determining global minima within the potential energy surfaces. Because of the highdimensional nature of PS, this method may cause a force field parameter to be trapped in an undesired local minimum reproducing some of the training data precisely while poorly fitting the rest, also called 'the curse of dimensionality' . This requires careful attention and skill of a force field developer during the optimization. Another drawback is that this method generates one local minimum because of the one-by-one parameter search, preventing the parallelization of this algorithm. This prompts the development of ML based optimization algorithms such as INitial-DEsign Enhanced Deep learning-based OPTimization (INDEEDopt) [57], JAX [58], iMOGA [59], etc, with a purpose of overcoming the issues caused by particularly high dimensional nature of the ReaxFF PS.
INDEEDopt is one of the ML-based algorithms recently developed by Sengul et al [57] and utilizes a deep learning-based optimization algorithm during the force field parameters' training (figures 2(a)-(c)). INDEEDopt produces improved accuracies in a shorter development time compared to the conventional optimization method. It is capable of finding several local minima in PS and produces multiple optimized force fields with low discrepancies that allow the precise fitting of whole reference properties defined in a force field training set, contrary to the conventional method [57]. Figure 2( Kaymak et al [58] developed another ML-based optimization algorithm, JAX-ReaxFF to speed up the optimization process of ReaxFF parameters. Implementing the auto differentiation tool of JAX into the ReaxFF force field training automates the calculation of the gradient of an error function in equation (2). This gradient-based optimization technique automatically generates forces from a given atomic position, thus, significantly reduces the computational time required for the force calculations. JAX-ReaxFF combined with a MultiStart approach allows the generation of several local minima, hence improving the efficiency of optimization and enabling its parallelization on architectures such as Graphical Processing Unit and Tensor Processing Unit. The Python-based JAX-ReaxFF code offers a computational platform for experimenting with the functional forms of the ReaxFF interactions or adding/removing interactions as desired. Figure 2(e) provides an overview of the JAX-ReaxFF workflow. In addition to single-parameter and ML-based methods, researchers also explored other optimization techniques for training ReaxFF parameters. For example, the Monte-Carlo (MC) simulated annealing approach enables to randomly scan high dimensional potential energy surface to determine global optimum parameters [60]. Still this approach depends on the initial conditions and is systemspecific, similar to the conventional approach. Guo et al [61] developed a gradient-based optimization method using ML, called intelligent-ReaxFF. This method allows fitting all parameters simultaneously, contrary to the sequential search of the conventional way; thus the search does not get trapped in a nondesirable local energy minimum. However, one of the drawbacks is that this method is sensitive to the training data size, and obtaining high-accuracy force fields becomes troublesome as the data size increases. Multiobjective optimization methods such as genetic algorithm [62][63][64][65][66] and MC simulated annealing [59] are among the other optimization algorithms adopted during ReaxFF parameterization.

Atomic scale simulations 2.2.1. Molecular dynamics
Molecular dynamics (MD) is a computational simulation technique initially developed in the 1950s with an attempt to solve the classical many-body problem numerically beyond the limits of conventional theoretical methods (i.e. lattice dynamics) [67]. By analyzing the motion of individual particles governed by empirical interatomic potentials, MD can reconstruct the collective dynamics of particles to explain, predict and design observable phenomena pertaining to the entire system [68][69][70]. Quantum behaviors can be incorporated in the interatomic potentials as semiclassical corrections [71][72][73].
MD simulations with the ReaxFF method provide fundamental insight into the discovering new reaction pathways and intermediates and predicting product compositions [23,74,75]. However, the short time scale, i.e. up to nanosecond, of ReaxFF MD simulations limits their applications to processes beyond the microsecond time scale. Therefore, previous studies have employed higher temperatures (>1800 K) to simulate the temperature-dependent processes, e.g. 2D materials growth from gas phase precursors in CVD systems [40,49,76,77], pyrolysis of 2D materials for hybrid materials design [78][79][80], and oxidation resistance tests of 2D materials at elevated temperatures [26,81,82]. The evident drawback of this approach is that entropically favored reactions are more likely to be picked [83][84][85], whereas it becomes a challenge when mapping a biased pathway obtained from high temperature simulations to an experimentally relevant process carried out at lower temperatures. The problems can be addressed using the accelerated methods in ReaxFF MD environments as discussed in the following section.

Common accelerated methods
The bias potential-based methods [86] such as metadynamics [87,88], hyperdynamics [89], and collective variable (CV) hyperdynamics (CVHD) [90,91] operate within the framework of free energy landscape by applying a bias potential to potential energy surface, filling energy minima and consequently lowering energy barrier for transitions. On the other hand, parallelization-based methods [92] such as parallel replica dynamics (PRD) [93] and parallel replica tempering [94], replicate the configurations with unique random velocities over multiple processors, such that each trajectory samples the basin of potential energy until a transition is inspected, hence providing a speedup proportional to the number of replicas. The combination of CVHD and PRD is identified to provide additional speedup to CVHD from milliseconds to seconds in a reasonable wall-clock time [95]. As an alternative approach, force-biased MC (fbMC) [96][97][98] and Grand Canonical MC (GCMC/MD) [99] utilize MC techniques, allowing for the system to reach the global minima faster.

CV-driven hyperdynamics
The CVHD which can be interpreted as a self-learning implementation of hyperdynamics, also incorporates the CV based feature of metadynamics, where a suitable history-dependent bias potential is slowly grown on-the-fly to trace the long time-scale evolution of a system. CVHD has been recently used in ReaxFF environments, particularly for the investigations of combustion mechanisms. Bal and Neyts [91] simulated oxidation of hydrocarbons at temperatures as low as 700 K. Ashraf et al [100] performed the pyrolysis simulations for single JP-10 (exo-tetrahydrodicyclopentadiene, C 10 H 16 ) and ndodecane molecules at an experimentally accessible temperature of 1200 K and density of 0.4 g cm −3 . Arash et al [101] investigated thermal degradation of amorphous polyamide 6,6 at low temperatures down to 1000 K and over a long timescale up to 0.17 µs. As there is growing interest in expanding the application of CVHD to ReaxFF MD simulations, we anticipate that CVHD could be extensively used to model the synthesis of 2D materials, such as the MOCVD growth of TMDs [102], in experimentally accessible conditions.

Parallel replica dynamics
PRD is a method that takes advantage of parallelism [92,93,103]. Unlike the other accelerated methods that are derived from transition state theory and are suitable for the systems with infrequent events, PRD requires no advanced knowledge of the available pathways. This is extremely useful to investigate system behaviors of high complexity. PRD as well as parallel replica hyperdynamics has recently experienced a coming of age, due to the increasing accessibility of parallel computing resources and the development of reactive potentials for MD. A number of studies that incorporate PRD to ReaxFF environments have been identified to extend the MD timescale and length scale effectively. Zamora et al [104] applied the extended Lagrangian Born-Oppenheimer MD in conjunction with PRD to study the shock compression behavior of liquid benzene on the order of nanoseconds. Joshi et al [93] employed PRD MD simulations for pyrolysis of a 1-Heptene system containing 40 molecules at temperatures as low as 1350 K for up to 1 µs. Hossain et al [105] explored the Li-electrolyte solvation, solvent exchange, and subsequent solvent decomposition reactions at the anode/electrolyte interface at room temperature over a couple of nanoseconds. Ganeshan et al [95] combined PRD with CVHD to analyze the speedup obtained during the pyrolysis of n-dodecane at low temperatures and timescales on the order of 2 µs.

Force-biased MC (fbMC)
As discussed in the previous sections, a class of methods were developed with the aim of extending the timescale of MD [106] and applied to ReaxFF environments. However, these methods are limited to infrequent event systems which could evolve through infrequent transitions from one metastable state to another [107]. A more general approach to achieve longer timescale is the coupling of MD with stochastic MC [108], where an MD cycle is employed to simulate fast processes, yet the following MC step is used for thermal relaxation at a longer timescale. In order to increase the acceptance rate of MC in a strongly interacting system, the fbMC was developed [109,110], through introducing deterministic forces into MC. The fbMC was later turned into a uniform acceptance version, i.e. uniform-acceptance fbMC (ufMC) [111], in which each MC step is accepted with unit probability [96] and the maximum allowed displacement needs to be sufficiently small [112].
The ufMC has shown great potential in revealing growth mechanisms of 2D materials in highly reactive environments. Neyts et al [27] performed the reactive ufMC simulations to investigate the growth of single layer graphene on a Ni (100) film in the temperature range of 300-1500 K, in far-from-equilibrium high precursor flux conditions. The nearly continuous graphene layers were obtained at 900 K and above. Sang et al [35] explored the Frank-van der Merwe growth of 2D h-TiC through the combination of experiments (in situ STEM) with DFT and ReaxFF MD using ufMC, reporting clustering of surface Ti and C adatoms into crystalline triangular islands. Lee et al [113] elucidated the nucleation and coalescence dynamics and 2D metallic growth of Ag adatoms on stoichiometric ZnO, O-deficient ZnO, and O-excessive ZnO surfaces from ufMC with ReaxFF. Jacquelín et al [114] employed ufMC with ReaxFF to study the dynamics of network structures of trimesic acid molecules on graphene at various temperatures, in response to different stabilities for honeycomb, filled honeycomb, flower, zigzag (ZZ), and closed packed 2D motifs.

Other computational approaches: from nanoscale to microscale
In addition to the atomistic simulation techniques used with ReaxFF methods discussed in the previous section, the following computational approaches can also be listed among the methods that model and simulate 2D materials.

Non-reactive reactive force fields
Empirical potentials are mainly categorized as reactive and nonreactive force fields based on the type of bond-connectivity adopted in simulations. Classical nonreactive potentials such as Lennard-Jones [13], AMBER [14] and CHARMM [15] assume fixed bond connectivity, meaning that bond rupture and reformation leading to chemical reactions does not account for during simulations. This allows faster simulations, thus enabling to model computationally demanding tasks. Although nonreactive force fields are useful to study the mechanical properties of 2D materials, such as elastic and thermal properties, they are not capable of simulating complex growth processes of 2D materials including different phases of materials, defect dynamics, and surface interactions, due to the static bond assumption. These potentials are system-specific like reactive potentials, therefore it is critical to choose one that can capture target interatomic interactions in a system. For example, Lennard-Jones and coulomb potentials describe nonbonded interactions, in other words, not-covalently bonded interactions such as van der Waals and electrostatic interactions in a molecule. Harmonic bond potentials such as AMBER are often used to describe the structural properties of a system containing covalent interactions.

Ab initio methods
Ab initio calculation is based on principles of QM, also called as first-principles calculations that can be used to model the electronic, optical and structural properties of periodic and non-periodic systems. The only input parameters required by these calculations are well-defined physical constants such as charges and masses of electrons and initial atomic positions. These methods are accurate and have been extensively utilized to study fundamental solid state, liquid and gasphase phenomena observed during the synthesis and characterization of 2D materials. However, ab initio methods are computationally expensive which limits their applications to small system size (<100 atoms).

Computational fluid dynamics
CFD model depends on the fundamental conservation laws of mass, energy, and momentum to describe a macroscopic behavior of fluids such as gas and liquids. In this method fluid is assumed as a continuous medium, thus enabling to simulate the fluid behaviors such as transport properties and helping with designing and optimizing fluid systems such as turbines and pipelines. This method has been extensively utilized to model the vapor-phase and colloidal synthesis of 2D materials (i.e. CVD/MOCVD growth) by simulating the transport profiles of volatile precursor molecules in a system. CFD models can provide insights into the effects of different parameters, such as flow rate and temperature, on the synthesis process, as well as the formation of gradients in precursor concentration and temperature that can affect the growth and morphology of the resulting 2D material. However, the mathematical description of these models relies on empirical correlations and assumptions deduced from experimental measurements and/or atomic scale calculations. Therefore, they are often used in conjunction with atomistic simulations such as ReaxFF and DFT, and/or experiments to provide a more complete understanding of fluid behaviors [25,35,113,115] as discussed in section 7 in detail.

Phase field
Phase field method is a coarse-grained approach aiding to model and simulate the microstructural evolution of a complex system with multiple phases or interfaces at mesoscale. Phase field equations governed by thermodynamics laws describe the spatial and temporal evolution of a phase. Therefore, this method is an effective numerical model to predict material science phenomena observed in systems such as interfacial problems, phase transformation, and topological changes during synthesis, solidification, and recrystallization. Depending on target applications, multiscale approaches can be further developed to improve the accuracy of phase field simulations by establishing an empirical correlation between atomistic simulations and phase field equations. In such models, atomistic simulations such ReaxFF and DFT are useful to generate input data on interactions and energetics that drive the dynamic evolution of GBs, interfaces, and defects in a system, for phase field models [116][117][118] as discussed in detail in section 7.

van der Waals layered 2D materials
van der Waals layered materials which are assembled layer-by-layer via weak van der Waals forces in their bulk forms, can be fabricated by bottom-up (e.g. CVD/MOCVD, molecular beam epitaxy (MBE), physical vapor deposition (PVD)) and top-down (e.g. mechanical exfoliation) growth techniques. Especially, van der Waals epitaxy with bottom-up techniques provides an effective way to mitigate lattice and thermal mismatches during the growth of highly lattice-mismatched systems. To date, extensive experimental studies have been devoted to the synthesis and characterization of novel 2D materials, representing great progress [11,12,41,76,[119][120][121][122][123]. However, thermodynamic and kinetic level understanding of growth mechanism, e.g. catalytic combustion of transition metal (TM) hydrocarbons by chalcogen precursors, surface reactions leading to the van der Waals epitaxy growth of TMD materials [49,53,76,119], morphological evolution of growthfronts as a function of the local chemical environment is not yet fully understood. In addition, the atomic-level influence of structural impurities such as point defects, adatom, GBs, and ripplocations on the material properties still remains elusive. As a complement, the theoretical synthesis and characterization of materials is a fundamental tool to deepen our understanding of the physical and chemical properties of materials, thus ensuring the highest quality in the design and manufacturing processes of 2D materials. Over the last decade, ReaxFF force fields have been used for large-scale (>1000 atoms) and longtime (>1 ns) dynamics simulations for 2D materials. In this section, we review the ReaxFF force fields and their applications for 2D layered materials summarized in figure 3.

Force field development history
Graphene has exceptional mechanical, electronical, chemical, and thermal properties, and therefore, is a promising material for a range of applications such as energy storage, electronics, optoelectronics, electrochemical sensing and biosensing. The ReaxFF development history of carbon-related parameters has initiated with the study of complex reactions in hydrocarbons [19]. Subsequently, this initial parameter set of C/H/O-2001 was modified to correctly predict transition states and initial chemical events associated with high-temperature gas phase oxidation of hydrocarbons, called as C/H/O-2008 [20]. It is noteworthy that both parameter sets do not target to reproduce any of the graphite or graphene characteristics. For example, Jensen et al [124] investigated the mechanical properties for carbon based materials including graphene with the use of C/H/O-2008. Despite a comprehensive discussion regarding a suitable choice for a time step value, type of the thermostat or value for a thermostat coupling parameter for the simulation of a fracture mode, authors acknowledged that since no data related to the mechanical response of the condensed phase carbon was used to derive the C/H/O-2008 parameter  set, one should not expect for a force field to correctly describe the mechanical properties of graphene. Srinivasan et al [125] further improved the ReaxFF parameter set to correctly reproduce the equations of state for graphite and diamond, as well as the formation energies of defects in graphene and amorphous phases of fullerenes. Based on the presented results on fullerene fragmentation, the improved parameter set, C-2013, is expected to be applicable for simulations of coal pyrolysis or soot incandescence. Ashraf and van Duin [126] further modified the C-2013 parameter set, and published as the C/H/O-2016, to be applicable also for syngas combustion, especially the conversion of CO 2 from CO. These force field development efforts were followed by extension of the C/H/O-2016 parameter set to carbonnitrogen interactions via introducing triple bond stabilization for nitrogen molecule [127]. The extended parameter set, C/H/O/N-2019 [127], subsequently, was applied to investigate the carbonization simulations of C/H/O/N-based polymers such as polyacrylonitrile, poly-p-phenylene-2,6-benzobisoxazole [127][128][129], as well as the polymer blends [31]. Further modification to the carbon-carbon interaction was introduced by Damirchi et al [130] to correctly reproduce the structure of flattened carbon nanotubes (CNTs), with a dumbbell-like shape characteristic of the curved part of CNTs. A summary of evolution of the carbon-carbon parameters to model carbonbased material is given in table 1.

Synthesis techniques 4.1.2.1. CVD growth
A number of synthetic techniques such as mechanical exfoliation [1,2], high temperature sublimation of SiC [131][132][133], chemical reduction of graphite oxide [134,135], and CVD on TMs [136][137][138][139][140][141][142] have been employed to produce high quality graphene. Among these methods, CVD is considered the most costeffective and promising approach to synthesize high quality and large area single or multilayer graphene [143,144] due to the tunable parameters in production. Although experimental research has long been dedicated to reveal the mechanisms of graphene growth, the lack of fundamental understanding of atomic details hinders its advancement. In recent N to allow for stable N2 production fCNTs-2020 [130] To correctly reproduce the structure of flattened carbon nanotubes (CNTs), with a dumbbell-like shape characteristic of the curved part of the CNTs. Optimized to represent graphite, diamond, defects in graphene and stable N2 production C to reproduce a dumbbell-like shape for the flatten CNTs decades, theoretical efforts, with the use of ReaxFF, have been mounted to address such a problem, so as to thoroughly understand the underlying phenomena and the reaction pathways of deposition of Cbased precursors, C nucleation, graphene growth, as well as the catalytic role of TMs. For example, the ReaxFF parameter set of C/H/Ni, developed for the chemisorption of hydrocarbon molecules into Ni (111), (100), and (110) surfaces by Mueller et al [145], has been widely used to simulate the CVD growth of graphene in the presence of Ni. Meng et al [146] studied the effects of C concentration and temperature on the kinetics of graphene growth on Ni (111) surface using elemental C precursors, from which the optimal defect annealing temperature of 1000 K was determined (figure 4(a)). Aside from the impact of precursor concentration and temperature, the role of the surface orientation on the graphene growth was also investigated by Lu and Yang [147]. Particularly, for the simulations using naphthalene and fluorene as precursors, Ni (100) and Ni (211) were found to have stronger surface interaction with volatile molecules, leading to the formation of discontinuous carbon networks, whereas Ni (111) surface exhibited moderate surface interaction and was beneficial to the formation of larger graphene nanosheets ( figure 4(b)). Chen et al [148] identified the CVD pathway of the CH 4 decomposition and the subsequent evolution from C segments to graphene. The optimal deposition rate of CH 4 and annealing temperature were found to be 0.037 ps −1 and 1200 K in theory (figure 4(c)). The ReaxFF C/H/Ni force field [145] has also been used for the CVD growth of graphene in the absence of Ni, such as the GB driven graphene growth on bi-crystal diamond [149] (figure 4(d)) and the graphene growth in a plasma-enhanced CVD environment sourcing from CH 3 , CH 2 , CH, and C radicals [150] (figure 5(a)).

Gliding arc plasma
Further development of ReaxFF force fields has paved the way for modeling the synthesis of graphene with novel methods and from new source of materials. For instance, Zhong and Hong [151] reported an ecofriendly and rapid method for preparing few-layer graphene through alternative-current rotating gliding arc plasma. The transformation mechanism of H 2 /CH 4 to graphene at atmospheric pressure was also investigated (figure 5(b)) by using ReaxFF including a London dispersion term developed by Liu et al [152]. Zhang and van Duin [77] developed ReaxFF C/Si/H-2020 force field to study the growth mechanisms of epitaxial graphene (EG) on C-and Si-terminated surfaces of SiC with C 2 H 2 and fluorene precursors (figure 5(c)).

Figure 5.
Representative ReaxFF studies of graphene growth (a) in a plasma-enhanced CVD environment, reproduced from [150] with permission from the Royal Society of Chemistry, (b) sourcing from H2/CH4, Reprinted from [151], © 2020 Elsevier B.V. All rights reserved. Copyright 2011 American Chemical Society, (c) on SiC surface, Reprinted with permission from [77]. Copyright (2020) American Chemical Society, and (d) using a laser-induced method. Reprinted with permission from [153]. Copyright (2020) American Chemical Society. this parameter set was also applied to study the changes in atomistic characteristics for the laserinduced graphene formation observed in experiment [153] (figure 5(d)). Moreover, with the use of the C/H/O/N-2019 parameter set [127], the initial edge growth of graphene during carbonization process was reported by Gao et al [154] and Rajabpour et al [155].

Mechanical properties
ReaxFF is a powerful tool for modeling mechanical properties, such as Young's modulus or fracture mode, due to its ability to treat covalent bond cleavage for any considered atom. For example, Jensen et al [156] tested the ability of the C-2013 parameter set [125] to predict the mechanical properties for several allotropes of carbon in the elastic regime as well as fracture mode. The reported simulation data indicates a better agreement between the DFT data and the one obtained with use of C-2013 parameter set rather than with use of C/H/O-2008 [20] (figures 6(a) and (b)). In 2021, Liu et al [152] also published a comparative study of a few bond order-based potentials suitable for carbon allotrope simulations. The results for lattice constant, thermal stability, cohesive and defect formation energies, mechanical properties as well as van der Waals interactions were taken into account in this assessment. However, the authors considered only C-2013 parameter set, not C/H/O-2016 [126] or fCNTs-2020 [130]. Nonetheless, the C-2013 parameter set placed as the second in overall accuracy measure, among all tested empirical force fields (figures 6(c) and (d)), with an accuracy of 85% for the adatom defect, 80% for cohesive/edge energies and thermal stability, 70% for the point defects and only 50% in predicting mechanical properties compared to the DFT target on average. Aghajamali et al [158] also compared the data obtained with use of the common carbon force fields, such as AIREBO [159] and Tersoff [16], including the C/H/O-2016 [126] and fCNTs-2020 [130], and identified these two ReaxFF force fields as the best to represent the interaction energies in layered carbon onions. Recently, Fthenakis et al [160] compared ten of the existing ReaxFF parameter sets to assess their applicability to correctly describe sp 2 carbon, and proposed a new parameter set with two different types of carbon atoms: one to describe sp 2 carbon and another one to represent carbon atoms in any other molecules. In their assessments, the authors considered the structural and mechanical properties as well as the energetics of graphene and CNTs. Based on this comparative study, the authors indicated that only C-2013 [125] and any other ReaxFF parameter sets published after 2013 (table 1) are suitable for studies of the elastic properties of materials with sp 2 carbon atoms. The ReaxFF method was also applied to investigate the mechanical response of a model of flattened CNT represented as layered graphene [161]. The presented results showed a significant impact of the irradiation-induced chemical crosslinking on the mechanical response despite a possible loss of the Reprinted from [163], Copyright © 2021 Elsevier Inc. All rights reserved. Schematic illustration of random distribution of the epoxide (-O-) functional group at the GB of bi-polycrystalline graphene in (g) the sandwich structure and with (h) the failure morphology with 50% epoxide coverage for the ZZ 21.8 • bi-crystalline graphene. Reprinted from [164], with the permission of AIP Publishing. structural integrity of the graphene layers due to the presence of vacancies.

Structural impurities
The quality and performance of as-grown graphene is sensitive to the growth conditions and precursors during CVD, which can be affected by the structural impurities such as point defects, dislocations, and GBs [162].

Point defect
Atomic-scale probing of the structural impurity induced thermal, mechanical, electrical, and chemical properties can provide valuable information about the coupling between the microstructure of graphene and its performance. Combined with studies on graphene growth [77] and as-grown microstructures [165], such investigations can contribute to the growth-performance mapping for the synthesis of graphene and other 2D materials. For example, Zheng et al [82] explored the defect driven oxygen adsorption to and the desorption of gasification products, i.e. CO, CO 2 , from graphene, respectively. As seen in figures 7(a) and (b), undercoordinated C surface atoms facilitate the O-functionalization of the surface by weakening the O=O bond. The adsorption of extra O atoms on the surface further induces the dissociation of CO 2 and/or CO through the formation of epoxy groups, resulting in the increase in the size of defect. In the CO 2 -rich and O 2 -poor conditions, the unpaired electrons on the C surface atoms can also catalyze the reduction of CO 2 to CO molecule through the O-functionalization of the surface.

Grain boundaries
Xu et al [166] also studied the mechanical properties of graphene oxide (GO) with GBs. They showed that GBs trigger the formation of defect clusters and provide zones for failures to propagate along. The degradation in tensile strength is obvious and can be worsened by defect coupling with functional groups with high oxygen content on GBs. Due to excellent mechanical, thermal, and electrical properties, graphene is also a promising reinforcing material added to nanocomposites to enhance relevant properties. Researchers have investigated how defects in graphene affect the performance of the reinforced nanoparticles. Verma et al [163] studied the interfacial properties of graphene/polyethylene (PE) based nanoparticles using ReaxFF and revealed that a bi-crystalline graphene containing GB outperforms the pristine graphene in enhancing the tensile strength, shear strength, and cohesive strength of reinforced PE nanoparticles. Figures 7(c)-(f) compares the shapes of pristine and bi-crystalline graphene in a graphene/PE system. Atoms at GB in bi-crystalline graphene are in higher energy states and induce stronger nonbonded interactions with PE molecules, bringing about wrinkling of the graphene sheet. The curved geometry further drives the local polarization and increases the number of adhesion points across the GB, which are the main mechanisms of improving tensile strength. The improvement of shear strength has been reported as due to higher adsorption capacity of defects sites to PE molecules at GB and improved π-π interactions at graphene/PE interface with the presence of sp 2 C atoms. Misorientation angles further amplify these effects and maximize the strengths. The thermal transport across graphene/PE interface was also studied by Verma et al [164] using ReaxFF combined with non-equilibrium MD. Bi-crystalline graphene with the presence of GBs yields higher thermal conductivity across the graphene/PE interface than a pristine graphene does. This is because high energy GBs can enhance the vibrational coupling between bicrystalline graphene and PE matrix, which improves phonon transport. These studies suggest that even though graphene with GBs has degraded strength and thermal conductivity, it serves as better reinforcing agent than a pristine graphene in polymer-based nanoparticles. On the other hand, the hydroxyl functional group only induces brittle failure for both of the GB functionalization case and sandwich structure. Figures 7(g) and (h) shows the sandwich structure of the system and its failure morphology. Since defects are almost inevitable, another worthwhile research direction is exploring the impact of the GB functionalization on the performances of graphene and its composites. Functional groups at GBs can alter the morphology, bonding, and stress distribution on nanosheets, and by consequence, the mechanical properties and transport behaviors for graphene.

Edge dynamics
Verma et al [167] studied the mechanical performance of ZZ and armchair (AC)-oriented bi-crystalline graphene containing different misorientation angles, and epoxide and hydroxyl functional groups. The overall tensile strength decreases with the increase in the number of epoxide groups at GBs, which is due to the hybridization transition from sp 2 to sp 3 . The offplane sp 3 bond imposes less restraint on bond rotation than the in-plane sp 2 and, thus, is more readily to distort and break. The functionalization of graphene with epoxy groups up to 50% coverage significantly improved ductility and toughness of functionalized GBs. This is because the sandwich structure introduces mixed failure modes, prolonging failure path and time.

Transitional metal dichalcogenides 4.2.1. Force field development history
2D TMDs with the chemical formula MX 2 are in a structural form of X-M-X where a hexagonally packed plane of M metal atoms (group 4-10B) is sandwiched between two planes of X chalcogen atoms (group 7A). Depending on the electron numbers in d orbitals of a TM, the TMD phase can be found in either semiconductor (i.e. 2H) or metallic (i.e. 1T, 1T ′ ) phase. The ReaxFF development history of TMDs has been started with MoS 2 by Ostadhossein et al [48] and continued with WSe 2 [49,50], WS 2 [24] and MoSe 2 [25] (figures 8(a) and (b)). The TMD force fields [24,25,48,50] that have been trained against the extensive QM data on both periodic and non-periodic systems describe the thermodynamic and structural properties of 2D single-layer and multi-layer TMD sheets. In particular, the force fields were shown to be able to accurately predict the formation energies and migration barriers of various vacancy defects, and the 1T → 2H phase transition of single-layer MX 2 . The energetic interplay between defects and sheet curvature was investigated by studying the formation energy and kinetics of a chalcogen vacancy migration along a moving ripplocation, suggesting a possibility to regulate defect distribution based on morphological changes in a 2D layer. Xuan et al [49] also developed a force field to model the MOCVD gas phase physiochemical processes of the crystal growth of TMD and combined it with CFD method for a multiscale modeling discussed in detailed in section 7 (figure 8(c)).
Hong et al [168] developed a force field for the CVD growth of MoS 2 by sulfidation of the MoO 3 surface deposited on Al 2 O 3 (figure 8(d)). This force field allows probing the early stages of the growth driven by the surface reactions of MoO 3 and S 2 from kinetic and thermodynamics aspects. It further predicts accurately the S 2 -assisted chemical reactions that lead to the formation of the byproducts O 2 , SO/SO 2 , the self-reduction and sulfidation of the MoO 3 surface, and finally, the Mo-S bond formation. Mishra et al [59] reparameterized the Hong et al's parameters [168] by adopting a dynamic approach based on the in-situ multiobjective genetic algorithm during the training. With this approach, the optimization of ReaxFF parameters of interest can be realized by fitting reactive MD trajectories of, in particular, farfrom-equilibrium chemical process (e.g. CVD growth of MoS 2 ) against MD trajectories on the fly. Mao et al [26] studied thermal degradation, oxidation, and hydrogenation of MoS 2 under O 2 and H 2 O gas phase exposures with and without Ti clusters as the compositing agent (figure 8(e)). This newly developed ReaxFF force field for Mo/Ti/Au/O/S/H interactions enables in-depth understanding of MoS 2 surface reaction mechanisms as well as the mitigating mechanisms of compositing agents on oxidation and hydrogenation. Ponomarev et al [169] developed another MoS 2 ReaxFF potential to improve the crystallization behavior of layered MoS 2 in melt-quench simulations (figure 8(f)). The training set consists of relative energies and structural features of both Mo x S y crystalline structures and amorphous MoS 2 models produced by melt-quench. This new force field accurately reproduces the structural parameters, the relative stability of Mo x S y phases and the energies of amorphous MoS 2 models. It further enables modeling the crystallization of MoS 2 from its amorphous form by melt-quencing method, outperforming the results from other state-of-the-art MoS 2 force fields (ReaxFF [48,76,168] and AIREBO [159]). Table 2  summarizes the force field development history for TMD materials.

Mechanical properties
MoS 2 is a typical 2D-TMD material, which has been extensively studied, due to its potential applications in electronics, catalysts, and solid lubricant, etc. ReaxFF, as an empirical interatomic force description, is especially suitable for systems involving van der Waals interactions and provides diverse research directions to probe TMD systems at atomic level. This method can illuminate mechanisms of dynamic processes of interest and provides guidelines for TMD growth and design. Islam and Haque [170] compared the mechanical performance of MoS 2 with an S-vacancy and O-doping, and discussed the deformation mechanism of MoS 2 containing GBs with different titling angles at the ReaxFF level. This research showed that both vacancy and O-functionalization influence the mechanical properties of MoS 2 -the increase in the density of the structural impurities in MoS 2 significantly weakens the ability of MoS 2 to resist the elastic deformation and improves the ductility of the materials (figures 9(a)-(f)).
Shi et al [171] showed that the tribological behavior of MoS 2 is closely related to the interlayer distance and temperature (figures 9(g) and (h)). Besides enhancing the thermal motion to compromise highsymmetry configurations, decreasing the interlayer distance can break high-symmetry configurations by deforming Mo-S bonds and deviating S atoms from equilibrium positions. This results in the increase in pressure and the friction force. Mao et al [26] also showed the MoS 2 surface is susceptible to the O 2 exposure which can degrade the tribological   properties of MoS 2 . Ti clusters were found to be an effective player in reducing the surface oxidation, by consequence, leading to the improved tribology behavior of MoS 2 . Rahman et al [172] explored the mechanical performance of MoS 2 affected by high temperature oxidation using ReaxFF ( figure 10(a)). The study reported the oxidation kinetics under different temperatures and oxidation mechanism that initates with the O 2 adsorption and is followed by the oxy-sulfide (MoS 2−x O x ) formation. The elastic behavior of MoS 2 is also contingent to the level of defect and surface oxidation which induces deformation in a monolayer and the associated mechanical properties such as fracture strength, fracture strain, and fracture toughness.

Synthesis techniques
ML assisted ReaxFF simulations [173,174] were conducted using a similar model with Hong et al [168], but with multimillion atoms, for the computational synthesis of MoS 2 on the MoO 3 surface deposited on Al 2 O 3 under the S 2 exposure. It was shown that the number of cooling/heating cycles on the CVD growth significantly influences the quality of the resultant structure, which offers a new way to control and repair structural impurities such as defects, GBs

Structural impurities 4.2.4.1. Defects and GBs
Sakib et al [54] reported the effect of GB defects with varying tilts on strain hardening of WSe 2 at different temperatures using Nayir et al's force field [50] (figures 11(a)-(d)). They showed that the fracture strength of WSe 2 depends on temperature, tilt and type of GBs. The increase in temperature and tilt yields decrease in the material toughness, inducing the plastic deformation due to the thermal activation of the bond breakage and dislocations blocked in GBs. Another work [43] compares the formation energies of various point defects predicted by Stilinger Weber [175], DFT and ReaxFF [50], showing that ReaxFF provides a much better agreement with DFT values, than those calculated by Stilinger Weber (figures 11(e) and (f)). Hickey et al [52] conducted a joint experiment and theory investigation on the underlying mechanism of invisible translational GBs in a WS 2 monolayer. Based on bright-field transmission electron microscopy (TEM) imaging and selected-area electron diffraction (SAED), the as-grown WS 2 films showed uniform and featureless nearly single-crystalline film across multimicrometer regions. However, dark-field TEM (DF-TEM) imaging uncovered unexpected, irregular, linear defects throughout the almost epitaxial monolayer (figures 11(g)-(i)). The ReaxFF MD simulations and scanning TEM (STEM) imaging revealed that, at GBs, two grains with inequivalent ZZ edges (chalcogen-terminated and metal-terminated) meet with geometry influenced both by the lattice offset and the dissimilar van der Waals interactions with the substrate during growth on a sapphire substrate. This result suggested that grains with metallic and chalcogen characters interact differently with electron beam, causing the out-of-plane distortion (figures 11(j) and (k)), as a consequence, the contrast differences observed in the DF-TEM image ( figure 11(h)). This work indicated the impact of the growth conditions and interface interactions on the crystal quality of epilayer. Another work [53] shed light on the underlying growth mechanism of WS 2 monolayer which contains linear arrays of W-vacancy defects. ReaxFF simulations revealed that the edge separation between two neighboring WS 2 grains drives the formation of W-vacancy defects. This work showed that the sterically hindrance of precursor/surface reactions may also induce the increase in the separation distance of grains and inhibit the complete coalescence of thin film, suggesting the deployment of precursors that are less willing to involve in sterically hindered reactions during growth.  [52]. Reprinted with permission from [52]. Copyright (2021) American Chemical Society.

Edge dynamics
Nayir et al [25] presented a theoretical approach that combines the kinetic Wulff construction (KWC) [176] and ReaxFF simulations to explore the morphological evolution of MoSe 2 during growth. The KWC model generated based on the energetics of various edge types with metal and chalcogen terminations showed that the two types of ZZ edges, called antenna and ordinary ZZ edges coexist in MoSe 2 domains during growth. The equilibrium morphology is expected to exchange between truncated triangular and hexagon shapes depending on the chalcogen chemical potential. The atomic force microscopy images further support this observation, showing the formation of truncated triangular MoSe 2 domains, which in some cases evolve to quasi hexagon shape.

Force field development history
MXenes (M n+1 X n T x ) [177,178] is a new class of 2D layered TM carbides and metal nitrides, which are synthesized by chemically etching bulk ceramics known as M n+1 AX n (MAX) phases, where M is an early TM (e.g. Ti, Zr, Hf), A is a group IIIA or IVA element and X is C or N with n = 1, 2 or 3. From etching process, the A element is removed from the MAX phase, leaving MX (MXene) layers bound by dispersion forces. MXenes exhibit a unique combination of properties, e.g. excellent conductivity, flexibility and volumetric capacitance, thus, can be candidates for a wide range of applications, e.g. lithium batteries [179], supercapacitors [180,181], sensing materials [182,183] or water purification [184]. ReaxFF force fields have been developed (table 3) and employed in the studies of various aspects of MXenes-from structural and mechanical properties to reactive chemical events on surface. Here, we highlight and discuss the recent progresses on MXene materials with an emphasis on the diverse applicability of ReaxFF.

Synthesis techniques 4.3.2.1. CVD growth
Bottom-up methods, e.g. CVD, allow 2D materials to grow on a large-scale and with high-quality compared to top-down methods (e.g. chemical exfoliation). Although there is a theoretical study [195] on the bottom-up synthesis of single-layer 2D TM carbides (TMCs) the assembly of ultra-thin layers (i.e. monolayer) have remained challenging in practice.
Since metal atoms of MXenes are provided from the metal substrate and bonded to gas atoms on the surface most MXenes prepared by CVD have non-layered structures with a thickness of tens of nanometers. In addition, residues and impurities are generated in the process of transferring the CVD sample attached to the substrate, and the production quality is significantly influenced by the growth conditions of the CVD ReaxFF parameter set Parameter optimization MX-2016 [36] The force field was generated by combining water force field [32] and Ti/O/H force field [185]. Ti/C/O parameters were optimized against Ti-C bond, C-Ti-C, and C-Ti-O angle distortions and the energy of condensed phases. MX-2016 [36] was employed in [35,37,81,178,[186][187][188][189][190] for ReaxFF simulations.
process [196]. Therefore, in order to develop largescale assembly methods for ultrathin TMC layers, it is required to search for suitable substrates and understand the interaction between the substrate and TMC layer.

E-beam induced growth
Sang et al [35] investigated the growth mechanism of a hexagonal-TiC (h-TiC) single adlayer synthesized on surfaces of Ti 3 C 2 substrates using in situ aberration corrected scanning transmission electronic microscopy. They described the process of homoepitaxial growth of a h-TiC single adlayer in two steps. The first step is that Ti and C atoms from the Ti 3 C 2 substrate were activated by electron beam irradiation and/or thermal exposure and migrated to the h-Ti surface, providing source materials for growth. In this process, pores were created on the substrate. The second step is the self-assembly of Ti and C atoms that have already migrated to the h-Ti surface, resulting in the island formation of a h-TiC adlayer. The crystal structure of the island is the new MXene Ti 4 C 3 and Ti 5 C 4 . Therefore, the growth of a h-TiC single adlayer was induced by a unique growth mechanism in which source atoms were provided through surface diffusion on the substrate. The dynamic growth process was observed using fbMC/ReaxFF MD simulations as shown in figures 12(a) and (b). DFT calculations suggested that the growth mechanism is related to the complex interplay between the low diffusion barrier and the high stepedge barrier. TiC dimers with low diffusion barrier migrate on the surface and clump together. This becomes the nucleation site for island formation and expansion. The step-edge barrier for the adatom to rise from the Ti 3 C 2 surface to the h-TiC adlayer through the edge is significantly large. Therefore, until the substrate is completely covered by the first h-TiC adlayer, a single layer grows as the adatom diffuses on the h-Ti surface rather than going up to the h-TiC adlayer.

Wafer-scale lateral self-assembly method
As another approach to fabricate ultrathin and largescale 2D materials using the bottom-up method, Mojtabavi et al [197] focused on controlling the material surroundings to be energetically favorable for the assembly process. They showed that the lateral self-assembly of MXenes at the interface between two immiscible solvents produces uniform mosaictype monolayers over large areas. To understand the interaction of MXenes with solvents, ReaxFF MD simulations were performed in various solvent environments. The results showed that the presence of methanol at the water/chloroform interface provides stability. Methanol molecules diffuse into MXenes and adsorb on the surface, forming hydrogen bonds. The formation of hydrogen bonds stabilizes MXenes even when solvated by chloroform. Therefore, the addition of methanol allows the liquid/liquid interface to be a favorable site for MXenes to self-assemble.

Mechanical properties of MXene
Plummer et al [186] studied the indentation of Ti n+1 C n T x MXenes using a ReaxFF interatomic potential to understand the effect of defects. Graphene has the highest tensile strength and was reported to have a Young's modulus of 1 TPa [187].
Recently, MXene has attracted great interest in applications such as energy storage, membranes and sensors that require robust materials that can withstand stress and pressure. The Young's modulus of MXene is 500-800 GPa, which is lower than that of graphene, but has advantages in processing. MXene can be produced in large quantities at a single time, whereas graphene can only be produced in small quantities from a singlelayer flake, and a typical graphene production process that requires mechanical exfoliation from graphite is inefficient. In addition, both MXene and graphene have hydrophilic surface that can be easily combined with polymer-based composites [198] after solutionprocessing. When graphene is solution processed, it results in GO with a greatly reduced Young's modulus of 200 GPa [199]. On the other hand, the Young's modulus of Ti 3 C 2 O 2 and Ti 2 CO 2 were calculated to be 466 GPa and 983 GPa, respectively, from the forcedisplacement curves, which are higher than that of GO. When Ti 3 C 2 O 2 contains 1% Ti vacancies and 10% C vacancies, the Young's modulus decreases, but still surpasses the Young's modulus of GO. Therefore, MXene is a promising material in terms of strength for solution-processed 2D materials. In addition, they presented that defects have a significant impact on fundamental mechanical behavior as an obstacle to cracks propagation [199].

Intercalation of ions in MXene
MXenes are one of the excellent platforms to be intercalated by diverse ions and molecules which can control the interlayer spacing of MXenes by occupying different positions with varying extents of hydration when intercalated between MXene layers. The dynamical response of Ti 3 C 2 T x MXene with different surface terminations to metal ion intercalants has been extensively studied with ReaxFF. Based on the hybrid GCMC/MD simulations [36,188], it was shown that cations can affect the ordering of the intercalated water due to the formation of the solvation shell. The ordering of the confined water increases with the ion size, with K+ > Li+. Prenger et al [193] investigated the deintercalation process of K + in MXene by performing ReaxFF MD simulations. Protonation of MXene surface and edges in H 2 SO 4 induces the removal of K + from an MXene interlayer. In addition, HSO − 4 ions promote the solvation of K + in the solution and its deintercalation. Reintercalation was found to be difficult due to the repulsive interaction between K + and the protonated MXene-edge.

Surface chemistry
Studies on the catalytic performance of MXene surfaces have been conducted using ReaxFF. As MXene is a redox active material, its surface can play a role as a cathode catalyst for oxygen-reduction and oxygenevolution reactions (ORR/OER) in the Li-O 2 batteries. For example, Ostadhossein et al [191] studied redox reactions on bare and O-and F-functionalized Ti 3 N 2 T x when the surface was exposed to Li/O 2 mixture. The ReaxFF MD simulations and first-principles calculations showed that the surface functionalization promotes OER reaction, whereas vacancies on the functionalized-MXene surfaces improve the ORR kinetics. Thus, the catalytic activity of surface was modulated by the surface terminations and defects of the MXene catalyst. Controlling the interlayer spacing of MXene is essential to achieve high capacitance. Urea intercalant between MXene layers can be readily decomposed on catalytically active MXene surface [37], yielding NH 4 + and CO 2 . Intercalation of NH 4 + as it remains between the MXene layers, leading to the interlayer expansion. Figures 12(c) and (d) show the catalytic effect of Ti 3 C 2 (OH) 2 MXene on the urea decomposition by comparing the rate of the NH 3 formation in three systems consisting of urea molecules (Urea), urea molecules immersed in water molecules (Urea_H 2 O), and urea molecules with a single Ti 3 C 2 T z MXene sheet immersed in water molecules (MXene_Urea_H 2 O) at the ReaxFF level. In ReaxFF, the energy barrier for the urea decomposition on MXene surface is 21.8 kcal mol −1 , which is lower than the experimental value for urea decomposition in water (32.4 kcal mol −1 ), implying that MXene surface can catalyze effectively in this reaction. Lastly, oxidation of Ti 3 C 2 MXene surface was used for generating hybrid carbon-supported nano titania structures. Lotfi et al [81] evaluated the effect of oxidizing agents and temperature on the oxidation of Ti 3 C 2 MXene by performing MD simulations at various temperatures under dry air (O 2 ), wet air (O 2 /H 2 O), and H 2 O 2 . It was clearly shown that as the oxidation proceeds under controlled temperature, Ti atoms migrate to the MXene surface, forming carbon-supported titania ( figure 12(e)). The crystal phase of titania was further investigated by Badawy et al [189]. In dry air, anatase and rutile TiO 2 were formed on Ti 3 C 2 T 2 (T=O, OH) MXene, whereas rock-salt TiO was formed under vacuum. In water vapor environment, the Ti-OH termination hinders the formation of continuous Ti-O bond, disrupting the formation of a distinct phase.

Nuclear quantum effects on proton transfer in MXene
To capture the nuclear quantum effects (NQEs) on proton transfer in the aqueous electrolytes confined between Ti 3 C 2 MXene layers, path integral MD (PIMD) in conjunction with ReaxFF [194] was adopted. PIMD is a method based on the Feynman path integral formulation that captures NQEs by sampling deviations from the classical path using several replicas. Although the structure and transport of intercalants were not severely affected by NQEs, it was found that NQEs enable the proton transfer from the surface OH group to confined water, which was classically not sampled in the high hydration system.

Friction of MXene
MXenes have gained great attention in tribology research as they have good physical and chemical properties such as low friction and good anticorrosion properties. A ReaxFF study [190] indicated that replacing a surface O atom with OH and OCH 3 in Ti 3 C 2 O 2 MXene significantly decreases the friction coefficients between the layers, compared to the Otermination. This can be explained by the strong van der Waals repulsion of hydrogen atoms between the layers. On the other hand, Ti-and O-vacancy defects increase the friction coefficients by increasing atomic-scale roughness and additional attractive forces between adjacent layers.

Hexagonal boron nitrides 4.4.1. Force field development history
2D h-BN (figure 3), with alternate B and N atoms in a honeycomb lattice, is isostructural and isoelectronic to graphene. It is an excellent electric insulator [200] while graphene is a zero-bandgap semimetal [1][2][3]. Various techniques have been exploited to synthesize 2D h-BN nanosheets. It can be achieved through mechanical exfoliation method, which is a simple technique of micromechanical cleavage of h-BN powders using adhesive tape [201]. Thin nanosheets obtained through this strategy process have high crystallinity and continuity over several microns; however, the thickness of h-BN layers is uncontrollable, the size is limited, and the production yield is relatively low [196]. Sonicating the h-BN powder in selective solvents (like dimethylformamide [202]) or water [203] can accelerate the dispersion and facilitate exfoliation in liquid phase (i.e. liquid exfoliation) [202][203][204], thus improving the production yield of h-BN nanosheets. Sputtering is the most used PVD methods for 2D h-BN nanosheets synthesis. By reactive magnetron sputtering of a solid B target in Ar/N 2 gas, a few-layer h-BN films with high crystallinity was achieved by Sutter et al [205]. This fabrication method paved the way for scalable synthesis of 2D layered h-BN for practical applications in devices. Epitaxial growth strategies, such as MBE [206], atomic layer epitaxy [207] and CVD epitaxy [208][209][210] can produce high quality h-BN nanosheets with excellent crystallinity and continuity, and controllable thickness. CVD epitaxy through thermal decomposition of gas phase borazine [208]/borazane [209,210] precursors on specific substrates such as Cu [209], Ni [208] and Pt [210] proved itself as a technique to fabricate thickness controllable, largesize and high-crystalline 2D h-BN layers. However, these synthesis methods are sensitive to the status of substrates; therefore, they are infeasible for large-scale fabrication of materials. This prompted the development of theoretical tools to have an in-depth understanding of the growth mechanism of h-BN layers. Benefiting from its capability and reliability to simulate the bond breakage and formation in chemical reactions, ReaxFF has been widely used in recent years to investigate the synthesis process for 2D nanomaterials as discussed in the following sections.

CVD nucleation and growth
Weismiller et al [215] developed a B/N ReaxFF force field for large-scale simulations of the ammonia borane dehydrogenation and combustion (table 4). Lele et al [211] recently optimized the ReaxFF force field developed by Weismiller et al [215] by fitting the parameters against DFT data for reaction energies and barriers for the formation of B-N species and hexagonal rings in the B-N gas-phase chemistry using B 2 H 6 /NH 3 as precursors. With this newly developed force field, Lele et al [211] further explored the CVD synthesis process of h-BN nanostructures from BN and HBNH precursors at elevated temperatures. They concluded that HBNH precursor forms H-terminated 2D h-BN structure with relatively small size and low quality and is less sensitive to temperature; while BN precursor evolves to a large cocoon-like BN nanostructure which converts to polymeric structures at elevated temperatures ( figure 13(a)).

Mechanical properties
Han et al [216] developed a force field for the interactions of hydrogen with single-walled BN nanotubes. Using this force field, Kumar et al [212] systematically explored the mechanical behaviors of partially and fully hydrogenated, and pristine h-BN nanosheets. The Young's modulus of pristine h-BN was predicted as ∼1.12 TPa which is higher than the experimental value, ∼0.865 TPa [217]. Kumar et al [212] showed that partially or fully hydrogenating h-BN lowers its Young's modulus down to ∼0.86-0.97 TPa but improves its integrity and stability which was demonstrated by smaller full width at half maximum values of the radial distribution function of those systems. The toughness of the partially hydrogenated h-BN with H atom on boron only was improved substantially ( figure 13(b)), which promises the potential application of h-BN for hydrogen storage. With the same ReaxFF force field [216], Kumar and Parashar [218] also predicted that the hydrogen passivation of the crack edge improves the fracture toughness up to 16%-23%, benefiting from the lowered crack-edge energy and blunting of the crack tip.

Defects, surface, and interfaces
The ReaxFF developed by Han et al [216] is also capable of capturing the interfacial interactions between h-BN nanosheets and other concerned materials.

Group III elements
Research interest in synthesizing group III elemental materials has rejuvenated over the past few years, primarily driven by research for atomically thin materials beyond graphene with unique properties [223]. Among group III elements, borophene [224,225] and gallenene [123], which are 2D analogue of graphene, have been experimentally synthesized to date, while the properties and applications of other group III elements can be derived by predictive modeling methods. Considering the lower number of valence electrons in group III elements, a number of elemental 2D-analogues have been theoretically predicted pertaining multiatomic unit cells [226]. Recent experimental measurement on gallenene confirmed such a predicted conformation, where two stable polymorphs in a buckled and planar 2D configurations are observed [123].

Borophene
The atomic structure of borophene ( figure 3) was, first, theoretically predicted by ab initio calculations [227]. Then, the first experimental demonstration of formation of borophene was achieved by Mannix et al [225] by growing borophene on an Ag substrate. Borophene attracted special attention due to its high strength along the AC direction and high thermal transport properties that can lead to their potential application to nanoscale heat transport. The buckled morphology in borophene opened doors for its use as energy storage material [228]. To date, various phases of borophene have been experimentally fabricated [224,225,229] in mainly hexagonal [230] and/or triangular [231] lattices in addition to the hydrogenated phase of borophene, called borophane [232]. The 2-Pmmn phase of borophene has risen significant interest due to its remarkable properties attributed to anisotropy in a material. The negative in-plane Poisson's ratio (∼−0.04 along the AC and ∼−0.02 along the ZZ direction) makes it an ideal candidate for flexible material. Low mass density, excellent tensile strength and intrinsic phonon mediated superconductivity makes it a possible candidate for flexible nanoelectronics and near-visible plasmonic applications. Borophene has a low vacancy formation energy and hence defects tend to form even at low temperatures during the synthesis of borophene [224]. Therefore, a wide range of studies have focused on understanding defect-induced thermal and mechanical properties. In 2010, Weismiller et al [215] developed a ReaxFF force field parameter set (table 4, [233] used this force field to study the mechanical behavior of borophene. Le et al [233] also compared the Weismiller et al [215] and Pai et al [221] force field parameters against the DFT data by Peng et al [234] and reported the Pai et al's parameters [221] to be more consistent with DFT. The Tersoff and AIREBO potentials were reported to be not suitable for the mechanical properties of borophene [235] while Stillinger-Weber (SW) potential [236] was used to simulate borophene's mechanical properties. Noroozi et al [237] used the Weismiller et al's force field [215] to study the defect-induced phase transition temperature, heat capacity, thermal expansion co-efficient and thermal conductivity of borophene, and the effect of chirality angle on these properties. At room temperature, the thermal conductivity was reported as 26.861 J (K·mol) −1 , in good agreement with an earlier DFT calculation (29 J (K·mol) −1 ) [238]. Their linear thermal expansion co-efficients along the AC and ZZ directions (−5.27 × 10 −6 K −1 and −2.01 × 10 −6 K −1 ) are also consistent with the previous DFT calculations (−7 × 10 −6 K −1 and −0.5 × 10 −6 K −1 ) [239]. The calculated melting point (2123 K) in their study is lower than that for borophene (2352 K) which was attributed to the size effect ( figure 14(a)). They reported a structural change (micro-phase transition) between 1525 and 1725 K as a sudden change in potential energy curve was observed, which was further confirmed by the mean square deviation from the initial structure ( figure 14(a)). They also investigated the effect of vacancy defects (0.5% and 1%) on the thermal conductivity and showed that the thermal conductivity is reduced by 36.5% with 0.5% and 50% on the presence of 1% defects. To study the effect of anisotropy, they reported the thermal conductivity along various chirality angles. The value of the thermal conductivity along ZZ is approximately half of the one along the AC direction. Mortazavi et al [240] used the same potential and reported the thermal conductivity of borophene to be approximately 76 W (m·K) −1 in the ZZ and 147 W (m·K) −1 in the AC direction. In 2017, Sadeghzadeh [235] used the Weismiller et al's force field parameters [215] to study the effect of chain-like boundary defects on the mechanical properties of borophene. Although the mechanical strength decreased with the increase in the defect size (line defect, figure 14(b)), periodic defects (figure 14(c)) with separation >6 Å did not affect the mechanical strength. In 2018, Sadeghzadeh and Khatibi [241] used the same force field parameters to study the vibrational frequency of graphene and borophene and concluded that borophene resonators are considerably more efficient than graphene resonators. Di Pierro et al [242] used classical MD simulations to calculate the thermal conductance between graphene/borophene and PDMS. Figures 14(d) and (e) show the graphene-PDMS and borophene (buckled)-PDMS systems, respectively. They used COMPASS force field [243] to define bonded and non-bonded interactions within PDMS. Tersoff potential [244] was used for graphene whereas the ReaxFF force field parameters developed by Weismiller et al [215] were used to define borophene. The heat transfer simulations resulted that the thermal conductance across the graphene-PDMS and borophene-PDMS interface are as 30 MWm −2 K −1 and 33 MWm −2 K −1 , respectively (figure 14(f)).
Arabha et al [245] used the Weismiller et al's force field [215] in non-equilibrium MD simulations to study the effect of defects/pores on the thermal and mechanical properties of borophene. They reported elastic modulus, toughness, and fracture strain by applying a constant strain rate of 1 × 10 9 s −1 on the borophene periodic structure with a box-length of 30 nm. While calculating the thermal properties, few layers of the structure were fixed to avoid slippage at the boundary. Their results showed the dependence of the thermal properties of pristine borophene on the length up to 60 nm, however, its mechanical properties became consistent after 30 nm. The 30% porosity (point defect) results in the reduction of more than 50% thermal conductivity in both directions. Upon increasing the porosity up to 30%, the elastic modulus, fracture strain and strength, and toughness decreases linearly. Also, the thermal conductivity and the elastic modulus anisotropy slightly decrease on adding layers (bilayer and trilayer). Mortazavi et al [240] used the Weismiller et al force field [215] to study the mechanical properties and the thermal conductivity of borophene and reported the Young's modulus of 188 N·m −2 and 403 N·m −2 and the thermal conductivity of 75.9 ± 5.0 W (m·K) −1 and 147 ± 7.3 W(m·K) −1 along the ZZ and AC directions, respectively. They reported an anisotropy ratio of 43% for thermal conductivity and 56% for elastic modulus of pristine borophene.
Dethan [44] used the Pai et al ReaxFF potential [221] to study the mechanical and thermal properties of borophane. Dethan reported the hexagonal borophene has similar conductivity in both AC and ZZ directions, which is contrary to the triangular borophene. The hydrogenation of borophene resulted in increase of the Young's modulus, especially along the AC direction. They also reported that a single vacancy defect only causes 0.11% and 0.09% decrease in the thermal conductivity in borophene and borophane, respectively.
Several studies have been devoted to investigate the mechanical properties of 2D layered group IV materials. For example, it was shown that the sp 2sp 3 hybridization and the buckled lattice structure lead to a much weaker mechanical strength for silicene, germanene and stanene compared with that for the sp 2 hybridized planar graphene [253]. Based on first-principles calculations, Mortazavi et al [254] also reported that the Young's modulus of silicene, germanene and stanene nanosheets along the AC and ZZ directions are 61.7 and 59 GPa, 44 and 43.4 GPa, 25.2 and 25.2 GPa, respectively, which are much lower than that for graphene. Roman and Cranford [255] explored the mechanical properties of silicene using the ReaxFF reactive force field Si/C/H/O-2012 [256] (table 5). Even though this force field was not specifically developed for silicene, it reproduces the lattice constant, buckling distance, and Si-Si bond length of silicene well comparable with the ab initio values [257]. With this force field, they obtained the in-plane stiffness of 50.44 and 62.31 N m −1 along the ZZ and AC directions, respectively, which are in  good agreement with those from DFT calculations [257]. The out-of-plane bending stiffness of silicene ( figure 15(a)) was calculated by rolling the 2D sheets to impose radii of curvature and was found to be ∼38 eV with R 2 value of 0.85, much larger than that for a monolayer graphene (2.1 eV) [258]. In contrast, Qian and Li [259] recently reported much smaller bending moduli for both monolayer and multilayer silicene ( figure 15(b)) (∼0.43 eV with R 2 value of 0.99 for monolayer silicene) than those for graphene [258] with the Si/O-2012 potential developed for oxygen interactions with silica surfaces [260]. The fracture patterns and edge reconstructions along the stretching process of silicene membranes were investigated by Botari et al [261] using Si/C/H/O-2012 [256]. They observed the continuous decrease of buckling for both AC and ZZ membranes but the edge reconstructions highly depend on the membrane type. For ZZ silicene membranes, only few pentagon and heptagon reconstructed rings were present once the rupture started (figures 15(c)-(e)) while for AC membranes, the edge was reconstructed significantly (figures 15(f)-(g)) [261]. Based on the MD simulation [262] with Si/O-2003 [263], the fracture strain of silicine was predicted to be further reduced when functionalized with hydrogen atoms. Baughman et al [264] predicted a new low-energy phase of carbon material, called graphyne, with planner sheet of the same symmetry as graphite, but part of the carbon-carbon bonds replaced by carboncarbon triple bond linkages. This new family of 2D carbon-based materials has some properties similar to graphene, such as the presence of the Dirac's cone or chemical stability, but this new carbon allotrope has uniformly distributed pores with tunable pores as well as a tunable electronic band gap. The schematic models of the graphyne-type structures are presented in figure 16(a). Raju et al [265] tested the desalination performance of a graphyne-based membrane using the C/H/O-2008 parameter set ( figure 16(b)), while Lin and Buehler [266] evaluated the selective water purification performance of the graphyne membranes as well as their mechanical response (figures 16(c) and (d)) using the same ReaxFF parameter set [20]. The C/H/O-2008 ReaxFF parameter set [20] was also used by Cranford and Buehler [267] to investigate also graphyne mechanical properties, and by Autreto and Galvao [268] to study the hydrogenation process of model graphyne membranes. Zhang et al [269] analyzed the thermal conductivity of an oxidized γ-graphyne via reverse nonequilibrium MD (figure 16(e)) based on the C/H/O-2008 ReaxFF parameter set, while and Paupitz et al [270] investigated an auxetic character of the considered graphyne-based membranes, with use of the C/H/O-2001 ReaxFF parameter set [19]. The uniqueness of the graphyne-type structures have resulted in significant emergence of fundamental research and a wide range of possible applications of this new material, e.g. in filtration/separation processes, energy storage, catalysis, nitrogen reduction, and also for photocatalysis or biomedical applications that are extensively discussed in the numerous recent reviews [271][272][273][274].

Group V elements 4.7.1. Phosphorene
Phosphorene can exist in many allotropic forms, α-P (black), β-P (blue), γ-P, δ-P, and θ-P. Black and blue phosphorene are the mostly studied and experimentally observed phases. Black phosphorene has a hexagonal crystal structure whereas blue phosphorene exists as buckled hexagonal sheets similar to silicene and germanene. Blue phosphorene was recently fabricated on an Au substrate by using black phosphorus as a precursor [277]. Black phosphorus is anisotropic whereas blue phosphorene is isotropic, but both show low band gap, making them a promising candidate for semiconductor application.
Moreover, the negative Poisson's ratio results in a flexible structure.
In 1982, a valence force field model was proposed to study the mechanical properties of black phosphorene, which was later used by Jiang and Zhou [278] to develop an SW potential. However, this potential underestimates the Young's modulus along the ZZ direction and the accuracy above near equilibrium region is not satisfactory. The ReaxFF force field parameters widely used for phosphorene were developed by Xiao et al [279] (table 6). They developed the P/H interatomic description to study the chemical and mechanical properties of phosphorene. Xiao et al [279] calculated the Young's modulus, Poisson's ratio, and adatom binding energy and compared them to the DFT and available SW potential-predicted values [278]. The results were in good agreement with DFT whereas the SW potential produced inconsistent results. Their ReaxFF force field predicted the temperature driven phase transition of AC/ZZ phosphorene nanotubes to faceted nanotubes. However, the SW potential was observed to underestimate the thermal stability of phosphorene nanotubes.
Xiao et al [279] further extended the P/H parameters to a preliminary P/H/O/C parameter set to study the oxidized phosphorenes and P-containing 2D materials. The P-O bond and several angles containing P/O/H/C atoms were trained for small molecules of P, H, and O. These parameters show good agreement with DFT data for binding energy of oxygen on phosphorene ( figure 17(a)). However, they did not study the oxidation of phosphorene in detail  and recommended further development of these parameters. Xiao et al [279] demonstrated the application of the developed P/C parameters by subjecting a graphene-phosphorene-graphene layered heterostructure to strain ( figure 17(b)). Chen et al [280] used this force field to study the mechanical performance of phosphorene in presence of defects (Stone-Wales) and elevated temperatures. Figures 17(c) and (d) show the α-P and its transition to β-P without and with Stone-Wales defect, respectively. In black phosphorene (α) with Stone-Wales-2 defect, a phase transition to mix β-P and γ-P was observed along the AC direction at high temperature under strain. The Stone-Wales defect was identified as a catalyst for this phase transition. Using the Xiao et al [279] force field, Le [281] studied the mechanical properties (Young's modulus and tensile strength) including the fracture mechanism of various allotropes of phosphorene at room temperature. Phosphorene was found to be more brittle (tensile strength < 25% of graphene) than other established 2D materials such as graphene, h-BN, SiC, AlN, and MoS 2 . Le [281] compared the ReaxFF optimized structural parameters and mechanical properties against various previously published DFT studies [282][283][284][285][286][287][288] for different allotropes of phosphorene. In 2016, Zhuo et al [283] used the SW potential developed by Dethan et al [44] to study the formation of fracture patterns along the AC and ZZ directions of phosphorene. In another work, Guan et al [284] used the Li et al [272] force field to study the mechanical properties of polycrystalline phosphorene. They reported an unusual linear elasticity characteristic to grain sizes ranging from 2 to 12 nm, where the elastic modulus increases by only 15.9%. Wei et al [285] studied the defect formation and evolution in black phosphorene under argon ion irradiation. Note that the phosphorus-phosphorus interactions were defined by Xiao et al ReaxFF parameters [279] and argon-phosphorus interactions were defined using the ZBL potential [286]. Figure 17(e) shows the computational model of the phosphorene used for irradiation. They reported that under low fluency (<2.86 × 10 13 ions cm −2 ), voids form in a coherent lattice, followed by healing (∼42% reduction in void area) upon annealing. This defect healing was attributed to re-organization of puckered P rings (figures 17(f)-(i)).

Bismuthene, arsenene, and antimonene
Theoretical calculations of bismuthene suggest it to be a good insulator. Although bismuthene has been fabricated on a substrate [290] a free standing bismuthene is yet to be experimentally observed. Thus far, the ReaxFF parameter descriptions for arsenic (As) and antimony (Sb) are not available and a force field for bismuthene is not developed. However, Bismuth's atomic parameters are available [291] which can be further developed to study bismuthene. Chowdhury et al [292] used the SW potential to study the mechanical properties and the phonon transport mechanisms in bismuthene. However, no force field based study is available for arsenene and antimonene to the best of our knowledge.

Non-layered 2D materials
Fabrication of atomically thin non-layered materials is quite challenging compared to van der Waals structures because inplane interactions within a nonlayered material are dominantly governed by chemical bonds which are difficult to break without sacrificing the crystal quality. Direct epitaxy of nonlayered materials (e.g. GaN, InN, Pb) on latticemismatched substrates is strain-driven and results in 3D islands as a consequence of the combined effect of the Volmer-Weber and Stranski-Krastanov growth modes. Alternatively, taking advantage of the Frankvan der Merwe growth mode, recent experimental studies offer novel techniques [11,12] such as CHeT that enables the growth of 2D-nonlayered materials atom by atom by encapsulating them at the half van der Waals interface of EG and substrate (e.g. silicon carbide (SiC)) where EG is used as a cap layer to stabilize intercalants at the interface.

Group III metals
Experimental identification of 2D materials has indeed made significant progress within the post graphene era; however, atomistic reaction pathways leading to the CVD growth of 2D group III elements cannot readily be obtained through experiment. First principles methods have been extensively used for predicting the structural characteristics [293,294] as well as the reaction kinetics of group III element growth in the presence of pristine or functionalized graphene [11,42,295,296]. Nonetheless, these DFT studies can only capture the reaction mechanisms of group III elements in small systems, and the need for a more computationally effective method, such as ReaxFF which allows for larger spatial scale and longer timescale MD simulations than the DFT calculations has been raised. Rajabpour et al [40] developed Ga/C/H-2020 and In/C/H-2020 ReaxFF parameter sets, aiming to reproduce the gas phase reactions of Ga and In having the 2D growth from trimethylgallium (TMGa) and trimethylindium (TMIn) precursors with a graphene substrate in the CVD environment. The force field was used to simulate the optimal processing conditions for the formation of Ga/In nanoclusters from TMGa/TMIn yet with low carbon contamination. In addition, Ga atoms growing on AC-edged graphene were found not only to exhibit a superior growth ratio to In atoms but also to produce a widely spread 2D thin layer between graphene edges ( figure 18(a)). Nayir et al [39] connected further ReaxFF and DFT based simulations to the experimental measurements done by Raman spectra, x-ray photoelectron spectroscopy, scanning tunnelling microscopy and spectroscopy that reveal defects in graphene which can act as pathways for Ga intercalation. The theoretical results showed that the binding strength of Ga to graphene is impacted by the size and parity (i.e. evenodd number) of defects, and the size of Ga clusters. Ga anchored at defects experiences a stronger interaction with graphene as the number of C-dangling bonds nearby defects is increased, but this comes at the expense of larger dissociation energy of Ga to intercalate from graphene at the SiC/graphene interface. This suggests that bare defects have a tendency to 'trap' Ga atoms in a graphene layer. Alternatively, as reported by Briggs et al [11], passivating C-dangling bond with functional group (e.g. epoxy, carbonyl, hydroxyl) via plasma treatment significantly weakens the Ga binding strength to the surface. Particularly, carbonyl and epoxy functional groups still serve as attractive centers but with lower dissociation energy, enabling Ga atoms to intercalate beneath graphene. However, the presence of hydroxyl groups can significantly lower the binding strength of Ga to the surface (even lower than the pristine case), indicating that hydroxylated graphene regions may inhibit the Ga binding to graphene surface as also reported by Bansal et al [41]. The Nayir et al's ReaxFF [39] and DFT [42] results further show that the sizes of vacancy defect and intercalant affect the thermodynamic and kinetic preference of intercalant between the adsorption on graphene surface or the intercalation at the interface. For example, the diffusion of a relatively larger-sized Ga atom through small sized defects (<DV) is kinetically hindered, where Ga encounters a kinetic barrier. But the potential energy of the system monotonically decreases as Ga approaches to graphene with larger sized defects (>DV) and the diffusion is nearly barrierless (figures 17(b)-(f)). This suggests that defect engineering provides an effective way to lower the temperature required for intercalation, by consequence, energy savings and cost reduction during the 2D material fabrication. A Ga atom intercalated at the interface is also found to thermodynamically more stable than adsorbed on the graphene surface, where the presence of SiC generates a lower local minimum for Ga to adsorb than the graphene surface (figures 17(g) and (h)), as reported in reference [42]. This signifies that the SiC substrate facilitates the Ga intercalation into the SiC/graphene gallery.

Heterostructures
2D materials with a dangling bond free surface provide a unique opportunity to fabricate high quality heterostructures with diverse multifunctionality (e.g. Janus [297], Z-scheme heterojunctions [298]) without concerning thermal and lattice mismatch between the abrupt interfaces of the joined components. 2D materials can be used further as contact for the fabrication of mixed dimensional devices which are the integration of materials with various dimensional ranges such as 0D/2D, 1D/2D, 2D/2D and 2D/3D architectures [299] ( figure 19). In this section, we will focus on the ReaxFF applications of mixed dimensional heterostructures.

Heterostructures containing different phases of one material
More rich and tunable TMD properties can be achieved by forming heterostructures. ReaxFF serves as an ideal investigation method for these systems due to the transferability of atomic interactions across interfaces. Mortazavi et al [300] investigated the mechanical performance of single-layer MoS 2 consisted of the semiconducting 2H and metallic 1T phases with varying domain sizes at room temperature (figures 20(a) and (b)). Such 2H/1T heterostructures exhibited reconciled elastic modulus, tensile strength and strain as a function of the concentration of 1T domains. Based on the fact that pristine semiconducting 2H phase has better mechanical performance than pristine metallic 1T phase, 2H/1T heterostructures simply exhibited negative correlations between the mechanical properties and the 1T phase concentration regardless of the domain size ( figure 20(d)). Except for the correlation that the elastic modulus deceases as 1T domain size decreases, the domain size did not affect the tensile strength and failure strain. Additionally, figure 20(c) shows that in all cases, lattice distortion initiates from the center of the 1T domains. The stress accumulates at triangular corners at first and propagates to the center, and then the deformation bands penetrate through the phase boundaries and propagate perpendicular to the loading direction across the entire sample. Such deformation behavior is different than that of single-phase polycrystalline MoS 2 , but still shares

0D/2D heterostructures
Zero-dimensional nanoparticles hold a great promise for next generation devices owing to their large surface area to volume ratio, super electrical conductivity and electrocatalytic activity. However, because of their high surface energy they are far from the equilibrium state and chemically highly reactive. They tend to form complex fluids by aggregating to each other in order to reach to energetically stable state. But, this also results in the decrease in their chemical activity. Alternatively, assembly of nanoparticles with 2D materials offers a unique opportunity to control the size and morphology of nanoparticles by inhibiting their aggregation into larger clusters and increasing their effective active surface, thus, enhancing catalytic activity and enabling the long-term storage for energy applications. However, the synthesis of such hybrid materials is quite challenging, requiring atomic-level precise control of each component joining in a heterostructure. This has prompted the development of multiscale frameworks that can enable engineering the morphology and interface of the 0D/2D composite. For example, ReaxFF force fields have been developed for various nanoparticles such as Ag, Pt, Cu nanoparticles supported on graphene. Zhang et al [301] used the Pt/C/O/H ReaxFF force field [302,303] to study the catalytic effect of a single Pt nanoparticle on the decomposition of PE (figures 21(a) and (b)). The presence of a Pt nanoparticle (with or without GO support) was shown to facile the PE depolymerization by lowering the reaction activation energy, compared to simple thermal degradation.
The catalytic behavior of a Pt nanoparticle also selectively yielded different distributions of product molecular mass, depending on the reaction temperature. This study highlights utilizing ReaxFF in catalyst design for plastic chemical recycling applications. Dabaghmanesh et al [304] used ReaxFF simulations to investigate the structure and morphology of Cr 2 O 3 nanoclusters over graphene and CNTs. The mobile Cr 2 O 3 nanoparticles were reported to selfassemble quickly and form worm-like large nanoclusters on graphene as well as on both inner and outer surfaces of CNTs due to weak van der Waals interactions between clusters and substrate ( figure  21(c)). Structural analysis of the clusters revealed that the short-range order resembles that of the bulk. In addition, the formation of large aggregation significantly deforms the graphene at low temperature, while increasing temperature tends to flatten the graphene surface. Such a finding enhances the basic understanding of decorating graphene and CNTs with magnetic Cr 2 O 3 nanoparticles with potential applications in spintronic devices. Tainter and Schatz [305] developed a ReaxFF force field to model the ZnO nanoparticle formation from diethyl zinc interacting with an epoxide-functionalized graphene surface. The mechanism of oxygen abstraction from the graphene surface was investigated using ALD simulations. Further, the condensation reactions that ultimately lead to the formation of a non-crystalline ZnO nanoparticle was modeled by an MC acceleration approach ( figure 21(d)). Time dependent DFT calculations on capped ZnO clusters were also conducted to reveal significant red-shift in the electronic absorption spectrum due to structural disorder. This study provides important insight in the formation mechanism and optoelectronic properties of ZnO grown by ALD technique. Lee et al [28] engineered silicon nanoparticles to improve the capacity and performance of lithium-ion batteries (LIBs). For this, they performed a comparative study of the lithiation process of a silicon nanoparticle and crystal systems using ReaxFF simulations. On the basis of their simulations, as the size of silicon nanoparticle approaches the bulk dimensions Si-bonding network becomes much weaker (figure 21(e)). This faciles the lithiation process of the silicon surface and thus improves the capacity of LIBs. Another study reported the correlation between the number of multilayered graphene and their ballistic penetration energy at the ReaxFF level [306]. In these simulations, Ni projectile was shot against graphene layers with varying thickness (from single to three layers) at various angles and velocities ( figure 21(f)). Their results showed that the penetration energy of graphene sheets decreases with the increase in the number of layers, indicating that with the increase in the number of layers, the multi-layered graphene sheets become less effective in absorbing the impact energy of projectiles per affected graphene mass.

1D/2D and 1D/3D heterostructures
Mixed dimensional 1D/2D and 1D/3D heterostructures offer new functionalities beyond the individual layer components. This has indeed prompted the synthesis and design of various combinations of heterostuctures [30,130,307,308]. 1D/3D heterostructures have been extensively studied at the ReaxFF level, with a particular interest in the nucleation and growth of CNTs on a 3D matrix such as Ni, Al [130,307,308]. For example, Khalilov et al [307] combined ReaxFF MD and time-stamped force-bias MC (tfMC) simulations to probe the nucleation and growth of single-walled CNTs (SWCNTs) on a Ni 55 catalyst from hydrocarbon precursors. Their results showed that the nucleation and growth depend on the absorption/desorption rates of hydrocarbon/H 2 molecules. The hydrocarbon degradation on the catalyst led to the surface carbonization, followed by the H 2 desorption. During the dehydrogenation of the surface, the carbon cap formation occurred via either the diffusion of the surface C atoms into a cluster in case of the full dehydrogenation or floating Catoms incorporated in a Ni catalyst toward the surface as a consequence of the partial dehydrogenation. Eventually these C atoms involved in a carbon network to form CNTs ( figure 22(a)). Another crucial result drawn from the simulations is that compared to the ReaxFF cycles to tfMC, ReaxFF potential utilized in this work [19,145] generated much reliable results for the nucleation and growth of CNTs from hydrocarbons such as the hydrocarbon degradation into carbon species and H 2 and the surface re/dehydrogenation. Another work [308] utilized the DFT calculations to benchmark the performances of empirical potentials (i.e. ReaxFF, AIREBO, EAMFF) in describing the nickel/CNT interactions by designing the pull-out simulations of Ni-decorated or coated CNTs on an Al matrix. The crucial highlights from this work are as follows: (i) On the basis of the ReaxFF simulations, the Ni coat significantly triggers the energy dissipation observed during the pull-out of CNT from the substrate. This may be the cause of the high pull-out force produced during the process, as a consequence, the rupture of CNT, showing a good agreement with earlier studies [309,310] (figures 22(b)-(e)).
Additionally, the Ni coating comes with benefits especially in increasing the resistance of the Al matrix to the plastic deformation by driving the recrystallization of Al. The comparative studies also showed that ReaxFF, developed particularly for C/Ni/Al interactions, offers an opportunity to investigate the complex surface interactions between CNT/Ni and Al/CNT, but it is noteworthy that this potential still needs improvement for better description of the mechanical properties of the Al matrix and of the Ni-Al alloy system. The aforementioned two investigations briefly exemplified the ReaxFF applications to 1D/3D heterostructures. However, the area of 1D/2D heterostructure at the ReaxFF level is still infant, there exist a few studies reported on 1D/2D heterostructure of same materials [30]. For instance, Zhang et al [30] modeled a SWCNT/carbon fiber to explore the catalytic impact of Ni and NiO nanoparticles on the growth of SWCNT on the carbon fiber surface ( figure 22(f)). Ni particles were found to hinder the growth of SWCNT on the surface by weakening the C-C bond interactions between the joint components in the heterojunction while the O-rich environment effectively saturates Ni nanoparticles, leading to the formation of NiO nanoparticles. Thus, this prevents Ni particles from weakening the C-C bonds formed between the SWCNT and carbon fiber. This result suggests the realization of the SWNT growth on the surface under the aerobic condition. It is noteworthy that although force field descriptions for the individual components of 1D/2D heterostructures of different materials may be available (e.g. the force fields for CNTs and MoS 2 exist but no force field description has been reported yet for the assembly), modeling and simulation of a mixture may need the extension of the individual force fields to the interface interactions.

2D/2D heterostructures
Two types of 2D/2D heterostructures exist: vertical heterostructures with different material layers stacked vertically, or lateral heterostructure with in-plane domains of different materials. The lattice continuity at lateral interfaces yield unique electronic properties, including outstanding current rectificaticon behavior and photocurrent generation characteristics [311,312]. On the other hand, vertical heterostructures are held together by van der Waals interactions, which may induce structural changes and charge redistribution between neighboring and even more distant layers. By combining the intrinsic merits of 2D components with the heterojunction effect, 2D/2D heterostructures have attracted significant interest over the recent years due to their unique properties and potential applications in electronics, photonics, and energy storage. By carefully engineering the composition and arrangement of layers, it is possible to create materials with enhanced properties such as high carrier mobility [312], tunable stiffness and strength [209], and improved catalytic activities [313]. For example, Berdiyorov et al [314] modeled a silicene/graphene heterostructure. As the crystallographic structures of silicene strongly depended on the atomic arrangement of substrates, silicene synthesized on metallic surfaces were structurally very different from the theoretically predicted buckled phase. Berdiyorov et al [314] proposed a novel method of using a weakly interacting surface, i.e. graphene bilayer, to stabilize silicene. From ReaxFF MD simulations, upon intercalating bilayer graphene, a flat hexagonal arrangement of Si atoms intercalating bilayer graphene transitioned to a buckled honeycomb structure, showing a close agreement with first principles predictions (figures 23(a)-(h)). Such quasi-2D structures have been found to be stable well above the room temperature, suggesting the feasibility of using graphene templates to grow pristine silicene.
Ganeshan et al [192] also evaluated the structure and dynamics of water in a heterostructure of 2D TiO 2 :Ti 3 C 2 T 2 (T=O, OH) MXene with the intercalation of alkali cations. According to the ReaxFF MD simulation results, 2D-TiO 2 decreased the water selfdiffusion in both cation-free and cation systems due to the water molecules entrapped on the notch oxygen sites. Among them, Li + ion hindered water diffusion the least, because Li + ion also occupied preferentially this oxygen site and prevented water molecules from being trapped, enabling water to be more mobile than the cation-free system. On the other hand, Na + ion hindered the water diffusion the most because of its position in the water layer similar to MXene layers ( figure 23(i)). In addition, 2D-TiO 2 has a more favorable adsorption surface for protons, similar to Li + , than the MXene surface and increased the accommodation of protons, compared to the MXene structure. Ibrahim et al [315] performed ReaxFF simulations on the water desalination performance of TiO 2 /MoS 2 -hexagonal and rhombohedral nanocomposites bilayer membranes. They revealed that adding one layer of TiO 2 onto the MoS 2 layer effectively enhanced both the water permeability and salt rejection percentage through the membrane ( figure 23(j)).  [192]. Copyright (2020) American Chemical Society. (j) Water permeability through membrane at different applied pressures at t = 4 ns and salt rejection percentage at different applied pressures at t = 2 ns. Reproduced from [315], Copyright © 2022, The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature.
Fan and Yao [316] applied ReaxFF to study the coupling effect of functional groups and GBs with different contents and mismatching angles, respectively, on the mechanical response of monolayer GO/h-BN heterostructures (figures 24(a) and (b)) by comparing the mechanical properties of graphene with GBs (Gr-GBs), graphene/h-BN with and without GBs (Gr-BN-GBs and Gr-BN), and GO/h-BN with GBs (GO-BN-GBs). This study revealed a reverse trend in the bending rigidity shown in figure 24(c), compared to the calculated Young's modulus , and considered it as an example of defected 2D materials which violates the classical plate theory. The combination of epoxy and hydroxyl groups improved the ductility of GO-BN-GBs compared to Gr-BN, as shown in figure 24(d). The coupling of hydroxyl and epoxy groups gave GO-BN-GBs moderate failure strength between the highest given by hydroxyl and the lowest given by epoxy for all cases. Increasing the density of functional groups decreased the failure strength for all cases (figure 24(e)).

2D/3D & 3D/2D heterostructures
Liu et al [222] developed a ReaxFF force field for B/N/Ni system to investigate the growth mechanism of h-BN monolayer from elemental B and N on a crystalline nickel substrate. The parameterization of Ni-B and Ni-N pair interactions is based on the first principles calculations carried out with DFT. The MD simulations based on this ReaxFF force field showed that B and N atoms initially form BN chains. Those chains evolved into longer branches with increase of B and N atoms, and then folded into hexagon rings serving as a nuclei for growing into larger h-BN monolayer ( figure 25(a)). They also concluded that growth of this atomic thin h-BN layer prefers higher temperature. With the well-developed ReaxFF force field for B/N/Ni system [222], Liu et al [34] further investigated the evolution of morphologies of individual h-BN domains on the crystalline Ni substrate with respect to the B:N ratio. They found that with elemental B and N as sources, excessive B species hinders the growth of h-BN domains while excessive N leads to larger domain (figures 25(b)-(i)). The results from MD simulations with ReaxFF are consistent with those from DFT calculations, which validate the capability and reliability of ReaxFF force field for investigating 2D h-BN growth with tuned chemical environments. Neyts et al [27] investigated the mechanisms of graphene growth on the Ni (100) surface under far-from-equilibrium high precursor flux conditions using hybrid reactive MD/ufMC simulations. In contrast to experimentally observed surface segregation mechanism at low flux conditions where the graphene formed through the segregation and precipitation of dissolved and diffused carbon atoms at low flux conditions, a combined deposition-segregation mechanism has been identified under high precursor flux conditions. The impinging carbon atoms that dissolve into the first and second subsurface layer became saturated, and the subsequent carbon atoms were mostly deposited on the surface.
In the simulations, nearly continuous graphene layers grew at temperatures of 900 K and above, suggesting a possibility to control the graphene growth mechanisms through adjusting the precursor flux and substrate temperature ( figure 25(j)). Other structural characteristics such as rippling of graphene sheets over catalysts have been explored with ReaxFF as well. Kowalik et al [29] validated the ReaxFF capability to simulate graphene bending and chemical interaction with catalyst metal surface. The research suggested that graphene is prone to ripple with external stress. Graphene geometries can be altered differently with different heating and annealing conditions. Comparisons between flat and rippled graphene on copper surface showed that the binding of copper atoms is more favored by rippled graphene. The examples showing the good agreement in describing graphene bending, ripple formation as well as binding characteristic for copper surface were recently reported [29,317]. With the use of previously published parameter set where C-Cu interactions were optimized for Cu-metal surface catalysis [318], the authors were able to reproduce experimentally measured graphene draping and rippling angle at the copper step edges and also show better binding properties of the rippled graphene.

Multiscale modeling of CVD growth of 2D materials
The elucidation of nucleation density, morphology, size, and orientation of growth fronts on a preferred substrate and their coalescence into a thin film is critical to improve the crystallinity of materials. Multiscale models are powerful tools to connect macroscopic growth parameters such as temperature, flow rate and inlet velocity to the morphological properties of 2D domains. Such models approaching actual experimental dimensions provide a direct way to communicate with experiments and can test experimental conditions. Thus, they provide insights into more effective protocols, shedding light on the significance of reactor geometry, and improving the reproducibility of reactor results across different experimental growth chambers.
For instance, the following two studies prove the power of multiscale models in developing and optimizing the growth parameters of 2D materials in a riskfree way, providing an alternative to exhaustive trialand-error experimentations. A multiscale framework developed by Xuan et al [49] bridges spatial scales that range from 10 −9 to 10 −3 m in space by combining DFT, ReaxFF and CFD (figures 26 and 27) to model the MOCVD gas phase physiochemical processes of the crystal growth of TMD. A ReaxFF based chemical kinetic model developed in this work ( figure 27(b)) provides a new and fundamental insight into kinetic and thermodynamic properties of elementary reaction pathways leading to the TMD growth from metal hexacarbonyl and chalcogen in a cold-horizontal MOCVD chamber. The full-scale CFD simulations of the MOCVD chamber showed that W(CO) 6 is rapidly consumed at early stage of the reaction before reaching the heated substrate by CO release reactions and H 2 Se addition reactions. CO-rich global reaction zones over the heated provide the basis for the formation of stable products such as CH 4 , H 2 O, CO 2 , CSe 2 , etc. Specific WCO(H x Se y ) related intermediate species like WCO(SeH) 2 Se 2 can reach and be located over the heated plate, which can result in the formation of a WSe 2 thin film with carbon-impurities. In the CFD simulations, W-species without CO (W(SeH) 2 Se 2 ) was observed in significant amount only towards the end of the heated substrate, indicating that depending on the local chemical environment all CO can be released from W-species before reaching the substrate ( figure 27(c)). Moreover, the increase in temperature and lumped species concentration resulted in more growth near the front edge of the growth substrate, in line with experimental observation ( figure 27(d)). All these results together provide a deep insight into the reactions prior to the growth.
In another work [116], the multiscale model development efforts further extended to phase field simulations. For this, ReaxFF simulations based on the force field developed by Nayir et al [50] were conducted to predict a thermodynamic model that describes the stability of growth fronts with various edge types and coverages as a function of local chemical environment (figure 26). Then this information was transferred to phase field models to predict the morphological evolution of 2D-domains as a function of growth conditions (e.g. T, P). This model revealed the impact of precursor concentration and homogenity of precursor distribution across the substrate on the coverage and continuity of 2D thin film. This model further showed that the substrate rotation is the key to the MOCVD/CVD growth of 2D materials, leading to a significant improvement in the uniformity and coverage of a 2D thin film.

Hybrid ML/ReaxFF applications
2D materials exhibit many unique properties and play a pivotal role in a wide range of research areas such as sensors [319,320], desalination membranes [34,317,321,322], energy storage [323,324] and biomedical applications [325][326][327]. However, the  mechanical properties of 2D materials cannot be accurately inferred from the chemical composition. Moreover, traditional methods only focus on one or two parameters for prediction, even with accelerated or multiscale computational efforts but it is still challenging to incorporate high-dimensional factors to models to fulfil the goal. ML algorithms, on the other hand, can readily process high-dimensional inputs and solve complex combinatorial optimization problems for predicting mechanical properties of 2D materials. For example, Yu et al [328] adopted the ConvLSTM models to predict the crack propagation of graphene with different orientations or defects. The datasets of single crystal and bicrystal graphene as well as defective graphene with three sets of single pore defects from ReaxFF simulations, which possessed the dimensional tensors of time sequences, rows and columns of images, and color channel, were input to the ConvLSTM architecture. The results indicated 99% and 98% accuracy for predicting the fracture path of graphene and different sets of defects. Xu et al [329] tested the four ML models, eXtreme Gradient Boosting (XGboost), multilayer perceptron, gradient boosting decision tree, and random forest, aiming to discover the wrinkling morphology and mechanical properties of nanocrystalline GOs (NCGOs) and GOs, with respect to the effects of grain size, oxidation, hydroxylation, epoxidation, GB hydroxylation, GB epoxidation, and GB oxidation by fitting the ReaxFF data [20]. The XGboost model revealed that the strength of NCGOs is greatly dictated by oxidation and grain size, and the hydroxyl group plays more prominent role in the strength of NCGOs than the epoxy group does. Zheng et al [330] maximized the toughness of GOs by strategically assigning functional groups by using deep reinforcement learning (RL) models. The RL approach reached optimized function group distribution within only 5000 rollouts, while the simplest designing task has 2 × 10 11 possibilities based on the datasets obtained from ReaxFF MD simulations using Chenoweth et al's force field [20].
Data driven approaches such as ML based multiscale models offer further a time-saving and cost-effective way to model the complex CVD reactions that lead to the growth of 2D material and enable to analyze and optimize the growth processes or parameters by the generation of a number of data using atomic scale simulations. For instance, two recent works [173,174] conducted ML assisted ReaxFF simulations using a similar model to Hong et al's [168] but with multimillion atoms, to investigate the nucleation, growth and coalescence of grains and defect healing during the CVD synthesis of MoS 2 growth on the MoO 3 surface deposited on Al 2 O 3 under the S 2 exposure. Taking the advantage of the data-driven approach with multimillion atoms, the authors were able to establish a correlation between the number of CVD heating/cooling cycles and defect healing. Some of the crucial takeaways from this work are summarized as follows: (i) first cooling/heating cycles drive the random nucleation and growth of grains with 2H, 1T, and disordered phases on a presulfurized MoO 3 surface. Upon the second heating/cooling treatment, 2H-phased grains coalescence into bigger grains and external heat further induces the phase transformation between 1T and 2H phases. As they grow, they coalesce and develop GBs. As the number of heating/cooling cycles increases, the size of grains becomes larger with reduced line defects and GBs. This result indicates the number of cooling/heating CVD cycles significantly influences the quality of the resultant structure, which offers a new way to control and repair structural impurities such as defects, GBs in a monolayer, thus, enhancing the crystal quality (figures 28(a)-(e)).
Patra et al [55] also employed ML to scale up their systems' size inaccessible at the MD level in order to investigate the kinetics and dynamics of spatial defect distribution and defect induced phase transition in an atomically thin material (e.g. MoS 2 ). In this approach, they combined genetic algorithm, ReaxFF simulations and TEM images. This ML-based model predicted that randomly distributed vacancy defects migrate into a thermodynamically more stable line defect, but requiring longer time duration over several minutes. These extended line defect formations were usually followed by the endothermic 2Hto-1T phase transitioning through the α-phase. The energy required for such a transition was typically supplied by heating and electron beam irradiation. The correlation of the size and shape of 1T domain to the intensity and alignment of line defects suggests that 2H-to-1T phase transition can be tuned in a controllable manner (figure 28(f)). Given the influence of defects on the atomic structure of 2D materials, similar to Patra et al [55], Banik et al [331] developed an effective model to accelerate the material search and discovery of impurities in a low-dimensional material, a decision tree reinforcement learning model, e.g. MC tree search (MCTS), combined with ReaxFF reactive simulations. This method cycles back and forth between MCTS and ReaxFF simulations, where atomic configurations of defective models are minimized at the ReaxFF level and the relaxation energies are passed to the MCTS as rewards. Then, MCTS update the defect configurations accordingly and feed them to ReaxFF simulations. This search continues until reaching the energetically optimal defect configurations. Thus, this model enabled to detect the energetically favorable pathways for defect migration and aggregation that include unstable intermediates and requires overcoming energy barrier, in a 2D material, e.g. MoS 2 after fewer search iterations than that of genetic algorithm based ML models ( figure 28(g)). Another work [332] addressed the acceleration of chemical reactions that are observed during the growth and characterization of 2D materials and require larger sizes and/or longer time scales than that are accessible to MD simulations, by developing ML assisted hybrid ReaxFF simulation method. Combining nonreactive (i.e. optimized potential for liquid simulations (OPLSs)) and ReaxFF reactive force fields with the use of dense neural network (DNN) algorithm enables to substantially decrease the computational time by accelerating the chemical reactions. As shown in figure 28(h) this ML-driven framework consists of three main steps: (i) the relaxation of topology based on OPLS, (ii) the identification/update of chemical reactions and bond topology of a system of interest at the ReaxFF level and (iii) cycling back to OPLS simulations by feeding the bond topology data from ReaxFF simulation using the DNN algorithm to predict the OPLS parameters and atomic charges. This approach cycles back and forth between OPLS and ReaxFF simulations until optimizing the parameters of interest. The case study of cross-linking reaction of dicumyl peroxide and decane reported in this work benchmarks the efficiency of this ML-driven approach in the acceleration of chemical reactions, by consequence, the substantial reduction of the computational time required for simulations (figures 28(j) and (k)).

Perspective and future directions
Despite significant advances, 2D materials' area still poses many interesting challenges and unexplored applications that need to be addressed through interdisciplinary approaches including both theoretical and experimental efforts. Here, we provide perspectives on knowledge gaps that we observe in this field and future research directions that we may pursue in this fascinating field.

Modeling and simulation of vapor-phase and colloidal synthesis of 2D materials
Synthesis of 2D materials (e.g. TMD, h-BN, graphene, MXenes) on substrates is quite challenging and timeand source-expensive, and there are still several open research problems that need to be addressed [21][22][23]. Some of the challenges that the Materials Science community faces in selective growth of large-area and high-quality 2D materials and their potential solutions using ReaxFF can be listed as follows: 9.1.1. Unidirectional alignment of 2D domains on substrate van der Waals epitaxy has emerged as a promising route for the synthesis of 2D materials on highly thermal-and lattice-mismatched substrates. However, controlling the unidirectional alignment of growth fronts on a substrate and their seamless stitching into a high-quality thin film is still challenging. Recent studies show that the nucleation and alignment of 2D domains are intimately tied to step edges, growth conditions that modulate the surface chemistry of substrate, and the underlying substrate symmetry [333][334][335][336]. Although some of these studies have used DFT calculations to model phenomena observed in experiments, a fundamental understanding of such chemically reactive growth process still requires multiscale approaches that should allow the dynamic evolution of 2D material growth process as a function of time under external stimuli (temperature, pressure, flow rate, etc).

Role of reaction intermediates
Synthesis of 2D materials from volatile and colloidal source is highly chemically reactive. This may lead to the formation of various reaction intermediates which indeed affect the resulting product's morphology and properties. Owing to the dynamic-bond connectivity adopted in the functional form, ReaxFF is extensively used to model chain of chemical reactions by large-scale simulations, thus enabling the tracking of the formation and behavior of these intermediates and providing insights into their role in the synthesis process.

Modeling solvent effect in colloidal synthesis
In colloidal synthesis, the solvent is another critical player in determining the course of the synthesis and morphology and properties of a resulting 2D material. Bond-order dependent ReaxFF method can be a good choice to accurately model (i) the solventmediated interactions between solvent and precursor molecules to understand the nucleation and growth of nanoparticles through acquiring the solubility and reactivity of precursor molecules from ReaxFF simulations, (ii) the interactions between solvent molecules and colloidal particles to gain an insight into the role of solvent in governing the morphology and structure of nanoparticles on the atomic scale, and (iii) the solvent's effects on the resulting properties of colloidal particles, such as their size, shape, and surface chemistry, which could be controlled by adjusting the concentrations of precursor molecules and reaction conditions using ReaxFF simulations.

Surface functionalization and reactions
As also emphasized in another review [337], surface functionalization of 2D materials with various polymers or nanoparticles makes them promising candidates for a broad range of applications, such as membrane-based separation, sensing, and biomedical applications, due to their outstanding features. However, the state-of-the-art fine-tuning of the interlayer spacing and stability in physiological conditions of 2D materials has remained a challenge. Furthermore, the mechanisms of various surface reactions between precursors and substrate, which are promoted by the active sites of pristine and functionalized 2D material surfaces, play a critical role in the synthesis of 2D materials; however they are not fully understood yet. ReaxFF can be used to accurately model these interactions and their impact on synthesis, thus, guiding experiments to grow highquality and large-area atomically thin materials.

Effect of external parameters
The morphology and materials properties of 2D materials are significantly affected by operating conditions such as temperature, pressure, and precursor concentration. To gain an atomic insight into physical and chemical behavior of 2D materials under external stimuli during the growth is essential for the development and improvement of 2D materialsbased devices. Large-scale atomic simulations based on ReaxFF come forth as an effective tool to study stimulus-driven dynamic evolution of 2D domains. 9.2. 2D materials and organic hybrids 2D materials, for example, MXenes combined with organic materials have recently attracted attention as they can provide unique thermal and mechanical properties. However, the fundamental understanding of interactions between MXenes and organic molecules has not been fully grasped. Theoretical approaches can provide insights into the covalent and non-covalent interactions between MXene and organic molecules to guide synthesis. Exploring optimal conditions that can control chemically active and functionalized MXene surfaces is also essential for designing fine-tuned hybrid materials.

Half van der Waals growth of group III metals and their nitrides
The real-world semiconductor industry calls for the implementation of group III nitrides, including AlN, GaN, InN, and their ternary and quaternary compounds into a myriad of next-generation high-frequency power electronic devices and optoelectronics, such as light emitting diodes, lasers, photodetectors, and solar cells [338,339]. The key to understanding and engineering the epitaxial growth process of group III nitrides in an MOCVD system is to acquire the knowledge of gas phase reaction mechanisms through predictive modeling [340]. Future ReaxFF force field development and the research can step in the following directions: (a) modeling the gas phase reaction chemistry during the MOCVD growth of group III nitrides, with an aim to understand the reactivity in group III precursor species with ammonia; (b) modeling crystal lattice mismatch and defect dynamics, along with the solutions to compensate mismatch or stress concentration on interfaces between group III nitrides and foreign substrates; and (c) modeling surface kinetics with varying mixing ratios of groups III and V materials, so as to mimic the compound layers grown in mass transport limited conditions.

Phonon dynamics with ReaxFF and ML
van der Waals layered materials exhibit superlative electronic and optical properties and promise applications in novel nanoelectronic devices [341,342]. The pre-existing or introduced structural defects, such as vacancies, dislocations, and GBs, govern the electron and thermal transport properties of these devices [343]. A fundamental understanding of the phonon dynamics is, thus, of vital importance for implementations of these materials in nanotechnology. Simulations for these defective systems require largescale (in nm scale) atomic models, first-principles methods are thus computationally expensive to be utilized. MD simulation is feasible to handle a largescale system, and a validated empirical potential is indispensable for characterizing the phonon properties of defective materials. Several ReaxFF reactive force fields have recently been developed for van 2D Mater. 10 (2023) [24,25,48,50]. Through training against first-principles data, one can tailor these force fields to well describe the phonon dynamics of TMD systems. Combining reactive MD simulations with ML techniques for exhausting the distributions of defective systems [55] and experimental measurements, such as Raman spectra and high-resolution transmission electron microscopy (HRTEM), would be a promising method to advance the understanding of phonon dynamics and structural evolutions of defective TMD systems and pave the way for their applications in nanodevices.

2D confinement by assembly of van der Waals heterostructures for molecular and ion transport
Despite significant progress made in the assembly of van der Waals 2D materials for ionic and molecular transport, there are still several knowledge gaps that need to be addressed, such as scalability, reproducibility, dynamics of transport, selectivity and control, as well as the stability and durability [344]. Particularly, large-scale modeling and simulation of such complex systems are still infant, and the development of new reactive potentials is sorely on demand. One of the major challenges is an accurate representation of the complex interactions of ion/molecules with interfaces due to the confinement and surface interactions. Therefore, there is a need for accurate and efficient models that can capture the complex dynamics of molecules/2D materials interfaces through sophisticated force fields such as ReaxFF. To date, ReaxFF force fields have been developed to address some of these challenges, as discussed in this review, to guide experimental studies to design novel materials and devices, but further refinement and validation are still necessary.

Optimizing ReaxFF parameters with guidance from ML
Handling the high dimensional PS of ReaxFF during the parameterization has been troublesome and requires well-developed optimization algorithms. The reliability and accuracy of a force field strongly relies on the communication between QM-based training data set and test data set which can be produced by MD simulations. Several studies show that ML algorithms are proven themselves as an effective means to automate such communication during the optimization. In future studies, we aim to develop a dynamic data-driven approach in a combination of ML, ab initio MD (AIMD) and classic MD simulations for the force field optimization. With the use of AIMD simulations, we will use forces that will set as reference values for force field parameters to be trained. This will immensely reduce the computational time and sources required for the force calculations, by consequence, the force field parameterization process, particularly for the CVD/MOCVD growth of 2D materials including complex gas-gas and gas-solid reaction networks. Furthermore, with the help of ML algorithm, this tandem approach will cycle back and forth between the training and test data sets until reaching an optimized parameter set of ReaxFF force field.

e-ReaxFF
Since standard ReaxFF description is lack of the explicit electron description it is incapable of accurately simulating the electron dynamics of reactions such as redox reactions and polarization of bulk material. In response to this need, e-ReaxFF [345] is developed for explicit treatment of electrons within the ReaxFF force field framework, which is computationally less demanding and several magnitudes of order faster than QM method. Comparative studies of e-ReaxFF and QM methods on hydrocarbon molecules show the capability of eReaxFF in simulating the electron transfer dynamics in inexpensive way. Another study [105] successfully applied eReaxFF to graphitic anodes of Li-batteries shows that (i) excess electrons are localized nearby defects in a material, in agreement with those generated by QM method and (ii) lithium metal formation around the graphitic anode of lithium batteries initiates with electron transfer from graphitic anode to Li + ions. Due to the high surface-volume ratios, moderate band gaps, and short carrier transport paths, 2D materials assembled into van der Waals heterostructures hold great promise for energy storage applications. For instance, a direct Z-scheme and Janus heterostructures with anisotropy shape and chemical composition present outstanding performance in catalyzing the hydrogen evolution reactions. In the future, we will extend e-ReaxFF applications to 2D materials for large-scale and cost-effective modeling of such as electron distribution, ion and/or cation exchange reactions in a 2D material.

Concluding remarks
In summary, ReaxFF simulations offer a computationally cost-effective way to connect atomic-scale impurities to larger morphologies, growth modes, and material properties. As discussed in the previous sections, the use of ReaxFF can provide a valuable tool for tackling open research problems encountered in the modeling and simulation of synthesis and characterization of 2D materials. It helps to guide experimental studies and the development of new materials and devices. However, it is important to note that 2D materials' field is an active area of research and is rapidly evolving. Therefore, there is a growing demand for the development and improvement of more sophisticated and accurate ReaxFF reactive force fields. These advancements will not only help the Materials Science community establish a theoretical understanding of the behavior of these materials but also offer the advantage of relatively low computational cost compared to QM-based methods.

Data availability statement
No new data were created or analysed in this study.

Acknowledgments
Financial support was provided by the National Science Foundation (NSF) through DMR-1808900 and The Two Dimensional Crystal Consortium−Materials Innovation Platform (2DCC-MIP) at The Pennsylvania State University under NSF cooperative Agreements DMR-1539916 and DMR-2039351.