# Prediction of Carbon Nanostructure Mechanical Properties and Role of Defects
Using Machine Learning

Qi Zhao  Department of Computer Science and Engineering, University of
California San Diego, La Jolla, CA 92093  Jordan J. Winetrout  Department of
Chemical and Biological Engineering, University of Colorado Boulder, Boulder,
CO 80309, and Materials Science and Engineering Program, University of
Colorado Boulder, Boulder, CO 80309  Yanxun Xu  Department of Applied
Mathematics and Statistics, Johns Hopkins University, Baltimore, MD 21218
Yusu Wang<EMAIL_ADDRESS>Department of Computer Science and Engineering,
University of California San Diego, La Jolla, CA 92093  Halıcıoğlu Data
Science Institute, University of California San Diego, La Jolla, CA 92093
Hendrik Heinz<EMAIL_ADDRESS>Department of Chemical and Biological
Engineering, University of Colorado Boulder, Boulder, CO 80309, and Materials
Science and Engineering Program, University of Colorado Boulder, Boulder, CO
80309

(September 2021)

###### Abstract

Carbon fiber and graphene-based nanostructures such as carbon nanotubes (CNTs)
and defective structures have extraordinary potential as strong and
lightweight materials. A longstanding bottleneck has been lack of
understanding and implementation of atomic-scale engineering to harness the
theoretical limits of modulus and tensile strength, of which only a fraction
is routinely reached today. Here we demonstrate accurate and fast predictions
of mechanical properties for CNTs and arbitrary 3D graphitic assemblies based
on a training set of over 1000 stress-strain curves from cutting-edge reactive
MD simulation and machine learning (ML). Several ML methods are compared and
show that our newly proposed hierarchically structured graph neural networks
with spatial information (HS-GNNs) achieve predictions in modulus and strength
for any 3D nanostructure with only 5-10% error across a wide range of possible
values. The reliability is sufficient for practical applications and a great
improvement over off-the shelf ML methods with up to 50% deviation, as well as
over earlier models for specific chemistry with 20% deviation. The algorithms
allow more than 10 times faster mechanical property predictions than
traditional molecular dynamics simulations, the identification of the role of
defects and random 3D morphology, and high-throughput screening of 3D
structures for enhanced mechanical properties. The algorithms can potentially
be scaled to morphologies up to 100 nm in size, expanded for chemically
similar compounds, and trained to predict a broader range of properties.

##  1 Introduction

Carbon fiber is among the most promising engineering materials for the 21st
century due to superior tensile strength and modulus relative to low weight
(Figure 1).[1] It is extensively used in the automotive, aviation, and
aerospace industry (Figure 1a). Carbon fiber is composed of nanostructured
carbon such as carbon nanotubes (CNTs), graphitic layers, and sometimes
polymer binders, for example, in high-strength carbon fiber yarns. In addition
to structural applications and fabricating ultra-strong composites, CNTs and
graphitic materials are also used in catalyst supports, ultra-small electrical
and thermal conductors, as well as membranes for water desalination that
benefit from controlled microstructure, conductivity, and mechanical
stability.[2][3]

One of the major challenges is closing the gap between the current mechanical
properties, which are already superior to most other materials, and the
theoretically possible performance, which is significantly higher (Figure
1b).[4] Specifically, the theoretically possible Young’s Modulus is 1 TPa
while current IM-7 fiber achieves 400 GPa, and the theoretically possible
tensile strength amounts to 100 GPa while current IM-7 fiber barely reaches 8
GPa.[1] There is a space of improvements by multiples to explore, which can
facilitate revolutionary changes in multiple industries including automotive,
air and space flight.

Since the discovery of graphene and CNTs, there has been a multidisciplinary
effort to understand and utilize their unique and remarkable properties. Loss
in mechanical performance, in particular, is related to defects, impurities,
random orientations as well as finite lengths and resulting gaps between
CNTs.[4] It remains difficult to quantify the impact of these features on
mechanical properties and tackle related synthesis and processing challenges
towards customized manufacturing from the atomic scale to the microscopic
scale. Characterization of the structure in experiments relies on X-ray
scattering, microscopy and 3D tomographic reconstruction, as well as indirect
methods such as Raman spectroscopy.[5] Current measurements remain tedious and
expensive, limited in the number of samples and in revealing failure
mechanisms down to the atomic and nanometer-scale. To complement such efforts,
simulation methods have been employed such as molecular dynamics simulations
and density functional theory. The advantages of these techniques include
insight at the scale of atoms, access to the large nanometer scale, and the
ability to simulate entire stress-strain experiments from equilibrium to
failure in high accuracy.[6][7] Simulations then allow to design and screen a
larger number of model structures with a variety of defects and nanoscale
features of interest, inspired by experimental data and theory, and examine
the failure mechanisms in depth. However, although more efficient than
experiments, the computational cost is still considerable.

Recently, there has been a major trend to explore properties of materials
using rapidly developing techniques from computer science and machine learning
(Figure 1c, d).[8] Machine learning has already achieved remarkable success on
a wide range of applications, such as image classification and segmentation,
language translation, and DNA sequence analysis [9, 10, 11, 12, 13, 14, 15].
Different from conventional approaches, machine learning (ML) algorithms build
models and infer knowledge about materials from sampled data called training
data. The trained models can then be used to predict materials properties or
changes in the dynamics of fictive materials structures. Researchers first
focused on using vectorized global representations of molecular structures to
predict their properties. Besides some simple features describing basic
properties of molecular structures like the number of atoms or bonds, a list
of molecule representation methods including Hamiltonians [16], Coulomb matrix
[17], Bag of Bonds (BoB) [18], electronic density [19], symmetry functions
[20], and fingerprints [21] have been suggested. Kernel methods and Gaussian
Models are commonly applied with those global representations in downstream
machine learning tasks [22, 17, 18]. Kernel ridge regression and Gaussian
process regression on such representations can be used to predict cohesive
energies of molecular structures and other properties [23]. Although kernel
methods are powerful when dealing with relatively small datasets of vectorized
representations, it still has two main drawbacks. First, when the size of
training data set grows increasingly large, the computation and storage kernel
methods become extremely expensive. For example, the computation of a millions
dimensional kernel matrix is a disaster. Second, as these vectorized
representations cannot extract all major information from molecules, hidden
features that are potentially important will be lost.

A more powerful machine learning framework in the age of data is deep learning
(DL) [24, 25]. Roughly speaking, deep learning uses multiple neural network
layers to progressively learn representations, and is often used to model the
map from input raw data to output target properties. Neural networks can learn
molecules representations and apply them to prediction tasks in a more
intelligent and data-driven way. Depending on how they represent molecular
structures, we can classify deep learning approaches as follows. First, SMILES
[26] is a traditional text-processing architecture that represents molecules
as sequences, on which it has been common for researchers to apply Recurrent
Neural Networks (RNN) or Long Short-term Memory (LSTM) techniques [27, 28,
29]. Second, to better encode the 3D structural information of molecules, one
could model a molecular structure as a 3D point cloud where each atom is
represented by a point in $\mathbb{R}^{3}$, and then apply the PointNet [30]
or RS-CNN [31] architecture to perform property prediction over such point
cloud representation [32]. However, simply representing molecules as point
clouds can lose significant structural information such as chemical bonds.
[33] and [34] then improved this approach by designing more expressive point
cloud convolution networks which align local neighbor information and symmetry
factorization of point distributions. Furthermore, fingerprints [21, 35, 36]
have the power to encode the chemical bond relations which can be viewed as
the start of applications of Graph Neural Networks for the prediction of
properties of molecular structures.

Finally, a more natural way to organize molecules is to represent them as
graphs (Figure 1c). Atoms are represented as nodes in graphs and two atoms are
connected when there exists a chemical bond between them. Graph Neural
Networks (GNNs) like Graph Convolution Networks (GCN) [37] and Graph Attention
Networks (GAT) [38] can perform convolution over graph-structured data. MPNN
[39] has built a unified framework on designing GNNs as well as a specific
architecture for predicting molecular properties. DimeNet [40] and D-MPNN [41]
further developed this approach by incorporating bonded interactions into
models. Although GNNs outperform other approaches mentioned above in tasks
like classification of small molecules, there are limitations in existing
approaches: (1) In addition to atom connectivity and bonded interactions, a
good model should also take non-bonded relations and features characterizing
spatial (geometric) patterns of atomic structures into consideration. (2) Due
to over-smoothing or over-squashing phenomena, long-range interactions /
signals are hard to capture, making it challenging to process large scale
molecular structures and assemblies.

##### Contribution of this paper.

In this paper, we develop a new machine learning framework to predict
mechanical properties for carbon nanostructures directly from their atomic
structure configurations (Figure 1d). Our ML framework, called HS-GNN, is an
enhanced graph neural network architecture, carefully designed to process CNT
atomic structures more effectively. In particular, our HS-GNN provides a much
more powerful way to encode the input microstructure configuration than
previous approaches: (i) The input is modeled as a heterogenuous graph so as
to capture both the covalent bonds and the short-range (non-bonded)
interactions among atoms. (ii) The geometric shape of the input structure
(e.g, how the CNT sheet is curved, and what geometric patterns it forms with
atoms from neighboring sheets) can affect the mechanical properties of the
structure. We leverage local geometric and topological features to better
encode the spatial geometric shape of the input structure. (iii) We deploy a
hierarchical neural network design to capture large-scale interaction among
different parts of CNT composite.

We train our HS-GNN on initial configurations of about 1000 independent CNTs.
Tensile properties used as ground truth in training set are generated from
their stress-strain curves of complex carbon nanostructures using molecular
dynamics simulations with the thoroughly validated IFF force field that
includes bond breaking (IFF-R) [42, 43, 7]. The resulting HS-GNN can predict
Young’s modulus and strength for new atomic arrangements (not used in training
set) with high accuracy (5-10% MSE error), including larger graphitic
structures. The time it takes is less than 10% of the computational cost of
traditional molecular simulations. We also compare our HS-GNN with six other
ML approaches including the state-of-the art DimeNet [40], and our new
approach achieves significantly better accuracy than these other ML
approaches. The ML source code and examples, as well as key data and run
scripts for the training set, will be shared publicly.

Figure 1: Hierarchical structure of carbon-based materials, challenges and new
machine learning workflow to predict mechanical properties at the nanoscale.
(a) Automotive, air and space vehicles (top) increasingly require strong and
lightweight materials, for which carbon fiber and related composites are a
promising choice. At the microscale, carbon fiber sheet and yarn consist of
aligned yet highly defective carbon fiber bundles, which can be composed of
carbon nanotubes and irregular graphitic structures (bottom left). 3D imaging
by tomography is typically limited to the coarse nanometer scale. In this
work, we examine stress-strain relationships up to failure for a library of
over 1000 structures at the small nanometer scale and feed the data into a ML
model (bottom right). (b) The current performance of carbon fiber is much
lower than the theoretical limits, which are ca. 1000 GPa in tensile modulus
and ca. 100 GPa in tensile strength. While efforts over the last 50 years have
brought steady improvements, the bottleneck is understanding and engineering
of defects in the materials that originate from the molecular scale to the
nanoscale. (c) We introduce a hierarchical spatial graph neural network (HS-
GNN) that shows excellent performance in mechanical property predictions for a
wide range of graphitic morphologies up to the large nanometer scale. (d) The
workflow includes data collection by reactive MD simulations (validated
relative to experiment), representation of the structures by the HS-GNN,
training of the ML model with the data for modulus and strength, and fast
screening (prediction) of mechanical properties of unknown structures to
quantify the role of defects and desirable new nanoscale designs.

##  2 Results and Discussion

###  2.1 Molecular structures, training data, feature definition, and causal
relationships.

The benchmark dataset used to train and test the ML models was obtained from
precise reactive molecular dynamics simulations of 1159 structures of carbon
nanotubes, graphite, defective and deformed graphitic assemblies, utilizing
available experimental data and images to the extent available (Figure
2).[5][44] The simulations used the reactive interface force field (IFF-R),
which quantitatively reproduces pi-pi stacking, surface and interfacial
energies, Young’s moduli and tensile strength of CNTs, graphene, and graphite
in agreement with experiments.[43][7] CNT morphologies in the training set
included individual and bundled single-wall (Figure 2a), double-wall (Figure
2b inner wall highlighted in blue), and triple-wall CNTs (Figure 2c inner most
wall highlighted in yellow). The nanotubes simulations were modeled with
experimentally determined structural and defect considerations such as
pristine (Figure 2a,d), with missing atoms (Figure 2b-d), with discontinuities
or fracture defects (Figure 2d,f), and with reconfiguration or Stone-Wales
defects (Figure 2d). Global features such as the number of nanotubes in a
bundle, the type of nanotube modeled, the diameter of the nanotubes, the
number of atoms and bonds in the simulation, and the cross-sectional area of
the simulated bundle were identified as training features that could
adequately capture the structural considerations that may influence CNT
mechanical properties (Figure 2e). Tensile strength and tensile modulus are
important parameters for designing safe materials for structural applications
because they represent critical loads before material failure. The
significance of strength and modulus for materials’ design motivated the
development of a machine learning model capable of predicting these two
quantities for carbon nanotube bundles. The stress-strain curves of all
structures were simulated up to failure to extract the tensile modulus and
tensile strength (Figure 2f). Tensile strength is defined as the maximum
stress before material failure. Tensile modulus is a metric that defines how
hard a material is, and is calculated as the slope of the linear stress-strain
response (typically 0.0-0.01 strain for hard materials). For a double-wall
carbon nanotube (DWCNT) with fractured inner walls the tensile strength and
tensile modulus were determined to be 19.2 GPa and 210 GPa respectively.

Figure 2: Overview of the molecular dynamics data set and methods that were
used to train the hierarchical spatial graph neural network (HS-GNN). Top and
bottom view of (a) single-walled carbon nanotube with no defects, (b) bundle
of double-walled carbon nanotubes with mono-vacancy defects, (c) dundle of
triple-walled carbon nanotubes with mono-vacancy, di-vacancy, and tri-vacancy
defects. (d) Overview of the structural features considered in the data set
including defectless (pristine), discontinuous nanotubes, reconfiguration
defects (Stone-Wales), and missing atoms (mono-, di-, and tri-vacancy).(e)
Top-down view of a bundle of 3 DWCNT demonstrating the global features
considered in the data set including the number of nanotubes in a bundle, the
type (double-wall) of nanotube simulated, the diameter of each nanotube,
number of atoms in the simulation, and cross-sectional area of the nanotube
bundle. (f) Stress-strain curve from a tensile simulation of two double-walled
carbon nanotubes with snapshots of the nanotubes at different strain points in
the tensile simulation. Highlighted in red are most of the mechanical
properties used for training HS-GNN (strength, modulus, and strain at break).
HS-GNN was used to predict strength and modulus, so the values calculated from
IFF-R are listed next to their labels in parentheses. (f) Relationship between
some simulation features and their mechanical performance (strength and
modulus highlighted in red), as determined by a causal additive model (CAM).
The relationship determined by the CAM was used to inform the training of HS-
GNN in which the prediction of strength was used to predict modulus.

To rank the importance of each structural features’ influence on CNT
mechanical performance, structural feature significance for strength and
modulus predictions were evaluated using a causal additive model (CAM), which
aims to infer causal relations among features. [45] The CAM estimates strength
as a direct cause of modulus (Figure 2g). The causal discovery information was
added to our feature list and improvements by adding this relations into the
ML predictions of modulus and strength will be discussed later.

###  2.2 Machine learning pipeline for property prediction

The high level framework of our machine learning pipeline, called HS-GNN
(Hierarchical Spatial Graph Neural Networks), for mechanical property
predictions is shown in Figure 3. We describe the key ideas below. More
details can be found in Section 4.6.

The input to HS-GNN is the initial atomic structure S of a CNT configuration.
Instead of using a PointNet-type [30] architecture to take the set of atoms,
viewed as a set of points each equipped with a radius, as input, we recognize
the importance of both covalent bonds in the structure, as well as short-range
interactions among atoms. Therefore we use a specially designed graph neural
network-based architecture to process S which we denote by _HS-GNN_. In
particular, novel features of our HS-GNN architecture include: a heterogeneous
graph representation of input, a hierarchical neural network model so as to
capture large-scale interaction among different parts of the CNT bundle, and
the injection of local geometric and topological features to better encode the
shape of CNT bundle into the neural network via attention mechanism. We
briefly describe these key components below.

Heterogeneous graph:

It is common to use a _bond-graph_ $G_{bond{\texttt{S}}}$ to represent S,
where each node corresponds to an atom, and there is an edge between two nodes
if the corresponding atoms form a covalent bond. We go beyond the bond-graph
and generate a heterogeneous graph $G_{\texttt{S}}$ with multiple types of
edges among nodes from S. In particular, we add edges between nodes whose
corresponding atoms are spatially close (i.e, atoms from different nanotubes
but within 6Å distance) or have small effective resistence distance. We also
connect nodes whose corresponding atoms form a dihedral angle, the angle of
two planes formed by four sequentially bonded atoms rotated about a central
bond. This heterogeneous graph is then fed to a hierarchical GCN.

Hierarchical GCN model:

A GCN (graph convolutional network) is a popular type of GNN (graph neural
network; see Supplementary Material for a brief introduction). On the high
level, a standard GCN takes a graph with initial node features as input. The
input graph then goes through $\ell$ convolutional layers, during each of
which it will update the feature vector stored at each graph node by
aggregating features from neighboring nodes and performing a learned
transformation of node features. See Figure 3 (a, c) for an illustration.
However, it is known that GNNs tend to have the oversmoothing issue [46] and
often cannot go very deep, thus limiting the aggregation of information from
long range interactions. To address this issue, we use a hierarchical GCN that
process the input graph $G_{\texttt{S}}$ at multiple (three in our current
experiments) resolutions: see Figure 3 (a, b). Different from hierarchical GNN
approaches proposed by [47, 48] etc., HS-GNN finds hierarchies by spatial
geometry information instead of graph topology or nodes features. The lowest
level $L_{1}$ operates on input graph $G_{\texttt{S}}$; while in a higher
level $L_{i}$, each node corresponds to a cluster of nodes of level $L_{i-1}$
and we can call such a node a _super-node_. The feature vector associated to
each super-node is obtained by a pooling layer on level $L_{i-1}$ graphs.
Within each resolution level $L_{i}$, we perform several GCN aggregation
layers, which we refer to as micro-layers. In the end, a last pooling layer is
applied so as to obtain graph-level prediction.

Encoding local geometry/topology:

The shape of nanotubes and the spatial relation among neighboring tubes/sheets
impact the final property. To encode local geometry, we use Principle
Component Analysis (PCA) to capture a “curvature"-like quantity for each
resolution level $L_{i}(i>1)$. To encode local interactions among spatially
close atoms (from potentially different nanotubes), we use the so-called
persistent homology (PH) to characterize the spatial distribution of points
within each cluster in each level. Persistent homology is one of the most
important developments in the field of topological data analysis in the past
two decades and have already been applied to characterize different types of
complex shapes [49, 50]. In our case, given a set of points (atom centers)
forming a local cluster, by tracking the creation and death of topological
features w.r.t. a growing sequence of space around these points, it can
provide a meaningful yet succinct summary of the shape formed by these points.
See Figure 3 (d) for a simple illustration, and see Section Persistent
Homology in the Supplement for a brief introduction.

Node features and edge attention:

The above information are incorporated into our HS-GNN via node features and
via the _edge attention mechanisms._ Given any graph node $v_{a}$
correspoinding to atom $a$, our initial node feature vector $\mu^{(0)}(v_{a})$
in level $L_{1}$ includes the degree of $v_{a}$ in the bond graph (measuring
local “defect"), the 3D coordinates of atom $a$, as well as random generated
features. For the higher level GNNs, the persistence summary of a cluster is
included in the features of corresponding super-node. Node features will be
transformed and aggregated through the neural networks. The edge-attention
mechanism intuitively allows one to compute “differential" of information at
two endpoints $v_{a},v_{b}$ of an edge ($v_{a},v_{b}$) ($v_{a}$ and $v_{b}$
will be super-nodes/clusters in $L_{j}(j>1)$ level HS-GNN), and use this to
weight this connection (edge) $(v_{a},v_{b})$ when aggregating information for
node $v_{a}$ from its neighbors. We use the two major principal vectors and
the norm vector in local Principal Component Analysis (PCA) of points in each
cluster for the edge-attention in $L_{j}$ level HS-GNN.

Figure 3: The machine learning pipeline is a carefully designed hierarchical
spatial GNN. (a) Overview of the pipeline to predict mechanical properties of
CNTs and carbon nanostructures. Line 1 illustrates the generation of a
heterogeneous graph for a given CNT structure and processing it with a
hierarchical GNN. In each level of the hierarchy, we create a specific GNN to
learn the graph representations and the list of vertical arrays within each
level indicates the micro-layers. In addition, we compute topological and
geometric features, more specifically, persistence summaries and PCA
information of clusters from CNT structure. (b) Constructing the hierarchical
graph series: In order to obtain a higher level graph, we group nodes in the
lower level graph that are spatially close into clusters. Each node in the
higher level graph represents a node cluster in the lower level, and we
connect nodes if their corresponding clusters share overlap. (c) Illustration
of the message-passing workflow of the GNN in Level-2 and Level-3. Each node
passes its features reweighted by geometric features of the corresponding
cluster relative to the neighborhood.

Figure 3 illustrates the entire machine learning pipeline. Message passing is
the scheme applied in GNN, while persistence summaries and geometric features
are used in message passing process.

###  2.3 Results on Modulus and Strength prediction

Using the data set of 1159 tensile simulations of CNT bundles and graphitic
structures, we trained HS-GNN and evaluated the predictions of mechanical
properties (Figure 4, 5).

##### Tensile properties distribution.

The distribution of tensile strengths within the training set is between 0-120
GPa, and the distribution of tensile moduli is between 0-1000 GPa. Few CNT
bundles had strength and modulus values at the upper and lower limits
attributing to strictly pristine or completely fractured cases respectively.
However, even in highly crystalline CNT bundles it is impossible to achieve
defect-less CNTs. To account for this, the training set is predominately
comprised of CNT bundles that are flawed through missing atoms, Stone-Wales
defects, discontinuities, or a combination of defects. The training set
distribution reflects this, demonstrating a large portion of tensile strengths
between 40-100 GPa and tensile moduli between 300-700 GPa (Figure 4 e,f). To
achieve values of 0 GPa for tensile strength and modulus, some structure files
generated after the tensile simulation were used as training structures.
Completely fractured CNT bundles are not expected to be the primary CNT bundle
morphology of interest, so a small population of completely fractured
structures was included in the training set. 115 initial 3D computer generated
CNT structures were used for testing of HS-GNN. Test structures were randomly
sampled to ensure that complex structural features such as straight and
compressed nanotubes, pristine and defective structures, single and multi-wall
nanotubes, etc. are captured in the test set.

##### Results.

The HS-GNN method has shown to predict the tensile properties with a relative
mean squared error (MSE) of 4.1% deviation from the IFF-R simulations for
strength, and 7.6% deviation for tensile modulus when the strength prediction
and pretrained HS-GNN model (see Section 4.7 for details) is used to predict
modulus (HS-GNN-C) as determined by the CAM (Table 1). A 0.6% MSE improvement
is observed for predicting tensile modulus when using the strength prediction
to predict modulus (HS-GNN-B) is compared to HS-GNN without using any property
predictions as a input feature (HS-GNN-A). In a more detailed analysis,
tensile strength and modulus prediction error is dependent on tensile
properties distribution in training set. (Figure 4). Tensile strength (Figure
4a) and tensile modulus (Figure 4b) prediction error was greater than or equal
to a 20 GPa difference between the prediction and simulated values only for
strength populations less than 3.7% of the training set, for example, cases in
which strength was simulated to be 0-40 GPa, or modulus between 0-200 GPa. The
number of those cases in training set is relatively smaller as shown in Figure
4 e, f. Despite having a smaller portion in the training set, the predictions
for CNT bundle cases close to the upper limit of performance (100-120 GPa and
800-1000 GPa for strength and modulus respectively) were below a 20 GPa
difference between the prediction and simulated values, albeit with a larger
prediction error deviation compared to the predominate sample range. The
accurate predictions for a smaller portion in training set may be accredited
to the less complex morphology of pristine CNT bundles that closely resemble
somewhat flawed cases (Figure 2). On the other hand, fractured CNT bundles (0
GPa strength and modulus cases) can have dramatically different structures
after material failure when compared to their unfractured counterparts (Figure
2f) leading to a higher property prediction error (Figure 4 a-d). What’s more,
HS-GNN takes a server with 24 CPU cores and 1 RTX A6000 GPU about 7 hours to
compute the persistence summaries and PCA information of CNT bundles in
training set and train the model, but the inference process only takes about 2
minutes to obtain the predicted tensile strength and modulus for a CNT case in
test set while IFF-R simulation takes 24 CPU cores about 30 minutes to compute
the tensile properties.

Figure 4: Prediction of mechanical properties. (a) Box plot of prediction
results of HS-GNN on tensile strength. We divide the CNT-based models into 6
groups based on their tensile strength from 0 to 120 GPa. (b) Box plot of
prediction results of HS-GNN on tensile modulus. We divide the CNT-based
models into 5 groups based on their tensile modulus from 0 to 1000 GPa. The
Box plots show the average and standard deviation inside the box, as well as
the maximum and minimum outside the box. (c, d) Comparison between HS-GNN and
6 baselines on prediction over tensile strength and tensile modulus. The
baselines are linear regression based on global features of CNT materials,
PointNet and RS-CNN over point cloud formation CNT structures, GNN over bond
graphs without the help of local topological and geometric features, MPNN and
dimeNet taking simple geometric features like bond angles. (e,f) The
distribution of tensile strengths and tensile moduli in the training dataset.

##### Comparison with other ML approaches.

We compare the performance of HS-GNN to 6 baseline machine learning models,
linear regression (LR) based on global features (i.e., number of nanotubes,
number of atoms, diameter of a nanotube) of CNT bundles, PointNet [30] and RS-
CNN [31] on point cloud formation CNT bundles, and three GNN models, vanilla-
GNN [37], MPNN [39] and DimeNet [40]. Among them, DimeNet is considered to be
the state-of-the-art in molecular structures properties prediction tasks. HS-
GNN outperforms all 6 baselines in predicting tensile strength and modulus of
CNT bundles (Figure 4 c,d and Table 1). In general, all models follow the same
prediction error trend by which a lower prediction error is observed for cases
in the property ranges of more training cases (i.e. between 40-100 GPa for
strength and 400-800 GPa for modulus). However, only while using HS-GNN is a
large upswing in prediction error not observed for the cases within the upper
limit of strength (Figure 4c). For nearly all methods of training HS-GNN, the
mechanical property prediction accuracy for all non-large cases was 2-5 times
better than baselines.

##### Generalization to large structures.

We also test our HS-GNN on 16 larger (>20,000 atoms) CNT structures who
consist of much more atoms than those in the training set, so as to see how
well our trained model generalize as the size of input scales up. HS-GNN-A
method predicts strength with a MSE deviation of 12.7% and modulus with a MSE
deviation of 17.6% from IFF-R simulations (Table 1). A 1.6% MSE improvement is
shown for HS-GNN-B and 2.3% MSE improvement for HS-GNN-C, exemplifying the
benefit of using property predictions as input features. For generalization of
the trained ML models, HS-GNN was nearly 3 times more accurate than vanilla-
GNN or LR, and about 1.3 times more accurate than DimeNet (Table 1). Figure 5
highlights 5 of the 16 larger cases analyzed in-depth. It was determined that
the tensile strength and tensile modulus predictions deviated less than 20%
from the IFF-R simulations in all larger cases except for the structures that
were compressed and had fracture defects (Figure 5 c, d), and for structures
that had low densities and diameters much larger (i.e. 10 times larger
diameters) than what was used in the training set. In Figure 5a, the MD
structure is modelled after an experimentally determined CNT/polymer cross-
section[5], and has CNT diameters of 10-12 nm (about 10 times larger than in
the training set). The accurate predictions (4.0% strength deviation and 7.5%
modulus deviation) for Figure 5a can be attributed to the close packing of the
CNTs, exemplified by the typical dog-bone formation of CNTs with diameters
larger than 4 Angstroms[51]. Figure 5b demonstrates another simulation modeled
after an experimentally determined cross-section of a double-wall CNT bundle
with hexagonal packing and diameters of 1.45 nm[44]. The structure in Figure
5b contains 64 CNTs in the bundle and 78080 atoms (about 6 times larger than
the majority of the training set). However, HS-GNN predicts a strength value
for this structure with a value of about 3% from the IFF-R prediction,
suggesting that there is limited influence on the size of the input structure
when predicting tensile strength for pristine CNT morphologies. When
predicting the strength of large (>64 nanotubes), compressed CNT structures
with fracture defects (Figure 5 c,d) the strength predictions of HS-GNN
deviate 20% or more. This suggests that complex morphologies with simultaneous
fracture defects greatly diminish the performance of HS-GNN more than any
other observed factor. In contrast, when analyzing a similar morphology to
Figures 5c,d without fracture defects (as shown in Figure 5e) it is shown that
HS-GNN is capable of predicting the tensile strength with less than 6%
deviation from IFF-R and tensile modulus with less than 6.5% deviation from
IFF-R. It should be mentioned, in the test set, cases 1051-1053 also had
diameters of 12 nm. The strength and modulus deviations were >30% for the
cases 1051-1053, suggesting that CNT spacing in the input structure has some
influence on prediction accuracy. In summary, it is shown that the number of
atoms is of little concern when generalizing to larger structures as long as
the nanotube diameters are within the range of the training set (0.271-2.98
nm), and if there few, if any, fracture type defects.

##### Ablation study.

Besides the 3 GNN based baselines and our HS-GNN model, we also evaluate two
GNN setups. In the first setup, we found hierarchies within a CNT structure as
HS-GNN, but we didn’t take topology and geometry information of each cluster
as super-node features or edge attention. Its prediction error for test set on
strength and modulus are 7.4% and 14.6%. In the second setup, we didn’t design
the hierarchical neural network anymore. But we computed persistence summaries
and PCA information of neighborhood of each point, and took them as the input
features of nodes or edge attention in GNN. Its prediction error on strength
and modulus are 10.7% and 19.5%. Compared to performances of vanilla-GNN and
HS-GNN, we can conclude that both hierarchical design capturing long range
interactions among atoms and local topological and geometric information
contribute to our HS-GNN model. The improvement from hierarchies plays a more
important role than local topology and geometry.

Table 1: Mean squared error in prediction of strength and modulus, using 90% of the data for training. Mechanical property prediction error is also provided for larger cases (>20,000 atoms) not included in the data set and distinguished here by the proceeding (L). Our ML method was trained in 3 ways: HS-GNN-A excludes using any mechanical property predictions as input features, HS-GNN-B is trained using predicted strength to predict modulus, and HS-GNN-C is trained using predicted strength and a pretrained model to predict modulus. Predictions using HS-GNN methods are then compared to 6 baselines (right of the double vertical line). The HS-GNN is 2 to 5 times more accurate even for larger structures not included in the training set. Uncertainties under 10% can be considered competitive with experimental measurements, which often have similar errors.  Prediction error | HS-GNN-A | HS-GNN-B | HS-GNN-C | LR | PointNet | RS-CNN | Vanilla-GNN | MPNN | DimeNet
---|---|---|---|---|---|---|---|---|---
Strength (%) | 4.1 | 4.1 | 4.1 | 21.3 | 15.18 | 11.60 | 14.5 | 12.02 | 8.6
Modulus (%) | 8.8 | 8.2 | 7.6 | 37.4 | 30.01 | 20.55 | 31.2 | 19.72 | 16.4
Strength (L) (%) | 12.7 | 12.7 | 12.7 | 41.5 | 32.3 | 24.5 | 35.0 | 27.2 | 20.9
Modulus (L) (%) | 17.6 | 15.8 | 15.3 | 42.1 | 37.4 | 30.5 | 40.9 | 34.8 | 23.5
Figure 5:  Comparison of 5 larger carbon nanostructures and calculated tensile
strength and tensile modulus from molecular dynamics simulation (IFF-R) with
the ML-prediction using the HS-GNN method (under MD image). (a) Transmission
electron microscopy (TEM) image of a CNT composite cross-section taken by
Jolowsky et al. [5] (left) and corresponding models of flattened DWCNTs in MD
simulation (right). The deviations are on the order of 5%. (b) Hexagonal
packing of a DWCNT bundle observed by TEM by Colomer et al. [44] (left) and
corresponding models in MD simulation (right). Differences between MD
simulation and ML are small for strength and larger for modulus. (c-e)
Defective structures derived from compression of SWCNT bundles. SWCNT bundle
of sizes 72 (c) and 64 (d) nanotubes, each containing some nanotubes with
fracture defects, and a SWCNT bundle of 64 nanotubes without defects. The
model size is larger than the training set and uncertainties in ML predictions
remain in the range of 4% to 25%.

##  3 Conclusion

A hierarchical spatial graph neural network (HS-GNN) has been developed and
trained to predict the tensile strength and modulus of computer-generated
models of pristine and defective carbon nanotube and graphitic assemblies. We
used an extensible data set of 1159 unique molecular dynamics (MD) simulations
of tensile stress-strain curves that were carried out using the IFF-R
molecular dynamics force field, which accurately represents chemical bonding,
physical properties, and bond dissociation. The average deviation using the
HS-GNN is 4% for strength and 9% for tensile modulus relative to the MD
simulation, and by a factor of 2 to 5 lower than mechanical property
predictions using other methods such as DimeNet, vanilla-GNN, PointNet, and
linear regression. The incorporation of spatial information and use of a
hierarchical method to process graph information is adequate for generalizing
the neural network property predictions to larger graphitic structures outside
the training set and still achieve less than 18% deviation in ML-predicted
mechanical property values relative to MD simulation. The ML model can be used
for rapid screening of the mechanical performance of CNT morphologies and
understanding the relationship between nanostructural features and carbon
nanotube bundle performance for carbon fiber yarns used in automotive and
aerospace applications. It is yet to be determined if the HS-GNN can be
generalized to accurately predict mechanical properties for nanostructures
that are structurally similar but chemically different (e.g. boron nitride,
macromolecular structures and composites), or structures that are chemically
similar but structurally different (e.g. diamond nanofibers).

##  4 Methods

###  4.1 Building the Model Structure Database

Models of carbon nanotube structures, assemblies, and defective graphitic
structures were generated using the Materials Studio Graphical User
Interface.[52] A total of 1159 different structures were created to build the
structure database of CNTs and graphitic models. The CNT systems of interest
include individual and bundled single-wall, double-wall, and triple-wall CNTs.
CNT bundles are defined as a group of more than one CNT. In general, the
bundles ranged from 2 to 12 nanotubes in a simulation. Some larger structures
were generated for generalizing our developed ML method, which contained up to
72 CNTs in a simulation. The diameter of the CNTs were modified by changing
the chiral index values (N, M) using Materials Studio’s "Build Nanostructure"
tool. For simplicity, all structures were generated such that N = M. CNTs with
N and M values from 2 to 21 were evaluated. The initial outer wall diameters
of these CNTs ranged from 0.271 nm to 2.98 nm. In general, the height of the
CNTs ranged from 2.42 nm to 2.52 nm. To introduce complex structural features,
some CNT systems were compressed to a height of 2.29 nm before being strained
in tension. The mechanical properties of pristine nanotubes and defective
nanotubes were of interest. A defective nanotube is defined as a CNT that has
one or more missing atoms such that some carbon atoms in the nanotube have 2
instead of 3 bonded neighbors, or the CNT lattice is structured such that
pentagons and heptagons are formed instead of the typical hexagonal
arrangement. Defects were implemented in three ways: firstly, by selectively
severing bonds along the radial direction of the CNT such that nanotube is
broken into two segments; secondly, by semi-randomly deleting atoms in the CNT
structures effectively mimicking mono-, di-, and tri-vacancy defects commonly
afflicting CNTs after synthesis, stress, or alignment; thirdly, by re-
configuring the bonding of carbon atoms in the CNT lattice to form 2 heptagons
and 2 pentagons, otherwise known as a Stone-Wales defect. In some cases,
structures with 0 modulus and 0 strength were of interest, and these
structures were obtained from the data file created after the CNT ruptured in
tension during the simulation. IFF-R parameters were used to define the CNT
atomic properties (supporting information S0). Each carbon atom’s force field
type was defined as cg1.

###  4.2 Choice of Force Field

In order to reproduce surface energies of graphitic materials that deviate
less 5 percent from experimentally determined values, IFF incorporates virtual
pi electrons on the corresponding carbon atoms. However, IFF is unable to
simulate bond dissociation because of the use of a harmonic bond energy
potential. To simulate bond dissociation IFF-R utilizes a Morse bond energy
potential, which offers a better description for atomic bond behavior at large
displacements. IFF-R maintains the use of virtual pi electrons and has been
shown to reproduce experimentally determined surface energies and mechanical
properties for graphite and CNTs, making it a suitable force field for
investigating the influence of various nanostructural features on CNT
mechanical performance (IFF-R reference). However, it should mentioned that
virtual pi electrons were not considered for the current developed ML model
for the sake of simplicity and ease of CNT structure generation.

###  4.3 Reactive Molecular Dynamics Simulation to Create Training and
Validation Data

All structures were exported as .car/.mdf files so that they could be
converted to a Large-scale Atomic/ Molecular Massively Parallel Simulator
(LAMMPS) readable data file using the msi2lmp tool (part of the LAMMPS
release).[53] The data file contains useful structural information such as the
height of each CNT, the number of atoms in an individual CNT or a CNT bundle,
and so on. The property predictions made from the developed machine learning
pipeline used the data files as inputs. A data file was created for the
initial CNT structures, equilibriated CNT structures, compressed CNT
structures (where applicable), and the ruptured CNT structures.

A series of 1159 molecular dynamics simulations was carried out using the
LAMMPS molecular dynamics simulation software. All simulations were run at 298
K using the canonical ensemble and Nose-Hoover thermostat and barostat. The
Morse bond parameters were used for the carbon-carbon graphitic type bonds
defined in IFF-R. The angles, dihedrals, and improper terms used class II
potentials with a global cutoff of 12.0 Å. Each simulation was minim zed with
an energy stopping criteria of 1e-4 and a force stopping criteria of 1e-6
kcal/molÅ. Before the tensile simulation began each simulation was allowed to
equilibrate for 10000 femtoseconds. When applicable, the simulation cell was
reduced in size after equilibration for 3500 femtoseconds and at a rate of
-20/ns (7% of the original CNT height) to form the compressed structures. The
tensile simulation was run until failure (no more than 100,000 femtoseconds
were needed) and at an engineering strain rate of 20/ns.

###  4.4 Generation of Features

###  4.5 Molecular Dynamics Mechanical Property Calculations

The mechanical properties were calculated using Microsoft Excel and an example
can be found in the supporting information. The stress, strength, and Young’s
Modulus were calculated in units of GPa. To calculate stress the pressure
tensor along the axial direction of the CNT (zz-direction) was normalized to
account for only the stress experienced by the CNTs. The normalization
involves computing the total cross-sectional area of the CNTs by using

$A=N\pi(r+0.19)^{2}$

where N is the number of CNTs and r is the CNT radius plus a van der Waals
correction of 0.19. The total area of the simulation cell is then divided by
the normalized area (A) to get the correction term, C. Stress can then be
defined as

$\sigma=|C\Delta P_{zz}|\;\;\;where\;\;\;\Delta P_{zz}=P_{t}-P_{0}$

here the subscript t is the value at any time step and the subscript 0 is the
initial value after equilibration and compression. Strain is calculated by

$\epsilon=\frac{\Delta L_{zz}}{L_{0}}\;\;\;where\;\;\;\Delta
L_{zz}=L_{t}-L_{0}$

The strength was determined to be the maximum corrected stress value and
modulus was calculated by taking the slope of the stress-strain curve where
the response is linear (e.g. between 0-0.01 strain). For the compressed
structures, the Young’s Modulus was calculated after the CNTs were stretched
to their original length (e.g. after 0.07 strain).

###  4.6 HS-GNN Neural network architecture

The overall HS-GNN pipeline is given in Figure 3 with a brief description in
Section 2.2. We now describe the components in more details.

####  4.6.1 Hierarchical heterogeneous graph formation

##### Heterogeneous graph representation of CNT structures.

A CNT structure can be modeled as a graph in different ways. The bond graph
$G_{bond{\texttt{S}}}$ captures the information of chemical bond connections
which plays a pivot role in molecular structures. However, it ignores the
spatial relations between nearby but non-bonded atoms (e.g, two nearby carbon
atoms from neighboring nanotubes). This issue can be addressed by connecting
nearby pairs of atoms, namely, two graph nodes are connected if their
corresponding atoms are within certain Euclidean distance $c_{e}$ to each
other. In addition, effective resistance distance [54] between two nodes in a
graph (network) can measure how well the two nodes are connected via paths in
the graph, which can reflect the graph topology. Thus we also connect pairs of
nodes if their effective resistance distance in the bound graph
$G_{bond{\texttt{S}}}$ is smaller than a threshold $c_{r}$. Finally, we also
connect nodes whose corresponding atoms form a dihedral angle. All these
together give rise to a _heterogeneous graph_ $G_{\texttt{S}}=(V,E)$ with four
different types of edges $E=E_{b}\cup E_{e}\cup E_{r}\cup E_{d}$, representing
the edge sets formed by chemical bonds, Euclidean distance, effective
resistance distance, and dihedrals, respectively.

##### Creating hierarchical representations.

Existing GNN models suffer the so-called over-smoothing issue, where signals
on graph are quickly smoothed out as the number of layers increase. Thus a GNN
usually has only a small number of layers, causing the receptive field of each
node limited to a relatively local region when a GNN model is applied on large
size graphs. We design a hierarchical architecture leveraging spatial
geometric information to tackle this issue.

An input CNT S can be viewed as a point cloud $V_{\texttt{S}}$ with each point
representing the center of an atom in S. We first apply a hierarchical
clustering algorithm on this 3D point cloud $V_{\texttt{S}}$, and obtain a
series of coarser and coarser point sets $V^{(1)}=V_{\texttt{S}}\subset
V^{(2)}\subset\ldots\subset V^{(k)}$. In particular, $V^{(i+1)}$ is obtained
by taking a so-called $\delta_{i}$-net of point set $V^{(i)}$, for each
$i\in[1,k)$. The computation of a $\delta$-net $Q$ of a point set $V$ is given
in Algorithm 1 in Supplement C. Note that each point (called a _super-node_)
in $V^{(i+1)}$ corresponds to a cluster of points in $V^{(i)}$.

The coarser point set $V^{(i)}$ will serve as the node set for a coarser graph
$G^{(i)}$ in level $L_{i}$ for any $i\in[1,k]$. Set $G^{(1)}=G_{\texttt{S}}$
as the graph representation in level $L_{1}$. For any $i>1$, we connect two
super-nodes in $V^{(i)}$ if there exists overlapping nodes from their
corresponding clusters (recall each cluster consists of a set of nodes in
$V^{(i-1)}$) in $G^{(i-1)}$.

####  4.6.2 Heterogeneous GNN in $L_{1}$ level hierarchy

The input to the GNN in level $L_{1}$ is the heterogeneous graphs
$G_{\texttt{S}}$ with 4 types of edges $E=E_{b}\cup E_{e}\cup E_{r}\cup
E_{d}$. The GNN used in level $L_{1}$ consists of $q$ GIN layers [55],
followed by $p$ GAT layers [38], and it will train different weights for
different types of edges. (See Supplement B for the message
passing/aggregation in the GIN or GAT layers.) We now describe this in details
using the edge set $E_{b}$ as an example:

Suppose there are $N_{\texttt{S}}$ nodes in $G_{\texttt{S}}$
($|V|=N_{\texttt{S}}$). We denote the node (feature) representation matrix in
the $l$-th ($l=1,2,...,q+p$) layer as
$H_{l}=[h_{1}^{(l)},h_{2}^{(l)},...,h_{N_{\texttt{S}}}^{(l)}]^{T}$ in which
$h^{(l)}_{i}$ is feature representation of the $i$-th node in the $\ell$-th
layer. The input node features are denoted as
$H_{0}=[h_{1}^{(0)},\ldots,h_{N_{\texttt{S}}}^{(0)}]^{T}$, where each
$h_{i}^{(0)}$ consists of the $i$-th atom’s 3D coordinates, degree, and
randomly generated features.

For $l=1,...,q$, $G_{\texttt{S}}$ is processed by GIN layers (Equation (10) in
Supplement), and we process $G_{\texttt{S}}$ (where nodes are equipped with
feature representations $H_{q}$ output from the $q$th GIN layer) by $p$ GAT
layer for $l=q+1,...,q+p$. From Equation (11) in Supplement we obtain the
final node representations $H_{b}=[h_{1;b},h_{2;b},...,h_{N_{s};b}]^{T}$
according to chemical bond edge set $E_{b}$.

Using the same GNN architectures as described in Equations (10) and (11) on
edge sets $E_{d}$, $E_{e}$ and $E_{r}$, respectively, we have node
representations $H_{*}=[h_{1;*},h_{2;*},...,h_{N_{s};*}]^{T}$ (* denotes
$d,e,r$). We get final node representations
$H=[h_{1},h_{2},...,h_{N_{s}}]^{T}$ by integrating the 4 kinds representations
into the convolution layer as follows:

$\bar{h}^{(1)}_{u}=\mathrm{ReLU}(W^{\prime}(h_{u;b}||h_{u;d}||h_{u;e}||h_{u;r}))$
(1)

where $\cdot||\cdot$ stands for concatenation.

####  4.6.3 Spatial information enhanced GNN in higher level hierarchies
(i.e., $L_{i}$ for $i>1$)

Level $L_{1}$ outputs node feature representations $\bar{h}^{(1)}_{u}$ as in
Eqn (1). In general, suppose we have already finished processing level
$L_{k-1}$ with the final node representations $\bar{h}^{(k-1)}_{u}$ for nodes
in $u\in V_{k-1}$ in graph $G^{(k-1)}$. In level $L_{k}$, the graph we will
process is the coarser graph $G^{(k)}=(V_{k},E_{k})$. In order to initialize
node features for nodes in $V_{k}\subseteq V_{k-1}$, recall that each node
$v\in V_{k}$ in fact corresponds to a cluster of nodes $C_{v}\subseteq
V_{k-1}$. Let $G^{(k-1)}(v)$ denote the subgraph of $G^{(k-1)}$ spanned by
nodes in $C_{v}$ – intuitively, this subgraph $G^{(k-1)}(v)$ from graph
$G^{(k-1})$ in level $L_{k-1}$ is collapsed into a single node $v$ in graph
$G^{(k)}$ in level $L_{k}$. We simply perform a max-pooling (see Equation (2))
of the node features of subgraph $G^{(k-1)}(v)$ to obtain a feature
representation $h_{G^{(k-1)}_{v}}$ for the entire subgraph $G^{(k-1)}(v)$.

$h_{G^{(k-1)}_{v}}=[\mathrm{max}(\\{\bar{h}^{(k-1)}_{u}[0]|u\in
V_{k-1}\\}),\mathrm{max}(\\{\bar{h}^{(k-1)}_{u}[1]|u\in
V_{k-1}\\}),\cdots,\mathrm{max}(\\{\bar{h}^{(k-1)}_{u}[d_{k-1}]|u\in
V_{k-1}\\})]$ (2)

where $d_{k-1}$ is dimension of $\bar{h}^{(k-1)}_{u}$. This representation
$h^{(k),0}_{v}=h_{G^{(k-1)}_{v}}$ is then used as the initial feature for node
$v\in V_{k}$.

Now we have the coarse graph $G^{(k)}$ with initial features
$H^{(k),0}=[h^{(k),0}_{1},\ldots,h^{(k),0}_{N_{k}}]^{T}$ for all
$N_{k}=|V_{k}|$ nodes in $V_{k}$. The GNN for level $L_{k}$ consists of $r$
GCN layers [37] as GCN layer is easier to add reweighting factors introduced
in the following than GIN layer. More precisely, in the $l$-th layer
($l=1,...,r$) in hierarchical level $L_{k}$, the message passing function a
node $u\in V_{k}$ is:

$m_{u}^{(l)}=\sigma(\sum_{v\in u}W^{(l)}h_{v}^{(l-1)})$ (3)

where $\sigma(\cdot)$ is a non-linear function like $\mathrm{ReLU}(\cdot)$,
and the parameter $W^{(l)}$ is a linear transformation matrix that will be
learned by training.

In our model, we add a reweighting factor $\tau_{uv}$ for the message between
two nodes, $u$ and $v$, during graph convolution:

$m_{u}^{(l)}=\sigma(\sum_{v\in N(u)}\tau_{uv}^{(l)}W^{(l)}h_{v}^{(l-1)})$ (4)

This reweighting factor can tell the difference and significance of any
message based on local topological and geometric features at nodes, or in the
corresponding subgraphs. We describe a method to learn this factor in
Supplement C. After the message aggregation, we take a 3 layer MLP as the
update function in a GCN.

$h_{u}^{(l)}=\mathbf{MLP}(h_{u}^{(l-1)},m_{u}^{(l)})$ (5)

After 3 layers convolution, we could obtain a graph representation for any
subgraph $G^{\prime}$ of $G^{(k)}$ based on the node representation of nodes
in $G^{\prime}$ through max-pooling the same as what used before for level
$L_{1}$.

In the highest level hierarchy, we add one regression layer predicting the
tensile properties taking both graph representations and the global geometric
and topological features as input in the final step. Specifically, we take the
global persistence summaries (see Supplement A) of all the atoms in the CNT
bundles in the final regression layer.

###  4.7 Property prediction

Dataset:

Our dataset has 1159 CNT bundles which we split into training and test dataset
as the ratio of 9:1. Each CNT bundle consists of one or more carbon nanotubes,
and preserves 3D positions of all carbon atoms ($(x,y,z)$ coordinates),
chemical bonds between atoms as well as chemical dihedrals among atoms. The
sizes of CNT bundles vary from 480 atoms to over 10,000 atoms. We denote those
features a CNT bundle as $(V,E_{b},E_{d})$, where
$V=\\{(x_{i},y_{i},z_{i})|i=1,\cdots,N\\}$ is the set of atoms, $E_{b}$ is the
set of chemical bonds and $E_{d}$ is the set of chemical dihedral relations.
Our task is to predict two mechanism properties, tensile strength and tensile
modulus, based on those input features.

Architecture:

We take a 3 level hierarchical HS-GNN. As we introduced in Section (4.6), we
take $\delta$-nets to get the coarser graphs. In particular, we take a
$\delta_{1}$-Net to obtain the level $L_{2}$ graph $G^{(2)}$ from the original
level $L_{1}$ graph $G^{(1)}=G_{\texttt{S}}$, and we take a $\delta_{2}$-Net
to obtain the level $L_{3}$ graph $G^{(3)}$ from $G^{(2)}$. To choose the
parameter $\delta_{1}$ and $\delta_{2}$, we randomly sample 100 CNT bundles
from training set, compute the average 3D Euclidean distance (denoted by
$d_{E}$) between every two nodes connected by chemical bonds, and set
$\delta_{1}=15d_{E}$ and $\delta_{2}=30d_{E}$.

We then design GNN architecture for each level. In the basic $L_{1}$ level, a
heterogeneous graph $G^{(1)}=(V,E_{b},E_{d},E_{e},E_{r})$ is formulated from
each point set $V$. If there exists a chemical bond between points $u$ and
$v$, we add an edge $(u,v)_{b}$ in $E_{b}$. If there two points $u$ and $v$
can formulate a chemical dihedral angle with a third point, we add an edge
$(u,v)_{d}$ in $E_{d}$. If Euclidean distance between two points $u$ and $v$
is smaller than a cut-off distance $c_{E}=6d_{E}$, we add an edge $(u,v)_{e}$
in $E_{e}$. In the final, we compute the average effective resistance distance
between every two nodes connected by chemical bonds in 100 sampled CNT bundles
and denote it as $d_{r}$. We add an edge $(u,v)_{r}$ in $E_{r}$ if effective
resistance distance between two points $u$ and $v$ are smaller than
$c_{r}=10d_{r}$. For each relation in $\\{E_{b},E_{d},E_{e},E_{r}\\}$, we take
$q$ GIN layers followed by $p$ GAT layers to process as we mentioned in
Section (4.6). In each higher level, a coarser graph is processed through $r$
reweighted GCN layers. Thus we have 3 hyperparameters, $q$, $p$ and $r$, to
tune. They are chosen from $\\{3,4,5,6\\}$, $\\{1,2,3\\}$, and $\\{3,4,5,6\\}$
respectively. We tune those hyperparamters by 5 fold cross-validation. In a
more explicit manner, we split the training set into 5 folds, each time we use
4 folds to train our model, and then compute the prediction error on the rest
1 fold as validation set. Repeat this operation 5 times so that each fold is
taken as validation set once. We take the average mean squared error over 5
iterations as the evaluation performance of a model with certain
hyperparameters. Finally we test the model and hyperparameters with best
evaluation performance on our test set and report the results.

As an exploratory study, we predict tensile strength and tensile modulus in
two ways. The first way is to apply generalized additive models [56] to
predict tensile strength and tensile modulus using only features (e.g., CNT
diameter, CNT height) separately. The second way is to first estimate causal
relations among tensile strength, tensile modulus, and features using a causal
discovery method, CAM [45], and then utilize the inferred causal relations
between tensile strength and tensile modulus in prediction. As shown in Figure
(2 g), the CAM finds that tensile strength is a direct cause of tensile
modulus. When we only use features to predict tensile modulus, the generalized
additive model yields a fitted $R^{2}$=0.931 (the higher, the better). When we
add tensile strength in predicting tensile modulus besides features, the
$R^{2}$ increases to 0.941. This exploration suggests that using tensile
strength helps the prediction of tensile modulus. We therefore design three
different final prediction setups: (A) We train two models to predict tensile
strength and tensile modulus independently. (B) We first train a model to
predict tensile strength. We then train a second model to predict tensile
modulus which takes predicted strength as part of input features in the final
layer. (C) We train a model to predict tensile strength first. We then take
the trained GNN parameters as initialization to train a second model to
predict tensile modulus which also takes the predicted strength as input
features in the final layer to predict modulus.

## Acknowledgments

The authors acknowledge support by the National Science Foundation
(OAC-1931587, CMMI-1940335) and the University of Colorado Boulder. Work by
Wang and Zhao are partially supported by NSF under grants OAC-2039794 and
CCF-2051197. The allocation of computational resources is acknowledged at the
Argonne Leadership Computing Facility, which is a DOE Office of Science User
Facility supported under contract DE-AC02-06CH11357, and at the Summit
supercomputer supported by the National Science Foundation (ACI-1532235 and
ACI-1532236).

## References

  * [1] H. Chang, J. Luo, P. V. Gulgunje, and S. Kumar.  Structural and functional fibers.  Ann. Rev. Mater. Res., 47(13):13.1–13.29, 2017.
  * [2] Céline Merlet, Benjamin Rotenberg, Paul A. Madden, Pierre-Louis Taberna, Patrice Simon, Yury Gogotsi, and Mathieu Salanne.  On the molecular origin of supercapacitance in nanoporous carbon electrodes.  Nat. Mater., 1(4):306–310, 2012.
  * [3] S. P. Koenig, L. D. Wang, J. Pellegrino, and J. S. Bunch.  Selective molecular sieving through porous graphene.  Nat. Nanotechnol., 7(11):728–732, 2012.
  * [4] Y. X. Bai, R. F. Zhang, X. Ye, Z. X. Zhu, H. H. Xie, B. Y. Shen, D. L. Cai, B. F. Liu, C. X. Zhang, Z. Jia, S. L. Zhang, X. D. Li, and F. Wei.  Carbon nanotube bundles with tensile strength over 80 gpa.  Nat. Nanotechnol., 13(7):589–593, 2018.
  * [5] C. Jolowsky, R. Sweat, J. G. Park, A. Hao, and R. Liang.  Microstructure evolution and self-assembling of cnt networks during mechanical stretching and mechanical properties of highly aligned cnt composites.  Compos. Sci. Technol., 166:125–130, 2018.
  * [6] C. Pramanik, T. Jamil, J. R. Gissinger, D. Guittet, P. J. Arias-Monje, S. Kumar, and H. Heinz.  Polyacrylonitrile interactions with carbon nanotubes in solution: Conformations and binding as a function of solvent, temperature, and concentration.  Adv. Funct. Mater., 29:1905247, 2019.
  * [7] Jordan J. Winetrout, Krishan Kanhaiya, Geeta Sachdeva, Ravindra Pandey, Behzad Damirchi, Adri van Duin, Gregory Odegard, and Hendrik Heinz.  Implementing reactivity in molecular dynamics simulations with the interface force field (iff-r) and other harmonic force fields.  arXiv preprint arXiv:2107.14418, 2021.
  * [8] G. R. Schleder, A. C. M. Padilha, C. M. Acosta, M. Costa, and A. Fazzio.  From dft to machine learning: Recent approaches to materials science-a review.  J. Physics: Materials, 2:032001, 2019.
  * [9] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25:1097–1105, 2012.
  * [10] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [11] Alex Graves, Abdel-rahman Mohamed, and Geoffrey Hinton.  Speech recognition with deep recurrent neural networks.  In 2013 IEEE international conference on acoustics, speech and signal processing, pages 6645–6649. Ieee, 2013.
  * [12] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, pages 6000–6010, 2017.
  * [13] Babak Alipanahi, Andrew Delong, Matthew T Weirauch, and Brendan J Frey.  Predicting the sequence specificities of dna-and rna-binding proteins by deep learning.  Nature biotechnology, 33(8):831–838, 2015.
  * [14] Ngoc Giang Nguyen, Vu Anh Tran, Duc Luu Ngo, Dau Phan, Favorisen Rosyking Lumbanraja, Mohammad Reza Faisal, Bahriddin Abapihi, Mamoru Kubo, Kenji Satou, et al.  Dna sequence classification by convolutional neural network.  Journal of Biomedical Science and Engineering, 9(05):280, 2016.
  * [15] Michael I Jordan and Tom M Mitchell.  Machine learning: Trends, perspectives, and prospects.  Science, 349(6245):255–260, 2015.
  * [16] Roberto Todeschini and Viviana Consonni.  Handbook of molecular descriptors, volume 11.  John Wiley & Sons, 2008.
  * [17] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller, and O Anatole Von Lilienfeld.  Fast and accurate modeling of molecular atomization energies with machine learning.  Physical review letters, 108(5):058301, 2012.
  * [18] Katja Hansen, Franziska Biegler, Raghunathan Ramakrishnan, Wiktor Pronobis, O Anatole Von Lilienfeld, Klaus-Robert Muller, and Alexandre Tkatchenko.  Machine learning predictions of molecular properties: Accurate many-body potentials and nonlocality in chemical space.  The journal of physical chemistry letters, 6(12):2326–2331, 2015\.
  * [19] Matthew Hirn, Stéphane Mallat, and Nicolas Poilvert.  Wavelet scattering regression of quantum chemical energies.  Multiscale Modeling & Simulation, 15(2):827–863, 2017.
  * [20] Albert P Bartók, Risi Kondor, and Gábor Csányi.  On representing chemical environments.  Physical Review B, 87(18):184115, 2013.
  * [21] David Duvenaud, Dougal Maclaurin, Jorge Aguilera-Iparraguirre, Rafael Gómez-Bombarelli, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams.  Convolutional networks on graphs for learning molecular fingerprints.  In Proceedings of the 28th International Conference on Neural Information Processing Systems-Volume 2, pages 2224–2232, 2015.
  * [22] Atsuto Seko, Hiroyuki Hayashi, Keita Nakayama, Akira Takahashi, and Isao Tanaka.  Representation of compounds for machine-learning prediction of physical properties.  Physical Review B, 95(14):144110, 2017.
  * [23] Sandip De, Albert P Bartók, Gábor Csányi, and Michele Ceriotti.  Comparing molecules and solids across structural and alchemical space.  Physical Chemistry Chemical Physics, 18(20):13754–13769, 2016.
  * [24] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [25] Ian Goodfellow, Yoshua Bengio, Aaron Courville, and Yoshua Bengio.  Deep learning, volume 1.  MIT press Cambridge, 2016.
  * [26] David Weininger, Arthur Weininger, and Joseph L Weininger.  Smiles. 2. algorithm for generation of unique smiles notation.  Journal of chemical information and computer sciences, 29(2):97–101, 1989.
  * [27] Stanisław Jastrzębski, Damian Leśniak, and Wojciech Marian Czarnecki.  Learning to smile (s).  arXiv preprint arXiv:1602.06289, 2016.
  * [28] Zheng Xu, Sheng Wang, Feiyun Zhu, and Junzhou Huang.  Seq2seq fingerprint: An unsupervised deep molecular embedding for drug discovery.  In Proceedings of the 8th ACM international conference on bioinformatics, computational biology, and health informatics, pages 285–294, 2017.
  * [29] Benjamin Sanchez-Lengeling and Alán Aspuru-Guzik.  Inverse molecular design using machine learning: Generative models for matter engineering.  Science, 361(6400):360–365, 2018.
  * [30] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [31] Yongcheng Liu, Bin Fan, Shiming Xiang, and Chunhong Pan.  Relation-shape convolutional neural network for point cloud analysis.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8895–8904, 2019.
  * [32] Ryan S DeFever, Colin Targonski, Steven W Hall, Melissa C Smith, and Sapna Sarupria.  A generalized deep learning approach for local structure identification in molecular simulations.  Chemical science, 10(32):7503–7515, 2019.
  * [33] Kui Xu, Zhe Wang, Jianping Shi, Hongsheng Li, and Qiangfeng Cliff Zhang.  A2-net: Molecular structure estimation from cryo-em density volumes.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 1230–1237, 2019.
  * [34] Niklas WA Gebauer, Michael Gastegger, and Kristof T Schütt.  Symmetry-adapted generation of 3d point sets for the targeted discovery of molecules.  arXiv preprint arXiv:1906.00957, 2019.
  * [35] Steven Kearnes, Kevin McCloskey, Marc Berndl, Vijay Pande, and Patrick Riley.  Molecular graph convolutions: moving beyond fingerprints.  Journal of computer-aided molecular design, 30(8):595–608, 2016\.
  * [36] Connor W Coley, Regina Barzilay, William H Green, Tommi S Jaakkola, and Klavs F Jensen.  Convolutional embedding of attributed molecular graphs for physical property prediction.  Journal of chemical information and modeling, 57(8):1757–1772, 2017\.
  * [37] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  International Conference on Learning Representations, 2017.
  * [38] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  International Conference on Learning Representations, 2018.
  * [39] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In Proceedings of the 34th International Conference on Machine Learning-Volume 70, pages 1263–1272. JMLR. org, 2017.
  * [40] Johannes Klicpera, Janek Groß, and Stephan Günnemann.  Directional message passing for molecular graphs.  arXiv preprint arXiv:2003.03123, 2020.
  * [41] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  Journal of chemical information and modeling, 59(8):3370–3388, 2019\.
  * [42] Hendrik Heinz, Tzu-Jen Lin, Ratan K. Mishra, and Fateme S. Emami.  Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: The interface force field.  Langmuir, 29(6):1754–1765, 2013.
  * [43] C. Pramanik, J. R. Gissinger, S. Kumar, and H. Heinz.  Carbon nanotube dispersion in solvents and polymer solutions: Mechanisms, assembly, and preferences.  ACS Nano, 11:12805–12816, 2017.
  * [44] J.F. Colomer, L. Henrard, G. Van Tendeloo, A. Lucas, and P. Lambin.  Study of the packing of double-walled carbon nanotubes into bundles by transmission electron microscopy and electron diffraction.  J. Mater. Chem., 14:603–606, 2004.
  * [45] Peter Bühlmann, Jonas Peters, and Jan Ernest.  Cam: Causal additive models, high-dimensional order search and penalized regression.  The Annals of Statistics, 42(6):2526–2556, 2014.
  * [46] Qimai Li, Zhichao Han, and Xiao-Ming Wu.  Deeper insights into graph convolutional networks for semi-supervised learning.  In Thirty-Second AAAI conference on artificial intelligence, 2018\.
  * [47] Rex Ying, Jiaxuan You, Christopher Morris, Xiang Ren, William L Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 4805–4815, 2018.
  * [48] Jingjia Huang, Zhangheng Li, Nannan Li, Shan Liu, and Ge Li.  Attpool: Towards hierarchical feature representation in graph convolutional networks via attention mechanism.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 6480–6489, 2019.
  * [49] Herbert Edelsbrunner, David Letscher, and Afra Zomorodian.  Topological persistence and simplification.  In Proceedings 41st annual symposium on foundations of computer science, pages 454–463. IEEE, 2000.
  * [50] Herbert Edelsbrunner and John Harer.  Computational topology: an introduction.  American Mathematical Soc., 2010.
  * [51] James A Elliott, Jan KW Sandler, Alan H Windle, Robert J Young, and Milo SP Shaffer.  Collapse of single-wall carbon nanotubes is diameter dependent.  Physical Review Letters, 92(9):095501, 2004.
  * [52] Biovia/Dassault Systemes.  Materials Studio 2019 Program Suite and User Guides.  Cambridge, UK, 2019.
  * [53] S. Plimpton.  Fast parallel algorithms for short-range molecular dynamicss.  J. Comput. Phys., 117:1–19, 1995.
  * [54] D Babić, DJ Klein, I Lukovits, Sonja Nikolić, and N Trinajstić.  Resistance-distance matrix: A computational algorithm and its application.  International Journal of Quantum Chemistry, 90(1):166–176, 2002\.
  * [55] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  International Conference on Learning Representations, 2019.
  * [56] Simon N Wood.  Generalized additive models: an introduction with R.  CRC press, 2017.
  * [57] Peter Bubenik.  Statistical topological data analysis using persistence landscapes.  Journal of Machine Learning Research, 16(1):77–102, 2015.
  * [58] Genki Kusano, Kenji Fukumizu, and Yasuaki Hiraoka.  Kernel method for persistence diagrams via kernel embedding and weight factor.  Journal of Machine Learning Research, 18(189):1–41, 2018.
  * [59] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Shipman, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier.  Persistence images: a stable vector representation of persistent homology.  Journal of Machine Learning Research, 18:218–252, 2017.
  * [60] Mathieu Carrière, Marco Cuturi, and Steve Oudot.  Sliced Wasserstein kernel for persistence diagrams.  International Conference on Machine Learning, pages 664–673, 2017\.
  * [61] Qi Zhao and Yusu Wang.  Learning metrics for persistence-based summaries and applications for graph classification.  In Advances in Neural Information Processing Systems, pages 9855–9866, 2019.
  * [62] Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume 2, pages 729–734. IEEE, 2005.
  * [63] F Scarselli, M Gori, Ah Chung Tsoi, M Hagenbuchner, and G Monfardini.  The graph neural network model.  IEEE Transactions on Neural Networks, 1(20):61–80, 2009.
  * [64] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  International Conference on Learning Representations, 2014.
  * [65] Mathias Niepert, Mohamed Ahmed, and Konstantin Kutzkov.  Learning convolutional neural networks for graphs.  In International conference on machine learning, pages 2014–2023, 2016.
  * [66] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In Advances in Neural Information Processing Systems, pages 1024–1034, 2017.

##  Appendix A Persistent Homology

Persistent homology [49, 50] is one of the most important developments in the
field of topological data analysis, and persistence summaries are often used
as features in statistics or machine learning tasks associated with graphs,
point clouds and 3-D shapes [57, 58, 59, 60, 61]. Suppose we are given a
topological space $X$ and a $filtration$, a sequence of growing subsets, of
$X$: $X_{1}\subseteq X_{2}\cdots\subseteq X_{n}=X$. As we inspect $X$ through
this filtration, sometimes a new topological feature like a void or a loop is
created when entering $X_{i}$, and destroyed in $X_{j}$. Persistent homology
can capture the birth and death of topological features in the form of a
$persistence$ $diagram$ $\mathrm{Dg}X$. Specifically, the $k$-dimensional
persistence diagram $\mathrm{Dg}_{k}X$ consists of a multi-set of persistence
points in the birth-death plane. Each persistence point $(b,d)$ indicates that
a $k$-dimensional topological feature appears when entering $X_{b}$ and
disappears upon entering $X_{d}$. The $persistence$ of this feature is its
lifespan $|d-b|$. In other words, the persistence diagram provides a simple
yet rich summary for the entire evolution of space $X$ through the lens of the
filtration $X_{1}\subseteq X_{2}\cdots X_{n}=X$, encoding multi-scale features
in $X$.

When $X$ is a point cloud, one common approach to obtain persistent summaries
is constructing a so-called _Vietoris-Rips filtration_. Here, a space is
modeled by a _simplicial complex_ spanned by a vertex set $V$: Roughly
speaking, a $k$-dimensional simplex is the $k$-dimensional generalization of
vertices ($0$-D), edges ($1$-D) and triangles ($2$-D simplices). A simplicial
complex is then simply a union of simplices with the condition that if a
simplex is contained in this complex, then any of its face will also be in the
complex.

Given a set of point $V\subseteq\mathbb{R}^{n}$, the Vietoris-Rips complex at
scale $r$ consists of all simplices with diameter less than $r$:

$\mathrm{VR}_{r}(V)=\\{\sigma\subset V\mid\forall u,v\in\sigma,||u-v||\leq
r\\}$ (6)

In particular, $\mathrm{VR}_{0}(V)=\\{\\{u\\}||u\in V\\}$, and
$\mathrm{VR}_{\infty}(V)$ consist of all simplices spanned by vertices in $V$.
By increasing $r$ from 0 to $\infty$, we obtain a filtration
$\mathrm{VR}_{0}(V)\subseteq\mathrm{VR}_{r_{1}}(V)\subseteq\mathrm{VR}_{r_{2}}(V)\subseteq...\subseteq\mathrm{VR}_{\infty}(V)$
($0\leq r_{1}\leq r_{2}\leq...$). We use this Vietoris-Rips filtration on
point clouds, see Fig 6, to track molecular structures topological features in
our experiments.

Figure 6: Vietoris-Rips filtration on point clouds. We create a ball centered
at each point in the 2D space, and connect two points when their corresponding
balls intersect. As radius of balls increase, two voids in the 2D plane appear
and then disappear. Their birth and death are recorded as two red persistence
points in persistence diagram.

In order to facilitate the downstream machine learning tasks, a persistence
diagram is usually converted a fixed dimensional vector [57, 58, 60]. We
deploy persistence image [59] in our approach, which is a stable and powerful
persistence diagrams vectorization approach developed by Adams et al. Set
$T:{\mathbb{R}}^{2}\to{\mathbb{R}}^{2}$ to be the linear transformation
$T(x,y)=(x,y-x)$. Given a persistence diagram ${\mathrm{D}}$, let
$T({\mathrm{D}})$ Denote the transformed diagram of a persistence diagram
${\mathrm{D}}$ as $T({\mathrm{D}})$. Let
${\phi}_{u}:{\mathbb{R}}^{2}\to\mathbb{R}$ be a differentiable distribution
function whose mean locates at $u\in{\mathbb{R}}^{2}$: For example, in our
implementation later, we will use the Gaussian distribution function (a common
choice), where for any $z\in{\mathbb{R}}^{2}$,
${\phi}_{u}(z)=\frac{1}{2\pi{\tau}^{2}}e^{-\frac{||z-u||^{2}}{2{\tau}^{2}}}$.

######  Definition 1 (Persistence images)

Let $\alpha:{\mathbb{R}}^{2}\to\mathbb{R}$ be a non-negative weight function
for the persistent plane ${\mathbb{R}}^{2}$. Given a persistence diagram
$\mathrm{Dg}X$, its _persistence surface_
${\rho}_{{\mathrm{D}}}:{\mathbb{R}}^{2}\to\mathbb{R}$ (w.r.t. $\alpha$) is
defined as: for any $z\in\mathbb{R}^{2}$,

$\rho_{{\mathrm{D}}}(z)=\sum_{u\in T({\mathrm{D}})}\alpha(u)\phi_{u}(z).$ (7)

The persistence image is a discretization of the persistence surface as
follows. Set a fixed grid within a rectangle in the plane with a collection
$\mathcal{P}$ of N pixels. The _persistence image_ for a persistence diagram
${\mathrm{D}}$ is
${\mathrm{PI}}_{{\mathrm{D}}}=\\{{\mathrm{PI}}_{\mathrm{D}}[p]\\}_{p\in\mathcal{P}}$,
where ${\mathrm{PI}}_{\mathrm{D}}[p]:=\int\int_{p}\rho_{\mathrm{Dg}(X)}dydx$.

Note that ${\mathrm{PI}}_{\mathrm{D}}$ can also be viewed as a vector in
${\mathbb{R}}^{N}$, and thus persistence images are naturally equipped with
the $L_{2}$-distance in $\mathbb{R}^{N}$.

##  Appendix B Graph Neural Networks

Graph Neural Networks (GNN) [62, 63, 64, 65] are the generation of neural
networks to the graph structured data. A node in a graph iteratively receives
information from its neighborhood and update its representation or features.
These node $\&$ edge representation information transferred between vertices
are called messages [39]. A transformation of the messages and updated
representations can be learned through training. This message passing scheme
can be formulated in a more explicit manner. Given an undirected graph
$G=(V,E)$ where $V$ is the node set and $E$ is the edge set, the input to a
GNN are node features $h_{u}^{0}$ and edge features $e_{uv}$ for every $u\in
V$ and $(u,v)\in E$. Then in the $t$-th hidden layer of the GNN, the forward
convolution consists of two functions:

$\displaystyle\mathrm{AGGREGATE}~{}~{}~{}~{}m^{t+1}_{u}$
$\displaystyle=f^{t}(h^{t}_{u},e_{uv},\\{h^{t}_{v}|v\in N(u)\\})$ (8)
$\displaystyle\mathrm{UPDATE}~{}~{}~{}~{}h^{t+1}_{u}$
$\displaystyle=g^{t}(m^{t+1}_{u},h^{t}_{u})$

where $N(u)$ is the neighborhood of node $u$, $f(\cdot)$ and $g(\cdot)$ are
message aggregation and update function respectively. Finally, if the task is
on the graph level, there is a readout function $r(\cdot)$ in the final layer
mapping node representations to a graph representation

$h_{G}=r(\\{h^{T}_{u}|u\in V\\})$ (9)

where $T$ is the total number of hidden layers.

Some popular GNNs include GCN [37], GraphSAGE [66], GIN [55], GAT [38], etc.
GIN first simply sums the aggregated messages over nodes’ neighborhood, and
then takes Multiple Layer Perceptrons (MLP) as the the update function. Its
forward functions under the formulation of Equation (8) are:

$\displaystyle m^{(t+1)}_{u}$ $\displaystyle=\sum_{v\in N(u)}h_{v}^{(t)}$ (10)
$\displaystyle h^{(t+1)}_{u}$
$\displaystyle=\mathrm{MLP}^{(t)}((1+{\epsilon}^{(t)})h_{u}^{(t)}+m_{u}^{(t+1)})$

GAT takes self-attention mechanism to re-weight the messages passed across
nodes, and takes a 1-layer MLP to update the aggregated messages. Its message
passing functions are:

$\displaystyle\alpha_{uv}^{(t+1)}$
$\displaystyle=\frac{\mathrm{exp}(\mathrm{LeakyReLU}(a^{T}[W^{(t)}h_{u}^{(t)}||W^{(t)}h_{v}^{(t)}]))}{\sum_{v^{\prime}\in
N(u)}\mathrm{exp}(\mathrm{LeakyReLU}(a^{T}[W^{(t)}h_{u}^{(t)}||W^{(t)}h_{v^{\prime}}^{(t)}]))};$
(11) $\displaystyle h_{u}^{(t+1)}$ $\displaystyle=\mathrm{ReLU}(\sum_{v\in
N(u)}\alpha_{uv}^{(t+1)}W^{(t)}h_{v}^{(t)}).$

Here, $\alpha_{uv}$ is the attention between nodes $u$ and $v$, $a$ and $W$
are parameters learned from the training process, $(\cdot||\cdot)$ denotes
vector concatenation operation, and the two non-linear functions processing
features are:

$\displaystyle\mathrm{LeakyReLU}(x)=\begin{cases}x&\text{x > 0}\\\ \lambda
x&\text{x $\leq$
0}\end{cases}~{}~{}\text{and}~{}~{}\mathrm{ReLU}(x)=\begin{cases}x&\text{x >
0}\\\ 0&\text{x $\leq$ 0}\end{cases}$ (12)

Our HS-GNN takes GIN and GAT layers in our experiments.

##  Appendix C More on HS-GNN

###  C.1 $\delta$-Net Clustering

See Algorithm 1 for the detailed $\delta$-Net clustering algorithm used to
find hierarchies in CNT bundles for HS-GNN.

Algorithm 1 $\delta$-Net clustering

0: A set of points $P=\\{p_{1},p_{2},...,p_{n}\\}$, radius $\delta>0$

0: A $\delta$-net $Q$ of $P$, and a set of clusters $\Pi$ where each point
$q\in Q$ corresponding to a cluster in $\Pi$

1: $P^{\prime}=\\{\\}$, $\Pi=\\{\\}$

2: while $P^{\prime}\neq$ do

3: Randomly pick point $p_{i}\in P\setminus P^{\prime}$

4: Obtain a cluster $C=\\{p_{j}:||p_{i}-p_{j}||_{2}\leq\delta|p_{j}\in P\\}$

5: $P^{\prime}=P^{\prime}\cup C$, $\Pi=\Pi\cup\\{C\\}$

6: end while

7: Set a node set $V=\\{\\}$ and an edge set $E=\\{\\}$

8: for $C_{i}\in\Pi$ do

9: Create a node $v_{i}$, add it to $V$

10: end for

11: for $(v_{i},v_{j})\in V\times V$ do

12: if $i\neq j$ and $C_{i},C_{j}$ have overlapping points in P then

13: Add edge $\\{v_{i},v_{j}\\}$ to $E$

14: end if

15: end for

16: Construct net $Q$ with node set $V$ and edge set $E$

###  C.2 Reweighting factors

We now describe how to compute the reweighting factor appearing in the message
passing process (4), $\tau_{uv}^{l}$, for two supernodes $u,v\in G^{(k)}$.
This $\tau_{uv}^{l}$ is computed by a function taking local geometric
information around $u$ and $v$ as input.

More precisely, the reweighting factors are learned by MLP taking principal
and norm vectors mentioned in Section (2.2). In particular, recall given any
supernode $a$ in $V_{k}$ corresponds to a cluster $C_{a}$ (resp. $C_{a}$) of
nodes from $V_{k-1}$, let $n_{a}$ denote the approximated normal vector and
$\nu_{a}^{i}(i=1,2)$ denote two principal vectors at $a$ computed by
performing PCA for points in $C_{a}$.

Now to compute $\tau_{uv}$, we concatenate information computed from PCA as
$\psi_{uv}=n_{u}||\nu_{u}^{1}||\nu_{u}^{2}||n_{v}||\nu_{v}^{1}||\nu_{v}^{2}$.
The rewieghting factor in the $l$-th convolution layer is:

$\tau^{l}_{uv}=S^{l}(f^{l}(\psi_{uv}))=\frac{e^{f^{l}(\psi_{uv})}}{\sum_{v^{\prime}\in
N(u)}e^{f^{l}(\psi_{uv^{\prime}})}}$ (13)

$f^{l}$ is a 3-layer MLP and $S^{l}$ is a softmax function.