Revisiting the origin of non-volatile resistive switching in MoS2 atomristor

Introduction

Atomristors represent the ultimate scaling of memristors in the vertical stacking configuration, i.e., Metal-Monolayer 2D Material-Metal¹. Recent experiments demonstrated NVRS in these monolayer 2D materials using various combination of Au/Ag/G^2,3,4 electrodes, however, the origin of NVRS remains debatable. On the experimental front, the nature of NVRS is extrinsic, caused by the electrically induced dynamic adsorption of Au atoms into sulfur vacancies in MoS₂ in Au-MoS₂-Au and Au-MoS₂-G atomristors^3,4,5. On the other hand, atomistic computations indicate the existence of both intrinsic^6,7 and extrinsic^7,8,9 NVRS in MoS₂ atomristors. The intrinsic NVRS is caused by the popping of a sulfur atom at a monosulfur vacancy (parent state) into the molybdenum plane (popped state) in MoS₂. For extrinsic NVRS, we showed in our previous work⁸ how single Au atom adsorption at a S vacancy modulates the local environment around the adsorption site thereby reducing exit barrier significantly at the vdW interface of Au-MoS₂-Au atomristor. However, apart from Au atom adsorption at S vacancies, the possibility of intrinsic non-volatile resistive states emerging from MoS₂ have not been investigated properly in light of an ab-initio framework such as Density Functional Theory (DFT). Since in an actual atomristor it is not possible to perform in-situ experiments for these 3-atom thick 2D monolayer TMD materials to observe atomic events taking place at Metal-2D material interface, accurate atomistic computations are required to investigate the role of defects in NVRS in 2D atomristors.

DFT as an indispensable tool has been used for studying static properties of various atomic systems hosting upto 1000 atoms. At the same time utilizing DFT for Ab-initio Molecular Dynamics (AIMD) is extremely time consuming and resource demanding for a system with even 100–200 atoms. To overcome this problem, researchers developed various simplified forcefields such as Tersoff¹⁰, Embedded Atom Method (EAM)¹¹, Lennard Jones¹² and Reactive (ReaxFF)¹³ etc. These forcefields are used to study the dynamics of a given system across nanosecond timescales with moderate computational resources. Among all these forcefields, only ReaxFF provides description of dynamic bond making and breaking and is often used for investigating accurate dynamics of a given system. The downside of such forcefields is that they are fitted/trained on a narrow range of dataset that can help describe certain but not all properties of a system. One such well known reactive forcefield for MoS₂ was proposed by Ostadhossein et al.¹⁴, which can describe covalent interactions between molybdenum (Mo), sulfur (S) and hydrogen (H). It has been utilized to study many MoS₂ based systems^{15,16,17,18,19} but unfortunately it gives an incorrect description of popped states in MoS₂. Apart from the above-mentioned conventional forcefields, recently Machine Learning (ML) and Deep Learning (DL) have been used to develop complex forcefields with multiple atomic compositions. These machine/deep learned forcefields also allow accessibility across nanosecond timescales and at the same time have an accuracy near that of ab-initio methods. Chen et al., developed one such deep learned universal interatomic forcefield – M3GNet²⁰. It is based on graph neural networks and includes three body interactions. M3GNet was trained on 10 years of massive database generated by Materials project²¹ from VASP while utilizing plane-wave basis set and can be used for all the elements in the periodic table.

In this work, for the first time, rigorous DFT and M3GNet based calculations are carried out to investigate the previously reported intrinsic resistive state in MoS₂ atomristors caused by the popping of S atom at monosulfur vacancy into the molybdenum plane. It is found that such an intrinsic state is neither stable nor non-volatile but is a false manifestation of underlying ReaxFF¹⁴ due to inaccuracies in its parameterization. It is also found that Au atoms adsorbed at sulfur vacancies represent stable and non-volatile resistive state as observed in the STM experiment⁵. Furthermore, local heating rates and local sulfur (S) vacancy diffusion are identified as leading causes of cycle-to-cycle V_RESET and V_SET variability observed in experiments, respectively.

Results and discussion

Revisiting intrinsic NVRS in MoS₂

Recently, Mitra et al.⁶, reported an intrinsic NVRS in G-MoS₂-G atomristor using reactive forcefield (ReaxFF) molecular dynamics. They reported that during the set process a vertical electric field pops S atom at a monosulfur vacancy (parent state) (Fig. 1a) into the molybdenum plane (popped state) (Fig. 1b) which then stays stable in that position at 300 K and 1 atm, thereby representing a non-volatile Low Resistive State (LRS) in G-MoS₂-G atomristors. Furthermore, they show that during reset process a high local temperature (Joule’s heating) alongwith an electric field is required to bring the popped state back to its parent state which represents a non-volatile High Resistive State (HRS).

**Fig. 1: Schematic of 1D MoS₂ chain showing.**

The ReaxFF used in their study describes chemical interaction between Mo, S and H atoms only, whereas the interaction with graphene electrode was introduced through an LJ potential. On the other hand, M3GNet is a universal interatomic potential which can describe interactions for all the atoms in the periodic table. It is important to note that both M3GNet and ReaxFF provide description of a system at an atomic level and do not contain any electronic information. However, DFT which is more fundamental theory provides description at an electronic level and has been used as a tool for generating training dataset to develop ReaxFF and M3GNet forcefields. Thus, DFT stands at the highest level of accuracy and any observations made from ReaxFF and M3GNet need to be consistent with DFT in-order to be valid and any inconsistencies arising between a forcefield and DFT is bound to be a parametric error in the design of forcefield.

To demystify the origin of intrinsic non-volatile resistive states in detail, an orthogonal supercell of MoS₂ with dimensions of 2.18 × 1.89 nm² was taken. Then, one monosulfur vacancy was created resulting in a defect density of 2.4 × 10¹³ cm⁻² as observed in experiments^5,22 (Fig. 2a). The molecular dynamics were carried out independently using two different forcefields, i.e., reactive forcefield- ReaxFF¹⁴ and deep learned forcefield- M3GNet²⁰. While in experiments the application of an electric field at a given value is of the order of micro-milliseconds, simulating such longer time scales is not possible in simulation setups due to increasing computational and/or time complexity. Therefore, electric field of a given value was applied for a short time-period which was enough to capture the findings from previous study^6,7 whilst demystifying the origin of intrinsic NVRS.

The electric field was increased in steps of 0.5 V/Å after every 0.1 ns (Fig. 2b). It was observed that the dynamics carried out using ReaxFF resulted in popping of S atom at monosulfur vacancy into the molybdenum plane (popped state) at a critical electric field of 2.5 V/Å in agreement with previous study⁶ (Fig. 1b, Supplementary Video 1- The red highlighted atom in supplementary videos denotes sulfur atom at monosulfur vacancy). The dynamics carried out using M3GNet resulted in S atom at a monosulfur vacancy to move from one S plane to another S plane (swapped state) at a critical electric field of 2.5 V/Å (Fig. 1c, Supplementary Video 2). Furthermore, after removing the electric field the popped S atom obtained from ReaxFF based molecular dynamics remained in the molybdenum plane indicating a possible intrinsic non-volatile resistive state. The M3GNet, however, did not reveal any such intrinsic non-volatile resistive state. It is important to note that in the previous study the popped state was observed with inert graphene electrodes⁶, nevertheless same popped state was observed here without any electrodes. This suggests that there is no role of graphene electrodes in observing the popped states and these states self-emerge from MoS₂ ReaxFF under out-of-plane electric field due to its limitation in properly describing Potential Energy Surface (PES) around popped states as discussed further in sections below.

Supercell optimization & total energy analysis

To settle the confusion between two contrasting results obtained from two different forcefields, further investigation about the popped state was carried out rigorously in light of DFT. MoS₂ supercells of varying sizes were taken, i.e., 7 × 7, 8 × 8, 9 × 9, 10 × 10, each hosting one monosulfur vacancy with S atom initialized to popped state. All these supercells were relaxed using three different methods- DFT, M3GNet and ReaxFF. For DFT, Linear Combination of Atomic Orbitals (LCAO) basis set was used since using planewave DFT becomes computationally expensive for full relaxation of large supercells. However, to test with planewave basis set, we used indirect approach by utilizing M3GNet. Finally, supercells were also relaxed using ReaxFF which was used in the earlier works^6,7. It was observed that all the relaxed supercells from ReaxFF had the S atom in popped state (Fig. 1b) (Supplementary Fig. 1). However, for M3GNet and DFT, the relaxed supercells had the S atom returned to parent state (Fig. 1a). Furthermore, we took eight different MoS₂ supercells, four of which had one S atom in popped state and other four had S atom in parent state. All these supercells were first relaxed using ReaxFF and then DFT was used to calculate total energies of these supercells (Table 1). Clearly, MoS₂ supercells with popped states have a higher energy by about ≈ 8.1 eV than parent states for all the supercell sizes. This huge energy difference indicates that popped states (if they exist) are less stable than the parent states and require higher external energies to switch. This translates to higher V_SET (voltage required for HRS to LRS transition) voltages in G-MoS₂-G atomristors. Besides providing such high energies through external electric fields, it is important that there should be an energy valley in the PES of MoS₂ for the popped state which will hold S atom in the molybdenum plane to represent a non-volatile resistive state.

Table 1 DFT calculated energies of parent and popped states for various supercell sizes

Full size table

Minimum Energy Path (MEP) from parent to popped state

To further study the popped state, Climbing Image Nudge Elastic Band (CINEB) calculations were carried out on a 7 × 7 MoS₂ supercell to investigate Minimum Energy Path (MEP) connecting the parent and popped state in the Potential Energy Surface (PES) of MoS₂ using all the three methods mentioned above (Fig. 3a). All computational settings were kept same as mentioned in the methods section except for DFT, the atomic forces and stresses were relaxed below 0.04 eV/Å and 0.1 GPa, respectively. It can be observed that both DFT and M3GNet are in a qualitative agreement with each-other, predicting the MEP connecting the two states as concave upwards with energy migration barriers of ≈ 2.6 and 1.2 eV, respectively. The initial/final state in MEP represents MoS₂ with parent/popped state. Clearly, there is no energy valley at top (popped state) to hold S atom in the molybdenum plane. Thus, a small perturbation due to thermal vibrations will steer it down to parent state as discussed further in molecular dynamics section. On the other hand, only ReaxFF predicts the MEP to have an energy valley for popped state with a migration barrier of ≈ 1.5 eV, similar to previous reporting⁶. Furthermore, it is observed that the popped state remains unstable under both low (10 × 10-MoS₂) and high (7×7-MoS₂) defect densities (Fig. 3b). Thus, even if an external electric field does provide sufficient energy required for transition to popped state, MoS₂ does not have an energy valley in its PES for popped state to represent a non-volatile resistive state. Therefore, the reason why M3GNet and DFT did not relax the supercells to popped state is due to the absence of energy valley therein and tendency of optimization algorithms to steer towards lower energy minimum, in this case the parent state. The underestimation of CINEB migration energies obtained from M3GNet against DFT is due to the absence of defect configurations in the training data of M3GNet causing PES softening²³. Nonetheless, M3GNet provides correct qualitative description of energy landscape around defects. Also, the radial and angular distribution functions obtained from M3GNet show good agreement with DFT in describing atomic environment of MoS₂ (Supplementary Fig. 2). Moreover, M3GNet generates molecular dynamics in agreement with DFT and experiment which suffices to test the nature of possible resistive states as explained in the following section. As for ReaxFF, it is in complete contrast with DFT and M3GNet, highlighting its limitations in predicting structural properties under large electric field induced deformations²⁴.

**Fig. 3: Minimum Energy Path (MEP) as predicted by CINEB calculations.**

Assessing the stability and non-volatility of popped states using molecular dynamics

To this point, popped state analysis carried out show that ReaxFF is in contrast with more accurate models (DFT & M3GNet). However, it is possible that under certain conditions and different interfaces, it might become stable and remain non-volatile. Therefore, molecular dynamics was carried out to assess the stability and non-volatility of popped states under ambient conditions in both MoS₂ and MoS₂ atomristors.

Molecular dynamics for MoS₂

A 7 × 7 supercell of MoS₂ was taken with S atom initialized to popped state and relaxed using ReaxFF. The molecular dynamics were then carried out using NPT ensemble since it describes most real systems, however, contrary to the claim⁶, it was found that the popped state relaxes back to parent state within few femtoseconds when using DFT (Supplementary Video 3) or M3GNet (Supplementary Video 4). The average S-Mo bond length in popped/parent states is 1.91/2.40 Å (Fig. 4a, b). As the system evolves, the S atom which was initialized to popped state relaxes back to parent state as indicated by a shoot in red/blue curves (Fig. 4c) from 1.91 Å to 2.82/2.62 Å after which it settles at an average bond length of 2.40 Å from neighboring Mo atoms. To further avoid any doubts of incorrect equilibration using NPT ensemble, popped states were tested using NVT Berendsen ensemble which is known to give best equilibration. To no surprise, it was again found the popped state relaxes back to the parent state using DFT and M3GNet (Fig. 4d) (DFT – Supplementary Video 5 & M3GNet – Supplementary Video 6). On the other hand, when molecular dynamics were carried out using ReaxFF, the S atom remained stable in the popped state in the molybdenum plane for both NPT (Supplementary Video 7) and NVT (Supplementary Video 8) ensembles (Fig. 4c, d), with an average S-Mo bond length of 1.91 Å. Given there is a consistent agreement between DFT and M3GNet so far, next we proceed to utilize M3GNet for investigating MoS₂ atomristors since carrying out ab-initio molecular dynamics using DFT is not practically feasible for an atomristor device hosting more than 700 atoms.

Molecular dynamics for MoS₂ atomristors

Three atomristors structures were created, i.e., Graphene-MoS₂-Graphene (GMG-hypothesised⁶), Gold-MoS₂-Graphene (AuMG-experimental^2,3) and Gold-MoS₂-Gold (AuMAu-experimental^2,4) with S atom initialized to popped state (Fig. 5a, c, e). The molecular dynamics were carried out using both NPT (Supplementary Videos 9–11) and NVT (Supplementary Videos 12–14) ensembles while utilizing M3GNet as forcefield. However, M3GNet lacks account for vdW interactions between materials. To overcome this issue, M3GNet was combined with Grimme’s D3Z vdW correction²⁵ using PBE exchange-correlation functional. It was again found in all the atomristors that the popped state always relaxes back to parent state within few femtoseconds in both NPT (Fig. 5b, d, f) and NVT ensembles (Supplementary Fig. 3). Fig. 5g provides time evolution of popped S atom which settles at an average Mo-S bond length of 2.4 Å in parent state. The MD simulations were also carried out with upto 9 Au layers while keeping outer few layers fixed and similar results were observed.

**Fig. 5: M3GNet molecular dynamic (NPT) snapshots taken for various atomristors.**

Revisiting intrinsic NVRS in multilayer MoS₂

More recently, intrinsic non-volatile resistive states were also theorized to exist in multilayer MoS₂⁷. Therefore, non-volatility of popped states was tested in multilayer MoS₂ using ReaxFF and M3GNet forcefields. A tri-layer orthogonal supercell of MoS₂ having dimensions of 1.93 × 2.22 nm² with multiple S atoms initialized in popped states (Fig. 6a–c) was taken. A vacuum padding of 30 Å was added in the out-of-plane direction to avoid spurious interactions. It was found that only ReaxFF based molecular dynamics holds S atoms in popped states (Supplementary Video 15) while M3GNet based molecular dynamics relaxes them back to parent state within few femtoseconds (Supplementary Video 16). This implies that irrespective of the number of MoS₂ layers interlayer vdW interactions also do not help stabilize these popped states.

**Fig. 6: Revisiting NVRS in multilayer MoS₂.**

Do adsorbed Au atoms at sulfur vacancies represent non-volatile resistive states?

Experimental observations show gold atom adsorption at a sulfur vacancy causes NVRS in gold contacted MoS₂ atomristors⁵. Here molecular dynamics were carried out using M3GNet combined with Grimme’s D3Z vdW correction to analyze the non-volatility of Au atom adsorbed at monosulfur and disulfur vacancy. Before testing the non-volatile nature of Au adsorbed atoms, DFT-CINEB calculations were carried out to calculate Au atom migration energies from Au electrode to adsorption at a monosulfur vacancy (V_S). It is important to note that in an experimental atomristor, the Au surface at Au-MoS₂ interface hosts Au ad-atoms⁹, vacancies (leading to atomic roughness) and random interface orientation separated by a vdW gap from the MoS₂. Moreover, in DFT CINEB calculations, the migration path of an Au atom from Au metal to sulfur vacancy is calculated at 0 K while in experiments the thermal vibrations and electric field profile at interface will introduce additional complexities in the migration path of an Au atom before it adsorbs at a sulfur vacancy resulting in modifications of energy barriers. But due to the computational complexities, the DFT-CINEB calculations in this work were carried out on a simplified interface with flat Au (111) surface next to MoS₂. The atomic relaxations were carried out until forces fell below 0.04 eV/Å. The CINEB calculations were carried out for two Au atoms (Au₁/Au₂) at a distance of 3.10 Å and 3.76 Å from a neighboring V_S, respectively. The migration barriers for Au₁V_S/ Au₂V_S are 0.22/0.49 eV (Fig. 7a–d) in accordance with previous reporting (0.378 eV)²⁶. It is observed that a slight difference in Au-V_S distances can affect migration barriers significantly. Now in experiments the observed resistance switching is not due to the adsorption of one Au atom at a sulfur vacancy but due to multiple Au atom adsorptions. All these Au atoms will have different distances from neighboring sulfur vacancies, thus, different energy barriers. The observed V_SET in experiments is then a reflection of average over all those energy barriers under realistic conditions. Therefore, V_SET will depend largely on the atomic description of Au surface above MoS₂ and migration path Au atoms take from Au metal to sulfur vacancies which will give out a different average over these energy barriers from device-to-device or cycle-to-cycle explaining the observed V_SET variability^2,3,4. Moreover, the energy barriers in our calculation might be low due to ideal conditions nevertheless they elucidate the impact of Au atom-V_S distances on energy barriers. A similar argument follows for disulfur vacancies though they have low defect density. Next, two atomristors of the type AuMG were considered having Au atom adsorption at monosulfur vacancy and disulfur vacancy, respectively. The molecular dynamics were carried out using NPT ensemble for a period of 10 ps (Supplementary Videos 17 and 18). The MD snapshots (Fig. 8a, b) taken at 10 ps show Au atoms remaining adsorbed at mono and disulfur vacancies, indicating a stable behavior. Furthermore, (Fig. 8c) shows time evolution of average Au-Mo (nearest 3 Mo atoms) bond length settling at 2.87 Å. Thus, Au atom adsorbed states are truly non-volatile in nature in agreement with experiments^2,3,5 and our previous AIMD study carried on MoS₂ only⁸.

**Fig. 7: DFT-CINEB calculation for Au-MoS₂ interface in atomristor (View from first Au interfacial layer).**

**Fig. 8: Snapshots taken from NPT molecular dynamics of AuMG atomristor at 10 ps with Au atom adsorbed at (as indicated by red dashed circle).**

Desorption and cycle to cycle variability

Evidences from direct experiment⁵ and atomistic computation⁷ show Au atom adsorption at sulfur vacancies under applied field result in LRS. As for desorption process, due to the high electrical conductivity across local filament (Au adsorbed site) local heating (joule’s heating) plays a critical role. Nonetheless, how this leads to cycle-to-cycle variability has remained an open question. To understand this, M3GNet based molecular dynamics were carried out under various local heating rates and electric field. Due to M3GNet being computationally resource expensive for long time simulations (Supplementary Table 1), we took only an MoS₂ supercell of size 2.52 × 2.73 nm² with an Au atom adsorbed at monosulfur vacancy. The local heating was introduced by increasing temperature at different rates for all the atoms in 4 Å radius around the Au adsorbed site. At the same time, the electric field was also gradually increased at a rate of 0.05 V/Å per 0.1 ns in the direction as shown in Fig. 9b. It was observed that when only electric field is present the Au atom gets desorbed at 0.5 V/Å, however, when the filament was heated locally at different rates the desorption occurred earlier. The electric field values for Au atom desorption at a local heating rate of 1 Kelvin/ps and 3 Kelvin/ps were 0.25 V/Å and 0.15 V/Å, respectively (Fig. 9a). This reveals how different local heating rates result in cycle-to-cycle V_RESET (voltage required for LRS to HRS transition) variability. Furthermore, post Au-atom desorption from a sulfur vacancy the generated local heating will take some time to dissipate through the atomristor. This depends on the thermal conductivity of in- plane MoS₂ and across Metal-MoS₂ vdW interfaces in an atomristor (Fig. 9b). While an accurate estimation of these factors is a matter of further research, nevertheless, it was observed that if the electric field and local heating are allowed to sustain for few picoseconds after desorption, then this kicks in S vacancy diffusion at the adsorption site. The sustaining of electric field can be attributed to the finite fall time of applied electric pulses in experiments.

**Fig. 9: Desorption process and cycle-to-cycle variability.**

For this, an electric field of 0.01 V/Å and local heating at the rate of 3 Kelvin/ps was applied. It was observed that after 42 ps of Au atom desorption the S vacancy diffused in the in-plane direction of MoS₂ (Fig. 10a, b) (Supplementary Video 19). While the desorbed Au atom can be seen moving randomly over MoS₂ but in an atomristor it will go back to Au electrode and sit at any random position. In the next SET cycle, both these phenomena (random position of Au atom and V_S diffusion) would lead to change of Au-V_S distances and energy migration barriers, eventually causing cycle-to-cycle V_SET variability.

**Fig. 10: Sulfur vacancy diffusion after desorption.**

Limitations of forcefields

It is noteworthy that ReaxFF was trained on pristine and defect configurations but fails to predict both qualitative and quantitative PES description around monosulfur vacancies. On the other hand, M3GNet was not trained on any defect configurations yet it did capture correct qualitative behavior in accordance with DFT. As such M3GNet not only unveiled reality of popped states but also revealed factors responsible for cycle-to-cycle variability. Nevertheless, this does not deem ReaxFF/M3GNet incorrect/correct for all applications, but this work highlights the importance of validating the results obtained from forcefields with DFT/experiments in-order to build a better picture of reality.

To conclude, this work investigated the nature of popped sulfur atom in the molybdenum plane of MoS₂ which was theorized to give rise to intrinsic non-volatile resistive state in MoS₂ atomristors. Based on DFT and M3GNet, the analysis from electric field coupled molecular dynamics, equilibrium molecular dynamics, supercell relaxations, total energy and CINEB calculations reveal that such a state is neither stable nor non-volatile but an artefact of ReaxFF. On the other hand, Au atoms adsorbed at sulfur vacancies were found to be stable and non-volatile. The CINEB calculations also revealed that the Au atom migration energies from gold electrode to adsorption at sulfur vacancies depend significantly on the Au atom – V_S distances. Additionally, it is observed that different local heating rates during reset operation cause Au atom desorption at different electric field values explaining the V_RESET cycle-to-cycle variability observed in experiments. Furthermore, the generated local heating causes local S vacancy diffusion post Au atom desorption. This local S vacancy diffusion together with the desorbed Au atom which will go back and sit randomly on Au electrode causes change of Au atom-V_S distances in the next set cycle. This results in a change of Au atom migration barriers, eventually leading to cycle-to-cycle V_SET variability in experiments.

Methods

DFT calculations

All DFT and MD calculations were carried out using QuantumATK²⁷. For structural relaxations carried out using DFT, we utilized spin polarized GGA-PBE exchange-correlation functional alongwith state of the art norm-conserving SG15 pseudopotential to describe the core electrons²⁸. The density mesh cutoff used was 120 Hartree with a k-point Monkhorst-Pack grid of 7 × 7 × 1. The structural relaxations using DFT were carried out until forces and stresses fell below 0.01 eV/Å and 0.1 GPa, respectively, and relaxations using ReaxFF¹⁴ and M3GNet²⁰ forcefields were carried out until atomic forces and stresses fell below 0.001 eV/Å and 0.05 GPa, respectively.

Structural details

For G-MoS₂-G (Au-MoS₂-Au) atomristors, initially an orthogonal MoS₂ supercell with an area of 2.21 × 2.24 nm² was created. The supercell had one S atom in popped state resulting in a defect density is 2.02 × 10¹³ cm⁻² which agrees with experiments^5,22. Then an interface was created between G (Au) and MoS₂ using interface builder²⁹ in QuantumATK. To accommodate lattice mismatch, a mean absolute strain of 1.28% (1.61%) was introduced in G (Au) with directional strains ϵ₁₁ = 1.4% (1.89%), ϵ₁₂ = −1.87% (2.6%), ϵ₂₂ = −0.56% (−0.34%). Then, G (Au) was replicated to other side of MoS₂ resulting in G-MoS₂-G (Au-MoS₂-Au) atomristor. For Au-MoS₂-G, after making Au-MoS₂ interface as mentioned above, G was introduced to other side. Also, a large vacuum of 70 Å was padded in the out-of-plane direction to avoid spurious interactions with periodic images. For Au contacted atomristors, each electrode on either side of MoS₂ was modelled using 4-9 Au layers. The initial distance between Au-MoS₂/G-MoS₂ interfaces was set at 2.7 Å³⁰ / 3.66 Å³¹ which post optimization was 2.87 Å / 4.0 Å. For Au-MoS₂ structures used in DFT-CINEB calculations, a mean absolute strain of 0.38% with ϵ₁₁ = 0.57%, ϵ₂₂ = 0.57%, ϵ₁₂ = 0% was introduced in Au. The resulting supercell had an angle of 60° between vectors, a relative surface rotation of 23.41° between Au and MoS₂ and an interfacial area of ~1.4 nm². The surface plane of Au used in all interface simulations was (111). For MoS₂ only and CINEB calculations, a vacuum of ≈ 35 Å was used in all simulations.

Equilibrium molecular dynamics

Molecular dynamic simulations were carried out while utilizing NPT-Martyna Tobias Klien and NVT-Berendsen thermodynamic ensembles using DFT (Ab-initio), M3GNet (Deep-learned) and ReaxFF forcefields. For Ab-Initio Molecular Dynamics (AIMD), similar DFT settings were utilized as mentioned in DFT details. All MD simulations were carried out at 300 K and 1 atm pressure (NPT) with a step size of 1 femtosecond for AIMD and 250 attoseconds for ReaxFF and M3GNet based MD. The initial velocities for atoms were taken from Maxwell–Boltzmann curve at 300 K. Equilibrium molecular dynamics was carried out for Figs. 4–6 and 8.

Electric field coupled molecular dynamics

Electric field coupled molecular dynamics were carried out independently using M3GNet and ReaxFF forcefields while utilizing NPT Martyna Tobias Klien barostat at 300 K and 1 atm pressure with a step size of 1 femtosecond. The initial velocities for atoms were taken from Maxwell–Boltzmann curve at 300 K. Each system was relaxed and equilibrated for 100 picoseconds before applying electric field. The partial charges (q) on all atoms were calculated dynamically at each step of MD using the charge equilibration (QEq) scheme of Rappe and Goddard³². Then force on each atom was calculated in the direction of electric field (vec{E}) as (overrightarrow{{F}_{{Ext}}}=qtimes overrightarrow{E}). Next, all atomic forces were updated at each step of MD using (overrightarrow{{F}_{{Net}}}=overrightarrow{{F}_{0}}+overrightarrow{{F}_{{Ext}}}), where (vec{{F}_{0}}) includes forces due to forcefield alone. The average partial charges obtained from QEq scheme for sulfur (S) and molybdenum (Mo) were 0.77 and −0.385, respectively, whereas those reported using DFT bader charges³³ are 1.72 and −0.86, respectively. Due to the underestimation of charges, the values of electric fields used in our simulation may appear higher. Moreover, in an actual system the charge at a site is influenced by various factors such as local potential, strain, defects, substrate, residues etc. Therefore, it is quite difficult to correlate experimental electric fields values with the simulation environment. Nonetheless, it does not affect the conclusion of findings as evident from various agreements with DFT, experiments and previous studies. The local heating around Au adsorbed site in MoS₂ and electric field was implemented using posthook method in QuantumATK. Electric field coupling with molecular dynamics was introduced for Figs. 1–2 and 9–10.