The empirically fitted set of parameters listed in Table III were used to carry out a series of MD simulations in the temperature range of 400 K to 600 K with 25 K intervals, and zero atmospheric pressure. The Parrinello-Rahman isobaric-isothermal (NPT) MD method , which allows for changes in the simulation box sizes while keeping angles fixed, was used. The temperatures and pressure in the system were controlled using thermostating and barostating techniques, whereby some dynamic variables are coupled with particle velocities and simulation box dimensions. Simulations commenced from the reported X-ray diffraction structures wherein the simulation supercell was constructed by 7 × 4 × 3 unit cells comprising 3696 atoms in an orthorhombic symme-try following P2 1 2 1 2 1 (No. ). Several partially occupied crystallographic sites were identified in the Na conduction plane. In MD simulations, we placed the Na + ions at the crystallographically reported highest occupied sites (Na2 and Na5) to avoid strong cation-cation repulsion. For clarity, all reported crystallographically distinct sites (Na1, Na2, Na3, Na4, Na5, and Na6) in the layered frameworks of Na 2 LiFeTeO 6 are shown in Figure . The inter-layer interaction is usually weaker in nature; thus, inter-layer sliding is expected unless the initial configuration is not set properly. To avoid inter-layer sliding, we performed a constrained MD simulation for 500 ps for each temperature case; the framework atoms (Li, Fe, and Te) were frozen at the crystallographic sites and the other atoms were allowed to move. Finally, we lifted the former constraint and performed the usual MD simulation. A time step of 2 fs was applied to the velocity Verlet algorithm to solve Newton equations of motion. A typical run-time of 6 ns or longer was used with trajectory samples stored at intervals of 200 fs for detailed analyses. To guarantee the thermodynamic convergence properties, a few longer run-time simulations of 100 ns were performed separately. Periodic boundary conditions in all three directions and the Ewald summation technique for the convergence of long-range Coulombic interactions were applied using the LAMMPS software package . A cut-off distance of 13 Å was used for both the short-range interactions and the short-range part of the Ewald summation. Micro-canonical MD (NVE-MD) simulations were further performed at 600 K using the final structure obtained from NPT-MD simulations.