Theoretical analysis of zirconium oxynitride/water interface using
  neural network potential
Zr oxides and oxynitrides are promising candidates to replace precious metal
cathodes in polymer electrolyte fuel cells. Oxygen reduction reaction activity
in this class of materials has been correlated with the amount of oxygen
vacancies, but a microscopic understanding of this correlation is still
lacking. To address this, we simulate a defective Zr$_7$O$_8$N$_4$/H$_2$O
interface model and compare it with a pristine ZrO$_2$/H$_2$O interface model.
First, ab initio replica exchange Monte Carlo sampling was performed to
determine defect segregation at the surface in the oxynitride slab model, then
molecular dynamics accelerated by neural network potentials was used to perform
1000 of 500 ps-long simulations to attain sufficient statistical accuracy of
the solid/liquid interface structure. The presence of oxygen vacancies on the
surface was found to clearly modify the local adsorption structure: water
molecules were found to adsorb preferentially on Zr atoms surrounding oxygen
vacancies, but not on the oxygen vacancies themselves. The fact that oxygen
vacancy sites are free from poisoning by water molecules may explain the
activity enhancement in defective systems. The layering of water molecules was
also modified considerably, which should influence the proton and O$_2$
transport near the interfaces which is another parameter that determines the
overall activity.
Akitaka Nakanishi, Shusuke Kasamatsu, Jun Haruyama, Osamu Sugino
# Structural analysis of zirconium oxynitride/water interface using neural network potential

###### Abstract

Zr oxides with oxygen-nitrogen substitutions and oxygen vacancies are promising candidates to replace Pt as electrocatalyst for the oxygen reduction reaction. To understand the microscopic structure of the catalyst/water interface, many nanosecond-long molecular dynamics simulations were performed using the interatomic force field constructed by machine learning the ab initio calculations. A defective Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\)/H\({}_{2}\)O interface model was simulated and compared with a pristine ZrO\({}_{2}\)/H\({}_{2}\)O interface model. Water was found to be adsorbed on the surface partly as a H\({}_{2}\)O molecule or as dissociated components, OH and H, on both surfaces. On the pristine ZrO\({}_{2}\) surface, H\({}_{2}\)O molecules show a monolayer adsorption structure and adsorb on all Zr atoms with equal probability. On defective Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) surface, however, H\({}_{2}\)O molecules show a bilayer adsorption structure and do not adsorb on the oxygen vacancies but only on some of the surrounding Zr atoms. The total number of the molecular and dissociative adsorptions are almost the same on both surfaces, but a larger number are dissociatively adsorbed on the pristine surface. Possible implications for the catalytic behaviour are discussed.
## Introduction

Polymer electrolyte fuel cells are attracting attention as energy conversion devices, and their performance depends on the catalytic activity of the oxygen reduction reaction (ORR) taking place on the cathode side of the electrode. Pt has been alloyed or coated to increase its catalytic activity. Although Pt has high catalytic activity, it has problems in terms of durability and cost.[1]

Zr oxides[2, 3, 4] and Ti oxides[5, 6, 7, 8] are expected to replace Pt as electrocatalysts due to their high durability, low cost and high catalytic activity obtained by introducing defects such as oxygen vacancies (V\({}_{\mathrm{O}}\)s) and oxygen-nitrogen substitutions (N\({}_{\mathrm{O}}\)s).[9, 10] Previous studies have shown that the catalytic activity correlates with the amount of V\({}_{\mathrm{O}}\)s and other factors in ZrO\({}_{2}\),[2, 11, 12] but the mechanisms of the reaction pathways, rate-determining steps and other aspects of the ORR remain unknown. To the best of our knowledge, there have not been many studies on the ORR activity of defective oxide surfaces compared to pristine ones.[13, 14, 15]

Muhammady et al. studied the activity-energy diagrams of ORR on tetragonal ZrO\({}_{2}\)(101) surfaces with V\({}_{\mathrm{O}}\) and N\({}_{\mathrm{O}}\) introduced from ab initio calculations.[15] The results suggest that the Zr atoms and the V\({}_{\mathrm{O}}\) do not differ significantly in activity. Their model had neither solvent nor adsorbed water molecules. In addition, their studies discussed only the static state, stable at 0 K, based on first-principles structural optimization, and did not go as far as to discuss the dynamic properties at finite temperatures.

Conventional theoretical studies of solid-water interfaces, including their work, have been mainly based on first-principles calculations, in particular ab initio molecular dynamics (AIMD) for dynamic properties, which is computationally expensive although accurate. The high computational cost is a major problem, especially for models with a large number of atoms in the unit cell, such as catalysts with introduced defects, or when dynamic properties such as proton diffusion are studied using molecular dynamics calculations of 1 ns or longer. In recent years, neural network potentials (NNPs), which allow more efficient molecular dynamics calculations while maintaining the same level of accuracy as first-principles calcu
lations, have been used to study adsorption structures and proton transfer at the solid-water interface [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31].

The aim of this study is to investigate the structure of water molecules at the interface between water and defective oxide catalysts, by using neural network potentials molecular dynamics (NNPMD) as a first step towards understanding the effects of defects and solvents on catalytic activity of ORR. Examples of such effects include the influence of water layering on proton transfer and the poisoning of O\({}_{2}\) adsorption sites by water. Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\)[32], a defective ZrO\({}_{2}\) with V\({}_{\rm O}\) and N\({}_{\rm O}\) introduced, was used as a research object and ZrO\({}_{2}\) as a comparison object. For these materials, NNPs reproducing first-principles calculations were constructed based on active learning. Based on the NNPMD trajectories, the structures of water molecule adsorption on the surfaces of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\) are investigated and the possibility of the V\({}_{\rm O}\)  as catalytically active sites is discussed by comparing them.

## 2 Methods

Figure 1 shows the interface models used in this study. The atomic configuration of the 5-layer Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is one of the stable structures obtained by Monte Carlo sampling. The V\({}_{\rm O}\) and N\({}_{\rm O}\) defects are found both on the surface and inside. The atomic configuration of ZrO\({}_{2}\) is that of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) with the V\({}_{\rm O}\) and N\({}_{\rm O}\) replaced by oxygen atoms. These atomic configurations are listed in the Supporting Information. Both the unit cells of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\) have a hexagonal lattice with lattice constants \(a=b=9.542394\) A, \(c=25\) A. The interface model is a random arrangement of water molecules on a unit cell with an additional 15 A vacuum layer. The number of water molecules was set to 21, which corresponds to about three layers of water. The xy-coordinates of the water molecules were set so that the oxygen of the water molecules was above one of the atoms of the slab and the dipole moment of the water molecules was parallel to the x-axis. The 21 water molecules were randomly assigned to three different groups of seven water molecules each and, depending
Figure 1: The initial structures of molecular dynamics for interfaces between H\({}_{2}\)O and (a) Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) or (b) ZrO\({}_{2}\). Green, blue, red and white indicate Zr, N, O and H atoms respectively. The figures have been generated using the VESTA software package [33].
on the group, the z-component of the atomic coordinates was set to a value 3, 6, 9 A higher than the z-coordinate of the largest atom in the slab. The oxides prepared in this way, with water in a random arrangement, were used as the initial structure for AIMDs and NNPMDs. This procedure was used to create different initial structures that were not correlated with each other. When NNPMD is performed in a \(2\times 2\times 1\) supercell instead of a unit cell, the initial structure is not a \(2\times 2\times 1\) supercell of randomly placed water on a unit cell oxide, but randomly placed water on a \(2\times 2\times 1\) supercell oxide. Therefore, the initial structure of the slab always has translational symmetry, but that of the water molecules does not in most cases.

Active learning is used to improve the accuracy of NNPs in complex systems such as interfaces.[22, 25, 27, 30] The following active learning procedure was used to generate the NNP of these systems. First, we obtained the total energy, interatomic forces and stress tensor data corresponding to the atomic configuration from the AIMD for the slab model. An NNP reproducing the data was generated and used for the NNPMDs. We performed ab initio self-consistent calculations on the structures appearing in the trajectory and obtained the total energy, interatomic forces and stress tensor data corresponding to the atomic configuration. These new data were added to the existing data to generate a new NNP. The accuracy of the NNP was improved by active learning, which repeated this process. The active learning was stopped when the structural properties (e.g. the number of adsorptions) calculated from the NNPMD trajectories converged so that the difference between Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\) could be discussed. In this paper, the first NNP trained with AIMD only is referred to as the loop 0 NNP, while subsequent NNPs are referred to as loop 1, 2,... NNPs, depending on the number of times the ab initio calculation results on the structures of the NNPMD were added.

The calculation software used for each procedure and the detailed calculation conditions are as follows. First-principles calculations were performed using the VASP code,[34, 35] which is based on the plane-wave pseudopotential method. Projector augmented wave pseudopo
tentials of the type [36] and the exchange correlation functional of the generalised gradient approximation of the Perdew-Burke-Ernzerhof type [37] were used. The projection operator was evaluated in real space, including the non-spherical contribution to the gradient of the density of potential spheres. Only \(\Gamma\) points were sampled for the Brillouin zone and partial occupancies were calculated by Gaussian smearing with a width of 0.03 eV. The cutoff energy for the plane-wave-basis set is 520 eV. Self-inconsistent calculations were performed for up to 100 loops and stopped when the energy difference was less than 0.1 meV. The initial magnetic moment was set to 0.6 \(\mu_{\mathrm{B}}\) for all atoms to account for spin polarisation. The mass of the hydrogen atoms was replaced by the mass of the deuterium atoms, allowing a larger time step without breaking the time evolution. AIMD used the NVT ensemble whose temperature was controlled at 300 K by the Langevin thermostat with a friction coefficient of 10 ps\({}^{-1}\). Four initial structures were created using the method described above, each with the different initial velocities and 1 ps AIMDs were performed for each with a 1 fs time step. The initial data set thus contains 4000 samples.

The NNPs were generated using SIMPLE-NN (ver. 2) [38, 39], which is based on the symmetry function [40]. The parameters of the radial symmetry function G\({}_{2}\) are \(R_{\mathrm{cut}}=6.0,\eta=0.003214,0.035711,0.071421,0.124987,0.214264,0.3571 06,0.714213,1.428426,R_{\mathrm{s}}=0.0\). Those of the angular symmetry function G\({}_{3}\) are \(R_{\mathrm{cut}}=6.0,\eta=0.000357,0.028569,0.089277,\lambda=-1,1,\zeta=1,2,4\). A neural network is constructed with these symmetry functions as input and the total energy, the interatomic forces and the stress tensor as output. The nodes of the neural network are 30-30 and the activation function is tanh, and its parameters are optimised up to 1000 epochs using the Adam method with a learning rate of 0.001. The initial dataset of 4000 samples was split into 400 and 3600, the latter being used as test data for NNP accuracy testing. The former was further split 9:1 into training data and validation data; the training data is used to optimise the NNP weights and the validation data is used to calculate the loss function for each epoch to determine which epoch's NNP is adopted.

The NNPMDs were performed using LAMMPS [41, 42]. The NVT ensemble was used with
the temperature controlled to 300 K by the Langevin thermostat with the damping factor of 0.1 ps. The mass of the hydrogen atom was replaced by the mass of the deuterium atom. Four initial structures were created using the method described above, each with a different initial velocity and a 1 ns NNPMDs was performed for each with a 1 fs time step.

Self-inconsistent field calculations were performed in VASP for 400 structures every 10 ps in four 1 ns trajectories, and the ab initio and NNP energy errors were checked. Samples with an energy error greater than 1 meV/atom were added to the training and validation data, but not to the test data. Thus, as active learning progresses, the total of training and validation data increases from 400 samples, while the test data remains at 3600 samples. When the energy error in the unit cell is less than 10 meV/atom, in addition to the unit cell, a \(2\times 2\times 1\) supercell was used to test the energy error in the ab initio calculation and the NNP under almost similar conditions. In the case of the \(2\times 2\times 1\) supercell, the difference is that 40 structures every 10 ps in four 1 ns trajectories are validated and the ab initio calculation results are not used to train the NNP.

When the energy error was less than 10 meV/atom for both the unit cell and the \(2\times 2\times 1\) supercell, the structural properties used for the convergence decision condition of active learning were calculated. For the \(2\times 2\times 1\) supercell, 1000 new initial structures were generated and NNPMD was performed on each of them for 500 ps with a time step of 1 fs. The mean and standard deviation of the number of molecular and dissociative adsorptions were calculated using 1000 structures at 500 ps. The conditions for O adsorption on the Zr site and OH covalent bonds were set to 3.09651 A and 1.2 A respectively. When it was determined that these properties had sufficiently converged, active learning was stopped and the other properties - number density distribution for the surface perpendicular and parallel, adsorption distance d, adsorption angle \(\alpha\), molecular orientation \(\beta\), OH orientation \(\gamma\), number of hydrogen bonds and radial distribution function - were also calculated. Definitions of angles and hydrogen bonds will be explained later. In summary, if the RMSE of the unit cell is below 10 meV/atom, the RMSE of the \(2\times 2\times 1\) supercell is calculated; if it is also
below 10 meV/atom, the number of adsorptions is calculated; if it converges, active learning is stopped and the other properties are calculated.

## Results and discussion

Table 1 shows the root mean squared error (RMSE) for the Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\) NNPs against existing data sets. The energy error is less than 10 meV/atom in each loop. In most cases the relative difference between loops is less than 10 %. Although the NNP appears to converge sufficiently without going to loop 4, the NNP accuracy should not be judged solely on the basis of the RMSEs of the existing data. This is because, as discussed below, the RMSEs of the newly acquired data in the molecular dynamics trajectories performed with these NNPs are not necessarily 10 meV/atom.

Table 2 shows the errors in the ab initio self-consistent calculations for the structures included in the NNPMD trajectories of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). Ab initio self-consistent calculations were also performed for \(2\times 2\times 1\) supercells only in loops 1-4, where the RMSE in the unit cell is less than 10 meV/atom. In loop 1, the RMSEs and structural properties were not

\begin{table}
\begin{tabular}{r r r r r r r r r} \hline Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) & Energy & (meV/atom) & Force & (eV/Å) & & Stress & (kbar) & \\ loop & train & valid & test & train & valid & test & train & valid & test \\ \hline
0 & 1.119 & 1.184 & 1.729 & 0.088 & 0.118 & 0.157 & 1.752 & 1.987 & 2.022 \\
1 & 1.294 & 1.521 & 1.548 & 0.109 & 0.133 & 0.160 & 2.298 & 2.263 & 2.178 \\
2 & 1.998 & 1.695 & 2.347 & 0.113 & 0.128 & 0.148 & 2.023 & 1.980 & 2.103 \\
3 & 2.122 & 2.587 & 2.302 & 0.133 & 0.162 & 0.158 & 2.144 & 2.241 & 2.292 \\
4 & 1.530 & 1.866 & 1.818 & 0.133 & 0.141 & 0.159 & 2.488 & 2.558 & 2.795 \\ \hline ZrO\({}_{2}\) & Energy & (meV/atom) & Force & (eV/Å) & & & Stress & (kbar) & \\ loop & train & valid & test & train & valid & test & train & valid & test \\ \hline
0 & 0.826 & 0.819 & 1.061 & 0.083 & 0.103 & 0.130 & 1.822 & 1.781 & 1.938 \\
1 & 3.865 & 2.135 & 2.117 & 0.137 & 0.164 & 0.147 & 2.866 & 4.953 & 1.927 \\
2 & 3.565 & 2.088 & 2.177 & 0.151 & 0.163 & 0.141 & 3.232 & 4.531 & 2.647 \\
3 & 3.043 & 4.622 & 1.826 & 0.153 & 0.155 & 0.150 & 3.011 & 3.020 & 2.092 \\
4 & 3.267 & 1.834 & 1.295 & 0.151 & 0.163 & 0.138 & 2.919 & 2.508 & 1.783 \\ \hline \end{tabular}
\end{table}
Table 1: The root mean squared error (RMSE) of the existing data set for the total energy, the interatomic force and the stress tensor for the NNP of each loop in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\) for training, validation and test data.
calculated because the self-consistent calculations of \(2\times 2\times 1\) supercells did not converge for some structures. The RMSEs of loops 2, 3, and 4 were \(<\) 10 meV/atom in both cells, so structural properties were calculated. Note that the total energies in 1 ns NNPMDs reach equilibrium after about 0.1 ns and are conserved without significant increase or decrease thereafter.

Whether sufficient time and the number of trajectories were used to calculate the structural properties was verified by the number of adsorptions calculated from the loop 4 NNPMD trajectories of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\). Figure 2 (a) shows the time dependence of the number of adsorptions calculated from the 10, 20,..., 500 ps structure of 1000 trajectories. It can be seen that it converges to 0.1 at 400 ps. Figure 2 (b) shows the dependence of the number of adsorptions on the number of trajectories calculated from the 500 ps structure of 2, 3,..., 1000 trajectories. It can be seen that the number converges to 0.1 for 200 trajectories. Thus the number of adsorptions calculated from the 500 ps structure of 1000 trajectories can be trusted to the order of 0.1.

Table 3 shows the number of molecular and dissociative adsorptions (N) of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and

\begin{table}
\begin{tabular}{r r r r r r r r} \hline Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) & unit cell & & & \multicolumn{4}{c}{\(2\times 2\times 1\) supercell} \\ loop & run = 1 & 2 & 3 & 4 & run = 1 & 2 & 3 & 4 \\ \hline
0 & 10.572 & 17.765 & 11.313 & 10.633 & & & & \\
1 & 4.833 & 7.105 & 5.323 & 5.113 & & & & \\
2 & 2.317 & 2.140 & 8.145 & 2.143 & 2.279 & 2.353 & 5.398 & 2.299 \\
3 & 1.626 & 1.541 & 1.603 & 1.442 & 1.931 & 6.112 & 6.553 & 6.647 \\
4 & 1.202 & 1.178 & 1.392 & 1.384 & 1.310 & 1.199 & 1.428 & 1.345 \\ \hline ZrO\({}_{2}\) & unit cell & & & & \multicolumn{4}{c}{\(2\times 2\times 1\) supercell} \\ loop & run = 1 & 2 & 3 & 4 & run = 1 & 2 & 3 & 4 \\ \hline
0 & 28.576 & 28.416 & 28.955 & 25.323 & & & & \\
1 & 1.900 & 1.939 & 1.899 & 2.100 & & & & \\
2 & 1.628 & 2.553 & 1.514 & 1.573 & 6.632 & 6.231 & 5.308 & 5.153 \\
3 & 1.469 & 1.094 & 1.088 & 1.152 & 4.936 & 5.024 & 5.159 & 5.116 \\
4 & 2.644 & 1.160 & 1.625 & 2.123 & 3.021 & 3.902 & 3.884 & 3.543 \\ \end{tabular}
\end{table}
Table 2: RMSE in the total energy of the newly acquired data, which are included in the NNPMD trajectories of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). For loops 0 and 1, RMSEs for \(2\times 2\times 1\) supercells were not calculated due to large RMSEs for unit cell and the non-convergence for self-consistent calculation, respectively.
\begin{table}
\begin{tabular}{r r r r r r r} \hline \multicolumn{6}{l}{Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\)} \\ loop & \multicolumn{1}{c}{N\({}_{\rm mol}\)} & \multicolumn{1}{c}{N\({}_{\rm dis}\)} & \multicolumn{1}{c}{M\({}_{\rm mol}\)} & \multicolumn{1}{c}{M\({}_{\rm dis}\)} & \multicolumn{1}{c}{C\({}_{\rm mol}\)} & \multicolumn{1}{c}{C\({}_{\rm dis}\)} \\ \hline
2 & 9.6 (1.2) & 5.8 (1.1) & 2.4 (0.3) & 1.5 (0.3) & 34 (4) & 21 (4) \\
3 & 10.2 (1.2) & 5.5 (1.0) & 2.5 (0.3) & 1.4 (0.2) & 36 (4) & 20 (3) \\
4 & 9.5 (1.1) & 5.5 (1.0) & 2.4 (0.3) & 1.4 (0.3) & 34 (4) & 20 (4) \\ \hline \multicolumn{6}{l}{ZrO\({}_{2}\)} \\ loop & \multicolumn{1}{c}{N\({}_{\rm mol}\)} & \multicolumn{1}{c}{N\({}_{\rm dis}\)} & \multicolumn{1}{c}{M\({}_{\rm mol}\)} & \multicolumn{1}{c}{M\({}_{\rm dis}\)} & \multicolumn{1}{c}{C\({}_{\rm mol}\)} & \multicolumn{1}{c}{C\({}_{\rm dis}\)} \\ \hline
2 & 4.1 (1.4) & 9.3 (1.3) & 1.0 (0.4) & 2.3 (0.3) & 15 (5) & 33 (5) \\
3 & 6.2 (1.4) & 9.1 (1.2) & 1.6 (0.4) & 2.3 (0.3) & 22 (5) & 33 (4) \\
4 & 5.6 (1.4) & 9.6 (1.2) & 1.4 (0.4) & 2.4 (0.3) & 20 (5) & 34 (4) \\ \hline \end{tabular}
\end{table}
Table 3: Mean values of the number of molecular (mol) and dissociative (dis) adsorptions in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). The values in brackets are standard deviations. For both materials, these were calculated from loop 2, 3 and 4 NNPMD trajectories. N = the number of adsorptions in a \(2\times 2\times 1\) supercell. M = N/4 = the number of adsorptions per unit cell. C = coverage (%) = N/28(= the number of Zr sites in \(2\times 2\times 1\) supercell surface)*100.

Figure 2: Dependence of the number of adsorptions on (a) timestep and (b) the number of trajectories. Blue and red represent molecular (mol) and dissociative (dis) adsorption, respectively. These were calculated from loop 4 NNPMD trajectories of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\).
ZrO\({}_{2}\) calculated by each loop NNP. For reference, the number of adsorptions per unit cell (M), converted from the number of adsorptions per \(2\times 2\times 1\) supercell, and the coverage (C) are also given. For the number of adsorptions N, the mean difference between loops 3 and 4 is up to 0.7 and the standard deviation between the materials is up to 2.5. Both are small compared to the difference between materials in loop 4, 2.9, so we consider that the convergence is not a problem for comparing these materials.

The relationship between the number of molecular and dissociative adsorptions is reversed between Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). Because Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) has less oxygen for the V\({}_{\rm O}\), the dissociated excess H\({}^{+}\) is less likely to be adsorbed on oxygen and more likely to re-bind to dissociatively adsorbed OH\({}^{-}\), which is the cause. This idea is also supported by the fact that the total number of adsorptions is almost the same for Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). In addition, Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is charge neutral due to the presence of the V\({}_{\rm O}\) and N\({}_{\rm O}\) defects, so it is unlikely that the V\({}_{\rm O}\) is responsible for attracting the O atom. According to a study by Muhammady et al. [15], where one V\({}_{\rm O}\) and two N\({}_{\rm O}\) defects were introduced into the tetragonal-ZrO\({}_{2}\) (101) surface, the V\({}_{\rm O}\) is not more easily adsorbed by the O\({}_{2}\) molecule than the Zr atom. Therefore, it is consistent with the present results to consider that the V\({}_{\rm O}\) does not promote dissociation by accepting the O of water molecules. The dissociative adsorption becomes molecular adsorption, which reduces the probability that intermediates and the H\({}^{+}\) that binds to them in the ORR will have an additional reaction with nearby adsorbates, i.e. the ORR is more likely to proceed in the presence of oxygen defects. Also, while the total number of adsorptions is almost the same and the number of vacant Zr sites is almost the same, Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) has two V\({}_{\rm O}\)s per unit cell, which could also be an active site for the ORR. In other words, Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) has more active sites and is more likely to promote the ORR than ZrO\({}_{2}\).

Figure 3 shows the surface perpendicular number density distributions of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). The densities of total, oxygen and hydrogen atoms in water at the interface and hydrogen atoms in bulk water are shown respectively. The surface perpendicular number density distribution shows that the oxygen peak in the adsorption layer (16-18 A) is one in
ZrO\({}_{2}\), whereas it is divided into two in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\), which is due to the fact that the height of the Zr sites is almost the same in all ZrO\({}_{2}\) and is divided into two types in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\). Differences in the surface-derived adsorption structures will affect their function as catalysts for the ORR, as well as the diffusion of molecules and protons.

Away from the surface, the number density oscillates and is not constant as in the bulk. In other words, the motion is limited compared to the bulk. According to a previous study of the Cu/H\({}_{2}\)O interface [16], water molecules move like the bulk when the thickness of their layer is more than 30 A. In the present study, the thickness of the water molecule layer is less than 15 A, which is consistent with the results that it is not constant like the bulk.

Figure 4 shows the surface parallel number density distributions of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\). Only the densities of the oxygen atoms in the H\({}_{2}\)O or OH\({}^{-}\) adsorbed on the Zr sites are shown. The distribution of ZrO\({}_{2}\) is uniform and translational symmetry is observed. As mentioned above, this is not due to the initial structure, as the translational symmetry of water is not guaranteed. ZrO\({}_{2}\) was assumed to have a uniform adsorption distribution due to the equivalent Zr sites on its surface. However, if the adsorption distribution does not

Figure 3: Surface perpendicular number density distributions (Å\({}^{-1}\)) of (a) Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and (b) ZrO\({}_{2}\). The blue, red and yellow solid lines represent total, oxygen and hydrogen atoms in water at the interface. The red and yellow dashed lines represent oxygen and hydrogen atoms in bulk water. These densities for both materials were calculated from loop 4 NNPMD trajectories.
Figure 4: Surface parallel number density distributions (Å\({}^{-3}\)) of O atoms in (a) molecular and (b) dissociative adsorption for Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and (c) molecular and (d) dissociative adsorption for ZrO\({}_{2}\). For visibility, the sign of the molecular adsorption number density is reversed. Two types of oxygen vacancies in a Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) unit cell are represented by v1 and v2.
change much in 500 ps MD, the adsorption distribution calculated from only one trajectory may not be uniform. In the present study, the adsorption distribution was calculated from 1000 independent trajectories, so that a uniform distribution was obtained, regardless of whether the adsorption distribution was likely to change in 500 ps MD. The distribution of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is non-uniform, but translational symmetry is established. This is also a result of obtaining a sufficiently large number of trajectories, rather than the initial structure, and is due to the low computational cost of NNPMD, which allows so many MDs to be performed.

Particles were found to adsorb on the Zr, but not on the V\({}_{\rm O}\). Zr atoms around V\({}_{\rm O}\) 2 are more likely to adsorb particles than those around V\({}_{\rm O}\) 1 These results suggest that the V\({}_{\rm O}\) act as active sites for the ORR and that differences in adsorption around them affect their function. In particular, we suggest that the V\({}_{\rm O}\) 1 is more likely to adsorb O\({}_{2}\) molecules as there is less adsorption on the surrounding Zr atoms.

For other properties, (adsorption distance d, adsorption angle \(\alpha\), molecular orientation \(\beta\), OH orientation \(\gamma\), and number of hydrogen bonds N\({}_{\rm HB}\)) their definiti

Figure 5: (a) Definition of adsorption angle \(\alpha\), molecular orientation \(\beta\) and OH orientation \(\gamma\). (b) Definition of hydrogen bonding. Oxygens donating and accepting hydrogen are denoted O\({}_{d}\) and O\({}_{a}\). Green, red and white indicate Zr, O and H atoms respectively. The figures have been generated using the VESTA software package.
results are presented. Figure 5 (a) shows the definition of the adsorption angle \(\alpha\), molecular orientation \(\beta\) and OH orientation \(\gamma\). The adsorption angle \(\alpha\) is the angle between the position vector of O with respect to Zr and the surface normal. The adsorption distance d is the absolute value of the vector. The water molecule orientation \(\beta\) is the angle between the dipole vector of the water molecule and the surface normal. The OH orientation \(\gamma\) is the angle between the position vector of H with respect to O and the surface normal. Figure 5 (b) shows the definition of a hydrogen bond. A hydrogen bond is defined when the angle H-O\({}_{d}\)-O\({}_{a}\) is less than 30\({}^{\circ}\), the distance O\({}_{d}\)-H is less than 1.2 A and the distance O\({}_{d}\)-O\({}_{a}\) is less than 3.5 A, where O\({}_{d}\) and O\({}_{a}\) are the oxygen donating and accepting hydrogen, respectively [43].

Table 4 shows the adsorption distance d, adsorption angle \(\alpha\), molecular orientation \(\beta\), OH orientation \(\gamma\), and the number of hydrogen bonds N\({}_{\rm HB}\) calculated based on the above definition. The adsorption distance d of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is 0.1 A shorter in both molecular and dissociation. The increase or decrease of the adsorption angle \(\alpha\) is opposite for molecular and dissociation, with the adsorption angles becoming unequal, whereas for ZrO\({}_{2}\) they were almost equal at 18\({}^{\circ}\). The molecular orientation \(\beta\) of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) decreases, i.e. it turns upwards. The increase and decrease of the OH orientation \(\gamma\) is opposite in molecular and dissociation, going from unequal in ZrO\({}_{2}\) to aligned at almost 62\({}^{\circ}\). The number of hydrogen bonds in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is reduced by one. We consider this to be the effect of a distortion in the distribution of water molecules due to the defect.

Figure 6 shows the Radial Distribution Function (RDF). The oxygens in the slab and in the water is distinguished as O\({}_{\rm s}\) and O\({}_{\rm w}\). Compared to ZrO\({}_{2}\), the RDF of Zr-O\({}_{\rm w}\) in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) has a lower peak at 2.1 A and is shifted to the right. This means that the average adsorption

\begin{table}
\begin{tabular}{l c c c c c c c} \hline  & d\({}_{\rm mol}\) & d\({}_{\rm dis}\) & \(\alpha_{\rm mol}\) & \(\alpha_{\rm dis}\) & \(\beta_{\rm mol}\) & \(\gamma_{\rm mol}\) & \(\gamma_{\rm dis}\) & N\({}_{\rm HB}\) \\ \hline Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) & 2.29 & 2.11 & 12.81 & 20.23 & 35.03 & 62.14 & 62.38 & 67.56 \\ ZrO\({}_{2}\) & 2.39 & 2.20 & 18.28 & 17.74 & 51.68 & 68.40 & 53.35 & 68.70 \\ \hline \end{tabular}
\end{table}
Table 4: Adsorption distance d (Å), adsorption angle \(\alpha\), molecular orientation \(\beta\), OH orientation \(\gamma\) and number of hydrogen bonds N\({}_{\rm HB}\) for Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\), for molecular (mol) and dissociative (dis) adsorption.
distance increases. Both molecular and dissociative adsorption distances are shorter for Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\), but the increase in the proportion of molecular adsorption at longer adsorption distances is thought to be the cause. The lower peak at 1 A in the O\({}_{\rm s}\)-H RDF of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is thought to be due to a decrease in H adsorbed on O\({}_{\rm s}\) as dissociative adsorption decreases.

In addition to structural properties such as those in this study, other studies of solid/water interface systems using NNPMD include the free energy surface for proton transfer [18, 22, 29, 44], that for water dissociation [30, 31], solvation dynamics [27], and vibrational spectra of OH bonds [19]. Systems with intermediates in the ORR could also be added to the data to learn NNPs, sampled by NNPMD and calculate the free energies on the reaction pathways, so that the influence of defects and the proportion of molecular and dissociative adsorption on catalytic activity could also be discussed. There is also much to be explored in studies based on conventional first-principles calculations. In contrast to the present study, there are previous studies that have discussed the electric double layer (EDL) by calculating the electrostatic potential and density when the system is about a system with ions in water [45, 46, 47, 48, 49]. For example, Ando et al. performed AIMD at the interface between Pt and water containing Na\({}^{+}\) ions and calculated the electrostatic potential due to the EDL based on the difference in

Figure 6: Radial distribution function \(g(r)\) of (a) Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and (b) ZrO\({}_{2}\). O\({}_{\rm s}\) and O\({}_{\rm w}\) represent oxygen in the slab and in the water molecule. Blue, green, yellow and red represent Zr-O\({}_{\rm w}\), Zr-H, O\({}_{\rm s}\)-O\({}_{\rm w}\) and O\({}_{\rm s}\)-H, respectively.
the potentials between the case with and without Na\({}^{+}\) ions [45]. For Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) and ZrO\({}_{2}\), the electrostatic potential due to the EDL could also be evaluated in NNPMD by training NNP with additional data for the case of the solvent containing ions. Alternatively, as in the work of Otani et al. [50], the electrostatic potential of H\({}_{3}\)O\({}^{+}\) in oxygen can be discussed by creating a system with excess H\({}^{+}\) and electrons in water and performing an AIMD based on the effective screening medium method. Other studies [51, 52, 53, 54], have performed AIMD and structural optimisation to discuss band offsets, which are important values in interface systems, taking into account electron affinity, charge neutrality level and pinning factor, and have discussed the band offset, which is an important value in interface systems. This multifaceted study will allow a more detailed investigation of the effects of defects at the oxide/water interface.

## 4 Conclusion

In this study, NNPs were constructed by active learning for interfaces between the water and defective Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) or pristine ZrO\({}_{2}\). Based on 1000 NNPMD trajectories, the structures of water molecule adsorption on these interfaces were investigated and it was found that Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) has a bilayer structure, which is different from a monolayer structure of ZrO\({}_{2}\). On the Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) surface, water molecules do not adsorb on the V\({}_{\rm O}\)  but only on some of the surrounding Zr atoms, suggesting that the V\({}_{\rm O}\) is expected to act as an active site for the ORR. While the proportion of dissociative adsorption in Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) is less than in ZrO\({}_{2}\) due to the less oxygen on which the H\({}^{+}\) can adsorb, the total number of molecular and dissociative adsorption is the same for both materials. Therefore, Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) may promote the ORR more than ZrO\({}_{2}\) because the V\({}_{\rm O}\)s act as additional active sites. In the future, the influence of the solvent on the ORR activity of Zr\({}_{7}\)O\({}_{8}\)N\({}_{4}\) could be discussed by updating the NNP with new data for systems with intermediate oxygen and hydrogen molecules in addition to the solvent, water molecules.
## Acknowledgement

This research was supported by the New Energy and Industrial Technology Development Organization (NEDO) project, MEXT as "Program for Promoting Researches on the Supercomputer Fugaku" (Fugaku battery and Fuel Cell Project) (Grant No. JPMXP1020200301, Project No.: hp220177, hp210173, hp200131), and the Japan Society for the Promotion of Science (JSPS) as Grants in Aid for Scientific Research on Innovative Area "Hydrogenomics" (Grant No. 18H05519). The calculations were conducted at the ISSP Supercomputers Center, The University of Tokyo.
