Insights into the Structure and Dynamics of Water at Co$_3$O$_4$(001) Using a High-Dimensional Neural Network Potential

TL;DR

This study employs a high-dimensional neural network potential trained on DFT data to simulate and analyze the atomistic structure, dynamics, and reactivity of Co$_3$O$_4$(001)-water interfaces over extended time scales, revealing key surface behaviors.

Contribution

The paper introduces a neural network potential that extends ab initio molecular dynamics simulations to larger scales, enabling detailed investigation of Co$_3$O$_4$-water interfaces.

Findings

Distinct hydration structures for A and B terminations.

B-termination stabilizes a well-organized hydration layer.

Proton transfer pathways and surface hydroxylation dynamics.

Abstract

Co$_3$O$_4$ is an important catalyst for the oxidation of organic molecules in the liquid phase. Still, understanding the atomistic details of Co$_3$O$_4$-water interfaces under operando conditions remains extremely challenging. While ab initio molecular dynamics have become an essential tool for investigating these dynamic interfaces in silico, they are limited to only a few picoseconds and a few hundred atoms. In this work, we overcome these limitations by training a high-dimensional neural network potential (HDNNP) on density functional theory data, which allows us to significantly extend the accessible time and length scales. Employing this HDNNP, we perform simulations to unravel the structure, dynamics, and reactivity of Co$_3$O$_4$(001)-water interfaces in detail. Our simulations reveal distinct characteristics of the two possible A and B terminations. The B-terminated surface stabilizes a compact, quasi-epitaxial hydration layer with strong templating effects, enhanced hydroxylation, and a well-organized hydrogen-bond network. In contrast, the A-termination forms a more diffuse contact layer with weaker templating, lower hydroxylation, and less ordered interfacial water. Extended simulations further uncover proton transfer pathways, including intermittent protonation of surface hydroxyls, migration of water molecules into the epitaxial layer, and rare hydronium-like configurations.

Figures

Equations

$$E=\sum _{\alpha =1}^{N_{\text{elements}}}\sum _{n=1}^{N_{\text{atoms}}^{\alpha }}{E}_n^{\alpha }.$$

Full text

I Introduction

Cobalt oxide (Co3O4) is a transition metal oxide with mixed valence states that has attracted significant interest in recent years, owing to its unique chemical, physical, and electronic characteristics. In particular, Co3O4 exhibits rich redox behavior due to the coexistence of Co2+ and Co3+ ions, remarkable catalytic activity, tunable electronic structure, and intriguing magnetic propertiesAnke et al. (2019); Doheim and El-Shobaky (2002); Yang et al. (2021); Omranpoor et al. (2023); Anke et al. (2020); Douma et al. (2023); Omranpoor et al. (2022); Kenmoe et al. (2022); Omranpoor and Kenmoe (2023). These properties make Co3O4 highly promising for diverse technological applications, including its role as a catalyst in alcohol oxidationFalk et al. (2021), water oxidationFeizi et al. (2019); Jiao and Frei (2009), methane combustionHu et al. (2008), and CO oxidationXie et al. (2009), as well as for use in lithium-ion batteries and gas sensorsLi et al. (2005). Among the different structural forms, the Co3O4 spinel has been particularly well studied for its effectiveness in promoting oxidation reactions, notably in the selective oxidation of hydrocarbonsWaidhas et al. (2020); Hill and Hartwig (2017); Finocchio et al. (1997).

As a catalyst, Co3O4 has been widely used for the oxidation of alcohols. Recently, there has been a growing interest in carrying out catalytic reactions such as oxidation of alcohols in the liquid phaseFalk et al. (2021), typically in aqueous solutions, which allows for milder conditions enhancing selectivity. However, for Co3O4 as a catalyst, the role of the aqueous environment remains unclear, as both beneficial and detrimental effects have been reported in both computationalWang et al. (2012) and experimentalFalk et al. (2021) studies.

A combination of density functional theory (DFT) calculations and high-resolution scanning transmission electron microscopy (HRSTEM) was employed by Zasada et al. to study the morphology of cobalt spinel nanostructuresZasada et al. (2011). Their surface energy calculations indicated the stability sequence (001) $>$ (111) $>$ (110) for the low-index facets in vacuum, which aligns with the findings of Montoya et al.Montoya and Haynes (2011) from an independent DFT investigation. Utilizing these surface energies, they performed a Wulff construction to predict the equilibrium shape of Co3O4 nanoparticles. The resulting model revealed rhombicuboctahedral grains predominantly exposing the (001) and (111) facets, accounting for 48% and 41% of the surface area, respectively, along with a smaller contribution from the (101) facet of about 11%. These theoretical predictions are consistent with HRSTEM images of the synthesized spinel nanocrystalsZasada et al. (2011).

Although extensive experimental studiesSchwarz et al. (2018); Budiyanto et al. (2023); Natarajan et al. (2021); Tran-Phu et al. (2022); Qiu et al. (2024); Wiegmann et al. (2022); Reikowski et al. (2019); Haunold et al. (2025); Zhang et al. (2018); Varhade et al. (2023) have been performed on Co3O4–water interfaces, they often lack the spatial and temporal resolution required to investigate these complex interfaces in operando. Thus, computational studies are crucial to achieve a comprehensive atomistic understanding. To date, investigations using static DFT calculationsChen and Selloni (2012); Kaptagay et al. (2015); Yan and Sautet (2019); Yan et al. (2019); Haunold et al. (2025); Huo et al. (2024); Peng et al. (2021) have been the dominant theoretical approach, albeit they do not explicitly account for the impact of solvation and the dynamic nature of the interface. Ab initio molecular dynamics (AIMD) simulations overcome this limitation by intrinsically incorporating dynamics and the role of finite temperatures. So far, there have been three main AIMD studies that focused purely on Co3O4–water interfaces: Kox et al. on (001)Kox et al. (2020), Gaigeot et al. on (110)Creazzo et al. (2019), and Kox et al. on (111)Kox and Kenmoe (2024). As these are computationally demanding AIMD studies, all three works were limited to simulation times of about 20 picoseconds and systems consisting of a few hundred atoms. One of the principal aims of the present work is to explicitly examine the effect of time and length scales in atomistic simulations and to extend the study of Kox et al.Kox et al. (2020) on the (001) surface to significantly larger scales. In their work, increasing water coverage enhanced Co2+ solvation and stabilized the A-terminated surface, while for the B-termination it primarily added diffuse outer water layers without affecting the rigid first layer or the dissociation degree. However, their simulations were restricted to up to 2 monolayers (32 water molecules), thus lacking the explicit and implicit influence of bulk water. In the present study, we extend these simulations to more than 20 monolayers in a large supercell containing about 2000 water molecules, thereby allowing to describe realistically the role of solvation.

There are two possible terminations of the Co3O4(001) surface: the A-termination, which is predominantly terminated by Co2+ ions, and the B-termination, which is predominantly terminated by Co3+ ions (for more details, see Ref. 38). Experimentally, determining which termination is more favorable for catalytic reactions is of utmost importance, as this holds the key to explaining why certain cation dopants, such as Fe3+, reduce the catalytic activity of Co3O4 Falk et al. (2021). One possible explanation is that the incorporation of Fe3+ preferentially occurs at Co3+ sites, which are only exposed in the B-termination, potentially altering its catalytic behavior. Therefore, a detailed characterization of these two surface terminations is a central focus of the present study.

The spatial and temporal limitations of AIMD simulations have largely been overcome by the development of machine learning potentials (MLPs) Behler (2016); Deringer et al. (2019); Unke et al. (2021); Friederich et al. (2021); Kocer et al. (2022), which have become a mainstay in atomistic simulations of condensed systems and are now extensively used to study solid–liquid interfacesHou et al. (2024); Schran et al. (2021); Omranpour et al. (2024, 2025). However, developing accurate MLPs for oxides of comparable complexity, particularly those containing transition metals in multiple oxidation states, remains a significant challenge for theoretical studies. To date, the only other examples include recent investigations on LiMn2O4 Eckhoff et al. (2020, 2020); Eckhoff and Behler (2021), hematite (Fe2O3) Schienbein and Blumberger (2022), and magnetite (Fe3O4) Romano et al. (2025, 2024), all of which employed high-dimensional neural network potentials (HDNNPs)Behler and Parrinello (2007); Behler (2021). For magnetic oxides a significant complication arises from their complex spin-dependent properties, which pose challenges even for DFT, often leading to convergence issues or metastability of high-energy configurations Meredig et al. (2010). This, in turn, complicates the construction of accurate and robust MLPs, as their quality fundamentally depends on consistent, high-quality reference data. These challenges are well documented in Ref. 54, where both Gaussian Approximation Potentials (GAP)Bartók et al. (2010) and HDNNPs initially struggled to produce stable potentials for the magnetite (Fe3O4)–water interface, primarily due to inconsistencies in the DFT+U reference dataset.

In the present work, we overcome the spatial and temporal limitations of previous ab initio molecular dynamics studies by developing and employing a HDNNP for Co3O4(001)–water interfaces. This enables us to perform nanosecond-scale simulations of realistic interfacial systems comprising thousands of atoms. Using this approach, we present a detailed atomistic investigation of the Co3O4(001)–water interface, with particular focus on comparing the A-terminated and B-terminated surfaces. We analyze their structural ordering, interfacial water dynamics, surface hydroxylation, and dynamic surface reconstruction, thereby revealing how surface termination governs the properties of interfacial water in case of this catalytically relevant oxide.

II Methods

Machine Learning Potentials are supervised regression models that establish a mapping from atomic structures to their associated potential energies, with labeled data consisting of structures and their corresponding energies and often forces. Unlike ab initio molecular dynamics, where energies and forces are explicitly computed on-the-fly at each timestep using electronic structure methods, MLPs such as HDNNPs provide a continuous function that represents the potential energy surface (PES). This function efficiently predicts the total energy as well as its derivatives like the forces for a given atomic configuration.

Developing an MLP involves two essential tasks: representation and regression. The representation problem is addressed by employing descriptors that encode the local atomic environments, transforming Cartesian coordinates into features that are invariant under translation, rotation, and permutation of identical atoms. In the case of HDNNPs, this is typically achieved through atom-centered symmetry functions (ACSF) Behler (2011), which act as structural fingerprints to ensure that equivalent configurations yield identical energy predictions. The regression task is then performed by training a multilayer feed-forward neural network to learn the relationship between these descriptors and the target energies (and optionally forces), thus constructing an accurate model of the PES.

In HDNNPs, the local chemical environment of an atom, defined by all neighboring atoms within a cutoff sphere of radius $R_{\mathrm{c}}$, is encoded by a vector of ACSF values. In second-generation HDNNPs Behler and Parrinello (2007), the cutoff $R_{\mathrm{c}}$ must be sufficiently large to capture all relevant interactions, and typically cutoffs between 5 and 10 Å are employed. In this work, two types of ACSFs are used: radial and angular symmetry functions, as defined in Ref. 60. Depending on the complexity of the system, between 30 and 150 ACSFs per atom are typically employed.

Beyond specifying the chemical elements of the atoms, no additional information such as atom types, fixed oxidation states, or predefined bonds is required. This makes HDNNPs inherently reactive, capable of accurately describing bond formation and breaking as well as changes in oxidation states, all governed by the underlying electronic structure. Since the dimensionality of the ACSF vectors is determined by the chosen set of symmetry functions and is independent of the particular local environment, these vectors can serve directly as input vectors to atomic neural networks that require a fixed dimensionality.

A separate atomic neural network is constructed for each chemical element $\alpha$, which processes the local structural environment of each atom $n$ to produce its atomic energy contribution $E_{n}^{\alpha}$. The total energy $E$ of the system, comprising $N_{\text{elements}}$ elements and $N_{\text{atoms}}^{\alpha}$ atoms of element $\alpha$, is then given by

[TABLE]

These atomic contributions are predicted by multilayer feed-forward neural networks, which use a linear activation function in the output layer and, in most cases, hyperbolic tangent functions in the hidden layers. For each atom, the ACSF vector is computed and passed through the appropriate neural network, with the resulting atomic energies summed to yield the total potential energy of the system. Forces, which can also be used for training, can then be obtained as analytical derivatives of the energy with respect to the atomic positions.

The weights of all atomic neural networks are optimized simultaneously using iterative gradient-based procedures. In this study, an adaptive, global, extended Kalman filter Kalman (1960); Blank and Brown (1994) is employed to minimize the errors of both total energies and atomic force components for a DFT training set of reference structures. For further details on the construction, training, and applications of HDNNPs, the reader is referred to several comprehensive reviews Behler (2021, 2017, 2014, 2015); Tokita and Behler (2023); Omranpour et al. (2024, 2025).

III Computational Details

III.1 Density Functional Theory Calculations

All electronic structure calculations were performed using the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996, 1996) version 6.3.2, employing spin-polarized DFT. Exchange–correlation effects and van der Waals interactions were treated using the optPBE-vdW functional Perdew et al. (1996); Klimeš et al. (2009, 2011). On-site Coulomb interactions were further accounted for via the DFT+U approach to properly treat the strongly correlated Co $3d$ electrons, with an effective Hubbard U parameter of 2.43 eV determined following Dudarev et al.Dudarev et al. (1998).

Ionic cores were represented using Projector Augmented Wave (PAW) potentials Blöchl (1994), as formulated by Kresse and Joubert Kresse and Joubert (1999); Kresse and Furthmüller (1996). Plane-wave expansions were performed with a kinetic energy cutoff of 500 eV. The Brillouin zone was sampled using a Monkhorst-Pack scheme with a $5\times 5\times 5$ k-point mesh for bulk Co3O4, a $3\times 3\times 3$ k-point mesh for bulk water, and a $3\times 3\times 1$ k-point mesh for the Co3O4–water interfaces. The typical bulk Co3O4 systems consisted of 56 atoms, with a simulation box of about $8\,\text{\AA }\times 8\,\text{\AA }\times 8\,\text{\AA }$ (see Ref. Omranpour et al. (2025)). The typical systems for bulk liquid water consisted of 192 atoms, with a simulation box of about $12.42\,\text{\AA }\times 12.42\,\text{\AA }\times 12.42\,\text{\AA }$, in addition to some ice structures and liquid–vacuum interfaces (see Ref. 75). The interface systems consisted of 176 atoms, and the typical simulation box dimensions were about $8\,\text{\AA }\times 8\,\text{\AA }\times 20\,\text{\AA }$. Gaussian smearing of 0.1 eV was applied to determine partial occupancies, and non-spherical contributions within the PAW spheres were included. Electronic self-consistency was considered achieved when total energy differences fell below $10^{-6}$ eV.

III.2 Construction of the Reference Data Set

When constructing the reference data for systems such as Co3O4–water interfaces, it is necessary to sample three regions of configuration space corresponding to three types of atomic environments: bulk Co3O4, bulk water, and the Co3O4–water interfaces Omranpour et al. (2025). For the bulk Co3O4 reference data, the set of 10,700 structures from Ref.76 was adopted without modification. This dataset spans temperatures from 0 K up to 700 K, thus exceeding the temperature range typically used in catalysis at solid-liquid interfaces. For bulk water, 814 structures from Ref. 75 were employed, and energies and forces were recomputed using the electronic structure settings of this work. This particular water dataset was chosen for its efficiency and diversity, as it was curated from a large pool of aqueous configurations to accurately reproduce the key properties of water.

To build an initial dataset for Co3O4–water interfaces for training a preliminary set of HDNNPs, AIMD simulations were conducted with VASP on a single unit cell of Co3O4 exposed to 40 water molecules, with the (001) surfaces consisting of both A- and B-terminations in contact with water. The electronic structure settings mirrored those described earlier, except that a $\Gamma$-centered $1\times 1\times 1$ k-point grid was used to enhance computational efficiency. Simulations were performed in the $NVT$ and $NPT$ ensembles at temperatures ranging from 300 K to 600 K to generate diverse interfacial structures, using a Nosé–Hoover thermostat to regulate temperature Evans (1985). Each simulation ran for approximately 15 ps with a time step of 0.5 fs. The initial 4000 configurations (corresponding to about 2 ps) were discarded to allow for equilibration. Subsequently, configurations were extracted every 50 steps, yielding an initial dataset of roughly 3000 structures. For these configurations, single-point electronic structure calculations were performed with a dense $3\times 3\times 1$ k-point mesh to obtain precise energies and forces.

Following the creation of this dataset and training of initial HDNNPs, an active learning strategy, following Eckhoff et al. Eckhoff and Behler (2021, 2019), was employed to identify additional structures associated with high prediction uncertainty across an ensemble of similarly trained HDNNPs. The HDNNPs differed by random seeds, resulting in varied train/test partitions and different weight initializations. Electronic structure calculations were then performed on these high-uncertainty structures, which were added to the dataset to further improve the HDNNPs. This active learning cycle was repeated multiple times until the variance across all sampled structures was comparable to the RMSEs of the models Tokita and Behler (2023).

III.3 Construction of the High-Dimensional Neural Network Potential

The HDNNPs were constructed using the RuNNer code Behler (2015, 2017) (version dated August 16, 2023). The atomic neural networks for each element employed two hidden layers with 25 and 20 neurons, respectively. The input nodes corresponded to the element-specific ACSFs, and a single output node provided the atomic energy. A cutoff radius of $R_{c}=12\,a_{0}$ (6.35 Å) was used for the ACSFs to ensure an accurate description of atomic interactions in the Co3O4–water interfaces. The Kalman filter parameters were set to $\lambda=0.98000$ and $\nu=0.99870$. Further details on the employed symmetry functions and their parameters, the RuNNer settings for constructing the HDNNP, as well as energy and force correlation plots, are provided in the supporting information.

After completing the iterative generation of the reference dataset via active learning, a total of 3540 Co3O4(001)–water interface structures (each containing 176 atoms) were obtained. Combined with the previously generated 10,766 bulk Co3O4 structures Omranpour and Behler (2024) and 812 bulk water structures Schran et al. (2020), this resulted in a dataset of 15,118 structures, providing 15,118 potential energies and their associated atomic force components. Of these, 90% were used for training the HDNNP and 10% for testing. In each epoch 100% of all energies and 5% of all forces were used in the training process. The root mean squared error (RMSE) for the energy is 1.201 meV/atom for the training set and 1.288 meV/atom for the test set, while the RMSE for the force components is 0.1277 eV/Å and 0.1284 eV/Å, respectively, which are typical for state-of-the-art MLPs.

III.4 Molecular Dynamics Simulations

Molecular dynamics simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) Plimpton (1995) (version from 2nd August 2023), incorporating the n2p2 library (version 2.2.0 from 23rd May 2022) Singraber et al. (2019) for employing HDNNPs. Simulations used for active learning Eckhoff and Behler (2021, 2019) were conducted in the $NVT$ and $NPT$ ensembles in a temperature range from 300 to 450 K, with a time step of $\Delta t=0.5\,\text{fs}$. The velocity Verlet scheme Swope et al. (1982) was used as the integration algorithm. Temperature and pressure were controlled using the Nosé–Hoover thermostat and barostat Nosé (1984); Hoover (1985); Evans (1985), with all $NPT$ simulations performed at $p=1$ bar. For the final simulations, systems were first equilibrated in the $NPT$ ensemble for a 1 ns run at 300 K and 400 K, and the resulting configurations were then used to perform production runs in the $NVT$ ensemble at 300 K and 400 K, respectively, for an additional 1 ns. The same simulation setup was previously applied to the ZnO–water interface using HDNNPs Quaranta et al. (2017). Figure 1 shows the $4\times 4\times 2$ Co3O4 supercell representing the slab in the MD simulations, with the system comprising a total of 7,936 atoms. The top surface of the slab corresponds to the A-termination, while the bottom surface corresponds to the B-termination.

IV Results and Discussion

The structure of this section is as follows: first, a few of the most important reaction mechanisms near the interfaces are discussed in detail (Section IV.1). This is followed by an analysis of the contact layer structure, where snapshots of different contact layers and radial distribution functions are provided to characterize the water structure near the surfaces (Section IV.2). Subsequently, we discuss the density profiles of O and H atoms to describe the overall water distribution throughout the system (Section IV.3). Thus, the section progresses from a detailed, zoomed-in view of the reaction mechanisms to the contact layers and finally to the entire system. In the end, we discuss the implications of our findings for catalysis (Section IV.4).

IV.1 Proton Transfer Mechanism

In liquid water and at oxide–water interfaces, protons can migrate via proton transfer events (i) in the bulk liquid through the hydrogen-bond network or (ii) near the surface through reactions involving adsorbed species (e.g., hydroxyls) and/or lattice oxygen atoms. To visualize such events in our trajectories, we monitor, for selected oxygen atoms (the “tagged” O), the time series of its distances to the nearest four hydrogen atoms (H1–H4) and, when relevant, to the nearest surface Co atom. In these plots, two persistent O–H distances near $\approx 1.0~\text{\AA }$ indicate a neutral water molecule; one near $1.0~\text{\AA }$ indicates a hydroxyl; and three near $1.0~\text{\AA }$ sustained over multiple frames is the hallmark of a hydronium-like configuration. H3/H4 values in the $1.5\!-\!2.4~\text{\AA }$ range report non-bonded, hydrogen-bond donor hydrogens. An O–Co distance near $2.0~\text{\AA }$ signals an inner-sphere Co–O bond (chemisorbed OH or H2O), whereas a distance of $3.2\!-\!3.6~\text{\AA }$ is typical for outer-sphere water in the epitaxial/solvation layer.

Outer hydration layer case (A-termination).

Figure 2a shows a tagged water molecule located relatively far from the surface. Its initial distance from the A-terminated surface is about $10~\text{\AA }$, and the closest it approaches is roughly $4~\text{\AA }$, as indicated by the black trace in the figure. During the brief hydronium-like episode discussed below, its distance from the surface is around $5~\text{\AA }$, implying that even at this range, the molecule remains under the influence of the surface. The red and blue traces (H1, H2) stay close to $1.0~\text{\AA }$ throughout, confirming two covalent O–H bonds and no persistent dissociation. The green/orange traces (H3, H4) fluctuate around $1.6\!-\!2.2~\text{\AA }$, consistent with typical hydrogen-bond donors in the first and second hydration shells. Near the end of the trajectory, H3 transiently approaches $\sim 1.0~\text{\AA }$ while H4 simultaneously recedes toward $\sim 2.5\!-\!2.7~\text{\AA }$. This concerted behavior corresponds to a short-lived hydronium-like configuration in which the tagged O momentarily binds three protons before rapidly relaxing back to water. In this example, the residence time is on the order of $\sim 50$ ps, underscoring that protonation events in classical MD can be rare and require long sampling to capture convincingly.

Chemisorbed hydroxyl on Co (B-termination).

In Fig. 2b the black O–Co curve indicates a distance of about $2.0~\text{\AA }$, identifying an inner-sphere species bound to a surface Co site. Only one O–H trace (red) stays at $\approx 1.0~\text{\AA }$, as expected for an OH group. The second-nearest H (blue) resides around $1.6~\text{\AA }$ but sporadically dips to $\approx 1.0~\text{\AA }$, indicating transient protonation–deprotonation of the surface OH to form/re-form molecular water at the same site. These excursions are more frequent in the early part of the run, suggesting that the interfacial acid–base equilibrium is still relaxing before settling into a steady population of bound hydroxyls and outer-sphere water molecules.

Water migrating from the outer hydration layer to the epitaxial layer (B-termination).

Figure 2c shows a tagged water molecule that initially resides in the outer hydration region, with an O–Co distance of $\sim 4.5\!-\!5.0~\text{\AA }$ during the first $\sim 100$ ps of the simulation. After this period, the molecule gradually moves closer to the surface and relocates into the epitaxial (outer-sphere) water layer, where its O–Co distance stabilizes around $\sim 3.3\!-\!3.6~\text{\AA }$. This relocation places the molecule in close interfacial proximity but still distinct from chemisorbed species, such as those in Figure 2b, which bind at $\sim 2.0~\text{\AA }$. The H1/H2 traces remain near $\approx 1.0~\text{\AA }$, confirming the molecule retains its molecular water character. Intermittent approaches of H2 toward $\approx 1.6~\text{\AA }$ indicate transiently strengthened hydrogen bonds with nearby OH/H2O groups; these become somewhat more frequent after the relocation, suggesting enhanced interfacial exchange once the molecule has joined the epitaxial layer (see Figure 4(a)).

Water in the epitaxial layer with a late protonation event (B-termination).

The tagged O–Co distance in Fig. 2d is again $\approx 3.3\!-\!3.6~\text{\AA }$, diagnosing an outer-sphere water molecule at the interface. The most striking feature is the progressive approach of H3 (green) toward $\approx 1.0~\text{\AA }$ in the final tenth of the trajectory, while H4 (orange) simultaneously recedes, mirroring the bulk case in panel (a). This pattern is consistent with the formation of a hydronium-like configuration in the epitaxial layer near the end of the simulation and is corroborated by our coordination-count check, which flags frames where the tagged oxygen binds three hydrogens within a tight covalent cutoff.

Across (b–d), all examples are taken from the B-terminated interface for consistency. Where hydronium-like behavior is inferred, we verify that three O–H distances are simultaneously near $1.0~\text{\AA }$ over multiple frames; isolated single-frame spikes can also arise from fluctuating hydrogen bonds and should not be over-interpreted. Likewise, an O–Co distance near $2.0~\text{\AA }$ is a robust fingerprint of a chemisorbed OH/H2O, whereas $\approx 3.3\!-\!3.6~\text{\AA }$ indicates outer-sphere binding in the epitaxial layer and is consistent with the second peak in the interfacial Co–Ow RDF (Fig. 4a). Finally, the observation of late-time protonation events highlights the importance of long trajectories: such rare events are easily missed in short AIMD runs but emerge in extended classical MD sampling.

Overall, it can be concluded that (i) Bulk and interfacial protonation show the same mechanistic signature (three near-$1.0~\text{\AA }$ O–H distances accompanied by the retreat of the fourth H), yet occur in different structural environments (no Co partner vs. outer-/inner-sphere binding). (ii) Surface OH groups undergo intermittent protonation/deprotonation while staying attached to Co, mediating charge/proton mobility at the interface. (iii) Waters entering the epitaxial layer can participate in enhanced hydrogen-bond rearrangements and, occasionally, hydronium-like episodes late in the trajectory, especially on the B-termination.

IV.2 Structure of the Contact Layers: A-termination vs. B-termination

Due to the complexity and diverse chemistry at the contact layers, this section first presents several analyses, including atomic snapshots of the contact layers, trajectory lines of the surfaces and the first interfacial layers, and radial distribution functions (RDFs). The discussion then follows, based on the combined insights from these analyses.

IV.2.1 Snapshots, Trajectory Lines, and Radial Distribution Functions

Figure 3 illustrates the contact layers of the equilibrated Co3O4(001)–water interfaces at 300 K for the A-termination on the left and the B-termination on the right. Row 1 shows side views of only the first interfacial layer with respect to the outermost surface plane. (Color code: tetrahedral Co2+ green, octahedral Co3+ purple, lattice O red, interfacial oxygen species blue, H white.) The main difference between the A- and B-terminated surfaces is that, in the A-termination, the tetrahedral Co2+ (green) are located at the surface and directly exposed to the interfacial layer, whereas in the B-termination, the octahedral Co3+ (purple) occupy the surface and directly interact with the interfacial water layer.

Row 2 shows similar side views but overlays the trajectory lines of the surface Co atoms (Co2+ green for the A-termination and Co3+ purple for the B-termination) and the oxygen species of the first interfacial water layers over the last 50 ps of the simulation. Due to the highly dynamic and reactive nature of the solvent layer, it is difficult to clearly distinguish between H2O molecules, OH- groups, and transient H3O+ species. The boundaries between these species are blurred near the interface; for example, a water molecule can exchange a proton with a surface oxygen, a neighboring OH-, or another H2O, and may even transiently share an excess proton. Therefore, in these snapshots, all oxygens of the interfacial layer are consistently shown in blue (in contrast to the red lattice oxygens). It should be kept in mind, however, that these blue oxygens may correspond to H2O, OH-, or H3O+. For a detailed mechanistic discussion, see Section IV.1.

Row 3 shows top views of the same first-layer region. Row 4 overlays the trajectories of the first-layer surface Co atoms and the interfacial oxygen species (similar to Row 3 but from the top view) to highlight their in-plane mobility. For visual clarity, only a small interfacial portion of the slab is shown in each panel. In addition, the trajectory lines of the interfacial oxygens are plotted thinner than the Co trajectories, so that they remain distinguishable.

To complement the above structural analysis and resolve local water reactivity and bonding, we analyze two key radial distribution functions (RDFs): (i) the RDF between surface cobalt atoms and oxygen species from any interfacial adsorbed species (Co*∗*–Ow), where Ow includes oxygens from molecular water, hydroxyl groups, and transient hydronium ions; and (ii) the RDF between lattice surface oxygen atoms and hydrogen atoms (Os–H), which quantifies surface protonation via proton transfer from the interfacial layer.

It should be emphasized that, in the first RDF (Co*∗–Ow), the Co ions for the A-termination are Co2+, while for the B-termination they are Co3+*. Thus, comparing the A- and B-terminations also reflects how different cobalt oxidation states influence the adsorption geometry and structuring of the interfacial layer.

Co∗–Ow (Fig. 4a). On the B-termination, the Co*∗–Ow RDF shows a sharp first maximum at $r\approx 2.0~\text{\AA }$, followed by a deep minimum and a second maximum near $3.3\!-\!3.6~\text{\AA }$. The A-termination instead exhibits a broader, weaker feature shifted to larger distances ($\approx$3 Å), reflecting weaker and less specific adsorption of interfacial species. For reference, bulk Co–O distances for octahedral Co3+* and tetrahedral Co2+ are $\approx$1.95–1.96 ÅOmranpour and Behler (2024), bracketing the B-termination first maximum and confirming its inner-sphere character. Importantly, Ow here refers to all oxygen species in the interfacial layer (H2O, OH-, and H3O+), excluding lattice oxygens, i.e., not only molecular water alone.

Os–H (Fig. 4b). The Os–H RDF clearly separates covalently bound surface hydroxyls from hydrogen-bonded water. A sharp peak near $\approx$1.0 Å corresponds to surface OH groups formed via protonation by interfacial species, while a broader feature at 1.6–1.8 Å originates from hydrogen bonding with water molecules in the interfacial layer. The B-termination exhibits a much stronger OH peak and significantly more pronounced hydrogen bonding, demonstrating its higher reactivity toward water dissociation and the formation of a dense, oriented interfacial network. In contrast, the A-termination shows considerably weaker protonation and a more disordered hydrogen-bond network, reflecting its substantially lower activity for water dissociation.

IV.2.2 Discussion

The first peak of the Co*∗–Ow RDF (Fig. 4a), located at $\approx$2 Å, coincides with the equilibrium Co–O distance in bulk Co3O4 and can thus be attributed to oxygens of chemisorbed species coordinated to surface Co3+* ions on the B-termination.

But what does the second peak represent? In Row 4 of Fig. 3, two distinct mobilities of interfacial oxygen species (blue trajectories) can be identified: (i) oxygens located directly on top of surface Co atoms (purple) show reduced lateral motion, whereas (ii) oxygens positioned between these sites display significantly higher mobility. However, as visible in Rows 1 and 2 for the B-termination, both groups of oxygens belong to the same first interfacial layer at similar distances from the surface. The second peak can therefore be assigned to an outer-sphere (epitaxial) interfacial layer, consisting of oxygens laterally displaced relative to surface Co sites and hydrogen-bonded to chemisorbed OH groups or other interfacial species, rather than to a distinct, vertically separated second layer.

For the A-termination, the Co*∗–Ow RDF (Fig. 4a) shows its first feature at a larger distance, $\approx$ 3.0 Å, indicating weaker binding of interfacial oxygen species to surface Co2+* ions. This weaker interaction is also evident from the trajectory lines in Rows 2 and 4 of Fig. 3, which show much higher mobility of the interfacial oxygen species (both laterally and vertically) compared to the B-termination. Additionally, Row 3 of Fig. 3 shows a more disordered interfacial layer for the A-termination, with no clear adsorption registry relative to the underlying lattice.

In summary, a direct comparison reveals a pronounced difference between the two terminations:

B-termination: the first interfacial water layer organizes into a compact, laterally ordered adlayer with well-defined adsorption sites exhibiting clear registry with the surface lattice. In the top view, interfacial oxygen species cluster into recurring motifs with relatively uniform nearest-neighbor spacing, while in the side view the layer remains geometrically tight with reduced out-of-plane roughness.

A-termination: the first interfacial water layer is more diffuse and spatially heterogeneous. Adsorption sites are more broadly distributed without a persistent repeating pattern, and the layer exhibits larger height fluctuations.

These qualitative differences indicate a stronger templating effect and more specific binding geometries on the B-termination, whereas the A-termination exhibits a more liquid-like and weakly structured hydration environment. Moreover, the B-terminated surface also shows significantly higher degree of protonation (see Section IV.3).

IV.3 Density Profiles

The structural characteristics of interfacial water are analyzed by computing the number density profiles of oxygen and hydrogen atoms, as shown in Fig. 5. The density profiles $\rho_{i}(z)$ ($i=\mathrm{H},\mathrm{O}$) provide a spatially resolved view of how liquid water arranges itself at an oxide interface, capturing the contact layer, subsequent layering, and convergence to the bulk structure on the Å scale. For Co3O4(001) slabs exposed to water on both faces, these profiles reveal how the two crystallographic terminations, A and B, template the adjacent liquid differently through the height, position, and sharpness of the first and second maxima, as well as how far this ordering propagates into the liquid phase. Comparing oxygen and hydrogen profiles also provides an orientational cue: H peaks that are closer to the surface than O peaks indicate water molecules donating H-bonds toward surface acceptors. It is important to note that according to the definition used in these plots, if a surface oxygen is protonated, it is also counted as an OH group.

Overall features and assignment. Because the Co3O4(001) slab is exposed to water on both sides, the profiles in Fig. 5 exhibit two interfacial regions separated by the slab interior (gray-shaded region). Reading from left to right, the first interfacial zone corresponds to the B-terminated surface and the second to the A-terminated surface. At each interface, both $\rho_{\mathrm{O}}(z)$ and $\rho_{\mathrm{H}}(z)$ display a sharp contact-layer maximum followed by damped oscillations that decay into a flat bulk plateau, characteristic of packing frustration and layering at solid–liquid interfaces.

A–termination vs. B–termination. The main structural differences between the two terminations are revealed most clearly in the region of the second maxima, which dominate the interfacial water density. The B-termination exhibits a more pronounced and narrower second maximum, a deeper inter-layer minimum, and a slower decay toward the bulk plateau compared to the A-termination, reflecting stronger templating and tighter lateral ordering of the interfacial bilayer. The smaller first O peak visible at the B-side interface arises predominantly from protonated surface oxygens (surface OH groups), which are by definition counted as part of the “water” density. In contrast, the first peak on the A side is mostly molecular water, which explains why the apparent O peak height at the B interface is lower than at the A interface but does not contradict the stronger templating effect of the B-termination overall.

Hydrogen vs. oxygen: interfacial orientation and dissociation. For molecular water at both terminations, the H maxima lie slightly closer to the surface than the corresponding O maxima, consistent with preferential donation of H-bonds toward surface oxygens. On the B-termination, however, an additional, smaller H peak appears even closer to the surface than the main molecular H maximum. This peak originates from dissociated protons bound to surface oxygens (forming OH groups), as confirmed by the strong Os–H RDF peak at $\sim 1.0$ Å (Fig. 4). Thus, for the B-termination, the innermost hydrogen density partly reflects surface protonation rather than orientational ordering of intact water molecules. This resolves the apparent mismatch between O and H peak positions: the smaller first O peak is dominated by protonated surface species, while the main second O/H peaks correspond to intact molecular water.

Bulk region and normalization. Far from either surface, both $\rho_{\mathrm{O}}(z)$ and $\rho_{\mathrm{H}}(z)$ approach constant plateaus with the expected stoichiometric ratio $\rho_{\mathrm{H}}^{\mathrm{bulk}}\!\approx\!2\,\rho_{\mathrm{O}}^{\mathrm{bulk}}$ of bulk water. The flatness of these plateaus and the correct 2:1 ratio confirm proper normalization, adequate sampling, and negligible finite-size artifacts.

Temperature dependence (300 K vs. 400 K). Raising the temperature from $300$ K to $400$ K reduces the height and sharpness of the interfacial peaks, consistent with enhanced thermal fluctuations and increased mobility. Before the production runs, the systems were equilibrated in the NPT ensemble at their respective temperatures, so the simulation cell volumes are slightly different at 300 K and 400 K. Consequently, small differences in the bulk density plateaus are expected and consistent with the setup described in Section III. Importantly, the qualitative ordering persists: the B-termination maintains stronger interfacial structuring than the A-termination at both temperatures.

In summary:

Water at Co3O4(001) forms a structured interfacial bilayer at both terminations, followed by damped oscillations that converge to bulk density within a few Å. 2. 2.

The B-terminated surface supports a denser, more laterally ordered interfacial structure than the A-termination, as seen from its sharper secondary maxima, deeper minima, and slower decay. 3. 3.

The innermost O and H peaks on the B side originate mainly from protonated surface oxygens (OH groups), while the main O/H peaks farther from the surface correspond to intact molecular water donating H-bonds. Thus, the additional close-in H peak on the B side reflects surface hydroxylation rather than molecular orientation. 4. 4.

Increasing the temperature slightly reduces interfacial ordering but preserves bulk stoichiometry and the strong B > A contrast in templating.

Together with the discussion of the previous section, these density profiles confirm a consistent picture: the B-terminated Co3O4(001) surface forms a more ordered, tightly bound interfacial bilayer of partially dissociated water molecules that both hydroxylates and hydrogen-bonds to the surface, while the A-termination supports a looser and more weakly structured hydration layer. This ordered, quasi-epitaxial water adlayer on the B-termination has been proposed to enhance oxidative catalytic activity Omranpoor et al. (2022); Kox et al. (2020).

IV.4 Discussion on the Implications for Catalysis

As mentioned in the Introduction, the broader context of this work lies in catalysis in the aqueous phase in contact with solid surfaces and aims to elucidate a deeper understanding of such complex interfaces. Performing reactions such as alcohol oxidation in water can improve catalytic selectivity and enable catalyst reusability, both of which are central to practical catalysis. The interaction between water and transition-metal oxides (TMOs) such as Co3O4 is, however, more complex than at metal/water interfaces Henderson (2002); Björneholm et al. (2016). Water plays a multifaceted role, with both beneficial and detrimental consequences for oxidation chemistry Davies (2016); Lin et al. (2021). On the beneficial side, it can facilitate proton-coupled steps by providing mobile proton carriers (transient OH-, H3O+, and surface hydroxyls) that temporarily store and relay protons Omranpoor et al. (2022); Stuchebrukhov and Hammes-Schiffer (2010); Agmon et al. (2016); it has also been shown to promote O2 activation in certain oxidation regimes Muñoz-Santiburcio et al. (2018); Tran et al. (2016). On the detrimental side, water competes with reactants (e.g., alcohols) for adsorption sites and can also passivate lattice oxygen functionality Omranpoor et al. (2023), thereby diminishing surface reactivity. Which of these tendencies dominates is highly dependent on the reaction and conditions (pH, potential, temperature, coverage, defects, termination) to be investigated in future investigations.

This work specifically sought to bridge the gap between previous ab initio simulations and experimental observations and contributes to the understanding of Co3O4(001)–water interfaces. It confirmed several earlier findings while elucidating additional features of the Co3O4(001) A- and B-terminations, in particular the higher activity of the B-termination for water dissociation. This conclusion is directly supported by the Os–H RDF (Sec. IV.2), which shows a significantly stronger $\approx$1.0 Å peak on the B side, indicating more extensive protonation of surface oxygens compared to the A-termination. Furthermore, the Co∗–Ow RDF clearly reveals how the oxidation state of the surface cobalt ions affects their interaction with water: surface Co3+ ions on the B-termination exhibit stronger inner-sphere coordination with interfacial oxygens (sharp $\approx$2.0 Å peak), whereas surface Co2+ ions on the A-termination bind more weakly, resulting in a broader, more distant $\approx$3.0 Å feature. However, it should be emphasized that higher interfacial water structuring or enhanced hydroxylation does not necessarily translate into superior catalytic performance for all reactions; for example, in alcohol oxidation, many competing factors are at play. As our analysis showed (Sec. IV.1), multiple microscopic pathways operate simultaneously near the interface (e.g., proton transfer between waters, protonation/deprotonation of chemisorbed hydroxyls, and exchange between inner- and outer-sphere species), and introducing specific reactants only adds further complexity. These results reinforce the need for a more complete atomistic-level description of the Co3O4–water interface as a prerequisite for predictive understanding and optimization of aqueous-phase oxidation catalysis.

V Conclusions

In this work, we developed and applied a high-dimensional neural network potential to investigate the structure, dynamics, and reactivity of Co3O4(001)–water interfaces. By combining density functional theory reference data with active learning, we constructed a robust potential that enables nanosecond-scale molecular dynamics simulations for systems comprising thousands of atoms, thereby overcoming the spatial and temporal limitations of previous ab initio studies.

Our simulations reveal pronounced differences between the two possible surface terminations. The B-terminated Co3O4(001) surface forms a highly ordered, quasi-epitaxial hydration layer with strong templating effects, enhanced hydroxylation, and well-defined hydrogen-bonding networks. In contrast, the A-termination exhibits a more diffuse contact layer with weaker templating, lower hydroxylation, and less structured interfacial water. These structural differences are directly linked to the distinct reactivity of the terminations: the B-terminated surface shows a higher propensity for water dissociation, consistent with experimental observations of its greater catalytic activity.

Furthermore, the extended timescales captured in our simulations allow us to identify rare events that are inaccessible to shorter ab initio molecular dynamics runs. These include late-time proton transfer episodes within the outer hydration layer and dynamic protonation–deprotonation of surface hydroxyl groups, which play an essential role in interfacial charge and proton transport. Temperature-dependent simulations demonstrate that increasing the temperature reduces interfacial ordering but preserves the stronger templating and enhanced reactivity of the B-terminated surface.

Overall, this study provides an atomistically detailed and dynamically resolved picture of Co3O4(001)–water interfaces under realistic solvation conditions. The insights gained into interfacial structure, hydroxylation, and proton transfer mechanisms bridge the gap between experimental findings (particularly regarding the crucial role of Co3+) Falk et al. (2021) and previous ab initio molecular dynamics simulations Kox et al. (2020).

Acknowledgments

We are grateful for funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in TRR/CRC 247 (A10, project-ID 388390466) and under Germany’s Excellence Strategy – EXC 2033 RESOLV (project-ID 390677874). The authors also acknowledge the computing time provided to them by the Paderborn Center for Parallel Computing (PC2).

Data Availability

The data generated and analyzed during the current study are available from the corresponding authors upon request.