A Review of Recent Developments in Molecular Dynamics Simulations of the Photoelectrochemical Water Splitting Process

Abstract

In this review, we provide a short overview of the Molecular Dynamics (MD) method and how it can be used to model the water splitting process in photoelectrochemical hydrogen production. We cover classical non-reactive and reactive MD techniques as well as multiscale extensions combining classical MD with quantum chemical and continuum methods. Selected examples of MD investigations of various aqueous semiconductor interfaces with a special focus on TiO₂ are discussed. Finally, we identify gaps in the current state-of-the-art where further developments will be needed for better utilization of MD techniques in the field of water splitting.

1. Introduction

It is needless to say that global climate change, due to the emission of vast amounts of CO₂ and other greenhouse gases from burning fossil fuels, is one of the main challenges mankind is facing today. Therefore, in order to reduce such emissions into the atmosphere, alternative energy sources are indispensable. Hence, a transition to renewable energy sources has gained strong momentum within the last decade or so. Various alternative renewable energy sources have been introduced in the past including solar, geothermal, wind and biomass. Although these renewable resources make a significant contribution to the modern world, they fluctuate with time and strongly depend on geographical locations. Over the past decade, investments in renewable energy have increased tremendously. Among the renewable sources, sunlight has a pivotal role since it is by far the world’s most abundant source of renewable energy. However, due to its unequal temporal and geographical availability, the storage of energy converted from solar radiation will become more and more important with an increasing share of solar energy in the world’s overall energy consumption. Here, the direct synthesis of hydrogen by photoelectrochemical water splitting is considered a promising technological pathway that overcomes the storage problem of our most abundant but strongly fluctuating energy source [1].

This review article was written as part of the European Cooperation in Science and Technology initiative “COST Action 18234—Computational materials sciences for efficient water splitting with nanocrystals from abundant elements,” funded by COST [2]. Up to now, different modeling techniques from ab initio methods up to the continuum scale are mostly employed independently in the field of photoelectrochemical water splitting. In order to overcome the limitations of the single modeling approaches, the COST action 18234 intends to combine different modeling scales and methods with the ultimate goal of generating an integrated multi-scale modeling environment for water splitting materials. The need to combine the competences from different areas and institutions to accelerate the development of viable and cost-effective solar water splitting solutions was also recognized by other research initiatives such as the HydroGEN consortium in the United States, which brings together various capabilities from modeling, synthesis, characterization and system integration [3]. Here, we review state-of-the-art applications of molecular dynamics simulations (MD) in the field with the aim of identifying the integration potential with other simulation methods such as ab initio and continuum calculations.

Water Splitting Application Background

Since the seminal article of Fujishima and Honda [4] reporting unbiased dissociation of water into hydrogen and oxygen with the help of an illuminated TiO₂ anode and a Pt cathode, photoelectrochemical (PEC) water splitting has become a main research field in photoelectrochemistry due to its great promise of the eco-friendly and renewable production of hydrogen as a practically inexhaustible source of energy [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. The solar splitting of water molecules into oxygen and hydrogen is based on a rather simple principle: (1) Photons excite electron-hole pairs in a semiconducting (SC) absorber material; (2) the excited hole in the valence band drives the oxidation of water to form O₂ (oxygen evolution reaction–OER) once it has reached the anode surface; (3) the excited electron in the conduction band reduces water to H₂ at the cathode (hydrogen evolution reaction–HER). Hence, the water splitting process can be summarized as follows [4]:

$$\mathrm{SC}+4 \mathrm{h}\mathsf{\nu }\to 4 {{\mathrm{e}}^{-}}_{\mathrm{exc}}+4 {\mathrm{p}}^{+}{\mathrm{exc}}_{.}$$

(1)

Oxidation (OER) at the anode:

$$4 {{\mathrm{p}}^{+}}_{\mathrm{exc}}+2{\mathrm{H}}_2\mathrm{O} \to {\mathrm{O}}_2+4{\mathrm{H}}^{+}$$

(2)

Reduction (HER) at the cathode:

$$4 {{\mathrm{e}}^{-}}_{\mathrm{exc}}+4{\mathrm{H}}^{+} \to 2{\mathrm{H}}_2$$

(3)

The net effect is the dissociation of water into the gases of its constituent elements:

$$2 {\mathrm{H}}_2\mathrm{O}+4 \mathrm{h}\mathsf{\nu }\to 2 {\mathrm{H}}_2+{\mathrm{O}}_2.$$

(4)

Thermodynamically, the overall water splitting reaction (4) requires photons with energies of at least 1.23 eV under standard conditions to run downhill in free energy. In practice, however, the oxidation and reduction half reactions (2) and (3) exhibit considerable overpotentials due to several factors such as energetically unfavorable intermediate states of the catalytic reactions and resistances within the PEC system at the electrolyte interfaces, solid-state junctions and bulk losses. Therefore, semiconductor systems with sufficiently large band gaps of at least 1.6 eV are required [1]. On the other hand, band gaps far beyond 2 eV only allow the absorption of a small part of the solar spectrum and thus also result in low so-called solar energy to hydrogen (STH) efficiencies. Moreover, the edges of the conduction and valence bands must be properly aligned with the water redox potentials so that the transfer of the excited charge carriers can occur at the solid/liquid electrolyte, see Figure 1.

Although PEC water splitting with the help of semiconducting materials is based on a very simple principle, and despite the enormous amount of research activity, so far no PEC system has been devised that allows for economically viable and scalable solar water splitting. The main reason for this is the manifold requirements to be met by the employed materials at the same time:

1.

High photoabsorption efficiency and efficient separation of excited charge carriers;

2.

High catalytic activity with fast reaction kinetics and low overpotentials;

3.

Proper alignment of conduction and valence band edges with respect to the water redox potentials;

4.

Corrosion resistance and long lifetime without performance decay;

5.

Earth abundance and low production costs.

An argument by Idriss [21] illustrates that it is not trivial to find a stable, earth-abundant semiconductor with the proper bandgap because, in nature, such a material prevented the formation of water on Earth at all.

The present paper reviews current approaches in the application of molecular dynamics techniques for water splitting. The structure of the paper is as follows: In Section 2, “Molecular dynamics,” we provide an overview of the molecular dynamics field, focusing on molecular dynamics algorithms including temperature coupling algorithms such as Berendsen temperature coupling, velocity-rescaling temperature coupling, and Andersen thermostat or Nose-Hoover thermostat. Next, we describe reactive force field simulations, which are necessary to describe the chemical reactions during the water splitting process. The focus is put on the popular and widespread ReaxFF models. This section ends with a review of how molecular dynamics simulations can be combined with other methods used in computational chemistry, such as quantum mechanical ab initio, Monte Carlo and continuum simulations. Section 3, “Review utilization of MD to water splitting,” is a review of how molecular dynamics simulations are used to improve artificial photosynthesis processes, highlighting characterization methods for the flow of reactants, the study of water splitting photocatalysts such as Pt/TiO₂ and solid–water interphases. Section 4, “Conclusion,” summarizes the main results and discusses the existing research gaps in the application of molecular dynamics to water splitting.

2. Molecular Dynamics

In this section, the Molecular Dynamics method (MD) is discussed. In Section 2.1 the key steps of the MD method are briefly explained. In Section 2.2, reactive MD (ReaxFF) is explained and discussed. Finally, multiple multiscale algorithms, which include MD, are discussed in Section 2.3.

2.1. A Brief Discussion about the Method—Formulation, Algorithm

Molecular Dynamics (MD) belongs to the realm of computational chemistry and is a deterministic approach used to simulate and analyze the dynamics of molecular systems. It is based on the Newtonian equations of motion. The particles interact through potentials—bonded and nonbonded—for a given number of simulation steps. In each simulation step, forces, velocities and positions of particles are updated. Thermodynamic quantities, such as temperature, pressure and number of particles, are controlled using statistical ensembles. Molecular dynamics is commonly used in material science, biophysics, chemistry, and so forth, for systems of 10⁵–10⁶ atoms/particles, and to cover time scales spanning nanoseconds to microseconds. On big supercomputers, however, these numbers can increase (e.g., systems of 10⁹ or 10¹² particles have been simulated [22,23,24,25,26]). The use of approaches, such as machine learning molecular dynamics (e.g., [27]) and coarse grain MD, can allow one to reach further extended scales (e.g., time scales beyond 1 microsecond can be reached by coarse grained MD).

There are many computational ingredients for the molecular dynamics simulations algorithm. However, at a high level we can describe a global molecular dynamics algorithm as a sequence of the following steps:

• Firstly, the initial conditions must be set and several inputs, such as the potential interaction V that is computed as a function of atom positions, the positions r of all atoms in the system and the velocities v of all atoms in the system, are required.
• In many subsequent computational steps, the following quantities are calculated:

  ◦

  new potentials—due to atomistic interactions;

  ◦

  corresponding forces.

• After that, an update configuration algorithm based on Newton’s equations of motion follows and new atomistic velocities, positions, system temperature and some other characteristics of interest are computed.
• After a given number of steps, configurations, forces, potential and other data of interest are saved repeatedly for post-processing.

Without being exhaustive, we will describe some important ingredients related to molecular dynamics algorithms. For further study, we can refer, for example, to the GROMACS user manual written by Spoel et al. [28].

2.1.1. Initial Step

First, the parameters of the force field must be loaded. A force field in molecular modeling consists of parameters that are used for calculating the potential energy of a system of atoms. They are deducted from quantum mechanical calculations, chemical/physical experiments or both. Next to the parameters of force fields, the size of the simulation box, and the coordinates and velocities of all particles are required. The size and shape of the box are determined by the corresponding basis vectors, and periodic boundary conditions along the different directions might be imposed.

2.1.2. Computation of Interaction, Forces and Update Configuration Algorithm

The potential energy of the interacting atoms and the corresponding forces are calculated. These can be conducted with a variety of methods, for example:

• (Reactive) force fields, in which case we deal with a classical molecular dynamics simulation. The use of a (reactive) force field will be discussed extensively in the next section;
• Quantum mechanical model derived from the Schrödinger equation, in which case, an ab-initio molecular dynamics simulation is performed;
• Using a combined quantum-classical approach, in which case a quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulation is performed.

We will first discuss the use of classical force fields to calculate the potential energy surfaces. Depending on the interactions to be described, various types of potentials are used within the force fields. Several examples of potentials include the Lennard–Jones potential [29,30], which can be used to describe the dispersion interactions; the Coulomb potential [31], which is used to describe the ionic interactions; the Tersoff potential [32], a bond order potential, which is used to describe covalent and variable bonding for elements like carbon, silicon, germanium, and so forth; Sutton–Chen potentials [33], a many-body potential used to describe metals (e.g., Cu, Ag, Au, Pt); and TIPxP potentials, used to describe water [34]. Often, a mixture of these potentials is used, for example, Coulomb-Lennard–Jones potentials, which are able to describe both ionic and dispersion interactions. Sometimes, in place of interaction potentials, constraints can be used.

The Lennard–Jones potential V_LJ between two non-bonded atoms can be written as follows:

$$V_{LJ}\left(r_{ij}\right)=\frac{C_{ij}{}^{12}}{r_{ij}{}^{12}}-\frac{C_{ij}{}^6}{r_{ij}{}^6},$$

(5)

where C_ij is a constant, and $r_{ij}$ is the distance vector between particles i and j. This potential is a difference between a repulsive term (due to electron–electron interactions occurring at short distances) and an attractive term (due to the long-range van der Waals interactions/dispersion).

Using the resulting force given by a force field (additional popular models are UFF [35], Dreiding [36], MMx [37], AMBER [38], CHARMM [39], GROMOS [40], OPLS [41], COMPASS [42], ReaxFF [43,44]), the update of positions and velocities can be computed through various algorithms: for example, the so-called leap-frog algorithm or Velocity Verlet, which integrate the Newtonian equation of motion. The Velocity Verlet integrator is preferable when extremely accurate integration with temperature and/or pressure coupling is required.

The leap-frog algorithm uses the positions $r$ at time $t$ and velocities $v$ at time $t$-dt/2 to update the positions at time $t+dt$ and the velocities at time $t+dt/2$ using the forces $F\left(t\right)$ that are determined by the positions and force field at time $t$. Explicitly, the time integration is carried out with the following relations:

$$v\left(t+1/2dt\right)=v\left(t-1/2dt\right)+dt\times F\left(t\right)/m$$

(6)

and

$$r\left(t+dt\right)=r\left(t\right)+dt\times v\left(t-1/2 dt\right).$$

(7)

The algorithm produces trajectories that are similar to the Verlet algorithm, whose position update relation is:

$$r\left(t+dt\right)=2r\left(t\right)-v\left(t-dt\right)\times dt+\frac{1}{2}\times \frac{F\left(t\right)}{m}\times d{t}^2+O\left(d{t}^4\right).$$

(8)

In the formulae above, $dt$ is the time step. This algorithm is time-reversible and of third order in $r$. In addition to the updates of velocities and positions, various algorithms for pressure or temperature coupling are used as explained in the following.

Many of the molecular dynamics simulations lead to the creation of trajectories within the NVE (constant number, constant volume, constant energy) thermodynamic ensemble, also called the microcanonical ensemble; or are conducted within the NVT (constant number, constant volume, constant temperature) ensemble, also named the canonical ensemble. To keep the temperature of the NVT ensemble at the desired reference temperature of the molecular system, temperature coupling algorithms are used. Among them, we can mention Berendsen temperature coupling [45], Velocity-rescaling temperature coupling [45], the Andersen thermostat [46] and the Nose-Hoover thermostat [47]. Below, as an example, we will present the Berendsen temperature coupling algorithm.

The Berendsen algorithm is based on a weak coupling algorithm with first-order kinetics. The coupling is made to an external bath at the desired reference temperature. The atomistic system temperature is corrected according to the following formulae:

$$r\frac{dT}{dt}=\frac{T_0-T}{\tau },$$

(9)

which means that the deviation of the temperature decays exponentially with a time constant $\tau$. In this way, the correction can be adapted to the user’s needs: for equilibration purposes, the coupling time can be short (e.g., $\tau =0.01$ ps); in general, however, the equilibrium will be reached in a much longer fashion (e.g., $\tau =0.5$ ps), in which case it hardly influences the dynamics of the atomistic system.

After updating the configuration, different quantities of interest are computed, such as temperature, pressure, and so forth. For example, the temperature is given by the total kinetic energy of the N-particle system:

$$E_{kin}=\frac{1}{2}{{\displaystyle \sum }}_{i=1}^N{m}_i{v}_i{}^2,$$

(10)

with $v_{i }$ being the initial velocity and $m_i$ being the mass of the particle. From this, the absolute temperature $T$ can be computed using:

$$\frac{1}{2}{N}_{df}kT=E_{kin},$$

(11)

where $k$ is Boltzmann’s constant and $N_{df}$ is the number of degrees of freedom of the molecular system. Pressure can be also controlled within MD simulations. We can calculate the pressure tensor $P$ as the difference between the kinetic energy and the virial $\Xi$. The following formula is used:

$$P=\frac{2}{V}\left(E_{kin}-\Xi \right),$$

(12)

where V is the volume of the computational box. The scalar pressure $P_s$, which can be used for pressure coupling in the case of isotropic systems, is computed as:

$$P_s=\frac{trace\left(P\right)}{3}.$$

(13)

The virial $\Xi$ tensor is defined as:

$$\Xi =-\frac{1}{2}{\Sigma }_{i<j} r_{ij}\otimes F_{ij}.$$

(14)

2.1.3. Output Step

Data of interest, such as a trajectory file, are saved at given intervals. Due to the large volume of data for the trajectory files, it is advisable to not save data at every step, but at certain desired intervals. These files contain frames that can include velocities and/or forces, and positions as well as valuable information about the dimensions of the simulation volume, integration step/time, and so forth. According to the integrator that has been chosen, the interpretation of the time may vary. From these, dataset quantities of interest, such as thermal conductivity, diffusion, and so forth of the simulated system, can be computed. This can happen after the complete simulation or for some characteristics on-the-fly.

2.2. Molecular Dynamics Simulations with a Reactive Force Field

Traditional MD simulations are based on non-reactive force fields. However, in different cases, such as water splitting, reactive force field simulations are needed [48,49,50,51]. There are multiple methods with empirical reactive force fields, such as ReaxFF [43,44], COMB [52,53] and the empirical valence bond (EVB) method [54]. For clarity, we will further limit the explanation to the (often used) ReaxFF method. ReaxFF reactive force fields were first developed by van Duin et al. [43,44] (shown in Figure 2). While most of the classical force fields cannot model chemical reactions because of the need for chemical reactions to break and form bonds, reactive force fields use bond order relations which allows continuous bond dissociation for all orders at the same time. In this sense, reactive force fields have been parameterized and tested for many complex applications and materials, for example, transition–metal–catalyzed nanotube formation, heterogeneous catalysis applications, alkoxysilane gelation, hydrocarbon reactions, batteries, high-energy materials, thermochemical heat storage, and polymers [55,56,57,58,59,60,61,62,63,64,65]. Furthermore, ReaxFF has successfully described water at saturation conditions [66] and has been applied to water splitting applications. For example, comparable results were found for the dissociation of water with different coverages on titania surfaces between ReaxFF and other theoretical and experimental studies [48]; dynamics related to the dissociation of water on aluminum nanoclusters were studied [49]; the molecular structure of water under thermal and electric field effects [50]; and water-electrolyte systems including cations Li⁺, Na⁺, K⁺, and Cs⁺, and anions F⁻, Cl⁻ and I⁻ [51].

Similar to the non-reactive classical force field, the potential used by ReaxFF is a summation of different energy contributions [19]. However, extra energy terms are included to capture the bond formation:

$$E_{system}=E_{coul}+E_{vdW}+E_{bond}+E_{over}+E_{under}+E_{val}+E_{pen}+E_{torsion}+E_{conj}+E_{others}.$$

(15)

The first two energy terms $\left(E_{Coul}, E_{vdW}\right)$ capture the non-bonded Coulomb and van der Waals interactions, respectively. The Coulomb interaction uses a shielded potential including atomic charges from the Electron Equilibration Method (EEM). The van der Waals interaction is described by a Morse-potential [44]. The term $E_{others}$ can be added to include interactions such as H-bonds or London dispersion interactions. The other terms are related to bonded interactions, in which $E_{bond}$ accounts for the σ-, π-, and ππ-bond energy. $E_{over}, E_{under}, E_{val}, E_{pen}$ and $E_{conj}$ are used to correct for over-, and under-coordination, valence and torsion angles, penalty energies, and conjugated systems, respectively [43,44]. The bond order between atoms is directly related to the distance between atoms and is described by:

$$B{O}^{\prime }{}_{ij}=B{O}^{\sigma }{}_{ij}+B{O}^{\pi }{}_{ij}+B{O}^{\pi \pi }{}_{ij}=exp[P_{bo,1}\cdot {(\frac{r_{ij}{}^{\sigma }}{r_0})}^{P_{bo,2}}\left]+exp\right[P_{bo,3}\cdot {(\frac{r_{ij}{}^{\pi }}{r_0})}^{P_{bo,4}}\left]+exp\right[P_{bo,5}\cdot \left(\frac{r_{ij}{}^{\pi \pi }}{r_0}{)}^{P_{bo,6}}\right].$$

(16)

The use of the bond order terms is a fundamental step in the ReaxFF concept and provides a direct relation between the bond order $B{O}^{\prime }{}_{ij}$ and the interatomic distance $r_{ij}$ of atoms $i$ and $j$. Functions that describe all the contributing energy terms, including empirical parameters, need to be trained with a set of known data (for example, from accurate DFT results). This will be explained in detail below.

2.2.1. ReaxFF Training

Reactive force fields (ReaxFF) can describe bond breaking and formation based on atomic parameters that only contain one single atom type for each element (e.g., the description of oxygen is the same for O in H₂O, O₂, and atomic O). To accomplish this, a complex force field is needed with many nonlinear terms and parameters. To obtain these parameters, an advanced and extensive training algorithm is required. Furthermore, the reference dataset should include relevant points such as reaction enthalpies, equations of state, bond dissociation energies, charges, surface energies, angle and torsion deformations, geometries, and so forth.

The data for the training can be obtained from all kinds of references. However, most often they are generated by DFT calculations [67] since these can provide plenty of accurate data. There are multiple optimization techniques to parameterize the ReaxFF force field such that it is able to replicate the reference data [68,69]. One of the proven training algorithms is the Metropolis Monte Carlo (MMC) force field optimizer developed by Iype et al. [69]. This optimization algorithm is based on the simulated annealing Metropolis algorithm devised by Kirkpatrick and others [70,71,72,73], and is a high-dimensional (treats multiple parameters simultaneously) and efficient training method. This allows one to search for the global minimum of the force field and prevents getting stuck in a local minimum. The goal of the MMC optimizer is to minimize the error between the dataset and the predicted results given by a ReaxFF force field. The error of optimization of a force field is computed via [69]:

$$Error\left(FF\right)={{\displaystyle \sum }}_{i=1}^n{\left[\frac{\left(X_i-X_{i,ReaxFF}\right)}{{\sigma }_i}\right]}^2,$$

(17)

which is the sum of the squared difference between the reference data $X_i$, and the calculated ReaxFF data $X_{i,ReaxFF}$; where the training set contains $n$ data points with a corresponding weight $\sigma$. This weight allows one to emphasize a specific part of the dataset. After each MMC iteration, a random set of selected parameters in the force field is modified in a random direction, and a new error is computed. With the difference between the new and old errors $\left(\mathsf{\Delta }Error\right)$, the new force field can be accepted with the MMC probability [69]:

$$P=min\left[1,exp\left(-\beta \mathsf{\Delta }Error\right)\right].$$

(18)

This is the mathematical formula that samples states according to the Boltzmann distribution. The acceptance rule results in the fact that when the new force field has a lower energy than the old force field (thus probably a more accurate ReaxFF), the new force field is accepted and becomes the old force field. When the error of the new force field is larger than that of the old force field, the new force field will only be accepted with a probability related to the size of the $\mathsf{\Delta }Error$ and the inverse of the thermodynamic temperature $\left(\beta =1/kT\right)$ such that the temperature $T$ acts as a control parameter for the acceptable $\mathsf{\Delta }Error$. A high $T$ allows larger error fluctuations, and a low $T$ only allows small $\mathsf{\Delta }Error$ fluctuations. Thus, the temperature $T$ is a measure that dictates the freedom of meandering above the energy surface of the force field. During the optimization, $T$ is lowered to mimic simulated annealing. A schematic view of the MMC force field optimizer algorithm is given in Figure 3 below.

2.2.2. Free Energy Sampling

Molecular dynamics (MD) is a powerful tool for studying the time evolution of many different chemical systems. In several cases, however, chemists are interested in understanding activated processes that lead to a chemical transformation between reactants and products. From a thermodynamic point of view, reactants and products represent free energy minima separated by high barriers. In most cases, these barriers are associated with extremely long time scales that are incomparably larger than those reachable by MD simulations. This fact has limited the application of molecular simulations to problems that are distant from the main interest of chemistry, that is, reactions.

Since the landmark paper of Torrie and Valleau, in which umbrella sampling was introduced, a large variety of methods have been proposed to study rare events within the framework of molecular simulations. Among many, the most relevant methods proposed over the years are local elevation [75], conformational flooding [76], hyperdynamics [77], metadynamics (MetaD) [78,79], and variationally enhanced sampling [80].

All these methods rely on the definition of a reaction coordinate, typically referred to as the collective variable (CV), along which a bias potential is added. The effect of such bias is to enhance the fluctuations within the basins (the free energy minima in which the system is trapped for long periods), accelerating the transitions between reactants and products. This type of sampling has the double advantage that exploring the conformational space allows the determination of the accurate free energy of the process of interest, thus allowing both thermodynamic and kinetic evaluations on the reactive system.

These methods have been largely used in molecular simulations to study a wide range of rare events ranging from protein folding, crystal nucleation and chemical reactions. In particular, and concerning the topic of this review, metadynamics has been used to study water splitting and oxidation by metal complexes in a solution using full ab initio molecular dynamics (see Refs. [81,82]) as molecular models for photosynthesis. Likewise, the adsorption and dissociation of water have been studied on the Rutile TiO₂ surface [83].

2.3. Multiscale MD Modeling—A Short Overview

In what follows, without being exhaustive, we will try to describe some ways in which molecular dynamics simulations can be combined with other methods used in computational chemistry for water splitting such as quantum mechanics, and Monte Carlo and quantum mechanics simulations.

2.3.1. Molecular Dynamics and Quantum Mechanics

There are several approaches that use a combination of molecular dynamics and quantum mechanics. One possibility is a combined molecular dynamics (MD) and quantum mechanics simulation, an approach that is appropriate for systems in which the electronically active part (i.e., the quantum mechanically treated part) of a large system is localized. Such an approach is implemented in different molecular dynamics tools such as Gromacs [28], Q-Chem [84] and NWChem [85], and so forth.

Several challenges appear when using a mixed quantum mechanics/molecular dynamics implementation: for example, setting of the boundary, the treatment of the bonds between the two regions, the nature of the QM/MM coupling, and the calculation of the total energy. The combined simulation is solved at the level of force fields: some particles are treated as having quantum mechanics potentials and others as having molecular dynamics potentials or both. The QM part of the system is described through the use of ab-initio (e.g., HF, CCSD(T), MP2), density functional (DFT), semi-empirical (e.g., AM1, PM6, PM7, DFTB) or ab-initio MD methods. The choice between them is made based on the size of the quantum region to be described and on the specific requirements for the description of the electronic part. Most often, methods such as DFT or DFTB are employed, as they provide a good balance between the size of the systems that can be described and accuracy. The MM part is described by (reactive) force fields like those presented in the previous sections of this work [86].

QM/MM methods have been used to study a variety of research topics, for example, zeolite chemistry [87,88,89], solvation studies [90] and solid electrolyte interphases [91], as well as various biomolecular studies [92,93].

Another possibility for the use of MD to describe phenomena with relevant electronic behavior is through the explicit inclusion of the treatment of the electron degrees of freedom in the classical force fields. Significant steps in this direction have been made through the development of eReaxFF [94], eFF [95,96] and LEWIS [97,98,99] Ref force fields.

Yet another possibility for combining quantum mechanics and molecular dynamics is through the use of quantum mechanics-based methods (most often DFT) in the training of atomistic potentials, which are further used in molecular dynamics simulations as described in detail in the review by Liang et al. [100]

Finally, MD and DFT can be used successively if one wants to first explore the overall potential energy surface of a complex system and later use the local minima found for further re-optimization and analysis at 0 K.

2.3.2. Molecular Dynamics and Monte Carlo

When Monte Carlo (MC) methods are applied to molecular simulations, we are dealing with Monte Carlo molecular modeling. Both the MD and Monte Carlo approaches can perform the same type of molecular simulations. The main difference is that the Monte Carlo method relies on equilibrium statistical mechanics rather than on molecular trajectories. As opposed to generating the dynamics of a system, it generates states in concordance with Boltzmann probabilities. Moreover, it is a subset of the more general Monte Carlo method, which is widely used in statistical physics.

In Monte Carlo, a new state of a system is generated from a previous one based on a Markov chain procedure. Given its stochastic nature, this new state is randomly accepted. Each trial counts as a move. While the freedom in choosing the moves makes this method very adjustable, the avoidance of dynamics restricts the method to studies of static quantities only. For equilibrium to be properly described, these moves must only satisfy a basic condition of balance. When designing new algorithms, however, a more detailed balance and a stronger condition are usually required. The Monte Carlo approach has given rise to a plethora of generalizations, such as the method of the already discussed simulated annealing for optimization, in which an artificial temperature is introduced and then gradually decreased.

2.3.3. Hybrid Molecular Dynamics

Hybrid Molecular Dynamics combines the two methods (MD and MC) and it is already implemented in some simulation packages such as the Abalone II code [101]. The main idea is that there are two steps for this algorithm: one in which the new momentums of the particles are generated (MC steps); the second when the MD step(s) are performed and the new configurations are accepted or rejected according to the MC criterion. Below we describe this method in more detail, taking the water splitting case as an example.

2.3.4. Hybrid Molecular Dynamics—Monte Carlo Simulations (Hybrid MD–MC)

The conditions at the top of the surface can influence the surface processes and reactions in a water splitting system. A combination of MD and MC can be used by performing MD near the surface boundaries for accuracy and MC in the gas bulk because of the low computational cost. In [102], this method is used to characterize the influence of different densities and molecule sizes on the equilibrium properties of the gas. The hybrid MD–MC simulation method is validated by comparing the results with those of pure MD and pure MC simulations. All three methods produced consistent results regarding density and temperature profiles both in the bulk and near the boundaries when hard-sphere interactions were used. When Lennard–Jones potentials were used to accurately model the interactions between the gas and wall molecules instead, the results of pure MD simulations differed significantly from the pure MC simulations near the boundaries. However, the results of the hybrid method compare well with the pure MD results near the wall and with the pure MC and pure MD results in the middle of the channel. The hybrid method also very accurately simulates the interface between the MD and MC simulation domains.

The speedup when using 50% of the domain for MD simulations and 50% for MC simulations is very small compared to pure MD simulation times, but this speedup increases drastically for more realistic situations where the region near the wall is small compared to the bulk region [102].

2.3.5. Reactive Molecular Dynamics—Monte Carlo Simulations (ReaxFF MD–MC)

MC simulations are based on the statistical sampling of an imposed probability distribution given by a specified ensemble. Thereby, contrary to MD, time integration by Newton’s laws is not needed and rare events or high energy barriers can be simulated or avoided without long simulation times. In this sense, numerous ReaxFF MD–MC combinations are made, to exploit the benefit of MC over MD in reactive simulations. For example, ReaxFF was combined with Gibbs ensemble Monte Carlo (ReaxFF–GEMC) to determine the vapor–liquid equilibrium (VLE) of a given force field [66]. Next to that, ReaxFF–grand canonical Monte Carlo (ReaxFF–GCMC) algorithms have been developed and used for monatomic catalytic oxidation, hydrogenation, or carbonation studies [103,104,105]. Related to water, the combination has also been used for gas phase water formation on platinum catalysts [106], and for water absorption in salt hydrates by some of the authors [107].

2.3.6. Molecular Dynamics–State-Space Modelling Simulations (MD–SSM)

SSM is a well-known method in control theory for simulating complex, interdependent systems [108]. The state–space representation is a mathematical model of a physical system as a set of input, output, and state variables related by first-order differential equations. The state variables can be reconstructed from the measured input–output data but are not themselves measured during the experiment. The general state–space representation is given as follows

$$\dot{x}\left(t,p\right)=\frac{dx\left(t,p\right)}{dt}=f\left[x\left(t,p\right),u\left(t\right),t,p\right],$$

(19)

in which $\dot{x}\left(t,p\right)$ represents the vector of the state variables depending on the time $t$ and the vector of unknown parameters, $p$. The vector $u\left(t\right)$ signifies the input variables that can be varied in the simulations. In addition to the differential equations, the state–space description contains the observation function $y\left(t,p\right)$, which denotes the observed quantities, and it is referred to as the model output

$$y\left(t,p\right)=g\left[x\left(t,p\right),p\right]$$

(20)

State–space modeling has been used in a few examples of electrochemical systems to calculate impedance spectra, current-voltage plots and Tafel plots in fuel cells and batteries [109,110,111,112]. In cited references [109,110,111,112], it was shown for the anodic solid oxide fuel cell system Ni, a H₂–H₂O|YSZ (yttria-stabilized zirconia) electrochemical model can be identified with such a state–space representation. If the mass and charge balances are formulated from an electrochemical model, the surface concentrations of the different adsorbed species are the vector of state variables $x\left(t,p\right)$. The concentration of the adsorbed species changes as a function of time $t$ and as a function of diverse other parameters $p,$ such as surface coverage and reaction-rate constants. The applied potential is the input variable $u\left(t\right)$, which can be changed in the experiments, and the Faraday current is the model output $y\left(t,p\right)$.

By implementing the electrochemical model as described in Refs. [109,110,111,112] in MATLAB and SIMULINK, the impedance at the Ni, H₂–H₂O|YSZ interface can be simulated and compared to experimental measurements under the same conditions. Required input values are, at first instance, derived with the help of the literature. Fitting the simulated data to experimental measurements under standard conditions allows the input values to be optimized. The simulated impedance yields values in the same order of magnitude as the experimental ones, and the spectrum exhibits one arc like in the experiment. The deviation between experiment and simulation is mainly attributed to the availability of kinetic input data. Such data can be determined by atomistic modeling, such as DFT. Diffusion effects can also be implemented in the SSM based on kinetic Monte Carlo simulations. In this sense, the SSM approach is a multiscale modeling and simulation approach. This potential, however, has never been fully tapped in the field of fuel-cell systems, as this approach was used for the first time for an electrochemical system. We note that it has never been applied for PEC systems. A detailed review of MD/KMC/SSM approaches can be found in reference [113].

2.3.7. Molecular Dynamics–Continuum Mechanics

Continuum mechanics (CM) is a sub-branch of mechanics that is concerned with the mechanical behavior of materials that are modeled as a continuous mass rather than as a set of discrete particles. There are many ways in which continuum methods can be combined with MD; several hybrid algorithms of this sort have been proposed in the literature [114,115,116,117,118,119,120]:

• Type 1: CM parametrization: run MD simulation to obtain tables of data [121,122], or a reduced model [123] to be included in the continuum simulation such as surface tension in surfactant-laden flows or the contact line dynamics. This kind of simulation includes the class of parameterizing continuum constants using MD, defining non-local viscosity kernels [124], defining slip boundary conditions [125,126,127] or any type of parameterization which allows the continuum model to be run for the entire space with no additional Molecular Dynamics simulation (e.g., in order to obtain surface tension values, you need to define contact angles and create a model for the moving contact line, which also includes fluctuations [128]);
• Type 2: spawn new MD simulations or call them dynamically during a continuum run to obtain parameters that are transferred to the continuum run, including the effect of complex molecules on viscosity [129,130] or slip-length [125];
• Type 3: link directly molecular and continuum solvers, each of them solving a portion of the same domain with momentum, mass and energy exchange at the interface while both are evolving together [114,115,118,119,120,125,131].

However, though these methods were applied to the molecular dynamics of reactive systems, they are still not used for water splitting applications. Possible extensions of the discussed multiscale approaches specifically tailored for application to water splitting will be discussed in the conclusions section.

3. Review–Utilization of MD to Water Splitting

Photoelectrochemical (PEC) water splitting is based on catalyzed charge transfer reactions at solid–liquid interfaces. Hence, a detailed understanding of charge transfer and the formation of intermediate and final reaction products, as well as the electronic charge carrier and ion distribution in this particular region of the PEC system, are of vital importance. However, the photoexcitation and transport of excited charge carriers take place within a semiconductor material and possibly involve solid–solid interfaces at semiconductor heterojunctions or semiconductor/co-catalyst interfaces. Moreover, the transport of ionic species within the liquid bulk phase is also highly relevant for optimizing the performance of PEC systems. Thus, to comprehensively describe and model PEC systems and reactions with numerical means, methods have to be employed and combined that are able to describe the electronic, optoelectronic and chemical properties of solid materials, the electrochemical reactions at solid–liquid interfaces, the structural properties of the electric double layer at solid–liquid interfaces, and ion transport within and the flow of liquids. Molecular dynamics simulations are able to tackle some, but not all, of these questions, and certain MD simulation methods and force fields are more or less suited to addressing specific problems. Hence, the applications of MD simulations in the field of water splitting are diverse and classical force field MD and ab initio MD, especially, have complementary applications.

In general, the description of electrochemical reactions requires the explicit incorporation of the electronic structure and high chemical accuracy such that quantum chemical ab initio based MD simulations are the natural choices for investigating chemical surface processes in water splitting processes. Moreover, to take into account effects from the liquid environment on the reaction energetics, free energy sampling methods, in combination with ab initio techniques such as DFT, are indispensable tools in this respect.

Classical Molecular Dynamics (MD), as described in detail in Section 2, belongs to the realm of computational chemistry and is one of the prominent techniques used for understanding the properties of the materials at the molecular scale level. It is based on the Newtonian equation of motion and at every simulation step, forces, velocities and positions of particles are updated. In particular, solid–water interphases are of vital importance for a wide range of systems, including photochemical conversion, corrosion and mineral activity, and so forth The catalytically reactive surface of these systems can be understood at the molecular scale. The particles interact through potentials—bonded and nonbonded—for a given number of simulation steps and mainly anticipate desorption and adsorption of reactants on a solid catalyst. Concerning water splitting on a solid surface, strong hydrogen bonding networks play a vital role in understanding and supporting fundamental phenomena such as the hydration of atoms/ions, electrochemical properties, diffusion, wetting and acid-base behavior. In particular, water splitting was investigated using molecular dynamics via photosynthesis, semiconductors and dye-synthesized PECs. First, we will discuss some studies related to the use of molecular dynamics simulations for photosynthesis, then molecular dynamics simulations for semiconducting materials for water splitting, and in the end molecular dynamics simulations of titanium-based materials for water splitting.

3.1. Molecular Dynamics Simulations for Photosynthesis

Classical MD simulations have proven to be useful for improving artificial photosynthesis processes [21] by simulating and understanding the flow of reactants over surfaces with active sites such as multi-scale interfacial phenomena at catalytic reactive surfaces. The authors of reference [1] have developed characterization methods for the flow of reactants induced by forcing along a specific dimension, over an active site treated as a metal ion. The dynamic and structural properties of the system are characterized using particle trajectory data. Additionally, repulsive LJ potentials were used in combination with Yukawa screened electrostatic potentials for modeling the absorption behavior at active sites. The simulations were run using LAMMPS [28]. Results were obtained from profiling the diffusion coefficient using coarse-grained schemes [4], the velocity in the vicinity of the active surface and radial distribution functions.

3.2. Molecular Dynamics Simulations Determining a Proper Level Alignment of Semiconducting Materials at Liquid/Solid Interfaces

One of the most important properties that determines whether a semiconducting material can be used for water splitting is the relative alignment of the valence and conduction band edges with respect to the water redox potentials. In a first approximation, the band alignment can be calculated by taking the vacuum level as a common energy reference. However, this approach neglects the formation of dipole layers at the solid–liquid interface and leads to inaccurate results. Recipes to correct vacuum calculations and to include electric dipole layers at the interface in a semi-empirical way were devised by some researchers around Zunger [132], but in general, such approaches do not take into account properties that are specific to the considered surface termination and material. One of the first attempts to include the effect of the solid–liquid interface on the band edge alignment was made in the works of Wu et al. [133,134], where classical MD simulations based on the TIP4P force field for water were combined with DFT calculations using the semilocal PBE functional. Here, the electrostatic potential profile across the interface was used to correct the band edge alignment. One of their main assumptions was that the error of the semilocal DFT calculations in calculating the conduction band edge of the investigated semiconductors TiO₂, WO₃, CdS, ZnSe, GaAs and GaP cancels with the error for the prediction of the water acceptor level, which was modeled as the lowest unoccupied molecular orbital of a solvated H₃O⁺ ion. While this error cancellation worked for the specific materials, it seems not to hold in general [132]. Similarly, ab initio MD simulations of TiO₂/water interfaces based again on the PBE functional, predicted the alignment of the conduction band edge with the normal hydrogen electrode (NHE) potential in reasonable agreement with experimental values, whereas the valence band edge position or ionization potential showed large deviations as presented by Cheng and Sprik [135]. This effect was traced back to the typical underestimation of the PBE functional of the TiO₂ bandgap. In their work, the free energy change of the reduction of TiO₂ at the NHE going along with the solvation of a proton was determined via free energy integration using a dummy proton with variable charge.

The importance of the exact nature of surface adsorbates for the position of conduction and valence band edges was investigated using DFT MD simulations by Zhou et al. [136]. These simulations predict that the relative positions of the band edges at {001} and {101} facets of anatase TiO₂ vary with the state of the adsorbed water molecules. Deprotonation of the water adsorbates favors electron accumulation on the {001} facets, whereas holes accumulate preferably at the {101} facets. Contrary to this, under acidic conditions, when no deprotonation occurs, the relative positions of the band edges close to the surfaces are reversed with respect to the anatase bulk. This means that the originally reductive {001} facets turn oxidative and likewise the oxidative {101} facets turn reductive.

In order to overcome the shortcomings of local and semilocal DFT to correctly describe the band gaps and absolute positions of the conduction and valence band edges while also taking into account the effects of the surface dipole layers, the many-body perturbation theory within the GW approximation in combination with ab initio MD have been employed. Early examples of this approach are found in the works of Pham et al. who investigated Si/water interfaces [137,138] containing aqueous interfaces of GaN and ZnO. In both cases, the overall simulation protocol consists of DFT–MD simulations of the solid–liquid interface to tackle the time-averaged interface-specific effects on the level alignment and separate GW calculations of both the bulk semiconductor and water systems. The connection between the different calculations is then made via the averaged electrostatic potential of the bulk regions in the interface model or, alternatively, via the averaged core potentials. Depending on the surface properties, typically the effect of the interfacial dipole was of the order of 0.5 eV but could even reach values of 1.5 eV [139].

3.3. Simulation of Water Splitting for Titanium Oxide Systems

In the field of semiconductor photocatalysis for water splitting [5,6], classical MD was employed in order to study titanium dioxide-based catalysts such as Pt/TiO₂ with two different configurations such as anatase (101) and rutile (110). These systems were investigated as water splitting photocatalysts [140,141,142] using LAMMPS. Universal force field (UFF) parameters were used to account for interactions between TiO₂ surfaces and H₂O. The authors extensively studied the behavior of water molecules over Pt-loaded titania surfaces and reported the following valuable characteristics: diffusion coefficient, radial distribution function, dipole angles, residence times and density profiles. The results revealed that H₂O dipoles do not have a preferred orientation and the Pt/rutile system shows a large residence time of over 30 ps.

Various studies have addressed the understanding of solid–water interphases [143,144,145,146,147,148]. It is important to analyze the interatomic potential that is used to describe the solid and the water–solid interface. In such a context, several materials and the corresponding potentials have been reviewed. In the case of TiO₂, various force fields have been proposed to elucidate the interaction between Ti and O atoms present in different morphologies and stoichiometries. Matsui and co-workers have developed a pair-wise potential for TiO₂ polymorphs; for instance, rutile, anatase, brookite, TiO₂. The authors reported a typical interatomic potential that takes into account the sum of Coulomb, dispersion and repulsion interactions, and the parameters can be found in reference [149]. On the other hand, Kim et al. proposed an alternative interatomic potential with an additional term to the initial Matsui potential composed of Coulomb, short-range repulsion, van der Waals and Morse interactions for rutile. The suggested potential reproduced the cell structures, bulk modulus and volume thermal expansivity; it can also be successfully applied to other polymorphs [150]. The reported potential was able to reproduce several properties that combine crystal structures, volume thermal expansion, volume compressibility and enthalpy relationships, and so forth. This typical potential is extensively used for MD simulations [151,152]. However, the anisotropic static relative permittivity of rutile is poorly calculated, and this is the main drawback of this potential.

In 2003, Bandura and Kubicki modified the potential for Ti-O with the Matsui and Kim potentials as starting points and quantum mechanical calculations as a reference [153]. In particular, non-electrostatic interactions are described in the form of the Buckingham potential (Equation (21)) and the interactions between Ti and O are given in

Table 1. Buckingham potential parameters of TiO₂ and atomic point charges.

i–j	Aij (kcal Mol−1)	ρij, Å	Cij (kcal Mol−1) Å6
Ti–O	315,480.8	0.194	290.3317
Ti–Ti	717,647.4	0.154	121.0676
O–O	271,716.3	0.234	696.8883
Atomic charges (e)
q(Ti)	2.196
q(O)	−1.098

.

$$E_{ij}=A_{ij}exp\left(-\frac{r_{ij}}{{\rho }_{ij}}\right)-C_{ij}/r_{ij}^6.$$

(21)

The majority of studies have focused on the development of potentials for perfectly crystalline TiO₂. However, in 2010, Schneider et al. proposed the formation of a native oxide layer on the Ti(0001) surface to elucidate the incorporation of O₂ on the surface as well as forming an amorphous interphase [154,155]. In these studies, a new potential for TiO₂ was introduced, called the Finnis–Sinclair many-body potential, which clearly illustrates the physical phenomena of the formed oxide layer on the surface and also explains the adsorption of water. It should be pointed out that the electrostatic interactions within the system were determined using the Electronegativity Equalization Method (EEM). Still, the electrostatic interactions between the oxide layer and the molecular species above the surface remain a problem in this EEM method.

Considering water models for classical MD simulation, many have been developed. Among them, simple three-site SPC/E, three-site transferable TIP3P and four-site transferable TIP4P are used [156,157,158,159]. These typical water models are able to accurately reproduce the properties of water at a molecular scale. Among the reported models for H₂O, the TIP4P potential shows great accuracy in representing the properties of water. However, the additional site in TIP4P demands significant computational costs, limiting extensive MD simulations with solid surfaces.

Bandura adopted the force field developed by Matsui for TiO₂ and additionally derived the interaction potential between relaxed TiO₂ surfaces and water using ab-initio calculations [153]. The fitted interaction potential was successfully applied to classical MD for rutile TiO₂ surfaces with SPC/E water. This force field was used to aid in the interpretation of quasielastic neutron scattering (QENS) spectra with water on rutile and anatase surfaces [141]. Predota and co-workers [160] studied the structural characteristics at the surfaces and the interfacial alignment of water for non-hydroxylated and hydroxylated rutile surfaces using ab-initio MD and classical MD and the studied Figure 4. The negatively charged surface was introduced in MD simulations to account for the effect of long-range interactions. The negatively charged surface can be counterbalanced by a surplus of cations in the solution; for instance Rb⁺ and Sr²⁺. The interesting results showed that the position of the water molecules at the interphase in both investigated systems (non-hydroxylated and hydroxylated) were found to be very similar in both terminal hydroxyl groups and terminal Ti atoms. The water molecules form two layers close to the surface and after that there is a transition to the bulk properties. Due to strong electric fields at hydroxylated surfaces, hydrogen-bonding and water orientation are more pronounced in those layers.

Skelton and Walsh made some minor modifications to the existing force field of Predota to describe the interaction between the rutile TiO₂ (110) surface and TIP3P water [161]. The optimized force field for non-hydroxylated surfaces was consistent with the final optimized structure obtained from periodic density functional calculations. Six different surfaces were considered in order to study the interaction with TIP3P water and the obtained results were in good agreement with the SPC/E model. Therefore, this typical water model can also be used for modeling with TiO₂. Furthermore, Nada et al. carried out MD simulations for 001 and 110 plane configurations of TiO₂ with TIP3P [162]. The results indicated that the water layers at the interphase play a dominant role in which it strongly affects the dynamics of ions. Moreover, water molecules are strongly bound to the TiO₂ (110) surface rather than to the (100) surface, which hindered the growth of the other ions. Similarly, the authors developed a new force field by modifying initial Bandura and Kubicki parameters for water with TiO₂ [163]. The suggested new set of Buckingham parameters was able to yield binding energies and conformations in good agreement with DFT results and experiments. The parameters also exclusively describe the adsorption of the hydroxyl groups that formed through dissociated water on the different surfaces of TiO₂. A recent MD review by YazdanYar et al. extensively investigates the interaction of ions and biological organic molecules on the rutile surface [148].

The van Duin group has developed a ReaxFF reactive force field, which is described in Section 2, to describe reactions in the Ti-O-H system with the help of quantum mechanical reference calculations. The developed ReaxFF reproduces the structures and energetics of clusters as accurately as the DFT calculations do, and also clearly describes the relative energies of different surfaces such as rutile, brookite and anatase. Water-binding energies, surface energies and H₂O dissociation energies also reasonably match DFT results [48]. In another study, they reported that a hydration layer close to the surface of anatase enhances water dissociation. Moreover, ReaxFF describes water binding and dissociation on various titania surfaces for water coverage within a large range of temperatures. Therefore, ReaxFF can be used for large-scale simulations of TiO₂/water interfaces [164]. The reported ReaxFF was successfully applied to TiO₂ surfaces (anatase and rutile) in order to understand the interaction of H₂O at the interphase [140]. The results revealed that interfacial water was heavily confined and the corresponding self-diffusion coefficient increased by two orders of magnitude. The steric hindrance has evolved on the adsorbed water molecules, and the hydrogen-bond life-time on the interphase was relatively shorter than that of bulk water. The calculated IR frequencies of the stretching vibrations of anatase and rutile have exhibited a stronger variation.

One more approach towards photocatalytic water splitting is to use dye-sensitized PEC [165,166,167,168,169]. In [170], the authors use constrained ab initio MD simulations to study the catalytic mechanism in water oxidation. In more detail, the authors studied a supramolecular complex integrating a mononuclear Ru-based water oxidation catalyst (WOC) with a fully organic naphthalene-diimide (NDI) dye for fast photoinduced electron injection into the conduction band of the titanium-dioxide semiconductor anode. The MD simulations have shown that the chosen dye has a good potential to be integrated into a dye-sensitized PEC device.

Viswanathan et al. [171] presented a multiscale model to simulate the linear sweep voltammetry of the electrochemical oxidation of water on Pt (111) and Pt3Ni (111). DFT is used to parameterize the reaction kinetics, and kMC is used to capture the kinetic steps of the electrochemical reactions. In this paper, the resulting voltammogram from the computations shows good agreement with the experimental result. The theoretical predictions show that OH adsorbs between 0.65 and 0.85 V_RHE and that, at 0.9 V_RHE, OH starts to become oxidized to O [113].

4. Conclusions

Efficient energy storage is a key requirement for the transition to cleaner and more renewable energy production technologies such as solar power. Energy storage in chemical fuels such as hydrogen can be obtained by water splitting, with photoelectrochemical water splitting being a prime candidate. However, up to this date, no photoelectrochemical system has managed to be both economically viable and scalable. This work reviews how molecular dynamics simulations can be employed in the study of the water splitting process in order to test and design materials empowering efficient photoelectrochemical systems. A general review of molecular dynamics is presented before the work focuses on specific use cases for molecular dynamics in water splitting. The first series of reviewed simulation works investigated the use of molecular dynamics simulations of processes at the solid–liquid interface that determine, amongst other things, the alignment of semiconductor electronic band edges with respect to the water redox potentials. Another reviewed series of molecular dynamics studies investigated titanium oxide systems for water splitting applications. While molecular dynamics simulations, especially for titanium dioxide, are well established, simulations of other water splitting materials are vastly lacking. An overall conclusion from this review is that there are few molecular dynamics studies for water splitting compared to, for example, the vast number of DFT simulations concerning the electronic structure of water splitting materials or of specific surface reactions.

Various directions for combining molecular dynamics with other methods that combine the atomistic scale with electronic structure, such as continuum or Monte Carlo methods, have been actively developed but are still not widely used in the water splitting community. From our review related to multiscaling molecular dynamics approaches, we identified the following gaps with respect to their application to water splitting modeling:

• Molecular Dynamics and Quantum mechanics: On the multiscaling side, the combination of these two methods is more mature as compared to Molecular Dynamics and Continuum modeling. One possible gap can be the use of artificial intelligence for the derivation of parameters for Molecular Dynamics and/or Quantum Mechanics in an automatic way;
• Molecular Dynamics and Monte Carlo: Again, in this area there are more mature methods available as compared to the coupling of Molecular Dynamics and Continuum Modeling. A gap can still be identified in the sense of developing tooling and methods to combine the two methods tailored for the application to water splitting;
• MD and state-space modelling (SSM): As in the case of MD and Monte Carlo, an identified gap is the development of more and easy-to-use tools and methods to combine the two methods specifically for the application to photoelectrochemical water splitting;
• Molecular Dynamics and Continuum modeling: The methods here are in a prototyping phase. As both continuum and molecular dynamics modeling techniques are quite diverse, there is sufficient room for new methods and ways to combine the two. A lot of research is still needed in this area.