Medium-range order in amorphous ices revealed by persistent homology

For MD simulations we used the TIP4P/2005 model [23] in the LAMMPS software [24]. We first prepared a liquid water structure as an isothermal-isobaric NPT ensemble of 2880 water molecules in a cubic box with periodic boundary conditions at $\newcommand{\K}{~\mathrm{K}} T=280\K$ and $\newcommand{\mpa}{~\mathrm{MPa}} P=0.1\mpa$. The isobaric cooling cycle on the liquid water to obtain hyperquenched glassy water, i.e. LDA ice, ran for 200 ns at a cooling rate of $\newcommand{\K}{~\mathrm{K}} 1\K$ ns⁻¹. Then it is followed by the isothermal compression of LDA ice to get HDA ice at 1.3 GPa, running for 18 ns with a compression rate of 0.1 GPa ns⁻¹. The time step for the integration algorithm was set to be 1 fs and the damping parameters for the thermostat and barostat are 0.5 and 5 ps, respectively. We cut off the Coulombic and Lennard-Jones interactions at distance 10 $\mathring{\rm A}$ and 10.3092 $\mathring{\rm A}$, respectively [18].

The structure factor $S(\mathbf{Q})$ measures the scattered intensity of an incident radiation for the scattering vector $\mathbf{Q}$ and is related to the radial distribution function (or pair correlation function) via the Fourier transform. For a perfect crystal, it exhibits infinitely many sharp peaks reflecting the periodicity of the system; for a liquid-like arrangement, a certain degree of fluctuations appears in the distribution function, yet the sharpness is weak. The numerical estimation of the structure factor of a multiatomic system of N atoms $\alpha_j\,(1\leqslant j\leqslant N)$ at positions $\mathbf{x}_j$ is as follows [25]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{SF} S(Q) = \frac{1}{N}\left|\sum_{j=1}^{N} f_j(Q)\, {\rm e}^{-{\rm i}\mathbf{Q}\mathbf{x}_j}\right|^2, \nonumber \end{align} \tag{ 1 }$$

where $Q=|\mathbf{Q}|=\frac{2\pi}{\lambda}$ is the wave number for the wavelength $\lambda$ of the plane wave and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_j(Q) = \sum_{k=1}^4 a_{j,k}\, {\rm e}^{-b_{j,k}\left(\frac{Q}{4\pi}\right)^2} + c_j \nonumber \end{align} \tag{ 2 }$$

is the (approximated) atomic form (or scattering) factor of each atom $\alpha_j$ with constants $a_{j}, b_{j}$ and c_j reflecting the electron density of atom $\alpha_j$ under the assumption that the electron density is spherically symmetric [26].

For structural analysis, we used persistent homology (PH) computations from an emerging research field, topological data analysis [27, 28]. PH has been actively utilized in a wide range of domains, including materials science and neuroscience [29–34], as an effective detector of topological features, i.e. holes embedded in a system. For each dimension n, we define the nth persistent homology group PH_n, which extracts qualitative features of n-dimensional holes, e.g. connected components for n  =  0, rings for n  =  1 and cavities for n  =  2, that persist over multiple scales. Since we aim to detect MRO ring structures in water that are responsible for the FSDP in the structure factor, our focus is on the case n  =  1.

As one of the most efficient representations of PH computations, a persistence diagram (PD) is used to visualize the distributions of size and shape of hole structures. PD can precisely reveal hierarchical structural units in various length scales from short- to long-range order in the system. In this sense, PD has great merit compared to other conventional methods for structural analyses, such as radial distribution functions (RDF), distributions of bond or dihedral angles and ring statistics [35]. The detailed PH computing mechanism for our study using HomCloud software [36, 37] follows in section 3.2 along with a schematic illustration.

We performed out-of-equilibrium MD simulations on the system of 2880 water molecules to obtain LDA and HDA ices, separated by two steps (figure 1(a)) as described in section 2.1. Figure 1(b) shows the change in density at each step, which is in good agreement with the experimental densities $\rho_{{\rm LDA}}\approx$ 0.95 $\mathrm{g~cm^{-3}}$ and $\rho_{{\rm HDA}}\approx$ 1.31 $\mathrm{g~cm^{-3}}$ at $\newcommand{\K}{~\mathrm{K}} T=77\K$ and $P=1\,$ GPa [5]. Note that we took LDA and HDA ice structures at P  =  0.1 and $\newcommand{\mpa}{~\mathrm{MPa}} 1300\mpa$, respectively, of which densities agree the best with the experimental results.

Next, the structure factor $S(Q)$ of each glassy water system is computed and compared with experimental results [22] to confirm the reliability of our MD-generated amorphous ice systems. (See section 2.2 for the computational details.) For an amorphous configuration, the FSDP at small $Q\neq 0$ in $S(Q)$, if it exists, is closely related to MRO within the system [14, 38]. Figure 2 compares the experimental structure factor (top) with the computed $S(Q)$ (bottom) of our water structures generated by MD simulations. The FSDP position Q^* in each case is well-matched to the experimental result at $Q^*_{{\rm LDA}}\approx 1.75$ and $Q^*_{{\rm HDA}}\approx 2.25$ within the error bound of $\sim 0.05$ $\mathring{\rm A}$⁻¹. Note that the position of the peak varies slightly from the experimental result due to the difference in external conditions of temperature and pressure in our simulation setup.

Zoom In Zoom Out Reset image size

The PH computing mechanism is as follows. Suppose we have an atomic configuration of N atoms $\alpha_j$ centered at $\mathbf{x}_j\in\mathbb{R}^3\,(1\leqslant j\leqslant N)$ and a set $\{r_j\in\mathbb{R}: 1\leqslant j\leqslant N\}$ of input radii. To compute PH₁, we place a spherical ball of radius r_j centered at $\mathbf{x}_j$, and enlarge the ball radius to $\sqrt{r_j^2+\varepsilon}$ with increasing $\varepsilon$, the scale parameter in the computation. When two balls meet, we add a line segment to connect the corresponding two atoms. As $\varepsilon$ increases, more line segments appear and eventually create a ring by connecting multiple atoms end to end. When a ring $\newcommand{\R}{\mathfrak{R}} \R$ is born, $\varepsilon$ is recorded as the birth scale $\newcommand{\R}{\mathfrak{R}} b_\R$ of the ring. In contrast, we define the death of ring $\newcommand{\R}{\mathfrak{R}} \R$ when $\newcommand{\R}{\mathfrak{R}} \R$ is completely triangulated and every three balls has nonempty intersection. The $\varepsilon$ at the death point is recorded as the death scale $\newcommand{\R}{\mathfrak{R}} d_\R$ of the ring. The very last triangulated ring is called the death simplex, and its maximum side length represents the size of the ring. [31, 39]

Figure 3(a) shows a schematic illustration for the PH computing process with the resulting PD of a 2D version of a water configuration. In our study, we use the following input radii of atoms: taking the first peak positions $r_{{\rm OH}}=0.95$ $\mathring{\rm A}$ and $r_{{\rm HH}}=1.55$ $\mathring{\rm A}$ of the RDF of OH and HH pairs into consideration (figure S1 in supplementary data (stacks.iop.org/JPhysCM/31/455403/mmedia)), we solve the equations $r_{{\rm HH}} = 2\,r_{{\rm H}}$ and $r_{{\rm OH}} = r_{{\rm O}} + r_{{\rm H}}$ to take $r_{{\rm H}}=0.775$ $\mathring{\rm A}$ and $r_{{\rm O}}=0.175$ $\mathring{\rm A}$ as the input radii of H and O atoms, respectively. As $\varepsilon$ increases, growing balls intersect, thereby line segments connect atoms, so that multiple rings $\newcommand{\R}{\mathfrak{R}} \R_j\, (\,j=1,\cdots,8)$ are created. Note that red triangles $\newcommand{\R}{\mathfrak{R}} \R_j\, (\,j\neq 4)$ are small enough to disappear shortly after their birth, while the blue ring $\newcommand{\R}{\mathfrak{R}} \R_4$ persists longest. In the PD, the longer a ring persists, the farther away it is placed from the diagonal. The largest ring $\newcommand{\R}{\mathfrak{R}} \R_4$ has $\newcommand{\R}{\mathfrak{R}} \R_j$ for $j=5,6,7$ as its subrings and the unfilled blue triangle at $\varepsilon=\varepsilon_5$ is the death simplex of $\newcommand{\R}{\mathfrak{R}} \R_4$.

Zoom In Zoom Out Reset image size

Next, the PD₁ for our MD-generated amorphous ices are aligned in figure 3(b). Note that the birth–death pairs of rings in PDs are plotted in a squared scale of $\mathring{\rm A}$². PDs of both LDA and HDA ices show a division of points into distinct areas, characterized by lines or band-like regions, where a large number of points (each representing a ring) are concentrated. We denote the separated characteristic regions by L_j and $H_j\, (\,j=1,2,3,4)$ for LDA and HDA, respectively, as depicted in the figure. Note that rings in vertical island regions L₁ and H₁ are composed of more than four member atoms ($m\geqslant 4$) and their triangular or quadrilateral subrings are recorded in other regions close to the diagonal. Since rings in L₄ and H₄ are small triangles whose length scales contribute little to the FSDP in the structure factor, our focus is rather on rings in other characteristic regions ($j=1,2,3$) that would suitably correspond to the MRO in glassy water and are likely to be associated to the sharp peaks in $S(Q)$ in figure 2.

The geometry of MRO ring structures of LDA and HDA ices substantially differs from each other in several aspects. First, stacked bar charts in figure 4(a) compare the distributions of the number of ring members ($m\leqslant 30$) in the vertical island regions L₁ and H₁ in PDs. We denote a ring with m members by m-ring. While L₁ has a sharp distribution of the number of members, mostly concentrated on 7–10 members, H₁ has a broader distribution, as many rings of H₁ involve a larger number of atoms as their components. The median number of ring members is 9 and 12 in L₁ and H₁, respectively. As mentioned earlier, rings in other characteristic regions for $j\neq1$ (close to the diagonal) are simple triangles or 4-rings. We observe in figure 4(a) that HH atom pairs, as opposed to OH or OO (none) pairs, dominantly determine death scales of MRO rings both in LDA and HDA ices. The inset figures show six types of representative m-rings in L₁ ($\mbox{\textcircled{\tiny 1}}$–$\mbox{\textcircled{\tiny 3}}$) and in H₁ ($\mbox{\textcircled{\tiny 4}}$–$\mbox{\textcircled{\tiny 6}}$) extracted from the marked spots in the PDs in figure 3(b) with red lines connecting the atom pairs responsible for the death of each ring.

Zoom In Zoom Out Reset image size

In addition to the number of ring members, MRO rings in LDA and HDA ices are differentiated by their size and shape. For this purpose, we consider the following length of ring $\newcommand{\R}{\mathfrak{R}} \R$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\R}{\mathfrak{R}} \displaystyle \label{ell_d} \delta(\R) = \sqrt{r_{\alpha}^2 + d_\R\vphantom{r_{\alpha'}^2}}+\sqrt{r_{\alpha'}^2 + d_\R}, \nonumber \end{align} \tag{ 3 }$$

which is the distance between the centers of atoms $\alpha$ and $\alpha'$, the pair of which determines the death scale $\newcommand{\R}{\mathfrak{R}} d_\R$ of ring $\newcommand{\R}{\mathfrak{R}} \R$ (i.e. depicted as red lines in figure 3(a)) [39]. Note that $\newcommand{\R}{\mathfrak{R}} \delta(\R)$ represents the size of $\newcommand{\R}{\mathfrak{R}} \R$, as $\newcommand{\R}{\mathfrak{R}} \delta(\R)$ is the maximum edge length within the triangulated ring, and is invariant under the choice of input radii, since it is determined by the given configuration itself. We plot red line graphs in figure 4(b) to show the mean value of $\newcommand{\R}{\mathfrak{R}} \delta(\R)$ of m-rings in L₁ and H₁ regions with respect to the number of ring members m, tending to increase as m becomes larger. The average sizes $\langle\delta\rangle$ of the rings in L₁ and H₁ contributory to the FSDPs are estimated to be notably different as $3.58\pm 0.23$ $\mathring{\rm A}$ and $2.99\pm 0.09$ $\mathring{\rm A}$, respectively.

It might be counterintuitive that rings with more members in H₁ are smaller than rings in L₁ with fewer members. This is also confirmed in PDs where the death scale (i.e. the size of the ring) exhibits a range of smaller length scales in HDA than in LDA ice. However, as water molecules become closer to each other in HDA, mainly by the pressure effect, the movement naturally leads to smaller voids in twisted ring bodies with more members, which is consistent with higher density of HDA. We define the following formula to study the degree of flatness (or planarity) of a detected MRO ring $\newcommand{\R}{\mathfrak{R}} \R$ of m members:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\R}{\mathfrak{R}} \displaystyle F(\R) = \frac{1}{m} \sum_{j=1}^m {\rm dist}(\mathbf{x}_j, P), \nonumber \end{align} \tag{ 4 }$$

where $\mathbf{x}_j$ is the position of atom $\alpha_j$ and P is the best fit plane to the ring geometry that minimizes the sum of the distances of member atoms from the plane. Note that $F\geqslant0$ with F  =  0 indicating that all the ring members are coplanar, hence a flat ring. MRO rings in L₁ display flatter shapes with $\langle F\rangle=0.42\pm 0.12$ and $F_{\max}=0.52$, while those in H₁ have $\langle F\rangle=0.56\pm 0.19$ and $F_{\max}=0.87$. The mean flatness $\langle F\rangle$ of rings in L₁ (left) and H₁ (right) are shown with blue dotted lines in figure 4(b). Two 9-rings, each found in L₁ and H₁ regions, are shown in figure 4(c). This geometry is compatible to the picture where water molecules in the second coordination shell in the basic tetrahedral network penetrate into the first shell by the increased pressure, resulting in the deformed tetrahedral network in HDA ice.

To reproduce the FSDP in $S(Q)$ from PD₁ results, we use the following equation [31]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\R}{\mathfrak{R}} \displaystyle S^*(Q) = \frac{1}{|\Pi|} \sum_{(b_\R,d_\R)\in \Pi} \exp{\left[-\frac{\left(Q-\frac{2\pi}{\delta(\R)}\right)^2}{2\sigma^2}\right]}, \nonumber \end{align} \tag{ 5 }$$

where $\Pi$ is a characteristic region L_j or $H_j\, (\,j=1,2,3)$ on PDs with $|\Pi|$ the number of rings in the region, $\newcommand{\R}{\mathfrak{R}} (b_\R,d_\R)$ is the birth–death pair of ring $\newcommand{\R}{\mathfrak{R}} \R$ in a squared scale of $\mathring{\rm A}$² and $\sigma=0.05$ applies a smoothing effect on the curve. Figure 5 shows the reproduced FSDPs of LDA (top) and HDA (bottom) by S^*(Q). We first consider each individual characteristic region $\Pi_j\,(\,j=1,2,3)$ on PD to plot the corresponding S^*(Q) and then their combination $\Pi_1\cup \Pi_2\cup \Pi_3$ for each LDA and HDA ice. Note that rings in L₁ and L₂ regions contribute much to reproducing the peak at $Q^*_{{\rm LDA}}$ apart from the rings in L₃ region, which are of smaller scale than those in other regions. This signifies that MRO rings in L₁ and L₂ are the main contributors to the FSDP in $S(Q)$. Also in the case of HDA, the reproduced peak from the combined characteristic region $H_1\cup H_2\cup H_3$ matches fairly well with the designated spot.

Zoom In Zoom Out Reset image size

In addition to the peak position, the width of the FSDP is also examined. We calculate the ratio $\Omega=w_{Q^*_\mathrm{HDA}}/w_{Q^*_\mathrm{LDA}}$, where $w_{Q^*_\mathrm{LDA}}$ and $w_{Q^*_\mathrm{HDA}}$ are the full widths at half maximum of the FSDPs of the amorphous ices. The ratio $\Omega_{S^*(Q)}$ of the reproduced peaks is estimated to be 1.63, which excellently agrees with $\Omega_{{\rm x-ray}} \approx$ 1.59 of the experimental data (top panel in figure 2) and $\Omega_{\mathrm{MD}} \approx$ 1.60 of our MD-generated ice configurations (bottom panel in figure 2). All these values consistently reflect that the FSDP width of HDA ice is larger than that of LDA ice. Overall, PH-extracted MRO rings successfully reproduce the quantities associated with the FSDPs, including the peak positions and widths, in the structure factor, which convincingly supports that these rings construct the hidden MRO in glassy water, leading to the appearance of the sharp peaks.

In a multiatomic system as water, the larger the electron density of an atom has, the more effect it adds to the diffraction pattern. This fact has been one of the reasons for higher attention to the oxygen arrangement only in the previous studies of water structures [17, 40]. As a comparison, we further carried out the FSDP reconstruction process based on the oxygen arrangement only (O-system hereafter) and compared the results to the case of the entire H₂O system (figure S2 in supplementary data). While the trend of the FSDP of HDA appearing at larger Q than that of LDA is preserved, O-system does not accurately reproduce the quantities associated with the FSDP (e.g. peak position and width). For the O-systems in the bottom panel in figure S2, the peak positions deviate by approximately 0.4–0.6 $\mathring{\rm A}$⁻¹, and also the notably larger peak width of HDA than that of LDA is not reproduced. Thus, it is critical to keep both oxygen and hydrogen atoms into consideration for an accurate PH analysis of water in ensemble. This is a hitherto unknown result, which would have been difficult to obtain in earlier approaches, as previously known structural information was mostly based on oxygen atoms in water, overlooking hydrogen atoms as a mere decoration.