Modeling the aqueous interface of amorphous TiO₂ using deep potential molecular dynamics Free

Amorphous metal oxides (AMOs) exhibit properties that are useful in a plethora of energy and electrochemical applications.¹ Compared to their crystalline counterparts, AMOs often have better mechanical properties, stronger corrosion resistance, and more desirable defect distribution.^2,3 For example, when used as anodes in lithium ion batteries, amorphous Fe₂O₃ offers small volume expansion and extra active sites from defects and vacancies as well as an ideal matrix for fast ion diffusion.¹ Similar advantages have been observed when V₂O₅, a cathode material for high energy Mg batteries, is amorphized.⁴ Moreover, non-precious AMOs have been found to have favorable mechanical properties, thermal stability, and corrosion resistance that enhance the activity and durability of fuel cells.^5,6 Generally speaking, AMOs do not have periodicity in their structures and therefore lack long range order. Furthermore, they tend to have a broad distribution of coordination numbers. The unique structural features of AMOs affect the structure and dynamics of water at their aqueous interfaces, relevant to many applications, in intricate and often unknown ways.

Among the AMOs, amorphous titanium dioxide (a-TiO₂) is of particular interest here because TiO₂ is one of the most widely used materials in photocatalysis, solar cells, self-cleaning windows, and hydrogen storage.^7–9 The functional properties of crystalline TiO₂ are largely preserved in a-TiO₂, which makes it a cost-effective coating material for applications ranging from antibacterial surfaces and surface passivation layers in biosensors to protective layers of photo-anodes for water splitting.^10–12 Despite its important role in such applications, the structure of the surface and aqueous interface of a-TiO₂ is not well known. The scarcity of detailed experimental information has been exacerbated by the dependence of the structure on preparation methods,¹³ hindering the understanding of the behavior of a-TiO₂ in aqueous environment.

The relevance of a-TiO₂ has also motivated numerous computational studies. The structural properties of bulk a-TiO₂ have been investigated using both classical force fields^14–18 and ab initio molecular dynamics (AIMD) simulations^19–22 with the amorphous structure generated via the melt-and-quench method.²³ While significantly more reliable than classical force fields, AIMD has a high computational cost so that relatively small models prepared by extremely fast quenching from the melt could be examined. To generate more realistic models of a-TiO₂ with first principles accuracy, recently, Calegari Andrade et al. developed deep neural network potentials (DPs) that reproduced the structural properties of bulk crystalline and disordered TiO₂ phases predicted by density functional theory (DFT) calculations at a much lower computational cost.²⁴ Molecular dynamics simulations based on such DPs (DPMD) could thus explore the role of the system size and quenching rate and further predict the effect of pressure on the structure of a-TiO₂ in good agreement with experiment. However, the above studies have all focused on bulk a-TiO₂, whereas little is known about a-TiO₂ surfaces and their interaction with an aqueous environment. To our knowledge, only one computational study considered the surface of a-TiO₂ using classical MD for structural generation and DFT for characterization.²⁵ 

In this work, we investigate the structure and dynamics of interfacial water on a-TiO₂ surfaces using a combination of AIMD and DPMD simulations. Employing DFT data at the level of the SCAN functional,²⁶ we trained the DPs via an active learning protocol²⁷ that has been utilized to generate accurate DPs for crystalline TiO₂/water interfaces.^28–30 We constructed models of the a-TiO₂ surface using a variation of the multi-step procedure previously employed to generate amorphous SiO₂ surfaces,³¹ which involves cycles of heating and quenching the surface region of the sample while keeping the bulk (central) part fixed, followed by relaxing the outermost surface region in order to reduce the strain at the surface–bulk boundary. Such procedure was accomplished using nanosecond timescale DPMD simulations with potentials that accurately reproduced the structure of all bulk crystalline and disordered TiO₂ phases given by DFT-SCAN.²⁴ 

In more detail, we here examine two relatively small (162- and 216-atom) slab models of the a-TiO₂ surface in contact with water for which short (picosecond) timescale AIMD simulations are also feasible. Our results show that the distribution of water at the interface of a-TiO₂ lacks distinct peaks and can extend below the outermost TiO₂ atoms. Due to the lack of compact and ordered water layers, water at the surface of a-TiO₂ diffuses at a rate about ten times faster than that on crystalline TiO₂ surfaces. Bridging hydroxyls (Ti₂–O_bH with the oxygen coordinated to two Ti atoms) resulting from water dissociation are relatively stable and decay more slowly than terminal hydroxyls (Ti–O_wH; note that on TiO₂ there are no geminal Ti–(O_wH)₂ species such as on silica surfaces³²). These findings provide new insights into the properties and behavior of amorphous TiO₂ in aqueous environments and more broadly demonstrate a computational procedure for studying AMO–water interfaces.

DPs were trained using a data-driven active learning protocol.^27,29 The three core components of active learning are training a committee of DPs, exploration of configuration space by DPMD, and appending the training data using high-error snapshots from the exploration step. The last component involves performing static DFT calculations on the high-error snapshots and generating energies and forces; this is called the labeling step in machine learning language. This procedure iterates until convergence is reached when the standard deviation in energy and forces among the committee of DPs falls below a certain threshold. A committee of four DPs was trained, and the threshold for convergence was set to the average standard deviation in forces being within 0.08 eV/Å for any 1 ns DPMD simulation.

The first step of the aforementioned protocol was carried out using the deepMD-kit.³³ The cutoff radius and smooth cutoff radius for the chemical descriptor were set to 6 and 3 Å, respectively. About 20% of the training data was used for validation during neural network training. Technical details on the neural network architecture are published in a previous study.³⁴ DPMD simulations were performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package.³⁵ 

The training data consist of five types of systems: bulk water, bulk a-TiO₂, a-TiO₂ surface in vacuum, a-TiO₂ surface with one adsorbed water in vacuum, and a-TiO₂–water interface. The training data on bulk water were taken from a previous study³⁰ while the initial training data for the other four systems were acquired from sufficiently separated snapshots of AIMD trajectories sampled in the canonical ensemble (see in the following). Subsequent training data for the latter four systems were acquired through the iterative protocol described above. The final training data are summarized in Table I.

All ab initio calculations were performed with the CP2K package.³⁶ Static density functional theory (DFT) calculations were conducted using the SCAN functional.²⁶ Analytic Goedecker–Teter–Hunter (GTH) pseudopotentials³⁷ and a double ζ-with-polarization (DZVP) basis set were used together with a cutoff of 1200 Ry for the plane wave expansion of the density. Only the Γ point was considered for k-sampling. Structural optimizations were carried out through the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, where forces on all atoms were minimized to within 4.5 × 10⁻⁴ hartree/bohr.

The main purpose of AIMD simulations was to generate initial representative configurations to be labeled by the energy and forces obtained from static DFT-SCAN calculations. The Born–Oppenheimer AIMD simulations were thus performed using the PBE functional,³⁸ which is computationally faster than SCAN and performs quite well for crystalline rutile and anatase TiO₂.¹⁹ Grimme’s D3 van der Waals dispersion (vdW) corrections³⁹ were applied to aqueous interfaces. The AIMD equations of motion were integrated using the Verlet algorithm.⁴⁰ Hydrogen atoms were replaced by deuterium, and a time step of 0.5 fs was used. For equilibrium sampling, a canonical ensemble at 330 K was used for aqueous interfaces and at 300 K for other systems. A Nosé–Hoover thermostat with four chains^41,42 was employed.

To construct models of the a-TiO₂ surface, we started from (atomic coordinates and velocities of) a 1296-atom bulk a-TiO₂ sample generated in a recent DPMD study²⁴ and utilized a “cut-melt-and-quench” procedure that allowed us to generate a large number of bulk models of different sizes rather easily. In detail, the 1296-atom unit cell was first “cut” to form smaller bulk a-TiO₂ cubic unit cells of 162- and 216-atoms (two sizes used in previous computational studies of a-TiO₂^20,24), which are amenable to relatively fast DFT calculations during an active learning procedure. To allow removal/relaxation of the broken bonds at the unit cell boundaries, the cut models were heated to 2250 K, cooled to 300 K at a rate of 1 K/ps, and further equilibrated for 1 ns in an isothermal–isobaric ensemble using DPMD with the DPs generated in Ref. 24. This procedure resulted in 20 different bulk a-TiO₂ models, which were characterized structurally (radial distribution functions, coordination numbers, etc.). Two of these models, one of 162 and the other of 216 atoms, were then chosen as starting points for the surface construction based on their agreement with the bulk a-TiO₂ structure reported in previous studies.

Starting from such small bulk a-TiO₂ models, the a-TiO₂ surface was generated based on the scheme in Ref. 31. We first added a vacuum region of about 18 Å along the z-direction of the bulk a-TiO₂ unit cell, which resulted in the formation of periodically repeated slabs separated by a vacuum region. Next, the two sides of each slab were equilibrated at 300 K for 500 ps while keeping the central (∼1/5) part fixed. The surface layers were annealed to 1000 K and subsequently cooled to 300 K at a rate of 1 K/ps, followed by equilibration at 300 K for 500 ps, while keeping the central part fixed throughout the process. We used a quench rate of 1 K/ps because this was found to yield an a-TiO₂ structure nearly identical to that obtained with a slower rate of 0.1 K/ps in Ref. 24. To alleviate the strain built up at the boundary of mobile–immobile atoms, the frozen portion at the center of the slab was increased (to ∼1/3) and the two surfaces were further equilibrated at 300 K for 500 ps. (Note that increasing the frozen portion of the slab decreases the strain in the surface region because it creates the boundary of mobile–immobile atoms in a region that had previously been equilibrated.³¹) The above steps were conducted in a canonical ensemble using DPMD with the DPs generated in Ref. 24. Finally, the two surfaces were optimized geometrically using DFT-SCAN.

Two water/a-TiO₂ interface models were constructed by confining 123 and 122 water molecules, respectively, at the experimental density between the periodically repeated 162-atom and 216-atom slabs; these will be called 162-atom and 216-atom a-TiO₂/water interface models in the following. Due to the irregular profile of the a-TiO₂ surface, water molecules could diffuse below the uppermost TiO₂ atoms [see Fig. 2(b)]. AIMD in an isothermal–isobaric ensemble was then carried out to allow equilibration of water along the z-direction. The resulting simulation boxes had dimensions of 13.6 × 13.6 × 29.77 and 15.6 × 15.6 × 26.49 ų for the 162- and 216-atom models, respectively. Finally, AIMD in the canonical ensemble at 330 K was carried out in order to generate initial representative configurations of the a-TiO₂/water interface for training (and validating) the DPs.

Structures were analyzed using cutoff distances of 2.75 and 4.1 Å, respectively, for Ti–O and Ti–Ti bonds. For surface structures, only the top and bottom 1/3 of each slab were considered because the middle 1/3 was frozen during surface construction and thus retained bulk-like properties. The surface roughness was evaluated as the average absolute deviation of the z-coordinates of surface atoms, i.e., $Ra=1N∑i|zi−z̄|$⁠. Here, N is the number of outermost (Ti or O) atoms, z_i is the z-coordinate of a surface (Ti or O) atom, and $z̄$ is the corresponding average z-coordinate.

The convergence of a committee of DPs in the active learning protocol implies self-consistency but not necessarily accuracy. It is therefore crucial to benchmark certain key quantities computed using the trained DPs against those reported in the literature. One 1 ns DPMD trajectory was carried out for bulk water and bulk a-TiO₂ at 330 and 300 K, respectively. Using an arbitrary DP out of the committee of four is sufficient for benchmarking purpose because the estimation of energy and forces is consistent within the committee, therefore, the standard deviation in the benchmarked quantities is negligible.

Table II reports the positions of the peaks in the radial distribution functions (RDFs) for bulk water, obtained from a 1 ns DPMD simulation at 330 K. The agreement with DFT-SCAN results (between parentheses) from Ref. 43 is excellent. The computed diffusion coefficient of 0.181 ± 0.002 Ų/ps is also very close to the experimental value of 0.187 Ų/ps.⁴⁴ 

The structures of the 162- and 216-atom models of bulk a-TiO₂ generated with DPMD were characterized in terms of RDFs [Figs. 1(a)–1(c)], distribution of Ti–O–Ti angles [Fig. 1(d)], coordination numbers of Ti- and O atoms (Table III), and polyhedra connectivity. The computed density was 3.64 g/cm³, which falls in the range of 3.62–4.09 g/cm³ obtained experimentally.⁴⁵ The trained DPs also reproduce the experimentally measured pair correlation and the mean coordination numbers of sputtered a-TiO₂.⁴⁶ The split-peak feature in the Ti–Ti RDF [Fig. 1(a)] corresponds to edge- and vertex-sharing (Ti-centered) polyhedra.²¹ This medium range order is also observed in the distribution of Ti–O–Ti angles [Fig. 1(d)] where the sharper peak around 100° is attributed to edge-sharing polyhedra whereas the smaller peak around 130° comes from vertex-sharing polyhedra. The computed angle distribution is consistent with previous findings.^20,24 Another characteristic metric, the fractions of Ti-centered polyhedra of various connectivities are found to be 0.71, 0.26, and 0.03 for vertex- (V), edge- (E), and face-sharing polyhedra, respectively. The resulting V/E ratio is 2.7 to be compared to the values of 2.1²⁰ and 2.6²⁴ reported in previous studies. Altogether, the above benchmarking of structural properties indicates that the cut-melt-and-quench procedure driven by DPMD produces reasonably good bulk a-TiO₂ structures, validating the trained DPs.

To characterize the structure of the a-TiO₂ surface, we performed a 10 ps AIMD (also used for DP training) and 1 ns DPMD (four independent trajectories with the different DPs of the committee) simulations at 300 K of the 162-atom a-TiO₂ slab model described in Sec. II C. Table IV reports the distribution of coordination numbers of Ti and O atoms at the a-TiO₂ surface obtained from these simulations. At the surface, both Ti and O have a higher fraction of smaller coordination numbers than the bulk structure, as can be seen by comparing Table IV with Table III for bulk a-TiO₂. Furthermore, the fraction of vertex-sharing (Ti-centered) polyhedra decreases whereas that of edge-sharing polyhedra increases (Table V), resulting in a slightly lower V/E ratio. We can see small differences between the distributions of coordination numbers given by AIMD and DPMD, which are likely related to both the different durations and (to a smaller extent) the different underlying density functionals (PBE vs SCAN) of the two simulations. However, both the average Ti and O coordinations and the polyhedra connectivity obtained from the DPMD and AIMD simulations are very similar. The surface energy of the a-TiO₂ slab was estimated to be 1.16 J/m² relative to bulk a-TiO₂.

In comparison to the surfaces of crystalline TiO₂, the surfaces of a-TiO₂ appear irregularly corrugated. We thus evaluated the surface roughness in terms of the deviation from the average of the outermost atoms along the z-direction. As shown in Table VI, our DPMD-generated a-TiO₂ models exhibit a substantially higher degree of surface roughness in comparison to crystalline surfaces such as anatase(001), anatase(101), and rutile(110).

Preliminary insight into the structure of our a-TiO₂/water interface models was obtained from the last 22 of a 27 ps AIMD trajectory that started from the final configuration of a previous simulation in the isothermal–isobaric ensemble (see Sec. II C). For comparison, four independent 27 ps DPMD simulations of the a-TiO₂/water interface were also performed, starting from the same configuration as the AIMD simulation. The water density profiles of the 162-atom a-TiO₂/water interface model obtained from the two types of simulations are presented in Fig. 2. Unlike at a crystalline TiO₂ interface, water at the interface of a-TiO₂ lacks distinct peaks. Moreover, the distribution of water is asymmetrical because the two sides of the slab are different and extend slightly below the outermost TiO₂ atoms. The agreement between the density profiles predicted by the short time scale DPMD and AIMD simulations is quite good, further supporting the reliability of our DPs.

Analysis of the dynamics of interfacial water on a-TiO₂ shows that the water diffusion coefficient is about ten times larger than that of water at crystalline TiO₂ surfaces, such as rutile (110)³⁰ and anatase (101)⁴⁷ (Table VII). This result can be attributed to the irregular surface morphology of a-TiO₂. Due to the relatively high surface roughness, many intra- and inter-surface hydrogen bonds are indeed disrupted, leading to the lack of compact and ordered water layers at the a-TiO₂ surface described above.

To investigate the equilibrium properties of the a-TiO₂/water interface, four independent 10 ns DPMD trajectories were performed, the last 5 ns of which were used for analysis. Results for the structural and dynamical properties of the 162-atom a-TiO₂/water model are presented in Figs. 3 and 4. The water density profile [Fig. 3(a)] shows a better defined bulk-like region at the center than that obtained from the 27 ps trajectories in Fig. 2. In comparison to the bulk-like water region, the water density near the a-TiO₂ surface exhibits a significantly higher standard deviation among the DPs as illustrated by the shaded orange region in Fig. 3(a). This is different from previous DPMD studies of aqueous crystalline TiO₂ interfaces^28,30 where the interfacial water density showed a relatively small standard deviation among the DPs for simulations on the nanosecond scale. This difference suggests an intrinsic higher degree of fluctuation in the water density near the a-TiO₂ surface in comparison to the crystalline TiO₂ surfaces. The relatively low peaks of the water density also suggest that a-TiO₂ is less hydrophilic than the crystalline surfaces, a result consistent with contact angle measurements for water on a-TiO₂ nanotubes.⁴⁸ Nonetheless, water dissociates to a significant extent on the a-TiO₂ surface. As shown by the radial distribution functions in Figs. 3(b) and 3(c), there is indeed a strong peak around 0.97 Å in the O_s–H_w RDF indicative of water dissociation and a split peak around 1.86 and 2.15 Å in the O_w–Ti_s RDF corresponding to dissociative and molecular adsorption of water, respectively. Remarkably, the same features, although less pronounced, are also present in the RDFS obtained from the short AIMD trajectory.

The time evolution of the fraction of dissociated water is plotted in Fig. 4. We find that H₂O–OH proton exchange between adsorbed water and hydroxyl at adjacent sites occurs frequently with an average water dissociation fraction of ∼15%. Moreover, the survival probabilities of bridging and terminal hydroxyls show that the terminal hydroxyls (O_wH) decay at a faster rate in comparison to their bridging counterparts (O_bH). The survival probability for the bridging hydroxyls exhibits a piece-wise linear behavior. In contrast, terminal hydroxyls display a fast decay followed by a slower one with characteristic times of 16 and 150 ps, both shorter than the lifetime of terminal OD species on anatase (101) reported in Ref. 28.

The computed diffusion coefficients of interfacial water obtained from DPMD are similar to those given by AIMD (Table VII) with diffusion perpendicular to the interface slightly slower than that parallel to the interface.

To explore the sensitivity of the results for the a-TiO₂/water interface to the specific model used for the simulations, we also conducted one 37 ps AIMD simulation and four independent 10 ns DPMD simulations of the 216-atom a-TiO₂/water interface model (see the  Appendix). The aqueous interfaces of the 162-atom and 216-atom a-TiO₂ models share many similarities even though the latter has a higher proportion of under-coordinated surface atoms than the former. At both interfaces, the distribution of water is asymmetrical and extends below the uppermost surface atoms (see Figs. 3 and 5). The equilibrium water dissociation ratio averages to 19.3% for the 216-atom a-TiO₂–water interface. This is higher than the 14.9% fraction found at the 162-atom interface, consistent with the higher proportion of under-coordinated atoms on the 216-atom surface.

Similarly, proton exchange is observed more frequently at the 216-atom interface than at the 162-atom a-TiO₂–water interface (compare top panels of Figs. 6 and 4). This results in a faster decay of both bridging and terminal hydroxyls on the former interface (compare bottom panels of Figs. 6 and 4). We calculated the survival probability of the two types of hydroxyls using a 100 fs sampling frequency. We observe a double-exponential decay behavior for both terminal and bridging hydroxyls with characteristic times of 13 and 121 ps for the terminal vs 24 and 568 ps for the bridging hydroxyls. The lifetime of the terminal hydroxyls on the 216-atom interface is indeed slightly shorter than that on the 162-atom interface.

We trained a committee of deep potentials capable of describing the potential energy surface of amorphous TiO₂/water interface with the accuracy of DFT at the level of the SCAN functional via a data driven active learning protocol. Through benchmarking key quantities of bulk water and bulk a-TiO₂, we validated the accuracy of the converged deep potentials. Using both AIMD and DPMD simulations, we then investigated the structure and dynamics of the interfacial water on two different models of a-TiO₂ surface, one of which is not in the training data, illustrating the transferability of the DPs.

The two interface models are qualitatively similar, despite their different equilibrium water dissociation ratios (14.9 and 19.3%) that result from the difference in the proportion of under-coordinated surface atoms. We find that the distribution of water lacks symmetrical and distinct layers normally found at the aqueous interface of crystalline TiO₂. As a result of that, water diffuses about ten times faster on the surface of amorphous TiO₂ than on crystalline TiO₂. Bridging hydroxyls resulted from water dissociation decay significantly more slowly than terminal hydroxyls due to frequent proton exchanges between H₂O and terminal OH species adsorbed at adjacent Ti sites. These results for the structure of the a-TiO₂–water interface will be useful for a better understanding of various phenomena involving this interface. Moreover, we expect that the computational procedure of constructing a-TiO₂/water interfaces employed here will be broadly applicable to studying other AMO–water interfaces as well.

We thank Marcos Calegari Andrade and Bo Wen for helpful discussions. This work was supported by DoE-BES, Division of Chemical Sciences, Geosciences and Biosciences under Award No. DE-SC0007347. We also acknowledge support from the Computational Chemical Center: Chemistry in Solution and at Interfaces, funded by the DoE under Award No. DE-SC0019394. We used resources of the National Energy Research Scientific Computing Center (DoE Contract No. DE-AC02- 05cH11231). We also acknowledge use of the TIGRESS High Performance Computer Center at Princeton University.

The authors have no conflicts to disclose.

Zhutian Ding: Formal analysis (equal); Investigation (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Annabella Selloni: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available at https://github.com/CSIprinceton/amor-TiO2-water.

We performed four independent 1 ns DPMD simulations for the 216-atom a-TiO₂ surface in the canonical ensemble at 300 K and analyzed the last 500 ps. The coordination numbers of the surface atoms are reported in Table VIII. The committee averages of polyhedra connectivity are 0.77, 0.23, and 0.00 for vertex-, edge-, and face-sharing polyhedra, respectively, giving rise to a V/E ratio of 3.42.

Results for the aqueous interface are summarized in Table IX (diffusion coefficients) and Fig. 6 (fraction of dissociated water vs simulation time and lifetimes of bridging and terminal hydroxyls).