Transferability of Molecular Potentials for 2D Molybdenum Disulphide

An ability of different molecular potentials to reproduce the properties of 2D molybdenum disulphide polymorphs is examined. Structural and mechanical properties, as well as phonon dispersion of the 1H, 1T and 1T’ single-layer MoS2 (SL MoS2) phases, were obtained using density functional theory (DFT) and molecular statics calculations (MS) with Stillinger-Weber, REBO, SNAP and ReaxFF interatomic potentials. Quantitative systematic comparison and discussion of the results obtained are reported.


Introduction
Group 6 transition metal dichalcogenide (G6-TMD) two-dimensional (2D) nanomaterials [1], and especially single-layer molybdenum disulphide (SL MoS 2 ), are probably the second most studied 2D materials following graphene [2]. The major disadvantages of graphene are the lack of a band gap in the electronic spectrum, its susceptibility to oxidative environments and that it has some toxic properties. That is why scientists and engineers, beyond ordinary human curiosity, have begun to look for materials free of these deficiencies [3,4].
Both synthetic and natural bulk transition metal dichalcogenides have layered structures with two primary distinguished allotropic forms, 2H and 3R, belonging to the hexagonal crystal family, but differing in a sequence of arrangement. Strong triple layers of metal-sulphur-metal are weakly bounded by the van der Waals forces, similar to graphene in graphite [1].
The most accurate methods of solid state physics are based on quantum mechanics, unfortunately, with the accuracy of the methods their cost increases. The number of atoms and the number of timesteps that can be analysed with the first-principles method using either energy minimisation or ab initio molecular dynamics (AIMD) is highly limited. For typical computational resources currently available, the use of these methods is limited to several hundreds of atoms for less than about several picoseconds. These restrictions justify the need for more approximate methods, such as molecular methods [7].
In general, there is a lack of perfectly transferable interatomic potential that would work with the various materials and systems we are interested in. Some are more transferable, others less [8]. It depends on the physics behind them, the mathematical flexibility of the model capable of describing the multimodal potential energy surface (PES) and the quality of the fitting process and, of course, on the "difficulty" of the material [7].
According to the author's best knowledge, there are no publications where the performance of different molecular potentials for molybdenum disulphide is analysed for all phases of SL MoS 2 , there are only partial comparisons, and so in [9] the results for 1H and 1T phases for potentials Stillinger-Weber, REBO and ReaxFF are only compared between each other. In [10], the geometric parameters and mechanical properties of 1H phase obtained from Stillinger-Weber and REBO potentials are compared with density functional theory (DFT) calculations. Thermal transport properties in 1H phase from molecular dynamics using Stillinger-Weber and REBO potentials were obtained in [11].
A partial comparison of different potentials for the 1H phase SL MoS 2 can be found in papers where new parametrisations are presented, e.g., [12][13][14]. There are also publications where using molecular simulations the authors try to determine certain SL MoS 2 properties that were not taken into account during the parametrization of potential, e.g., [15][16][17][18].
The paper is organised as follows. Following the above Section 1, Section 2.1 presents the computational methodology used in ab initio calculations of analysed structures and Section 2.2 describes the computational methodology used in molecular calculations and molecular potentials examined: four Stillinger-Weber (SW) potentials [19], the reactive many-body (REBO) potential, the spectral neighbour analysis potential (SNAP) and the reactive force-field (ReaxFF). Section 3 presents the structural and mechanical properties of SL MoS 2 and phonon spectra obtained from the ab initio and molecular calculations and evaluates the quality of the analysed potentials. The last Section 4 summarises and concludes the results obtained.

Computational Methodology
Analysing the available literature concerning phases of SL MoS 2 , it is not feasible to find all structural, mechanical and phonon data obtained in one consistent way. The availability of experimental data is actually limited to phase 1H only and therefore we must use ab initio calculations. Unfortunately, also ab initio calculations, most often DFT, differ in the calculation methodology, i.e., they use different functional bases, different pseudopotential or exchange-correlation (XC) functionals, and such a parameter as cohesive energy is not accessible at all. For this reason, structural and mechanical data-lattice parameters, average cohesive energy, average bond length, average height, 2D elastic constants as well as phonon data-are determined using a single consistent first-principle approach as described in the next Section 2.1. These data will be further considered as reference data and marked as Value DFT . Then the same data were determined, as described in Section 2.2, using the analysed molecular potentials Section 2.2 and will be marked as Value potential . Having both data, we can simply define mean absolute percentage error (MAPE): that will allow us to quantify the potentials under examination. Phonons were determined only for the three best, having the lowest ∑MAPE, molecular potentials.

Ab Initio Calculations
Ab initio computations by means of the density functional theory (DFT) [20,21] and the pseudopotential plane-wave approximation (PP-PW) programmed in ABINIT [22,23] code were done in the present study. Optimised norm-conserving Vanderbilt pseudopotentials (ONCVPSP) [24] were utilised to describe the interactions of non-valence electrons and ionic core. ONCVPSP pseudopotentials used were taken from PseudoDojo project [25].
To strengthen the reliability of the calculations as an exchange-correlation (XC) functional, three approximations were initially checked for their ability to reproduce the geometry of 1H-MoS 2 : local density approximation (LDA) [26,27], classical Perdew-Burke-Ernzerhof (PBE) generalised gradient approximation (GGA) [28] and modified Perdew-Burke-Ernzerhof GGA for solids (PBEsol) [29]. To provide access to all XC functionals used a library of exchange-correlation functionals for density functional theory, LibXC [30] was utilised.
All the computations were done by a proper adjustment of their precision, what was achieved by automatically set up the variables at accuracy level 4 (accuracy = 4 matches the default settings of ABINIT). The cut-off energy in line with ONCVPSP pseudopotentials of the plane-wave basis set was fixed at 35 Ha (952.4 eV) with 4d 5 5s 1 valence electrons for Mo and 3s 2 3p 4 valence electrons for S. K-PoinTs grids were derived with kptrlen = 35.0 (grids that specify a length of the smallest vector LARGER than kptrlen). In all the present ABINIT computations, the metallic occupation of levels with the Fermi-Dirac smearing occupation scheme and tsmear (Ha) = 0.02 was applied.
Initial data defining unit cells of SL 1H-MoS 2 , 1T-MoS 2 and 1T'-MoS 2 were taken from [31] and then all structures were relaxed by applying the BFGS minimisation scheme with full optimisation of cell geometry and atomic coordinates (a two-stage scheme was used here: in the first one, the ionic positions without cell shape and size optimization, and in the second, the full optimization of cell geometry). Tolerance for maximum stress (GPa) was specified as 1 × 10 −4 .
The cohesive energy the E c (MoS 2 ) was calculated, taking into account stoichiometry, as the total energy E total (MoS 2 ) difference of 2D molybdenum disulphide and a single Mo atom energy E iso (Mo) in a sufficiently large box and a single S atom energy E iso (S) in a similar large box [32]: The theoretical ground state elasticity tensor, C ij , of all the structures analysed, was identified with the metric tensor formulation of strain in density functional perturbation theory (DFPT) [33].
In order to examine the elastic (mechanical, Born) stability of all the structures, positive definiteness of the elasticity tensor was verified [34] by computing Kelvin moduli, i.e., eigenvalues of the elasticity tensor represented in second-rank tensor notation [35,36].
To get the elastic constants, C ij , for all pre-relaxed structures, the stress-strain method with the maximum strain amplitude of 10 −6 was utilised [41,43].
For the phonon calculations, phonoLAMMPS (LAMMPS interface for phonon calculations using Phonopy code [44]) [45] was utilised. Supercell and finite displacement approaches were used with 3×3×1 supercell of the unit cell and the atomic displacement distance of 0.01 Å. The cohesive energy, E c (MoS 2 ), (Equation (2)) in molecular calculations is simply potential energy.
• REBO [48]: the reactive many-body potential (REBO) fitted to structure and energetics of Mo molecules, three-dimensional Mo crystals, two-dimensional Mo structures, small S molecules and binary Mo-S crystal structures. • SNAP [49]: the machine-learning-based spectral neighbour analysis potential (SNAP) fitted to total energies and interatomic forces in SL 1H-MoS 2 obtained from firstprinciples density functional theory (DFT) calculations. • ReaxFF [50]: the reactive force-field (ReaxFF) parameters fitted to a training set of geometries, energies, and charges derived from DFT calculations for both clusters and periodic Mo x S y systems.

Results
The first step of the ab initio calculation was to select the exchange-correlation (XC) functional that most accurately reproduces the experimental geometry of 1H-MoS 2 . The measured lattice constant for SL 1H-MoS 2 a = 3.157 Å and average height (vertical separation between S atoms) h = 3.116 Å [50], while that calculated with the local density approximation (LDA) a = 3.144 Å, h = 3.111 Å, with the classical Perdew-Burke-Ernzerhof (PBE) generalised gradient approximation (GGA) a = 3.220 Å, h = 3.121 Å and with the modified Perdew-Burke-Ernzerhof GGA for solids (PBEsol) a = 3.165 Å, h = 3.120 Å, a similar trend can also be observed in other papers [12,51]. Once again, it was confirmed that the PBEsol is the overall best performing XC functional for identifying the structural and mechanical properties [4,52,53] and thus all subsequent calculations will use PBEsol XC functional.

Structural and Mechanical Properties
The basic cell for the SL 1H-MoS 2 polymorph is depicted in Figure 1 (hP3 in Pearson notation, P6m2-space group in Hermann-Mauguin notation, no.187-space group in the International Union of Crystallography (IUCr) notation), SL 1T-MoS 2 polymorph is shown in Figure 2 (hP3, P3m1, no.164) and SL 1T'-MoS 2 is depicted in Figure 3 (oP6, P2 1 /m, no.11), respectively [54]. Although 2D structures were studied, as is commonly practised, the 3D notation is used here. The crystallographic data for all calculated phases are additionally stored in crystallographic information files (CIFs) in Appendix A. Determined from DFT calculations structural and mechanical properties, namely, lattice parameters, average cohesive energy, average bond length, average height, 2D elastic constants and 2D Kelvin moduli, of the three analysed SL MoS 2 allotropes are gathered in Table 1. It can be seen that the calculated values match well the available experimental data [50] as well as those from other calculations [55][56][57]. This can be regarded as a confirmation of the correctness of the applied methodology. It is worth observing that the trend in the calculated cohesive energy matches the stability of the analysed phases and adding that all the calculated values are obtained using one consistent methodological approach. All calculated 2D Kelvin moduli for the three analysed phases are positive, which translates into mechanical stability. Calculated with the use of molecular statics and different molecular potentials, twelve structural and mechanical properties, namely, lattice parameters, average cohesive energy, average bond length, average height, 2D elastic constants and 2D Kelvin moduli, of the SL 1H-MoS 2 phase are collected in Table 2. The results obtained are then compared with those from DFT and quantified by calculating the mean absolute percentage error (MAPE) using the Equation (1). What follows from the results obtained? Overall, analysing the MAPE 1H for 1H-MoS 2 , the three most accurate potentials are: SW2017, SNAP, REBO, and the least are: ReaxFF and SW2015. A detailed look shows that only two potentials correctly reproduce cohesive energy. Mechanical stability is correctly reproduced by all potentials, i.e., K i > 0. Potential ReaxFF catastrophically badly reproduces mechanical properties of 1H-MoS 2 , even the symmetry of the elasticity tensor is not correct. The computed twelve structural and mechanical properties of the SL 1T-MoS 2 phase are summarised in Table 3. In general, analysing the MAPE 1T for 1T-MoS 2 , the three most accurate potentials are: SW2015, SW2016 and SW2017, and the least are: SNAP and ReaxFF. A detailed look shows that only two potentials correctly reproduce cohesive energy. Mechanical stability is correctly reproduced by all potentials, i.e., K i > 0. Potential ReaxFF again catastrophically badly reproduces mechanical properties of 1T-MoS 2 , even the symmetry of the stiffness tensor is again not correct. Unfortunately, the three potentials: SW2013, SW2015, SW2016, do not correctly reproduce the symmetry of the 1T-MoS 2 phase, i.e., during pre-relaxation input 1T-MoS 2 converges to 1H-MoS 2 phase.
The identified thirteen structural and mechanical properties of the SL 1T'-MoS 2 phase are summarised in Table 4. In general, analysing the MAPE 1T' for 1T'-MoS 2 , the three most accurate potentials are: SW2016, REBO and SW2017, and the least are: SNAP and ReaxFF. Once again, only two potentials correctly reproduce cohesive energy. Mechanical stability is reproduced in the right way by all potentials, i.e., K i > 0. Unfortunately, the two potentials, SW2017 and SNAP, do not properly restore the symmetry of the 1T'-MoS 2 phase, i.e., during pre-relaxation input 1T'-MoS 2 basic cell converges to 1T-MoS 2 .
Let us now analyse the cumulative performance of the analysed potentials for all SL MoS 2 phases. We see that ∑MAPE in Table 4 is the lowest, and almost the same, for three potentials: SW2017, SW2016 and REBO. However, only the REBO potential distinguishes three different phases, the other two potentials degenerate phases, i.e., instead of three they produce two.

Phonon Spectra
Phonon spectra along the appropriate high symmetry q-points [37], calculated by applying the PBEsol XC functional, for SL 1H-MoS 2 phase are depicted in Figure 4a), for SL 1T-MoS 2 phase are depicted in Figure 4b) and for SL 1T'-MoS 2 phase are shown in Figure 4c), respectively. Experimental data for single-layer molybdenum disulphide are very scarce and concern only Γ point in 1H-MoS 2 phase, see [58]. When we compare the results obtained here with those calculated by other authors, we see agreement typical for different DFT calculations, see [31,[59][60][61][62].
Analysis of the computed curves in Figure 4a-c allows us to conclude that phases 1H-MoS 2 and 1T'-MoS 2 are not only mechanically but also dynamically stable, i.e., phonon modes everywhere have positive frequencies. Phase 1T-MoS 2 is mechanically stable, but not dynamically stable, i.e., phonon modes also have negative frequencies. Similar observations can be found in [31,60]. Let us now compare the phonon spectra for SL MoS 2 phases calculated with DFT, Figure 4, and those calculated with LAMMPS and the three best potentials, i.e., SW2017, Figure 5, SW2016, Figure 6 and REBO, Figure 7. As a result that only the REBO potential distinguishes three different phases of SL MoS 2 , the molecular calculations of phonons utilise basic cells derived from DFT calculations. At first glance, we can see that it only makes sense to compare molecular phonons with DFT phonons merely qualitatively, not quantitatively. All three potentials are qualitatively well reproducing the phonon spectra for SL 1H-MoS 2 phase, see

Conclusions
A systematic quantitative comparison of Stillinger-Weber, REBO, SNAP and ReaxFF potentials for the reproduction of the properties of 2D molybdenum disulphide polymorphs was presented. To compare the potentials, the structural and mechanical properties and phonon dispersion of single-layer phases 1H, 1T and 1T' MoS 2 (SL MoS 2 ) obtained from the functional density theory (DFT) and molecular static (MS) calculations were used.
We can conclude that: • The transferability of analysed molecular potentials leaves much to be desired. To increase the transferability of potentials, the number of configurations to be taken into account in the parameter optimisation process should be significantly increased.