High-throughput first-principles-calculations based estimation of lithium ion storage in monolayer rhenium disulfide

Two-dimensional materials are promising candidates for lithium ion battery anodes due to their large surface to volume ratio. The distorted T′ phase of the rhenium disulfide crystal makes the study of lithium binding more complex than for other two-dimensional materials with symmetric crystal structures. Here we explore the lithium ion storage capacity of monolayer rhenium disulfide by first-principles based calculations. We employ hardware-accelerator-assisted high-throughput calculations, using a van der Waals density-functional-theory based ‘structure search’ technique, to emulate the lithiation process. Exploring 2000 structures, each containing 49 to 98 atoms, we find the most stable lithiated structures for various lithium concentrations. We then design a delithiation algorithm and apply it to those lithiated structures for the estimation of the reversible specific capacity. Despite possessing high molar mass, a reasonably high specific capacity (214.13 mAh/g) and open-circuit voltage (0.8 V), in agreement with experimental results, make rhenium disulfide a promising alternative anode material. Rhenium disulfide is a promising lithium ion battery material but its distorted structure makes computational modelling challenging. Here hardware-accelerator-assisted high-throughput DFT based structure searching is used to model the reversible lithiation of ReS2 including metal–sulphur bond cleavage.

R henium dichalcogenide, a recent addition to the twodimensional (2D) transitional metal dichalcogenides (TMDs), has warranted special attention due to its unique crystal structure. It exists in naturally stable distorted octahedral (T′) phase, where the metal (Re) atoms tend to cluster together giving rise to a triclinic symmetry with two in-plane distinct regions: wide inter-cluster and narrow intra-cluster. Specifically, the metal atoms arrange themselves in such a way that four clustered metal atoms can be visualized as a diamond (rhombus), which has coined up the term 'diamond-shaped phase' in literature 1 . Such distorted crystal structure has prompted strongly anisotropic electronic 1 and thermal transports 2 and even anisotropic 1/f noise 3 in rhenium disulfide (ReS 2 ), which are not observed in other TMDs. Monolayer ReS 2 is available with commercial vendor and has also been synthesized recently by molten-salt-assisted chemical vapour deposition 4 . Moreover, the interlayer coupling energy is much weaker in ReS 2 compared to other TMDs, which makes the exfoliation of the monolayer much easier 5 . Potential application of ReS 2 has been demonstrated for digital electronics 6 , optoelectronics 7 and flexible electronics 8 . At the same time, 2D materials are being explored extensively for possible application in next-generation lithium-ion battery (LIB) anode due to their large surface-to-volume and surface-to-mass ratio 9,10 . Furthermore, the weak van der Waals forces between atomic layers allow Li ions to diffuse easily and without a significant volume expansion. Experiments conducted for several such 2D materials and their composites [11][12][13][14][15][16][17] have yielded encouraging results paving way to intensifying research in this domain. Recent work by Zhang et al. 17 has demonstrated a reversible specific capacity of about 200 mAh/g and a corresponding open-circuit voltage (OCV) of around 0.9 V in vertical ReS 2 nanowalls anode-based LIB, which is nearly similar to the metrics of MoS 2 anode-based LIBs 12 . This is indeed remarkable since the atomic mass of Re is almost double than that of Mo.
In parallel to these experimental efforts, first-principles based investigations using density functional theory (DFT) were applied to probe the physical and chemical origin of Li ion binding and diffusion on the surface of 2D materials. Since the pursuit for appropriate anode material has become a major challenge for the development of high-performance LIBs, such computational methods allow exploration of large set of 2D materials [18][19][20][21][22][23][24][25] , even which are not yet synthesized or have been synthesized under very stringent laboratory conditions. Theoretically estimated figure-of-merit parameters (e.g. specific capacity, diffusion barrier, open-circuit voltage etc.) for 2D material anode-based LIB could act as a guideline for the experimentalists. It also paves the way for appropriate functionalization [26][27][28][29] to improve the performance. Plethora of materials [18][19][20][21][22][23][24][25][26][27] investigated by firstprinciples based methods demonstrate the great potential of 2D materials for realization of high-performance LIBs. It is, however, be noted that all such studies deal with materials possessing highly symmetric crystal structure and 'uniform-adsorption' model is mostly used to study the Li storage capacity. In a recent study, Putungan et al. 30 demonstrated the limitation of 'uniformadsorption' model and its overestimation of the specific capacity for VS 2 based Na ion battery. To the best of our knowledge, no investigations on Li ion binding in ReS 2 based on first-principles have been reported yet, perhaps owing to the following reasons. (i) Due to distorted crystal structure, the number of atoms in ReS 2 unit cell is much higher 5,31,32 than any other 2D material with symmetric crystal structure. This imposes higher computation budget to study Li ion adsorption on ReS 2 . Geometry optimization also takes much longer time if the crystal lacks symmetry. In addition, the higher atomic number of Re makes computation with associated pseudopotential expensive. (ii) It is difficult to find single Li binding sites on distorted crystals, which further inhibits the calculation of the diffusion barriers. Furthermore, recent in situ 33 and ex situ 34 measurements demonstrate metal-sulphur bond breaking during lithiation and restoration during delithiation process in TMD anode-based LIBs. None of those earlier studies consider such bond breaking in their firstprinciples based calculations.
In this work we employ ab-initio random structure searching (AIRSS) technique 35 to investigate the lithiation process in monolayer ReS 2 . Hardware accelerator-assisted supercomputing facility (http://www.serc.iisc.in/facilities/cray-xc40-named-assahasrat/) is used to support the high-throughput ab-initio calculations required for large ReS 2 supercell. In a systematic approach we first identify single Li binding sites using AIRSS. This is followed by the calculation of the binding energies for each site and diffusion barriers between the sites by climbing image nudged elastic band (CINEB) 36,37 approach. We infer that for single Li the diffusion in ReS 2 surface is indeed anisotropic with widely varying diffusion barriers. We further use structuresearch technique to find most energetically stable structures for higher Li concentrations. We observe the typical increasing nature in adsorption energy characteristics for low Li concentrations; however, a monotonically decreasing adsorption energy for higher Li concentrations is seen once the Re-S bonds start breaking. At the same time, we notice a sharp decline in cohesive energy of the lithiated structures with increasing Li concentration. In contrast to previous studies, we propose an algorithm to emulate delithiation process in order to examine the reversibility of Re-S bond breaking, and to estimate the reversible specific capacity and average open-circuit voltage. We find that our theoretically predicted values of those parameters are in agreement with experimental results 17 . To draw a comparison, using the same method we calculate the reversible specific capacity of MoS 2 and find it to be very close to that of ReS 2 , despite Mo atoms being much lighter than Re. We use crystal orbital hamilton population (COHP) [38][39][40][41] technique to explain this observation. Our study predicts ReS 2 to be a reasonable candidate for LIB anode application. The proposed methods, which combines structure-search-based lithiation processes followed by delithiation, may be applied to other 2D materials to assess their suitability as LIB anodes.

Results
Study of single Li adsorption and diffusion barriers. Figure 1a and b illustrate the top and side view of a fully relaxed 2×2 unit cell of ReS 2 with lattice parameters a = 6.42 Å, b = 6.52 Å, α = 91.88, β = 103.49̊and γ = 118.85̊(see Supplementary Note 1, Supplementary Figure 1), which agree with the previous reports 1, 31,32,42 . The pink shaded region indicates the intracluster area where clustering of Re atoms can be observed. The white space in between the shaded regions indicates the wider inter-cluster region. Also, in the intra-cluster region four Re atoms can be seen making a rhombus-like shape, which ultimately forms 'diamond chains' in this region.
For TMDs with symmetric H or T structures, one can easily find the (x, y) co-ordinates of the adsorption sites in a unit cell by simple intuition. When adsorbed, the Li ion usually is triplecoordinated by three chalcogen atoms, giving rise to only two kinds of adsorption sites on a unit cell: on top of the metal atom ('top' from now on) and on the centre of the hexagon ('hexagon' hereafter). The top site is energetically more favourable since the proximity to the metal atom allows the Li ion to bond with it apart from the three neighbouring chalcogens. However, in the case of ReS 2 , due to the highly distorted T′ structure, it is indeed difficult to locate the adsorption sites simply by intuition. Thus, it becomes essential to use 'structure-search' technique to study Li binding on the ReS 2 surface. Among different techniques reported for 'structure-search' (e.g. particle swarm optimization 43 , evolutionary algorithm 44 etc.), here, we use ab-initio random structure searching (AIRSS), since this method is highly scalable (allows two levels of parallelization using hardware accelerator as discussed in Methods) and free of any bias (i.e. minimizes any local minima trapping). AIRSS has already been widely used to predict high-pressure phases of SiH 4 45 , dense hydrogen phases 46 , incommensurate phases in aluminium 47 , ionic phases of ammonia 48 and also for metal ion binding on symmetric TMDs 30,49 .
For this study, we use a 2×2 ReS 2 supercell containing 16 Re and 32 S atoms to eliminate all possible Li-Li interactions from periodically repeated neighbouring cells. All the possible single Li adsorption sites as obtained from AIRSS are shown in Fig. 1c. Similar to previous studies on symmetric TMD crystals, we classify these Li binding sites into two categories: α and β, which are somewhat equivalent to top and hexagon sites, respectively. The equilibrium distances between the adsorbed Li ions and the nearest Re atoms are listed in column 4 of Table 1, which agree with this classification. Interestingly, we find only six stable Li adsorption sites in an ReS 2 unit cell, while an equal sized symmetric TMD supercell would contain eight such spots: four top and four hexagon sites. We find that due to Re atom clusterization there is no β site in the intra-cluster regions. However, the wider inter-cluster regions do contain β sites in the hollow space of the distorted 'hexagons'. It is worth noting that the adsorbed Li are indeed triply coordinated by S atoms, albeit with widely different equilibrium distances as listed in column 2 of Table 1. However, not all triply S coordinated sites are stable adsorption sites as we can see in case of intra-cluster regions.
Column 9 of Table 1 lists the amount of charge transfer from Li to ReS 2 for different binding sites obtained from Bader charge analysis [50][51][52] . We see that for all adsorption sites the charge transfer is very close to one electron, indicating that the Li atom would transfer its entire 2s 1 electron to the substrate and will exist in an ionized state after adsorption (see Supplementary Note 2, Supplementary Figure 2). We define the adsorption energy of single Li atom in ReS 2 , E ads as where E ReS 2 þLi is the energy of Li adsorbed ReS 2 system, E ReS 2 is the energy of the pristine monolayer ReS 2 substrate, and μ Li is the chemical potential of Li. As per this definition, more negative adsorption energy indicates tighter Li binding while positive adsorption energy is an indication of the possibility of Li clustering and phase separation with the substrate. To be consistent with literature, we calculate the adsorption energies using both bulk and atomic values of μ Li , which are listed in columns 7 and 8 of Table 1. We see that the adsorption energy for all the binding sites are well below zero even when referenced with bulk Li, indicating thermodynamically stable Li adsorption on ReS 2 . The Li binding energies of ReS 2 is comparable or superior (more negative) to other widely studied TMDs like MoS 2 , MoSe 2 , WS 2 , WSe 2 and their heterostructures 24 , and other 2D materials like phosphorene 20 and borophene 22 . It is very interesting to note that all the six adsorption spots in ReS 2 have widely different adsorption energies. To explain this observation, we use two techniques: first, we measure the amount of distortion induced by lithiation as the sum of the displacement of all Re and S atoms (see Supplementary Note 3). The binding site with least adsorption energy would have the least distorted structure, which here would be the α4 site. We take this as the reference structure, since the lattice parameters of all the lithiated cells with different binding sites are almost equal. The computed values of total displacement of Re and S atoms relative to the lithiated cell where Li is adsorbed at the α4 site are given in column 6 of Table 1. We observe that indeed there is a monotonic relation between adsorption energy and the amount of distortion introduced in the material by the lithiation process. We further use recently  Table 1, respectively, which confirms the theory that stiffer binding sites will result in loose binding. For instance, at α1 site, the Li faces least stiffly attached S atoms where the Li binding is greatest; similarly, at α4 site, the Li sees most stiff S atoms resulting in the loosest binding. At the α1 and β2 spots, the Li faces exactly equally stiff S atoms, but at the α1 site, the Li binding is much stronger, because the proximity to Re atom at this spot allows the Li to bind with it, which is not the case at β2. Same argument is applicable in case of α2 and β1 sites. Diffusion barrier height faced by Li ions while diffusing on the surface of the material is another extremely important property for an anode material since it decides how fast the LIB can be charged or discharged. We employ CINEB calculations to determine the activation energies for the minimum energy pathways (MEPs) between stable adsorption sites. Since each of these six binding sites have different binding energies, estimation of diffusion barrier is more complex than other symmetric 2D materials. In case of ReS 2 , we found 10 unique and several nonunique MEPs (see Supplementary Note 4, Supplementary  Figures 3-11) between neighbouring adsorption sites, which make the CINEB calculations very expensive. We overcome this hurdle using highly parallel computation in hardware accelerators. Figure 2a shows the relative energy profiles of all possible unique MEPs between α and β sites, which again varies widely due to the highly anisotropic crystal structure of ReS 2 . Figure 2b illustrates the energy barriers of all possible unique MEPs between two α sites. Note that in symmetric TMD materials there is no direct MEP between two 'top' sites, because each hexagon always contains a 'hexagon' site which takes the MEP through it. Absence of β sites (equivalent to 'hexagon' sites) in the intra-cluster region allows direct MEPs between two α sites in that region of ReS 2 . It should further be noted that, the barrier height between two α sites in the intra-cluster region is usually larger than that of α-β pathways, making Li diffusion much less likely in that direction since the diffusion coefficient in a specific direction depends on the negative exponential of the activation energy in that direction. All these complex MEPs can be summarized using a directional graph-like schematic, which is illustrated in Fig. 2c. Due to the periodic boundary condition, there could be two completely different MEPs between two binding spots. For instance, the α1-α2 MEP has a barrier height of 0.4 eV as the Li is diffuses to the right, but the MEP between the same sites has an activation energy of 0.42 eV while the diffusion happens to the left. Because of the absence of β sites and much higher diffusion barrier, the possibility of long range diffusion in the direction of b through the intra-cluster space is expected to be very low. The Li ion will mostly diffuse through the inter-cluster region along the direction of a in a sort of zig-zag motion, which indicates anisotropic Li diffusion in ReS 2 . However, this anisotropic diffusion will follow a zig-zag path instead of a relatively straight line as in the case of phosphorene 20 or borophene 22 , where, in one direction the diffusion barrier is almost 1-3 order less than other directions. The α-β MEPs, where most of the Li diffusion is expected to occur, have diffusion barrier heights comparable to or less than those of other widely studied TMDs like MoS 2 , MoSe 2 , WS 2 , WSe 2 , and their heterostructures 24 .  30 observed in their study of Na adsorption on VS 2 . For x = 2.5 and 3.125, the amount of Re-S bond breaking becomes much more severe, turning the lowestenergy lithiated phases almost amorphous. To the best of our knowledge, so far, no first-principle-based study of Li ion storage on TMD has included the effect of metal-chalcogen bond breaking. However, recent experimental studies have established this to be a real-life phenomenon. Su et al. 33 observed in their in situ transmission electron microscopic (TEM) study that single-crystalline MoS 2 nanosheets convert to Mo nanograins with broken Mo-S bonds which helps in forming an Li 2 S matrix for the nanograins after the first full lithiation. After the delithiation, they observed the disappearance of the Li 2 S phase and We define the adsorption energy (E ads ) for the multi-Li adsorbed ReS 2 phases as where E ReS 2 þLi is the energy of n number of Li ion adsorbed ReS 2 system, E ReS 2 is the energy of monolayer ReS 2 substrate, and μ Li is the chemical potential of Li (bulk). Figure 4a illustrates the variation of the adsorption energy with Li content per formula unit (x). Within x 2 [0, 1], we mostly observe a typical increment towards zero, and at x = 1.0 the adsorption energy becomes 0.099 eV/Li, closest to zero. This behaviour can be attributed to the fact that with increasing adsorbed Li concentration the repulsion between the Li ions increases. A local minima at x = 0.25 can be explained as follows. Increasing Li ion concentration introduces more and more distortions in ReS 2 . As discussed in the previous section, such deformation in the substrate can lead to stronger Li ion binding. Due to the large size of the supercell and equal Li ion distribution on both sides of ReS 2 , at very low values of x (0.0625, 0.125 and 0.25), there is insignificant coulombic repulsion between the Li ions. Thus, we observe decreasing adsorption energy (i.e., stronger binding) with increasing x. However, for x ! 0.5, the repulsion between the Li ions becomes significant and the adsorption energy starts to increase towards zero in a typical fashion. From x = 2, the adsorption energy becomes more negative with increasing Li content, which is a completely opposite trend of what was observed for x 1.0. It is obvious that the breaking of Re-S bonds and formation of Li-S compound aids in stronger adsorption of Li ions into the material. It is probably because , which can be written as where E Li x ReS 2 is the energy of the lithiated ReS 2 phase with x Li ion per formula unit, E Li , E Re and E S are energies of free Li, Re and S atoms respectively. By this convention, more negative cohesive energy indicates better cohesion of the structure. Figure 4b shows the variation of the cohesive energy with increasing Li content per formula unit. Unlike adsorption energy, the cohesive energy steadily declines (increases towards zero) with increasing Li concentration, which is not surprising considering the fact that the lowest-energy lithiated phases tend to move towards amorphous phases with increasing Li content.  Figure 14) and find that the latter is inferior in terms of both adsorption and cohesive energy.
The maximum amount of Li that can be adsorbed into a material (x max ) is a very important quantity, which explicitly decides metrics like specific capacity and OCV. In most previous studies, x max is determined as the value of x for which the value of adsorption energy is closest to zero, i.e., further increase in the value of x would result in positive adsorption energy. Unlike uniform adsorption model where the adsorption energy profile is monotonically increasing, our study finds a combination of two completely opposite trends. However, despite being more stable in terms of adsorption energy the almost-amorphous lithiated phases with high amount of Li are very weakly bound together as the cohesive energy indicates. At x = 3.125 the cohesive energy is only −3.94 eV/atom, which suggests that such phase might disintegrate on its own at slightly high temperatures. Clearly, trade-off between the adsorption and cohesive energy needs to be considered while deciding the value of x max . Since the experimental studies suggest that the lithiation and successive delithiation process must be reversible, we define x max as the maximum value of x for which the lithiated phase will return close to the pristine substrate crystal structure after delithiation, both visually and energetically.
Study of delithiation. Inspired by previously proposed delithiation algorithms for silicon 54-56 , we develop a simple but rigorous delithiation technique for our case (see Methods). Figure 5 illustrates the delithiated structure of the phase Li 2.5 ReS 2 , which clearly is very far from the pristine ReS 2 structure. On the other hand, delithiated structure of the phase Li 2 ReS 2 is visually as well as energetically very close to the pristine ReS 2 structure where all but two broken Re-S bonds have been restored. Supplementary  Movies 1-4 show the gradual structural evolution during the delithiation process. From the above studies, we thus take x max = 2.0 in case of ReS 2 . The reversible specific capacity (C) of ReS 2 has been calculated from the following equation: where MW ReS 2 is the molecular weight of ReS 2 , F ¼ 26:801 Ah/ mol is the Faraday's constant, and v ¼ 1 is the number of valence electrons for Li. Note that the former is calculated with two references: bulk value of chemical potential which is taken as cohesive enery of metallic BCC Li, and atomic value which is taken as the energy of an isolated Li atom We also calculate the average OCV using the following equation: where x min = 0 and x max = 2.0 are the amount of Li per formula unit of ReS 2 in initial and final structure after lithiation, respectively, and e is the electronic charge.
We find, using above equations, the value of C to be 214.13 mAh/g and the OCV turns out to be 0.8 V. These values of reversible specific capacity and average OCV are in very close agreement with experimentally reported values 17 . Note that, the value of C obtained here is the reversible specific capacity that can be much lesser than the specific capacity observed in the first few cycles.
Comparison with MoS 2 . Since MoS 2 probably is the most studied among all TMDs, it would be interesting to compare the metrics obtained for ReS 2 with those of MoS 2 using the same methodology. Uniform adsorption model predicts x max > 2 in case of MoS 2 25 . However, with the proposed methodology we find that the delithiated lowest-energy Li 2 MoS 2 structure is indeed far from the pristine MoS 2 structure. However, the lowest energy Li 1.625 MoS 2 phase seems to return very close to the pristine structure of MoS 2 after the delithiation (see Fig. 6). Thus, for MoS 2 , we take x max = 1.625, which yields reversible specific capacity, C = 272.08 mAh/g, and average OCV = 0.74 V. Again, these values match the experimentally reported values 12 closely, and are comparable to those of ReS 2 despite MoS 2 being more than 1.5 times lighter. To investigate the origin of lower x max for MoS 2 we conduct COHP analysis.
It is observed for both materials that most of the metal-sulphur bonds broken due to lithiation are restored after delithiation. However, the breaking of the 'metal backbone' of the TMDs due to lithiation is completely irreversible. The susceptibility to breakage of the metal backbone due to lithiation is determined by the strength of the metal-metal bonds in the material. This is a very important parameter as according to the model used in this study, the x max and associated specific capacity and OCV are directly dependent on it. Table 2 lists all the metal-S as well as metal-metal atomic distances and their average ICOHP values as well as cohesive energies for both materials. In case of ReS 2 , intra and inter-cluster values are computed separately because of the significant structural difference in these two regions.
It turns out that in the intra-cluster region of ReS 2 , Re-Re bonds are more than twice as strong compared to Mo-Mo bonds in MoS 2 . However, in the inter-cluster region there is almost no bonding between the Re atoms. Apparently, the clusters are then held together by only inter-cluster Re-S bonds. In the intracluster region of ReS 2 , the Re-S bonds are again stronger than Mo-S bonds of MoS 2 , but in the inter-cluster region the Re-S bonds are slightly weaker than Mo-S bonds. Thus, any mechanical failure like crack propagation will occur at the weaker inter-cluster region and separation of the clusters of ReS 2 will be easy. In fact, this can be observed in the case of delithiated Li 2.5 ReS 2 structure where the clusters have been separated from each other due to the deformation introduced by the lithiation process. All this bond-specific information obtained from COHP analyses can be roughly summarized into one quantity, the cohesive energy of the TMD. Put simply, since ReS 2 has slightly stronger cohesion than MoS 2 , its reversible Li storage capacity is higher.

Discussion
We propose a first-principle-based method for the estimation of reversible specific capacity and average open-circuit voltage of 2D materials-based Li battery anode. In contrast to earlier studies, here the maximum Li storage capacity is estimated by randomstructure-search-based lithiation followed by random-Liremoval-based gradual delithiation algorithm. Though computationally expensive, proposed method can capture the bondbreaking phenomenon during lithiation and successive restoration during the delithiation process. Thus, it can predict realistic values of the reversible specific capacity, which can be further improved with the inclusion of solvent and solid-electrolyte

Methods
DFT calculation. DFT calculations have been carried out using generalized gradient approximation (GGA) as implemented in the code VASP 57-60 with PAW 61 method using the Perdew-Burke-Ernzenhof (PBE) 62 exchange-correlation functional. The following electrons have been treated as valence electrons and are expanded in the plane wave basis set: Li: 1s 2 2s 1 , Re: 5p 6 6s 2 5d 5 , Mo: 4p 6 5s 1 4d 5 and S: 3s 2 3p 4 . A high cut-off energy of 700 eV was required for this expansion to avoid Pulay stress. A Γ-centred 5×5×1 k-mesh was found to be suitable to sample the Brillouin zone for structural relaxations and a denser k-mesh of 11×11×1 has been used for static energy calculations. Electronic convergence is achieved when the difference in energy of successive electronic steps becomes less than 10 −6 eV, whereas the structural geometry is optimized until the maximum force on every atom falls below 0.01 eV/Å. A large vacuum space of more than 25 Å in the direction of c is applied to avoid any interaction between successive layers. Semiempirical dispersion corrections with DFT-D3 method as developed by Grimme 63 is used in all the calculations. For CINEB calculations 3 to 12 images were used depending on the convergence speed. Calculation of CINEB was found to be much faster in hardware accelerators than in conventional CPU. All crystal structure images were generated using the tool VESTA 64 .
Hybrid processor system based AIRSS calculation. n/2 Li ions were placed on each side of the pristine ReS 2 supercell randomly using the AIRSS initial configuration generation engine and then relaxed using DFT. The (x, y) fractional coordinates of the Li ions were allowed to be completely random while keeping the minimum separation between any two atoms 1.8 Å to minimize atomic overlaps. The vertical distances (along c) of the Li ions from the surface was chosen randomly from the range of 1.5-4.5 Å. For each value of n, 200 random structures that could be relaxed successfully were explored. A detailed description of this AIRSSlowest-energy-phase-finding algorithm along with schematics can be found in Supplementary Note 8 and Supplementary Figures 15-16. AIRSS algorithm can be highly parallelized as every random initial configuration is generated independently of other structures. The two-level loop in our AIRSS algorithm can be completely 'flattened' (i.e., performed parallelly), provided enough processors are available. To take advantage of the parallel computing feature of AIRSS algorithm, we implemented it in a hybrid master-slave CPU-accelerator architecture that resulted in high-throughput ab-initio computations as illustrated in Fig. 7. To accelerate submodules of DFT calculations (e.g. fast Fourier transforms), GPU version of VASP 65,66 code is used, and regular CPU code is adapted for execution in Xeon-Phi 67 . Depending on the number of atoms in the supercell, for AIRSS structure relaxation, 2×2×1 or gamma-only k-points sampling was used with Li: 2s 1 electron and reduced cut-off energy requirement (340 eV). However, top five lowest-energy structures from every rank-list was processed further with the above-mentioned stringent DFT calculation setup. It is worth noting that previous efforts on structure discoveries 47 were limited to few atoms per structure. Use of hardware accelerators thus extends the structure-search technique to supercells with very large number of atoms. A comparison of computation time between ReS 2 and MoS 2 structures can be found in Supplementary Note 9.
COHP. COHP analyses have been performed using LOBSTER 68 . The wavefunctions generated by VASP deploying an 11×11×1 (for pristine TMD) or an 8×8×1 kmesh with the stringent DFT conditions mentioned above is used as an input to the tool. Unlike structural relaxation, the COHP calculation demands much more computer memory and hence it is conducted using all-CPU cluster at nano-scale device research laboratory.
Delithiation algorithm. For a given lithiated structure, single Li ion was randomly removed from each of the top and bottom surface at a time, and then the entire structure was relaxed where atomic positions as well as the lattice parameters were allowed to change. This process was repeated until zero Li atoms were left in the structure. Unlike AIRSS based lithiation process, the delithiation algorithm is sequential in nature, and thus computationally expensive. There could be numerous delithiation paths depending on which Li atoms are being removed, however, we only explored three such delithiation paths for each structure due to computational budget. In Figs 5 and 6 we show the final structures from one such path and the remaining are shown in Supplementary Figure 17. All the three random delithiation paths provide similar results. We further show that from a lithiated structure, if all the Li atoms are removed at once and then relaxed using a constrained volume, which we call 'quick delithiation' (see Supplementary Note 10 and Supplementary Figure 18), we get quite similar structure as obtained from the previous algorithm. A combination of quick and gradual delithiation algorithm could act as a rigorous technique for the estimation of specific reversible capacity of any material.

Data availability
The authors declare that the main data supporting the findings of this study are available within the paper and its Supplementary Information file. Other relevant data are available from the corresponding author upon reasonable request.