††thanks: A.R.T. and S.U. contributed equally to this work.††thanks: A.R.T.
and S.U. contributed equally to this work.

# Graph theory-based structural analysis on density anomaly of silica glass

Aik Rui Tan  Department of Materials Science and Engineering, Massachusetts
Institute of Technology  Shingo Urata<EMAIL_ADDRESS>Technology General
Division, Planning Division, AGC Inc.  Masatsugu Yamada  Innovative Technology
Laboratories, AGC Inc.  Rafael Gómez-Bombarelli<EMAIL_ADDRESS>Department of
Materials Science and Engineering, Massachusetts Institute of Technology

###### Abstract

Understanding the structure of glassy materials represents a tremendous
challenge for both experiments and computations. Despite decades of scientific
research, for instance, the structural origin of the density anomaly in silica
glasses is still not well understood. Atomistic simulations based on molecular
dynamics (MD) produce atomically resolved structure, but extracting insights
about the role of disorder in the density anomaly is challenging. Here, we
propose to quantify the topological differences between structural
arrangements from MD trajectories using a graph-theoretical approach, such
that structural differences in silica glasses that exhibit density anomaly can
be captured. To balance the accuracy and speed of the MD simulations, we
utilized force matching potentials to generate the silica glass structures.
This approach involves casting all-atom glass configurations as networks, and
subsequently applying a graph-similarity metric (D-measure). Calculated
D-measure values are then taken as the topological distances between two
configurations. By measuring the topological distances of configurations in
silica glass simulated structures across a range of temperatures, distinct
structural features could be observed at temperatures higher than the fictive
temperature. In addition to quantifying structural changes in the full
simulation box, we compared topological distances between local atomic
environments in the glass and crystalline silica phases. Evidence from this
approach suggests that more coesite-like local structures, which are less
symmetric, emerge in silica glasses when density is at a minimum during the
heating process.

##  I Introduction

Research in the physical sciences is increasingly leveraging data-driven
approaches. Materials Informatics (MI), for instance, applies the principles
of data science to materials research and discovery and is expected to greatly
reduce the risk and time for materials deployment[1, 2, 3, 4]. Despite this
promising outlook, experimental data in materials science is inconsistently
reported and typically sparse and scarce because of the high costs of
experiments. High-throughput simulations can be used to replace experiments or
to enhance the robustness and accuracy of MI applications by creating
informative descriptors[5, 6, 7, 8, 9]. A key part of MI is the generation of
quantitative structure-property relationships, where structural features that
encode atomic environments are related to observed properties. A wide variety
of approaches to process atomistic structure have been proposed.[10, 11, 12,
13, 14, 15].

Unlike small molecules, molecular liquids or crystalline materials, amorphous
materials are structurally non-unique, which is a standing bottleneck to
define descriptors that are capable of encoding the topological and chemical
landscapes in glass structures. In the case of oxide glasses, common global
descriptors include bond order of network formers, ratio of bonding oxygen and
non-bonding oxygen (BO/NBO), distributions of Qn speciation (where $n$ is
number of bridging oxygen connected to a network-former cation), ring sizes,
radial distribution functions and structure factors. More recently, persistent
homology has attracted attention as a novel algorithm to analyze amorphous
topology based on the vacancy distribution [16, 17, 18]. In terms of local
structural environments, descriptors such as tetrahedral order of SiO4 [19],
Qn speciation, and bond order for atoms are used to represent the short-range
coordination, while medium-range descriptors are still mostly limited to local
strain [20] and potential energy [21, 22].

However, tools are lacking to describe glass structure at the intermediate-
range. Here, we propose a graph-based approach to study the order/disorder of
glasses beyond local environments by comparing connectivity networks including
second and higher order neighbors. Expressing the atomic environment as a
discrete graph provides a rich signal that can easily distinguish dissimilar
(non-isomorphic) environments. However, quantifying the distance between two
arbitrary graphs is not computationally trivial. Here, we apply a continuous
graph dissimilarity metric known as D-measure [23] which ranges from 0 to 1
(isomorphic). D-measure has been applied to analyze social networks [24, 25],
community structures [26], brain networks [27], and in the realms of materials
science to characterize structural transition of zeolite nanoporous materials
[28].

This work applies D-measure to amorphous materials for the first time. We
choose silica glass as a representative oxide glass since it is in principle
simple, but exhibits non-trivial behaviors. Sharing the same tetrahedral order
and dynamic anomalies as water, silica exhibits a phenomenon known as density
anomaly, where its density reaches a minimum at temperatures higher than the
melting point, and then increases with temperature [29]. Silica also exhibits
typical nonlinear behaviors such as mechanical anomaly [30] and fragile-to-
strong transition of diffusion behavior [31]. These peculiar properties of
silica have been studied using molecular dynamics (MD) simulations [30, 31,
32, 33, 34, 35, 36, 37] using a plethora of interatomic potentials such as the
BKS potential[38], TTAM potential[39], and others[22, 30, 34, 35, 36, 40].
Since fluorine doping is known to reduce the density anomaly of silica glass
we also analyze MD simulations of fluorine-doped silica glasses.

Computational details on D-measure and MD simulations are provided in Section
II. Section III examines applicability of D-measure to analyze the global and
local structures in silica and F-doped silica glasses. Section IV summarizes
the concluding remarks.

##  II Computational Method

###  II.1 D-measure

D-measure, $D(G,G^{\prime})$ evaluates the dissimilarity between two
undirected graphs, $G$ and $G^{\prime}$. First, the probability distribution
function (PDF) of node-distances for every node is generated and then PDFs are
compared through the following definitions [23].

$\displaystyle D(G,G^{\prime})$ $\displaystyle=$ $\displaystyle
w_{1}\sqrt{{J(\mu_{G},\mu_{G^{\prime}})\over
log2}}+w_{2}\left|\sqrt{\textrm{NND(G)}}-\sqrt{\textrm{NND(G')}}\right|$
$\displaystyle+$ $\displaystyle{w_{3}\over 2}\biggl{(}\sqrt{{J(P_{\alpha
G},P_{\alpha G^{\prime}})\over\log 2}}+\sqrt{{J(P_{\alpha G^{C}},P_{\alpha
G^{C^{\prime}}})\over\log 2}}\biggr{)},$

where

$\displaystyle J({\bf P}_{1},...,{\bf P}_{N})$ $\displaystyle=$
$\displaystyle{1\over
N}\sum_{i,j}p_{i}(j)\log\biggl{(}{p_{i}(j)\over\mu_{j}}\biggr{)}$ (2)
$\displaystyle\mu_{j}$ $\displaystyle=$ $\displaystyle{1\over
N}\sum_{i=1}^{N}p_{i}(j)$ (3) $\displaystyle\textrm{NND}(G)$ $\displaystyle=$
$\displaystyle{J({\bf P}_{1},...,{\bf P}_{N})\over\log(d+1)}$ (4)

$G^{C}$ and $G^{C^{\prime}}$ indicate the complement of $G$ and $G^{\prime}$.
The distance distribution is defined as ${\bf P}_{i}=\\{p_{i}(j)\\}$, where
${p_{i}(j)}$ is the fraction of nodes connected to node $i$ at distance $j$.
Accordingly, a graph composed of $N$ nodes is written as {${\bf P}_{1}$, …
,${\bf P}_{N}$}. The first term of Eq. (II.1) compares graph distributions
representing averaged node’s connectivity patterns using the Jensen-Shannon
divergence (Eq. 2). The second term evaluates heterogeneity of the nodes in
terms of their connectivity profiles via the network node dispersion term, NND
(Eq. 4). The third term measures centrality of each node by considering the
span of nodes both directly and indirectly connected.

$D(G,G^{\prime})$ is constrained within $(0,1)$ under the limitation that
$w_{1}+w_{2}+w_{3}=1$. Two isomorphic networks return $D(G,G^{\prime})=0$, but
two non-isomorphic networks will not necessary return non-zero D-measure.
Furthermore, graph sizes of $G$ and $G^{\prime}$ need not be the same.

In this work, the weights are set to be $w_{1}=w_{2}=0.45$ and $w_{3}=0.1$ as
the original study suggested [23]. To simplify the graphs, only silicon atoms
are considered as nodes. Two silicons are defined to be connected when they
are linked through a bridging oxygen atom, with a Si-O cutoff distance of 1.9
Å.

###  II.2 Molecular Dynamics simulations

Table 1 summarizes all the models constructed by MD simulations. The LAMMPS
package [41] was employed for the MD simulations. All simulations were
conducted in an isothermal-isobaric ensemble (NPT) using Nosé–Hoover
thermostat [42] and barostat [43]. The F-doped silica glass models were
constructed by replacing 1, 3, 5 and 10 % of oxygen of the silica glass model
with fluorine. They are denoted as SiO2-Fx (x = 1, 3, 5, 10). The glass models
were obtained using the melt-quenching method [44] where initial glass
structures were melted at 3500 K for 2.5 ns, and 5 different configurations
were extracted every 500 ps to act as the starting configurations for 5
different MD trajectories. Each of the 5 configurations were cooled down from
3500 K to 300 K at a cooling rate of 1 K ps-1, followed by equilibration at
300 K for 500 ps each. They were then heated up again to 3500 K at a rate of
0.5 K ps-1. Individual glass models were recorded every 5 K (10 ps) and each
model was then equilibrated for 10 ps at the corresponding temperatures. Only
data from the last 3 ps was used to average the properties, such as density
and potential energy. These five independent replicas were examined for all
the glass models to verify reproducibility. In all cases, equations of motion
were integrated with a time step of 1.0 fs.

Two force-matching potentials (FMP) [45, 46] were utilized. These have been
reported elsewhere to reproduce forces from density functional theory (DFT)
calculations performed on a-SiO2 and F-doped a-SiO2. The configurations used
for DFT calculations were extracted from trajectories of long-time classical
MD simulations performed using the Buckingham-type functional in the BKS
potential. The first FMP (named as FMP-v1 hereafter) reproduces density minima
of silica at a lower temperature than the BKS potential, but overestimates
density of silica at room temperature[45]. The second FMP (FMP-v2) was
developed by considering both forces of atoms and energies of configurations
to remedy overestimation of density in FMP-v1 [46]. Interestingly, FMP-v2 is
able reproduce the room temperature density of silica glass accurately, but
does not exhibit density anomaly of silica glass. This study aims to identify
some structural variations relating to the silica density anomaly by
quantifying topological differences of silica glass configurations obtained
using FMP-v1 and FMP-v2 force fields since the former reproduces said anomaly,
but the latter does not. Hereafter, pure silica glass models obtained by
FMP-v1 and FMP-v2 are named glass-v1 and glass-v2, respectively. All F-doped
silica glasses shown were modeled using only FMP-v1, and the effect of
F-doping on the reduction of the density anomaly was also investigated using
D-measure.

Table 1: Glass models for molecular dynamics simulations.  Abbreviation | No of atoms | F content | F/O
---|---|---|---
| Si | O | F | Total | [wt%] | ratio [%]
a-SiO2 | 3333 | 6666 | - | 9999 | - | -
SiO2-F1 | 3333 | 6600 | 132 | 10065 | 1.2 | 2.0
SiO2-F3 | 3333 | 6467 | 398 | 10298 | 3.7 | 6.2
SiO2-F5 | 3333 | 6333 | 666 | 10332 | 6.1 | 10.5
SiO2-F10 | 3333 | 6000 | 1332 | 10665 | 11.8 | 22.2

##  III Results and Discussion

###  III.1 Global graph of glass network

The densities of a-SiO2 and F-doped a-SiO2 simulated using FMP-v1 and FMP-v2
force fields are shown in Fig 1. The density of glass-v1 (a-SiO2 obtained with
FMP-v1) decreases with increasing temperature and reaches a minimum just above
2000 K. The density then increases again until it reaches a plateau around
3500 K. In general, the density of a material monotonically decreases with
increasing temperature due to positive thermal expansion coefficient. However,
silica, like water, exhibits negative thermal expansion coefficient at some
range of temperatures higher than its melting point. Contrary to glass-v1, the
glass-v2 (a-SiO2 obtained with FMP-v2) shows monotonic decrease in density
with increasing temperature even though the two force-fields have been
optimized based on similar DFT data [45, 46]. Hence, using D-measure, this
work tries to identify structural changes of silica networks that may relate
to the density anomaly phenomenon by comparing and analyzing topological
differences between these two glass models exhibiting different temperature
dependence of density.

Figure 1: Change of density as a function of temperature for pure and
fluorine-doped silica glass systems. Left axis indicates density for systems
generated from FMP-v1 (glass-v1 and SiO2-Fx) while right axis indicates
density for systems obtained from FMP-v2 (glass-v2). Glass-v1 and glass-v2
both correspond to pure amorphous silica glass systems.

To understand differences in the structures of glass-v1 and glass-v2, two
heatmaps of D-measure values are plotted in Figs 2a and 2b. The D-measure
heatmaps are obtained by first constructing graphs from configurations taken
every 50 K from the long-time MD trajectories and calculating the pairwise
dissimilarity between graphs from each system. Connectivity of all graphs
account for periodic boundary conditions. The bright color around the diagonal
suggests that similar temperatures have similar global structures, and for
both potentials a sharp transition is seen at 2300 K and 2700K respectively.
Fig 2c shows the temperature dependence of the D-values, where graphs
constructed at all temperatures from each system are compared against the
initial configuration graphs at 300 K. An interesting finding is that the
D-value monotonically increases with the temperature even for glass-v1,
implying that the FMP-v1-based system continuously varies the glass structure
from room temperature to high temperature without reversal to a more low-
temperature-like global structure. This suggests that high density structure
at higher temperatures is structurally dissimilar from the high density
structure at lower temperatures.

Figure 2: Heatmaps of pairwise D-measure values (D-values) when comparing
graphs obtained at different temperatures from silica glass systems (a) with
FMP-v1, (b) with FMP-v2. Temperature dependence of D-values for graphs of (c)
silica glass systems, and (d) F-doped silica glass systems, where the
reference graphs being compared to are the initial configuration graphs of the
respective systems. Points shown are averages of 5 independent MD trajectories
with error-bar indicating the standard deviation of the values.

Figs 2a and 2c show that the structural transition captured by D-measure of
glass-v1 begins to increase at a lower temperature than glass-v2. This is
consistent with the fact that fictive temperature of glass-v1, which was
determined as a folding point of temperature-potential energy curve, is lower
than that of glass-v2 (Table S1 in the Supporting Information). Therefore, it
is inferable that the glass structure changes identified by D-measure may be
related to the glass transition of silica glass. To confirm this hypothesis,
we evaluated the temperature dependence of D-values for the F-doped silica
glasses (Fig. 2d), since higher fluorine doping is known to lower the fictive
temperature of silica glass further[47, 48, 45]. This figure demonstrates that
the D-value of the silica glass containing more fluorine begins to increase at
lower temperature, in keeping with the experimental trends. According to these
observations, it is well supported that D-measure appropriately captures the
global network changes related to the glass transition in silica glass.

Next, we quantified the variability of D-values across five independent models
to ensure reproducibility of the structural trends described above. Fig 3
tracks averaged D-values among the same temperature frames of both glass-v1
and glass-v2 trajectories. Each point shown in the figure corresponds to mean
value of D-values between any 2 frames of the 5 replicas (10 unique pairs) of
either glass-v1 or glass-v2. Before reaching the fictive temperature, D-values
are very small, indicating that glass structures constructed from the
different initial configurations are almost identical from the view point of
graph networks. After crossing the fictive temperature, D-values begin to
increase and more diverse structures have begun to appear at higher
temperatures. This also demonstrates that D-measure is applicable to discuss
entropy change at high temperature by analyzing diversity of graphs
corresponding to structures of glass materials. The scale of dissimilarity
between replicas, even at higher temperature, is much lower than the scale of
changes described in Fig. 2. D-values representing dissimilarity between
glass-v1 and glass-v2, are also shown in Fig 3 (square markers). Here,
dissimilarity between 5 x 5 = 25 pairs of graphs were taken into account to
evaluate both the means and standard deviations. As can be seen from the
figure, configurations of the two glasses obtained from different force-fields
are more dissimilar than those modeled by the same force-field, indicating
that FMP-v1 and FMP-v2 build very different network morphologies.
Additionally, the small magnitude of the D-value standard deviation at low
temperatures confirms good reproducibility. Note that, contrary to our
expectation, dissimilarity between glass-v1 and glass-v2 decreases with an
increase of temperature to reach a minimum around 3000 K but interestingly
increases again at temperatures higher than 3000 K. This nonlinear behavior
seems like a reverse trend of temperature-density curve of the glass-v1;
however, more detailed analysis is required to unravel the reason behind the
non-monotonic change.

In summary, the comparisons across replicate simulations confirm that
D-measure consistently quantifies (dis)similarity of an entire glass network,
but the reason of density anomaly of glass-v1 was still unclear. We thus apply
D-measure to analyze local structures of the glass models in the next
subsection.

Figure 3:  Temperature dependences of D-measure values of five independent
silica glass systems obtained from FMP-v1 and FMP-v2 (glass-v1 and glass-v2)
against their respective initial configuration graphs. D-value measuring
dissimilarity between glass-v1 and glass-v2 are also constructed (25 pairs).

###  III.2 Local graph of glass network

To identify intermediate-range structures in the silica glasses, a variety of
templates were prepared by taking local structures from the crystalline silica
phases $\alpha$-quartz, tridymite, cristobalite, coesite, and stishovite (see
Fig 4). Fig 5 depicts the partial graphs obtained from these phases with a
cutoff distance ranging from 4 Å to 7 Å from a randomly chosen central silicon
atom. Similarly, local graphs were constructed from local structures of all
the silicon atoms in both glass-v1 and glass-v2 systems with the same cutoff
distances. We then evaluated the distributions of D-measure dissimilarity
between the crystal-phase local graphs and the glass local graphs to
characterize the glass structures. In addition, the temperature dependence of
these dissimilarities was examined in search for structural signatures density
anomaly of glass-v1. The cutoff distance to construct the local networks is
arbitrary, and the values of 5 Å and 7 Å were selected to target short-range
structures including the second coordination shell, and intermediate-range
structures considering higher coordinations respectively (Fig 5).

Figure 4: Unit cells of crystalline silica phases used as templates/baselines
for D-measure calculations of local graph structures of glass-v1 and glass-v2
systems. Numbers in the parentheses are approximate densities in unit
$g/cm^{3}$.  Figure 5: Sub-graphs extracted from the crystalline silica phases
(Fig 4) with cutoff distances ranging from 4 Å to 7 Å. Red spheres are the
central silicon atom, while yellow spheres show silicon atoms within the
cutoffs. Note that central silicon atoms are chosen randomly since choice of
central atoms would not yield different structures under periodicity
conditions.

The distributions of D-measure dissimilarity between local structures in our
silica glass systems and crystal phase graphs, within cutoff distances of 5 Å
and 7 Å are shown in Fig 6. Since it was shown in the previous section that
trajectories obtained using the same force-field exhibit good reproducibility,
the histograms were obtained using configurations from a single MD trajectory.
The D-value distribution with 5 Å cutoff distance highlights several clear
differences between the glass-v1 and glass-v2. In the case of glass-v1, the
distributions of structural dissimilarity with $\alpha$-quartz and coesite
exhibit peaks around $D=0.1$, implying that the local glass structures are
highly similar to these crystals. As the temperature increases, the height of
the first peak decreases and another broader peak appears at around $D=0.3$.
Emergence of significant bimodality in these distributions demonstrates
increasing structural disorder and suggests that it is localized to a fraction
of highly distorted environments rather than homogeneously distributed. On the
other hand, the first peak positions when compared to cristobalite and
stishovite occur around $D=0.2$, and even higher for tridymite (peak around
$D=0.4$). Increasing temperature shifts the distribution of dissimilarity to
stishovite to higher D-value, whereas the distribution against tridymite is
slightly shifted towards lower D-value. Contrary to glass-v1, D-value
distributions of glass-v2 exhibit relatively small temperature dependence and
higher multi modality. Interestingly, glass-v2 displays larger changes in
density with temperature than glass-v1 (Fig 1), while our results suggest that
the local structures within glass-v2 stayed relatively constant even at high
temperatures aside from small changes in interatomic distances responsible for
the density change.

Figure 6: Histograms of D-measure values when comparing local glass structures
and the crystalline-based templates with cutoff distances of (a) 5 Å and (b) 7
Å (Fig 5). The silica glass-v1 and glass-v2 models are constructed with FMP-v1
and FMP-v2, respectively. Histograms of different colors indicate local
environment of glass extracted from trajectories at different temperatures.

For local structures defined at a cutoff distance of 7 Å (Fig 6b), the D-value
distributions are narrower for both glass-v1 and glass-v2. Still, the
distributions of glass-v1 are broader than those of the glass-v2 at low
temperatures, suggesting that the glass-v1 possesses more diverse middle-range
structures in comparison with the glass-v2. Similar to short-range local
structures at 5 Å cutoff (Fig 6a), the distributions of dissimilarity between
glass-v1 and all crystalline phases shift to larger D-values at higher
temperatures, but less bimodality is observed. Similar observations can also
be made for glass-v2. However, different from short-range local structures,
distributions against $\alpha$-quartz and cristobalite crystals maintain only
one distinct peak each and the peak positions remains constant. These
observations imply that glass-v2 is more likely to maintain intermediate-range
structure even at high temperatures, which might be the cause of monotonic
decrease in density with temperature.

When the intermediate-range glass structure (7 Å cutoff) is compared with
coesite, the distribution exhibits a peak shift to lower D-value at higher
temperature. This trend is common in both glass-v1 and glass-v2, but glass-v1
exhibits a larger shift and eventually the distribution becomes bimodal at
2500 K. Interestingly, this may have indicated that more coesite-like
structures appeared as glass melt density increased. Since the coesite
structure is less symmetric, as seen in Fig 5, more disordered structures have
appeared and are expected to appear in higher temperatures. In addition to
that, such a local structure change might be related to the density anomaly of
glass-v1 since coesite has a higher density than vitreous silica [29] and
$\alpha$-quartz.

To evaluate reliability of D-measure in identifying topological differences,
Smooth Overlap of Atomic Positions (SOAP) descriptors were computed for the
same local glass structures used in Fig 6 and compared against SOAP of
crystalline-based silica templates[49]. Again, cutoff distances of 5 Å and 7 Å
were used and comparisons were done using cosine distance metric on the SOAP
descriptors. The results are shown in Fig 7. With SOAP, we can establish that
glass-v1 bears higher similarity to $\alpha$-quartz and coesite for local
structures defined at both 5 Å and 7 Å , while glass-v2 have local structures
which are more similar to tridymite and cristobalite. This observation aligns
with that seen with D-measure comparison (Fig 6). In addition, since the
tridymite structure is much less denser than the stishovite structure, greater
(dis)similarity to these structures reflects density of local glass
structures. In this case, the fact that glass-v1 is denser than glass-v2 is
well reflected in both SOAP and D-measure comparisons. However, there seems to
be no considerable change in SOAP of local structures with temperature, which
is rather different from what we observed in D-measure comparison. Based
solely on SOAP comparison, we can see that for local structures of glass-v1 at
2500 K with a cutoff distance of 5 Å, cosine distance distributions of SOAP
against coesite and stishovite have become slightly wider, indicating that the
local structures have become less similar to these crystalline phases. In
addition, peak of distribution against cristobalite at 2500 K has become
closer to 0.0. Furthermore, no multi modality can be observed in any
distribution. Contrasting this with D-measure comparison, we argue that
D-measure is better than SOAP at capturing slight changes in the topological
landscape.

Figure 7: Histograms showing cosine distances of SOAP between local glass
structures and crystalline-based templates of silica using cutoff distances of
(a) 5 Å and (b) 7 Å . Glass-v1 and glass-v2 are constructed with FMP-v1 and
FMP-v2 respectively. Histograms of different colors indicate local environment
of glass extracted from trajectories at different temperatures.

The degree of disorder of silica glass is usually measured by evaluating
tetrahedral order, $q$ (see SI), or distance between the fifth neighbor
silicon atom and the fourth neighbor oxygen atom from a silicon atom
(modified-d5) [37]. The temperature dependences of these order parameters for
the glass-v1 and glass-v2 are compared in Fig S1. Both glass-v1 and glass-v2
exhibit similar temperature dependence for both parameters, demonstrating
difficulty in finding structural origin of the density anomaly (or the lack
thereof) in glass-v1 and glass-v2 using the conventional short-range
structural parameters. However, it has to be noted that while D-measure
provides an effective way to identify structural changes and packing order in
amorphous materials, this graph-based method lacks capability to distinguish
bond angle changes directly, which could be important in many applications,
although it captures coordination number effectively.

###  III.3 F-doped silica glass

The magnitude of the density anomaly in silica glass decreases with increasing
of fluorine content. When more than 10 wt% of fluorine was added in our
simulations, the density anomaly could no longer be observed, as shown in Fig
1. Thereby, F-doping effect on the local structure of silica glass was
examined using similar local D-measure evaluation method shown in Section 3.2.
D-value distribution of local structures in F-doped silica glass against
coesite is shown in Fig 8. Here cutoff distances of 5, 6 and 7 Å are used. At
both 300 K and 2500 K, D-value distributions have become wider as more
fluorine is added. A possible explanation is that added fluorine atoms
terminate connections between silicon atoms by forming Si-F…Si instead of Si-
O-Si. Such formation increases the interatomic distances between silicon and
oxygen atoms, and due to the graph construction method used in this study,
terminated connections could be captured within the graphs, causing the local
structures to be more dissimilar to the crystalline phase. This trend is
consistently seen for all cutoff distances used.

A noticeable observation is that when compared to coesite, the D-value
distributions exhibit similar bimodality shown in the case of glass-v1 at 2500
K. However, peak heights at $D=0.1$ decrease as more fluorine atoms are added
to the system. It is thus expected that emergence of intermediate-range
structures similar to coesite at high temperatures is suppressed by F-doping,
which might be the reason why the density anomaly of silica glass is reduced
by the doping.

Figure 8: D-value distributions of local structures in pure and F-doped silica
glasses against the coesite phase at 300 K and 2500 K. Different cutoff
distances are used to show difference between short-range and intermediate-
range structures. Note that all F-doped silica glass systems are obtained from
FMP-v1.

##  IV Conclusions

This work proposes to apply D-measure, a metric to quantify graph
dissimilarities, as a novel approach to understand the global and local
structures of amorphous materials. This method involves casting silica glass
networks as graphs and to quantify the similarities between full all-atom
simulations or among intermediate-rage local environments. This approach we
applied to investigate the structural origin of the density anomaly in silica
glass. Atomistic structures from MD simulations with two silica glass models
were compared, only one of which exhibits a density anomaly

First, “global” graphs were defined by analyzing connectivity of silicon atoms
bridged by oxygen in the entire systems. Since global D-value was found to
monotonically vary with increasing temperature without showing any abnormal
behavior, it was concluded that the density anomaly is independent of the
global network change of silica glass. However, several advantages of the
D-measure for analyzing entire glass structure were found as follows:

  1. 1.

D-measure appropriately distinguishes similar glass structures encountered
during heating progression of the systems and identifies structural difference
of the silica glass systems modeled with different force-fields.

  2. 2.

The temperature dependence of D-value when comparing with the cooled glass
structure at 300 K resembles potential energy change, and the changing-point
relates to the fictive temperature.

  3. 3.

Increases in D-value when comparing different trajectories at temperatures
higher than the melting point are related to diversity of the melt structure
and could be related to the entropy of glass.

Next, using D-measure, the change in local structural with temperature and its
relation to the density anomaly was explored by comparing with local sub-
graphs of five crystalline silica phases. Accordingly, the silica glass system
that does not exhibit density anomaly was found to maintain a relatively
constant profile of short-range structures even at high temperatures. On the
other hand, our approach method was able to identify specific local structural
changes in the silica glass system showing density anomaly, where coesite-like
intermediate-range structures have begun to increase at high temperatures when
density of silica increases. The emergence of such lower-symmetry and high-
density environments in the intermediate-range may be the reason of density
anomaly in silica according to the graph-based structural analysis.

In summary, this study demonstrated the applicability of graph theory to find
structural signatures in an amorphous material. In general, global and local
D-measure comparisons are expected to help with better understanding of
structural changes within systems such that design of novel amorphous
materials could be accelerated.

##  V Competing Interests

The authors declare no competing interests.

##  VI Author Contributions

A.R.T. and S.U. contributed equally. A.R.T. and S.U. designed the experiments,
performed data analysis and wrote the manuscript. M.Y. wrote the code for
D-measure. R.G.-B. supervised the research, assisted in results interpretation
and contributed to manuscript writing.

##  VII Materials & Correspondence

Correspondence to Rafael Gómez-Bombarelli and Shingo Urata.

## References

## References

  * Rajan [2005] K. Rajan, Materials informatics, Materials Today 8, 38 (2005).
  * Ward and Wolverton [2017] L. Ward and C. Wolverton, Atomistic calculations and materials informatics: A review, Current Opinion in Solid State and Materials Science 21, 167 (2017).
  * Agrawal and Choudhary [2016] A. Agrawal and A. Choudhary, Perspective: Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science, Apl Materials 4, 053208 (2016).
  * Jablonka  _et al._ [2020] K. M. Jablonka, D. Ongari, S. M. Moosavi, and B. Smit, Big-data science in porous materials: materials genomics and machine learning, Chemical reviews 120, 8066 (2020).
  * Urata  _et al._ [2002] S. Urata, A. Takada, T. Uchimaru, A. K. Chandra, and A. Sekiya, Artificial neural network study for the estimation of the c–h bond dissociation enthalpies, Journal of fluorine chemistry 116, 163 (2002).
  * Urata  _et al._ [2003] S. Urata, A. Takada, T. Uchimaru, and A. K. Chandra, Rate constants estimation for the reaction of hydrofluorocarbons and hydrofluoroethers with oh radicals, Chemical physics letters 368, 215 (2003).
  * Curtarolo  _et al._ [2013] S. Curtarolo, G. L. Hart, M. B. Nardelli, N. Mingo, S. Sanvito, and O. Levy, The high-throughput highway to computational materials design, Nature materials 12, 191 (2013).
  * Schleder  _et al._ [2019] G. R. Schleder, A. C. Padilha, C. M. Acosta, M. Costa, and A. Fazzio, From dft to machine learning: recent approaches to materials science–a review, Journal of Physics: Materials 2, 032001 (2019).
  * Rickman  _et al._ [2019] J. M. Rickman, T. Lookman, and S. V. Kalinin, Materials informatics: From the atomic-level to the continuum, Acta Materialia 168, 473 (2019).
  * Morgan [1965] H. L. Morgan, The generation of a unique machine description for chemical structures-a technique developed at chemical abstracts service., Journal of Chemical Documentation 5, 107 (1965).
  * Rogers and Hahn [2010] D. Rogers and M. Hahn, Extended-connectivity fingerprints, Journal of chemical information and modeling 50, 742 (2010).
  * Dong  _et al._ [2015] J. Dong, D.-S. Cao, H.-Y. Miao, S. Liu, B.-C. Deng, Y.-H. Yun, N.-N. Wang, A.-P. Lu, W.-B. Zeng, and A. F. Chen, Chemdes: an integrated web-based platform for molecular descriptor and fingerprint computation, Journal of cheminformatics 7, 1 (2015).
  * Valle and Oganov [2010] M. Valle and A. R. Oganov, Crystal fingerprint space–a novel paradigm for studying crystal-structure sets, Acta Crystallographica Section A: Foundations of Crystallography 66, 507 (2010).
  * Montáns  _et al._ [2019] F. J. Montáns, F. Chinesta, R. Gómez-Bombarelli, and J. N. Kutz, Data-driven modeling and learning in science and engineering, Comptes Rendus Mécanique 347, 845 (2019), data-Based Engineering Science and Technology.
  * Gómez-Bombarelli and Aspuru-Guzik [2020] R. Gómez-Bombarelli and A. Aspuru-Guzik, Machine learning and big-data in computational chemistry, in  _Handbook of Materials Modeling: Methods: Theory and Modeling_, edited by W. Andreoni and S. Yip (Springer International Publishing, Cham, 2020) pp. 1939–1962.
  * Hiraoka  _et al._ [2016] Y. Hiraoka, T. Nakamura, A. Hirata, E. G. Escolar, K. Matsue, and Y. Nishiura, Hierarchical structures of amorphous solids characterized by persistent homology, Proceedings of the National Academy of Sciences 113, 7035 (2016).
  * Onodera  _et al._ [2019] Y. Onodera, Y. Takimoto, H. Hijiya, T. Taniguchi, S. Urata, S. Inaba, S. Fujita, I. Obayashi, Y. Hiraoka, and S. Kohara, Origin of the mixed alkali effect in silicate glass, NPG Asia Materials 11, 1 (2019).
  * Lee  _et al._ [2017] Y. Lee, S. D. Barthel, P. Dłotko, S. M. Moosavi, K. Hess, and B. Smit, Quantifying similarity of pore-geometry in nanoporous materials, Nature communications 8, 1 (2017).
  * Shell  _et al._ [2002] M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Molecular structural order and anomalies in liquid silica, Physical Review E 66, 011202 (2002).
  * Urata and Li [2018] S. Urata and S. Li, A multiscale shear-transformation-zone (stz) model and simulation of plasticity in amorphous solids, Acta Materialia 155, 153 (2018).
  * Takada  _et al._ [2009] A. Takada, P. Richet, and T. Atake, New description of structural disorder in silica glass, Journal of non-crystalline solids 355, 694 (2009).
  * Takada  _et al._ [2010] A. Takada, P. Richet, and T. Atake, New description of structural disorder in silica glass. ii. application of atomistic energy distribution analysis to pressure effects, Journal of non-crystalline solids 356, 2486 (2010).
  * Schieber  _et al._ [2017] T. A. Schieber, L. Carpi, A. Díaz-Guilera, P. M. Pardalos, C. Masoller, and M. G. Ravetti, Quantification of network structural dissimilarities, Nature communications 8, 1 (2017).
  * Kolomeets  _et al._ [2019] M. Kolomeets, A. Chechulin, and I. V. Kotenko, Social networks analysis by graph algorithms on the example of the vkontakte social network., J. Wirel. Mob. Networks Ubiquitous Comput. Dependable Appl. 10, 55 (2019).
  * Horawalavithana and Iamnitchi [2019] S. Horawalavithana and A. Iamnitchi, On the privacy of dk-random graphs, arXiv preprint arXiv:1907.01695 (2019).
  * Jiang  _et al._ [2017] H. Jiang, Q. Li, X. Zhou, Y. Chen, S. Yi, H. Wang, and Z. Lu, A collective human mobility analysis method based on data usage detail records, International Journal of Geographical Information Science 31, 2359 (2017).
  * Mheich  _et al._ [2020] A. Mheich, F. Wendling, and M. Hassan, Brain network similarity: methods and applications, Network Neuroscience 4, 507 (2020).
  * Schwalbe-Koda  _et al._ [2019] D. Schwalbe-Koda, Z. Jensen, E. Olivetti, and R. Gómez-Bombarelli, Graph similarity drives zeolite diffusionless transformations and intergrowth, Nature materials 18, 1177 (2019).
  * Brueckner [1970] R. Brueckner, Properties and structure of vitreous silica. i, Journal of non-crystalline solids 5, 123 (1970).
  * Huang  _et al._ [2004] L. Huang, L. Duffrene, and J. Kieffer, Structural transitions in silica glass: thermo-mechanical anomalies and polyamorphism, Journal of non-crystalline solids 349, 1 (2004).
  * Geske  _et al._ [2016] J. Geske, B. Drossel, and M. Vogel, Fragile-to-strong transition in liquid silica, AIP Advances 6, 035131 (2016).
  * Vollmayr  _et al._ [1996] K. Vollmayr, W. Kob, and K. Binder, Cooling-rate effects in amorphous silica: A computer-simulation study, Physical Review B 54, 15808 (1996).
  * Horbach and Kob [1999] J. Horbach and W. Kob, Static and dynamic properties of a viscous silica melt, Physical Review B 60, 3169 (1999).
  * De Oliveira  _et al._ [2009] A. B. De Oliveira, P. A. Netz, and M. C. Barbosa, An ubiquitous mechanism for water-like anomalies, EPL (Europhysics Letters) 85, 36001 (2009).
  * Soules  _et al._ [2011] T. F. Soules, G. H. Gilmer, M. J. Matthews, J. S. Stolken, and M. D. Feit, Silica molecular dynamic force fields—a practical assessment, Journal of non-crystalline solids 357, 1564 (2011).
  * Jabes  _et al._ [2012] B. S. Jabes, D. Nayar, D. Dhabal, V. Molinero, and C. Chakravarty, Water and other tetrahedral liquids: order, anomalies and solvation, Journal of Physics: Condensed Matter 24, 284116 (2012).
  * Shi and Tanaka [2018] R. Shi and H. Tanaka, Impact of local symmetry breaking on the physical properties of tetrahedral liquids, Proceedings of the National Academy of Sciences 115, 1980 (2018).
  * Van Beest  _et al._ [1990] B. Van Beest, G. J. Kramer, and R. Van Santen, Force fields for silicas and aluminophosphates based on ab initio calculations, Physical review letters 64, 1955 (1990).
  * Tsuneyuki  _et al._ [1988] S. Tsuneyuki, M. Tsukada, H. Aoki, and Y. Matsui, First-principles interatomic potential of silica applied to molecular dynamics, Phys. Rev. Lett. 61, 869 (1988).
  * Bertolazzo  _et al._ [2016] A. A. Bertolazzo, A. Kumar, C. Chakravarty, and V. Molinero, Water-like anomalies and phase behavior of a pair potential that stabilizes diamond, The Journal of Physical Chemistry B 120, 1649 (2016).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computational physics 117, 1 (1995).
  * Nosé [1984] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, The Journal of chemical physics 81, 511 (1984).
  * Tuckerman  _et al._ [2006] M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim, and G. J. Martyna, A liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble, Journal of Physics A: Mathematical and General 39, 5629 (2006).
  * Urata and Sato [2017] S. Urata and Y. Sato, A study on the plasticity of soda-lime silica glass via molecular dynamics simulations, The Journal of chemical physics 147, 174501 (2017).
  * Urata  _et al._ [2021a] S. Urata, N. Nakamura, K. Aiba, T. Tada, and H. Hosono, How fluorine minimizes density fluctuations of silica glass: Molecular dynamics study with machine-learning assisted force-matching potential, Materials & Design 197, 109210 (2021a).
  * Urata  _et al._ [2021b] S. Urata, N. Nakamura, T. Tada, and H. Hosono, Molecular dynamics study on the co-doping effect of al2o3 and fluorine to reduce rayleigh scattering of silica glass, Journal of the American Ceramic Society (2021b).
  * Shiraki and Ohashi [1992] K. Shiraki and M. Ohashi, Scattering property of fluorine-doped silica glasses, Electronics letters 28, 1565 (1992).
  * Kirchhof  _et al._ [2018] J. Kirchhof, S. Unger, and J. Dellith, Viscosity of fluorine-doped silica glasses, Optical Materials Express 8, 2559 (2018).
  * Bartók  _et al._ [2013] A. P. Bartók, R. Kondor, and G. Csányi, On representing chemical environments, Physical Review B 87, 10.1103/physrevb.87.184115 (2013).

See pages 1 of SupportingInformation.pdf See pages 2 of
SupportingInformation.pdf See pages 3 of SupportingInformation.pdf