Free energy of proton transfer at the water–TiO₂ interface from ab initio deep potential molecular dynamics

Abstract

TiO₂ is a widely used photocatalyst in science and technology and its interface with water is important in fields ranging from geochemistry to biomedicine. Yet, it is still unclear whether water adsorbs in molecular or dissociated form on TiO₂ even for the case of well-defined crystalline surfaces. To address this issue, we simulated the TiO₂–water interface using molecular dynamics with an ab initio-based deep neural network potential. Our simulations show a dynamical equilibrium of molecular and dissociative adsorption of water on TiO₂. Water dissociates through a solvent-assisted concerted proton transfer to form a pair of short-lived hydroxyl groups on the TiO₂ surface. Molecular adsorption of water is ΔF = 8.0 ± 0.9 kJ mol⁻¹ lower in free energy than the dissociative adsorption, giving rise to a 5.6 ± 0.5% equilibrium water dissociation fraction at room temperature. Due to the relevance of surface hydroxyl groups to the surface chemistry of TiO₂, our model might be key to understanding phenomena ranging from surface functionalization to photocatalytic mechanisms.

---

1 Introduction

The photocatalytic activity of TiO₂ has attracted the attention of researchers for decades. Due to its natural abundance, chemical stability, and environmental compatibility, TiO₂ is one of the most widely used photocatalysts for scientific and technological applications to date.^1–3 Of particular relevance is the anatase phase of TiO₂, which predominates at the nanoscale. Since TiO₂ photocatalysis usually takes place in humid or aqueous environment, the interface of anatase TiO₂ with water is of fundamental importance, e.g. for elucidating the not yet fully understood mechanisms of photochemical water splitting,^4–6 UV-induced hydrophilicity,^7–9 and for improving the performance of TiO₂ nanomaterials in biomedical applications.^10,11 As a specific example, we focus here on the aqueous interface with the majority (101) surface of anatase. Our main goal is to solve an apparently simple but longstanding question: does water spontaneously dissociate on the defect-free anatase (101) surface? Hydroxyl groups are known to trap charges on TiO₂ surfaces,^12–14 hence the relevance of this particular question extends over the broad field of TiO₂ photocatalysis in aqueous environments.

Water adsorption on anatase (101) takes place at the under-coordinated five-fold Ti (Ti_5c) and two-fold oxygen (O_2c) surface sites (Fig. 1a).^15,16 Water oxygen atoms can form dative bonds with Ti_5c atoms (Fig. 1b), while surface O_2c's can either accept hydrogen bonds (Fig. 1b) or protons from water (Fig. 1c). Both Ti_5c and O_2c sites participate in water dissociation on the anatase (101) surface. As shown in Fig. 1c, water dissociation involves the transfer of a proton from a molecule adsorbed at a Ti_5c site to an O_2c. This results in the formation of two surface hydroxyl groups: a terminal hydroxyl on Ti_5c (OH_t, Fig. 1c, left) and a bridging hydroxyl, corresponding to a protonated O_2c atom (OH_br, Fig. 1c, right).

Fig. 1 (a) Surface undercoordinated Ti_5c and O_2c sites, the most relevant TiO₂ atoms for water adsorption on the anatase (101) surface. (b) Molecular adsorption of water on anatase TiO₂. (c) Dissociative water adsorption on anatase TiO₂. Terminal and bridging hydroxyls are indicated by OH_t and OH_br, respectively. (d) Model of the anatase (101)–water interface used to train the DNN potential. The water region is highlighted by the pale blue shading.

The presence/absence of dissociated water plays a major role in the surface chemistry of anatase TiO₂. Although both experiments and computer simulations agree on the nature of water adsorption on anatase (101) under vacuum conditions, the picture remains controversial at the interface of anatase (101) with liquid water. In vacuum, Temperature Programmed Desorption (TPD),¹⁷ Scanning Tunneling Microscopy¹⁶ (STM) and Density Functional Theory (DFT)^15,18 all agree that water dissociation can only take place at defect sites on the anatase (101) surface. At the aqueous interface, however, only few experiments can selectively probe water right at the TiO₂ surface, and the experimental data currently available for the anatase (101)–water interface were either obtained using reduced anatase,¹⁹ or could not unambiguously define the adsorption state of water.^20,21

Computational studies based on DFT have provided valuable insights into water adsorption on the surface of metals,^22,23 semi-conductors^24–26 and insulators.^27,28 When coupled with molecular dynamics (MD), DFT has also provided atomic-level interpretation of experimental surface sensitive vibrational spectrocopies.^29–31 DFT-based ab initio molecular dynamics (AIMD) employs forces derived from the quantum mechanical ground state of the electrons and could in principle be used to access the stability difference between molecular and dissociative adsorption of water on the anatase (101) surface. However, water dissociation on anatase TiO₂ is kinetically hindered and cannot be observed in the short-time dynamics currently accessible to AIMD. Neither the dissociation of water nor the recombination of hydroxyls were indeed observed in AIMD simulations of the anatase (101)–water interface of approximately 40 ps duration.²¹

In this work we extend both the time and length scales sampled by AIMD using a deep neural network (DNN) to represent the potential energy surface (PES), gaining a substantial speed-up of MD simulations. Several authors have demonstrated the DNNs' ability to reproduce the complex PES of ab initio potentials,^32–34 thus allowing accurate long-time simulations of several condensed phase systems^35–37 including metal oxide–water interfaces.^38,39 In the present study, the ab initio-based DNN potential reproduces energy and atomic forces obtained from DFT at a five orders of magnitude lower computational cost. This allows us to carry out MD simulations at the nanosecond timescale on systems of thousands of atoms and to use enhanced sampling techniques to obtain converged free energy surfaces of proton transfer at the aqueous anatase (101) interface as a function of suitable reaction coordinates. All our simulations are carried out assuming no net electric charge on the surface, a condition that can be experimentally realized at the point of zero charge and flat band potential of TiO₂ in contact with water. Altogether, the present results enable us to elucidate both the relative stability of molecular versus dissociated water and the mechanism of proton transport at the aqueous anatase (101) interface, which are two fundamental features for understanding the reactivity of this interface.

2 Methods

2.1 DNN training

Our DNN potential was constructed following the method recently proposed by Zhang et al.^34,40 In this scheme, the potential energy E of each atomic configuration is a sum of “atomic energies” , where E_i is determined by the local environment of atom i within a smooth cutoff radius R_c. E_i is constructed in two steps. First, for each atom a set of symmetry-preserving descriptors is constructed. Next, this information is given as input for a DNN, which returns E_i as the output. The additive form of E naturally preserves the extensive character of the potential energy. We refer the reader to ref. 34 and 40 for further details. With this method, the DNN potential has been shown to accurately reproduce interatomic forces and energies predicted by DFT.⁴¹

We trained our DNN potential using an active learning approach.⁴¹ In this scheme, training data is collected in an iterative scheme consisting of a computationally efficient exploration of the configurational space and evaluation of ab initio energies and atomic forces of a small set of selected atomic configurations. We performed such iterative procedure using three different systems: bulk anatase TiO₂ (108 atoms), bulk liquid water (192 atoms), and the anatase (101)–water interface (426 atoms), for which we included configurations with one or few interfacial water molecules either dissociated or near the transition state. The latter were generated with enhanced sampling techniques during the iterative refinement of the DNN and their inclusion has been essential to predict dissociation barriers close to those of DFT.

In more detail, our model of the anatase–water interface consists of 82 water molecules in contact with an anatase (101) slab (five TiO₂ layers exposing a (1 × 3) surface supercell) within a periodically repeated unit cell of size 10.4 Å × 11.4 Å × 36.8 Å along the three orthogonal directions [1̄01], [010] and [101] of the anatase TiO₂ crystal lattice (Fig. 1d). Forces and energies were computed with the SCAN functional,⁴² as implemented in the PWscf code of Quantum ESPRESSO.^43,44 The SCAN functional has been shown to accurately predict both the structure of water⁴⁵ and TiO₂.²¹ We performed all DNN potential-based MD simulations using the Lammps code⁴⁶ interfaced with the DeePMD-kit.⁴⁷ Enhanced sampling techniques were provided by PLUMED.⁴⁸ The DeePMD-kit code⁴⁷ was used to train the DNN potential. Additional details on the DNN training method and DFT calculations can be found in the ESI.† MD performed with the DNN potential is called DPMD in the following.

2.2 DPMD simulations

DPMD simulations of the anatase (101)–water interface sampled the canonical ensemble (constant volume and temperature), while the isobaric–isothermal ensemble (constant pressure and temperature) was sampled for liquid water simulations. Liquid water simulations were performed with a periodically repeated cell containing 64 water molecules. For the anatase (101)–water system, comparison between DPMD and AIMD involved simulations of a 1 × 3 anatase (101) surface supercell exposing 6Ti_5c and 6O_2c sites. This supercell was used also for umbrella sampling simulations used to evaluate the free energy barrier to transfer a proton from water to a surface O_2c atom. All other DPMD simulations of the anatase (101)–water interface were carried out with a larger 4 × 12 anatase (101) surface supercell exposing 96Ti_5c and 96O_2c sites. In all cases, the system was composed by a slab of 5TiO₂ layers in contact with a 20 Å water slab at the experimental density; essentially identical interfacial water properties were obtained in a simulation using a water thickness of 40 Å (see ESI†). Temperature and pressure were controlled by a Nosé–Hoover thermostat^49,50 at 330 K and Parrinello–Rahman barostat⁵¹ under 1 bar, respectively. A temperature 30 K higher than room temperature was used to roughly mimic nuclear quantum effects on the oxygen distribution in water.⁵² The classical equations of motion were numerically integrated using Verlet's method⁵³ with a time step of 0.5 fs for D₂O (0.25 fs for H₂O). Equilibrium properties were computed from simulations with D₂O, since it allows a larger time step in the numerical integration of the equations of motion without affecting the equilibrium statistical properties of the investigated systems. All error bars were evaluated from simulations ran with independently trained DNN potentials.

3 Results and discussion

3.1 Validation

The ability of DPMD to reproduce the DFT-predicted structures of pure water, pure anatase TiO₂, and anatase (101)–water interface is confirmed by the following results. (i) The anatase TiO₂ lattice constants given by DPMD agree well with the values predicted by SCAN²¹ (DPMD: a = 3.80 Å, c = 9.52 Å; SCAN: a = 3.77 Å, c = 9.52 Å). (ii) The radial distribution functions and average water density profile obtained from 4 independent 40 ps DPMD simulations of the anatase (101)–water interface agree, within their standard deviation, with analogous quantities obtained from a 40 ps AIMD trajectory of this interface (Fig. 2a and b, respectively). (iii) Our DNN potential reproduces the oxygen radial distribution function and density of bulk liquid water obtained from SCAN AIMD simulations, which are in good agreement with experiment⁴⁵ (see Fig. 2b, inset). DPMD also yields a water diffusion coefficient that compares well with that given by AIMD (DPMD: 0.17 ± 0.01 Ų ps⁻¹, AIMD:⁴⁵ 0.19 ± 0.03 Ų ps⁻¹). (iv) The water density profile along the direction perpendicular to the TiO₂ surface derived from DPMD simulations of our large (4 × 12) anatase (101) surface model (surface area of 18.5 nm²) essentially coincides with that obtained from analogous simulations of the smaller (1 × 3) anatase (101) model (surface area of 1.16 nm²) used in our training data (Fig. 2c).

Fig. 2 (a) Radial distribution function of selected atomic type pairs at the anatase (101)–water interface given from DPMD (lines) and AIMD (points). The definition of Ti_5c and O_2c is given in Fig. 1a, and O_w represents water oxygen atoms. (b) Density profile of water confined between two TiO₂ surfaces as predicted by DPMD and AIMD.²¹ In the inset of (b) we compare the oxygen radial distribution function in water, g(r), as obtained from AIMD⁴⁵ and DPMD. (c) Comparison of the water density profiles obtained from DPMD simulations of water in contact with 1 × 3 (1.16 nm²) and 4 × 12 (18.5 nm²) anatase (101) surface supercells. Distributions in (a) and (b) were obtained from 40 ps trajectories, while 2.5 ns of dynamics was sampled in (c). Shaded areas indicate the standard deviation obtained from four independent DPMD simulations.

Additional validation tests are reported in the ESI,† where we compare results from DPMD and AIMD for the vibrational densities of states of water and TiO₂ at the anatase–water interface and water dissociation at both the dry anatase surface and the aqueous interface. In particular, from Fig. S4 of the ESI,† one can see that the DFT result for the work to move a H⁺ from a surface O_2c to a nearby OH_t is reproduced very well by DPMD.

3.2 Equilibrium sampling of aqueous anatase (101) via DPMD

The equilibrium structure of the aqueous anatase (101) interface was determined by averaging the results of 4 independent DPMD trajectories of 25 ns duration each, from which error bars of computed quantities were also estimated. All simulations started from configurations with purely molecular water at the interface and a surface model exposing 96Ti_5c and 96O_2c sites was used, which allows a minimum non-zero hydroxyl surface coverage of 1.04%. The water density profile predicted by such DPMD simulations is perfectly symmetric with respect to its center (Fig. 2c), consistent with the symmetry of the anatase (101) slab used for the simulations. Since the water density profile does not depend on the surface supercell size (Fig. 2c), the lack of symmetry of the profile in Fig. 2b can be attributed to the short duration of that trajectory. In fact, AIMD simulations have shown a substantial slowdown of the rotational and translational dynamics of water close to the anatase (101) surface relative to bulk liquid water.²¹ Hence, a longer simulation time is needed to obtain the equilibrium density distribution of interfacial water on the anatase (101) surface relative to bulk liquid water.

Analysis of the DPMD trajectories shows that the equilibrium adsorption configuration of water at the aqueous anatase (101) interface consists mostly of intact molecules at Ti_5c sites, but also includes an average 5.6 ± 0.5% of water dissociated into OH_t and OH_br hydroxyl groups (see Fig. 3a; note that the average surface hydroxyl coverage is calculated excluding the first 10 ns of this curve). Hydroxyls have a finite lifetime before they recombine to form molecular water adsorbed at Ti_5c sites. This can be important in photo-oxidation processes, e.g. photocatalytic water splitting, where photoexcited holes are transferred from TiO₂ to adsorbed hydroxyls to form OH radicals.^54–56 The hydroxyl lifetime can be estimated from its survival probability distribution P(t), which is the probability to continuously observe water dissociated during a time interval t. We computed P(t) for OH_t and OD_t separately (Fig. 3b), and defined the average lifetime (τ) as . We found the average lifetime of OD_t (τ = 0.6 ± 0.2 ns) to be twice as long as the one of OH_t (τ = 0.3 ± 0.1 ns). For comparison, the hydroxyl lifetimes observed in our simulation are shorter than the lifetimes of several ns for electron–hole recombination in anatase,^57,58 and on the same timescale as the characteristic times of hole transfer to organic adsorbates on TiO₂ surfaces.³

Fig. 3 (a) Time evolution of the surface hydroxyl coverage on the anatase (101)–water interface, as obtained from an average of 4 DPMD simulations. (b) Un-normalized survival probability of hydroxyl groups on the anatase (101)–water (H₂O and D₂O) interface as a function of time. In both plots τ denotes the average lifetime of terminal hydroxyl groups. (c) Mechanisms of proton transfer reaction at the anatase (101)–water interface. I: molecular adsorption; II: transition state; III and IV: equivalent configurations of the dissociated water. Water dissociation mechanism follows the path: I → II → III (or IV), while proton transport follows: III → II → IV or IV → II → III.

The mechanism of hydroxyl formation that emerges from our simulations involves a concerted proton transfer from molecular water to surface O_2c sites. As shown in Fig. 3c (path I → II → III or I → II → IV), the initial state is a stable hydrogen-bond complex, in which the path from the water coordinated to a Ti_5c site to two adjacent O_2c sites is fully connected with hydrogen bonds. The proton transfer proceeds through a Grotthuss-like mechanism,⁵⁹ with a transient hydronium ion formed along the reaction path. The dissociation reaction produces a terminal (bonded to Ti_5c) and a bridging hydroxyl (protonated O_2c) connected through an intermediate water molecule that is H-bonded to both of them. The TiO₂ surface remains charge neutral along the proton transfer reaction, and no proton transfer to liquid water is observed.

A similar mechanism, involving an intermediate water molecule, is observed also for the proton transport along O_2c sites (path III → II → IV in Fig. 3c). A proton at a bridging O_2c site can either recombine with a terminal hydroxyl to form molecular water or transfer to a neighboring non-protonated O_2c site. The latter proton transport process requires a fully connected hydrogen-bond path between initial and final states similar to that of water dissociation/recombination, and has also the same transition state. For this reason, the recombination of bridging and terminal hydroxyls to form molecular water competes with the proton transfer from one bridging hydroxyl to a non-protonated O_2c site.

3.3 Enhanced sampling

We estimated the free energy difference between molecular and dissociated water on aqueous TiO₂ from enhanced sampling techniques. Due to the complex nature of proton transfer reactions in water,^60–62 we devised a two-step procedure to reconstruct the free energy surface of water dissociation on titania. First, we estimated the free energy barrier to transfer a proton from water to a surface O_2c atom using umbrella sampling.⁶³ With this method, the free energy F(S_o) of the system was obtained as a function of the minimum distance S_o between a particular surface O_2c (index o) and any hydrogen atom. A continuous version of the minimum distance was used, (λ = 500 Å), and quadratic (or umbrella) potentials were applied at different values of S_o, allowing the sampling of configurational space regions that would be rarely visited from unbiased simulations. The free energy surfaces obtained from umbrella sampling are given in the ESI.†

The second step in our procedure involved a new set of enhanced sampling simulations with a bias potential V = Σ_o − F(S_o) applied to the O_2c sites in the system, in which F(S_o) is the free energy previously obtained from umbrella sampling. This set of bias potentials enhances the probability of proton transfer to every O_2c on the anatase (101) surface. A probability distribution of the two collective variables given in Fig. 4a was then constructed in the biased ensemble and re-weighted by the Boltzmann factor e^βV, where β is the inverse temperature. The probability distribution was obtained from the average of 12 independent 3 ns long simulations (viz. 3 simulations with each of the four independent DNN potentials) and the average over all water adsorption sites in the system. The choice of the two collective variables in Fig. 4a was inspired by earlier computer simulations of proton diffusion in basic and acidic aqueous solutions.⁶⁴

Fig. 4 (a) Interatomic distances used to define the two collective variables (CVs) describing proton transfer reactions. The first CV, d_ooo = (d₁ + d₂)/2, represents the average distance between hydrogen-bonded oxygen atoms connecting proton-donor and proton-acceptor sites. The second CV, (v₁ + v₂)/2 (v_i = b_i − h_i), describes the progress of the proton transfer reaction, with positive and negative values corresponding to the dissociated and the molecular states, respectively. (b) Free energy surface of water dissociation (left) and proton transport (right) at the anatase (101)–water interface. Roman numerals indicate the adsorption state of water at specific regions in the free energy plot, as shown in Fig. 3c. Molecular water is more stable than the dissociated state by 8.0 ± 0.9 kJ mol⁻¹, with a free energy barrier of 32 ± 4 kJ mol⁻¹ separating these states.

From the computed free energy surface in Fig. 4b (left), two main features emerge: (i) water adsorbed in molecular form is thermodynamically more stable than dissociatively adsorbed water on the anatase (101) surface by ΔF = 8.0 ± 0.9 kJ mol⁻¹; (ii) the free energy barrier for water dissociation, ΔF^‡_d = 32 ± 4 kJ mol⁻¹, is around 13 times larger than the thermal energy at room temperature and would unlikely be crossed during the typical timescale (tens of ps) of AIMD simulations. Analysis of the computed trajectories in the presence of the bias potential confirms the mechanism of water dissociation observed in the unbiased simulations: water dissociation always involves a water molecule mediating the proton transfer between water adsorbed on a Ti_5c site and a neighboring O_2c atom (Fig. 3c). The enhanced sampling simulations also captured the mechanism of proton transport at the aqueous anatase (101) interface observed in the unbiased simulations. As shown in Fig. 4b (right), proton transport produces an approximately symmetric free energy surface with a free energy barrier of 25 ± 4 kJ mol⁻¹ separating the two equivalent states.

To verify the consistency of the free energy surfaces in Fig. 4b with the results of our unbiased DPMD simulations (Fig. 3), we can estimate the equilibrium fraction of dissociated water on the anatase (101)–water interface from the Boltzmann factor e^−βΔF, where β is the inverse temperature and ΔF is the free energy difference between molecular and dissociated water. Using the computed value ΔF = 8.0 kJ mol⁻¹, we obtain a dissociation fraction of 5.0% at our simulation temperature of 330 K, in very good agreement with the DPMD simulations. This result is also confirmed by Monte Carlo sampling of a 2-D lattice gas model that includes the constraint present in the real system: adsorbed water can only dissociate if at least one neighboring O_2c site is non-protonated (see ESI†). Using the lattice gas model, we can further estimate the free energy of the anatase (101)–water interface as a function of the surface hydroxyl coverage, which allows us to rule out the existence of additional metastable states at other surface hydroxyl fractions.

Another relevant comparison is that between the hydroxyl lifetime determined from our equilibrium simulations and that estimated from the computed free energy barrier. Using transition-state theory, the recombination rate constant can be written as , where ν_c is the attempt frequency and ΔF^‡_r = ΔF^‡_d − ΔF is the barrier for a terminal and bridging hydroxyl to recombine into a water molecule. Approximating ν_c with the H-bonded OD or OH stretching frequency of heavy (ν_c ≈ 71 THz) and light (ν_c ≈ 101 THz) water, and taking the recombination barrier ΔF^‡_r = 25 kJ mol⁻¹, we obtain a lifetime τ = 1/k of the terminal hydroxyl group equal to 0.09 ns and 0.13 ns for OH and OD, respectively. Thus, this simplified approach predicts lifetimes on the same order of magnitude as the ones evaluated from our DPMD simulations, although with slightly underestimated values. However, the fact that the simulations give τ_OD ∼ 2τ_OH, rather than τ_OD ∼ 1.4τ_OH as in the simplified approach, indicates a process more complex than a simple exponential relaxation dynamics. This is apparent from the curves in Fig. 3b, that cannot be fitted by a single exponential.

4 Conclusions

In summary, the ability of neural networks to reproduce the ab initio potential energy surface of the anatase (101)–water interface together with their scalability has allowed us to substantially extend the size and time scales of our first principles simulations for this system. By combining the computationally efficient neural network potential with enhanced sampling techniques, we have also been able to reconstruct the free energy surface of water dissociation and proton transport at the anatase interface with liquid water. The present simulations support a picture of predominantly molecular adsorption on the defect-free anatase (101) in contact with liquid water at ambient conditions, with an entropically driven 5.6% fraction of water dissociation at the interface. They also provide detailed mechanisms of proton transfer at the interface, free energy barriers and lifetimes that are all features to be taken into account when investigating the reactivity of the TiO₂–water interface, either chemical or photo-induced, which is key to many practical applications.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was conducted within the Computational Chemical Center: Chemistry in Solution and at Interfaces funded by the DoE under Award DE-SC0019394. A. S. acknowledge the support of DoE-BES, Division of Chemical Sciences, Geosciences and Biosciences under Award DE-SC0007347. M. C. A. acknowledges partial financial support from the CNPq-Brazil. 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.

Notes and references

1. A. Fujishima, X. Zhang and D. A. Tryk, Surf. Sci. Rep., 2008, 63, 515–582.
2. A. L. Linsebigler, G. Lu and J. T. Yates, Chem. Rev., 1995, 95, 735–758.
3. J. Schneider, M. Matsuoka, M. Takeuchi, J. Zhang, Y. Horiuchi, M. Anpo and D. W. Bahnemann, Chem. Rev., 2014, 114, 9919–9986.
4. A. Fujishima and K. Honda, Nature, 1972, 238, 37–38.
5. S. U. M. Khan, M. Al-Shahry and W. B. Ingler Jr, Science, 2002, 297, 2243–2245.
6. M. Gratzel, Acc. Chem. Res., 1981, 14, 376–384.
7. R. Wang, K. Hashimoto, A. Fujishima, M. Chikuni, E. Kojima, A. Kitamura, M. Shimohigoshi and T. Watanabe, Nature, 1997, 388, 431–432.
8. L. Zhang, R. Dillert, D. Bahnemann and M. Vormoor, Energy Environ. Sci., 2012, 5, 7491.
9. T. Zubkoy, D. Stahl, T. L. Thompson, D. Panayotov, O. Diwald and J. T. Yates, J. Phys. Chem. B, 2005, 109, 15454–15462.
10. A. Fujishima, T. N. Rao and D. A. Tryk, J. Photochem. Photobiol., C, 2000, 1, 1–21.
11. R. Cai, K. Hashimoto, K. Itoh, Y. Kubota and A. Fujishima, Bull. Chem. Soc. Jpn., 1991, 64, 1268–1273.
12. M. Anpo, T. Shima and Y. Kubokawa, Chem. Lett., 1985, 14, 1799–1802.
13. O. I. Micic, Y. Zhang, K. R. Cromack, A. D. Trifunac and M. C. Thurnauer, J. Phys. Chem., 1993, 97, 7277–7283.
14. S. H. Szczepankiewicz, A. J. Colussi and M. R. Hoffmann, J. Phys. Chem. B, 2000, 104, 9842–9850.
15. A. Vittadini, A. Selloni, F. P. Rotzinger and M. Grätzel, Phys. Rev. Lett., 1998, 81, 2954–2957.
16. Y. He, A. Tilocca, O. Dulub, A. Selloni and U. Diebold, Nat. Mater., 2009, 8, 585–589.
17. G. S. Herman, Z. Dohnálek, N. Ruzycki and U. Diebold, J. Phys. Chem. B, 2003, 107, 2788–2795.
18. A. Tilocca and A. Selloni, Langmuir, 2004, 20, 8379–8384.
19. I. M. Nadeem, J. P. W. Treacy, S. Selcuk, X. Torrelles, H. Hussain, A. Wilson, D. C. Grinter, G. Cabailh, O. Bikondoa, C. Nicklin, A. Selloni, J. Zegenhagen, R. Lindsay and G. Thornton, J. Phys. Chem. Lett., 2018, 9, 3131–3136.
20. S. Hosseinpour, F. Tang, F. Wang, R. A. Livingstone, S. J. Schlegel, T. Ohto, M. Bonn, Y. Nagata and E. H. Backus, J. Phys. Chem. Lett., 2017, 8, 2195–2199.
21. M. F. Calegari Andrade, H.-Y. Ko, R. Car and A. Selloni, J. Phys. Chem. Lett., 2018, 9, 6716–6721.
22. J. Cerdá, A. Michaelides, M. L. Bocquet, P. J. Feibelman, T. Mitsui, M. Rose, E. Fomin and M. Salmeron, Phys. Rev. Lett., 2004, 93, 116101.
23. J. Carrasco, B. Santra, J. Klimeš and A. Michaelides, Phys. Rev. Lett., 2011, 106, 026101.
24. B. C. Wood, E. Schwegler, W. I. Choi and T. Ogitsu, J. Am. Chem. Soc., 2013, 135, 15774–15783.
25. T. A. Pham, Y. Ping and G. Galli, Nat. Mater., 2017, 16, 401–408.
26. P. Liao, J. a. Keith and E. a. Carter, J. Am. Chem. Soc., 2012, 134, 13296–13309.
27. D. Ceresoli, M. Bernasconi, S. Iarlori, M. Parrinello and E. Tosatti, Phys. Rev. Lett., 2000, 84, 3887–3890.
28. G. Melani, Y. Nagata, J. Wirth and P. Saalfrank, J. Chem. Phys., 2018, 149, 014707.
29. M. Sulpizi, M. Salanne, M. Sprik and M.-P. Gaigeot, J. Phys. Chem. Lett., 2013, 4, 83–87.
30. R. Khatib, E. H. G. Backus, M. Bonn, M.-J. Perez-Haro, M.-P. Gaigeot and M. Sulpizi, Water orientation and hydrogen-bond structure at the water–fluorite interface, Nature Publishing Group, 2016, vol. 6, p. 24287.
31. M. Dellostritto, S. M. Piontek, M. L. Klein and E. Borguet, J. Phys. Chem. C, 2018, 122, 21284–21294.
32. J. Behler and M. Parrinello, Phys. Rev. Lett., 2007, 98, 146401.
33. S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt and K.-R. Müller, Sci. Adv., 2017, 3, e1603015.
34. L. Zhang, J. Han, H. Wang, R. Car and W. E, Phys. Rev. Lett., 2018, 120, 143001.
35. B. Cheng, E. A. Engel, J. Behler, C. Dellago and M. Ceriotti, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1110–1115.
36. L. Bonati and M. Parrinello, Phys. Rev. Lett., 2018, 121, 265701.
37. H.-Y. Ko, L. Zhang, B. Santra, H. Wang, W. E, R. A. DiStasio and R. Car, Mol. Phys., 2019, 117, 3269–3281.
38. V. Quaranta, M. Hellström and J. Behler, J. Phys. Chem. Lett., 2017, 8, 1476–1483.
39. M. Hellström, V. Quaranta and J. Behler, Chem. Sci., 2019, 10, 1232–1243.
40. L. Zhang, J. Han, H. Wang, W. A. Saidi, R. Car and W. E, Adv. Neural Inf. Process. Syst., 2018, 31, 4436–4446.
41. L. Zhang, D.-Y. Lin, H. Wang, R. Car and W. E, Phys. Rev. Mater., 2019, 3, 023804.
42. J. Sun, A. Ruzsinszky and J. P. Perdew, Phys. Rev. Lett., 2015, 115, 036402.
43. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
44. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. De Gironcoli, P. Delugas, R. A. Distasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H. Y. Ko, A. Kokalj, E. Kücükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H. V. Nguyen, A. Otero-De-La-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
45. M. Chen, H.-Y. Ko, R. Remsing, M. F. Calegari Andrade, B. Santra, Z. Sun, A. Selloni, R. Car, M. Klein, J. Perdew and X. Wu, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 10846–10851.
46. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
47. H. Wang, L. Zhang, J. Han and W. E, Comput. Phys. Commun., 2018, 228, 178–184.
48. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
49. S. Nosé, Mol. Phys., 1984, 52, 255–268.
50. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
51. M. Parrinello, J. Appl. Phys., 1981, 52, 7182–7190.
52. J. A. Morrone and R. Car, Phys. Rev. Lett., 2008, 101, 1–4.
53. L. Verlet, Phys. Rev., 1967, 159, 98–103.
54. J. Cheng, M. Sulpizi, J. Vandevondele and M. Sprik, ChemCatChem, 2012, 4, 636–640.
55. J. Chen, Y.-F. Li, P. Sit and A. Selloni, J. Am. Chem. Soc., 2013, 135, 18774–18777.
56. J. Cheng, J. Vandevondele and M. Sprik, J. Phys. Chem. C, 2014, 118, 5437–5444.
57. J. Tang, J. R. Durrant and D. R. Klug, J. Am. Chem. Soc., 2008, 130, 13885–13891.
58. M. Xu, Y. Gao, E. M. Moreno, M. Kunst, M. Muhler, Y. Wang, H. Idriss and C. Wöll, Phys. Rev. Lett., 2011, 106, 1–4.
59. D. Marx, ChemPhysChem, 2006, 7, 1848–1870.
60. P. L. Geissler, C. Dellago, D. Chandler, J. Hutter and M. Parrinello, Science, 2001, 291, 2121–2124.
61. A. Hassanali, M. K. Prakash, H. Eshet and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 20410–20415.
62. M. Chen, L. Zheng, B. Santra, H. Y. Ko, R. A. Distasio, M. L. Klein, R. Car and X. Wu, Nat. Chem., 2018, 10, 413–419.
63. G. Torrie and J. Valleau, J. Comput. Phys., 1977, 23, 187–199.
64. A. Hassanali, F. Giberti, J. Cuny, T. D. Kühne and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 13723–13728.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sc05116c

---