Next Article in Journal
Multifunctional Hollow Porous Fe3O4@N-C Nanocomposites as Anodes of Lithium-Ion Battery, Adsorbents and Surface-Enhanced Raman Scattering Substrates
Next Article in Special Issue
Theoretical Insights on ORR Activity of Sn-N-C Single-Atom Catalysts
Previous Article in Journal
Antioxidant Activities, Anticancer Activity, Physico-Chemistry Characteristics, and Acute Toxicity of Alginate/Lignin Polymer
Previous Article in Special Issue
Experimental and Simulation Studies on Hematite Interaction with Na-Metasilicate Pentahydrate
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations

1
Department of Theoretical, Atomic and Optical Physics, University of Valladolid, 47011 Valladolid, Spain
2
The Abdus Salam International Centre for Theoretical Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy
*
Author to whom correspondence should be addressed.
Molecules 2023, 28(13), 5182; https://doi.org/10.3390/molecules28135182
Submission received: 5 June 2023 / Revised: 23 June 2023 / Accepted: 28 June 2023 / Published: 3 July 2023
(This article belongs to the Special Issue Computational and Theoretical Chemistry for Material Research)

Abstract

:
Density functional theory (DFT) calculations are employed to study the oxygen evolution reaction (OER) on the edges of stripes of monolayer molybdenum disulfide. Experimentally, this material has been shown to evolve oxygen, albeit with low efficiency. Previous DFT studies have traced this low catalytic performance to the unfavourable adsorption energies of some reaction intermediates on the MoS2 edge sites. In this work, we study the effects of the aqueous liquid surrounding the active sites. A computational approach is used, where the solvent is modeled as a continuous medium providing a dielectric embedding of the catalyst and the reaction intermediates. A description at this level of theory can have a profound impact on the studied reactions: the calculated overpotential for the OER is lowered from 1.15 eV to 0.77 eV. It is shown that such variations in the reaction energetics are linked to the polar nature of the adsorbed intermediates, which leads to changes in the calculated electronic charge density when surrounded by water. These results underline the necessity to computationally account for solvation effects, especially in aqueous environments and when highly polar intermediates are present.

1. Introduction

In the face of the urgent global need for sustainable and clean energy sources, solar and wind energy have emerged as viable alternatives to traditional fossil fuels. These renewable resources are not only abundant but also capable of significantly reducing the carbon footprint associated with energy production. However, a major challenge associated with these intermittent energy sources is their inability to provide a consistent supply of power in the absence of sunlight or wind. Hence, it is crucial to have effective energy storage solutions to harness the full potential of these renewable energy sources.
One promising approach to address this intermittency issue involves converting surplus energy from solar and wind sources into hydrogen via water splitting [1]. Hydrogen is a high-energy, non-polluting fuel that can be stored and transported with relative ease, offering a potential solution to energy storage and delivery [2]. Nevertheless, the water-splitting process, which involves separating water into hydrogen and oxygen, requires the presence of an efficient catalyst to be energetically feasible.
In the pursuit of viable catalysts for this purpose, a key focus has been the use of earth-abundant materials. Such materials are not only widely available, but they also circumvent the problems associated with the cost, rarity, and potential toxicity of more traditionally used precious-metal-based catalysts. As such, research has intensified to identify earth-abundant materials that can act as efficient catalysts for water-splitting.
Layered Transition-Metal Dichalcogenides (TMDs) have shown significant promise as catalysts in the water-splitting process [3,4,5]. These compounds, which are characterized by layers of transition metal atoms sandwiched between chalcogen atoms [6], have demonstrated unique catalytic properties, including favorable reaction kinetics and long-term stability. Moreover, their layered structure provides ample active sites for the reaction, which may further enhance their catalytic efficiency [7,8].
Among the layered TMDs, Molybdenum Disulfide (MoS2) is a prominent representative, which is attracting considerable attention in scientific research due to its exceptional properties and potential applications [3].
Apart from its emerging role as a catalyst in water-splitting, MoS2 has found many other applications in various fields. For instance, owing to its semi-conductive nature and impressive electron mobility, MoS2 has been utilized in the fabrication of transistors [9,10,11], making it a potential contender for next-generation electronics. Moreover, its unique photoelectric properties make it an ideal candidate for photovoltaic applications and optoelectronic devices [12]. The compound has also been used in the realm of energy storage, serving as an electrode material in lithium-ion batteries due to its high capacity for lithium intercalation [13]. Furthermore, the high surface area and strong adsorption capacity of MoS2 have led to its application in environmental fields, such as water- and air-purification technologies [14,15].
Recently, there has been particular interest in towards the catalytic properties of MoS2 monolayers, especially regarding the Oxygen Evolution Reaction (OER), a critical half-reaction in water-splitting. It is well-established that the edge sites of MoS2 are highly catalytically active [16,17]. For example, the HER is known to be catalysed by edges, not by the inert basal plane of MoS2 [18]. For the case of the OER, we have recently shown [19] that the pristine monolayer exhibits a strong over-binding of the O * reaction intermediate, which leads to an unfavourable alignment of the energy levels [19]. Bulk MoS2 provides limited edge sites, restricting its overall catalytic performance. A monolayer of MoS2, with its two-dimensional structure, presents a significantly higher proportion of edge sites compared to its bulk counterpart. This maximization of active sites potentially enhances its catalytic efficiency, making monolayer MoS2 a promising catalyst for the OER. Indeed, layers of MoS2 and MoS2-WS2 heterojunctions have been experimentally observed to split water [20].
In the field of catalysis, theoretical and computational strategies are being increasingly recognized for their instrumental role in supplementing experimental observations [21]. These approaches not only provide a computational framework to validate empirical findings but also offer an interpretative lens to discern underlying catalytic mechanisms that are otherwise challenging to investigate experimentally. In this regard, Density Functional Theory (DFT) frequently emerges as the methodology of choice, due to its potency in calculating electronic structures and predicting chemical properties [22,23]. In the particular case of the OER, the key processes to understand are proton-coupled electron transfers (PCETs). The production of one molecule of O2 from water involves four distinct PCET steps. Even though these processes happen in a liquid environment, it is common practice to perform the corresponding DFT calculations in a vacuum, i.e., in a setup where the reaction intermediates are bonded to the catalyst surface, but the (water) environment surrounding the catalyst and reactants is neglected.
Approximate methods to account for the presence of a liquid environment have existed for many years. They are typically focused on mimicking the effects of the dipolar nature of polar solvents such as water. Most commonly, one resorts to so-called continuum solvation models where a dielectric is defined as a smooth self-consistent function of the electronic density of the solute. Such methods have very successfully been applied in the realm of theoretical chemistry, while applications in condensed matter or solid-state physics are more rare.
In this paper, we apply the revised self-consistent continuum solvation (SCCS) methodology to the edges of monolayer MoS2. We show how the strong dipoles present in typical OER reaction intermediates lead to important corrections of the adsorption energies and hence to a non-negligible correction of the estimated overpotential for the OER in this system.

2. Methodologies

2.1. The Implicit Solvent Model

Implicit solvent models, also known as continuum models, have been an integral part of computational chemistry and molecular simulations for several decades [24,25,26]. These models aim to describe the solvation environment around a molecule by approximating the solvent as a continuous medium, thereby bypassing the need to explicitly represent each solvent molecule. This simplification significantly reduces the computational demand and enables the study of larger systems or longer timescales than can typically be achieved with explicit solvent models. Many different continuum solvation models have been developed, predominantly within the realm of chemistry, since Onsager’s seminal contributions [27]. Among the multitude of such methodologies, one of the most well-received is the Polarizable Continuum Model (PCM) advanced by Tomasi and his colleagues [25,26]. The latest iteration of PCM incorporates a broad spectrum of similar techniques, such as the Conductor-like Screening Model (COSMO) approach [24]. Despite its popularity, PCM’s utilization has largely been confined to the chemistry discipline, and it remains relatively unexplored within the spheres of condensed matter and solid-state simulations. Specifically, the capacity to handle metallic systems within the PCM framework was subsequently introduced in an implicit manner [28]. Moreover, the corresponding algorithm was not inherently structured to interface with periodic systems and ab initio molecular dynamics (MD) simulations, further limiting its application in these domains.
Here, we employ an implicit solvation methodology that was first proposed by Fattebert and Gygi [29,30] and is well-suited to plane-wave pseudopotential calculations in periodic systems. In this approach, the dielectric function is defined as a smooth self-consistent function of the electronic density of the solute. A modified version of this procedure was suggested by Andreussi, Darbi and Marzari [31]. It overcomes several of the limitations of previous approaches. Their methodology is publicly available in the quantum-ENVIRON software, which is interfaced with, among others, the  quantum-ESPRESSO suite of programs [32].
This methodology, which is also employed here, involves defining an interface function, s ( r ) , with respect to certain degrees of freedom within the system. The design of this interface function is such that it demonstrates a smooth transition from a value of 1 within the volume where the system’s degrees of freedom exist to a value of 0 in the region of embedding. This progression ensures a comprehensive representation of the embedded system while also capturing the nuances of its interaction with the embedding environment. In the case of SCCS, which are used here, the interface function is defined in terms of the solute’s electron density ρ ( r ) in the following way:
s ( r ) = 0 ρ ( r ) ρ m i n 1 t log ( ρ m a x ) log ( ρ ( r ) ) log ( ρ m a x ) log ( ρ m i n ) ρ m i n < ρ ( r ) < ρ m a x 1 ρ ( r ) ρ m a x ,
where t ( x ) is a smooth step function going from 0 to 1 and is defined as t ( x ) = x sin 2 π x / 2 π .
The limiting values ρ m i n and ρ m a x determine how close or distant the surface s ( r ) is from the solute. In terms of the interface function s ( r ) , the dielectric function ϵ ( r ) is defined as
ϵ ( r ) = exp log ϵ 0 · 1 s ( r ) .
Once the dielectric function is defined as a function of s ( r ) , the electrostatic potential ϕ ( r ) of the system in the presence of the continuum solvent is given by a generalized Poisson equation:
· ϵ ( r ) ϕ ( r ) = 4 π ρ ( r ) .
The modified electrostatic potential ϕ ( r ) stabilizes the electronic densities of localized charges, and thus provides an electrostatic stabilization, especially in charged or dipolar systems like those present in the OER. The implicit solvent, which represents the effect of polar water molecules, can, in this way, play an important role in the stabilization of various reaction intermediates.

2.2. Model System of the MoS2 Edge

Molybdenum disulfide (MoS2), analogously to most TMDs from groups 4–7, exhibits a graphite-analogous layered lattice structure. Each stratified unit of this structure comprises a triadic arrangement of atom planes: a central plane of molybdenum (Mo) atoms lies between two exterior planes of sulfur (S) atoms. The spatial separation of these atomic planes is approximately 3.14 Å.
Intralayer bonding within this structure predominantly involves covalent interactions between the Mo and S atoms in a hexagonal lattice. This lattice exhibits a trigonal prismatic ( D 3 h point group) configuration, frequently designated as 2H [33,34].
In macroscopic manifestations of MoS2, these layers are bonded primarily through weak, non-directional van der Waals (vdW) forces. This form of interaction underlies the facile cleavage of the bulk material into individual monolayers. Figure 1a,b graphically depicts such a representative monolayer of MoS2 [35,36].
As we are interested in the catalytic properties of a single monolayer, we model this within periodic boundary conditions by including about 15 Å  of vacuum between repeated images of the layer. This distance is enough to effectively decouple adjacent layers.
The edge of the layer is represented by an infinite stripe model. A  4 × 4 MoS2 supercell is periodically repeated in one direction only, which leads to a stripe that exposes two inequivalent edges; see Figure 1c. One edge, known as the sulfur edge, exposes the outermost sulfur atoms of the hexagonal lattice. As it is the case inside the MoS2 layer, there are two sufur atoms at the edge, one above and one below the Mo layer. Each S-atom is bonded to two of the outermost Mo atoms. The sulfur edge is therefore obtained from a simple cut of the pristine MoS2 layer.
The second kind of edge is called molybdenum edge. It is created by a cut from the pristine layer with Mo atoms at the edge, followed by decorating the edge with one sulfur atom per Mo atom. In the minimum energy configuration of this edge, the terminating S atoms form pairwise bonds with one of their edge neighbours, with a long bonding distance of 2.19 Å. The molybdenum edges are thus terminated by pairs of weakly bound sulfur atoms, as shown in Figure 1c.
Such models for the edges of MoS2 layers have frequently been studied in the context of catalysis. In Ref. [37], the stability and structures of the considered Mo and S edges are examined by DFT. In Ref. [38], DFT calculations and experiments are used to study the role of transition metal atoms on the catalytic activities of the edges. Very recently, a joint experimental and computational study of MoS2 edges highlighted the role of noble metals at such edges [39].

2.3. Water Splitting and PCETs

The splitting of water into oxygen and hydrogen,
2 H 2 O 2 H 2 + O 2 ,
can be written as two separate half-reactions:
2 H + + 2 e H 2 ,
known as the hydrogen evolution reaction (HER) and
2 H 2 O 4 H + + 4 e + O 2 ,
which is the oxygen evolution reaction (OER). As  is clear from Equation (6), the OER involves four distinct deprotonations, each one coupled with the transfer of an electron, i.e., four PCETs. In order to assess the energy needed for the OER in Equation (6), one should have access to the (free) energy of a proton and an electron, as they appear on the right-hand side of Equation (6). It is very difficult to evaluate those two energies separately, but one can use the fact that the relevant quantity is the sum of both energies. Using the fact that the free energy of (H+ + e) is equal to half the energy of a hydrogen molecule under standard conditions, it is possible to avoid the calculation of a proton’s or an electron’s energy by replacing their sum with 1 2 E ( H 2 ) . This is known as the “computational hydrogen electrode” (CHE), and was first used by Nørskow and collaborators [40].
The CHE only provides the energy difference between subsequent reaction intermediates, with one difference for each PCET. It does not provide kinetic information, for which the barriers between the initial and final state would also need to be known.
One particularly important barrier is the formation of the O−O bond. This formation is non-trivial because of the spin polarization of the oxygen molecule [41]. As a direct simulation of this step goes beyond the CHE approach used here; we point the interested reader to a recent review about modeling of elementary steps in electrocatalysis [42]. Work along the lines described in Refs. [42,43] is currently under way to study the O−O bond formation at MoS2 edges.
Here, we consider the following four PCETs, which are established in the OER:
2 H 2 O + ( ) * A H O * + H 2 O + ( H + + e )
                                      B O * + H 2 O + 2 ( H + + e )
                          C H O O * + 3 ( H + + e )
                                      D ( ) * + O 2 + 4 ( H + + e ) ,
where the *-symbol indicates binding to the catalyst. From steps A to D above, it can be seen that, apart from the pristine MoS2 edge, three more systems with reaction intermediates bound to the edge are needed, namely O * , HO * , and O*HO * . The MoS2 infinite stripe with the reaction intermediates is shown in Figure 2a–c in their minimum energy configuration.
In order to evaluate the energetics of the four PCET steps above, the relevant quantity is the Gibbs free energy. Apart from the internal energy, calculated as the DFT total energy, also the pressure-volume work and entropic contributions are part of the Gibbs free energy. By far the largest corrections to the internal energy are linked to the free energy of the gas-phase O2 and H2 molecules. In this work, we follow the usual practice [44] to neglect entropic contributions to the reaction intermediates, which are bound to the catalyst, and consider the gas-phase molecules at standard conditions.

2.4. Density Functional Theory Calculations

Using spin-polarized Density Functional Theory (DFT), all computations are carried out with the quantum-ESPRESSO software, a comprehensive suite designed for the exploration and post-processing of electronic structures and molecular and condensed systems. The exchange and correlation functional is approximated with the semi-local gradient-corrected form, as parametrized by Perdew, Burke, and Ernzerhof (PBE). We use a plane-wave basis set throughout. The kinetic energy cutoffs for the orbitals and the charge density have been set to 45 Ry and 360 Ry, respectively. It was carefully verified that these cutoffs lead to well-converged energies and forces. The interaction between electrons and ions is modeled by employing Vanderbilt’s ultrasoft pseudopotentials (USPP). Concretely, we use the pseudopotentials Mo_ONCV_PBE-1.0.oncvpsp.upf (valence charge Z v a l = 14 ), s_pbe_v1.4.uspp.F.UPF ( Z v a l = 6 ), O.pbe-n-kjpaw_psl.0.1.UPF ( Z v a l = 6 ), and H.pbe-rrkjus_psl.1.0.0.UPF ( Z v a l = 1 ) for molybdenum, sulfur, oxygen, and hydrogen, respectively. These pseudopotentials can be found in the Standard Solid-State Pseudopotentials (SSSP) database [45,46].
The Brillouin zone (BZ) of this quasi-one-dimensional edge system is sampled using a 4 × 1 × 1 regular grid of k-points with a shift [47]. For the occupation numbers of the KS orbitals, we use the “cold” smearing method of Marzari [48] with an electronic temperature of 136 meV.
In order to account for the weak dispersion interactions, which can play an impotant role in weak binding regimes, we also include the correction known as Grimme-D3 [49].

3. Results and Discussion

3.1. Electronic and Magnetic Structure

In a first step, we consider the density of states (DOS) for the MoS2 stripes with and without the reaction intermediates adsorbed. In Figure 3a, we depict the pristine stripe. It is characterized by a band gap of 0.28 eV. The valence band edge is dominated by majority-spin (spin-up) states that are mainly localized on sulfur atoms. The conduction band edge, in contrast, is dominated by Mo states. Clearly, the main exchange splitting in the spin channels can be observed for the sufur states. This is in agreement with the observation of the main magnetic polarization, which is localised on S atoms at the sufur edges of the stripe.
When the reaction intermediates are bonded to the stripe, it loses its semi-conducting nature. In all cases, the majority spin channel is characterized by a partially occupied sulfur-state close to the Fermi level. In the cases of the HO * and O*HO * adsorbates, the minority spin channel features a partially occupied state, mostly of Mo character. The case of the O * adsorbate is semi-metallic, with a band gap of 0.4 eV in the minority spin channel.
Interestingly, the adsorbate atoms themselves, i.e., the oxygen atoms, contribute to the DOS close to the Fermi energy only in the cases of the HO * and O*HO * adsorbates, but not when a single oxygen atom is adsorbed on the stripe.
In none of the considered cases is any magnetic moment localized on the adsorbed atoms of the reaction intermediates. Instead, the magnetization is localized at the sulfur edge of the stripe.

3.2. Energetics of the OER

In order to evaluate the effect of the liquid solvent on the OER, we consider the particular case of the edge of a MoS2 stripe as catalytic site. Using the common practice of considering catalyst and adsorbed reaction intermediates in a vacuum environment, it has previously been shown that the molybdenum edge is the more advantageous site for the OER on MoS2, with an overpotential of 1.15 eV [19]. In the following, we will study the OER at such Mo-edges when the water environment is taken into account and compare the key results with the standard computations in vacuum.
The geometry of Mo-edges in the minimum energy configuration is unaltered by the inclusion of an implicit solvent model: weakly bonded pairs of sulfur atoms decorate the molybdenum atoms at the edge, as shown in Figure 1c above. When it comes to adsorbing the three reaction intermediates, their bonding geometries are shown in Figure 2, and the corresponding bond lengths are summarized in Table 1. The general picture of the adsorption remains unaltered by the inclusion of the solvent: one of the two pairs of edging sulfur atoms is broken, and the intermediates bind to one (for the case of O * and HO * ) or both (for the case of O*HO * ) of these sulfur atoms.
Looking at Table 1, it is notable that the bonding distances are little affected by the presence of the water solvent. In fact, the bond lengths between the edging sulfur and the O * adsorbate are essentially the same in both cases. The same is true for the O−H bond, both in the HO * adsorbate and in the O*HO * cases, where the bond length varies by less than 1 percent. The only significant change in bonding is observed for the bond length of the edging sulfur atom with the oxygen of the HO * group, both for the case of the HO * and the O*HO * intermediates. In this case, the inclusion of the water solvent results in a reduction in the bond length of about 2.5%.
In order to understand the influence of the water solvent on the bonding of the hydroxyl group, we plot, in Figure 4, the difference in charge density ρ w a t e r ρ v a c u u m for the case of the HO * intermediate, where ρ w a t e r is the electronic charge density obtained with the implicit solvent model and ρ v a c u u m the charge density without any solvent. Both charge densities were calculated in the same geometry, namely that of the system relaxed with water solvent. Clearly, the presence of the (polar) solvent leads to a reduction in the electronic charge around the proton and an increase around the oxygen adatom, with the shape of a d-orbital. It is this population of the polarization orbital on O that goes along with the reduction in the bonding distance to the edging sulfur atom. In Figure 4, the only places at which where there is a significant change in electronic density due to the solvent are around the adsorbed hydroxyl group.
Considering only the adsorption geometry, the influence of the solvent thus seems negligible. This picture changes when the adsorption energies are considered. We define the adsorption energies E a d s of the various intermediates in the following way:
E a d s = E M o S 2 + intermediate E M o S 2 N O μ O N H μ H ,
where E M o S 2 + intermediate is the DFT energy of the MoS2 edge system with the reaction intermediate, as shown in Figure 2a–c. E M o S 2 is the DFT energy of the pristine MoS2 stripe of Figure 1c. N O and N H are the numbers of oxygen and hydrogen atoms in the considered intermediate, respectively. The chemical potentials μ O and μ H of oxygen and hydrogen are obtained as half the energies of O2 and H2 molecules. For the systems considered in water or vacuum, all energies and chemical potentials appearing in Equation (11) are calculated with or without implicit solvent contributions, as detailed in Section 2.1.
The resulting adsorption energies are summarized in Table 2. The presence of the water environment strongly reduces the adsorption strength of all considered reaction intermediates. The strongest reduction is observed for the O*HO * intermediate, where the solvent leads to a reduction in the adsorption energy of 0.37 eV.
Linked to the adsorption energies of the intermediates are the energies involved in the four PCET steps A–D from Equations (7)–(10). The energy required for each step is defined as the energy difference between the products and reactants. In Table 3, we report the height of each one of these steps, both for the case of vacuum (from Ref. [19]) and water solvation (this work). The over-potential of the overall OER reactions is given by the difference between the highest step (step D in both cases considered here) and the ideal step height of 1.23 eV.
The data from Ref. [19] in Table 3 show that the molybdenum edge of MoS2 stripes has an overpotential of 1.15 eV for the OER, which seriously reduces the efficiency of this catalyst. This overpotential can be traced back to the too-strong binding of all the reaction intermediates to the stripe, which results in the huge overpotential for the PCET step D in Equation (10). It is therefore clear that the reduction in adsorption energy due to the solvent can play a role in reducing the overpotential.
Figure 5 shows the energy-level diagram for the four PCET steps, calculated using the methodology summarized in Section 2.3, corresponding to the data in Table 3. The reduction in the adsorption energy leads to the increased step heights of the PCET steps A to C (Equations (7)–(9)), and hence to a net reduction in the overpotential determining step D. Considering the water solvent, PCET step D has now a height of 2.0 eV, compared with 2.38 eV in the vacuum case. The calculated reaction overpotential is therefore reduced to 0.77 eV by the introduction of the implicit solvent. This difference shows how important it is to properly account for the liquid environment in which the OER takes place. Calculations in vacuum lead to similar adsorption structures, but the relative energies of the reaction intermediates may change considerably, as is the case in this system.

4. Conclusions

We performed DFT-based computer simulations of the OER at Mo-edges of MoS2 sheets. This TMD is intensely studied in the context of various applications in catalysis, and experimentally shown to be able to evolve O2 [20], but with low efficiency. In a previous DFT study [19], the binding strength of some of the reaction intermediates was found to be unfavorable for the OER, leading to the observed large overpotential and low activity for the O2 evolution.
These previous DFT simulations were based on a computational model where the catalyst and the adsorbed reaction intermediates are placed in a vacuum, at variance with the experimental setup, where the catalyst is in aqueous solution. The key point of the present work is to go beyond this common approximation and estimate the effect of a water solvent surrounding the MoS2 edges. This is achieved by applying an implicit solvent model, which can account for some of the main effects that the solvent has on the system by embedding the system in a realistic dielectric medium.
Our calculations, which include such solvent effects, show that the adsorption geometries and reaction mechanics of the OER are only very weakly altered by the presence of the solvent. Comparing the system’s electronic charge density with and without the dielectric medium reveals that the main effect of the solvent is localized in space, close to the adsorbed reaction intermediates. These changes in the charge density, and hence in the electrostatic interactions, have a strong effect on the energies of adsorption. In the case of the Mo-edges of MoS2 sheets studied here, the variance in adsorption energies leads to a reduction in the OER’s overpotential from 1.15 eV to 0.77 eV. This considerable change in the overpotential exemplifies that DFT-based computer simulations of reactions in liquid environments should study solvation effects whenever possible.
In a future study, we plan to extend the current computational approach by going beyond a continuum solvent model and explicitly including water molecules. Such an approach [50], which combines molecular dynamics simulations with the continuum models employed here, will allow for an even more detailed view of the energetics and reaction mechanics of the OER in a real application.

Author Contributions

Conceptualization, R.G.; methodology, R.G.; calculations, E.G. and R.G.; writing, E.G. and R.G. All authors have read and agreed to the published version of the manuscript.

Funding

Work supported by the Ministerio de Ciencia e Innovación of Spain (Grant PID2019-104924RB-I00 funded by MCIN/AEI/10.13039/501100011033) and the University of Valladolid (GIR Nanostructure Physics). E.G. acknowledges a postdoctoral contract with the University of Valladolid.

Data Availability Statement

Not applicable.

Acknowledgments

E.G. would like to thank the ICTP and the Simons Foundation for providing the opportunity to undertake this research by supporting her visits to the ICTP as an ICTP-Simons Associate.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hota, P.; Das, A.; Maiti, D.K. A short review on generation of green fuel hydrogen through water splitting. Int. J. Hydrogen Energy 2023, 48, 523–541. [Google Scholar] [CrossRef]
  2. Zhang, H.; Wang, J.; Zhao, X.; Yang, J.; Bu sinnah, Z.A. Modeling a hydrogen-based sustainable multi-carrier energy system using a multi-objective optimization considering embedded joint chance constraints. Energy 2023, 278, 127643. [Google Scholar] [CrossRef]
  3. Solati, N.; Karakaya, C.; Kaya, S. Advancing the Understanding of the Structure–Activity–Durability Relation of 2D MoS 2 for the Hydrogen Evolution Reaction. ACS Catal. 2023, 13, 342–354. [Google Scholar] [CrossRef]
  4. de Oliveira, C.C.; Autreto, P.A.d.S. Optimized 2D nanostructures for catalysis of hydrogen evolution reactions. MRS Adv. 2023, 8, 307–310. [Google Scholar] [CrossRef]
  5. Wang, F.; Shifa, T.A.; Zhan, X.; Huang, Y.; Liu, K.; Cheng, Z.; Jiang, C.; He, J. Recent advances in transition-metal dichalcogenide based nanomaterials for water splitting. Nanoscale 2015, 7, 19764–19788. [Google Scholar] [CrossRef] [PubMed]
  6. Manzeli, S.; Ovchinnikov, D.; Pasquier, D.; Yazyev, O.V.; Kis, A. 2D transition metal dichalcogenides. Nat. Rev. Mater. 2017, 2, 17033. [Google Scholar] [CrossRef]
  7. Chhowalla, M.; Shin, H.S.; Eda, G.; Li, L.J.; Loh, K.P.; Zhang, H. The chemistry of two-dimensional layered transition metal dichalcogenide nanosheets. Nat. Chem. 2013, 5, 263–275. [Google Scholar] [CrossRef]
  8. Politano, A.; Chiarello, G.; Kuo, C.N.; Lue, C.S.; Edla, R.; Torelli, P.; Pellegrini, V.; Boukhvalov, D.W. Tailoring the Surface Chemical Reactivity of Transition-Metal Dichalcogenide PtTe2 Crystals. Adv. Funct. Mater. 2018, 28, 1706504. [Google Scholar] [CrossRef]
  9. Ye, Z.; Tan, C.; Huang, X.; Ouyang, Y.; Yang, L.; Wang, Z.; Dong, M. Emerging MoS2 Wafer-Scale Technique for Integrated Circuits. Nano-Micro Lett. 2023, 15, 38. [Google Scholar] [CrossRef]
  10. Li, X.; Wei, Y.; Wang, Z.; Kong, Y.; Su, Y.; Lu, G.; Mei, Z.; Su, Y.; Zhang, G.; Xiao, J.; et al. One-dimensional semimetal contacts to two-dimensional semiconductors. Nat. Commun. 2023, 14, 111. [Google Scholar] [CrossRef]
  11. Liu, L.; Kumar, S.B.; Ouyang, Y.; Guo, J. Performance Limits of Monolayer Transition Metal Dichalcogenide Transistors. IEEE Trans. Electron Devices 2011, 58, 3042–3047. [Google Scholar] [CrossRef] [Green Version]
  12. Kumar, A.; Faella, E.; Durante, O.; Giubileo, F.; Pelella, A.; Viscardi, L.; Intonti, K.; Sleziona, S.; Schleberger, M.; Di Bartolomeo, A. Optoelectronic memory in 2D MoS2 field effect transistor. J. Phys. Chem. Solids 2023, 179, 111406. [Google Scholar] [CrossRef]
  13. Li, X.; Wang, R.; Wu, Q.; Yu, Y.; Gao, T.; Yao, T.; Wang, X.; Han, J.; Song, B. Synergistically Designed Dual Interfaces to Enhance the Electrochemical Performance of MoO2/MoS2 in Na- and Li-Ion Batteries. Small 2023, 19, 2206940. [Google Scholar] [CrossRef]
  14. Bai, Z.; Xu, H.; Li, G.; Yang, B.; Yao, J.; Guo, K.; Wang, N. MoS2 Nanosheets Decorated with Fe3O4 Nanoparticles for Highly Efficient Solar Steam Generation and Water Treatment. Molecules 2023, 28, 1719. [Google Scholar] [CrossRef] [PubMed]
  15. Guo, M.; Yuan, B.; Sui, Y.; Xiao, Y.; Dong, J.; Yang, L.; Bai, L.; Yang, H.; Wei, D.; Wang, W.; et al. Rational design of molybdenum sulfide/tungsten oxide solar absorber with enhanced photocatalytic degradation toward dye wastewater purification. J. Colloid Interface Sci. 2023, 631, 33–43. [Google Scholar] [CrossRef]
  16. Li, G.L.; Miao, Y.Y.; Qiao, X.Y.; Wang, T.Y.; Deng, F. Engineering edge sites based on NiS2/MoS2/CNTs heterojunction catalyst for overall water splitting. Appl. Surf. Sci. 2023, 615, 156309. [Google Scholar] [CrossRef]
  17. Chen, B.; Hu, P.; Yang, F.; Hua, X.; Yang, F.F.; Zhu, F.; Sun, R.; Hao, K.; Wang, K.; Yin, Z. In Situ Porousized MoS2 Nano Islands Enhance HER/OER Bifunctional Electrocatalysis. Small 2023, 19, 2207177. [Google Scholar] [CrossRef] [PubMed]
  18. Ouyang, Y.; Ling, C.; Chen, Q.; Wang, Z.; Shi, L.; Wang, J. Activating Inert Basal Planes of MoS2 for Hydrogen Evolution Reaction through the Formation of Different Intrinsic Defects. Chem. Mater. 2016, 28, 4390–4396. [Google Scholar] [CrossRef]
  19. German, E.; Gebauer, R. Why are MoS2 monolayers not a good catalyst for the oxygen evolution reaction? Appl. Surf. Sci. 2020, 528, 146591. [Google Scholar] [CrossRef]
  20. Pesci, F.M.; Sokolikova, M.S.; Grotta, C.; Sherrell, P.C.; Reale, F.; Sharda, K.; Ni, N.; Palczynski, P.; Mattevi, C. MoS2/WS2 Heterojunction for Photoelectrochemical Water Oxidation. ACS Catal. 2017, 7, 4990–4998. [Google Scholar] [CrossRef] [Green Version]
  21. Wang, Z.; Tang, M.T.; Cao, A.; Chan, K.; Nørskov, J.K. Insights into the Hydrogen Evolution Reaction on 2D Transition-Metal Dichalcogenides. J. Phys. Chem. C 2022, 126, 5151–5158. [Google Scholar] [CrossRef]
  22. Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864–B871. [Google Scholar] [CrossRef] [Green Version]
  23. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. [Google Scholar] [CrossRef] [Green Version]
  24. Klamt, A.; Schüürmann, G. COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc. Perkin Trans. 1993, 2, 799–805. [Google Scholar] [CrossRef]
  25. Tomasi, J.; Persico, M. Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent. Chem. Rev. 1994, 94, 2027–2094. [Google Scholar] [CrossRef]
  26. Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999–3093. [Google Scholar] [CrossRef] [PubMed]
  27. Onsager, L. Electric Moments of Molecules in Liquids. J. Am. Chem. Soc. 1936, 58, 1486–1493. [Google Scholar] [CrossRef]
  28. Corni, S.; Tomasi, J. Enhanced response properties of a chromophore physisorbed on a metal particle. J. Chem. Phys. 2001, 114, 3739–3751. [Google Scholar] [CrossRef]
  29. Fattebert, J.L.; Gygi, F. Density functional theory for efficientab initio molecular dynamics simulations in solution. J. Comput. Chem. 2002, 23, 662–666. [Google Scholar] [CrossRef] [PubMed]
  30. Fattebert, J.L.; Gygi, F. First-principles molecular dynamics simulations in a continuum solvent. Int. J. Quantum Chem. 2003, 93, 139–147. [Google Scholar] [CrossRef]
  31. Andreussi, O.; Dabo, I.; Marzari, N. Revised self-consistent continuum solvation in electronic-structure calculations. J. Chem. Phys. 2012, 136, 064102. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Buongiorno Nardelli, M.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys. Condens. Matter 2017, 29, 465901. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Chen, Y.; Deng, J.J.; Yao, W.W.; Gurti, J.I.; Li, W.; Wang, W.J.; Yao, J.X.; Ding, X.L. Non-stoichiometric molybdenum sulfide clusters and their reactions with the hydrogen molecule. Phys. Chem. Chem. Phys. 2021, 23, 347–355. [Google Scholar] [CrossRef]
  34. Wang, Y.Y.; Deng, J.J.; Wang, X.; Che, J.T.; Ding, X.L. Small stoichiometric (MoS2) n clusters with the 1T phase. Phys. Chem. Chem. Phys. 2018, 20, 6365–6373. [Google Scholar] [CrossRef] [PubMed]
  35. Ulman, K.; Nguyen, M.T.; Seriani, N.; Gebauer, R. Passivation of surface states of α-Fe2O3(0001) surface by deposition of Ga2O3 overlayers: A density functional theory study. J. Chem. Phys. 2016, 144, 094701. [Google Scholar] [CrossRef]
  36. Nguyen, M.T.; Farnesi Camellone, M.; Gebauer, R. On the electronic, structural, and thermodynamic properties of Au supported on α-MoS2 surfaces and their interaction with CO. J. Chem. Phys. 2015, 143, 034704. [Google Scholar] [CrossRef] [PubMed]
  37. Byskov, L.S.; Nørskov, J.K.; Clausen, B.S.; Topsøe, H. Edge termination of MoS2 and CoMos catalyst particles. Catal. Lett. 2000, 64, 95–99. [Google Scholar] [CrossRef]
  38. Wang, H.; Tsai, C.; Kong, D.; Chan, K.; Abild-Pedersen, F.; Nørskov, J.K.; Cui, Y. Transition-metal doped edge sites in vertically aligned MoS2 catalysts for enhanced hydrogen evolution. Nano Res. 2015, 8, 566–575. [Google Scholar] [CrossRef]
  39. Kim, T.; Lee, T.H.; Park, S.Y.; Eom, T.H.; Cho, I.; Kim, Y.; Kim, C.; Lee, S.A.; Choi, M.J.; Suh, J.M.; et al. Drastic Gas Sensing Selectivity in 2-Dimensional MoS2 Nanoflakes by Noble Metal Decoration. ACS Nano 2023, 17, 4404–4413. [Google Scholar] [CrossRef]
  40. Norskov, J.K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L.; Kitchin, J.R.; Bligaard, T.; Jónsson, H. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108, 17886–17892. [Google Scholar] [CrossRef]
  41. Gebauer, R.; Serra, S.; Chiarotti, G.; Scandolo, S.; Baroni, S.; Tosatti, E. Noncolinear spin polarization from frustrated antiferromagnetism: A possible scenario for molecular oxygen at high pressure. Phys. Rev. B 2000, 61, 6145–6149. [Google Scholar] [CrossRef]
  42. Groß, A. Challenges in the modeling of elementary steps in electrocatalysis. Curr. Opin. Electrochem. 2023, 37, 101170. [Google Scholar] [CrossRef]
  43. Chan, K.; Norskov, J.K. Electrochemical Barriers Made Simple. J. Phys. Chem. Lett. 2015, 6, 2663–2668. [Google Scholar] [CrossRef]
  44. Righi, G.; Fabris, S.; Piccinin, S. Oxygen Evolution Reaction on the Fe3O4(001) Surface: Theoretical Insights into the Role of Terminal and Bridging Oxygen Atoms. J. Phys. Chem. C 2021, 125, 18752–18761. [Google Scholar] [CrossRef]
  45. Prandini, G.; Marrazzo, A.; Castelli, I.E.; Mounet, N.; Marzari, N. Precision and efficiency in solid-state pseudopotential calculations. npj Comput. Mater. 2018, 4, 72. [Google Scholar] [CrossRef] [Green Version]
  46. Lejaeghere, K.; Bihlmayer, G.; Bjorkman, T.; Blaha, P.; Blugel, S.; Blum, V.; Caliste, D.; Castelli, I.E.; Clark, S.J.; Dal Corso, A.; et al. Reproducibility in density functional theory calculations of solids. Science 2016, 351, aad3000. [Google Scholar] [CrossRef] [Green Version]
  47. Monkhorst, H.J.; Pack, J.D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192. [Google Scholar] [CrossRef]
  48. Marzari, N.; Vanderbilt, D.; De Vita, A.; Payne, M.C. Thermal Contraction and Disordering of the Al(110) Surface. Phys. Rev. Lett. 1999, 82, 3296–3299. [Google Scholar] [CrossRef] [Green Version]
  49. 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. 2010, 132, 154104. [Google Scholar] [CrossRef] [Green Version]
  50. Zhao, W.N.; Liu, Z.P. Mechanism and active site of photocatalytic water splitting on titania in aqueous surroundings. Chem. Sci. 2014, 5, 2256–2264. [Google Scholar] [CrossRef]
Figure 1. Structure of the MoS2 model systems. (a) Pristine MoS2 monolayer, top view. (b) Side view of a pristine monolayer. (c) Infinite stripe model in top view. The two inequivalent edges (sulfur edge and molybdenum edge, respectively) can be discerned. See the main text for a description of the edges. The supercell along the axis of the stripe is indicated by the dashed lines. Colour code: Sulfur (yellow) and molybdenum (dark blue).
Figure 1. Structure of the MoS2 model systems. (a) Pristine MoS2 monolayer, top view. (b) Side view of a pristine monolayer. (c) Infinite stripe model in top view. The two inequivalent edges (sulfur edge and molybdenum edge, respectively) can be discerned. See the main text for a description of the edges. The supercell along the axis of the stripe is indicated by the dashed lines. Colour code: Sulfur (yellow) and molybdenum (dark blue).
Molecules 28 05182 g001
Figure 2. Infinite stripe model of MoS2 with key reaction intermediates: (a) Oxygen O * . (b) Hydroxyl group HO * . (c) HOO * , which binds in the lowest energy configuration as HO * + O * . Colour code: Sulfur (yellow), molybdenum (dark blue), oxygen (red), hydrogen (grey).
Figure 2. Infinite stripe model of MoS2 with key reaction intermediates: (a) Oxygen O * . (b) Hydroxyl group HO * . (c) HOO * , which binds in the lowest energy configuration as HO * + O * . Colour code: Sulfur (yellow), molybdenum (dark blue), oxygen (red), hydrogen (grey).
Molecules 28 05182 g002
Figure 3. Spin resolved density of states (DOS) and projected DOS (pDOS) of the systems considered. The spin-up and spin-down DOS is shown as positive and negative DOS, respectively. (a) Pristine MoS2 stripe without adsorbates. (b) Stripe with HO * adsorbate. (c) Stripe with O * adsorbate. (d) Stripe with O*HO * adsorbate. The pDOS of the adsorbed oxygen atoms is multiplied by a factor of 10 for better visual clarity.
Figure 3. Spin resolved density of states (DOS) and projected DOS (pDOS) of the systems considered. The spin-up and spin-down DOS is shown as positive and negative DOS, respectively. (a) Pristine MoS2 stripe without adsorbates. (b) Stripe with HO * adsorbate. (c) Stripe with O * adsorbate. (d) Stripe with O*HO * adsorbate. The pDOS of the adsorbed oxygen atoms is multiplied by a factor of 10 for better visual clarity.
Molecules 28 05182 g003
Figure 4. Reaction intermediate HO * adsorbed on the molybdenum edge of the MoS2 stripe. A contour plot of the charge density difference between the solvated system and the system in vacuum is presented. The geometry is that of the solvated system. The blue-coloured isosurface corresponds to a reduction in the charge density in the solvated system, while the red colour corresponds to an increase in charge density due to solvation.
Figure 4. Reaction intermediate HO * adsorbed on the molybdenum edge of the MoS2 stripe. A contour plot of the charge density difference between the solvated system and the system in vacuum is presented. The geometry is that of the solvated system. The blue-coloured isosurface corresponds to a reduction in the charge density in the solvated system, while the red colour corresponds to an increase in charge density due to solvation.
Molecules 28 05182 g004
Figure 5. Energy diagram for the four PCET steps of the OER, Equations (7)–(10). The step heights are strongly dependent on the solvation of the catalyst.
Figure 5. Energy diagram for the four PCET steps of the OER, Equations (7)–(10). The step heights are strongly dependent on the solvation of the catalyst.
Molecules 28 05182 g005
Table 1. Adsorption geometries of reaction intermediates at the MoS2 molybdenum edge with and without solvation in water. Distances in Å.
Table 1. Adsorption geometries of reaction intermediates at the MoS2 molybdenum edge with and without solvation in water. Distances in Å.
SolvationIntermediated(S−O)d(O−H)d(S−O(H))
O*HO * 1.4810.9791.704
VacuumHO * -0.9801.702
O * 1.480--
O*HO * 1.4850.9861.661
WaterHO * -0.9861.661
O * 1.482--
Table 2. Adsorption energies of reaction intermediates at the MoS2 molybdenum edge with and without solvation in water. Energies in eV.
Table 2. Adsorption energies of reaction intermediates at the MoS2 molybdenum edge with and without solvation in water. Energies in eV.
SolvationIntermediate E ads
O*HO * −2.375
VacuumHO * −1.651
O * −1.225
O*HO * −2.003
WaterHO * −1.530
O * −0.946
Table 3. Step heights for PCET steps A–D, Equations (7)–(10). Energies in eV.
Table 3. Step heights for PCET steps A–D, Equations (7)–(10). Energies in eV.
SolvationPCET Step Δ E
A0.809
Vacuum [19]B0.426
C1.311
D2.375
A0.929
WaterB0.585
C1.403
D2.004
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

German, E.; Gebauer, R. The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations. Molecules 2023, 28, 5182. https://doi.org/10.3390/molecules28135182

AMA Style

German E, Gebauer R. The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations. Molecules. 2023; 28(13):5182. https://doi.org/10.3390/molecules28135182

Chicago/Turabian Style

German, Estefania, and Ralph Gebauer. 2023. "The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations" Molecules 28, no. 13: 5182. https://doi.org/10.3390/molecules28135182

APA Style

German, E., & Gebauer, R. (2023). The Oxygen Evolution Reaction at MoS2 Edge Sites: The Role of a Solvent Environment in DFT-Based Molecular Simulations. Molecules, 28(13), 5182. https://doi.org/10.3390/molecules28135182

Article Metrics

Back to TopTop