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.

Absorption and reactive properties of fluids in porous media are key to the design and improvement of numerous energy related applications. Molecular simulations of these systems require accurate force fields that capture the involved chemical reactions and have the ability to describe the vapor–liquid equilibrium (VLE). Two new reactive force fields (ReaxFF) for CO₂ and H₂O are developed, which are capable of not only modeling bond breaking and formation in reactive environments but also predicting their VLEs at saturation conditions. These new force fields include extra terms (ReaxFF-lg) to improve the long-range interactions between the molecules. For validation, we have developed a new Gibbs ensemble Monte Carlo (GEMC–ReaxFF) approach to predict the VLE. Computed VLE data show good agreement with National Institute of Standards and Technology reference data as well as existing nonreactive force fields. This validation proves the applicability of the GEMC–ReaxFF method to test new reactive force fields, and simultaneously it proves the applicability to extend newly developed ReaxFF force fields to other more complex reactive systems.

Thermochemical heat-storage applications, based on the reversible endo-/exothermic hydration reaction of salts, are intensively investigated to search for compact heat-storage devices. To achieve a truly valuable storage system, progressively complex salts are investigated. For these salts, the equilibrium temperature and pressure conditions are not always easy to predict. However, these conditions are crucial for the design of thermochemical heat-storage systems. A biased grand-canonical Monte Carlo (GCMC) tool is developed, enabling the study of equilibrium conditions at the molecular level. The GCMC algorithm is combined with reactive force field molecular dynamics (ReaxFF), which allows bond formation within the simulation. The Weeks–Chandler–Andersen (WCA) potential is used to scan multiple trial positions for the GCMC algorithm at a small cost. The most promising trial positions can be selected for recomputation with the more expensive ReaxFF. The developed WCA–ReaxFF–GCMC tool was used to study the hydration of MgCl₂·nH₂O. The simulation results show a good agreement with experimental and thermodynamic equilibriums for multiple hydration levels. The hydration shows that water, present at the surface of crystalline salt, deforms the surface layers and promotes further hydration of these deformed layers. Additionally, the WCA–ReaxFF–GCMC algorithm can be used to study other, non-TCM-related, reactive sorption processes.

The interface of two solids in contact introduces a thermal boundary resistance (TBR), which is challenging to measure from experiments. Besides, if the interface is reactive, it can form an intermediate recrystallized or amorphous region, and extra influencing phenomena are introduced. Reactive force field Molecular Dynamics (ReaxFF MD) is used to study these interfacial phenomena at the (non-)reactive interface. The non-reactive interfaces are compared using a phenomenological theory (PT), predicting the temperature discontinuity at the interface. By connecting ReaxFF MD and PT we confirm a continuous temperature profile for the homogeneous non-reactive interface and a temperature jump in case of the heterogeneous non-reactive interface. ReaxFF MD is further used to understand the effect of chemical activity of two solids in contact. The selected Si/SiO₂ materials showed that the TBR of the reacted interface is two times larger than the non-reactive, going from $1.65\times {10}^{-9}$ to $3.38\times {10}^{-9}$ m²K/W. This is linked to the formation of an intermediate amorphous layer induced by heating, which remains stable when the system is cooled again. This provides the possibility to design multi-layered structures with a desired TBR.