High-entropy materials for electrocatalytic applications: a review of first principles modeling and simulations

Abstract

High-entropy materials, for both complexity in structure and superiority in performance, have been widely confirmed to be one possible kind of advanced electrocatalyst. Significant efforts have been dedicated to modeling the atomic-level details of high-entropy catalysts to improve the viability for bottom-up design of advanced electrocatalysts. In this review, first, we survey developments in various modeling methods that are based on density functional theory. We review progress in density functional theory simulations for emulating different high-entropy electrocatalysts. Then, we review the advancements in simulations of high-entropy materials for electrocatalytic applications. Finally, we present prospects in this field.

Abbreviations: HEMs: high-entropy materials; CCMs: compositionally complex materials; DFT: density functional theory; LDA: local density approximation; GGA: generalized gradient approximation; VASP: Vienna Ab initio simulation package; ECP: effective core potential; PAW: projector-augmented wave potential; VCA: virtual crystal approximation; CPA: coherent potential approximation; SQS: special quasi-random structures; SSOS: small set of ordered structures; SLAE: similar local atomic environment; HEAs: high-entropy alloys; FCC: face-centered cubic; BCC: body-centered cubic; HCP: hexagonal close-packed; ORR: oxygen reduction reaction; OER: oxide evolution reaction; HER: hydrogen evolution reaction; RDS: rate-limiting step; AEM: adsorbate evolution mechanism; LOM: lattice oxygen oxidation mechanism; HEOs: high-entropy oxides; OVs: oxygen vacancies; PDOS: projected densities of states; ADR: ammonia decomposition reaction; NRR: nitrogen reduction reaction; CO₂RR: CO₂ reduction reaction; TMDC: transition metal dichalcogenide; TM: transition metal; AOR: alcohol oxidation reaction; GOR: glycerol oxidation reaction; UOR: urea oxidation reaction; HEI: high-entropy intermetallic.

GRAPHICAL ABSTRACT

IMPACT STATEMENT

This paper reviews recent developments in the field of atomistic simulations of high-entropy electrocatalysts, one of emerging state-of-the-art catalytic materials.

KEYWORDS:

• High-entropy materials
• first principles
• modeling
• density functional theory
• electrocatalysis

1. Introduction

HEMs form a quintessential class of CCMs. They can form simple structures with chemical compositions of multiple elements [1–3]. The combination of properties that arise from both the structure and composition, has led to interesting applications [4]. HEMs show remarkable mechanical properties [5], corrosion resistance [6], and radiation resistance [7]. Recently, nanostructured HEMs have been developed and have shown promise for superior catalytic properties [8]. The birth of the design concept based on multiple elements, HEMs, provides a promising path to break the shackles of electrocatalyst composition design. It overcomes the strict limit of the chemical composition for conventional electrocatalysts. HEMs have shown great potential for electrocatalytic applications and have developed into a rapidly expanding field during this period. Various kinds of HEMs have been reported, e.g. HEAs and high entropy ceramics. High entropy metallic glasses [9], oxides [10], carbides [11], borides [12], nitrides [13], sulfides [14], silicides [15], fluorides [16], phosphides [17], and intermetallics [18] have also been explored. Atomic and electronic structures play vital roles in the properties of HEMs. It provides a variety of potential applications for different high-entropy materials. However, the underlying mechanisms of intrinsic structural stability, mechanical capacity and catalytic applications remain unknown.

Combined with experiments, theoretical models and simulations can serve as effective methods to improve the basic understanding of the relationship between the structure and properties and enhance the performance by optimizing the material composition and structure. For high-entropy catalysts, the multielemental synergy enables a diverse range of adsorption sites that could accelerate multistep electrocatalytic reactions that require multifunctional catalysts. However, the wide range of constitutive compositions and complex atomic arrangements also create grand challenges in the understanding of real active sites and the underlying mechanism [8]. The computational method, especially combined with the experimental results could provide a new platform to disclose the obscure composition-structure-performance relationships. Molecular dynamics simulation and DFT calculation are the most widespread methods for the prediction of the structure and the fundamental understanding of various (structural/electronic/thermal) properties [19]. Molecular dynamics simulation is usually used at a larger scale (over 100∼10,000 atoms), which is appropriate for covering materials with many elements. However, the molecular dynamics method cannot consider the interatomic potentials to solve the complex chemistry of multi-component ceramics and requires a large number of computing resources. Thus, the microscopic electronic features cannot be precisely analyzed by molecular dynamics. On the other hand, DFT could be used as a useful computational tool to investigate the relationship between material properties and physical/chemical properties, complementing experiments and acting as a helpful predictive tool for identifying and characterizing materials with HEMs [20]. These detailed simulations are related to the theory of quantum mechanics and inherent constants and thus the prediction and exploration of intrinsic material properties without extra empirical input [21].

Furthermore, there are still significant challenges to the theoretical modeling and simulations [4]. The first is the difficulty in constructing empirical atomic interaction models. Due to the large number of chemical interactions involved, high mixing entropy, serious lattice distortion, and short-range order in HEMs, the physical model and simulation for traditional binary or ternary materials are no longer suitable for studying the enhanced mechanism of HEMs. The second challenge occurs in the first-principles simulations. DFT simulations are generally computationally expensive. However, they are more prohibitive for HEMs due to the necessity of using large supercells (up to more than 1000 atoms) to represent the non-stoichiometric compositions and complex microstructures. Therefore, to study the unique characteristics of their microstructures and properties by DFT calculations, it is necessary to develop a more efficient method of atomistic simulations for various HEMs.

In this review, we introduce the current development of atomistic simulations of HEMs for catalytic applications. The following major topics are presented. Section 2 introduces the basic development and concepts of DFT calculations. The most used methodologies of modeling for HEMs and DFT simulations for HEM electrocatalysts are reported in Sections 3 and 4. The applications of HEMs in electrocatalytic field are discussed in Section 5. This review closes with the consideration of the suggested outlooks (Section 6).

2. Theoretical basis of DFT

The DFT concept can be first traced back to the Thomas–Fermi model in the 1920s [22]. The introduction of the Hohenberg–Kohn theorem in 1964 accelerated the establishment of DFT [23]. Then, the Kohn–Sham equation made DFT methods a reality [24]. DFT originates in the observation that it is prohibitively expensive to numerically estimate the many-electron wave function Ψ. For N electrons, Ψ is a function of 3N variables such as the Cartesian coordinates x₁, y₁, z₁, x₂  … , z_N. The amounts of required bits to store such a function scale exponentially with the number of electrons N. This unavoidable fact critically limits the applicability of numerical methods for calculations that are dependent on wave functions [23], as shown in several applications by the Hartree–Fock method [25].

The use of the electron density ρ could overcome the aforementioned exponential barrier. Due to the nature related to the 3 spatial coordinates, the use of the electron density averted the limitation of exponential scaling. As a real function in a drastically smaller number of dimensions, electron density might contain less information than a wave function. However, the Hohenberg-Kohn theorem pointed out that electron density (ρ) implicitly carries all necessary information about the properties of the ground state [23]. In principle, there exists a ρ-functional for most physically observable measurements. However, these functionals are uncertain, and we can only perform approximate DFT calculations.

The Kohn–Sham DFT depends on the solution of the noninteracting problem, by defining effective Coulomb potentials derived from external (ionic) potentials V and electrons, (1) $V_C(\overrightarrow{r})=V(\overrightarrow{r})+\int \frac{\rho ({\overrightarrow{r}}^{\mathrm{\prime }})}{|\overrightarrow{r}-{\overrightarrow{r}}^{\mathrm{\prime }}|}d{\overline{r}}^{\mathrm{\prime }}$(1) and the exchange correlation potential V_xc. The ground-state density is the solution of the self-consistent equations (2) $\left(-\frac{1}{2}{\mathrm{\nabla }}^2+V_C(\overrightarrow{r})+V_{XC}(\overrightarrow{r})-E_i\right){\psi }_i(\overrightarrow{r})=0$(2) and (3) $\rho (\overrightarrow{r})=\sum _i|{\psi }_i(\overrightarrow{r}){|}^2$(3) ${\psi }_i$ are the wave functions for the noninteracting Kohn–Sham quasiparticles. Thus, the complex N electron problem is simplified to N one-electron problems. The equation in this form is accurate, and the Kohn–Sham eigenvalue E_i can be used to calculate the total ground state energy. However, the form of the exchange correlation functional V_xc is still unclear. Although V_xc is often written as a local function, the value of V_xc at $\overrightarrow{r}$ relies on the distribution of electron density in the whole space. Compared to the Hartree potential, the whole magnitude of V_xc is much smaller. However, errors in approximating V_xc can cause inevitable qualitative errors. The disadvantages of the present approximations to exchange correlation functionals lead to the unreliability of DFT in some materials systems and the choice of the functional requires the experience of researchers.

Exchange–correlation functionals can be classified on the basis of features. To make clear classification of the exchange correlation functional, Jacob's ladder was proposed [26], based on the descriptions of the Hartree world without any exchange–correlation and guaranteed much higher chemical accuracy. The first step of Jacob's ladder is LDA, and the following are the GGA, meta-GGA, and finally, hybrid functionals. Each rung on the ladder should contain the former functionality and perform better than the previous ones. The LDA and the GGA [27,28] have been widely used. For single materials and ordered compounds, LDA and GGA have been demonstrated to be reliable for experimental material properties [29]. For HEMs, GGA has been widely used in most studies. Some of the target alloys are magnetic and need to be treated with spin polarization [30]. Hybrid functionals can also be used and incorporated into the Hartree–Fock theory which is widely applied for studying electronic properties of materials such as PBE0 due to Perdew-Burke-Ernzerhog and B3LYP by Becke-Lee–Yang–Parr [31], especially in the modeling of water molecules and other complex compounds.

Over the past few decades, many quantum-chemical software packages have been explored for the prediction of material properties. Most of the popular DFT computer programs are summarized in Table 1. Among them, the Gaussian program is the most extensively used for chemical systems, and the VASP is widely used for physics and materials science. Moreover, the DMol3 code based on a numerical basis is applicable for isolated molecules, solids and surfaces [32], which exhibits low computational cost. To effectively lower the computational cost, the ECP [33], pseudopotential [34], and PAW methods [35] are introduced for relatively inert core electrons.

Table 1. Various developed software and the related basis set, functionals, and core potential.

Software	Basis set	Functionals	Core treatment
Gaussian^1	Gaussian basis	All functionals	Effective Core Potential
ORCA^2	Gaussian basis	GGA, meta-GGA, hybrid functionals	Effective Core Potential
Q-chem^3	Gaussian basis	LDA, GGA, and global hybrid functionals	Effective Core Potential
DALTON 2.0^4	Gaussian basis	LDA, GGA, hybrid functionals	Effective Core Potential
VASP^5	Plane wave	LDA, GGA, meta-GGA, and hybrid functionals	Ultrasoft pseudopotentials or projector-augmented wave potentials
CASTEP^6	Plane wave	LDA, GGA	Ultrasoft pseudopotentials / Nor-conserving pseudopotentials
ABINIT^7	Plane wave	LDA, GGA	Norm-conserving pseudopotentials
CP2K^8	Mixed Gaussian and plane waves	LDA, GGA, meta-GGA, hybrid functionals, and double-hybrid XC functionals	Pseudopotential, all-electron
WIEN 2K^9	The augmented plane wave plus local orbitals (APW + lo) method	LDA, GGA, meta-GGA, and hybrid functionals	Full potentials
DMol3^10	Numerical basis	LDA, GGA, B3LYP, meta-GGA, Hartree–Fock exact exchange	All-electron, elective Core potentials, all-electron relativistic, and semi-core pseudopotential
ADF^11	Slater basis	Numerical integration scheme	Frozen core approximation
Octopus^12	Numerical basis set, plane waves, and atomic orbitals	LDA, GGA	Optimized effective potential method (OEP)
LMTART^13	Atom centered local muffin-tin orbitals	LDA, GGA	Full potential
GPAW^14	Localized atomic orbitals and plane waves	LDA, GGA	Pseudopotential, all-electron

Notes: The websites of the software are listed as follows: 1. http://www.gaussian.com/; 2. https://orcaforum.kofo.mpg.de; 3. https://www.q-chem.com/; 4. https://www.daltonprogram.org/; 5. https://www.vasp.at/; 6. http://www.castep.org/; 7. https://www.abinit.org/; 8. https://www.cp2k.org/; 9. http://www.wien2k.at/; 10. https://wiki.fysik.dtu.dk/ase/ase/calculators/dmol.html; 11. https://www.scm.com/amsterdam-modeling-suite/; 12. https://octopus-code.org/wiki/Main_Page; 13. http://www.fkf.mpg.de/andersen/docs/lmtoart_programs.html; 14. https://wiki.fysik.dtu.dk/gpaw/.

3. Atomistic modeling of HEMs

DFT calculations are a significant method to investigate and predict the various properties of targeted materials. The limitations of the first-principles method arise from the approximations required to deal with the high computational cost that increases exponentially with the size and complexity of simulated HEMs. One of the major limitations of DFT is the limited number of atoms it can handle, and the chemical disorder in HEMs poses an especially difficult problem. Due to the chemical disorder of HEMs, it is a huge challenge to build appropriate HEMs models (the number of atoms lower than 1000) to meet the requirement of chemical randomness [36]. Recently, two strategies have been proposed to solve the intractable problem: (a) the ‘effective medium’ method, including the VCA method [37] and CPA method [38]; (b) supercell methods, including SQS [39], SSOS [40], and SLAE structures [41]. These methods have their own advantages and shortcomings. The practical use relies on the specific application.

3.1. VCA method

VCA is usually suitable and useful for materials in which the constituent elements are adjacent in the periodic table of elements [37]. The core of the VCA is the replacement of the true metal atoms by a ‘virtual’ metal, which is composed of a weighted average between the heterogeneous elements in the parent materials. Generally, the potentials that represent atoms of two or more elements are averaged into a composite atomic potential. It is an overly simplified approach to substitutional solid solutions describing the elastic properties of a specific material. The supercell is not necessarily needed, and the calculation time can be significantly saved. For instance, as shown in Figure 1(a), Wang et al. calculated the lattice constants of various HEAs, FeCoNiCr, FeCoNiCrMn, and FeCoNiCrMnGe, via the VCA methods [38]. It was demonstrated that for HEAs composed of atoms with similar properties, the VCA method exhibits more efficiency than SQS methods for the estimation of LCs. Arun et al. introduced a facile HEA design with engineered configuration entropy and a band tuning approach [43]. With the help of DFT calculations via the VCA model, the Na atoms were verified to locate Pb/Sn substitutional doping sites of the original PbSnTeSe HEA (Figure 1(b)). The tiny doping of Na (up to 0.01 at.%) in the Pb sites exhibits a prominent effect on the band tuning. Previous reports have also demonstrated the reliability of the VCA method for the simulation of alloys containing refractory elements [44,45].

Figure 1. (a) VCA models of FeCoNiCr, FeCoNiCrMn, and FeCoNiCrMnGe, respectively. Reprinted with permission from Ref. [42] Copyright 2011, Elsevier. (b) Crystal structure of PbSnTeSe HEA and Pb_0.99SnTeSe-Na_0.01 HEA. Reprinted with permission from Ref. [43] Copyright 2022, Elsevier. (c) CPA models for the equimolar ABCDE high-entropy alloys (HEAs). Reprinted with permission from Ref. [46] Copyright 2017, Creative Commons Attribution License (CC BY).

However, the VCA model is not suitable for calculations of the full material properties of HEMs, such as lattice distortion. The VCA methods cannot be used for HEMs systems consisting of heterogeneous atoms with large differences in properties. Thus, more atomistic modeling methods have been developed to solve this problem.

3.2. CPA method

The CPA method was proposed to simulate random media, which depends on the mean field approach for the consideration of chemical disorder in high entropy material systems. Figure 1(c) shows an illustration of the CPA method for HEMs [46]. The quinary HEA composed of the various elements A, B, C, D, and E is substituted by an effective medium. There are two kinds of approximation methods. The first one is based on the assumption that the local potentials (P_A, P_B, P_C, P_D, P_E) around a certain kind of atom from the HEAs are equal, which ignores the local environment effect. Another is that the target system is replaced by a monatomic medium determined by the site-independent coherent potential $\tilde{P}$. For the CPA method, the elements in the HEMs are introduced into an equal medium, which is determined self-consistently by the mutual effect of the whole atoms, where information on the charge distribution due to the presence of defects is not provided. Based on this facile mean field approach, CPA methods could well simulate chemical disorder in a single primitive cell, which significantly lowers the computational requirements. In addition, CPA-based methods are suitable for various composition designs. The CPA also shows high efficiency in the calculation of the ground-state properties of HEMs, such as lattice parameters, bulk moduli, and mixing enthalpies. Due to the easy extraction of the energetics of HEMs systems, it could also be used for multiscale simulations. These unique merits ensure a high-efficiency computational exploration of a wide range of HEMs. CPA is also coupled with other methods, such as exact muffin-tin orbitals (EMTO) [47] or Korringa-Kohn-Rostoker (KKR) methods [48].

3.3. SQS method

Compared with the previous two methods, the difference in the supercell method is reflected in the simulation of the distribution of the local chemical environment in HEMs and how to reduce the computational demands. As shown in Figure 2(a), a disordered chemical configuration is mimicked on a large and finite supercell with periodic boundary conditions [49].

Figure 2. (a) Schematics for the supercell method. Reprinted with permission from Ref. [49] Copyright 2019, Elsevier. (b) 64-atom SQS used for DFT calculation with atomic species labeled. Reprinted with permission from Ref. [52] Copyright 2022, Elsevier. (c) Schematic representations of SQS and SSOS models for an HEA with FCC lattice. Reprinted with permission from Ref. [56] Copyright 2020, American Physical Society. (d) SLAE models for equiatomic BCC, FCC, and HCP HEAs. Reprinted with permission from Ref. [46] Copyright 2017, Creative Commons Attribution License (CC BY).

The SQS method is a widely used supercell method to build a complex HEMs structure. The SQS concept is based on the set of physical correlation pairs and multilocus correlation functions of totally disordered HEMs. Via the SQS approach, dozens of atoms are needed to realize the chemical randomness of HEMs with a supercell where the Monte Carlo method and genetic algorithms are commonly utilized to create the structures. In this way, supercells built by the SQS method have been applied for HEMs simulations with FCC [46], BCC [50], and HCP [51] structures.

Kim et al. analyzed the interatomic distance distributions with the help of DFT calculations in the SQS model (Figure 2(b)) [52]. It was found that the HEA exhibited large lattice distortion. In the exact consideration of the lattice distortion, the obtained elastic properties well supported the experimental measurements, and the calculated elastic constant values were within 5% of the neutron-diffraction measurements. With the help of DFT calculations in the SQS model, Yao et al. explored the thermodynamic properties of W_x(TaTiVCr)_1-x HEA, and the effect of W content was systematically investigated. The obtained DFT results showed that the bulk modulus of W_x(TaTiVCr)_1-x decreased with increasing temperature, and the softening tendency coincided with the presence of W [53].

DFT simulations based on SQS supercells could obtain exact results under the condition of short-range interatomic interactions in HEMs. In addition, the properties of produced materials are usually dependent on the local atomic arrangements (local lattice distortions, equilibrium bond lengths, density of states, band gap, charge-transfer, mixing, and formation enthalpies) [54]. However, this method is not suitable for systems dominated by long-range interactions.

3.4. SSOS method

For HEMs, the incremental amount of metal elements could result in a significant increase in the correlation functions. Under this circumstance, a small-sized SQS model could not meet the requirements to effectively simulate the statistics of a random materials system [39]. On the other hand, for polymetallic HEMs, the number of correlation functions grows with the increased metallic component. To solve this problem, Jiang and Uberuaga developed a facile SSOS method in 2016 [55]. Via the SSOS method, the actual HEMs system is replaced by a set of small ordered structures, simultaneously combined with approximate weight-averaged properties, leading to a quite different structure from the typical single supercell, to simulate the random state.

Sorkin et al. proposed that the SSOS method has enormous potential for high-throughput DFT calculations of solid-solution HEAs [56]. As shown in Figure 2(c), an SQS model with 125 atoms of a quinary HEA was first built. The equivalent HEA model was simplified by the SSOS method, which could match the equal atomic pair-correlation functions of the SQS model up to the second-nearest neighbor range. Thus, the computational cost for calculating the three 5-atom SSOS model is much lower than that of the SQS model (N = 125). This is beneficial for the practical applications of the SSOS method for high-throughput DFT calculations.

3.5. SLAE method

Different from previously mentioned models, the SLAE model for HEMs is proposed on the basis of local atomic environment settings of all lattice sites, which are similar to each other [52].

To intuitively examine the radial distribution and correlation, standard deviations may be adopted to describe the disorder degree. The SLAE-based supercell for various lattice sites could be formed by random search. The number of atomic pairs on the boundary is equally divided by the two neighboring supercells to treat the periodic boundary condition. As shown in Figure 2(d), the three types of SLAE base models (FCC, BCC and HCP structures) are illuminated, in which the radial distribution is a function of the neighboring pairs of equimolar HEA atoms [57].

4. DFT simulations for HEMs electrocatalysts

In the past decade, the field of atomistic simulations for HEMs with multiple components has progressed rapidly. Most atomistic simulations of HEMs are focused on the explanations of the observed structural, mechanical, and thermodynamic properties. In this section, we present several concepts and descriptors derived from DFT calculations for the prediction of electrocatalyst performance indicators.

4.1. Thermodynamic phase stabilities

Due to the unique features of HEMs, the thermodynamic property and the underlying structure-performance relations is also vital for the further exploitation of potential electrocatalytic applications.

Defect. In particular, structure defects, such as point defects (e.g. vacancies and self-interstitials), line (e.g. dislocations) and planar (e.g. twin, grain and interface boundaries) intrinsic defects, could have a crucial effect on the improvement of electrocatalytic activity on the nanoscale [56]. For an ideal crystal, the atoms should be arranged periodically within the specific structure. However, due to the complex chemical composition and heavy lattice distortion, it is impossible for HEMs, (mostly nanoscale) to fit perfectly into a complete crystal during the synthesis progress, usually under extreme conditions. Defects are produced when the original atoms are arranged incorrectly, which inevitably affects the electronic structure of the materials. It provides a potential design platform to improve the electrocatalytic performance by controllably adjusting the types and/or locations of defects through appropriate strategies. The recently proposed defect engineering has become an effective strategy for tuning the electronic structure and adsorption behavior of catalysts to enhance their activity and stability. For example, Jia et al. reported novel high-entropy metallic glass materials with five equal atomic elements for the alkaline and acidic HER [9]. Furthermore, with the help of a simple dealloying treatment, they obtained a nano-sponge-like architecture. The obtained samples showed superior HER activity, due to the enhanced large specific surface area. To further understand the origin of the catalytic activity, DFT calculations have shown that the interstitial solid solution of P atoms in the Pt–Pd sublattice produces lattice distortion, stabilizes the adsorption/desorption of hydrogen protons, and significantly promotes the HER.

In addition, precise DFT calculations of thermodynamic properties are necessary to understand the formation behaviors of HEMs. The most relevant features include the vacancy/self-interstitial formation energy, vacancy/interstitial migration energy barriers, dislocation simulations, and stacking fault energy [56]. For example, Zhang’s group studied the vacancy formation energies in equiatomic HEAs by first-principles calculations [58]. They reported that the d electrons could affect the distributions of formation and migration energies for point defects in NiCoCr and NiCoFeCr alloys. The formation energies of interstitials in HEAs are lower than those in pure Ni, whereas the formation energies of vacancies are higher. A detailed analysis of electronic properties reveals that the electronic charge deformation flexibility regarding the e_g to t2_g transition has a dominant effect on defect energetics for different elements in HEAs.

Segregation free energy. Owing to the vast compositional space of HEAs, the concentration of the elements can be tuned to obtain specific surface structures to design highly active, stable, and selective high-entropy catalysts. In general, surface segregation in HEAs is in most cases not negligible. Surface segregation in HEAs is not conceptually different from conventional alloys and segregation patterns follow canonical relations [59]. However, magnetism and configurational entropy counteract chemical driving forces, and in some cases, segregation is prevented even at low temperature. The mechanism of surface segregation in HEMs is usually helpful for understanding of the surficial behaviors under reaction conditions, where most electrocatalysis occurs on the surface [58]. A strong segregation of passivating elements could increase corrosion resistance on the one hand but could compromise the long-term stability of potential high-entropy catalysts on the other hand. For example, Körmann et al. calculated the segregation free energy of FCC HEA (Cr–Mn–Fe–Co–Ni) in vacuum and in the presence of 1 monolayer (ML) of atomic O on the surface. The results indicated that Ni segregated to the surface in vacuum for Cr-, Mn-, and Fe-rich alloys, whereas Cr and Mn segregated to the surface when O was present for every alloy composition. This proposed methodology and approach can be straightforwardly applied to other alloys, such as noble-element HEAs or refractory HEAs.

4.2. Basic descriptors for electrocatalytic reactions

Formation Energy. The formation energy is an important physical property that determines the stability of specific catalysts. More detailed information about DFT calculations of the formation energy is summarized by Persson et al. [57]. The formation energy is first used to estimate the potential of newly designed electrocatalysts [60,61], which is commonly necessary for the prediction of electrocatalysts where the dynamics of oxygen play an important role and need to be studied in detail. When coupled with the obtained experimental results, the preliminary structure of the electrocatalyst can be deduced, and the energetically favorable structure of the catalytic reaction can be screened out [62]. Furthermore, it has also been used to direct the design of next-generation electrocatalysts. For example, Nemani et al. first synthesized two high-entropy MXenes (TiVNbMoC₃T_x and TiVCrMoC₃T_x) and the corresponding precursor high-entropy MAX phases (TiVNbMoAlC₃ and TiVCrMoAlC₃) [63]. They used DFT tools to calculate the formation energies and estimate the formability of the obtained high-entropy MAX phases. It was found that the huge compositional space combined with the potential of entropic stabilization led to various applied directions and challenges.

Adsorption Energy. According to the generally accepted Sabatier principle [64], the appropriate adsorption strength of bonding adsorbates endows superior electrocatalytic activities. The weaker binding interaction may not promise the propitious electrochemical reaction. On the other hand, adsorption that is too strong increases the difficulties in the conversion/desorption of intermediates and final products. A negative adsorption free energy indicates thermally favorable adsorption, and a positive adsorption free energy suggests an energetically favorable desorption process. Generally, the range of suitable adsorption energy varies for different electrocatalytic reactions, and the intermediate free energy plays a key role in the various catalytic reactions [29]. Jan Rossmeisl et al. coupled DFT calculations with supervised machine learning [65]. They successfully predicted the adsorption energies of CO and H on all surface sites over the (111) surfaces of high-entropy alloys (CoCuGaNiZn and AgAuCuPdPt) [65]. Figure 3(a) illustrates that the adjustment of compositions in HEAs could tune the distribution of adsorption energies, which could further change the selectivity and activity of HEAs.

Figure 3. (a) Adsorption energy distributions. Reprinted with permission from Ref. [65] Copyright 2020, American Chemical Society. (b) Scaling relations on HEAs. Reprinted with permission from Ref. [70] Copyright 2020, Elsevier. (c) The heatmap for comparison of the strain effect on d-band center depth variations for the three low-index surfaces. Reprinted with permission from Ref. [73] Copyright 2020, Wiley-VCH GmbH. (d) Calculated charge-density difference of the P1 site for Co_0.6(VMnNiZn)_0.4PS₃. The red and blue regions refer to electron accumulation and depletion, respectively. Reprinted with permission from Ref. [83] Copyright 2022, American Chemical Society.

Scaling Relations. According to previous reports, the scaling relation of adsorption energies between similar adsorbates has already been widely investigated [66,67], and it could act as a representative descriptor for the correlation between different reaction species. The equation is as follows [68]: $E_A=\text{a}{E}_B+b$. Here, E_A and E_B are the corresponding adsorption energies of species A and B, respectively; a and b represent the slope and intercept of the linear relation. Complex electrocatalytic reactions could build a close correlation of the adsorption energies with one/two important intermediates. Meanwhile, the scaling relations reduce the degrees of freedom during the electrocatalytic process and infer the mechanism of relationships between adsorbates and the surface [69]. As an example, Jan Rossmeisl et al. explored the scaling relations between the adsorption energies of key intermediates (*OH and O*) over IrPdPtRhRu (111) surfaces (Figure 3(b)) [70]. The energy of *OH has been demonstrated to act as an effective descriptor with molecular power for the ORR. On the other hand, the difference in adsorption energies between O* and *OH describes the OER.

d-Band Model. The d-band model was first established and complemented by Nørskov and Hammer [71,72], which has been successfully used in understanding the mechanism of the electroactivity of various catalysts. To evaluate the effect of the d-band on the adsorption interactions, the corresponding d-band center (${\epsilon }_d$) is described as the local average of the d-electron energies: ${\epsilon }_d=\frac{{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\rho (x)\text{d}x}{{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\rho (x)\text{d}x}$. Here, x represents the energy level, and $\rho (x)$ is the density of states of the corresponding d-orbitals [71]. The ${\epsilon }_d$ center was described as a scaling relationship for the adsorption energy of reactants/intermediates. In general, a positive shift in ${\epsilon }_d$ results in an enhanced metal—site-adsorbate interaction. On the other hand, a negative energy shift leads to weaker binding. Due to the various multi-metallic sites in HEMs, ${\epsilon }_d$ exhibits a huge adjustable space and thus is beneficial for controlling the locations of ${\epsilon }_d$—active sites. Under these circumstances, it could easily regulate the binding energy of intermediates and further improve the electivity for the reaction products. The d-band center tends to be adjusted by introducing strains caused by lattice distortions. Huang et al. investigated the strain effects of five kinds of transition metals (Fe, Co, Ni, Pd, and Pt) in facet engineering, surface atom density, and ${\epsilon }_d$ [73]. They introduced a mechanism of Pt, Pd, Fe, Co, and Ni within three low index surfaces (111), (110) and (100) over HEAs. By convention, when the strain property changes from compression to tension, the center position of the d-band is expected to move up linearly. In this work, the variation tendency of the surfaces under different strain levels (ranging from −4% to 8%) is summarized. The ${\epsilon }_d$ levels of each metal surface are given by the scale bars (Figure 3(c)). As the color changes from blue to red, the represented level increases monotonically. Hence, the relative variation trend of the data of each element/surface under a given strain sequence is provided.

Orbital Occupation. Although d-band theory has been successfully used, it is not suitable for TM-based materials with split d-band orbitals [74], especially for metals coordinated by six identical ligands in an octahedral shape with O_h symmetry. The d orbitals of the center TM atoms can be split into low-energy t_2g orbitals and high-energy antibonding e_g orbitals [75,76]. The filling state of the e_g orbitals could be effectively used to describe the activity of the TM catalyst, which acts as a key factor for altering the binding strength of adsorbed/desorbed intermediates over the surface of catalysts. Due to this efficient descriptor, huge efforts have been devoted to regulating the catalytic activity of various kinds of metal oxides [77,78]

Charge distribution. Due to the different work functions of the active metal site, the redistribution of charge over the catalyst surface results in alternating accumulation and scarcity of electrons across the entire HEMs. The specific charge distribution could also regulate the adsorption of key intermediates and improve the selectivity of the reaction product. For alloys (especially HEAs), electron transfer over different metal surface atoms easily occurs due to the various work functions, which could result in an obvious charge redistribution on the surface [79]. The alteration of the local charge density is favorable to the generation of adsorptive active sites and reactive reactant molecules [80], promoting electrochemical kinetics. Moreover, the charge redistribution could also affect the d-band centers, which regulate the adsorption methods of reagents and the corresponding electrocatalytic activity and selectivity [80]. Due to the nature of large heterogeneous metal components, HEMs possibly endow significant redistribution of surface charges [81,82]. Thus, for HEMs, the charge density of various surface metallic active sites is different from that of the coordinated atoms, which provides the effective active center for chemical transformation. Wang et al. successfully introduced a typical 2D high-entropy nanosheet catalyst, Co_0.6(VMnNiZn)_0.4PS₃, with a high concentration of active sites [83], which exhibited superior HER performance to most reported HER catalysts. With the help of DFT calculations, the effect of the charge redistribution on hydrogen adsorption was demonstrated by the charge density differences of Co_0.6(VMnNiZn)_0.4PS₃ and CoPS₃ structures around the adsorbed P1 sites (Figure 3(d)). The different states of electron accumulation/depletion over P sites and S sites suggested fast electron transfer between S and P atoms, which was in favor of the optimization of the P–H bond and thus reduced ΔG_H*.

5. Applications of HEMs simulations in electrocatalysis

HEMs exhibit enhanced catalytic activity, selectivity, and stability relative to their mono-metallic and phase-separated counterparts. Serving as a discovery platform for constructing highly efficient catalysts with unexpected performance, HEMs have attracted increasing attention.

5.1. Electrocatalytic HER

During the past decades, the excessive consumption of fossil energy has caused the urgent development of clean and green fuel. Hydrogen has been regarded as a potential next-generation energy source due to its high gravimetric energy density and zero-carbon emissions [84,85]. Electrochemical water splitting is a particular method for the large-scale production of hydrogen fuel [86], which consists of two half-reactions. The HER is one of the semi-reactions during the electrochemical water splitting process, including the two-electron transfer reaction [84]. The reaction processes of the acidic/alkaline HER are different [87]. However, for the two kinds of HER, the Volmer reaction is usually the first step, which involves one-electron transfer to the active surface sites of electrocatalysts, followed by the interaction with H⁺ for the formation of adsorbed H*. The Volmer reaction is usually regarded as the RDS for the HER [85].

Wang et al. developed a facile magnetron sputtering method to fabricate a good-crystallinity high entropy alloy (FeCoNiCuPd) thin film catalyst on carbon fiber substrates [88]. The obtained HEA electrocatalyst shows superior alkaline HER and OER performance with ultralow overpotentials as well as Tafel slopes, showing how efficiently the electrode can produce current by applying a potential. DFT simulations were carried out to identify the roles of individual metal sites and their synergistic effects. As shown in Figure 4, the site-dependent d-PDOS of various active metal sites are calculated to verify their roles in the improvement of alkaline HER activity. The electron density difference result demonstrated that Fe sites favor the adsorption of water molecules. In addition, the Gibbs free energy of the two steps was evaluated. It was found that the synergistic effects of numerous metallic active site configurations lowered the reaction barrier to accelerate the HER process.

Figure 4. DFT calculation of FeCoNiCuPd HEA catalyst for alkaline HER. Reprinted with permission from Ref. [88] Copyright 2022, Elsevier.

5.2. Electrocatalytic OER

As another semi-reaction, OER is also important for the overall electrochemical water splitting. The OER process involves a more complex four-electron transfer process, which seriously hampers the rate of water splitting [89,90]. The key point in the development of OER electrocatalysts relies on improving the efficiency of electron transportation over catalytic active sites and advancing the step of adsorption/desorption of the reactive intermediates produced during the OER process [90]. For OER catalysis, the different local atomic environments result in a distribution of binding energies for the catalytic intermediates. Svane et al. outlined a strategy for the theoretical optimization of the composition by using the rutile (110) surface of HEO catalysts as an example [89]. DFT calculations predicted the reaction energies for all possible local atomic environments. The different local environments around the catalytically active sites result in a distribution of adsorption energies for the catalytic intermediates and therefore also in a range of overpotentials, with some sites having lower overpotentials than any of the pure oxides.

Zhang et al. creatively fabricated a novel high entropy (oxy)-hydroxide (CoCuFeAgMoOOH) to achieve a highly efficient OER [91]. First, they performed a detailed DFT calculation to predict the real high-efficiency active sites for OER and thus guide the design of high-entropy electrocatalysts. It was found that effective M-O d–p orbital hybridization could lower the energy barrier for RDS (Figure 5(a)). The ubiquitous multi-metallic nature is also beneficial for the adjustment of the electron structure to accelerate adsorption/desorption of key intermediates. Overall, this work provides a deep understanding of the influences caused by the poly-metallic material system on each metal reactive site, which is always related to the mature AEM.

Figure 5. (a) Free energy landscape and TDOS/PDOS plots for high entropy (oxy)hydroxides. Reprinted with permission from Ref. [91] Copyright 2022, Wiley-VCH GmbH. (b) DFT calculation of HEO catalyst for alkaline OER. Reprinted with permission from Ref. [92] Copyright 2022, Wiley-VCH GmbH.

Tang et al. successfully synthesized a high entropy perovskite cobaltate with five equimolar metal elements in the B-sites (Mg, Mn, Fe, Co, and Ni), which exhibited superior electrocatalytic activity for OER (with a low overpotential of 320 mV at a current density of 10 mA cm⁻²) compared to commercial IrO₂ [92]. Systematic experiments and DFT calculations are carried out to uncover the underlying catalytic mechanism and the unknown relationship between configuration entropy. As shown in Figure 5(b), the reaction free energies by the AEM and LOM paths were both taken into consideration, which could be helpful for understanding the whole thermodynamics and primary steps of the OER process on HEOs with different configuration entropies. It is found that the configuration entropy could play a key role in improving the inherent activity of the Co sites and improving the OER kinetics. In addition, the enhanced OER performance could be partially caused by the incorporation of OVs and the accompanying formation of high valence metal active sites.

5.3. Electrocatalytic ORR

The ORR is the reverse reaction with respect to OER. The ORR is also a complex four-electron transfer process, thus inevitably resulting in sluggish kinetics [93–96]. The ORR is important for the further development and practical applications of efficient fuel cells and has attracted considerable attention from researchers worldwide. Based on theoretical explanations and predictions, similar to OER, the binding energies of key intermediates could be convincingly used as descriptors for catalytic activity according to the well-established Sabatier principle [64]. As previously mentioned, the adsorption energy of reactive intermediates should be modest (neither too strong nor too weak). For HEMs, the intrinsic single solid solution phase favors a homogeneous distribution of all metal sites, prospectively providing numerous efficient active sites for the ORR, thus further advancing the reaction process. Batchelor et al. reported that the multidimensionality challenge caused by high entropy electrocatalysts can be overcome by a data-driven discovery cycle method [97]. Iteratively refined computational models precisely predicted the activity trends of HEA films. Furthermore, high-throughput characterization datasets are used as input for refinement of the model. The refined model correctly predicts the activity maxima of the exemplary model system Ag–Ir–Pd–Pt–Ru. This method can identify optimal HEA catalysts for electrocatalytic reactions in an unprecedented manner.

Jin et al. developed a facile dealloying method to fabricate various predetermined rugged HEA nanowire catalysts for efficient trifunctional HER, OER, and ORR [97]. It was found that the high entropy design strategy could effectively reduce the use of noble metals and significantly improve the adjustability of specific electronic structures for the ORR. Detailed DFT calculations were performed to investigate the origins of the enhanced ORR kinetics with the introduction of the fifth non-noble transition metal. DFT results indicated that the octahedral site is more energetically favorable and stable for potential substitution than the tetrahedral site. Under these circumstances, the e_g occupancies of Co/Ru atoms are calculated, which significantly simplifies the computational process (Figure 6(a)). In addition, the corresponding PDOS are provided in Figure 6(a). In the NiCo₂O₄ model, the e_g occupancy of Co is 0.71. The obtained e_g occupancies of Co within the doping of heterogeneous atoms (Mo, Cu, V, and Fe) are 1.10, 0.77, 0.79, and 0.83, respectively. This thus caused a decrease in the absolute differences from unity in the sequence of Cu > V > Fe > Mo. For Ru sites, the calculated values for substitution by Mo, Cu, V, and Fe doping are 0.98, 0.86, 0.76, and 0.76, respectively. The e_g occupancy trends of Co/Ru sites well supported the experimental OER/ORR activity trend.

Figure 6. (a) DFT calculation of HEAs for ORR. Reprinted with permission from Ref. [97] Copyright 2020, American Chemical Society. (b) Parameterization of the surface configurations and the activities of re-engineered compositions of the HEA IrPdPtRhRu. Reprinted with permission from Ref. [98] Copyright 2018, Elsevier.

Batchelor et al. found that the adsorption energies of HEAs exhibited a near-continuous distribution, ultimately promoting HEAs as a well-posed platform for next-generation alloys for various catalytic applications [98]. The authors built a systematic database of adsorption energy as well as the construction function models by detailed DFT calculations of the adsorption energies of OH* and O* over the surface of HEAs (Figure 6(b)). Benefiting from the easy adjustment of the metal compositions in HEAs, it is generally achievable to design and fabricate outstanding ORR electrocatalysts with the help of a volcano map construction of highly active sites on the surface. Optimization of the HEA composition results in Ir_10.2Pd_32.0Pt_9.30Rh_19.6Ru_28.9 HEAs. For example, Ir_10.2Pd_32.0Pt_9.30Rh_19.6Ru_28.9 HEAs, which possessed the optimized adsorption energies on the volcanic map, showed much lower overpotentials than commercial Pt/C catalysts.

5.4. Electrocatalytic ADR

Ammonia (NH₃) is a significant raw material for the industrial production of fertilizers and acts as an outstanding hydrogen storage carrier [99,100]. The catalytic decomposition of ammonia into hydrogen has attracted increasing attention due to its ability to liquify easily at a mild pressure of approximately 8 bar at room temperature. Xie et al. reported a novel HEA catalyst consisting of earth-abundant elements to realize efficient ammonia decomposition [99]. The obtained CoMoFeNiCu nanoparticles exhibited a single solid-solution phase with a controllable Co/Mo atomic ratio. These HEA electrocatalysts exhibited substantially improved catalytic activities and stability, even superior to those of Ru catalysts. According to the DFT calculation results of the adsorption properties of various metallic sites, Co−Mo-based alloy nanocrystals are also demonstrated to be potential candidates for ammonia decomposition (Figure 7(a)). The increase in compositional complexity may be a promising method to solve this problem due to the single-phase stabilization caused by increasing the configurational entropy. Saidi et al. developed a novel HEA catalyst, Co₂₅Mo₄₅Fe₁₀Ni₁₀Cu₁₀, for ammonia decomposition [100], which rivals that of state-of-the-art, but prohibitively expensive, ruthenium catalysts. They built a model to rapidly compute the adsorption energy of H, N, and NHx (x = 1, 2, 3) species on HEA surfaces with varied alloy compositions and atomic arrangements with the help of first-principles calculations in conjunction with data analytics and machine learning. The results indicated that the Co/Mo ratio (25/45) identified experimentally as the most active composition for ammonia decomposition increases the likelihood that the surface adsorbs nitrogen equivalently to that of ruthenium while at the same time interacting moderately strongly with intermediates.

Figure 7. (a) Schematic illustration of the rate-limiting factors in NH₃ decomposition. Reprinted with permission from Ref. [99] Copyright 2019, Springer Nature, Creative Commons Attribution License (CC BY). (b) Prediction of cathodic NRR activities for HEOs. Reprinted with permission from Ref. [104] Copyright 2022, Wiley-VCH GmbH.

5.5. Electrocatalytic NRR

As a typical technique for nitrogen reduction, the Haber–Bosch method always requires a high energy demand with inevitably numerous CO₂ emissions. Under such circumstances, a clean electrochemical NRR was developed [101]. Many researchers have relied on DFT tools to verify the electrochemical NRR mechanism. It was concluded that the N₂ electrolysis process is hampered by the ultra-high overpotentials caused by the two sluggish kinetics semi-reactions (OER and NRR), which require an external impressed voltage to actuate the multielectron-transfer reaction efficiently. To date, much effort has been devoted to developing outstanding electrocatalysts to improve NRR kinetics due to enhanced electron transfer and flexible coordination [101,102]. These various redox capabilities are beneficial to improve the adsorption and activation of reactants and intermediates, leading to superior electrocatalytic performance, especially multielectron-transfer processes [103]. These cation compounds can generate a synergistic effect to regulate the adsorption energy of the intermediate to improve the catalytic activity.

With the help of DFT calculations (density of state and difference charge density of HEO models), Sun et al. predicted that HEOs are promising for effectively advancing the NRR and OER progress (Figure 7(b)) [104]. It was found that Fe is the most effective active site of the NRR. The high-density valence electron orbitals also demonstrated that HEOs exhibited fast electron transfer, which could promote N₂ activation. Then, they successfully proposed a facile method to fabricate hollow HEO nanospheres with surface ultrathin nanosheets, which showed superior NRR performance compared to commercial catalysts.

5.6. Electrocatalytic CO₂RR

Converting carbon dioxide into valuable chemical fuels not only provides a sustainable way to produce clean and renewable energy but also lowers carbon dioxide emissions [105]. However, it depends on highly active and selective catalysts. The electrocatalytic reduction of CO₂ is a potential strategy that generates various products during complex CO₂RR pathways.

Cavin et al. reported 2D high-entropy TMDC alloys with four/five transition metals via systematic prediction, synthesis, and multiscale characterization [106]. Among the obtained samples, the quinary alloy catalyst (MoWVNbTa)S₂ with the higher configurational entropy exhibited outstanding performance. The experimental test displayed a superior current density of 0.51 A cm⁻² and a high turnover frequency of 58.3 s⁻¹ at ≈ −0.8 V (Vs RHE). As shown in Figure 8(a), DFT simulations indicated that the improved CO₂RR activity was related to a multisite synergistic effect and atomic-scale disorder, which activated the RDS (CO desorption). In addition, 2D high-entropy TMDC alloys serve as a materials platform to guide the design of outstanding catalysts for various applications.

Figure 8. (a) Theoretical catalysis results for the HEA, select pure TMDCs and silver. Reprinted with permission from Ref. [106] Copyright 2021, Wiley-VCH GmbH. (b) Comparison of adsorption characteristics of HEA. Reprinted with permission from Ref. [107] Copyright 2021, Springer Nature, Creative Commons Attribution License (CC BY).

Mori et al. proposed that TiO₂ could act as an effective substrate to support the synthesis of CoNiCuRuPd HEA nanoparticles at a relatively low temperature (400°C) [107]. The frequency of CO adsorbed on FCC sites over HEA (111) facets was calculated (Figure 8(b)), with an average ν_CO of 2079 cm⁻¹. The CO H adsorption energies (E_ad) on FCC, FCC hollow and HCP hollow configurations were calculated in detail on the HEA/pure metal surface. The average E_ad values for CO and H adsorption on the HEA model were −37.5 and −50.3 kcal/mol, respectively. However, the corresponding values for the Pd model are −26.2 kcal/mol and −54.2 kcal/mol, respectively. This result suggested that the interaction between the CO/HEA interface was stronger than that of the CO/Pd interface, which indicated the easier formation of CH₄ and CO on the HEA surface. Furthermore, the obtained results supplemented solid evidence that the unique cocktail effect of high entropy catalysts related to the synergistic effect of the polymetallic sites could effectively tune the electronic properties for the improvement of catalytic activity.

5.7. Other electrocatalytic applications

In addition to the abovementioned electrocatalytic reactions, HEMs could also exhibit the potential for other complicated electrocatalytic reactions that involve more electrons/protons transfer steps, such as the AOR [108], GOR [109], UOR [110], and decomposition of organic pollutants [111]. The tools of DFT calculations and the combination of multimethod theoretical simulations were proven to be effective for understanding the underlying mechanism of these intricate processes.

Chen et al. developed a facile method to synthesize HCP PtRhBiSnSb HEI, which could effectively accelerate the AOR progress in the alkaline environment [108]. Benefiting from the ensemble effect of those highly active elements, the PtRhBiSnSb HEI samples showed superior mass activities of 19.529, 15.558, and 7.535 A mg⁻¹_Pt + Rh for the electrooxidation of methanol, ethanol, and glycerol, respectively, outperforming most reported multifunctional electrocatalysts for the AOR. DFT calculations indicated that the addition of Rh could optimize the intrinsic electronic structures and enhance the electron transfer efficiency of the HEI samples, which could boost the oxidation capability. The robust electronic structures achieved by the synergistic protections from Bi, Sn, and Sb sites lead to stronger selectivity and durability for the CO₂ reaction pathways of alcohol oxidations.

Fan et al. successfully fabricated well-designed CoNiCuMnMo HEA nanoparticles and the as-prepared electrodes exhibited outstanding GOR electrocatalytic performance with low overpotential and high selectivity toward mate products [109]. The self-developed machine learning-based Monte Carlo simulation and DFT calculation were synergistically employed to simulate the surficial atomic configurations and further explore the underlying mechanism. The scaling relationship between different intermediates’ adsorption energies and the adsorption Gibbs energies of C₃H₆O₃ was established and acted as the descriptor for the electrocatalytic performance. The Mo sites coordinated by Mn, Mo, and Ni showed superior catalytic activities, supported by the theoretical volcano plot of the GOR.

6. Outlooks

First principles modeling and simulations have been widely confirmed to be helpful for understanding both the atomic and electronic structures of HEMs. The relationship between structure and remarkable electrocatalytic performance has been studied extensively with the aid of numerous atomistic simulation methods. In the current review, we analyze recent advances in first principles modeling and simulations of HEMs. We review the advancements in modeling and simulations of HEMs for electrocatalytic applications, including the HER, OER, ORR, ADR, NRR and CO₂RR. However, simulation research on HEMs continues to grow rapidly. We conclude by summarizing critical future prospects at the current time as follows.

1. Enhanced ensemble effect in HEMs: Electrocatalysts with specific surface compositions can provide highly catalytically active sites because of the ensemble effect [112,113]. Due to the heterogeneous atomic distribution, it is difficult to differentiate the effect on surfaces experimentally. There is an enhanced ensemble effect in the case of HEMs due to their multi-component surfaces [8]. There are many possible atom arrangements on the surface of HEMs, inducing different adsorption modes of reactants and intermediates. The march of Moore’s law has led to unprecedented computational power that enables such high-performance modeling and simulation [114]. Identifying the active center accurately is important in atomistic simulations of HEMs for electrocatalytic applications.

2. High-throughput simulations of HEMs: Effective design of HEMs with remarkable properties is a prerequisite for applications. However, their vast space of chemical composition makes the design an overwhelming task. During the design, most of the current HEMs simulations have been case-by-case constructions, i.e. researchers model the structure and perform atomistic simulations only for their specific electrocatalytic application. Using high-throughput DFT simulations [115,116] to screen for new HEMs and perform corresponding experimental studies can show a hopeful opportunity for designing advanced electrocatalysts. High-throughput simulations of HEMs require computations on a vast number of high-entropy structures, and it requires developing new DFT calculation methods and data management technologies.

3. Machine learning based on DFT simulations of HEMs: With the mentioned high-throughput approach, up to tens of thousands of simulation results of HEMs can be obtained. Among the developed tools for extracting knowledge, machine learning has been coined in the field of materials informatics. This makes it possible to handle the data analysis of atomistic simulations in a scalable way. In so doing, it overcomes several inherent disadvantages of simple DFT methods, e.g. low efficiency, long development cycle, and high cost [117,118]. After machine learning training, the prediction models can become powerful surrogates of the atomistic simulations and show remarkable advantages over case-by-case simulations of high-entropy electrocatalysts. It can be used to obtain stoichiometry effectively and phase-structure–property relationships efficiently for material-property database establishment [117].

4. Combination of standard theoretical calculations: taking HEOs as an example, although empirical predictors, e.g. mean ionic radius and size mismatch between the involved cations, are widely used in HEO modeling [119], there are still several novel combinations of standard theoretical calculations for HEO modeling. For example, the degree of cation order/disorder in HEOs can be evaluated using ab initio molecular dynamics using SQS modeling and Reverse and Metropolis Monte Carlo simulations [120]. Short-range order has been confirmed to have a significant impact on the mechanical properties of HEMs [121,122]. The relationship between structure and electrocatalytic performance in HEMs should be deeply revealed and understood using different combinations of standard theoretical calculations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The work was supported by Natural Science Foundation of Jiangsu Province, China [BK20220428], National Natural Science Foundation of China [52171110], Natural Science Foundation of the Jiangsu Higher Education Institutions of China [21KJB430014]. Also, support was acknowledged from the European Union Horizon 2020 Research and Innovation Program [857470], the European Regional Development Fund via the Foundation for Polish Science, International Research Agenda PLUS program [MAB PLUS/2018/8], and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) [NRF-2021R1A2C3006662, NRF-2022R1A5A1030054].