Persistent homology-based descriptor for machine-learning potential of amorphous structures

High-accuracy prediction of the physical properties of disordered systems such as amorphous materials is significantly challenging in condensed-matter physics. Ab initio calculations, based on quantum mechanics, are versatile and highly accurate. However, the substantial computational cost limits the application of these calculations to realistic amorphous structures. In recent years, machine-learning potentials have emerged as an attractive method to resolve the trade-off between computational cost and accuracy.

Machine-learning potential involves the construction of a surrogate model of the relationship between atomic coordinates and physical properties, which provides an alternative to computationally intensive ab initio calculations. Various machine-learning models have been developed over the past decade, including Gaussian process regression,^1,2 artificial neural networks,^3–6 and their derivatives.^7–12 A key requirement of machine-learning potential models is the development of a suitable method for converting atomic coordinates into vector data.

Predictions derived from machine-learning potentials should adhere to physical laws. For example, the total energy of a system is invariant to the spatial translation, rotation, reflection, and permutation of chemically equivalent atoms. A direct approach to achieve this is to ensure that the vectorial inputs of the machine-learning potential remain invariant under these symmetry operations. Currently, two methods exist for creating vectorial representations of atomic coordinates, which are commonly known as descriptors. The first method involves using handcrafted feature representations such as the smooth overlap of atomic positions (SOAP)¹³ and atom-centered symmetry functions.¹⁴ The second method uses graph neural networks (GNN),^5,6,11,12 where essential structural characteristics are extracted by neural networks during training. Using the former descriptors, even a simple machine-learning model can adequately perform predictions. However, the optimal values of numerous hyperparameters in the descriptor, such as the cutoff, shape, and size of the basis functions, should be determined. The latter method has the advantage of not requiring descriptor hyperparameter tuning, but it requires increased model complexity.

This study proposes a descriptor that differs from the aforementioned approaches, which eliminates the need for hyperparameter tuning while predicting the energies of amorphous materials using a simple machine-learning model. In this study, we focus on persistent homology (PH),^15,16 a recently developed computational topology technique. PH effectively extracts topological and geometrical characteristics and has been applied to structural analyses in wide-ranging fields, from glass to bio-related materials.^17–22 In addition, recent studies have demonstrated that PH effectively identifies correlations between the physical and structural properties of materials.^23–29 

Herein, we provide a brief overview of the fundamental concepts of PH. In mathematics, a shape is considered a topological space, X, consisting of a set of points (in our case, atoms) and a topology. The topological characteristics of X are represented by holes in space that are encoded in an algebraic structure called the k-th homology group H_k(X). A nested sequence of topological space, X₁ ⊆ X₂ ⊆… X_n, results in the series $HkX1,HkX2⋯HkXn$⁠. The concept of PH involves monitoring the elements of H_k(X_i) as i increases, where i represents the scale and is often referred to as “time.” The procedure used to obtain this sequence is called filtration and is illustrated in Fig. 1.

During the filtration process for the first homology in a two-dimensional system, circles with a radius r are placed at the respective atoms. Subsequently, r is increased gradually. Consequently, the circles begin to intersect, and then, edges are set between the centers of the intersecting circles. The growing network formed by these edges defines the sequence of the topological space, X_i. A closed path formed by these edges corresponds to a topological feature called a cycle. As the radius continues to increase, the closed path becomes completely covered by circles, which is interpreted as the conversion of the cycle into another type of topological feature called a boundary. The features of PH are represented by pairs of birth and death times at which cycles appear and are converted into boundaries. Filtration can be generalized into three or more dimensions using spheres or hyperspheres. A two-dimensional visualization of birth–death time pairs is known as a persistence diagram (PD).

Based on this definition, the PD captures information regarding the bonding state and distribution of atoms and remains invariant to the spatial translation, rotation, reflection, and permutations of chemically equivalent atoms. Moreover, PD does not require any hyperparameter to be adjusted manually.

In this study, we demonstrate that PH can be used as an input for the machine-learning potential of amorphous structures. The specific target is amorphous carbon (aC). Initially, we demonstrate that the normalized two-dimensional histogram derived from the PD can predict the average energy per atom of aC at various densities, even when using a Ridge model. Moreover, we demonstrate that the prediction accuracy is improved by utilizing a convolutional neural network (CNN) model. Subsequently, we compare the PH descriptors, conventional handcrafted descriptors, and latent space of the GNN-based architectures. SOAP descriptors are selected for the conventional handcrafted descriptor, whereas SchNet is used for the GNN-based architecture. We reveal the differences in how each descriptor captures structural characteristics. The 2D maps obtained from the dimensionality reduction of the PH descriptors resemble those of the SchNet descriptors but differ from those of the SOAP descriptors. Based on these results, we conclude that PH can construct descriptors with characteristics similar to those of a latent space in a GNN, without deep learning and hyperparameter tuning.

The aC dataset was generated by ab initio calculations using the VASP software.^30–33 We utilized the local-density approximation (LDA) exchange-correlation functional and PAW potential for carbon. Melt-quench simulations were performed to generate amorphous and liquid-state structures. A simple cubic lattice with 216 carbon atoms was selected as the initial state. Simulations were conducted at densities of 1.5, 1.7, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, and 3.5 g/cm³ to generate different structures. The NVT ensemble^34,35 was used for all the melt-quench simulations, and density was adjusted by varying the size of the simulation cell. A time step of 1 fs was used for the simulations. For all densities, only the Γ points were sampled in k-space. To increase the structural diversity, six independent simulations were performed.

During the melt-quench simulations, the temperature was increased from 300 to 9000 K in 2 ps to melt the carbon. Equilibrium molecular dynamics (MD) was conducted at 9000 K for 3 ps to form the liquid state, followed by a decrease in temperature to 5000 K in 2 ps and equilibration for 2 ps. Finally, the temperature was decreased from 5000 to 300 K in 2 ps to generate the amorphous structure. During this process, 30 snapshots were obtained from the equilibrium MD trajectory at 9000 K, 100 from the cooling process between 9000 and 5000 K, 25 from the equilibrium MD trajectory at 5000 K, and 100 from the cooling process between 5000 and 300 K, yielding 16 830 data points. The coordination number analysis of this dataset confirms that the carbon atom neighbor environments are adequately diverse (Fig. S2).

Additionally, the data for diamond structures containing 216 atoms at densities of 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, and 3.5 g/cm³ were prepared. Further data on the diamond structures were obtained from the 80 snapshots of the 2 ps equilibrium MD trajectory at 300 K, yielding 560 data points.

To validate the predictions for larger structures, we generated the data for 512-atom systems using the same procedure as that for the 216-atom systems. A single simulation was conducted for each density. The number of data points was 2805 for the amorphous and liquid states. These datasets can be downloaded from Zenodo (DOI: https://doi.org/10.5281/zenodo.7905585).

In this study, we utilized the HomCloud code^36,37 to analyze the PH of the amorphous models. Herein, we focused on the ring structures that correspond to the cycles from one-dimensional homology. The birth and death times of the respective cycles were determined by alpha complex filtration. The basic procedure of alpha complex filtration is the same as that depicted in Fig. 1, although the area covered by each sphere is restricted within the Voronoi cell. The radii of the spheres were continuously and uniformly increased from zero until the filtration was completed (i.e., all cycles were converted to boundaries). In accordance with the convention, we denote the birth and death times using the squared values of the radii of the spheres employed during the filtration process. Consequently, the results of these measurements are given in the units of Ų.

To create the input for the machine-learning potential, we converted the PDs into two-dimensional histograms using a 128 × 128 mesh. The square region of [0.0, 8.0] _birth × [0.0, 8.0] _death in the histogram was targeted. Unlike the previously proposed methods,^38,39 our method does not either apply a small weight to or remove the birth–death time pairs close to the diagonal line in PD. This is because we expect that the importance of each birth–death time pair will be automatically recognized during the machine learning process. Because the PD represents the number and size of ring structures, the histogram values depend on the system size. Consequently, the histogram was transformed into a probability distribution, and z-score was standardized to ensure that the descriptors were independent of model size. Because of the small size of the model, periodic boundary conditions (PBC) were applied when calculating the PDs. We discuss the limitations originating from the use of PBC in Sec. 3 of the supplementary material.

By tracing the filtration procedure, we can identify the candidates of the cycle corresponding to each birth–death pair. This enabled inverse analysis of PH, which was used to determine the atomic arrangements corresponding to the low- and high-energy regions. In some cases, several candidates exhibit equivalent topological characteristics. Various approaches^40,41 exist for selecting a representative cycle for these cases. In this study, the stable volume method⁴¹ was selected, which provides the tightest cycles with robustness against noise.

The SOAP descriptors were calculated using DScribe.⁴² The hyperparameters for the SOAP descriptor include the cutoff radius (rcut), number of radial basis functions (nmax), maximum degree of spherical harmonics (lmax), and standard deviation of the Gaussian function used to expand the atomic density (sigma). Optuna⁴³ was used to search for the optimal hyperparameter values, which were determined to be rcut = 4.1, nmax = 8, lmax = 9, and sigma = 0.16 (see Fig. S1 in the supplementary material). These hyperparameters resulted in a SOAP descriptor dimension of 360. The average value across the axis representing the number of atoms was used to evaluate the descriptor space.

By updating the learnable weight parameters in the architecture with a combination of the four components during the training procedure, the network generated descriptors that capture the structural characteristics to determine the total energy. In this study, the output from the last interaction block of the trained SchNet model was considered the descriptor generated by SchNet. Because the descriptors were computed for each atom, the resulting descriptors formed an n × D matrix, where D represents the descriptor dimension that was set to 100. To evaluate the descriptor space, the average value of the D-dimensional vector along the axis of the number of atoms was used.

Ridge regression for the total energies, based on the PD-generated descriptors, was conducted using the Scikit-learn package.⁴⁵ An intercept term was included in all the regression models. Because all the input values were non-negative, the regression coefficient was directly interpreted. The regularization parameter was set to 200.0.

The CNN model was implemented using PyTorch.⁴⁶ The model comprises three convolution layers followed by max-pooling layers and two fully connected layers. The ReLU function was used as the activation function in the model. The channel size of the first two convolution layers was 64, whereas that of the last layer was 32. The filter size was 3 × 3 for all the convolution layers.

During the training process, the loss function was evaluated based on the mean squared error between the predicted and actual mean energy per atom. The weights and biases in the convolution and fully connected layers were optimized using stochastic gradient descent (SGD) with a Nesterov momentum of 0.9 and a weight decay of 0.001. The milestones were set at 200, 400, 800, and 1200, and the number of epochs was set to 1600. When the epoch reached the milestones, the learning rate was reduced by a factor of 0.5. The initial value of the learning rate was set to 0.0015.

Neural network potentials based on SOAP descriptors were implemented using PyTorch. The model consists of two fully connected hidden layers with 40 nodes in each layer. During the training process, the loss function was evaluated based on the mean squared error between the predicted and actual mean energy per atom. The weights and biases in the model were optimized using the Adam optimizer. The milestones were set at 200, 500, and 1000, and the number of epochs was set to 2000. When the epoch reached the milestones, the learning rate was reduced by a factor of 0.5. The initial value of the learning rate was set to 0.001.

First, we present the typical PDs for liquid and amorphous carbon. Figures 2(a) and 2(d) show samples of liquid and amorphous structures with a density of 2.2 g/cm³. The PDs derived from these structures are presented in Figs. 2(b) and 2(e). The colors and symbols in these PDs signify the number of vertices of the cycle corresponding to each point, which is determined by stable volume analysis.

The three-vertex cycle with a long birth–death time that is observed only in the amorphous state originates from the pore structure present under low- and medium-density conditions. The local structure corresponding to this cycle is depicted in the inset of Fig. 2(e).

A remarkable difference between the liquid and amorphous states emerges in the birth–death time distributions for cycles with more than five vertices. We present the typical local structure of the five-, six-, and seven-vertex cycles in both liquid and amorphous states in the insets of Figs. 2(b) and 2(e), along with their birth–death times. Compared with the cycles in the amorphous state, those in the liquid state exhibit structures with greater distortion.

In the PD of amorphous structures, these cycles form a vertical line in the birth time range of 0.5–0.7. Such a structure is not observed in the case of liquid. This corresponds to the fact that medium-range order (MRO) exists in amorphous structures but not in the liquid state. This correspondence can be clearly observed in the radial distribution function using ring structures corresponding to five-vertex and six-vertex cycles [Figs. 2(c) and 2(f)].

Although similar differences in the PDs between amorphous and liquid structures also appear in the cases of lower and higher densities, it was found that the area of birth–death times distribution varies with the density (Figs. S5 and S6). These analysis results demonstrate that it is possible to detect the presence or absence of MRO and pore structures from the PD.

Next, we demonstrated that the mean energy per atom in amorphous/liquid carbon can be predicted using PDs. As described in Sec. II, we constructed descriptors based on PH by converting the PDs into two-dimensional histograms. These PH descriptors were used as inputs for the Ridge regression and CNN models. Figures 3(a) and 3(b) show the training and test results of the respective models for the 216-atom aC system. In this study, 80% of the data were used as training data, and the remaining 20% were used as test data. Even using a simple Ridge regression model, the mean energy per atom in both amorphous and liquid carbon was predicted by the PH descriptors. The root mean squared error (RMSE) for the test data was 59.2 meV/atom. The prediction accuracy improved by utilizing the CNN model, and the RMSE for the test data decreased to 40.7 meV/atom.

As shown in Fig. S7, both ridge regression and CNN predictions were valid when the diamond structures were included in the dataset. This result highlights that the PH descriptors effectively capture the correlation between the atomic structures and the energies of various structures.

Machine-learning potentials should maintain their prediction accuracy when applied to a system larger than that of the training data. Therefore, we predicted the mean energy per atom in aC for a system with 512 atoms. As shown in Figs. 3(c) and 3(d), both the Ridge regression and CNN models accurately predicted the energies based on the PH descriptor.

Subsequently, we compared these results with those obtained using other methods. SOAP and SchNet were used as examples of a conventional handcrafted representation and GNN model, respectively. Figure 4 shows the results of the training and test for the 216-atom aC system using these methods. Among the investigated models, the combination of the SOAP descriptor and neural network exhibited the highest accuracy, and the RMSE for the test data was 16.4 meV/atom. The prediction accuracy is sensitive to the selection of hyperparameters. As shown in Fig. S8, different hyperparameters can reduce prediction accuracy. For the SchNet model, the RMSE for the test data was 30.5 meV/atom, which falls between those obtained from the SOAP+NN and PH descriptor+CNN models. Note that all these prediction accuracies achieved by machine-learning potentials are superior to those obtained from empirical potentials (Fig. S9).

Beyond the comparison of the prediction accuracy, we investigated how these descriptors capture the structural characteristics of amorphous structures using dimensionality reduction. As described in the Methods section, the average values of the descriptors for the atoms were used to evaluate the descriptor space constructed using SOAP and SchNet. We note that this procedure did not cause a significant loss of information. As shown in Fig. S10, the mean energies per atom were accurately predicted by the averaged descriptors using a neural network. Figure 5 shows the t-SNE analysis results. The colors in Figs. 5(a)–5(c) correspond to the density of aC, whereas those in Figs. 5(d)–5(f) correspond to the energy of aC. In the t-SNE 2D maps of the SOAP descriptor space, descriptors were grouped into discrete island-like structures based on their densities [Fig. 5(a)].

Similarly, the t-SNE 2D maps for SchNet exhibited island-like structures; however, the distribution of each island was broader than that in the case of SOAP, and several islands partially overlapped [Figs. 5(b) and 5(e)]. The 2D visualizations of the PH descriptor space exhibited the most continuous distributions with respect to both energy and density, as shown in Figs. 5(c) and 5(f). These characteristics of the 2D visualizations were valid when different hyperparameters values were used in the t-SNE analysis (Figs. S11 and S12). Moreover, principal component analysis (PCA) exhibited similar results. As shown in Fig. 6, the 2D visualization of the SOAP descriptor using PCA revealed a discrete distribution, whereas those of the PH and SchNet descriptor spaces were continuous.

Notably, the above-mentioned differences between the SOAP and PH/SchNet descriptors were independent of the hyperparameters used in the SOAP descriptors. As shown in Fig. S13, the dimensionality reduction of the SOAP descriptor used in Fig. S8 had a discrete distribution similar to those shown in Figs. 5 and 6. This indicates that the continuous distribution originated from the additional information captured by the PH and SchNet descriptors but not by the SOAP descriptor. One promising candidate for this type of information is the topological structure, which is intrinsically incorporated into PH and GNN.

Based on these results, we concluded that the extraction of structural characteristics by PH is similar to that by GNN, but different from that by SOAP. Thus, PH can construct a descriptor similar to that of a GNN without deep learning and hyperparameter tuning.

These characteristics of the PH descriptor can be used to reduce the complexity of machine-learning models for predicting the energy of amorphous materials. Another advantage of this descriptor is its prediction interpretability. In the case of Ridge regression, we can estimate the contribution of each birth–death pair from the value of the Ridge coefficient. Figure 7(a) shows the projection of the Ridge coefficient onto the mesh grid used to generate the 2D histogram of the PD. The birth–death pairs in the dark blue region (the large-negative-coefficient region) correspond to the local structures that reduce the energy. By contrast, those in the dark red region (the large-positive-coefficient region) correspond to the local structures that increase the energy. In this study, we define the birth–death pairs in the region where the Ridge coefficient is larger (lower) than 0.003 (−0.0045) as the local structures with high (low) energy.

We then determined the local structures that correspond to the high- and low-energy regions obtained by inverse analysis based on stable volumes. Figure 7(b) shows a visualization of the high- and low-energy regions in the structure obtained from the equilibrium MD trajectory at 5000 K with a mass density of 2.4 g/cm³.

The same analysis was conducted on the remaining 24 samples obtained from the same trajectory. Figures 7(c) and 7(d) show the distributions of the birth–death times of the cycles assigned to the high- and low-energy regions. Examples of the ring structures are shown in the insets.

The distribution of the birth–death pairs for the low-energy region overlaps with the vertical line structure in Fig. 2(e), which characterizes the evolution of the MRO in the amorphous phase. By contrast, the distribution of the birth–death pairs for the high-energy region overlaps with that obtained from the liquid state, characterized by more distorted ring structures.

We also analyzed the histograms of the bond lengths and angles in the low- and high-energy regions, as shown in Figs. 7(e) and 7(f). The histograms clearly show differences between the local structures assigned to the low- and high-energy regions. In the bond length histogram, the low-energy structures exhibited a sharp peak at ∼1.5 Å, whereas the peak for the high-energy structures was significantly broader. The histograms of the bond angles for the low- and high-energy structures differ in the range between 50° and 80°. These differences in the histograms indicate that high-energy structures exhibit greater randomness in terms of bonds and angles.

These analyses indicate that the regression model for energy prediction based on PDs focuses on the following two local structures as descriptors of the low- and high-energy regions: (1) large ring structures with low variation in bond length and angles, corresponding to the evolution of MRO, and (2) ring structures with greater randomness. These results verify the ability of PH to effectively extract the correlation between geometric and topological structures and their energies.

In this study, we demonstrated that the descriptors obtained from PH can predict the mean energy per atom in aC using a simple ridge regression model. The prediction accuracy was improved by utilizing a CNN model. Furthermore, we compared the results with those obtained using other methods. As examples of a conventional handcrafted representation and a GNN model, we used SOAP and SchNet, respectively. We demonstrated that the PH descriptor and latent space in the trained SchNet model exhibited similar characteristics. We visualized this point using the dimensionality reduction of the descriptor spaces obtained from PH, SOAP, and SchNet. The t-SNE 2D maps of the SOAP and SchNet descriptor spaces exhibited island-like structures. However, the distribution of each island for SchNet was significantly broader than that for SOAP, and several islands partially overlapped. The 2D visualizations of the PH descriptor space exhibited continuous distributions. A similar result was observed for the 2D maps obtained by PCA. The difference between the SOAP and PH/SchNet descriptor spaces originates from the information on the topological structure, which is intrinsically incorporated into PH and GNN.

Another advantage of the PH descriptor is the interpretability of prediction. Local structures corresponding to the high- and low-energy regions can be determined by utilizing the inverse analysis technique for PH.

In summary, the PH descriptor provides several advantages including a simple model architecture for predicting energy, the absence of hyperparameter tuning in descriptor construction, and the interpretability of prediction. However, several of these aspects require further investigation. First, a method to improve prediction accuracy is required because the accuracy of prediction by PH is not superior to that of the other methods. One possible reason for this is that PH may lack detailed geometrical information. The complementary use of PH and conventional descriptors may be a promising approach. Second, to apply this approach to large-scale molecular dynamics simulations, an extension to the prediction of the forces acting on the atoms is necessary. This would necessitate the development of PH-based techniques for extracting local topological characteristics of atom environments and the utilization of the PD derivatives with respect to atomic coordinates.

Another aspect is the construction of PH descriptors for multielement systems. In the case of monoelement systems, such as the amorphous carbon studied here, we can obtain descriptors suitable for energy prediction by applying a filtration process in which the radii of all spheres are uniformly increased. However, when the system includes different elements, special measures are required to extract structural features that reflect the variations in the interactions between atoms of the same and different species. This problem may be solved simply by varying the initial radii of the spheres used in filtration according to the atomic species. For our proposed CNN-based model, it is also possible to combine the PD data for the partial atomic coordinates of each element, similar to the RGB channels in the image data. The exploration of these possibilities is expected to improve machine-learning potentials based on descriptors that depict both topological and geometrical information.

1.Hyperparameter optimization for SOAP descriptor2.Diversity of neighboring condition of carbon atoms in the dataset3.Limitations originating from the use of periodic boundary conditions4.Structures and PDs for liquid and amorphous carbon with low and high densities5.Prediction results by PH descriptor for the dataset including crystal data6.Prediction results by SOAP+NN model with hyperparameters different from optimal values7.Calculation results of the energies using empirical potentials8.Results of prediction using average values of SOAP and SchNet descriptors9.Results of t-SNE with perplexity = 10.010.Results of t-SNE with perplexity = 50.011.Dimensional reduction analyses for SOAP descriptor space with different hyperparameters

This study was supported by PRESTO Grants (Grant Nos. JPMJPR1923 and JPMJPR2198), JST, and KAKENHI Grant (Nos. 23H04470, 23H04100, 22H05106, 21H01816, 21H05552, 19H02544, 19H00834, and 20H05884), MEXT, Japan. The calculations were performed using a computer facility at the Research Center for Computational Science (Okazaki, Japan).

The authors have no conflicts to disclose.

Emi Minamitani: Conceptualization (lead); Funding acquisition (equal); Investigation (lead); Methodology (lead); Software (equal); Validation (lead); Writing – original draft (lead); Writing – review & editing (equal). Ippei Obayashi: Conceptualization (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal). Koji Shimizu: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Satoshi Watanabe: Conceptualization (supporting); Funding acquisition (equal); Investigation (supporting); Methodology (supporting); Writing – original draft (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.7905585. Raw data were generated at Research Center for Computational Science (Okazaki, Japan). Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.