Hierarchical structures of amorphous solids characterized by persistent homology

Edited by Giorgio Parisi, University of Rome, Rome, Italy, and approved April 22, 2016 (received for review October 22, 2015)
June 13, 2016
113 (26) 7035-7040

Significance

Persistent homology is an emerging mathematical concept for characterizing shapes of data. In particular, it provides a tool called the persistence diagram that extracts multiscale topological features such as rings and cavities embedded in atomic configurations. This article presents a unified method using persistence diagrams for studying the geometry of atomic configurations in amorphous solids. The method highlights hierarchical structures that conventional techniques could not have treated appropriately.

Abstract

This article proposes a topological method that extracts hierarchical structures of various amorphous solids. The method is based on the persistence diagram (PD), a mathematical tool for capturing shapes of multiscale data. The input to the PDs is given by an atomic configuration and the output is expressed as 2D histograms. Then, specific distributions such as curves and islands in the PDs identify meaningful shape characteristics of the atomic configuration. Although the method can be applied to a wide variety of disordered systems, it is applied here to silica glass, the Lennard-Jones system, and Cu-Zr metallic glass as standard examples of continuous random network and random packing structures. In silica glass, the method classified the atomic rings as short-range and medium-range orders and unveiled hierarchical ring structures among them. These detailed geometric characterizations clarified a real space origin of the first sharp diffraction peak and also indicated that PDs contain information on elastic response. Even in the Lennard-Jones system and Cu-Zr metallic glass, the hierarchical structures in the atomic configurations were derived in a similar way using PDs, although the glass structures and properties substantially differ from silica glass. These results suggest that the PDs provide a unified method that extracts greater depth of geometric information in amorphous solids than conventional methods.
The atomic configurations of amorphous solids are difficult to characterize. Because they have no periodicity as found in crystalline solids, only local structures have been analyzed in detail. Although short-range order (SRO) defined by the nearest neighbor is thoroughly studied, it is not sufficient to fully understand the atomic structures of amorphous solids. Therefore, medium-range order (MRO) has been discussed to properly characterize amorphous solids (13). Many experimental and simulation studies (47) have suggested signatures of MRO such as a first sharp diffraction peak (FSDP) in the structure factor of the continuous random network structure, and a split second peak in the radial distribution function of the random packing structure. However, in contrast to SRO, the geometric interpretation of MRO and the hierarchical structures among different ranges are not yet clear.
Among the available methods, the distributions of bond angle and dihedral angle are often used to identify the geometry beyond the scale of SRO. They cannot, however, provide a complete description of MRO because they only deal with the atomic configuration up to the third nearest neighbors. Alternatively, ring statistics are also applied as a conventional combinatorial topological method (2, 8, 9). However, this method is applicable only for the continuous random network or crystalline structures, and furthermore it cannot describe length scale. Therefore, methodologies that precisely characterize hierarchical structures beyond SRO and are applicable to a wide variety of amorphous solids are highly desired.
In recent years, topological data analysis (10, 11) has rapidly grown and has provided several tools for studying multiscale data arising in physical and biological fields (1116). A particularly important tool in the topological data analysis is persistence diagram (PD), a visualization of persistent homology as a 2D histogram (e.g., see Fig. 2). The input to the PD is given by an atomic configuration with scale parameters, and the output consists of various multiscale information about topological features such as rings and cavities embedded in the atomic configuration. Here, the atomic configurations are generated by molecular dynamics simulations in this article. Importantly, in contrast to other topological tools, PDs not only count topological features but also provide the scales of these features. Hence, PDs can be used to classify topological features by their scales and clarify geometric relationships among them; this is presumably the most desired function for deeper analysis of amorphous structures.
This article proposes a method using PDs for various amorphous solids in a unified framework. The method is applied to atomic configurations and enables one to study hierarchical geometry embedded in amorphous structures that cannot be treated by conventional methods. We first applied the method to silica glass as an example of the continuous random network structure and obtained the following results. (i) We found three characteristic curves in the PD of silica glass. These curves classify the SRO rings in the SiO4 tetrahedra and the MRO rings constructed by those tetrahedra. Furthermore, a hierarchical relationship among the SRO and MRO rings was elucidated. (ii) The PD reproduced the wavelength of the FSDP and clarified a real space origin of the FSDP. (iii) Each curve in the PD represents a geometric constraint on the ring shapes and, as an example, an MRO constraint on rings consisting of three oxygen atoms was explicitly derived as a surface in a parameter space of the triangles. Moreover, we verified that these curves are preserved under strain, indicating that the PD properly encodes the material property of elastic response. Next, as examples of the random packing structure, the Lennard-Jones (LJ) system and Cu-Zr metallic glass were studied by the PDs, and we clarified the following. (iv) These amorphous solids were also characterized well by the distributions of curves and islands in the PDs. (v) In the LJ system, the global connectivity of dense packing regions was revealed by dualizing octahedral arrangements. (vi) In Cu-Zr alloys, we found that the pair-distribution function defined by the octahedral region in the PD shows the split second peak. Furthermore, a relationship between the hierarchical ring structure and high glass-forming ability was discovered in Cu-Zr alloys.

PDs of Atomic Configuration

The input to PDs is a pair A=(Q,R) of an atomic configuration Q=(x1,,xN) and a parameter set R=(r1,,rN). Here, xi and ri are the position in 3 and the input radius for the ith atom, respectively. To characterize the multiscale properties in Q, we introduce a parameter α, which controls resolution, and generate a family of atomic balls Bi(α)={x3|xxiri(α)} having the radius ri(α)=α+ri2. We vary the radii of the atomic balls by α and detect rings and cavities at each α, where ααminmin{r12,,rN2}.
Let ck be a ring or cavity constructed in the atomic ball model B(α)=i=1NBi(α) at a parameter α. To be more precise, a ring [respectively (resp.) cavity] here means a generator of the homology H1(B(α)) [resp. H2(B(α))] with a field coefficient (17). Then, we observe that there is a value α=bk (resp. α=dk) at which ck first appears (resp. disappears) in the atomic ball model. The values bk and dk are called the birth and death scales of ck, respectively. The collection of all of the (bk,dk)2 of rings (resp. cavities) is the PD denoted by D1(A) [resp. D2(A)] for A (Fig. 1). It follows from the structure theorem of persistent homology (11) that the PDs are uniquely defined from the input. From this construction, (bk,dk) encodes certain scales of each ck. For example, in D1(A), bk indicates the maximum distance between two adjacent atoms in the ring ck, whereas dk indicates the size of ck.
Fig. 1.
Atomic balls (Top) and the PD D1(A) (Bottom). New rings appear at αi(i=2,3,4,5), and the dashed rings express the disappearance. This is a schematic illustration showing the rings on CP (red), CT (blue), CO (yellow), and BO (green) in silica glass. The large and small balls correspond to oxygen and silicon atoms, respectively.
In this article, our basic strategy is that we transform a complicated atomic configuration into PDs and try to identify meaningful shape information from specific distributions such as curves or islands in the PDs. Namely, we reconstruct characteristic atomic subsets from each distribution. To this aim, we compute the optimal cycle for each point (bk,dk)D(A) on the distribution. Mathematically speaking, for a given homology generator ck of (bk,dk), the optimal cycle is obtained by solving a minimizing problem in the representatives of ck under 1-norm (see refs. 18 and 19). Our method combining PDs with optimal cycles provides a tool to study inverse problems of PDs and effectively works in the geometric analysis of glass structures, as we will see shortly.
For a mathematically rigorous introduction of these concepts see Supporting Information or refs. 10 and 11. In this article, the PDs are computed by CGAL (20) and PHAT (21).

PDs for Continuous Random Network Structure

Fig. 2 shows the PDs D1(Aliq), D1(Aamo), and D1(Acry) of a liquid Aliq=(Qliq,Rliq), an amorphous Aamo=(Qamo,Ramo), and a crystalline Acry=(Qcry,Rcry) state of silica, respectively. Here, the horizontal and vertical axes are the birth and death scales, respectively, and the multiplicity of the PDs is plotted on a logarithmic scale. The configurations Qliq, Qamo, and Qcry are acquired by molecular dynamics simulations using the Beest–Kramer–Santen (BKS) model (22). The input radii R are set to be rO=1.275 Å and rSi=0.375 Å for each type of the atom (O or Si), which were determined from the first peak positions of the partial radial distribution functions of the amorphous configuration Qamo. The details of the molecular dynamics simulations and input radii are given in Supporting Information.
Fig. 2.
PDs of the liquid (Left), amorphous (Middle), and crystalline (Right) states with the multiplicity on the logarithmic scale. In the amorphous state, the three characteristic curves and one band region are labeled CP,CT,CO, and BO, respectively. The insets in D1(Aamo) show rings in the hierarchical relationship, where the red and blue spheres represent oxygen and silicon atoms, respectively.
We discovered that the PDs in Fig. 2 distinguish these three states. The liquid, amorphous, and crystalline states are characterized by planar (2-dim), curvilinear (1-dim), and island (0-dim) regions of the distributions, respectively. Here, the 0 and 2 dimensionality of the PDs result from the periodic and random atomic configurations of the crystalline and liquid states, respectively. Furthermore, we emphasize that the presence of the curves in D1(Aamo) clearly distinguishes the amorphous state from the others. This implies that specific geometric features of the rings generating the curves in D1(Aamo) play a significant role for elucidating amorphous states.
As shown in Fig. 2, D1(Aamo) contains three characteristic curves CP, CT, and CO and one band region BO, which are precisely characterized by using the invariance property with respect to the initial radius (15). These particular distributions, especially CO, start to become isolated near the glass transition temperature T=Tg (Supporting Information). Through further analysis of the persistent homology using optimal cycles, we found the following three geometric characterizations. (i) The rings on CP generate secondary rings on CT,CO, and BO (P is named after “primary”). That is, by increasing the parameter α, each ring on CP becomes thicker and starts to create new rings by pinching itself, and these newly generated rings appear on CT,CO, and BO (a schematic illustration of the pinching process is described in Fig. 1). This indicates a hierarchical structure from CP to CT,CO, and BO in the continuous random network. An example of the rings in this hierarchical relationship is depicted in the inset of D1(Aamo). (ii) The rings on CT are constructed by tetrahedra consisting of four oxygen atoms at the vertices with one silicon atom at the center. (iii) The rings on CO and BO are constructed only by the oxygen atoms (three and more, respectively). Recalling that the death scale indicates the size of rings, the rings on CT are classified as SRO, whereas those on CP, CO, and BO are classified as MRO.

Decomposition of FSDP

The FSDP observed in the structure factor S(q) (q1.5–2 Å1) has been supposed to be a signature of MRO, but its real space origin is still controversial (4). Here, we found that the distributions CP, CO, and BO of the MRO rings reproduce the q values of the FSDP fairly well. Moreover, we classified the MRO rings as a real space origin of the FSDP.
We first note that the death scales of the rings on CP, CO, and BO are determined only by the oxygen atoms, and this is directly verified by the PD computation. In addition, recall that α is the parameter controlling the radius ri(α)=α+ri2 of the ith atomic ball Bi(α), and the death scale α=dk indicates the size of the individual ring ck. More precisely, the ring ck disappears at α=dk by being covered up in the atomic ball model i=1NBi(α). Hence (dk)=2dk+rO2 measures the size of ck on CP, CO, and BO, where rO is the input radius of the oxygen atom.
From the aforementioned argument, we define a distribution
MAi(q)=1|Ai|(bk,dk)Aiδ(q2π(dk)),
where A1=CPCOBO, A2=CP, A3=CO, and A4=BO, and |Ai| is the number of the elements in Ai. Here, δ is the Dirac delta function, which is used to count the contribution of each MRO ring in q-space.
Fig. 3 shows the plots of MAi(q) and the structure factor S(q) around the FSDP. We found good agreement between the q values of the FSDP and the peak of MA1(q). This implies that the MRO rings composed of CO, CP, and BO are the real space origin of the FSDP. We also note that the distributions MA3(q),MA2(q),MA4(q) are located on the large, medium, and small q values and, hence, the rings in CO, CP, and BO provide a decomposition of the FSDP into those q values, respectively. It should be emphasized that the MAi(q) are derived from the configuration of the oxygen atoms only. This shows that the configuration of oxygen atoms plays a significant role in the FSDP. Furthermore, the invariant property (15) of (dk) for the rings on CP, CO, and BO induces that of MAi(q) under the choice of the input radius rO. This means that our analysis using PDs does not contain any artificial ambiguity of the input radii.
Fig. 3.
The distributions MAi(q) (i=1,2,3,4) and the structure factor S(q) around the FSDP.

Curves and Shape Constraints

The presence of curves CP, CT, and CO in D1(Aamo) clearly distinguishes the amorphous state from the crystalline and liquid states. We emphasize here that this characteristic property shows the constraints on the shapes of the rings induced by the normal directions of these curves.
For example, the shape of a ring on CO, which consists of three oxygen atoms (Fig. 4, Right), is determined by specifying the first and second minimum edge lengths d1 and d2 (d1<d2) and the angle θ between them, and hence is realized in a 3D parameter space. Then, the constraint for CO to be the curve requires that these three variables (d1,d2,θ) satisfy a certain relation f(d1,d2,θ)=0, and hence provides a restriction on the shape of the O-O-O triangles. Fig. 4 shows a plot of (d1,d2,θ) for the O-O-O triangles on CO in D1(Aamo), and we found a surface corresponding to this constraint f(d1,d2,θ)=0 in the parameter space. This demonstrates one of the medium-range geometric structures in the amorphous state. It is worth noting that this previously unidentified characterization of MRO of O-O-O triangles in 3D parameter space cannot be derived by separately analyzing the conventional distributions of lengths or angles, because each of them can only deal with the single parameter.
Fig. 4.
Plot of (d1,d2,θ) for the O-O-O triangles on CO.

Response Under Strain

The presence of the curves also indicates variations in the shapes of the rings. That is, by following each curve along its tangential direction and studying its rings, we can observe the deformation of the rings. It is reasonable to suppose that these variations are due to thermal fluctuations, and hence the deformations of the ring configurations along the curves are probably softer than those in the normal directions. Consequently, the response under strain is expected to follow the same shape constraints.
To verify this mechanical response of PDs, we performed simulations of isotropic compression for our amorphous state and computed the PD of the strained state. The strain is set to be 1% of the original volume, which is sufficiently small to satisfy a linear response relation with the stress. The three panels on the top in Fig. 5 show the contours of the histograms restricted on CP, CO, and CT for the original configuration Qamo of the amorphous state (black), for the strained state Qamostrain (red), and for the artificial 1% linear shrink Qamolinear of the original coordinates (blue), respectively. The contours of CO are depicted on the coordinates along the tangential direction. The bottom three panels in Fig. 5 show the projections of the histograms on their normal axes.
Fig. 5.
(Top) Contour plots of CP, CO, and CT for the original configuration Qamo of the amorphous state (black), the 1% strained state Qamostrain (red), and the artificial 1% linear shrink of the original coordinates Qamolinear (blue), respectively. (Bottom) Projections of the histograms on the normal directions.
Fig. 5 shows that the contours of the strained state shift along the original curves CP, CO, and CT in D1(Aamo) (i.e., downward in CP, leftward in CO, and almost fixed in CT, respectively). This is in contrast to a linear shrink, in which the contours simply move in the direction of decreasing both birth and death scales because of the uniform shrink of the system size. This strongly suggests that the configuration Qamostrain reflects the shape constraints and supports the expected mechanical response of PDs. We also note from the figures of Qamostrain that the contour of CT is almost fixed compared with those of CP and CO. Because CP and CO (CT) are classified as MRO (SRO), this implies that MRO is mechanically softer than SRO. Here, we remark that we can observe a similar behavior of the PDs for a shear deformation.
We have revealed that PDs encode information about elastic response, similar to how the radial distribution function encodes volume compressibility (23). It should also be emphasized that the curves or islands appear in the PDs of only the solids. This evidence suggests that these isolated distributions are related to the rigidity of the materials. This hypothesis follows from the fact that isolations represent geometric constrains reflecting mechanical responses. Future numerical and theoretical studies to unravel this relationship would be of great value.

PDs for Random Packing Structures

We next study the geometry of amorphous states close to random packing structures. In this case, we found that both D1 and D2 capture characteristic features of the amorphous structures. Fig. 6 shows the PDs Di(AcryLJ) and Di(AamoLJ) for i=1,2 of the monodisperse LJ system in the crystalline and the amorphous states, respectively. The input radii r=r1==rN in R are set to be zero, because changing r only causes translations of the PDs for the single component system. The details of the simulation are explained in Supporting Information.
Fig. 6.
PDs for the LJ system with the multiplicity on the logarithmic scale. Left and right panels correspond to D1 and D2, respectively. Top panels correspond to FCC crystal at T=0.1 and bottom panels correspond to amorphous state at T=103.
Similar to the case of silica, the crystalline structure is characterized by the island distributions in the PDs (top panels in Fig. 6). They correspond to the regular triangles in D1(AcryLJ) and the regular octahedra, the regular tetrahedra, and the quartoctahedra (24) in D2(AcryLJ) of the face-centered cubic (FCC) configuration. For the amorphous structure, the curves in D1(AamoLJ) and D2(AamoLJ) represent variations of triangles and tetrahedra, respectively. We also note that the isolation of the octahedral distribution is preserved well even in D2(AamoLJ). Its peak is separated from the curve of the tetrahedra (bottom right panel in Fig. 6), demonstrating the quantitative classification into two typical local structures of different-sized cavities.
In random packings, it is known that the atomic configuration can be divided into dense packing regions built from tetrahedra and the complement that patches those regions together (25). In particular, the network structure of the dense packing regions characterizes the global connectivity beyond MRO, which has not yet been investigated in detail. Note that D2(AamoLJ) clearly separates the dense packing regions as the deformation curve of the tetrahedra and identifies the octahedral island as the main component of the complement structure. From the Alexander duality in 3 (e.g., ref. 17), the connectivity of the dense packing regions can be studied by the cavities of the complement. Namely, we first extract all octahedra from the octahedral island in D2(AamoLJ) and construct a new set of points AoctLJ=(QoctLJ,RoctLJ) by putting a point at the center of each octahedron (Fig. 7, Left). Here, we set RoctLJ to be zero, because the SD of the octahedral sizes is very small.
Fig. 7.
(Left) Red balls express the octahedra in the amorphous structure of the LJ system, and the empty part corresponds to the dense packing regions. (Right) PD D2(AoctLJ) of the left figure with the multiplicity on the logarithmic scale.
Fig. 7, Right shows D2(AoctLJ), and we clearly observe that almost all cavities are located close to the diagonal. This means that the octahedral distribution rarely generates persistent cavities, and hence by the duality we can conclude that the dense packing regions mostly construct a giant connected network. It should be remarked that the dual treatment is much easier and computationally more efficient than directly studying the intricate huge network structures. We also emphasize that the analysis here of treating the octahedra as a new input is an iterative use of the PD method and can be an effective approach for studying multiscale geometry.
As an example of multicomponent systems, we also studied metallic glasses composed of Cu and Zr (26), in particular, focusing on Cu50Zr50 and Cu15Zr85, which display the different glass-forming ability. The PDs for these alloys are shown in Fig. 8. The input radii for Cu and Zr are set to be 1.30 Å and 1.55 Å, respectively, for the multicomponent PDs (Fig. 8, Top). These values are obtained by the same procedure as for the silica. Even in the multicomponent system, the PDs basically show similar behaviors to those of the LJ system. Specifically, the island distribution corresponding to octahedra appears in D2(AamoCu50Zr50) and the characteristic curves are also observed in D1(AamoCu50Zr50) and D2(AamoCu50Zr50).
Fig. 8.
PDs for Cu-Zr alloys with the multiplicity on the logarithmic scale. The left and right panels correspond to D1 and D2, respectively. In the top panels, Di(AamoCu50Zr50) are described. The middle and bottom panels show the PDs of the atomic configurations of only Zr atoms in the alloys. The middle panels correspond to Cu50Zr50 alloy, and the bottom panels correspond to Cu15Zr85 alloy. The blue and green spheres represent copper and zirconium atoms, respectively.
In the random packing structure, the split second peak of the radial distribution function has been supposed to be a signature of MRO (27). The shaded region of the radial distribution function in Fig. 9, Bottom shows the split second peak of Cu50Zr50. Then, we found that the pair-distance distribution of the atoms of generators in B[0.92,1.05]×[1.76,2.01]D2(AamoCu50Zr50) also shows a clear splitting of the distribution (black line in Fig. 9, Top) in the same length scale. Here, B is chosen to be a region around the octahedral distribution. Meanwhile the pair-distance distribution of generators other than B shows a slight change there (pink line in Fig. 9). This means that the generators around the octahedral distribution play a significant role for the split second peak. Therefore, similar to the FSDP in the silica, this result demonstrates that the PDs classify the length scale of MRO from other scales.
Fig. 9.
The normalized pair-distance distributions (Top) computed from D2(AamoCu50Zr50) and the radial distribution function for Cu50Zr50 alloy (Bottom) around the split second peak. The black line in the top panel was obtained by the pair-distances of the atoms in BD2(AamoCu50Zr50), whereas the pink line corresponds to those in the complement of B.
We also studied the PDs of only the Zr component. The PD D1 for Zr in Cu50Zr50 (Fig. 8, Middle Left) represents the existence of a hierarchical MRO structure similar to that of the silica in Fig. 2, whereas the PD D1 for Zr in Cu15Zr85 (Fig. 8, Bottom Left) does not show any hierarchical curves. An example of the hierarchical rings in Cu50Zr50 is depicted as the insets in the PD.
We here remark that Cu50Zr50 has higher glass-forming ability than Cu15Zr85 (28). Relationships between the glass-forming ability and the geometric structure are now actively studied (e.g., ref. 29). Then, this result suggests another possibility that the presence of the hierarchical MRO structure extracted from PDs is also related to the glass-forming ability of the alloy. In view of the results that the hierarchical MRO structure in PDs plays an important role for characterizing the glass states in SiO2, this statement seems to be reasonable. To understand the geometric and topological formation of the glass, it will be a great challenge to clarify this new perspective.

Conclusion

We have presented that PD is a powerful tool for geometric characterizations of various amorphous solids in the short, medium, and even further ranges. In this work, we have addressed two different types of amorphous systems: continuous random network and random packing structures. Both types of amorphous systems are characterized well by the existence of the curve and island distributions in the PDs. These specific distributions characterize the shapes of rings and cavities in multiranges, and the analysis using optimal cycles explicitly captures hierarchical structures of these shapes. We have shown that these shape characteristics successfully reproduce the FSDP for the continuous random network and the split second peak for the random packing and provide further geometric insights to them. Furthermore, the global connectivity of dense packing regions in the LJ system is revealed by the iterative application of the PD. For the binary random packing of the Cu-Zr metallic glass, we have also shown that the presence of the hierarchical MRO rings in the single component suggests the relationship with the glass-forming ability.
The methodology presented here can be applied to a wide variety of disordered systems and enables one to survey the geometric features and constraints in seemingly random configurations. Furthermore, because we investigated the mechanical response of the PDs, studying dynamical properties of materials using the PD method would be of great importance to understand the relationship between hierarchical structures and mechanical properties. We believe that further developments and applications of topological data analysis will accelerate the understanding of amorphous solids.

Computational Homology

Here we summarize the computational homological tools used in this article. For further mathematical and computational details, please refer to ref. 10.
Given a pair A=(Q,R) of an atomic configuration Q=(x1,x2,xN) and a set R=(r1,r2,,rN) of input radii, the persistence diagram captures hierarchical structures of the topological features such as rings and cavities in a one-parameter family of alpha shapes A(Q,R;α). Here, the alpha shape A(Q,R;α) is constructed as follows. Let 3=i=1NVi be the Voronoi decomposition for Q, where Vi is the Voronoi cell of xi. For each parameter α, we assign a ball Bi(α)={x3|xxiri(α)} with the radius ri(α)=α+ri2 for the ith atom. Then, the alpha shape A(Q,R;α) is defined by the dual of the cut balls Ci(α)=ViBi(α), i=1,,N. Namely, this is a polyhedron with the vertex set Q such that an edge |xixj| (triangle and tetrahedron, respectively) is assigned to A(Q,R;α) if and only if the intersection Ci(α)Cj(α) of the corresponding double (triple and quadruple, respectively) cut balls is nonempty (Fig. S1). In this article, the alpha shapes were computed using CGAL (20).
Fig. S1.
Alpha shape.
One of the important properties of the alpha shape is that the atomic ball model B(α)=i=1NBi(α) and the alpha shape A(Q,R;α) can be continuously deformed to each other. Hence, there is no loss of topological information in the use of the alpha shape model, which is much easier to analyze computationally than the atomic ball model. Furthermore, we can trace the hierarchical structures in B(α) by changing the parameter ααmin, where αmin=min{r12,rN2}. In other words, the parameter α controls the resolution from fine (ααmin) to crude (ααmin) geometric features.
Fig. S2, Left shows a schematic illustration of atomic ball models B(α) with the rings on CP (red), CT (blue), CO (yellow), and BO (green). Here, the large (small) balls correspond to oxygen (silicon). At each αi (i=2,3,4,5), a new ring appears, and the dashed rings express the rings that have disappeared. The hierarchical ring structures are observed at α3, α4, and α5, at which new rings are generated by pinching the primary ring.
Fig. S2.
Schematic illustration of atomic ball models B(α) showing the rings on CP (red), CT (blue), CO (yellow), and BO (green). Here, the large (small) balls correspond to oxygen (silicon). As we increase the parameter α, the primary ring generates the secondary rings in the hierarchy.
The persistence diagram is a 2D histogram recording the birth and death scales of the topological features in a one-parameter family {B(α)}α of atomic ball models. For example, let c be the yellow ring in Fig. S2. We observe that c first appears at α=α4 and first disappears at some value α*, where α4<α*<α5. Then, the birth and death scales b and d, respectively, of the ring c are determined by b=α4 and d=α*. The histogram of the birth and death pairs (b,d)2 for all of the rings appearing in the one-parameter family {B(α)}α is called the persistence diagram D1(A). Fig. S2, Right shows the persistence diagram of the left. From the definition of the persistence diagram, the birth scale indicates the maximum distance between adjacent atoms in the ring, whereas the death scale indicates the size of the ring. It should also be mentioned that persistence diagrams Di(A) can be similarly defined for any nonnegative integer i{0,1,2,}, and they capture higher-dimensional topological features (e.g., cavities for i=2). The persistence diagrams in this article were computed by PHAT (21).

MD Simulation for the Silica System

The atomic configurations of silica in a liquid Qliq, an amorphous solid Qamo, and a crystalline solid Qcry state were generated by the MD simulation. The system was composed of 2,700 silicon atoms and 5,400 oxygen atoms interacting via the BKS potential (22), and the 24-6 LJ potential is assigned to correct for the unphysical behavior during equilibration in the liquid state (30). The masses of the silicon and oxygen atoms were 27.98 and 15.99 g/mol, respectively. Despite its simplicity, the BKS model reasonably reproduces the static and dynamic properties of amorphous silica (31).
Starting from the beta-cristobalite structure, the crystalline configuration Qcry was obtained after the equilibration for 5 ps at 1 bar and 10 K. The MD simulation was performed with Parrinello–Raman barostat and Nosé–Hoover thermostat with the damping parameters τbar=0.5 ps and τth=0.8 ps, respectively. The time step was set to be 1.0 fs.
The liquid configuration Qliq was obtained after the equilibrating the system from a random initial configuration for 50 ps at 1 bar and 7,000 K. Here, the MD simulation was performed with the Andersen barostat and the Nosé–Hoover thermostat.
The amorphous configuration Qamo was obtained after the cooling MD simulation (31). Starting from the liquid configuration at the temperature 7,000 K, the heat-bath temperature was decreased linearly to the final temperature 10 K with dT/dt=1012 K/s. Throughout the cooling process the pressure was kept at 1 bar and we used the Nosé–Hoover thermostat with τth=1.0 ps. In the high-temperatures region 1,500 K T 6,000 K, we used the Anderson barostat with τbar=1.0 ps. In the low-temperatures region 10 K T 1,500 K, we used the Parrinello–Raman barostat with τbar=0.5 ps. Then, Qamo was obtained as the final configuration of the cooling process. We tested five cooling rates, dT/dt=1015,1014,1013,1012, and 1011 K/s, and checked that these choices do not affect our results qualitatively.
The set R=(r1,,rN) of input radii can be adjusted for each ith atom (i=1,,N). In our analysis, we chose the radii rSi and rO of silicon and oxygen using the first peak positions of the partial radial distribution functions of Qamo. The first peaks for SiO, OO, and SiSi appeared at rSiO=1.65 Å, rOO=2.55 Å, and rSiSi=3.15 Å, respectively. Then, we selected the smallest two peaks and determined rSi and rO by solving rSi+rO=rSiO and 2rO=rOO. This led to rO=1.275 Å and rSi=0.375 Å.
The glass transition temperature Tg and the melting temperature Tm are estimated through the enthalpy temperature curves. Starting from beta-cristobalite configuration, the system is heated with a heating rate dT/dt=0.1 K/ps from T=1 K to T=7,000 K (black curve in Fig. S3). Then, Tm is evaluated to be around 3,900 K (32). Subsequently, the cooling simulation is performed with a cooling rate dT/dt=10.0 K/ps (red curve). Then the glass transition temperature is evaluated as an intersection point between tangential lines of low (blue line) and high (green line) temperature regions of the curve and found to be Tg2,700 K (31). Fig. S4 shows PDs for 2,500 K (Left) and T=3,500 K (Right), and the curve CO starts to appear below Tg.
Fig. S3.
Enthalpy as a function of temperature for cooling and heating of silica system.
Fig. S4.
D1(A) below Tg (Left) and in the supercooled liquid region (Right).

MD Simulation for the LJ System and the Cu-Zr Alloy

The LJ system is composed of 4,000 monatomic particles interacting via LJ potential u(r)=4ε((σ/r)12(σ/r)6) with ε, σ, and mass equal to unity. The time step is set to be 0.005 in LJ units.
Starting from the FCC configuration, QcryLJ is obtained after equilibration at temperature T=0.1 and pressure p=1 with τth=1.0 and τbar=1.0. To obtain an amorphous structure QamoLJ, the equilibration at T=2 in the normal liquid phase has been performed for 2,000 steps at the constant number density 1.015 (33) with τth=0.5. Then, the cooling MD simulations have been performed for several cooling ratios from dT/dt=105 to 101. In Fig. S5, Left the enthalpy normalized by the number of particles is described as a function of T. As is observed, there is no glass transition in the LJ system. Therefore, we use an inherent structure for the liquid state as QamoLJ (black points in Fig. S5, Left).
Fig. S5.
Enthalpy as a function of temperature during cooling processes for the LJ system (Left) and Cu-Zr alloy (Right).
For the Cu-Zr alloy system, MD simulations have been performed using the embedded-atom method potential (34). The masses of Cu and Zr atoms were set to be 63.54 and 91.22 g/mol, respectively. The number of particles was 16,000, in which there were 8,000 of both copper and zirconium atoms for Cu50Zr50, 2,400 were copper atoms, and 13,600 were zirconium atoms for Cu15Zr85.
In this case, we can observe the glass transition (Fig. S5, Right) during a cooling simulation with the quench rate 1013 K/s (26) down to T=1 K keeping the pressure p=0 atm after the equilibration at T=17,000 K. Then, we use the final configuration at T=1 K as the amorphous configuration.
All MD simulations in this article were performed using LAMMPS (35).

Acknowledgments

We thank Mingwei Chen, Hajime Tanaka, Masakazu Matsumoto, and Daniel Miles Packwood for valuable discussions and comments. This work was sponsored by World Premier International Research Center Initiative, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also supported by Japan Science and Technology Agency (JST) CREST Mathematics Grant 15656429 (to Y.H.), Structural Materials for Innovation Strategic Innovation Promotion Program D72 (Y.H., A.H., and Y.N.), MEXT Coop with Math Program (K.M.), Japan Society for the Promotion of Science (JSPS) Grant 26310205 (to Y.N.), and JSPS Grant 15K13530 and JST PRESTO (to T.N.).

Supporting Information

Supporting Information (PDF)
Supporting Information

References

1
JD Bernal, The Bakerian Lecture, 1962 The structure of liquids. Proc R Soc Lond A Math Phys Sci 280, 299–322 (1964).
2
DE Polk, Structural model for amorphous silicon and germanium. J Non-Cryst Solids 5, 365–376 (1971).
3
JL Finney, Random packings and the structure of simple liquids. I. The geometry of random close packing. Proc R Soc Lond A Math Phys Sci 319, 479–493 (1970).
4
SR Elliott, Origin of the first sharp diffraction peak in the structure factor of covalent glasses. Phys Rev Lett 67, 711–714 (1991).
5
SR Elliott Physics of Amorphous Materials (Longman, 2nd Ed, New York, 1990).
6
S Susman, et al., Intermediate-range order in permanently densified vitreous SiO2: A neutron-diffraction and molecular-dynamics study. Phys Rev B Condens Matter 43, 1194–1197 (1991).
7
GN Greaves, S Sen, Inorganic glasses, glass-forming liquids and amorphizing solids. Adv Phys 56, 1–166 (2007).
8
SL Roux, V Petkov, ISAACS – Interactive structure analysis of amorphous and crystalline systems. J Appl Cryst 43, 181–185 (2010).
9
WH Zachariasen, The atomic arrangement in glass. J Am Chem Soc 54, 3841–3851 (1932).
10
H Edelsbrunner, J Harer, Computational Topology: An Introduction (Am Mathematical Soc, Providence, RI). (2010).
11
G Carlsson, Topology and data. Bull Am Math Soc 46, 255–308 (2009).
12
M Nicolau, AJ Levine, G Carlsson, Topology based data analysis identifies a subgroup of breast cancers with a unique mutational profile and excellent survival. Proc Natl Acad Sci USA 108, 7265–7270 (2011).
13
A Hirata, et al., Geometric frustration of icosahedron in metallic glasses. Science 341, 376–379 (2013).
14
H Kurtuldu, K Mischaikow, MF Schatz, Extensive scaling from computational homology and Karhunen-Loève decomposition analysis of Rayleigh-Bénard convection experiments. Phys Rev Lett 107, 034503 (2011).
15
T Nakamura, Y Hiraoka, A Hirata, EG Escolar, Y Nishiura, Persistent homology and many-body atomic structure for medium-range order in the glass. Nanotechnology 26, 304001 (2015).
16
L Kondic, et al., Topology of force networks in compressed granular media. Europhys Lett 97, 54001 (2012).
17
A Hatcher Algebraic Topology (Cambridge Univ Press, Cambridge, UK, 2001).
18
T Dey, et al., Optimal homologous cycles, total unimodularity, and linear programming. SIAM J Comput 40, 1026–1044 (2011).
19
EG Escolar, Y Hiraoka, Optimal cycles for persistent homology via linear programming. Optimization in the Real World –Towards Solving Real-World Optimization Problems. Mathematics for Industry, eds Fujisawa K, Shinano Y, Waki H (Springer, Tokyo), pp 79–96. (2015).
20
CGAL (2016) The Computational Geometry Algorithms Library. Available at www.cgal.org/.
21
PHAT (2014) Persistent Homology Algorithm Toolbox. Available at https://bitbucket.org/phat-code/phat.
22
GJ Kramer, Force fields for silicas and aluminophosphates based on ab initio calculations. Phys Rev Lett; van Beest BW; van Santen RA 64, 1955–1958 (1990).
23
JP Hansen, IR McDonald Theory of Simple Liquids (Academic, 3rd Ed, New York, 2006).
24
AV Anikeenko, et al., A novel Delaunay simplex technique for detection of crystalline neuclei in dense packings of spheres. Computational Science and Its Applications – ICCSA 2005. Lecture Notes in Computer Science, eds Gervasi O, Gavrilova ML, Kumar V, Laganà A, Lee HP, Mun Y, Taniar D, Tan CJK (Springer, New York), Vol 3480, pp 816–826. (2005).
25
R Yamamoto, M Doyama, The polyhedron and cavity analyses of a structural model of amorphous iron. J Phys F Metal Phys 9, 617–627 (1979).
26
SP Pan, et al., Origin of splitting of the second peak in the pair-distribution function for metallic glasses. Phys Rev B 84, 092201 (2011).
27
JA Barker, et al., Relaxation of the Bernal model. Nature 257, 120–122 (1975).
28
E Ma, Tuning order in disorder. Nat Mater 14, 547–552 (2015).
29
YQ Cheng, E Ma, Atomic-level structure and structure-property relationship in metallic glasses. Prog Mater Sci 56, 379–473 (2011).
30
Y Guissani, B Guillot, A numerical investigation of the liquid-vapor coexistence curve of silica. J Chem Phys 104, 7633 (1996).
31
K Vollmayr, W Kob, K Binder, Cooling-rate effects in amorphous silica: A computer-simulation study. Phys Rev B Condens Matter 54, 15808–15827 (1996).
32
TF Soules, et al., Silica molecular dynamic force fields – a practical assessment. J Non-Cryst Solids 357, 1564–1573 (2011).
33
G Monaco, S Mossa, Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. Proc Natl Acad Sci USA 106, 16907–16912 (2009).
34
MI Mendelev, et al., Using atomistic computer simulations to analyze x-ray diffraction data from metallic glasses. J Appl Phys 102, 043501 (2007).
35
S Plimpton, Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 117, 1–19 (1995).

Information & Authors

Information

Published in

Go to Proceedings of the National Academy of Sciences
Go to Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Vol. 113 | No. 26
June 28, 2016
PubMed: 27298351

Classifications

Submission history

Published online: June 13, 2016
Published in issue: June 28, 2016

Keywords

  1. amorphous solid
  2. hierarchical structure
  3. persistent homology
  4. persistence diagram
  5. topological data analysis

Acknowledgments

We thank Mingwei Chen, Hajime Tanaka, Masakazu Matsumoto, and Daniel Miles Packwood for valuable discussions and comments. This work was sponsored by World Premier International Research Center Initiative, Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also supported by Japan Science and Technology Agency (JST) CREST Mathematics Grant 15656429 (to Y.H.), Structural Materials for Innovation Strategic Innovation Promotion Program D72 (Y.H., A.H., and Y.N.), MEXT Coop with Math Program (K.M.), Japan Society for the Promotion of Science (JSPS) Grant 26310205 (to Y.N.), and JSPS Grant 15K13530 and JST PRESTO (to T.N.).

Notes

This article is a PNAS Direct Submission.

Authors

Affiliations

Yasuaki Hiraoka2,1 [email protected]
World Premier International Research Center Initiative–Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi 980-8577, Japan;
Takenobu Nakamura1
World Premier International Research Center Initiative–Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi 980-8577, Japan;
Akihiko Hirata
World Premier International Research Center Initiative–Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi 980-8577, Japan;
Emerson G. Escolar
World Premier International Research Center Initiative–Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi 980-8577, Japan;
Kaname Matsue
The Institute of Statistical Mathematics, Tachikawa, Tokyo 190-8562, Japan
Yasumasa Nishiura
World Premier International Research Center Initiative–Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi 980-8577, Japan;

Notes

2
To whom correspondence should be addressed. Email: [email protected].
Author contributions: Y.H., T.N., A.H., K.M., and Y.N. designed research; Y.H., T.N., A.H., E.G.E., K.M., and Y.N. performed research; Y.H. contributed new reagents/analytic tools; Y.H., T.N., and E.G.E. analyzed data; and Y.H., T.N., A.H., and E.G.E. wrote the paper.
1
Y.H. and T.N. contributed equally to this work.

Competing Interests

The authors declare no conflict of interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.


Citation statements




Altmetrics

Citations

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.

Cited by

    Loading...

    View Options

    View options

    PDF format

    Download this article as a PDF file

    DOWNLOAD PDF

    Get Access

    Login options

    Check if you have access through your login credentials or your institution to get full access on this article.

    Personal login Institutional Login

    Recommend to a librarian

    Recommend PNAS to a Librarian

    Purchase options

    Purchase this article to get full access to it.

    Single Article Purchase

    Hierarchical structures of amorphous solids characterized by persistent homology
    Proceedings of the National Academy of Sciences
    • Vol. 113
    • No. 26
    • pp. 7003-E3810

    Media

    Figures

    Tables

    Other

    Share

    Share

    Share article link

    Share on social media