Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles and JavaScript.

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

## Abstract

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.

## Introduction

Rhenium dichalcogenide, a recent addition to the two-dimensional (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 literature1. Such distorted crystal structure has prompted strongly anisotropic electronic1 and thermal transports2 and even anisotropic 1/f noise3 in rhenium disulfide (ReS2), which are not observed in other TMDs. Monolayer ReS2 is available with commercial vendor and has also been synthesized recently by molten-salt-assisted chemical vapour deposition4. Moreover, the interlayer coupling energy is much weaker in ReS2 compared to other TMDs, which makes the exfoliation of the monolayer much easier5. Potential application of ReS2 has been demonstrated for digital electronics6, optoelectronics7 and flexible electronics8. 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 ratio9,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 composites11,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 ReS2 nanowalls anode-based LIB, which is nearly similar to the metrics of MoS2 anode-based LIBs12. 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 materials18,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 functionalization26,27,28,29 to improve the performance. Plethora of materials18,19,20,21,22,23,24,25,26,27 investigated by first-principles 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 ‘uniform-adsorption’ model and its overestimation of the specific capacity for VS2 based Na ion battery. To the best of our knowledge, no investigations on Li ion binding in ReS2 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 ReS2 unit cell is much higher5,31,32 than any other 2D material with symmetric crystal structure. This imposes higher computation budget to study Li ion adsorption on ReS2. 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 situ33 and ex situ34 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 first-principles based calculations.

In this work we employ ab‒initio random structure searching (AIRSS) technique35 to investigate the lithiation process in monolayer ReS2. Hardware accelerator-assisted supercomputing facility (http://www.serc.iisc.in/facilities/cray-xc40-named-as-sahasrat/) is used to support the high-throughput ab‒initio calculations required for large ReS2 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 ReS2 surface is indeed anisotropic with widely varying diffusion barriers. We further use structure-search 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 results17. To draw a comparison, using the same method we calculate the reversible specific capacity of MoS2 and find it to be very close to that of ReS2, 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 ReS2 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 ReS2 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 reports1,31,32,42. The pink shaded region indicates the intra-cluster 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 triple-coordinated 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 ReS2, 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 ReS2 surface. Among different techniques reported for ‘structure-search’ (e.g. particle swarm optimization43, evolutionary algorithm44 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 SiH445, dense hydrogen phases46, incommensurate phases in aluminium47, ionic phases of ammonia48 and also for metal ion binding on symmetric TMDs30,49.

For this study, we use a 2×2 ReS2 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 ReS2 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 ReS2 for different binding sites obtained from Bader charge analysis50,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 2s1 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 ReS2, Eads as

$$E_{{\mathrm{ads}}} = E_{{{{\rm{ReS}}_2} + {\mathrm{Li}}}} - E_{{\mathrm{ReS}_{2}}} - \mu _{{\mathrm{Li}}}$$
(1)

where $$E_{{{{\rm{ReS}}_2 + {\rm{Li}}}}}$$ is the energy of Li adsorbed ReS2 system, $$E_{{{{\rm{ReS}}_2}}}$$ is the energy of the pristine monolayer ReS2 substrate, and $$\mu _{{\mathrm{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 ReS2. The Li binding energies of ReS2 is comparable or superior (more negative) to other widely studied TMDs like MoS2, MoSe2, WS2, WSe2 and their heterostructures24, and other 2D materials like phosphorene20 and borophene22. It is very interesting to note that all the six adsorption spots in ReS2 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 proposed COHP technique, which partitions the density of states (DOS) into bonding, nonbonding and antibonding contributions, and the presence of more antibonding states at the Fermi level are considered destabilizing53. The integration of the COHP up to the Fermi level (ICOHP) gives a reasonable estimate of the strength of an individual bond (i.e., more negative ICOHP values by convention indicate stronger bonds). We compute the average ICOHP values for the S and Re atoms with their nearest neighbours for different adsorption sites in the pristine ReS2 structure, i.e. before lithiation. These values tell us how stiffly the S and Re atoms are bound in the crystal at different possible Li binding sites. Average ICOHPS and ICOHPRe values are tabulated in columns 3 and 5 of 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 ReS2, we found 10 unique and several non-unique MEPs (see Supplementary Note 4, Supplementary Figures 311) 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 ReS2. 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 ReS2. 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 ReS2. However, this anisotropic diffusion will follow a zig-zag path instead of a relatively straight line as in the case of phosphorene20 or borophene22, 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 MoS2, MoSe2, WS2, WSe2, and their heterostructures24.

### Adsorption and storage of higher concentration of Li

Using AIRSS we find lowest-energy lithiated phases for the compound LixReS2, for x = 0.125, 0.25, 0.5, 1.0, 2.0, 2.5 and 3.125. Again, a 2×2 ReS2 unit cell is chosen as the substrate which allows large degrees of freedom while keeping the computational cost manageable even with a very high number (like n = 40 or 50; n = 16x) of Li ions in the supercell. Figure 3 shows lowest-energy lithiated ReS2 structures for x = 0.5, 1.0, 2.0 and 2.5 as found by AIRSS calculations, while lowest-energy structures with x = 0.125, 0.25 and 3.125 are illustrated in Supplementary Note 5 and Supplementary Figure 12. As expected, increasing the concentration of Li introduces distortion in the crystal structure of ReS2. It can be observed that till x = 1.0, there is no explicit Re–S bond breaking in the lowest-energy phases though at x = 1.0 an Re atom can be seen moving upwards along with the attached S atoms. However, at x = 2.0 significant amount of Re–S bond breaking, and Li2S-like compound forming can be seen. Also, a few floating and clustered Li ions can be observed at this concentration. This is similar to what Putungan et al.30 observed in their study of Na adsorption on VS2. For x = 2.5 and 3.125, the amount of Re–S bond breaking becomes much more severe, turning the lowest-energy 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 MoS2 nanosheets convert to Mo nanograins with broken Mo–S bonds which helps in forming an Li2S matrix for the nanograins after the first full lithiation. After the delithiation, they observed the disappearance of the Li2S phase and restoration of the broken Mo–S bonds which led to the re-formation of MoS2 nanocrystals. In other ex situ TEM study, Zhang et al.34 observed co-existence of NbS2, NbS and Li2S after charging the LIB with a codoped NbS2 electrode.

$$E_{{\mathrm{ads}}} = \frac{{E_{{\mathrm{ReS}_{2} + {\mathrm{Li}}}} - E_{{\mathrm{ReS}_{2}}} - n\mu _{{\mathrm{Li}}}}}{n}$$
(2)

where $$E_{{{{\rm{ReS}}_2 + {\rm{Li}}}}}$$ is the energy of $$n$$ number of Li ion adsorbed ReS2 system, $$E_{{{{\rm{ReS}}_2}}}$$ is the energy of monolayer ReS2 substrate, and $$\mu _{{\mathrm{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 \in$$ [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 ReS2. 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 ReS2, 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 \ge$$ 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\, \le$$ 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 of the Li–S compound formation, where the S atoms bonded with the Li ions screen the ionic charges and reduce the columbic repulsions between them. It also gives Li ions the opportunity to bond with the Re ions with dangling bonds which further stabilizes the system. With increasing Li content more and more S atoms are liberated from Re atoms, screening the positive charges of the Li ions entirely. In fact, at x = 3.125 the adsorption energy takes its lowest value, −1.17 eV/Li. Another very important metric for a structure’s stability is the cohesive energy (Ecohesive), which can be written as

$$E_{{\mathrm{cohesive}}} = \frac{{E_{{\mathrm{Li}}_x{\mathrm{ReS}}_2} - xE_{{\mathrm{Li}}} - E_{{\mathrm{Re}}} - 2E_{\mathrm{S}}}}{{x + 1 + 2}}$$
(3)

where $$E_{{{{\rm{Li}}_x{\rm{ReS}}_2}}}$$ is the energy of the lithiated ReS2 phase with $$x$$ Li ion per formula unit, $$E_{{\mathrm{Li}}}$$, $$E_{{\mathrm{Re}}}$$ and $$E_{\mathrm{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. COHP analyses (see Supplementary Note 6 and Supplementary Figure 13) of lithiated structures support this observation. From this analysis it can also be observed that at around x = 1.0 the band gap of the composite material reduces to almost zero, indicating metallic conductivity. This happens because the 2s1 electrons are entirely transferred from the Li atoms to the substrate as our Bader charge analysis revealed earlier. For completeness, we compare the atomic model obtained from AIRSS with conventional ‘uniformly adsorbed model’ (see Supplementary Note 7 and Supplementary 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 (xmax) is a very important quantity, which explicitly decides metrics like specific capacity and OCV. In most previous studies, xmax 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 xmax. Since the experimental studies suggest that the lithiation and successive delithiation process must be reversible, we define xmax 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 silicon54,55,56, we develop a simple but rigorous delithiation technique for our case (see Methods). Figure 5 illustrates the delithiated structure of the phase Li2.5ReS2, which clearly is very far from the pristine ReS2 structure. On the other hand, delithiated structure of the phase Li2ReS2 is visually as well as energetically very close to the pristine ReS2 structure where all but two broken Re–S bonds have been restored. Supplementary Movies 14 show the gradual structural evolution during the delithiation process. From the above studies, we thus take xmax = 2.0 in case of ReS2. The reversible specific capacity (C) of ReS2 has been calculated from the following equation:

$$C = \frac{1}{{MW_{{{{\rm{ReS}}_2}}}}}(x_{{\mathrm{max}}}vF\times 10^3)$$
(4)

where $$MW_{{{{\rm{ReS}}_2}}}$$ is the molecular weight of ReS2, $$F = 26.801$$ Ah/mol is the Faraday’s constant, and $$v = 1$$ is the number of valence electrons for Li.

We also calculate the average OCV using the following equation:

$${\mathrm{OCV}} \approx - \frac{{\left[ {E_{{\mathrm{Li}}_{x_{{\mathrm{max}}}}{{{\rm{ReS}}_2}}} - E_{{\mathrm{Li}}_{x_{{\mathrm{min}}}}{{{\rm{ReS}}_2}}} - \left( {x_{{\mathrm{max}}} - x_{{\mathrm{min}}}} \right)\mu _{{\mathrm{Li}}}} \right]}}{{\left( {x_{{\mathrm{max}}} - x_{{\mathrm{min}}}} \right)e}}$$
(5)

where xmin = 0 and xmax = 2.0 are the amount of Li per formula unit of ReS2 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 values17. 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 MoS2

Since MoS2 probably is the most studied among all TMDs, it would be interesting to compare the metrics obtained for ReS2 with those of MoS2 using the same methodology. Uniform adsorption model predicts xmax > 2 in case of MoS225. However, with the proposed methodology we find that the delithiated lowest-energy Li2MoS2 structure is indeed far from the pristine MoS2 structure. However, the lowest energy Li1.625MoS2 phase seems to return very close to the pristine structure of MoS2 after the delithiation (see Fig. 6). Thus, for MoS2, we take xmax = 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 values12 closely, and are comparable to those of ReS2 despite MoS2 being more than 1.5 times lighter. To investigate the origin of lower xmax for MoS2 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 xmax 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 ReS2, 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 ReS2, Re–Re bonds are more than twice as strong compared to Mo–Mo bonds in MoS2. 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 intra-cluster region of ReS2, the Re–S bonds are again stronger than Mo–S bonds of MoS2, 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 ReS2 will be easy. In fact, this can be observed in the case of delithiated Li2.5ReS2 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 ReS2 has slightly stronger cohesion than MoS2, 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 random-structure-search-based lithiation followed by random-Li-removal-based gradual delithiation algorithm. Though computationally expensive, proposed method can capture the bond-breaking 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 interphase effect. Computational techniques based on cutting edge processor technology is demonstrated to realize the high-throughput ab‒initio calculations efficiently. Potential of monolayer ReS2, a recent addition to TMD family with distorted T′ phase, as Li ion battery anode is estimated using the proposed methodology. In spite of possessing high molar mass, our study predicts a reasonably high reversible specific capacity and open-circuit voltage, which are in very close agreement with the experimental observation.

## Methods

### DFT calculation

DFT calculations have been carried out using generalized gradient approximation (GGA) as implemented in the code VASP57,58,59,60 with PAW61 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: 1s22s1, Re: 5p66s25d5, Mo: 4p65s14d5 and S: 3s23p4. 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 Grimme63 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 VESTA64.

### Hybrid processor system based AIRSS calculation

n/2 Li ions were placed on each side of the pristine ReS2 supercell randomly using the AIRSS initial configuration generation engine and then relaxed using DFT. The (x, y) fractional co-ordinates 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 AIRSS-lowest-energy-phase-finding algorithm along with schematics can be found in Supplementary Note 8 and Supplementary Figures 1516. 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 sub-modules of DFT calculations (e.g. fast Fourier transforms), GPU version of VASP65,66 code is used, and regular CPU code is adapted for execution in Xeon-Phi67. 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: 2s1 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 discoveries47 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 ReS2 and MoS2 structures can be found in Supplementary Note 9.

### COHP

COHP analyses have been performed using LOBSTER68. The wavefunctions generated by VASP deploying an 11×11×1 (for pristine TMD) or an 8×8×1 k-mesh 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.

## References

1. Lin, Y. C. et al. Single-layer ReS2: two-dimensional semiconductor with tunable in-plane anisotropy. ACS Nano 9, 11249–11257 (2015).

2. Jang, H., Ryder, C. R., Wood, J. D., Hersam, M. C. & Cahill, D. G. 3D anisotropic thermal conductivity of exfoliated rhenium disulfide. Adv. Mater. 29, 1700650 (2017).

3. Mitra, R., B. J., A. B. & A.D. Probing in-plane anisotropy in few-layer ReS2 using low frequency noise measurement. Nanotechnology 29, 145706 (2018).

4. Zhou, J. et al. A library of atomically thin metal chalcogenides. Nature 556, 355–359 (2018).

5. Tongay, S. et al. Monolayer behaviour in bulk ReS2 due to electronic and vibrational decoupling. Nat. Commun. 5, 3252 (2014).

6. Liu, E. et al. Integrated digital inverters based on two-dimensional anisotropic ReS2 field-effect transistors. Nat. Commun. 6, 6991 (2015).

7. Liu, F. et al. Photoresponse: highly sensitive detection of polarized light using anisotropic 2D ReS2. Adv. Funct. Mater. 26, 1146–1146 (2016).

8. Kim, Y., B. K., Y. C., J. H. C. & C. L. Direct synthesis of large-area continuous ReS2 films on a flexible glass at low temperature. 2D Mater. 4, 25057 (2017).

9. Chen, K.-S., Balla, I., Luu, N. S. & Hersam, M. C. Emerging opportunities for two-dimensional materials in lithium-ion batteries. ACS Energy Lett. 2, 2026–2034 (2017).

10. Shi, L. & Zhao, T. Recent advances in inorganic 2D materials and their applications in lithium and sodium batteries. J. Mater. Chem. A 5, 3735–3758 (2017).

11. Yoo, E. et al. Large reversible Li storage of graphene nanosheet families for use in rechargeable lithium ion batteries. Nano. Lett. 8, 2277–2282 (2008).

12. Chang, K. & Chen, W. L-cysteine-assisted synthesis of layered MoS2/graphene composites with excellent electrochemical performances for lithium ion batteries. ACS Nano 5, 4720–4728 (2011).

13. Hwang, H., Kim, H. & Cho, J. MoS2 nanoplates consisting of disordered graphene-like layers for high rate lithium battery anode materials. Nano. Lett. 11, 4826–4830 (2011).

14. Liu, Y., Zhao, Y., Jiao, L. & Chen, J. A graphene-like MoS2/graphene nanocomposite as a high performance anode for lithium ion batteries. J. Mater. Chem. A 2, 13109–13115 (2014).

15. Ai, W. et al. Supramolecular polymerization promoted in situ fabrication of nitrogen-doped porous graphene sheets as anode materials for Li-ion batteries. Adv. Energy Mater. 5, 1500559 (2015).

16. Xie, Y. et al. SnS2 nanoplates as stable anodes for sodium ion and lithium ion batteries. Mater. Technol. 31, 646–652 (2016).

17. Zhang, Q. et al. Extremely weak van der Waals coupling in vertical ReS2 nanowalls for high-current-density lithium-ion batteries. Adv. Mater. 28, 2616–2623 (2016).

18. Fan, X., Zheng, W., Kuo, J. & Singh, D. J. Adsorption of single Li and the formation of small Li clusters on graphene for the anode of lithium-ion batteries. ACS Appl. Mater. Interfaces 5, 7793–7797 (2013).

19. Xiao, B., Li, Y., chun, Yu, X., fang & Cheng, J. -bo Penta-graphene: a promising anode material as the Li/Na-ion battery with both extremely high theoretical capacity and fast charge/discharge rate. ACS Appl. Mater. Interfaces 8, 35342–35352 (2016).

20. Li, W., Yang, Y., Zhang, G. & Zhang, Y. W. Ultrafast and directional diffusion of lithium in phosphorene for high-performance lithium-ion battery. Nano. Lett. 15, 1691–1697 (2015).

21. Guo, G. C. et al. First-principles study of phosphorene and graphene heterostructure as anode materials for rechargeable Li batteries. J. Phys. Chem. Lett. 6, 5002–5008 (2015).

22. Jiang, H. R., Lu, Z., Wu, M. C., Ciucci, F. & Zhao, T. S. Borophene: a promising anode material offering high specific capacity and high rate capability for lithium-ion batteries. Nano Energy 23, 97–104 (2016).

23. Mortazavi, B., Dianat, A., Cuniberti, G. & Rabczuk, T. Application of silicene, germanene and stanene for Na or Li ion storage: a theoretical investigation. Electrochim. Acta 213, 865–870 (2016).

24. Wang, D., Liu, L.-M., Zhao, S.-J., Hu, Z.-Y. & Liu, H. Potential application of metal dichalcogenides double-layered heterostructures as anode materials for Li-ion batteries. J. Phys. Chem. C. 120, 4779–4788 (2016).

25. Jing, Y., Zhou, Z., Cabrera, C. R. & Chen, Z. Metallic VS2 monolayer: a promising 2D anode material for lithium ion batteries. J. Phys. Chem. C. 117, 25409–25413 (2013).

26. Tang, Q., Zhou, Z. & Shen, P. Are MXenes promising anode materials for Li ion batteries? Computational studies on electronic properties and Li storage capability of Ti3C2 and Ti3C2X2 (X = F, OH) monolayer. J. Am. Chem. Soc. 134, 16909–16916 (2012).

27. Hu, J., Xu, B., Ouyang, C., Yang, S. A. & Yao, Y. Investigations on V2C and V2CX2(X = F, OH) monolayer as a promising anode material for Li Ion batteries from first-principles calculations. J. Phys. Chem. C. 118, 24274–24281 (2014).

28. Schwingenschlögl, J. Z. & U. P and Si functionalized MXenes for metal-ion battery applications. 2D Mater. 4, 25073 (2017).

29. Bhandavat, R., David, L. & Singh, G. Synthesis of surface-functionalized WS2 nanosheets and performance as Li-ion battery anodes. J. Phys. Chem. Lett. 3, 1523–1530 (2012).

30. Putungan, D. B., Lin, S. H. & Kuo, J. L. Metallic VS2 monolayer polytypes as potential sodium-ion battery anode via ab‒initio random structure searching. ACS Appl. Mater. Interfaces 8, 18754–18762 (2016).

31. Biswas, D. et al. Narrow-band anisotropic electronic structure of ReS2. Phys. Rev. B 96, 085205 (2017).

32. Yu, Z. G., Cai, Y. & Zhang, Y. W. Robust direct bandgap characteristics of one-and two-dimensional ReS2. Sci. Rep. 5, 13783 (2015).

33. Su, Q., Wang, S., Feng, M., Du, G. & Xu, B. Direct studies on the lithium-storage mechanism of molybdenum disulfide. Sci. Rep. 7, 7275 (2017).

34. Zhang, J. et al. NbS2 nanosheets with M/Se (M = Fe, Co, Ni) codopants for Li+ and Na+ storage. ACS Nano 11, 10599–10607 (2017).

35. Needs, C. J. P. & R. J. Ab initio random structure searching. J. Phys. Condens. Matter 23, 53201 (2011).

36. Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 113, 9901–9904 (2000).

37. Henkelman, G. & Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 113, 9978–9985 (2000).

38. Dronskowski, R. & Bloechl, P. E. Crystal orbital Hamilton populations (COHP): energy-resolved visualization of chemical bonding in solids based on density-functional calculations. J. Phys. Chem. 97, 8617–8624 (1993).

39. Deringer, V. L., Tchougréeff, A. L. & Dronskowski, R. Crystal orbital Hamilton population (COHP) analysis as projected from plane-wave basis sets. J. Phys. Chem. A 115, 5461–5466 (2011).

40. Stefan, M., L., D. V., L., T. A. & Richard, D. Analytic projection from plane-wave and PAW wavefunctions and application to chemical-bonding analysis in solids. J. Comput. Chem. 34, 2557–2567 (2013).

41. Maintz, S., Esser, M. & Dronskowski, R. Efficient rotation of local basis functions using real spherical harmonics. Acta Phys. Pol. B 47, 1165–1175 (2016).

42. Luo, M., Xu, Y. E. & Song, Y. X. Tunable electronic properties of MoS2/ReS2 van der Waals heterostructure from first-principles study. Opt. - Int. J. Light Electron Opt. 144, 334–339 (2017).

43. Wang, Y., Lv, J., Zhu, L. & Ma, Y. Crystal structure prediction via particle-swarm optimization. Phys. Rev. B 82, 94116 (2010).

44. Oganov, A. R. & Glass, C. W. Crystal structure prediction using ab‒initio evolutionary techniques: Principles and applications. J. Chem. Phys. 124, 244704 (2006).

45. Pickard, C. J. & Needs, R. J. High-pressure phases of silane. Phys. Rev. Lett. 97, 45504 (2006).

46. Pickard, C. J. & Needs, R. J. Structure of phase III of solid hydrogen. Nat. Phys. 3, 473–476 (2007).

47. Pickard, C. J. & Needs, R. J. Aluminium at terapascal pressures. Nat. Mater. 9, 624–627 (2010).

48. Pickard, C. J. & Needs, R. J. Highly compressed ammonia forms an ionic crystal. Nat. Mater. 7, 775–779 (2008).

49. George, C., Morris, A. J., Modarres, M. H. & De Volder, M. Structural evolution of electrochemically lithiated MoS2 nanosheets and the role of carbon additive in Li-ion batteries. Chem. Mater. 28, 7304–7310 (2016).

50. Henkelman, W. T. and E. S. and G. A grid-based Bader analysis algorithm without lattice bias. J. Phys. Condens. Matter 21, 84204 (2009).

51. Sanville, E., Kenny, S. D., Smith, R. & Henkelman, G. Improved grid-based algorithm for Bader charge allocation. J. Comput. Chem. 28, 899–908 (2007).

52. Yu, M. & Trinkle, D. R. Accurate and efficient algorithm for Bader charge integration. J. Chem. Phys. 134, 64111 (2011).

53. Rao, F. et al. Reducing the stochasticity of crystal nucleation to enable subnanosecond memory writing. Science 358, 1423–1427 (2017).

54. Chevrier, V. L. & Dahn, J. R. First principles model of amorphous silicon lithiation. J. Electrochem. Soc. 156, A454–A458 (2009).

55. Chevrier, V. L. & Dahn, J. R. First principles studies of disordered lithiated silicon. J. Electrochem. Soc. 157, A392–A398 (2010).

56. Chan, M. K. Y., Wolverton, C. & Greeley, P. First principles simulations of the electrochemical lithiation and delithiation of faceted crystalline silicon. J. Am. Chem. Soc. 134, 14362–14374 (2012).

57. Kresse, G. & Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B 47, 558–561 (1993).

58. Kresse, G. & Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal--amorphous-semiconductor transition in germanium. Phys. Rev. B 49, 14251–14269 (1994).

59. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 6, 15–50 (1996).

60. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab‒initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 54, 11169–11186 (1996).

61. Kresse, G. & Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B 59, 1758–1775 (1999).

62. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996).

63. Grimme, S., Antony, J., Ehrlich, S. & Krieg, H. A consistent and accurate ab‒initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132, 154104 (2010).

64. Momma, K. & Izumi, F. VESTA3 for three-dimensional visualization of crystal, volumetric and morphology data. J. Appl. Crystallogr. 44, 1272–1276 (2011).

65. Hacene, M. et al. Accelerating VASP electronic structure calculations using graphic processing units. J. Comput. Chem. 33, 2581–2589 (2012).

66. Hutchinson, M. & Widom, M. VASP on a GPU: application to exact-exchange calculations of the stability of elemental boron. Comput. Phys. Commun. 183, 1422–1426 (2012).

67. Zhao, Z., Marsman, M., Wende, F. & Kim, J. Performance of hybrid MPI/OpenMP VASP on Cray XC 40 based on Intel Knights landing many integrated core architecture CUG Conference proceedings (2017).

68. Stefan, M., L., D. V., L., T. A. & Richard, D. LOBSTER: A tool to extract chemical bonding from plane-wave based DFT. J. Comput. Chem. 37, 1030–1035 (2016).

## Acknowledgements

This work is partially supported by Space Technology Cell, Indian Space Research Organisation (ISRO), under Grant ISTS/EED/SM/392. Authors thank Supercomputer Education and Research Centre (SERC), Indian Institute of Science (IISc) Bangalore, for GPU and Xeon-Phi based computations.

## Author information

Authors

### Contributions

A.K. performed the DFT calculations, implemented hybrid processor-based structure search and delithiation algorithms, and analysed the results. S.M. conceived the problem statement and overall supervised the work. All authors contributed in the writing.

### Corresponding author

Correspondence to Santanu Mahapatra.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

Reprints and Permissions

Kabiraj, A., Mahapatra, S. High-throughput first-principles-calculations based estimation of lithium ion storage in monolayer rhenium disulfide. Commun Chem 1, 81 (2018). https://doi.org/10.1038/s42004-018-0082-3

• Accepted:

• Published:

• DOI: https://doi.org/10.1038/s42004-018-0082-3

• ### Theory of nonvolatile resistive switching in monolayer molybdenum disulfide with passive electrodes

• Sanchali Mitra
• Arnab Kabiraj
• Santanu Mahapatra

npj 2D Materials and Applications (2021)