Revealing the role of ripples in phonon modes for MoS₂ and MoSe₂: insights from molecular dynamics and machine learning

Abstract

Interatomic potentials for single-layer MoS₂ and MoSe₂ were developed by training an artificial neural network with a reference data set generated using density functional theory. High accuracy was obtained for the phonon dispersion as well as a direct correspondence of the predicted Raman and infrared (IR) active modes (between 100 cm⁻¹ and 600 cm⁻¹) with molecular dynamics results. The ability to perform simulations with thousands of atoms allowed our system to accommodate spontaneous fluctuations in height and ripples propagating in the material. As a result of these ripples, our simulations suggest the existence of IR activity in the E′ and modes with inverted polarization, as well as the presence of IR activity at frequencies close to E′′ and , which are generally considered IR-inactive. The rippling of the system is expected to be responsible for breaking the 3-fold rotational symmetry around the z-axis of the material, resulting in non-trivial IR modes.

---

1 Introduction

Transition metal dichalcogenides (TMDs) have the stoichiometry MX₂ with M being a transition metal (e.g., M = Mo, W…) and X belonging to the oxygen family (e.g., X = S, Se, Te). Single-layer (SL) TMDs present unique behavior,^1,2 with potential optoelectronic applications. Notably, SL-MoS₂ is known for having a direct bandgap suitable for light absorption in the visible range of 1.88 eV while its bulk counterpart has an indirect bandgap of 1.29 eV.³ Furthermore, excitonic effects with strong binding energies⁴ play a major role in its optical properties. This system with its unique optical response can be leveraged in nanoscale transistors⁵ and several photoactive devices. MoS₂ is also known to have promising thermal properties due to carrier transport.⁶ Similarly, MoSe₂ is known to share such optical properties with MoS₂ by having a direct bandgap of 1.63 ± 0.01 eV.⁷

Molecular dynamics (MD) simulation is a well-known tool for studying structural and dynamical properties of TMDs that which often strongly depend on the temperature of the system^8–10 such as thermal conductivity and expansion coefficients which can be compared with experimental results. However, one of the main limitations of using MD simulations to explain real experiments is the accuracy of the interaction potentials/force fields between atoms in the simulation. Potentials with a complex functional form are usually developed in order to close the gap between experiment and simulation, e.g., Tersoff for carbon,¹¹ TIP3P for water,^12–14 Vashishta for Silica¹⁵ and the Embedded Atom Method for metallic systems.¹⁶ A widely employed potential for MoS₂ is the Stillinger–Weber proposed by Jiang et al.^17,18 to reproduce phonon dispersion but the breaking and formation of bonds is better described by the ReaxFF potential,¹⁹ which makes the work of simulating SL-MoS₂ a challenging prospect as one needs to select which potential is most suitable for the desired property, an indication of the non-transferability of these potentials such that their lack of accuracy can raise the aforementioned problems. Ab initio molecular dynamics (aiMD) could theoretically solve this but the calculations are so computationally expensive that only systems with at most a few hundred atoms could be explored. One way to overcome these issues is to use Machine Learning (ML) methods to develop potentials based on data generated using Density Functional Theory (DFT), also known as Machine Learning Force Fields (MLFFs).²⁰ The main idea is to generate a reference data set with DFT calculations of several structures and train a machine learning model (e.g. kernel based or neural networks) to reproduce its results within the desired accuracy. After training, the model will not have a specific functional form like the conventional force fields but will be able to reflect the DFT effective motion patterns with an arbitrary accuracy and the reliability of the ML model will be related with the structures included in the data set. Therefore, simulations using MLFFs will have an accuracy similar to aiMD with much less computational cost. To some extent, lattice dynamics offer a solid theoretical background to several optical experiments such as infrared (IR) terahertz (THz) absorption as well as Raman scattering.²¹ On a crystalline system, the lattice dynamics of a system is well described by its phonon dispersion and a common criterion to test the reliability of a force field (conventional or ML based) is to calculate its phonon dispersion and compare with ab initio calculations or experimental results.²²

Moreover, because SL-MoS₂ and MoSe₂ are ordered periodic 2D systems, one would expect them to be thermodynamically unstable at finite temperatures according to the Mermin Theorem,²³ and that thermal fluctuations would inhibit any kind of long-range order necessary in a 2D crystal. However, subsequent studies have shown that such fluctuations can be suppressed by an anharmonic coupling between the “bend” and “stretch” acoustic modes^24,25 and that such materials could be stable at finite temperature as long as they present large ripples. In graphene such ripples tend to be between 50 and 100 Å in lateral size and around 5 to 10 Å in height²⁶ while in MoS₂ those numbers vary between 60 to 100 Å and 6 to 10 Å respectively.²⁷ Also, ripples have an important role on the thermal expansion proprieties²⁸ and several efforts are underway to determine the structure of ripples through image reconstruction.²⁹ Here, we present two MLFFs developed for SL-MoS₂ and MoSe₂, we validate our models by reproducing experimental dielectric spectra results at low frequencies such as THz absorption with MD simulations and present results on how the ripples propagating through the system affect the spectra according to our MLFFs. In our simulations, we observed that these ripples are responsible for new peaks in the dielectric spectra beyond the usual peaks pertaining to E′ and modes, including ones next to phonon modes known to be IR inactive. This paper is structured as follows: in Section 2 we detail all models and methods used in this study, Section 3 show our results and have the main discussion, and we conclude in Section 4. Section 2 will be divided into three different subsections: in 2.1 we will describe the DFT calculations used to generate the reference data set, on 2.2 we will show the machine learning method selected to produce the MLFF and on 2.3 we will present the MD simulations and data analysis used on this study. Furthermore, Section 3 will also be divided into three subsections: the discussions about the development, accuracy and testing of the MLFFs for MoS₂ and MoSe₂ will be on Sections 3.1 and 3.2 respectively while Section 3.3 contains all the discussions about the role of thermal ripples on the THz spectra and the phonon modes for both compounds.

2 Methodology

2.1 DFT calculations & data set

In order to generate our data set, we performed self-consistent first-principle calculations using the plane wave pseudopotential method within the density functional formalism as implemented in the Quantum Espresso code.^30–32 We used optimized norm-conserving vanderbilt (ONCV) pseudopotentials^33,34 to treat the interaction between the atomic core and the valence electrons and the PBE generalized gradient approximation for the exchange–correlation potentials.³⁵ A Monkhorst–Pack grid of 6 × 6 × 1 was used for the integration on the Brillouin-zone and a kinetic energy cutoff of 70 Ry and 280 Ry was used for the wavefunctions and charge density respectively (see Fig. S1 in ESI†). The hexagonal primitive cell (see: red diamond on Fig. 1a) relaxation was achieved when the total force on the unit cell was less than 5.0 × 10⁻³ eV Å⁻¹. A vacuum region of 20 Å is set to ensure negligible spurious interactions between neighboring replicas. As a result of structure relaxation, the optimized lattice constant (indicated as a₀ on Fig. 1a) of SL-MoS₂ was 3.186 Å, which has a 0.79% error with respect to the experimental value of 3.16 Å.³⁶ The Mo–S relaxed distance (shown as d_Mo–X at Fig. 1b) was calculated to be 2.414 Å and the S–Mo–S angle as 80.76° (shown as θ at Fig. 1b). For MoSe₂, we adopted a cutoff radius of 40 Ry for the wave functions and 160 Ry for the charge density (see Fig. S1 in ESI†). The optimized lattice constant of SL-MoSe₂ was 3.32 Å, which has a 0.91% error with respect to the experimental value of 3.29 Å.³⁷ The Mo–Se relaxed distance was calculated to be 2.541 Å and the Se–Mo–Se angle was calculated to be of 82.12°. Those are the main parameters of the cell optimization and small displacements around those equilibrium positions are expected to have restoring forces.

Fig. 1 Top (a) and side (b) views of SL MoX₂ (X = S, Se). Mo is indicated in purple while X is in yellow. The red diamond indicates the hexagonal primitive cell with the lattice parameter a₀ indicated as an arrow. The green rectangle shows the rectangular unit cell with the arrows indicating the parameters |a₁| = |a₀| and . The distance between Mo and X is shown by d_Mo–X and the angle between X–Mo–X is shown by θ.

The initial set of reference structures for the MLFF construction was generated by systematically distorting the ideal structure. The lattice constants of the crystal structures were scaled within a range of ±45% around the equilibrium value on the xy and z directions‡. To get a sample of the atomic forces, for each pair (xy,z) of scale percentages we ran 3 to 5 different DFT calculations with atomic positions randomly displaced about their equilibrium positions with displacement amplitudes in the range of 0.01–0.1 Å. This initial set is used to generate a preliminary MLFF which is used to calculate the phonon dispersion and run a preliminary MD simulation to check the stability of the system. As we proceed to satisfy both criterion, additional structures were added to the initial data set. We shift atoms along lines corresponding to normal modes on various points of the Brillouin zone to ensure that our database contains sufficient information to reproduce the phonon dispersion, especially in the vicinity of the Γ point where the preliminary calculations often presented dynamical instabilities. This data collection method is known as normal mode sampling²⁰ and the overall workflow is illustrated at Fig. 2. In total, the reference data set comprises a total of 5259 structures for MoS₂ and 6187 for MoSe₂, each of which with has the number of atoms in the unit cell ranging from 3 to 48.

Fig. 2 Flowchart of the adopted methodology for construction, improvement and utilization of the machine learning force fields (MLFFs) generated in this work.

2.2 Machine learning methods

Here we used an Artificial Neural Network (ANN) as implemented in the Atomistic Energy Network (ænet)^38,39 package based on descriptors represented on the Chebyshev basis.⁴⁰ This package has tools for generating, training and testing ANNs. In our model, we considered interactions up to a certain cutoff radius r_cut, which are represented by an ANN as shown schematically in Fig. 3. The ANN returns the potential energy E_ξ(r) of one atom in the system (atom ξ) based on the position of itself and its neighbors within the cutoff. These positions must be represented in a set of coordinates {x₀ⁱ} that is invariant with respect to translation, rotation and exchange of atoms, which are the descriptors indicated at Fig. 3 and are used as the input layer for the MLFF.

Fig. 3 Schematic representation of a MLFF that consists of an input layer, three hidden layers and an output layer. The configuration of neighbors around atom ξ within a cutoff r_cut is transformed to a descriptor x₀ⁱ (i = 1,…,6) and used as the input to the MLFF. Arrows connecting nodes in adjacent layers represents weights. Biases and activation functions are not shown for simplicity. The selection of number of nodes and layers in this figure are only for illustration purposes.

In addition to the parameters that are determined during the training process of a machine learning model, e.g., the weights w^i,j_I in ANNs present in Fig. 3, most of models contain hyperparameters that are chosen before training and do not change during it. Their choice can affect the accuracy of model and its behavior during simulations since a model with too few parameters may not capture all desired features in the database. For our ANNs, there are two different types of hyperparameters. The first set is related to the representation of the local structure environment of the system, also known as “fingerprints”, which must be defined for both Mo and S (Se) atom types. We have used the same hyperparameters to define the descriptors for all atom types: a radial r_cut of 9.0 Å, a radial order expansion of 16, an angular r_cut of 8.0 Å and an angular order expansion of 8. The second set of hyperparameters is the number of hidden layers that connects the input layer to the output layer as well as the number of nodes in each. Moreover, an activation function is defined between two successive layers. Here, for MoS₂ we built an ANN composed by 2 layers with 24 nodes followed by 2 layers with 20 nodes (a total of 4 hidden layers), using hyperbolic tangent as the activation function. For MoSe₂ we built an ANN composed by 2 layers with 20 nodes followed by 2 layers with 16 nodes and the same activation function. We also tested neural networks with 2 and 3 hidden layers. However, none of them were able to reproduce the desired dynamical properties of the studied systems. Of the diverse number of nodes and layers that were tried, those architectures showed the best results. On the training process, the cost function chosen to be minimized is the root mean squared error (RMSE):

with N_Ref being the number of reference structures and σ indicating each structure in the data set. We employed the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method for weight optimization.^41,42 The reference data set was split into a training set of 80% and a testing set of 20% to check the prediction capability. A test of the accuracy and validation of the MLFFs can be found in Section 4.2 of ESI.†

2.3 Molecular dynamics simulations

All MD simulations in this study were made using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)⁴³ with the assistance of a module that allows Neural Networks generated with the æ net package to be used on LAMMPS simulations as “pair styles”.⁴⁴ The procedure to obtain the complex dielectric spectra at THz frequencies is based on expressing the dielectric spectra as a Fourier transform of the average total dipole moment correlation function on time domain:⁴⁵where M = M(t) indicates total dipole moment of system, ε₀ is the vacuum permittivity, V is the volume and ω the angular frequency. For these simulations, we used a timestep of dt = 1 fs within the following scheme: the first 250 000 steps were performed with the Nosé–Hoover Thermostat (NVT ensemble) at a chosen temperature to equilibrate the system and, after that, the thermostat is turned off and the simulations were carried out in the NVE ensemble, which is required to obtain dynamical properties. Since the Nosé–Hoover Thermostat is based on a extended Lagrangian responsible for coupling the system to an external heat-bath,⁴⁶ oscillations of the coupled heat-bath could affect the correlation functions needed in eqn (2). By switching from the NVT to the NVE ensemble, with the appropriate amount of equilibration steps, we were able to prevent these “virtual” oscillations from influencing the final result. After 8000 steps to equilibrate the system in the NVE ensemble, the subsequent 3 992 000 additional NVE steps were conducted to obtain individual M(t)·M(0) functions where t ranges from 0 to 22 ps. Then the overlap approach was used to increase correlation between successive samples, reducing the noise in the Fourier transform.^10,47 These individual functions were averaged to obtain the main correlation function of dipole moment which was used to obtain the dielectric spectra of the material. To obtain the error bars in our calculations, all M(t)·M(0) functions were partitioned in 10 different groups, each one of them giving one estimate to the dielectric function. Then, those estimates were averaged for the final results. The simulation cell was a rectangular supercell (see green rectangle at Fig. 1) with a total of 16 320 atoms and dimensions 68 a₁ × 40 a₂ (a₁ = a₀ and ). This size was adopted to minimize finite size effects and simulations with larger lattices did not improve the results.

3 Results & discussion

3.1 MoS₂

3.1.1 Phonon properties. Fig. 4 shows the results for the phonon dispersion calculated with the developed MLFF for MoS₂ in this study. The comparison with DFT in Fig. 4a has the purpose of verify how accurately the MLFF reproduces the data on the reference set. We can see that all the main features of the DFT dispersion are reproduced by our MLFF for both the optical and acoustic branches, indicating that our model “learned” the main motion patterns predicted by DFT. Fig. 4b shows the same MLFF phonon dispersion as Fig. 4a but compares it with Stillinger–Weber (referred in this section as SW-MoS₂ the parametrization on ref. 18) and with inelastic X-ray Scattering measurements⁴⁸ shown as the black dots. Both MLFF and Stillinger–Weber (SW) dispersion curves were obtained with the assistance of the Alamode package.⁴⁹ The first thing to be noticed here is that, although MLFF had a similar accuracy to SW-MoS₂ for the acoustic branches when compared with experimental data, MLFF follows the experimental results for optical branches closer than SW-MoS₂. When studying the lattice dynamics of a given system, special attention must be paid to the Γ modes, as they can be determined experimentally by IR and Raman measurements. Table 1 shows the predicted values for the four modes at Γ, comparing it with SW-MoS₂ and experimental results, each mode being drawn in Fig. 5. It is worth mentioning that of all the modes, only E′ and are IR active, since their oscillation generates an instantaneous dipole moment and produces instantaneous polarization. All the modes are known to be Raman active except , and both E′ and being used widely in literature to characterize MoS₂ samples.

Fig. 4 Phonon dispersion of MoS₂ along high symmetry directions Γ–M–K–Γ for MLFF, DFT calculations, Stillinger–Weber (SW-MoS₂) potential estimation and Inelastic X-ray scattering data.

Table 1 Phonon modes at Γ of SL-MoS₂. The unit is (cm⁻¹)

Mode	MLFF	SW-MoS2	Exp. results
a Ref. 50. b Ref. 51. c Ref. 48. d Ref. 52. e Ref. 53.
E′′	291.5	267.3	286^c
E′	372.3	360.5	384.3^a, 384.0^b, 383^c, 383.7^e
	390.8	429.2	403.1^a, 408^c
	459.0	506.8	470.0^b, 466^d, 468.2^e

---

Fig. 5 Sketch of the optical phonon modes of SL MoX₂ (X = S, Se) at Γ. Mo and X atoms were colored with turquoise and orange respectively.

3.1.2 THz dielectric spectra. Fig. 6 shows the MD results for the dielectric spectra for MoS₂ at 300 K at THz frequencies (from 100 cm⁻¹ to 600 cm⁻¹). The spectra are divided into two parts: in-plane (considering only the components of the dipole moment that are parallel to the layer) and out-of-plane (considering only the component of dipole moment that is perpendicular to the layer), represented by the blue and red curves respectively. Fig. 6b has two narrow absorption peaks, hence by fitting each peak to a Lorentzian oscillator,²¹ the peak positions were calculated to be at 368.3 ± 0.1 cm⁻¹ and 452.8 ± 0.3 cm⁻¹ respectively. Fig. 5 shows that the eigenvectors of the E′ mode only alter the dipole moment in the planar directions while eigenvectors of the only alter the dipole moment in the perpendicular ones. Thus, the peak illustrated by the blue curve can be compared to the E′ value in Table 1 while the peak denoted by the red curve can be compared with . Its noteworthy the proximity of the calculated values in Fig. 6b and Table 1, with the presented redshift being attributed to temperature effects of the MD simulation. We can compare our simulation results with IR reflectivity measurements by Wieting et al.⁵⁴ (more specifically with Fig. 6 of the paper). The same narrow peak found in our simulations are present in the reflectivity measurements while their peak positions are located at 384 cm⁻¹ for the in-plane mode and 472 cm⁻¹ for the out-of-plane mode, being very close to the values obtained by the simulations using our MLFF.

Fig. 6 Real (a) and imaginary (b) part of the dielectric function for MoS₂ as a result of our MD simulations at 300 K.

Fig. 7 shows our MD results for the THz dielectric spectra of MoS₂ obtained with the MLFF at 150 K, 350 K and 550 K. From it, we can discern several trends: with increasing temperature, both peaks in Fig. 7c and d broaden and redshift, a behavior in accordance with experiments.⁵⁵Fig. 8 shows the temperature dependence of the E′ and modes based on the peak positions of the imaginary part of the spectra. The numerical values were obtained by fitting Lorentzian oscillators to the peaks presented in the data. To obtain the SW-MoS₂ data points, the same simulation procedure described in Section 2.3 was applied. Both simulations with MLFF and SW-MoS₂ produced a redshift with increasing of temperature. For each mode and for each potential, we performed linear fits (represented by dashed lines) and extracted the slope, which contains information about the average redshift of the system due to temperature. Table 2 shows results for the slopes and compares it to experimental results, indicating that our calculations are consistent with literature. Fig. 9 shows our MD results for the THz dielectric spectra for MoS₂ obtained with the MLFF at 300 K using a unstrained lattice (0.0%), a lattice with 0.8% strain of the original and a lattice with 1.8% strain of the original. Qualitatively, the frequency peaks in Fig. 9(c and d) become redshifted with the application of strain, similarly to the reported temperature effects. However, for more strained lattices the peaks are narrower than for unstrained ones. The redshift can be attributed to the weakening of the bonds caused by the tensile strain⁵⁶ while the narrowing of the peaks can be attributed to the fact that in a strained lattice, acoustic modes do not propagate as freely as in an unstrained lattice and, therefore, cannot interfere with the propagation of optical modes. Fig. 10 shows the strain dependence of the of the E′ and modes. Similar to Fig. 8, we performed linear fits and extracted the slope for each mode and for each potential. Table 3 shows the results for the temperature coefficients and compares it to experimental and other theoretical results. We notice a qualitative difference between the behavior of the mode for simulations with MLFF and SW-MoS₂, the latter presenting a slight blueshift with the application of strain, a result significantly different from the DFT prediction. For both modes, the MLFF reproduces DFT results more closely than SW-MoS₂. When comparing with experimental results, both potentials had similar performance for the E′ mode. However, the diverse range of results for the E′ mode prevents us from making any claims about the accuracy of our calculations regarding the temperature coefficient. Moreover, since the mode is not Raman active, to our knowledge, there are not any experimental results that studies the variation of that mode with strain. In that sense, our results are a prediction of the mode redshift including temperature effects.

Fig. 7 THz dielectric spectra for MoS₂ at 150 K, 350 K and 550 K. Results in Fig. 8 for E′ and modes are related to the peaks in (c) and (d) respectively. We included fewer temperatures than in Fig. 8 for clarity.

Fig. 8 Temperature dependence of IR-active modes at MoS₂ for MLFF and SW-MoS₂ with error bars. Dashed lines indicates the linear fit performed to obtain the coefficients at Table 2.

Table 2 Temperature coefficient for E′ and modes of MoS₂ (10⁻² cm⁻¹ K⁻¹)

Mode	MLFF	SW-MoS2	Exp. results
a Ref. 55. b Ref. 57. c Ref. 58. d Ref. 59. e Ref. 60. f Ref. 61. g Ref. 62. h Ref. 63.
E′	−2.2	−2.2	−1.32^a, −1.47^b, −1.41^c, −1.79^d, −1.61^e, −1.36^f, −1.6^g, −1.3^h
	−3.1	−2.7	−1.99^e

---

Fig. 9 THz dielectric spectra for MoS₂ at 300 K for the strain percentages 0.0%, 0.8% and 1.4%. Results in Fig. 10 for E′ and modes are related to the peaks in (c) and (d) respectively. We included fewer strains than in Fig. 10 for the sake of clarity.

Fig. 10 Strain dependence of IR-active phonon modes at MoS₂ for MLFF and SW-MoS₂ with error bars. Dashed lines indicates the linear fit performed to obtain the coefficients at Table 3.

Table 3 Strain coefficient for E′ and modes of MoS₂ (cm⁻¹/%)

Mode	MLFF	SW-MoS2	Exp. results	DFT
a Ref. 64. b Ref. 65. c Ref. 66. d Ref. 67. e Ref. 68. f Spin Coated PVA encapsulation. g Conventional exfoliation.
E′	−3.49	−4.98	−7.4^a^f, −3.2^a^g, −4.5^b, −2.1^c, −4.67^d	−4.28^e
	−4.87	0.26		−4.19^e

---

3.2 MoSe₂

3.2.1 Phonon properties. Fig. 11 shows the results for the phonon dispersion calculated with the developed MLFF for MoSe₂ in this study. The purpose for the comparison with DFT had already been described in Section 3.1.1 but, to our knowledge, there are not any experimental results in the literature on the modes for MoSe₂ at several points in the Brillouin zone (such as in ref. 48). We compare our results with Stillinger–Weber (referred in this section as SW-MoSe₂) as described in the parametrization proposed by Kandemir et al.⁶⁹ For the optical phonon branches, we observe that the MLFF fits the DFT dispersion more closely than SW-MoSe₂ as expected. Table 4 shows the MLFF predicted values for the four modes at Γ, comparing it to SW-MoSe₂ and experimental results. From Table 4 we observe that SW-MoSe₂ predicts frequencies that are considerably redshifted when compared to experiment and that the MLFF predicts more accurate values. It is worth mentioning that the E′ mode has a lower frequency than the mode in MoS₂ while the opposite is observed in MoSe₂.

Fig. 11 Phonon dispersion for MoSe₂ along high symmetry directions Γ–M–K–Γ for the adjusted MLFF, DFT calculations and Stillinger–Weber (SW-MoSe₂) potential.

Table 4 Phonon modes at Γ of SL-MoSe₂. The unit is (cm⁻¹)

Mode	MLFF	SW-MoSe2	Exp. results
a Ref. 70. b Ref. 71. c Ref. 72. d Ref. 73. e Ref. 74. f Ref. 75. g Ref. 76.
E′′	165.1	153.8	168^a, 170^b, 175.8^f
	234.9	213.8	242^a, 243^b, 241.2^c
E′	275.3	252.2	286.9^d, 289^a, 287.3^c, 288^e, 286^f, 283^g
	342.7	311.5	352^a, 350^e, 352^f

---

3.2.2 THz dielectric spectra. Fig. 12 shows the MD results for the MoSe₂ THz dielectric spectra at long wavelengths similar to Fig. 6 in Section 3.1.1. In Fig. 12b, the peak for the blue and red curves respectively are localized at 269.8 ± 0.1 cm⁻¹ and 334.9 ± 0.1 cm⁻¹ respectively, which can be compared directly with the values from Table 4 for the E′ and for the same reasons given before. Similar to what was observed on MoS₂, we can assign those modes to the peaks, and the difference of the peak positions and the tabled values also can be attributed to a redshift effect due to temperature in our simulations. We can compare our results to infrared reflectance results from Sekine et. al.⁷⁰ (see Fig. 6 of the referred work), as they assigned the values 289 cm⁻¹ and 352 cm⁻¹ to the modes E′ and respectively. Fig. 13 shows our MD results for the THz dielectric spectra of MoSe₂ obtained with the MLFF at 150 K, 300 K and 450 K. The same qualitative trends observed on our MLFF simulations for MoS₂ can be observed on MoSe₂, i.e., the redshift and broadening of the peaks in Fig. 13c and d. However, the broadening is less pronounced than those in Fig. 7c and d. Fig. 14 and Table 5 shows our results for the temperature dependence of the infrared active phonon modes, comparing it with simulations with SW-MoSe₂ and experimental results. We notice that, when comparing with experiments, the simulations with our MLFF predicts more accurate coefficients than SW-MoSe₂ and, although the MLFF prediction for the temperature coefficient of the mode is considerably different than the experimental values, it is still a better prediction than those given by SW-MoSe₂. Also, SW-MoSe₂ predicts similar temperature coefficients for the two distinct modes while in the MLFF simulations the temperature coefficient seems to be sensitive to the oscillation mode. Fig. 15 shows results for the THz dielectric spectra of MoSe₂ at 300 K for the unstrained lattice (0.0%), and lattices with 0.8% and 1.8% strain respectively. Although a redshift of the peak is noticed, we did not observe any noticeable changes on the breadth of the peaks. Fig. 16 and Table 6 shows our results for the strain dependence of the infrared active phonon modes. The lack of other reported experimental or theoretical results made it difficult to rate the performance our MLFF. To our knowledge, the references cited at Table 6 were the only studies that reported the strain coefficients for the E′ mode and there are no reports for the mode. When comparing to these reported data, neither our MLFF nor SW-MoSe₂ presented accurate results.

Fig. 12 Real (a) and imaginary (b) part of the dielectric function for MoSe₂ as a result of our MD simulations at 300 K.

Fig. 13 THz dielectric spectra for MoSe₂ at 150 K, 300 K and 450 K. Results in Fig. 14 for E′ and modes are related to the peaks in (c) and (d) respectively. We included fewer temperatures than in Fig. 14 for the sake of clarity.

Fig. 14 Temperature dependence of IR-active phonon modes in MoSe₂ for MLFF and SW-MoSe₂. Dashed lines indicates the linear fit performed to obtain the coefficients in Table 5.

Table 5 Temperature coefficient for E′ and modes of MoSe₂ (10⁻² cm⁻¹ K⁻¹)

Mode	MLFF	SW-MoSe2	Exp. results
a Ref. 77. b Ref. 58. c Ref. 78. d Ref. 79.
E′	−0.79	−3.9	−1.18,^a −1.21,^b −0.69^c
	−2.1	−3.7	−0.64,^c −0.86^d

---

Fig. 15 THz dielectric spectra for MoSe₂ for the strain percentages 0.0%, 0.8% and 1.4%, Results in Fig. 16 for E′ and modes are related to the peaks in (c) and (d) respectively. We included fewer temperatures than in Fig. 16 for the sake of clarity.

Fig. 16 Strain dependence of IR-active phonon modes at MoSe₂ for MLFF and SW-MoSe₂. Dashed lines indicates the linear fit performed to obtain the coefficients at Table 6.

Table 6 Strain coefficient for E′ and modes of MoSe₂ (cm⁻¹/%)

Mode	MLFF	SW-MoSe2	Exp. results	DFT
a Ref. 80. b Ref. 81.
E′	−0.95	−1.59	−4.93^a	−4.28^b
	−2.54	−5.6

---

3.3 Rippling effects

Fig. 17 illustrates the imaginary part of the dielectric spectra for both systems (akin to Fig. 6(b) and 12(b)), but specifically focusing on regions with low absorption intensities, with a version of the same data plotted on a semi-log scale that is included in Fig. S7 in Section 4.3 of ESI.† Within the figure, it is evident that each system exhibits four distinct new peaks. Noticeably, two of these peaks align precisely with the frequencies of the E′ and modes, however with reversed polarization directions. Additionally, MoS₂ displays peaks at approximately 310 cm⁻¹ and 393 cm⁻¹, whereas MoSe₂ features peaks around 163 cm⁻¹ and 225 cm⁻¹, values closely matching the frequencies of the E′′ and modes as listed in Tables 1 and 4, respectively. Observing these new peaks allows us to identify two key characteristics in our simulations. Firstly, the E′ and modes exhibit infrared (IR) activity in a direction perpendicular to their usual oscillation. Secondly, there is a manifestation of some IR activity, although minimal, by the E′′ and modes, traditionally considered IR inactive. Moreover, we noticed that all of the newly observed peaks were extinguished under the application of tensile strain (see Fig. S8 in ESI†).

Fig. 17 Imaginary part of dielectric spectra for MoS₂ and MoSe₂ at 300 K with insets in the small intensities region. Arrows indicate peaks centered at 310(163), 368(225), 393(271) and 452(334) for MoS₂ (MoSe₂) in cm⁻¹.

Applying strain to the lattice reduces fluctuations in its height, thereby diminishing the amplitude of ripples, as depicted in Fig. 18. This figure demonstrates that in both systems, the average height (defined as the difference in the z coordinate between the highest and lowest Mo atom, as illustrated in the figure) decreases with increasing strain. For MoS₂ the average height decreases from 9 Å to 2 Å while for MoSe₂ it decreases from 6 Å to 2 Å when a strain 2% above the equilibrium lattice constant is applied. Therefore, considering the data presented in Fig. 17 and 18, it can be inferred that the observed low-intensity peaks observed correlate with the ripples present in SL-MoS₂ and MoSe₂. In the vicinity of the peaks corresponding to the frequencies adjacent to the E′ and modes, the mechanism by which ripples induce IR activity with altered polarizations becomes clearer. By allowing the lattice to ripple, the E′ mode is no longer confined strictly to the plane, and similarly, the mode is not restricted solely to the axis perpendicular to the plane. This results in the E′ mode acquiring an out-of-plane component and the mode acquiring a planar component, as depicted in Fig. 19. However, this explanation alone does not suffice to elucidate the observed IR activity near the frequencies associated with the E′′ and modes in our simulations. To gain further insight into the influence of ripples on the dielectric spectra and to explore potential explanations for the peaks observed at E′′ and , we constructed a supercell composed of 1020 atoms. Within this supercell, we introduced a sinusoidal wave characterized by a wave vector k parallel to the x̂-axis, with a wavelength determined by the total length of the cell in the x-direction. This supercell is depicted in Fig. S9 of the ESI† and maintains translational symmetry along the y direction and this simplified model for the ripple was used to qualitatively investigate the origins of the new peaks observed in the spectra. Using the MLFFs to model atomic interactions, we employed the Alamode package to compute the phonons of the primitive cell (see Fig. S10 at ESI†) corresponding to the rippled supercell at the Γ point. The calculations also provided the eigenvectors U_j(κ) associated with each frequency eigenmode ω_j. These eigenvectors enable the determination of the IR intensities for each mode, calculated using the equations defined in eqn (2.112) and (2.114) of ref. 82:

where e_κ and M_κ denote the charge and mass of atom κ, respectively§. An oscillation mode j contributes to Im[ε(ω)] only if f(ω_j) is non-zero, serving as a selection rule for IR absorption.^83,84 Consequently, it becomes feasible to attribute each peak in Fig. 17 to those specific oscillation modes with non-zero f(ω_j).

Fig. 18 Mean lattice height as function of strain for MoS₂ and MoSe₂ during a simulation. The sketch of the rippled single-layer serves to outline the definition of the height Δz of the lattice.

Fig. 19 Sketch of the decomposition of the components of the oscillation modes E′ and in a rippled lattice.

The IR intensities results for MoS₂ and MoSe₂ are shown in Fig. 20 and 21, respectively and the calculated frequencies for all modes are listed in Tables S1 and S2 in ESI† In each, the calculation for the unrippled (flat) cell was included at Fig. 20a and 21a for comparison purposes and the insets were included to facilitate visualization of all modes with low intensity f(ω_j). As expected for the flat cell in both compounds, only the E′ and modes showed a non-zero IR response. The E′ mode presents an in-plane response only while the mode presents an out-of-plane response only. For the rippled cell, a degeneracy breaking in frequencies is observed due to the breaking of some symmetries of the system caused by the introduction of the ripple. This symmetry breaking causes the only two modes present in the flat cell to split into several bars located at frequencies close to the E′ and modes in the flat cell. However, it is still expected that the new modes introduced by the ripple must have some similarity to those of the original cell due to the fact that they share similar frequencies. In both compounds the two largest bars, one blue and one red, correspond to the vibrations that best represent the E′ and modes of the flat cell in the rippled cell, the characteristic vibration of the E′ in the planar directions and the characteristic vibration of the in the out-of-plane direction. These bars would correspond to the main peaks depicted in Fig. 6b and 12b. Looking at the low IR intensity modes, it is noticeable that in both Fig. 20b and 21b the presence of a mode with a frequency close to E′ that presents out-of-plane IR activity (indicated by the red bar) as well as a mode with a frequency close to that presents in-plane IR activity (indicated by the blue bar) in both figures. Furthermore, both figures indicate the presence of modes with IR activity at frequencies close to those of the E′′ and modes which, in the flat cell, are known to be IR inactive. Thus, the degeneracy breaking caused by the introduction of the ripple is also responsible for the splitting of the E′′ and modes into several modes, some of which have some IR activity. The results shown in Fig. 20 and 21 strengths the hypothesis that lattice ripples would be responsible for the appearance of new peaks present in the simulations and shown in Fig. 17.

Fig. 20 IR intensities for the MoS₂. (a) and (b) indicate the calculated intensities for the cell without and with the ripple, respectively. The dotted lines represent the frequencies listed in Table 1.

Fig. 21 IR intensities for the MoSe₂. (a) and (b) indicate the calculated intensities for the cell without and with the ripple, respectively. The dotted lines represent the frequencies listed in Table 4.

The main advantage of calculating the eigenvectors of the ripples lattice lies in the ability to observe the patterns of new vibrational modes and understand how these modes affect the total electric dipole moment of the cell, which will be shown on the following figures. Here, we denote the modes of the rippled cell with the same notation of Tables 1 and 4 to indicate the mode of the flat cell that generated the new mode, but we will use a “*” to indicate that this mode belongs to the rippled cell and we will indicate the direction of polarization of the mode in parenthesis. Fig. 22 shows the modes *E′(z) and , found to have IR activity next to E′ and with inverted polarizations at Fig. 20(b-ii and iv) and 21(b-iii and iv). Upon comparing the modes depicted in Fig. 5 and 22, a discernible similarity emerges between modes *E′(z) and in relation to modes E′ and respectively. This similarity arises due to certain MoS(Se)₂ units within the rippled supercell exhibiting identical vibrational pattern as those observed in the primitive flat cell, leading to similar oscillation frequencies in these modes. However, Fig. 22 reveals how the oscillation of the whole rippled cell lead to alterations in its total dipole and in which direction. In mode *E′(z), the cell exhibits a oscillation pattern resembling E′, with atoms on either side of the ripple crest moving in opposite directions, which results in the cancellation of the horizontal component of total dipole moment of the mode, leaving only a resultant vertical component. Thus, mode *E′(z) manifests an out-of-plane infrared (IR) response with a frequency close to E′, consistent with the results shown in Fig. 17. Similarly, in mode , the cell displays an oscillation pattern akin to , also with atoms on either side of the ripple crest moving oppositely similar to *E′(z). This leads to the cancellation of the vertical component, leaving predominantly only a resultant horizontal component. Consequently, mode exhibits an in-plane IR response with a frequency close to , which is also consistent with the results in Fig. 17. The analysis of modes *E′(z) and on Fig. 22 strengths the hypothesis depicted in Fig. 19, wherein the rippled lattice causes modes E′ and to acquire out-of-plane and in-plane components, respectively. Fig. 23 shows the modes *E′′(x) and *E′′(y) related to the two highest blue bars at Fig. 20(b-i) for MoS₂ and Fig. 21(b-i) for MoSe₂, both of them presenting IR activity next to E′′. These modes exhibit an oscillation pattern similar to E′′, where atoms on opposite sides of the crest move in opposite directions. However, both modes appear to curve next to the ripple, with the *E′′(x) curving horizontally and the *E′′(y) curving vertically. Thus, the oscillation pattern depicted in Fig. 23, as predicted by the MLFFs, suggests that the ripple induces curvature in the displacement vector fields of S(Se) atoms, particularly near its inflection point. This curvature, visible on the top view of the sketch of both modes, is responsible for the emergence of a dipole moment perpendicular to the axis along which the atoms predominantly move. In the *E′′(x) mode, for instance, S(Se) atoms primarily move in the y direction and the curvature introduces a component of the displacement patterns in the x direction while in the *E′′(y) mode, the atoms primarily move in the x direction and the curvature is responsible for the component in the y direction. Therefore, the dipole moment component on the primary axis of movement of the atoms is canceled, and the resultant component on its perpendicular axis is due to curvature induced by the ripple. The emergence of a resultant dipole moment of the cell is apparent on the front view of the *E′′(x) mode, with all arrows slightly pointing to the right, and on the side view of the *E′′(y), with all arrows slightly pointing to the right and other modes with similar displacement patterns but with atoms on opposite sides of the crest move in the direction were found to have no IR activity. Our results for the eigenmodes suggest that this mechanism described for the *E′′(x) and *E′′(y) modes is consistent with the appearance of the peaks close to the E′′ frequency at Fig. 17 and provides a feasible explanation for them. Fig. 24 illustrates mode associated with the red bars in Fig. 20(b-iii) for MoS₂ and Fig. 21(b-ii) for MoSe₂, indicating their IR activity near . Displacement vectors are color-coded according to maximum displacement along the z-axis over half oscillation period, with red indicating greater displacement. This mode exhibit an oscillation pattern similar to , with the Mo-layer remaining stationary while the top and bottom S(Se)-layers moving symmetrically. However, upon examining the color of the displacement vectors, it is evident that the bottom layer exhibits larger displacements on average compared to the top layer, resulting in a greater distance from the Mo-layer. This disparity in displacement between the layers can induce a slight modification in the total dipole moment of the cell, thereby generating the observed out-of-plane IR activity depicted in Fig. 20(b-iii) for MoS₂ and Fig. 21(b-ii) for MoSe₂. Our results for the eigenmodes suggest that this mechanism described for the mode is consistent with the appearance of the peaks close to the frequency at Fig. 17 and provides a feasible explanation for it. It is noteworthy that several modes exhibiting planar IR activity were identified near the E′ mode in Fig. 20(b-ii) and near the mode in Fig. 21(b-ii). However, since the frequency is close to the E′ in both compounds, we believe that these peaks were suppressed by the broadening of the E′ mode main blue peak in our simulations as shown in Fig. 17.

Fig. 22 Sketch of the oscillation patterns of the modes *E′(z) and in the rippled cell. The directions of displacements of the Mo and S(Se) atoms are indicated by blue and red vectors respectively.

Fig. 23 Sketch of the oscillation patterns of the modes *E′′(x) and *E′′(y) in the rippled cell. The directions of displacements of the S(Se) atoms are indicated by red vectors.

Fig. 24 Sketch of the oscillation patterns of the mode in the rippled cell. The directions of displacements along the z-axis of the S(Se) atoms are indicated by vectors and its colors shows the maximum displacement of each during the first half period of oscillation, with the scale increasing from blue to red.

4 Conclusions

In this work, using machine learning, we developed two MLFFs for MoS₂ and MoSe₂ that captures their main dynamical properties. We tested its validity by calculating its phonon dispersion and making an accurate correspondence between the structures included in the reference data set and the modes desired to model on MD simulations. With that procedure, we were able to reproduce known experimental dielectric spectra at THz frequencies. To test the performance of the force fields generated in this work, we focused on reproducing results for the peak-frequencies and their redshift due to temperature and strain, comparing it with well known Stillinger–Weber parametrizations of the respective materials. In general, the phonon dispersion obtained using our MLFF were more accurate than Stillinger–Weber as well the MD results for the IR-active phonon frequencies. For the temperature coefficient, MLFF and SW presented similar capability of recreating experimental results for MoS₂ while for MoSe₂ the MLFF has a better performance. For the MoS₂ strain coefficient, SW predicted that the mode shows a blueshift, which contradicts other ab initio results previously reported. Simulations with MLFF predicted a redshift for the strain coefficient of both infrared active modes, which is in accordance to previous experimental and ab initio results. Furthermore, our simulations reveal the appearance of new low-intensity peaks in the THz dielectric spectra, especially those in proximity to the IR-inactive frequencies E′′ and . Molecular dynamics (MD) simulations investigating the height Δz of strained lattices establish a correlation between these emerging peaks and the propagation of thermal ripples, as evidenced by their significant reduction on lattices with smaller Δz. To explore potential explanations for the new observed peaks, we have calculated the eigenfrequencies and eigenmodes of a sine-rippled supercell with wave vector parallel to the zigzag direction using our MLFFs for each compound. The eigenmodes results for the frequencies next to E′ and with IR activities with inverted polarizations strengths our initial hypothesis suggesting that permitting the layer to undulate frees the modes from oscillating exclusively along their original axis, enabling a component along the perpendicular axis. Moreover, the eigenmodes results for frequencies next to E′′ suggests that, according the MLFFs modeled in this work, the ripples induces a curvature on the displacement vector field of the oscillation pattern and this curvature is responsible for the appearance of the peaks related to E′′ on the spectra. The reasons for the appearance of such curvature on a rippled layer have not yet been discovered. Finally, the eigenmodes results for frequencies next to suggests that the resultant dipole moment at z is originated by the disparity of the distance of the top and the bottom layer with respect to the Mo-layer. These results may have far-reaching implications on the understanding of thermal ripples on 2D materials and, after testing the reliability and functionality of the MLFFs, they can be applied to study similar problems related to the phonon properties of MoS₂ and MoSe₂. In future work, we plan to extend both MLFFs to study the role of defects, such as vacancies and dopants, on the layer dynamics.⁸⁵

Data availability

All data supporting the conclusions of this study, including the MLFFs presented in this work, the database used for its development and animations for the discussed modes are available on the following Github repository: https://github.com/gbsouza1997/TMD-MLFF-repo.git.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

This work received financial support from Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brazil (CNPq) under grants CNPq 402091/2012-4, CNPq 163465/2021-5, Fundação de Amparo à Pesquisa do Estado de Minas Gerais - Brazil (FAPEMIG) under grants FAPEMIG RED-00458-16 and FAPEMIG APQ-01414-22, Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) under grant CAPES 88887.840135/2023-00. G. B. G was also partially supported by the Simons Foundation and the Vice Provost for Graduate Education at The University of Georgia. S. B. H. and Y. A. acknowledge support from the Air Force Office of Scientific Research (AFOSR) grant numbers FA9550-19-0252 and FA9550-23-1-0375. This study was supported in part by resources and technical expertise from the Georgia Advanced Computing Resource Center (GACRC), a partnership between the University of Georgias Office of the Vice President for Research and Office of the Vice President for Information Technology. G. B. G and V. B also thank Dr Atanu Samanta from University of Pennsylvania (UPENN) for fruitful discussions.

References

1. A. Molina-Sánchez, K. Hummer and L. Wirtz, Surf. Sci. Rep., 2015, 70, 554–586.
2. M. Chhowalla, Z. Liu and H. Zhang, Chem. Soc. Rev., 2015, 44, 2584–2586.
3. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Phys. Rev. Lett., 2010, 105, 136805.
4. A. Ramasubramaniam, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 115409.
5. V. K. Sangwan, H.-S. Lee, H. Bergeron, I. Balla, M. E. Beck, K.-S. Chen and M. C. Hersam, Nature, 2018, 554, 500–504.
6. X. Wang and A. Tabarraei, Appl. Phys. Lett., 2016, 108, 191905.
7. M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. Da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain and Z.-X. Shen, et al. , Nat. Mater., 2014, 13, 1091–1095.
8. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017.
9. S. B. Hancock, D. P. Landau, N. A. Aghamiri and Y. Abate, Phys. Rev. Mater., 2022, 6, 076001.
10. D. C. Rapaport, The Art of Molecular Dynamics Simulation, Cambridge University Press, 2nd edn, 2004.
11. J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 6991–7000.
12. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
13. D. Braun, S. Boresch and O. Steinhauser, J. Chem. Phys., 2014, 140, 064107.
14. P. Höchtl, S. Boresch, W. Bitomsky and O. Steinhauser, J. Chem. Phys., 1998, 109, 4927–4937.
15. P. Vashishta, R. K. Kalia, J. P. Rino and I. Ebbsjö, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 12197–12209.
16. M. S. Daw and M. I. Baskes, Phys. Rev. B: Condens. Matter Mater. Phys., 1984, 29, 6443.
17. J.-W. Jiang, H. S. Park and T. Rabczuk, J. Appl. Phys., 2013, 114, 064307.
18. J.-W. Jiang, Nanotechnology, 2015, 26, 315706.
19. A. Ostadhossein, A. Rahnamoun, Y. Wang, P. Zhao, S. Zhang, V. H. Crespi and A. C. T. van Duin, J. Phys. Chem. Lett., 2017, 8, 631–640.
20. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142–10186.
21. M. Fox, Optical Properties of Solids, OUP Oxford, 2010.
22. G. Srivastava, The Physics of Phonons, Routledge, 1990.
23. N. D. Mermin, Phys. Rev., 1968, 176, 250–254.
24. D. Nelson, T. Piran and S. Weinberg, Statistical Mechanics Of Membranes And Surfaces, World Scientific Publishing Company, 2nd edn, 2004.
25. A. Fasolino, J. Los and M. I. Katsnelson, Nat. Mater., 2007, 6, 858–861.
26. J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth and S. Roth, Nature, 2007, 446, 60–63.
27. J. Brivio, D. T. L. Alexander and A. Kis, Nano Lett., 2011, 11, 5148–5153.
28. P. Anees, M. C. Valsakumar and B. K. Panigrahi, Phys. Chem. Chem. Phys., 2017, 19, 10518–10526.
29. S. Li, Y.-P. Wang, S. Ning, K. Xu, S. T. Pantelides, W. Zhou and J. Lin, Nano Lett., 2023, 23, 1298–1305.
30. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
31. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. De-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter., 2017, 29, 465901.
32. P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru and S. Baroni, J. Chem. Phys., 2020, 152, 154105.
33. D. R. Hamann, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 085117.
34. M. Schlipf and F. Gygi, Comput. Phys. Commun., 2015, 196, 36–44.
35. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
36. R. G. Dickinson and L. Pauling, J. Am. Chem. Soc., 1923, 45, 1466–1471.
37. P. B. James and M. Lavik, Acta Crystallogr., 1963, 16, 1183.
38. N. Artrith and A. Urban, Comp. Mater. Sci., 2016, 114, 135–150.
39. A. M. Cooper, J. Kästner, A. Urban and N. Artrith, npj Comput. Mater., 2020, 6, 54.
40. N. Artrith, A. Urban and G. Ceder, Phys. Rev. B, 2017, 96, 014112.
41. C. G. Broyden, J. Inst. Maths. Appl., 1970, 6, 76–90.
42. C. Zhu, R. H. Byrd, P. Lu and J. Nocedal, ACM Trans. Math. Softw., 1997, 23, 550–560.
43. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comp. Phys. Commun., 2022, 271, 108171.
44. M. S. Chen, T. Morawietz, H. Mori, T. E. Markland and N. Artrith, J. Chem. Phys., 2021, 155, 074801.
45. J. M. Caillol, D. Levesque and J. J. Weis, J. Chem. Phys., 1986, 85, 6645–6657.
46. S. Nosé, Mol. Phys., 1984, 52, 255–268.
47. I. Ponomareva, L. Bellaiche, T. Ostapchuk, J. Hlinka and J. Petzelt, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 012102.
48. H. Tornatzky, R. Gillen, H. Uchiyama and J. Maultzsch, Phys. Rev. B, 2019, 99, 144309.
49. T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys.: Condens. Matter, 2014, 26, 225402.
50. C. Lee, H. Yan, L. E. Brus, T. F. Heinz, J. Hone and S. Ryu, ACS Nano, 2010, 4, 2695–2700.
51. T. J. Wieting and J. L. Verble, Phys. Rev. B, 1971, 3, 4286–4292.
52. A. Garg, H. Sehgal and O. Agnihotri, Sol. State Commun., 1973, 12, 1261–1263.
53. Q.-C. Sun, X. S. Xu, L. I. Vergara, R. Rosentsveig and J. L. Musfeldt, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 205405.
54. T. J. Wieting and J. L. Verble, Phys. Rev. B, 1971, 3, 4286–4292.
55. S. Sahoo, A. P. S. Gaur, M. Ahmadi, M. J.-F. Guinel and R. S. Katiyar, J. Phys. Chem. C, 2013, 117, 9042–9047.
56. Q. Wang, L. Han, L. Wu, T. Zhang, S. Li and P. Lu, Nanoscale Res. Lett., 2019, 14, 1–9.
57. J. Wilson and A. Yoffe, Adv. Phys., 1969, 18, 193–335.
58. S. V. Bhatt, M. P. Deshpande, V. Sathe, R. Rao and S. H. Chaki, J. Raman Spectrosc., 2014, 45, 971–979.
59. S. Najmaei, P. M. Ajayan and J. Lou, Nanoscale, 2013, 5, 9758–9763.
60. S. Paul, S. Karak, A. Mathew, A. Ram and S. Saha, Phys. Rev. B, 2021, 104, 075418.
61. A. S. Pawbake, M. S. Pawar, S. R. Jadkar and D. J. Late, Nanoscale, 2016, 8, 3008–3018.
62. M. Yang, X. Cheng, Y. Li, Y. Ren, M. Liu and Z. Qi, Appl. Phys. Lett., 2017, 110, 093108.
63. N. A. Lanzillo, A. Glen Birdwell, M. Amani, F. J. Crowne, P. B. Shah, S. Najmaei, Z. Liu, P. M. Ajayan, J. Lou, M. Dubey, S. K. Nayak and T. P. O’Regan, Appl. Phys. Lett., 2013, 103, 093102.
64. Z. Li, Y. Lv, L. Ren, J. Li, L. Kong, Y. Zeng, Q. Tao, R. Wu, H. Ma and B. Zhao, et al. , Nat. Commun., 2020, 11, 1151.
65. H. J. Conley, B. Wang, J. I. Ziegler, R. F. J. Haglund, S. T. Pantelides and K. I. Bolotin, Nano Lett., 2013, 13, 3626–3630.
66. C. Rice, R. J. Young, R. Zan, U. Bangert, D. Wolverson, T. Georgiou, R. Jalil and K. S. Novoselov, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 081307.
67. N. Basu, R. Kumar, D. Manikandan, M. Ghosh Dastidar, P. Hedge, P. K. Nayak and V. P. Bhallamudi, RSC Adv., 2023, 13, 16241–16247.
68. Y. Cai, J. Lan, G. Zhang and Y.-W. Zhang, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 035438.
69. A. Kandemir, H. Yapicioglu, A. Kinaci, T. Çağin and C. Sevik, Nanotechnology, 2016, 27, 055703.
70. T. Sekine, M. Izumi, T. Nakashizu, K. Uchinokura and E. Matsuura, J. Phys. Soc. Jpn., 1980, 49, 1069–1077.
71. M. Jin, W. Zheng, Y. Ding, Y. Zhu, W. Wang and F. Huang, J. Phys. Chem. Lett., 2020, 11, 4311–4316.
72. S. Tongay, J. Zhou, C. Ataca, K. Lo, T. S. Matthews, J. Li, J. C. Grossman and J. Wu, Nano Lett., 2012, 12, 5576–5580.
73. S.-i Uchida and S. Tanaka, J. Phys. Soc. Jpn., 1978, 45, 153–161.
74. G. Lucovsky, R. M. White, J. A. Benda and J. F. Revelli, Phys. Rev. B, 1973, 7, 3859–3870.
75. T. J. Wieting and J. L. Verble, in Infrared and Raman Investigations of Long-Wavelength Phonons in Layered Materials, ed. T. J. Wieting and M. Schlüter, Springer Netherlands, Dordrecht, 1979, p. 385.
76. A. Garg, H. Sehgal and O. Agnihotri, Sol. State Commun., 1973, 12, 1261–1263.
77. M. Öper, Y. Shehu and N. K. Perkgöz, Semicond. Sci. Technol., 2020, 35, 115020.
78. D. Kumar, V. Kumar, R. Kumar, M. Kumar and P. Kumar, Phys. Rev. B, 2022, 105, 085419.
79. D. J. Late, S. N. Shirodkar, U. V. Waghmare, V. P. Dravid and C. N. R. Rao, ChemPhysChem, 2014, 15, 1592–1598.
80. S. Horzum, H. Sahin, S. Cahangirov, P. Cudazzo, A. Rubio, T. Serin and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 125415.
81. X. Cheng, L. Jiang, Y. Li, H. Zhang, C. Hu, S. Xie, M. Liu and Z. Qi, Appl. Surf. Sci., 2020, 521, 146398.
82. W. Bührer and P. Brüesch, Phonons: Theory and Experiments II: Experiments and Interpretation of Experimental Results, Springer Berlin Heidelberg, 2012.
83. P. Giannozzi and S. Baroni, J. Chem. Phys., 1994, 100, 8537–8539.
84. X. Gonze and C. Lee, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 10355.
85. H.-P. Komsa, S. Kurasch, O. Lehtinen, U. Kaiser and A. V. Krasheninnikov, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 035301.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5cp00163c

‡ The justification for the large variation of the lattice constant is that our database needs to contain information for when atoms are very close (r → 0) and very far (r → ∞) from each other to make our MLFFs reproduce those conditions.

§ The index κ denotes the atom within the unit cell. The vector U_j(κ) can be expressed in terms of the Cartesian unit vectors x̂, ŷ, and ẑ as U_j(κ) = U_j(κ, x)x̂ + U_j(κ, y)ŷ + U_j(κ, z)ẑ. The components referenced in eqn (3) are defined accordingly in this footnote.

---