{'fraction_non_alphanumeric': 0.026814111660224103, 'fraction_numerical': 0.01365171013358125, 'mean_word_length': 4.854482956173016, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The dynamics of copper intercalated molybdenum ditelluride

Layered transition metal dichalcogenides are emerging as key materials in nanoelectronics and energy applications. Predictive models to understand their growth, thermomechanical properties and interactions with metals are needed in order to accelerate their incorporation into commercial products. Interatomic potentials enable large-scale atomistic simulations at the device level, beyond the range of applications of first principle methods. We present a ReaxFF reactive force field to describe molybdenum ditelluride and its interactions with copper. We optimized the force field parameters to describe the properties of layered MoTe2 in various phases, the intercalation of Cu atoms and clusters within its van der Waals gap, including a proper description of energetics, charges and mechanical properties. The training set consists of an extensive set of first principle calculations computed from density functional theory. We use the force field to study the adhesion of a single layer MoTe2 on a Cu(111) surface and the results are in good agreement with density functional theory, even though such structures were not part of the training set. We characterized the mobility of the Cu ions intercalated into MoTe2 under the presence of an external electric fields via molecular dynamics simulations. The results show a significant increase in drift velocity for electric fields of approximately 0.4 V/A and that mobility increases with Cu ion concentration.


Introduction
The integration of single-layer and few-layers transition metal dichalcogenides (TMDs) into nanoelectronics [1] is paving the way for future developments of low power, ultrathin electronic devices [2,3]. carrier injection at the interface with the metal contact. Various metals have been investigated in order to improve contacts [19] as well as local TMD processing such as intercalation, phase engineering [20] and chemical doping [10]. Due of the wide compositional flexibility of TMDs and the various metals of interest, brute force experimental optimizations are likely to be ineffective and predictive computational techniques can help narrow down experimental studies to the most promising candidates. Electronic structure calculations using density functional theory (DFT), including corrections for the nonlocal dispersive forces part of vdW interactions [21,22], have played and continue to play a key role in this field [11,23,24,25].
However, the computational intensity of DFT calculations has limited its application to relatively small size, model systems composed of few hundreds of atoms. Thus, most of the DFT studies have focused on defect free TMD/metal interfaces (limited by the small size of the supercell) and little is known about the nanostructure of the interface.
Direct synthesis of single to few layers TMDs have been demonstrated however, even with sophisticated chemical vapor deposition (CVD) techniques [26,27], TMDs present high defect density [28].
Understanding defect formation in TMDs would not only enable the synthesis of large area, defect-free layers but would also provide the knowledge for defect-controlled TMD growth. Structural defects in TMDs, including point defects and grain boundary influence its electronic, thermal and mechanical properties [29,30,31]. Unfortunately, the atomistic mechanisms of nucleation and growth during the CVD process remain unclear.
Large-scale atomistic simulations would be critical to study the interactions between defects and interfaces and to provide atomistic insight into the process of synthesis and growth of these materials.
This paper contributes to bridge the gap between first principle calculations and experiments via the development of the first reactive force field capable of describing a TMD, a metal and their interactions, specifically molybdenum ditelluride and Cu. We use the ReaxFF formalism that has shown the flexibility to describe a wide range of materials including metals, semiconductors, their oxides [32,33] and molecular materials [34,35,36]. The force field is parameterized to describe various phases of Cu, Mo and MoTe 2 as well as the corresponding interfaces, the intercalation of Cu atoms and small clusters in the vdW gap of multiple layers of MoTe 2 including a proper description of the relative energies, charges and mechanical properties.
The rest of the paper is organized as follows. In Section 2 we present details about the simulations and methods. Section 3 describes the optimized force field and presents its comparison with DFT calculations.
In Sections 4 and 5 we study the interface between monolayer MoTe 2 and a Copper slab and the dynamics of Cu-intercalated MoTe 2 under various external voltages, respectively. We finally draw conclusions in Section 6.

Force field optimization approach and methods
ReaxFF is a bond order interatomic potential and naturally allows bond breaking and forming during molecular dynamics (MD) simulations [37]. The total energy is defined as the sum of various contributions including covalent (bond stretch, angles, torsions, and over/under coordination) and non-bonded interactions (van der Waals and Coulomb); additional terms to account for lone pair and bond conjugation. Electrostatic interactions are described using environment-dependent partial atomic charges which are computed at each step of the simulation using the charge equilibration method [38,39] (QEq). All contributions to the covalent energy depend on the partial bond order between atoms, a function that smoothly decays to zero when the bond between atoms is broken. Details on the mathematical form of ReaxFF can be found in the original paper, Ref. [37] and a recent review of applications in Ref. [40]. In order to determine the parameters that describe the functional form of the interaction, the force field has to be optimized against a training set consisting of first principle calculations or a combination between first principles calculations and experiments. The parameterization of such a force field is complex due to the large number of parameters and the need for extensive training sets and time-consuming optimization cycles. Thus, despite the enormous progress to date [40], several materials of technological importance lack accurate force field descriptions.

DFT simulation details
All DFT calculations have been performed with the Vienna Ab initio Simulation Package (VASP) [41,42] within the generalized gradient approximation as proposed by Perdew, Burke, and Ernzerhof (PBE) [43].
In addition, we use the non-local vdW density functional of Langreth and Lundqvist [44,45] (vdW-DF) throughout the paper to correct for the London dispersion, poorly described by PBE functionals. Atomic structures have been relaxed until energies and forces reach a tolerance of 1×10 −4 eV and 1×10 −2 eV/Å, respectively. The kinetic energy cutoff for the plane-wave basis was set to 500 eV and we integrate the k-mesh with various number of points depending on the size of the supercell. We use a dense 8×8×8

Parameter optimization
In order to optimize the force field we use an in-house code written in Python that implements a Monte Carlo-based (MC) simulated annealing algorithm [47] and coupled to the LAMMPS simulator [48] (via its Python interface) to evaluate energies, forces, pressures and charges of molecules and crystals included in the training set. This method has been used previously to optimize force fields for oxides [49] and metals [50]. The training set is included as an XML database of DFT calculations describing the dissociation of small molecules and equation of states for crystals. The goal of the simulated annealing optimization is to minimize a total error function ( tot ) defined as the sum of individual discrepancies between the force field and DFT data, corresponding to energy (E), force (f), pressure (P) and charge (Q) as: with N X a set of normalization constants chosen such that each individual error appears with equivalent weight; these parameters are also used to establish the relative importance between different quan-tities in the training set.
Before starting the simulated annealing run we perform an initial sensitivity analysis for each parameter in the force field in order to determine the steps by which each parameter will be stochastically modified during the MC run. A Metropolis criterion allows the random selection of some "non-optimum" parameters in order to expand the exploration space and avoid the force field function becoming trapped in local minima. The temperature used for the Metropolis acceptance is decreased during the simulation in order to converge towards the global minimum. We note that the QEq parameters are optimized first in order to reproduce partial atomic charges of the structures in the training set obtained from Bader analysis [51] of the DFT runs. This is done by setting E , f and P to zero in the energy expression so that only partial charges are compared. The simulated annealing code distributes MD simulations in parralel (using IPython for parallel computing) and at its current stage of development can perform 1,000 single point calculations in approximatively 5 seconds over 64 cores.
3 Cu-MoTe 2 reactive force field MoTe 2 has recently received significant attention because of the discovery of a metallic T' phase [52,53] energetically close [17,18] to its ground state semiconducting H phase, common to group VI TMDs.
Therefore, we decided to develop parameters for MoTe 2 rather than the more studied molybdenum disulfide (MoS 2 ), already explored by MD [54]. Copper is the metal of choice for electrodes and metallic interconnects; hence, we studied the system Cu-MoTe 2 . The following subsections describe the training data used in the force field parameterization and compare the optimized ReaxFF with the DFT training data.

Initial parameters and charge equilibration
We built an initial force field with parameters for Mo and Te extracted from Ref. [55] and Cu from Ref. [56].
The general parameters were re-optimized in order to improve the atomic features including equation of state and formation energy, as described in details in the following sections. Therefore, the parameters presented here are not transferable to the initial ReaxFF they belong to and the new force field we propose corresponds to a new branch of the ReaxFF tree, as defined in Ref. [40]. Additionally, we introduce three anglar terms corresponding to Te-Mo-Te, Te-Cu-Te and Cu-Te-Cu in order to describe the subtle energy difference between MoTe 2 phases and the potential energy landscape of Cu diffusion in MoTe 2 .
We first optimized the QEq parameters, traditionally named γ EEM , χ EEM and η EEM representing the shielding distance for Coulomb interactions, electronegativity and hardness parameters, respectively.
Partial atomic charges were optimized on each atom against Bader charges. We show in Table 1 the root mean square error of charges per atom computed between the optimized ReaxFF and DFT. ReaxFF describes accurately charges with an average error of 0.037 ± 0.015 e per atom compared to Bader charges.

Bulk phases and bond dissociation
A key feature of ReaxFF lies in its ability to describe atoms in various chemical environments, corresponding to different coordinations. Therefore, we fit ReaxFF against various bulk phases of Cu and Mo including simple cubic (SC), face centered cubic (FCC) and body centered cubic (BCC). Figure 1 shows various equations of state (EOS) computed with DFT and ReaxFF, with minimum energy shifted to their corresponding enthalpies of formation (see Table 1 for details). For both Mo and Cu, ReaxFF describes with accuracy the ground state structures, only the simple cubic phases present some discrepancy with respect to DFT. However, the simple cubic phases lie tens of kcal/mol higher than the ground state and therefore do not represent important contribution to their chemistry; it is only important that they remain energetically prohibited.
MoTe 2 has been observed in the stable H and T' phases as well as a metastable T phase, corresponding to trigonal prismatic, distorted octahedral and octahedral geometries, respectively. The semiconducting H phase has been shown to be the ground state configuration for Mo-based TMDs while the metallic T' phase lies few meV above the H for MoTe 2 . A transition from H to T' has been predicted from DFT calculations by mechanical strain in the range of 0.3 to 3% [17] or an applied electrostatic gate voltage [18]. Figure 2 shows the EOS corresponding to various phases of bulk MoTe 2 computed with DFT and ReaxFF, and shifted to their binding energy. The binding energy is defined as the total energy of the MoTe 2 phase minus the total energies of bulk Mo in its ground state BCC phase and the diatomic molecule Te 2 . The force field predicts the H phase to be the ground state however, the T' phase appears higher than the T phase. Additionally, we computed the potential energy surface corresponding to transition between H→T' within the rectangular unit cell constrained to the initial H phase, as depicted in Figure 3. We found that ReaxFF predicts the energy barrier between H and T' phases E A =32.3 kcal/mol in good agreement compared to 35.7 kcal/mol computed from DFT calculations. We note that ReaxFF is able to capture subtle energy differences between the various phases of MoTe 2 .
The bond order nature of ReaxFF enables dynamical connectivity during molecular dynamics simulations and the ability to describe dissociation (and formation) of chemical bonds. Therefore, the training set includes various bond dissociation curves, such as Te-Te and Mo-Te bonds in Te 2 molecule and mono-  layer MoTe 2 , respectively. The corresponding energies versus bond distance plots are shown in Figure   4. We found a good agreement between ReaxFF and DFT even though the force field predicts slightly lower dissociation energies for both Te-Te and Mo-Te dissociations.
Finally, in order to validate the stability of the force field we performed molecular dynamics simulations of bulk and monolayer MoTe 2 in the microcanonical ensemble with a timestep of 0.5 fs. We found little drift in the total energy (∼ 1×10 −6 kcal/mol/ps/atom) suggesting a good stability of the force field.

Cu-intercalated MoTe 2
The intercalation of Cu atoms within the vdW gaps of ground state H-MoTe 2 can occur either at octahedral (or hollow, noted "h" ) or tetrahedral ("t") sites. Figure 5 shows the energetics of h- (Fig. 5a) and t- (Fig. 5b) intercalated Cu x MoTe 2 for Cu concentrations x=0.125, 0.25, 0.375 and 0.5. We show energy as a function of volume upon uniaxial deformation along the direction perpendicular to the plane of the TMD in order to better test the interactions between Cu and the TMDs. The energy curves in Figure 5 are referenced to the perfect TMD and Cu as:   To summarize, we report on Table 1 the root mean square error per atom corresponding to charges, forces and pressure computed between DFT and ReaxFF as well as the various formation energies.
We acknowledge that the force field presents some uncertainties, and the main deviations, with respect to DFT calculations, are related to the energy barriers for Cu diffusion and cluster dissociation. The overestimation of diffusion barriers suggest a slower mobility of Cu atoms described by ReaxFF than predicted from DFT. We provide the optimized ReaxFF parameters in Supplementary Materials.

Mechanical properties
In addition to the energetics analysis, we evaluated the capacity of ReaxFF to describe the mechanical properties of the structures included in the training set. Equilibrium volumes, bulk moduli and elastic constants are reported on Table 1, computed with both DFT and ReaxFF. We found between 15 and

Force field validation: Cu/MoTe 2 interface
We used the force field to compute the interaction between a single layer of MoTe 2 and a (111) Cu surface. We follow the strategy developed in Ref. [57] to define supercells and minimize strain, chosen to be applied to the metal slab because of the corresponding small effect on the electronic structure of the interface [25]. We develop the in-plane supercell lattice vectors T 1 = n 1 a 1 + n 2 a 2 and T 1 =    Table 2: Parameters corresponding to the supercell construction of a Cu/MoTe 2 interface with minimal lattice mismatch (δ). The two sets (n 1 , n 2 ) and (m 1 , m 2 ) correspond to the coefficients of the in-plane supercells for MoTe 2 and Cu, respectively. We note that the two supercells have been rotated by the angle α.  parameterization. This indicates the type of regimes the force field is expected to be extrapolated to.
However, as with any force field, care should be exercised when using it away for the conditions for which it was trained.

Diffusion of Cu-intercalated MoTe 2
The relatively small computational intensity of ReaxFF enables the simulations of large scale systems; providing a unique tool for exploring ionic transport, fracture and the effect of defects on mechanical properties and growth. In this section we focus on the first of these examples.
We know from the first principle calculations presented above that the activation energy for thermal hop of Cu intercalated MoTe 2 varies between 9-13 kcal/mol (0.39-0.56 eV) range, depending on the intercalation site. These values are intermediate between the extreme cases of a good ion conductor like highly metal doped chalcogenide glasses with E A ∼ 0.2 eV [58] and amorphous silica with in average E A ∼ 0.8 eV [59]. In order to confirm this last statement, we performed simulations of Cu-intercalated bulk MoTe 2 at various temperature ranging from 340 to 380K and evaluated the diffusion coefficient D ≈ We now explore the effect of Cu concentration on ion mobility in Cu-intercalated bulk MoTe 2 under the action of an external electric field. We performed MD simulations of Cu-intercalated bulk H-MoTe 2 in the canonical ensemble at 300 K for up to 2 ns and applied various external electric fields (i.e. a constant force F i = Eq i to each atom i) ranging from 0.0 to 0.6 V/Å. The direction of the electric field is chosen to be aligned with the in-plane lattice parameter of the hexagonal cell of MoTe 2 , as shown in the inset of Figure 9. We varied the concentration of Cu ions from 1 to 4% with respect to MoTe 2 and study their mobility. In order to estimate uncertainties in the mobility of the Cu ions, we performed simulations over an ensemble of 4 independent samples differing by the initial position and velocity of copper atoms.
Since thermal diffusion leads to no net motion of the atoms, the average position of the Cu ions as a function of time < X(t) >, can be described as: The mobility u can then be obtained from the applied electric field E and drift velocity u = v d /E [60]. Figure 9 shows the computed drift velocity of Cu-intercalated MoTe 2 as a function of the electric field (left) and the Cu fraction (right). We found that a minimum electric field of approximately 0.4 V/Å is required to displace the Cu ions in the vdW gap of MoTe 2 within MD timescales. More importantly, we find that the ionic mobility increases with concentration. Finally, since ReaxFF predicts slightly higher diffusion barriers, with respect to DFT calculations, we expect the electric field threshold for drift diffusion to be lower than the 0.4 V/Å predicted here.

Conclusions
We presented a parameterization of ReaxFF to describe the interactions between MoTe 2 and copper and showed that such force field can be used to study large scale atomistic simulations of interest for emerging electronics. The optimization of the force field against first principle calculations shows that ReaxFF reproduces with accuracy the energetics, charges and mechanical properties of various systems composed of MoTe 2 and Cu and slightly overestimates the energy barrier for Cu diffusion in the vdW gap of the TMD. Additionally, we demonstrated the good transferability of the force field to describe Cu/MoTe 2 interfaces. As for any force field, care must be taken when used outside of the conditions used during parameterization. For example, the force field has not been optimized to describe defect formation energies and extensive validations must be performed in order to study growth, as we envisage in the future.