2D 1T′-MoS2-WS2 van der Waals heterostructure for hydrogen evolution reaction: dispersion-corrected density functional theory calculations

Dispersion-corrected density functional theory calculations were implemented to investigate structural characteristics, as well as the hydrogen evolution reaction (HER) capability of 2D 1T′ phase MoS2-WS2 van der Waals heterostructures. Two van der Waals corrections were utilized in the study, namely DFT-D3 (semi-empirical-based) and vdW-DF2-B86R (ab-initio-based) corrections. Results show that the DFT-D3 correction stabilized the binding of the monolayers consistent with experimental observations, with binding energy per unit cell of -0.54 eV/cell. The Gibbs free energy of hydrogen adsorption ΔGads,H, which is the lone descriptor of HER, were calculated for the two known adsorption sites in the 1T′ phase, termed S1 (sulfur site with elongated bonds, more active for HER) and S2 (sulfur site with compressed bonds, less active for HER). It is revealed that at the van der Waals region, the S1 and S2 sites, acting as a single adsorption site, become active for HER, with significantly lowered value of ΔGads,H at 0.20-0.24 eV. This is linked to the synergistic interaction of the two sites in adsorbing hydrogen. In terms of electronic structure, the enhanced states in the vicinity of the Fermi level for the interacting S1 and S2 sites at the van der Waals region resulted from orbital hybridization among 3p states of the sulfur sites from the inner top and bottom surfaces. The merging of the two sites at the van der Waals region would result to HER efficiency that is expected to be higher by a factor of 2 compared to that on the top and bottom surfaces. This work has showed that 2D heterostructures could be of importance in catalysis, particularly in HER. Furthermore, it is showed that building a 2D heterostructure could be a good alternative to the application of strain in improving HER capability of 1T′ 2D materials without compromising the adsorption properties of other sites.


Introduction
Due to its zero greenhouse emission property, hydrogen is currently the first in line when it comes to clean alternatives to fossil fuels [1][2][3]. Although hydrogen is not naturally available in large quantities, it can be produced via the hydrogen evolution reaction (HER) [3][4][5]. Here, a catalyst having an optimum hydrogen adsorption activity is employed, as both experimental and theoretical studies showed that the strength of hydrogen adsorption on the catalyst is a lone descriptor of HER activity [6,7]. In calculations, this quantity is computed as the Gibbs free energy of adsorption (ΔG), and for HER to be optimum this should be close to the thermoneutral value (ΔG ≈ 0) [8,9]. The noble metal platinum is currently the best catalyst for HER. However its limited availability and high cost makes it impractical for large scale hydrogen production [10][11][12]. As such, it is of immediate research concern to find cheap and abundant alternatives to platinum.
Upon the discovery of the existence of 2D materials, extensive research has been conducted to efficiently and effectively utilize these materials in improving current technologies and improving known active surfaces [13][14][15]. Specifically, 2D transition metal dichalcogenides (TMDCs) became a subject of interest due to their Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. excellent properties that can be utilized in a variety of applications such as optoelectronics, FETs, metal-ion battery electrodes and catalysis [16][17][18][19][20][21][22]. These 2D materials, although not single atom thick, have the chemical structure MX 2 which consists of one transition metal, M, and two chalcogen atoms, X [13,23]. In addition, several structural phases have been identified for 2D TMDCs, the most studied ones being the 1H (trigonal prismatic) 1T (octahedral) and 1T′ (distorted octahedral [23]. In the field of catalysis, it has been established that most TMDCs are potentially effective HER catalysts [5,[24][25][26]. Being cheap and easy to synthesize, TMDCs are gaining profound attention as potential alternatives to noble metals.
Recent studies showed that the 1T′ phase of most TMDCs are superior in HER compared to its 1H or 2H counterparts [27][28][29]. Moreover, it is revealed that the 1T′ phase has an active basal plane, whereas the 1H phase has only the edge sites as its active sites [24]. This is a very important result, since having an active basal plane translates to a very effective and efficient HER capability as there is no need to expose any edge sites. Among the considered TMDCs in [27][28][29], it was found that 1T′ phase MoS 2 and WS 2 are among the best for HER. Recently, there has been a great interest in stacking up 2D materials to form a 2D van der Waals heterostructure [30]. According to [30], there are several advantages of constructing a 2D van der Waals heterostructure. Combining different 2D materials could effectively result to a new material with unusual properties that could be used as platforms in the investigation of various phenomena. For instance, it was mentioned in [30] that a van der Waals heterostructure consisting of graphene and a dielectric crystal gave rise to properties similar to that of a superconducting copper oxide. Such observation serves as a motivation in studying other possible van der Waals heterostuctures, particulary those that could be of potential application in catalysis. In line with this, it is of great importance and of interest to investigate the properties of 1T′ phase MoS 2 and WS 2 as a 2D van der Waals heterostructure for HER.
In dealing with 2D van der Waals heterostructures, it is necessary to include what are known as van der Waals or dispersion corrections within DFT to account for the dispersion interactions between the monolayers. These corrections are typically accurate enough but are not computationally expensive. Two major implementations of these corrections are usually employed in calculations: those where the dispersion corrections are based on interatomic potentials of the form C 6 ·R −6 that is added to the total Kohn-Sham energy, and those where an exchange-correlation functional is formulated in which the dispersion corrections are incorporated through non-local correlation terms One of the most commonly utilized form of the former are the so-called DFT-D corrections [31][32][33][34], in which DFT-D3 is the latest version. For the latter, there are several methods typically used in calculations, such as the vdW-DF [35,36] and the opt methods (optPBE-vdW, optB86b-vdW, among others) [37,38].
In this work, we carried out van der Waals-corrected density functional theory calculations utilizing two different implementations: the DFT-D3 and vdW-DF2-B86R [39] (a revised version of vdW-DF), to probe the structural, energetics and HER capability of 1T′-MoS 2 -WS 2 2D van der Waals heterostructure. It is found that in terms of structural properties, the inter-monolayer distance is the most sensitive to the type of van der Waals corrections. Moreover, energetics wise, it is revealed that conventional PBE-DFT did not predict a stable 2D van der Waals heterostructure for 1T′-MoS 2 and WS 2 , while the DFT-D3-corrected DFT resulted to a stable binding energy consistent with experimental observations. It is revealed that at the van der Waals region of the heterostructure, both the active S1 and non-active S2 sites become a single active site for HER. Such phenomenon is linked to the synergistic interaction between the two sites at the van der Waals region. This important result implies that at the van der Waals region, the HER efficiency is doubled compared to the lone basal plane of single layer 1T′ MoS 2 and WS 2 since the non-active S2 sites are effectively 'turned on' and become active. Our results point to the importance of 2D van der Waals heterostructure, especially in catalysis, and it is expected that different combinations of 2D materials could be utilized in the control and fine tuning of HER activity.

Methodology
DFT calculations were carried out via the Quantum Espresso package [40,41] utilizing the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional within the Generalized Gradient Approximation (GGA) [42]. Projected-augmented wave (PAW) plane wave basis set was used, with an energy cut-off of 80Ry. The system was relaxed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization method. The reciprocal space is sampled with a 5×5×1 Monkhorst-Packk-point mesh, and the repeating slabs are separated by the vacuum space of 20Å. Electronic and ionic convergence was set to 10 −6 a.u. and 10 −3 a.u. respectively.
The 2D hterostructures were modelled by considering three different stackings of 1T′ phase MoS 2 -WS 2 monolayers (figures 1(a)-(c)). We define the A staking as a direct metal-metal and chalcogen-chalcogen alignment, while B stacking is a direct metal-chalcogen alignment. Finally, due to the asymmetry of the 1T′ phase, a third stacking was defined, stacking C, wherein there is an exact alignment of the chalcogens, but partial alignment of the metal atoms. We constructed 2×2 supercells consisting of four Mo atoms, four W atoms, and 16 sulfur atoms for the considered stackings. XCrysDen and Vesta were used in visualization and manipulation of the structures [43,44]. The top layer of the 2D heterostructure consist of the WS 2 monolayer while the bottom layer is composed of the MoS 2 monolayer. As these 2D heterostuctures are primarily held by disperion forces, Van der Waals corrections to the total Kohn-Sham energy were incorporated into the calculations. We utilized two flavors of van der Waals corrections, the semi-empirical DFT-D3 [32] and the ab-initio-based vdW-DF2-B86R [39].
The stability of the 2D heterostructures is effectively described by the binding energy: where E heterostructure is the total energy of the 2D heterostructure, and E MoS 2 and E WS 2 are the energies of the MoS 2 and WS 2 monolayers respectively. The binding energy per unit cell (eV/cell) is obtained by dividing equation (1) with the size of the supercell. Adsorption sites for hydrogen adsorption were identified consistent with a previous work utilizing 1T′-MoS 2 [27]. As the 1T′ phase possesses a level of asymmetry due to the dimerization of the metal atoms, the two adsorption sites consisting of a sulfur site with elongated bond (S1) and a sulfur site with compressed bond (S2) were chosen. Hydrogen atom was placed initially at 1.5 Å from the adsorption site and the system was optimized. Hydrogen adsorption energy (ΔE ads, H ) was calculated using where E substrate+H is the total energy of the surface with the adsorbed hydrogen, E substrate is the energy of the pristine 2D heterostructure and H 2 is the energy of a hydrogen molecule. From the adsorption energy the Gibbs free energy of adsorption, ΔG ads,H , was calculated via where DE ZPE is the difference in zero-point energy of hydrogen in the adsorbed state and the gas phase, and DS H is the entropy difference between the adsorbed state and the gas phase of hydrogen. This term is usually approximated as D » S S H H 1 2 [26,[45][46][47]. Here, S H is the entropy of the gas phase H 2 at standard conditions (300 K). The zero-point energy corrections for adsorbed hydrogen were obtained from the frequencies of hydrogen vibrational modes as calculated via the density functional perturbation theory (DFPT) [48]. As mentioned earlier, for optimum HER performance, ΔG ads,H must be close to zero and thus the lower ΔG ads,H is, the greater is the potential for HER.

Results and discussions
3.1. Structural properties and energetics of 1T′-MoS 2 -WS 2 2D heterostructure We first discuss the physical structure of the 2D heterostructures by calculating lattice parameters with reference to the model structure shown in figure 2. Table 1 shows the obtained lattice constants, lattice vectors and interlayer distance for plain DFT, DFT with D3 correction, and DFT with vdW-DF2-B86R corrections. Due to the  In general, there are no significant differences in the lattice parameters of the vdW heterostructures, except for the inter-planar distance between the constituent monolayers. This is expected as the interaction between the monolayers is through van der Waals forces. It is noted that addition of van der Waals corrections tends to reduce inter-layer distance, with the vdW-DF2-B86R correction giving the smallest inter-planer distance among the considered correction schemes. In terms of stacking, it can be observed that the B stacking has the least variation in inter-layer spacing, whereas A stacking has the most dramatic decrease in inter-layer distance (from 6.84 Å for no correction to 6.05 Å for the vdW-DF2-B86R correction). As for the heterostructures' stability, the binding energies per monolayer calculated for plain DFT, DFT with D3 correction and DFT with vdW-DF2-B86R correction are presented in table 2. As 2D MoS 2 -WS 2 heterostructures have been reported in experimental studies [50][51][52][53], it is expected that MoS 2 and WS 2 would bind and the interaction between the two monolayers is attractive.
It is interesting to note that using plain DFT resulted to a binding energy of effectively 0eV, implying that MoS 2 and WS 2 have no interaction and the 2D heterostructure is unstable. Surprisingly, the addition of vdW-DF2-B86R correction resulted to a more unstable heterostructure, with a very large positive binding energy. These results are inconsistent with experimental and other theoretical calculation results [50][51][52][53][54] that reported the existence of stable MoS 2 -WS 2 heterostructure. The only van der Waals correction scheme where a stable binding was obtained was found to be the DFT-D3. Thus, it can be inferred that for this particular 1T′-MoS 2 -WS 2 2D heterostructure, the DFT-D3 van der Waals correction appears to be the best correction that captures the inherent binding between the two monolayers. All of the considered three stacking configurations are stable under DFT-D3, with effectively equal binding energies.
We proceeded with hydrogen adsorption calculation, and subsequently ΔG ads,H by utilizing the B stacking configuration. The B stacking was selected due to it being one the most stable stackings energetically. Moreover, previous experimental studies suggest that for 1T′-MoS 2 -WS 2 2D heterostucture, the B stacking appears to be the most stable [53].

ΔG ads,H on B stacking 1T′-MoS 2 -WS 2 heterostructure
HER is a two-step reaction, which can be written in general as [55] The initial step of the reaction is the reduction of H + and its subsequent adsorption on to a site (H + + e − + * →H * , H * implies adsorbed hydrogen), followed by the combination of two adsorbed hydrogen atoms to form hydrogen gas (2H * →H 2 ) or the capture of H + and an electron by H * to form hydrogen gas. Here, the ΔG ads,H can be interpreted as a kind of energy barrier needed to cross the reaction from an adsorbed hydrogen state to a free hydrogen gas molecule (see figure 4).
The S1(longer bond) and S2 (shorter bond) adsorption sites are produced from the distorted symmetry of the 1T′ phase TMDC This is expected to induce asymmetric adsorption between the two sites [27,56]. Due to Table 2. Binding energies per cell of different 2D 1T′-MoS 2 -WS 2 heterostructure stackings. The suffix -GD3 implies binding energies with the Grimme D3 vdW corrections, while -rev denotes binding energies with the vdW-DF2-B86R corrections.

Stacking
Binding Energy per unit cell (eV/cell) C-rev unstable this asymmetry, hydrogen adsorption properties on the S1 and S2 sites are different, and it was showed that hydrogen adsorption is much more stable on the S1 site compared to the S2 site. This translates to a lower ΔG ads,H on the S1 site compared to that on the S2 site, which implies better HER capability. This scenario makes the basal plane of a single monolayer 1T′-MoS 2 and 1T′-WS 2 only 50% effective for HER. Application of biaxial strain helps in 'turning on' the S2 site and reduce its ΔG ads,H to significantly lower values [27]. However, one drawback of such process is that too large biaxial strain could make the S1 site's affinity to hydrogen too strong, making its ΔG ads,H shift away from the desired thermoneutral value. As such, it is of great interest to look into the potential advantages of combining these two monolayers to form a 2D van der Waals heterostructure, especially at the interface where the two monolayers meet (referred to as the van der Waals region). Figure 3 shows the optimized structures of adsorbed hydrogen on the top sites (WS 2 ), bottom sites (MoS 2 ) and on sites at the van der Waals region. For the considered B stacking 1T′-MoS 2 -WS 2 2D heterostructure, it was found that for the top sites (WS 2 ) and bottom sites (MoS 2 ), the S1 sites exhibited significantly lower ΔG ads,H compared to that on the S2 sites, consistent with the findings of [27,56] (table 3 and figure 4). At the van der Waals region, the calculated ΔG ads,H for both the S1 and S2 sites are comparable to that on the bottom S1 sites (c) On an S1 site located at the inner center surface (van der Waals region), WS 2 side; (d) On an S1 site (initially adsorbed on an S2 site, WS 2 side) located at the inner center surface (van der Waals region), MoS 2 side; (e) On an S1 site located at the top surface, WS 2 side; and (f) On an S2 site located at the top surface, WS 2 side.
(MoS 2 side), and are surprisingly lower than that on the top S1 sites (WS 2 side) of the 2D heterostructure. Moreover, within the van der Waals region, it is found that in the case of hydrogen adsorption on an S2 site, the hydrogen atom tends to be attracted to an opposite, nearest-neighbor S1 site. As such, a hydrogen atom initially adsorbed on an S2 site coming from the WS 2 side migrates to an S1 site on the MoS 2 side, and vice versa. This makes sense as the S1 site is known to have a higher affinity to hydrogen compared to an S2 site [27,56]. Further investigation revealed that hydrogen appears to be bonded to both sites, with bond lengths 1.36Åfrom the S1 site and 2.60Åfrom the S2 site. Surprisingly, these bond lengths are essentially the same for hydrogen adsorption on an S1 site where the hydrogen-S1 bond length is 1.36Åand that for hydrogen-S2 is 2.42Å. These imply that at the van der Waals region, hydrogen adsorption is mediated by an S1 site on one side and an S2 site on the other side of the van der Waals region. And even though the two monolayers comprising the 2D Figure 4. Graphical representation of the calculated ΔG ads,H with the thermoneutral value (0 eV) as the reference point for the different adsorption sites of the 2D 1T′-MoS 2 -WS 2 heterostructure. Note that the ΔG ads,H can be interpreted as an energy barrier that is needed to cross the reaction path from hydrogen reduction, passing through the adsorbed hydrogen state (denoted by H * ), to a free hydrogen atom that will subsequently form a H 2 gas molecule. heterostructure are different, the calculated ΔG ads,H at the van der Waals region are effectively the same for both S1 an S2 sites, as these two sites (one coming from MoS 2 side and the other from WS 2 side) synergistically adsorb hydrogen, making ΔG ads,H comparable to the bottom MoS 2 and much lower compared to top WS 2 S1 sites. The merging of the S1 and S2 sites at the van der Waals region is effectively equivalent to the 'turning on' of the S2 site for HER, similar to the effect of applying biaxial strain. However, one important advantage of this synergistic interaction between the S1 and S2 sites at the van der Waals region is that the HER ability of the S1 site is not compromised with the 'turning on' of the S2 site. Both the S1 and S2 sites, as a synergistically combined adsorption site, is converted to a single adsorption site that is effective for HER. And since there are effectively twice the number of merged S1 and S2 sites at the van der Waals region compared to that on the top and bottom surfaces of the 2D heterostructure, higher hydrogen coverages is expected here and thereby HER activity is expected to be higher by a factor of 2.
To make further inferences from these calculated ΔG ads,H , we plotted the projected density of states (PDOS) of the S1 and S2 sites at the top, van der Waals region (center) and bottom surfaces of the 2D heterostructure in figures 5(a)-(f). It has been established previously [27,56] that hydrogen's affinity to a sulfur site on 1T′ TMDCs is related to the availability of electronic states in the vicinity of the Fermi level, and that in the case of the 1T′ phase S1 sites typically possess more states near the level compared to S2 sites. In this regard, the abundance of states near the Fermi level is a signature property, especially of 2D 1T′ TMDs, that translates to a strong hydrogen atom adsorption. Our results are consistent with these findings. In figures 5(a) and (b), for the bottom MoS 2 surface layer, it is observable that there are more available sulfur 3p electronic states near the Fermi level for the S1 site compared to that for the S2 site. The same is true for the top WS 2 surface layer (figures 5(e) and (f). Since the bottom and top surface layers of the 2D heterostructures are effectively free layers with no interaction with other surfaces, the observed consistency in the ΔG ads,H trend is justified, wherein S2 sites for both surface layers exhibited much higher ΔG ads,H values compared to S1 sites. At the van der Waals region of the 2D heterostructure, the S1 and S2 are essentially interacting, and as such, the observed trend of ΔG ads,H for S1 and S2 differ compared to that on the bottom and top surface layers. As such, it can be inferred that at the van der Waals region, S1 and S2 sites 3p states hybridized, modifying the density of states in the vicinity of the Fermi level. Such interaction effectively created a new site composed of S1 and S2 with PDOS that is a superposition of S1 and S2 3p states (figures 5(c) and (d)). Furthermore, the interaction between the S1 and S2 sites of MoS 2 and WS 2 in the van der Waals region resulted to enhanced hybridization of 3p electronic states at the Fermi level. This hybridization between the 3p states of the S1 and S2 sites at the van der Waals region is essentially a reinforcing mechanism that maintains the abundance of states near the Fermi level. As such, at the van der Waals region, the merging of the S1 and S2 sites resulted to a new site that has hydrogen adsorption properties comparable to that of a lone S1 site and possesses very good potential for HER. Furthermore, as mentioned earlier, there are twice the number of merged S1 and S2 sites at the van der Waals region making HER at the region twice as efficient compared to that of the top and bottom surfaces. With these, it is inferred that van der Waals interaction could alter the adsorption properties of specific sites, which could lead to fine tuning of certain catalytic processes like HER by controlling the van der Waals interaction between the two monolayers.

Conclusions
Density functional theory calculations were carried out to probe the structural properties of 2D 1T′-MoS 2 -WS 2 vdW heterostructure, as well as the electronic structure of its sulfur sites to gauge their HER potential. It was found that for all the considered bilayer stackings, vdW corrections did not appreciably affect their structural properties except for the monolayer-monolayer distance within the vdW region. Energetics wise, the addition of vdW corrections to the Kohn-Sham total energy, in particular the DFT-D3 correction, resulted to good binding that is consistent with experimental findings. In terms of hydrogen adsorption, it was confirmed that for both the WS 2 top surface and MoS 2 bottom surface, the S1 site is the more active site compared to S2, with much lower ΔG ads,H that is consistent with previous results. Interestingly, it was revealed that within the vdW region of the 2D heterostructure, both S1 and S2 sites exhibited low ΔG ads,H resulting from the synergistic mediation of the two sites in adsorbing hydrogen. This is turn, resulted to a merged S1 and S2 site which has twice the efficiency of the top and bottom surfaces of the 2D heterostructure. Density of states plots showed that within the vdW region, the hybridization of the 3p states of both the S1 and S2 sites from the different monolayers (MoS 2 and WS 2 ) gave rise to enhanced density of states within the vicinity of the Fermi level, which is an indicator of good HER performance. This work has demonstrated the potential applications of 2D heterostructures in catalysis, particularly in HER. It is expected that other combinations of 2D materials would result to new HER capabilities, and HER fine tuning and control could be achieved through their utilization. Furthermore, building a 2D heterostructure could be one of the better alternatives to the application of strain especially for 1T′ 2D materials as there are no adsorption sites that are compromised.