# Efficient Parallel Multi-Hop Reasoning: A Scalable Approach for Knowledge
Graph Analysis

Jesmin Jahan Tithi  Parallel Computing Labs, Intel CorporationSanta
ClaraCAUSA95054<EMAIL_ADDRESS>,  Fabio Checconi  Parallel
Computing Labs, Intel CorporationSanta ClaraCAUSA95054
<EMAIL_ADDRESS>and  Fabrizio Petrini  Parallel Computing Labs, Intel
CorporationSanta ClaraCAUSA95054<EMAIL_ADDRESS>
(2024)

###### Abstract.

Multi-hop reasoning (MHR) is a process in artificial intelligence and natural
language processing where a system needs to make multiple inferential steps to
arrive at a conclusion or answer. In the context of knowledge graphs or
databases, it involves traversing multiple linked entities and relationships
to understand complex queries or perform tasks requiring a deeper
understanding. Multi-hop reasoning is a critical function in various
applications, including question answering, knowledge base completion, and
link prediction. It has garnered significant interest in artificial
intelligence, machine learning, and graph analytics.

This paper focuses on optimizing MHR for time efficiency on large-scale
graphs, diverging from the traditional emphasis on accuracy which is an
orthogonal goal. We introduce a novel parallel algorithm that harnesses
domain-specific learned embeddings to efficiently identify the top K paths
between vertices in a knowledge graph to find the best answers to a three-hop
query. Our contributions are: (1) We present a new parallel algorithm to
enhance MHR performance, scalability and efficiency. (2) We demonstrate the
algorithm’s superior performance on leading-edge Intel and AMD architectures
through empirical results.

We showcase the algorithm’s practicality through a case study on identifying
academic affiliations of potential Turing Award laureates in Deep Learning,
highlighting its capability to handle intricate entity relationships. This
demonstrates the potential of our approach to enabling high-performance MHR,
useful to navigate the growing complexity of modern knowledge graphs.

Multihop reasoning, Knowledge graph, Knowledge Graph Reasoning, Complex Query
Answering, Question Answering, Query answering

††copyright: acmlicensed††journalyear: 2024††doi: XXXXXXX.XXXXXXX††conference:
Irregular Applications; Nov, 2024; Atlanta, GA††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

Multi-hop reasoning (MHR) is a process in artificial intelligence and natural
language processing where a system needs to make multiple inferential steps to
arrive at a conclusion or answer. In the context of knowledge graphs or
databases, it involves traversing multiple linked entities and relationships
to understand complex queries or perform tasks that require a deeper level of
understanding.

For example, in a question-answering system, a single-hop query might ask,
”Where was Albert Einstein born?” The answer can be found directly by looking
up the birthplace of Albert Einstein in a knowledge graph. In contrast, a
multi-hop query might be, ”Which university did the physicist who proposed the
theory of relativity work for?” To answer this, the system must first identify
Albert Einstein as the physicist in question, then find the theory of
relativity associated with him, and finally, determine the university where he
worked.

Multi-hop reasoning is crucial for tasks that require connecting disparate
pieces of information that are not immediately adjacent in a knowledge graph.
It enables systems to handle more complex queries and provides a framework for
more sophisticated and contextually aware AI applications.

Encyclopedic KGs like Wikidata serve as repositories of structured world
knowledge, connecting entities through relational triples (head, relation,
tail) that are instrumental for downstream knowledge-driven applications.
Despite their breadth, these KGs are incomplete, often lacking critical
relational data. MHR emerges as a solution to infer these missing links,
thereby enriching the KG and enhancing its utility for applications that
demand a comprehensive knowledge base.

The essence of MHR lies in its ability to traverse expansive networks,
extracting salient insights by exploring multiple intermediate vertices or
paths. The MHR problem can be framed as the task of finding the ”best” $K$
paths between two vertices, $s$ and $t$, within a graph $G$ where vertices and
edges could be of different types and the edge weights can be computed based
on some learned embeddings (Guo and Berkhahn, 2016) of the vertices and edge
types reflecting domain-specific knowledge, guide the selection of paths that
reveal subtle and often hidden relationships. These embeddings can be learned
from a prior ML model such as a Graph Neural Network (Gilmer et al., 2020;
Khemani et al., 2024) or BERT (Reimers and Gurevych, 2019) training.

While existing research predominantly concentrates on enhancing the accuracy
of MHR processes, we assert the equal importance of performance optimization.
An unoptimized algorithm can be significantly slower, which is untenable for
time-sensitive applications such as question answering systems in Healthcare
Decision Support Systems, Financial Fraud Detection, Real-Time Recommendation
Engines, Interactive Voice Assistants, Emergency Response Systems, High-
Frequency Trading (HFT), etc.

Obtaining high performance on MHR is challenging since it involves large graph
or database traversals and hashtable accesses which requires extensive random
memory accesses and global synchronizations, which can pose significant
barriers to achieving high performance. Scaling MHR is also challenging due to
the large size of the input graphs and the associated learned embedding data.

Embedding vector tables can reach terabyte-scale dimensions, which could
exceed the memory capacity of a single server node.

A practical example of the MHR problem is the identification of academic
affiliations of potential Turing Award laureates within the Deep Learning (DL)
discipline (see Figure 1). This problem is an example of complex logical query
answering in a neural graph database and can be mapped to a three-hop
reasoning problem as discussed later in the paper. We introduce an efficient
parallel algorithm to address this MHR task. Our proposed model transcends
academic applications, providing valuable insights into complex networks where
MHR can uncover missing links, patterns, and predictions of substantial
importance.

Figure 1. An example snippet of Knowledge Graph (KG) and the composition of a
multi-hop reasoning query on the KG.

This work presents an efficient parallel algorithm that speeds up the MHR
problem. We begin with a basic parallel algorithm that a data scientist would
typically use to solve the problem. We explain the potential performance
improvement ideas and then propose an optimized algorithm to solve the same
problem significantly faster compared to the basic one on any machine. Our
approach uses concurrent hashmaps, k-heaps, and tree-based reduction to merge
k-heaps to improve computational throughput. The primary contribution of our
work is this optimized algorithm.

The contributions of this paper are threefolds:

  * •

We frame the MHR problem as determining the best $K$ paths between two
vertices in a graph, informed by a learned model. Given an input graph $G$
with entities $h,r$ connected with different types of edges/relationships $r$,
and given the learned embedding vectors for all the entities and relationships
in the graph, we try to answer, what are the $K$ best paths from $h$ to $r$.
This framework can be generalized to diverse domain-specific applications, can
be used to answer questions or suggest missing links in knowledge graphs.

  * •

We present a simple parallel algorithm for MHR that would be easy to implement
but may not be performant. Then introduce an improved algorithm that utilizes
concurrent hash tables, thread private k-heaps and tree merging of those heaps
that would result in significantly faster solutions.

  * •

We show performance analysis of these algorithms on cutting-edge AMD EPYC and
Intel SPR architectures, showcasing the practical effectiveness of our method.

Through these contributions, we aim to set a new precedent for performance-
oriented MHR, enabling faster and more efficient reasoning over complex
knowledge graphs.

##  2\. Related Work

Multi-hop reasoning (MHR) is a concept in natural language processing (NLP)
and knowledge graph (KG) analysis that involves answering complex queries by
traversing multiple connected entities in a sequence. Formally, multi-hop
queries, also known as path queries, are conjunctive queries (CQ) that form a
linear chain. In such a chain, the tail of one projection serves as the head
of the subsequent projection. Path queries can be resolved iteratively by
retrieving the neighbors of the current node. Multi-hop queries can also be
extended to include negation, allowing for more complex query formulations
(Hongyu et al., 2023).

MHR approaches often employ iterative methods like depth-first search (DFS) or
breadth-first search (BFS) to navigate paths within a knowledge graph (KG)
(Lao and Cohen, 2010). These methods may not consider the inherent meaning of
the entity and relationships and regardless can become computationally
intensive as the size of $K$ and the complexity of KGs increase.

Graph neural networks (GNNs) (Gilmer et al., 2020; Khemani et al., 2024) have
revolutionized learning over graph-structured data, offering sophisticated
means to deduce relationships. These networks can be trained to generate
entity and relationship embeddings (Guo and Berkhahn, 2016), which are dense
vector representations of categorical variables in a reduced-dimensional
space, encapsulating the relationships between categories.

TransE (Bordes et al., 2013) embedding model posits that the embedding of a
tail entity should be near the sum of the head entity and relation embeddings
if the triplet (head, relation, tail) is valid. A number of query reasoning
approaches have been developed to address multi-hop path queries through the
use of sequence models that leverage chainable knowledge graph embeddings,
such as TransE as demonstrated by Guu et al. (Guu et al., 2015), or employ
recurrent neural networks like LSTM, as utilized by Das et al. (Das et al.,
2016). These efforts represent some of the initial attempts to apply
embeddings and neural network methodologies to solve multi-hop path queries.

The Hop, Union, Generate (HUG) (Zhao et al., 2023) framework introduces an
explainable MHR framework that uses a pre-trained model to score paths for
their relevance to the reasoning task. The Tree-of-Mixed-Thought (ToMT) (Hu et
al., 2023) model combines rapid, one-stop reasoning with iterative refinement
from a learned model, striving for a balance between efficiency and accuracy
in MHR tasks.

None of the prior work is focused on improving runtime performance of MHR on
an HPC system. Our work introduces an efficient parallel algorithm for MHR. It
concurrently traverses the relavant knowledge graph nodes, using TransE-style
relations to find the most promising next paths to follow. Our algorithm
accelerates the reasoning process and scales well with the number of cores on
modern servers.

##  3\. Proposed Algorithm

In this section, we discuss the problem and the algorithm, and provide steps
to optimize the algorithm on a shared memory system.

###  3.1. The Problem and Its Mapping to MHR

Given a knowledge graph, we are trying to answer the following question: who
are the top K persons that are Winners of Turing Award, and also work in the
deep learning field and then find the top K universities these top K persons
are likely to be affiliated with. This problem can be rephrased as: 1) find
top K $X$ entities, who are most likely to be related to $A$ with a relation
$R1$, 2) out of those Top K, find those who are likely to be related to $B$
with relation $R2$ and 3) then find top K $C$s entities that are related to
those top $X$ with some other relationship $R3$. We can solve this problem
using a three-hop reasoning process - a particular type of Multi-hop Reasoning
(MHR). Note that this process can be extended to other problems without losing
generality.

For this particular problem, the first hop involves identifying nodes in the
knowledge graph that are award winners and has a connection to the ”TURING
AWARD” entity. We score these nodes based on the TransE relation (h + r - t),
where ”h” represents the TURING AWARD embedding, ”r” represents the award-
winning relationship, and ”t” represents the embedding for each person on the
list. Based on these scores, we select the top individuals who are likely to
receive the Turing Award.

In the second hop, we determine which of these individuals work in the field
of deep learning. We use the TransE relationship again, where ”h” represents
”DEEP LEARNING”, ”r” represents the ”works in” relationship, and ”t”
represents the person of interest. We rank them based on the (h + r - t)
TransE score.

In the final hop, we aim to determine the most probable university
affiliations of these individuals. In this case, we use the ”affiliated with”
relationship, where ”h” represents the person, and ”t” represents a
university.

Our problem requires three different inputs. Firstly, we have a knowledge
graph represented as edge lists in the format (h, r, t). Here, ”h” represents
the source node, ”t” represents the destination node, and ”r” represents the
relationship between them. These nodes represent different entities each with
their pre-learned embedding vectors. Secondly, we have a relationship
embedding table that contains learned meanings of each relationship from a
machine learning model. Finally, we have an embedding table that contains
embedding vectors for each entity, which were also pre-learned using a
machine-learning model.

###  3.2. Simple Algorithm

Listing 1 shows a simple algorithm to solve the multi-hop reasoning problem to
discover potential academic affiliations for Turing Award prospects in the
field of deep learning, utilizing pre-learned entity and relation embeddings.
Here are the steps in the algorithm:

Initialization of Graph Data: The algorithm begins by parsing triples (head
entity $h$, relation $r$, tail entity $t$) from an edge list file $E$. These
triples represent the edges of the knowledge graph and are used to initialize
the graph structure. Edge tables $E_{r}$ are populated in parallel for each
relation $r$, indexing the head entities $h$ as keys. This step is
synchronized using locks to prevent concurrent write conflicts.

Loading Embeddings: Entity embeddings are loaded from a file V into an entity
embedding hash table $V_{e}$, with each entity uniquely identified by its ID.
This operation is also performed in parallel with locks ensuring data
integrity. Similarly, relation embeddings are loaded from a file R into a
relation embedding hash table $R_{e}$, indexed by the relation type.

Extraction of Entities and Universities: Person entities are extracted from
the edge table $E_{AWARDWINNER}$ and stored in a hashtable or array P. This is
done in parallel using locks. A set of universities U is extracted from the
edge table $E_{AFFILIATEDWITH}$ in parallel, forming the university table.

Scoring and Ranking Process: A parameter K is set to define the number of top
paths to be considered. A ranked list of persons $L_{p}$ is initialized as
empty. In parallel, each person p related to the ”Turing Award” is scored
using the TransE model, and the pair (person p, score $S_{p}$) is added to the
list $L_{p}$ atomically. A TopK heap $H_{p}$ is created from the list $L_{p}$
based on the scores, and the heap is sorted in descending order to produce a
ranked list of persons $TP$.

Deep Learning Association: For each top-ranked person ($p_{k}$, $S_{p_{k}}$)
in $TP$, a score $S_{p_{k},DL}$ is computed for their association with ”deep
learning,” and the pair is added to the list $L_{p}$ using locks. A TopK heap
$H_{p}$ is created from the list $L_{p}$ based on the scores, and the heap is
sorted in descending order in $TP$.

University Affiliation Scoring: The process is repeated for each person in
$TP$ to compute affiliation scores $S_{p_{k},u}$ for each university $u$ in
the set U. These scores are added to the university list $L_{u}$ using locks.
A TopK heap $H_{u}$ is created from the list $L_{u}$ using the affiliation
scores, and the heap is sorted in descending order to rank the universities.
The algorithm concludes by returning a list of potential affiliations A for
the top-scored individuals, representing the final output of the process.

Algorithm 1 Parallel Multi-Hop Reasoning for Academic Affiliation Discovery

1:Edge list $E$, entity embeddings $V$, relation embeddings $R$

2:Top academic affiliations $A$ for Turing Award prospects in deep learning

3:Parallel Parse $(h,r,t)$ from $E$ to initialize graph grouped by relation
types $r$ $\triangleright$ Initialize graph

4:Parallel Populate $E_{r}$ using $h$ as key with locks $\triangleright$
Populate edge tables

5:Parallel Populate $V_{e}$ from $V$ using entity ID as key with locks
$\triangleright$ Load entity embeddings

6:Parallel Populate $R_{e}$ from $R$ using relation type as key with locks
$\triangleright$ Load relation embeddings

7:Parallel Extract $P$ from $E_{\text{AWARD WINNER}}$ with locks
$\triangleright$ Extract persons

8:Parallel Extract $U$ from $E_{\text{AFFILIATED WITH}}$ with locks
$\triangleright$ Extract universities

9:Set $K$ for the number of top paths $\triangleright$ Define top paths

10:Initialize $L_{p}\leftarrow\emptyset$ $\triangleright$ Initialize person
list

11:Parallel

12:for all $p\in P$ do

13: Compute $S_{p}$ using TransE model

14: Atomically insert $(p,S_{p})$ to $L_{p}$

15:end for

16:Create TopK heap $H_{p}$ from $L_{p}$ and sort it by score in
$TP$$\triangleright$ Rank persons

17:Parallel

18:for all $(p_{k},S_{p_{k}})\in TP$ do

19: Compute $S_{p_{k},\text{DL}}$ for ”deep learning” association

20: Add $(p_{k},S_{p_{k}})$ to $L_{p}$ with locks

21:end for

22:Create TopK heap $H_{p}$ from $L_{p}$ and sort it by score in
$TP$$\triangleright$ Rank persons

23:for all $(p_{k},S_{p_{k}})\in TP$ do

24: Parallel

25: for all $u\in U$ do

26: Compute $S_{p_{k},u}$ for ”AFFILIATED WITH”

27: Add $(u,S_{p_{k},u})$ to $L_{u}$ with locks

28: end for

29: Create TopK heap $H_{u}$ from $L_{u}$ and sort it by score in
$A$$\triangleright$ Rank universities

30: Print $A$ $\triangleright$ Output affiliations

31:end for

This could be a standard algorithm to solve the multi hop reasoning problem in
Python or C++. Although it is a good starting point, the algorithm has several
inadequacies and scalability issues. In this paper, we introduce an improved
algorithm that solves the same problem that can lead to a significantly lower
runtime cost.

###  3.3. Algorithmic Optimizations

For performance enhancement, we introduce a custom thread-safe hashtable
tailored for high-concurrency insertions. This approach mitigates the
inefficiencies inherent in standard STL hashtables that rely on user-level
locking. Our specialized hashtable is designed to handle the storage demands
of embedding vectors, which may scale to terabytes, thereby reducing
synchronization overhead and increasing the throughput of our MHR algorithm.
While hardware-accelerated hashing presents a further opportunity for
optimization, it is beyond the scope of this paper.

In the context of scoring and ranking individuals, we address the issue of
thread conflicts in shared hashmaps by utilizing thread-private k-heaps for
storing (person, score) pairs. These k-heaps, being priority queues, maintain
only the top K scores, offering several advantages:

  * •

Synchronization across threads during insertions is eliminated, enhancing
scalability.

  * •

Storage complexity is reduced from O(Persons) to O(threads * K), which is more
practical for large datasets.

  * •

K-heaps typically reside in thread-local storage, such as cache, leading to
faster insertions compared to global hashmaps that require locks and are
stored in main memory.

  * •

By tracking only the top K elements per thread, we conserve memory and
simplify the identification of the global top K elements in subsequent steps.

function ReduceTopKTree($localTopK$, $globalTopK$, $numThreads$, $threadId$)

Parallel

for $stride\leftarrow 1$ to $numThreads$ by $stride\times 2$ do

if $threadId$ is multiple of $2\times stride$ and $threadId+stride<numThreads$
then

$otherThreadId\leftarrow threadId+stride$

Merge heaps $localTopK[threadId]$ and $localTopK[otherThreadId]$

end if

Perform barrier synchronization

end for

single thread execution:

Update $globalTopK$ with $localTopK[0]$

end function

Algorithm 2 Tree-based reduction Algorithm to merge K-heaps from all threads
into one single K-heap.

Upon completion of individual scoring by all threads, we merge the thread-
local k-heaps into a unified k-heap to determine the global top K individuals.
For this purpose, we use a tree-based heap merging algorithm, ReduceTopKTree,
which employs a binary tree structure for efficient heap consolidation. The
pseudocode for this algorithm is detailed in Algorithm 2. This can be further
improved by using a quaternary tree merging. Note that, one could use locks to
merge the K-heaps. However, based on our experiences, we anticipate that such
approch would be less scalable on parallel machines with many threads.

The computeScorePerPerson algorithm, detailed in Algorithm 3, systematically
identifies individuals most likely associated with the Turing Award by
leveraging their entity and relation embeddings. The process unfolds as
follows:

Initially, the algorithm retrieves the Turing Award’s entity and relation
embeddings from their respective tables. These are then combined into a
composite embedding, encapsulating the award’s characteristics within the
knowledge graph.

A normalization constant, gamma, is set to 1.0, and the algorithm is
configured for parallel execution by initializing a local top-k heap for each
thread. This setup facilitates the simultaneous processing of individual
scores.

During the parallel computation phase, threads independently compute
similarity scores for individuals by contrasting their embeddings with the
composite Turing Award embedding. The score reflects the degree of association
with the award, calculated as the normalized difference between the
embeddings.

Each thread maintains a local top-k heap to store the highest scores, ensuring
that only the most pertinent candidates are considered. Upon completion, these
local heaps are consolidated into a global top-k heap through the reduction
function descrived before.

The outcome is a global top-k heap that ranks individuals according to their
association with the Turing Award, as determined by their embeddings. This
parallel algorithm enables scalable and efficient identification of the most
relevant candidates.

Algorithm 3 Compute Top-k Persons related to the Turing Award

1:$EntityEmbeddingTable$, $RelationEmbeddingTable$

2:$top\\_k\\_elements$ (MaxHeap)

3:function computeScorePerPerson($EntityEmbeddingTable$,
$RelationEmbeddingTable$, $top\\_k\\_elements$)

4: $turing\\_emb\\_ptr\leftarrow EntityEmbeddingTable.find(TURING\\_AWARD)$

5: $TURING\\_AWARD\\_EMB\leftarrow turing\\_emb\\_ptr.embedding$

6: $AWARD\\_WINNER\\_TABLE\\_EMB\leftarrow
RelationEmbeddingTable[AWARD\\_WINNER\\_TABLE\\_IDX].embedding$

7: $turing\\_arr\leftarrow$ embeddingAggregation($TURING\\_AWARD\\_EMB$,
$AWARD\\_WINNER\\_TABLE\\_EMB$)

8: $gamma\leftarrow 1.0$

9: $num\\_persons\leftarrow$ size of $person\\_ids$

10: $localTopK[num\\_threads]$ $\triangleright$ Array of MaxHeaps, one per
thread Parallel for: $id\leftarrow 0$ to $num\\_threads-1$

11: $localTopK\\_t\leftarrow localTopK[id]$

12: for $i\leftarrow 0$ to $num\\_persons-1$ do

13: $person\\_id\leftarrow person\\_ids[i]$

14: $record\leftarrow EntityScoreInfo(person\\_id,-\infty)$

15: $e,emb\leftarrow EntityEmbeddingTable.find(person\\_id)$

16: if $e\neq HashmapStatus.Fail$ then

17: $total\\_score\leftarrow 0.0$

18: for $j\leftarrow 0$ to $EMBEDDING\\_DIMENSION-1$ do

19: $total\\_score\leftarrow
total\\_score+|\text{turing\\_arr[j]}-emb.embedding[j]|$

20: end for

21: $record.score\leftarrow gamma-total\\_score$

22: end if

23: $localTopK\\_t.push(record)$

24: if $localTopK\\_t.size()>TOP\\_K$ then

25: $localTopK\\_t.pop()$

26: end if

27: end for

28: reduceTopKTree($localTopK$, $top\\_k\\_elements$, $num\\_threads$, $id$)
End Parallel for

29:end function

The next step in the process is to rank or rearrange these top-K Turing Award
winning people based on their probability of working in the Deep Learning
field. The computeScoreBasedOnWorksInDL function calculates a score for
individuals potentially associated with the Turing Award and active in the
Deep Learning field. It accepts entity and relation embeddings, a max heap for
top-k elements, and identifiers for the head entity, relation, and individual.

The algorithm retrieves the relevant embeddings for the head entity and
relation, then computes the L1 norm as the sum of absolute differences between
the head entity plus relation embeddings and the individual’s embedding. This
score reflects the individual’s alignment with the characteristics of the
Turing Award and their contributions to Deep Learning.

An adjustment is made to the score by subtracting it from a constant, gamma.
The resulting TransE score is used to rank the individual’s relevance. The
function ensures thread-safe operations while updating the max heap with the
individual’s score, maintaining only the top-k scores. This process identifies
the most relevant candidates for the award within the Deep Learning context.

Again, in this case, the simple algorithm might use a global hashmap to keep
the key (person) and score (association probability based on TransE) and
instead, we use K-heaps.

The next step of the MHR process is to find the most likely university
affiliations of these top-k people found in the prior step. The function
computeAffiliationScore calculates the top-k affiliations for a given a
person, based on their association with various universities. The function
takes the entity and relation embedding tables, a max heap to store the top-k
elements, and identifiers for the head entity and the relation of interest.

First, the function retrieves the embeddings for the head entity (person) and
the specified relation (AFFILIATED_WITH). These embeddings are aggregated to
form a composite embedding array, affiliated_arr, which represents the
combined characteristics of the head entity and the relation.

The algorithm sets a constant, gamma, to 1.0, which will be used to adjust the
final scores. In a parallelized block, each thread processes a subset of
universities. The threads each maintain a local max heap, localTopK_t, to
store the top-k affiliation scores for the universities they process. For each
university, the function computes a score by summing the absolute differences
between the composite embedding array and the university’s embedding. This
score is adjusted by subtracting it from gamma, yielding the final affiliation
score for the university.

The local top-k heaps are managed to ensure they contain the top-k scores. If
a heap exceeds the top-k limit, the lowest scoring element is removed. After
processing, the local heaps are merged into the global top-k heap using a
reduction function, reduceTopKTree.

The final output is the global top-k heap, which contains the top-k university
affiliations for the head entity based on their embeddings. This heap
represents the most relevant affiliations, identifying the universities most
closely associated with the person of interest. Algorithm 4 shows a pseudocode
for the process.

Algorithm 4 Compute Top-k Affiliations for a Person

1:$EntityEmbeddingTable$, $RelationEmbeddingTable$, $headId$, $relationId$

2:$top\\_k\\_elements$ (MaxHeap)

3:function computeAffiliationScore($EntityEmbeddingTable$,
$RelationEmbeddingTable$, $top\\_k\\_elements$, $headId$, $relationId$)

4: $head\\_emb\\_ptr\leftarrow EntityEmbeddingTable.find(headId)$

5: $head\\_emb\leftarrow head\\_emb\\_ptr.embedding$

6: $relation\\_emb\leftarrow RelationEmbeddingTable[relationId].embedding$

7: $affiliated\\_arr\leftarrow$ embeddingAggregation($head\\_emb$,
$relation\\_emb$)

8: $gamma\leftarrow 1.0$

9: $num\\_universitys\leftarrow$ size of $university\\_ids$

10: $localTopK[num\\_threads]$ $\triangleright$ Array of MaxHeaps, one per
thread Parallel for: $id\leftarrow 0$ to $num\\_threads-1$

11: $localTopK\\_t\leftarrow localTopK[id]$

12: for $i\leftarrow 0$ to $num\\_universitys-1$ do

13: $university\\_id\leftarrow university\\_ids[i]$

14: $record\leftarrow EntityScoreInfo(university\\_id,-\infty)$

15: $e,uni\\_emb\\_ptr\leftarrow EntityEmbeddingTable.find(university\\_id)$

16: if $e\neq HashmapStatus.Fail$ then

17: $total\\_score\leftarrow 0.0$

18: for $j\leftarrow 0$ to $EMBEDDING\\_DIMENSION-1$ do

19: $total\\_score\leftarrow
total\\_score+|\text{affiliated\\_arr[j]}-uni\\_emb\\_ptr.embedding[j]|$

20: end for

21: $record.score\leftarrow gamma-total\\_score$

22: end if

23: $localTopK\\_t.push(record)$

24: if $localTopK\\_t.size()>TOP\\_K$ then

25: $localTopK\\_t.pop()$

26: end if

27: end for

28: reduceTopKTree($localTopK$, $top\\_k\\_elements$, $num\\_threads$, $id$)
End Parallel for

29:end function

##  4\. Performance Result

###  4.1. Dataset

To demonstrate the performance improvement of our new algorithm, we used the
WikiKG90Mv2(Hu et al., 2021) dataset from the Stanford Open Graph Benchmark
(OGB) dataset (ogb.stanford.edu). The WikiKG90Mv2 dataset is a subset of
Wikidata, aimed at enhancing the graph’s coverage by predicting missing
triples. It contains 91,230,610 entities, 1,387 relations, and 601,062,811
triples (h, r, t), extracted from Wikidata dumps as of May 17th, June 7th, and
June 28th, 2021, for training, validation, and testing. Each entity and
relation is accompanied by text features, including titles and descriptions,
with MPNet(Song et al., 2020) sentence embeddings (Reimers and Gurevych, 2019)
provided to capture semantic meanings. This dataset is utilized to assess the
performance of our multi-hop reasoning algorithm. However, the dataset has
been trimmed (embedding size = 8) so that we can fit them in a single
machine’s memory.

The process of constructing a knowledge graph begins with reading an
WikiKG90Mv2 edgelist file, where each line represents a triple consisting of a
head, a relation, and a tail. These triples are then sorted into 1,387
distinct relation tables, with each (head, tail) pair being placed in the
table corresponding to its relation. Following this, two separate files
containing embeddings for approximately 91 million entities and all relations
are read to construct the entity and relation embedding tables, respectively.

A specialized table is created to store person entities, initially populated
by examining the ”AWARD WINNER” relation table. Similarly, a
university/institution table is assembled, primarily using the ”AFFILIATED
WITH” relation table. Table 1 shows different statistics about the data.

Table 1. Property of the Dataset used to test MHR Parameter for Edge Data
---
Num_edges | 6.01E+08 |
File size (B) | 1.44E+11 |
Record size (B) | 256 |
LINE SIZE | 256 |
Parameter for EntityEmbedding
num_entity | 91000000 |
File Size | 4.168E+12 |
Record Size | 3072 |
EMD size | 768 |
LINE SIZE | 3072 |
Parameter for RelationEmbedding
Num_relations | 1387 |
File Size | 6.48E+10 |
Record Size | 3072 |
EMD size | 768 |
LINE SIZE | 3072 |
Parameter for PersonTable
Num_person | 539603 |
person_record_size | 16 |
Parameter for UniversityTable
Num_university | 89689 |
University record size | 16 |
Parameter for TopK
Top_k | 50 |
Multihop Reasoning
EMB Used | 8 |
Elements | 6.92E+08 |

###  4.2. Systems Specifications

Table 2 presents the system parameters of the Intel and AMD platforms we used
to evaluate the performance of our algorithm implementations. Both systems
have similar capabilities, but the EPYC system has a higher average clock
speed and slightly higher core counts while having slightly less bandwidth
(about 30%). Additionally, the EPYC system has 8 NUMA dies, whereas SPR has 2.
Based on these specifications, we expect SPR to be faster for the bandwidth-
bound portions of the code, while EPYC could better in handling latency due to
its more cores. However, the EPYC system might be slightly more affected by
NUMAness than SPR. Overall, we anticipate that these platforms will perform
similarly. The EPYC system belongs to the Perlmutter Supercomputing cluster
hosted at Berkeley, while SPR is an in-house system. We collected the
bandwidth information using Memory Latency Checker which runs on Intel systems
but does not run on EPYC. For EPYC, we used the stream triad benchmark to
measure the bandwidth.

We utilized OpenMP and C/C++ languages to develop parallel implementations for
both the simple and the optimized algorithms. Afterwards, we compiled the
program using g++ compiler while applying standard optimization flags like
”-O3 -fopenmp -Ofast” on both platforms.

Table 2. Systems Specifications Used to Test The Algorithms System Parameters | SPR | EPYC
---|---|---
#cores/sockets | 56 | 64
Model | Intel(R) Xeon(R) Platinum 8480L | AMD EPYC 7763 64-Core Processor
#threads | 112 | 128
#hyperthreads | 224 | 256
#sockets | 2 | 2
Clock Speed (GHz) | 2 | 2.45
Memory Size (GB) | 527 total, 257 per socket | 527 total, 65 per socket
#Memory Nodes/Dies | 2 | 8
Total Memory BW (GB/s) | 527 (all read), 474 (stream triad) | 379 (stream triad)
Cache Size | 107 MB/core | 512 KB/core
Compiler Version (g++) | 11.4.1 20230605 | 12.3.0

###  4.3. Performance Trends

In this section, we discuss the performance results of the MHR code on the SPR
platform and compare it with the optimized version’s performance on both SPR
and EPYC platforms. We utilized the command ”KMP_AFFINITY=compact
OMP_PLACES=cores numactl -a -b -i 0” to affinitize the threads to cores, which
helped in creating a balanced and interleaved memory allocation across the
memory nodes. This setting was crucial for achieving satisfactory performance
on the EPYC machine.

Figure 2 illustrates the difference in the performance of the simple algorithm
and the optimized algorithm running on a single socket of the SPR node with 56
cores. The results show that the optimized algorithm is 100 times faster than
the simple algorithm when executed on a single socket of the SPR node. The
speedup difference between 28 cores and 56 cores indicates that the optimized
algorithm is more scalable than the simple algorithm. This is because the
speedup increases as the number of cores increases. The optimized algorithm
uses a customized hashtable with fine-grained locking and private K-heaps
instead of a global hashtable, as explained in the previous sections all of
which contribute to such a speedup.

Figure 2. Performance Benefit of the old algorithm compared on SPR Platform.

The graph depicted in Figure 3 demonstrates the strong scaling of an optimized
algorithm on the SPR platform. The multHopReasoning task represents the entire
process, while computeScorePerPerson, computeScoreBasedOnWorksInDL, and
computeAffiliationScore are the three sub-tasks of multHopReasoning, as
described in the algorithm description. One important point to note is that
among these three sub-tasks, computeAffiliationScore takes the most time since
it is executed O(Persons * Universities) times, whereas the other two sub-
tasks are executed O(Persons) and O(top-K) times, respectively. The graph in
Figure 3 shows that the algorithm performs well at different core counts and
scales well on 56 cores of one socket of SPR. The algorithm’s performance
improves by approximately 7 times on 56 cores. It is worth mentioning that the
MHR algorithm is very sparse and involves frequent hash lookup mixed with
embedding gathering, making it a challenging workload to obtain high
performance.

(a) Runtime scaling

(b) Speedup scaling

Figure 3. Strong Scaling of Runtime of the MHR Algorithm on SPR Platform.

Figure 4 illustrates the optimized algorithm’s performance on the EPYC
platform across varying core counts. The algorithm exhibits robust scaling up
to the utilization of one socket, equivalent to four NUMA nodes. However,
beyond three NUMA nodes, no further performance gains are noted. Notably, the
algorithm achieves a 5.6x speedup on 64 cores and a 7.2x speedup on 32 cores.
Critical to achieving this scaling is the use of thread pinning and
interleaved memory allocation, as specified by ”OMP_PLACES=cores numactl -a -b
-I 0-3”. Without such measures, performance is markedly reduced, underscoring
the importance of data interleaving in mitigating NUMA-related slowdowns.
Future work will focus on addressing NUMA latencies to further enhance
algorithm performance on NUMA architectures.

(a) Runtime scaling

(b) Speedup scaling

Figure 4. Strong Scaling of Runtime of the MHR Algorithm on EPYC Platform.

In evaluating the performance disparity between SPR and EPYC platforms, Figure
5 reveals that SPR surpasses EPYC by 40% when operating on a single socket
configuration. This performance advantage is primarily due to SPR’s superior
memory bandwidth and EPYC’s fasters clock speed. The algorithm experiences a
greater impact from the Non-Uniform Memory Access (NUMA) effect on the EPYC
system compared to SPR, indicating that SPR’s architecture mitigates NUMA-
related performance degradation more effectively.

Figure 5. Performance difference between SPR and EPYC.

Finally, we also experimented with a version of our optimized algorithm that
uses locks to merge in the global top-K heap instead of using the tree-based
reduction and as expected this version performed slightly worse than the tree
merge algorithm.

##  5\. Generalization of the MHR Algorithm

The versatility of our MHR algorithm extends beyond the specific use case of
identifying academic affiliations for Turing Award candidates. Its underlying
structure is designed to generalize across various domains where multi-hop
queries are essential. By abstracting the core components of the algorithm, we
can apply it to different types of entities and relationships within any
knowledge graph.

The generalization process involves redefining the entities and relationships
of interest and then applying the same reasoning steps. For instance, the
algorithm can be adapted to identify connections between medical conditions
and treatments, or to trace supply chain networks from manufacturers to end
consumers.

Below is a high-level pseudocode that illustrates the generalization of the
algorithm.

Algorithm 5 Generalized Three-Hop Reasoning Algorithm

1:GeneralEntityEmbeddingTable, GeneralRelationEmbeddingTable

2:GeneralTopKResults (MaxHeap)

3:function GeneralizedMultiHopReasoning(GeneralEntityEmbeddingTable,
GeneralRelationEmbeddingTable, GeneralTopKResults)

4: Define general entities and relationships of interest

5: GeneralEntityEmbeddingTable $\leftarrow$ LoadEntityEmbeddings()

6: GeneralRelationEmbeddingTable $\leftarrow$ LoadRelationEmbeddings()

7: FirstHopResults $\leftarrow$
computeRelation1ScorePerEntity1(GeneralEntityEmbeddingTable,
GeneralRelationEmbeddingTable, GeneralTopKResults)

8: SecondHopResults $\leftarrow$
computeScoreBasedOnRelation2Entity2(FirstHopResults, GeneralTopKResults)

9: FinalResults $\leftarrow$ computeRelation3ScorePerEntity1(SecondHopResults,
GeneralTopKResults)

10: GeneralTopKResults $\leftarrow$ MergeTopKHeaps(FinalResults)

11: return GeneralTopKResults

12:end function

In this pseudocode, LoadEntityEmbeddings() and LoadRelationEmbeddings() are
placeholders for functions that would load the relevant embeddings based on
the domain-specific problem. The computeRelation1ScorePerEntity1(maps to
computeScorePerPerson), computeScoreBasedOnRelation2Entity2
(maps to computeScoreBasedOnWorksInDL), and
computeRelation3ScorePerEntity1(maps to computeAffiliationScorePerLocal)
functions are reused from the Turing Award problem, demonstrating the
algorithm’s adaptability. The MergeTopKHeaps function represents the final
step of consolidating results from multiple threads.

Algorithm 6 Generic Multi-Hop Reasoning Algorithm

1:Graph $G$, Embeddings $E$, Source $S$, Target $T$, Number of Hops
$NUM\\_HOPS$, Top $K$

2:Results heap $results$

3:function Multihop_reasoning_generic($G$, $E$, $S$, $T$, $NUM\\_HOPS$,
$topK$, $results$)

4: $path\leftarrow$ vector of int

5: $path.insert(S)$ $\triangleright$ Initialize path with source vertex

6: $total\\_paths\leftarrow(\textsc{pow}(topK,NUM\\_HOPS-1)-1)/(topK-1)$

7: Define $paths[total\\_paths\times 2]$ $\triangleright$ Array of paths

8: $curr,next\leftarrow$ pointers to path arrays

9: $curr\\_size\leftarrow 0$

10: $next\\_size\leftarrow 0$

11: $curr[curr\\_size++]\leftarrow path$ $\triangleright$ Start with initial
path

12: for $level\leftarrow 0$ to $NUM\\_HOPS-1$ do

13: #pragma omp parallel for num_threads(min($curr\\_size$,
$total\\_threads$))

14: for $p\leftarrow 0$ to $curr\\_size-1$ do

15: expand_path($curr[p]$, $next$, $next\\_size$, $G$, $E$, $T$, $topK$,
$results$)

16: end for

17: swap($curr$, $next$)

18: $curr\\_size\leftarrow next\\_size$

19: $next\\_size\leftarrow 0$

20: end for

21:end function

Algorithm 7 Expand Path in Multi-Hop Reasoning

1:Current path $path$, Next paths array $next$, Size of next paths array
$next\\_size$, Graph $G$, Embeddings $E$, Target node $T$, Number of hops
$NUM\\_HOPS$, Top $K$ paths, Results heap $results$

2:Updated next paths array $next$, Updated size $next\\_size$, Updated results
heap $results$

3:function expand_path($path$, $next$, $next\\_size$, $G$, $E$, $T$,
$NUM\\_HOPS$, $topK$, $results$)

4: $horizon\\_node\leftarrow path[\text{end}]$

5: $neighbors\leftarrow G[horizon\\_node].neighbors$

6: if $neighbors$ is not empty then

7: $path\\_emb\leftarrow$ EmbeddingForPath($path$, $E$)

8: $score\\_heap\leftarrow$ new max heap of size $topK$

9: Parallel

10: for $i\leftarrow 0$ to $|neighbors|-1$ do

11: $neighbor\leftarrow neighbors[i]$

12: if not $neighbor$ in $path$ then

13: $extended\\_path\leftarrow path+[neighbor]$

14: $extended\\_emb\leftarrow$ ExtendEmbedding($path\\_emb$, $E[neighbor]$)

15: $score\leftarrow$ ComputeScore($extended\\_emb$)

16: if $neighbor=T$ then

17: Atomically $results$.insert($\\{extended\\_path,score\\}$)

18: else

19: Atomically $score\\_heap$.insert($\\{extended\\_path,score\\}$)

20: end if

21: end if

22: end for

23: Atomically UpdateNextPaths($score\\_heap$, $next$, $next\\_size$, $topK$)

24: end if

25:end function

Algorithms 6, 7 show another example of generic multi-hop reasoning algorithm
where the number of hops is also a parameter. The Multihop_reasoning_generic
algorithm is designed to perform multi-hop reasoning on a given knowledge
graph. It takes as input a graph $G$, a set of embeddings $E$, a source node
$S$, a target node $T$, the number of hops $NUM\\_HOPS$, and the number of top
paths $K$ to find. The algorithm outputs a heap of results that contains the
top paths from the source to the target.

The algorithm begins by initializing a path with the source node. It then
calculates the total number of paths to explore based on the number of hops
and the top $K$ paths. Two arrays of paths, curr and next, are defined to keep
track of the current and next set of paths to explore. The algorithm starts
with the initial path and iteratively expands each path in parallel using the
expand_path function.

The expand_path function takes a path and expands it by exploring its
neighboring nodes. For each neighbor, it checks if the neighbor is already in
the path to avoid cycles. If the neighbor is not in the path, it is added, and
the path embedding is extended. The score of the extended path is computed,
and if the neighbor is the target node, the path and its score are added to
the results heap. Otherwise, the path and its score are added to a local max
heap of scores. After exploring all neighbors, the top paths from the local
max heap are added to the next array of paths.

The algorithm continues this process for the specified number of hops,
swapping the curr and next arrays at each iteration. The final results heap
contains the top paths from the source to the target, as determined by the
computed scores.

This generalized approach allows the MHR algorithm to be applied to a wide
range of problems, making it a powerful tool for extracting insights from
complex, interconnected data.

The three-hop reasoning problem presented in this paper is a particular
instance of the generalized problem and the algorithmic optimizations
presented for the three-hop problem are applicable to the general problem as
well.

##  6\. Conclusions

In this study, we have presented novel parallel algorithms designed to
optimize multi-hop reasoning (MHR) tasks on large-scale knowledge graphs. Our
approach has been tested using the WikiKG90Mv2 dataset, a subset of Wikidata,
which serves as a challenging benchmark due to its extensive size and the rich
semantic information it encapsulates. Despite the necessity to trim the
dataset to fit within the memory constraints of a single machine, our
algorithm demonstrated significant performance improvements.

We evaluated the effectiveness of the algorithm on two high-performance
computing platforms, SPR and EPYC, each with unique architectural features.
Our findings show that the optimized algorithm achieves significant speed
improvements, surpassing the simple algorithm by a factor of 100 on a single
socket in the SPR system. This increase in efficiency is due to the
algorithm’s ability to better manage access to the hashmap, resulting in fewer
access conflicts, less synchronization overhead, and overall higher work
efficiency compared to the simple algorithm due to the usage of private
K-heaps and tree-based merging.

On the EPYC platform, the algorithm exhibited strong scaling up to the limits
of a single socket, with diminishing returns observed beyond 3 NUMA nodes.
Thread pinning and balanced memory allocation were crucial in realizing the
full potential of the EPYC system, underscoring the importance of memory
access patterns in optimizing MHR tasks.

Our findings underscore the significance of memory bandwidth and the impact of
NUMA effects on the performance of MHR algorithms. The SPR platform’s
architecture provided a notable advantage in memory-intensive operations,
leading to a 40% performance improvement over the EPYC system in single-socket
configurations.

The insights gained from this research contribute to the broader understanding
of MHR optimization strategies, highlighting the critical role of system
architecture in achieving efficient parallel processing. This study
contributes to the field of MHR by demonstrating the potential of algorithmic
optimization to leverage high-performance computing architectures effectively.
Our findings provide a foundation for future advancements in MHR algorithms,
aiming for enhanced scalability and efficiency in processing complex knowledge
graphs.

##  7\. Acknowledgement

This research is based upon work supported by the Office of the Director of
National Intelligence (ODNI), Intelligence Advanced Research Projects Activity
(IARPA), through the Advanced Graphical Intelligence Logical Computing
Environment (AGILE) research program, under Army Research Office (ARO)
contract number ¡W911NF22C0081¿. The views and conclusions contained herein
are those of the authors and should not be interpreted as necessarily
representing the official policies or endorsements, either expressed or
implied, of the ODNI, IARPA, ARO, or the U.S. Government.

## References

  * (1)
  * Bordes et al. (2013) Antoine Bordes, Nicolas Usunier, Alberto Garcia-Duran, Jason Weston, and Oksana Yakhnenko. 2013\.  Translating embeddings for modeling multi-relational data.  _Advances in neural information processing systems_ 26 (2013).
  * Das et al. (2016) Rajarshi Das, Arvind Neelakantan, David Belanger, and Andrew McCallum. 2016.  Chains of reasoning over entities, relations, and text using recurrent neural networks.  _arXiv preprint arXiv:1607.01426_ (2016).
  * Gilmer et al. (2020) Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. 2020\.  Message passing neural networks.  _Machine learning meets quantum physics_ (2020), 199–214.
  * Guo and Berkhahn (2016) Cheng Guo and Felix Berkhahn. 2016.  Entity embeddings of categorical variables.  _arXiv preprint arXiv:1604.06737_ (2016).
  * Guu et al. (2015) Kelvin Guu, John Miller, and Percy Liang. 2015.  Traversing knowledge graphs in vector space.  _arXiv preprint arXiv:1506.01094_ (2015).
  * Hongyu et al. (2023) Ren Hongyu, Galkin Mikhail, Cochez3 Michael, Zhu4 Zhaocheng, and Leskovec1 Jure. 2023\.  Neural Graph Reasoning: Complex Logical Query Answering Meets Graph Databases.  _arXiv preprint arXiv:2303.14617_ (2023).
  * Hu et al. (2023) Pengbo Hu, Ji Qi, Xingyu Li, Hong Li, Xinqi Wang, Bing Quan, Ruiyu Wang, and Yi Zhou. 2023\.  Tree-of-Mixed-Thought: Combining Fast and Slow Thinking for Multi-hop Visual Reasoning.  arXiv:2308.09658 [cs.CL]
  * Hu et al. (2021) Weihua Hu, Matthias Fey, Hongyu Ren, Maho Nakata, Yuxiao Dong, and Jure Leskovec. 2021\.  OGB Large-Scale Challenge – WikiKG90Mv2.  https://ogb.stanford.edu/docs/lsc/wikikg90mv2/.  Accessed: 2024-04-16.
  * Khemani et al. (2024) Bharti Khemani, Shruti Patil, Ketan Kotecha, and Sudeep Tanwar. 2024.  A review of graph neural networks: concepts, architectures, techniques, challenges, datasets, applications, and future directions.  _Journal of Big Data_ 11, 1 (2024), 18.
  * Lao and Cohen (2010) Ni Lao and William W Cohen. 2010.  Relational retrieval using a combination of path-constrained random walks.  _Machine learning_ 81 (2010), 53–67.
  * Reimers and Gurevych (2019) Nils Reimers and Iryna Gurevych. 2019.  Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks.  _arXiv preprint arXiv:1908.10084_ (2019).
  * Song et al. (2020) Kaitao Song, Xu Tan, Tao Qin, Jian Lu, and Tie-Yan Liu. 2020.  MPNet: Masked and Permuted Pre-training for Language Understanding.  _arXiv preprint arXiv:2004.09297_ (2020).
  * Zhao et al. (2023) Wenting Zhao, Justin Chiu, Claire Cardie, and Alexander Rush. 2023.  Hop, Union, Generate: Explainable Multi-hop Reasoning without Rationale Supervision. In _Proceedings of the 2023 Conference on Empirical Methods in Natural Language Processing_ , Houda Bouamor, Juan Pino, and Kalika Bali (Eds.). Association for Computational Linguistics, Singapore, 16119–16130.  https://doi.org/10.18653/v1/2023.emnlp-main.1001


# Exploring the Elastic Properties and Fracture Patterns of Me-Graphene
Monolayers and Nanotubes through Reactive Molecular Dynamics Simulations

Marcelo L. Pereira Junior1,†, José. M. De Sousa2,3,‡, Wjefferson H. S.
Brandão4,§,
Douglas. S. Galvão3,5,$, Alexandre F. Fonseca3,∥, and Luiz A. Ribeiro
Junior6,✠

1University of Brasília, Faculty of Technology, Department of Electrical
Engineering, Brasília, Brazil
2Federal Institute of Education, Science and Technology of Piauí – IFPI, São
Raimundo Notato, 64770-000, Piauí, Brazil.
3Applied Physics Department, “Gleb Wataghin” Institute of Physics - IFGW,
University of Campinas – UNICAMP, Rua Sérgio Buarque de Holanda, 777 - Cidade
Universitária, 13083-859, Campinas, São Paulo, Brazil.
4Department of Physics, Federal University of Piauí, Terezina, Piauí, Brazil
5Center for Computing in Engineering & Sciences, Unicamp, Campinas, SP, Brazil
6University of Brasília, Institute of Physics, Brasília, Brazil
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>.

Abstract: Me-graphene (MeG) is a novel two-dimensional (2D) carbon allotrope.
It is composed of sp2\- and sp3-hybridized carbon atoms assembled
topologically from C$-$(C3H2)4 molecules. Due to its attractive electronic and
structural properties, it is important to study the mechanical behavior of MeG
in its monolayer and nanotube topologies. In this work, we carried out fully
atomistic reactive molecular dynamics simulations using the ReaxFF force field
to investigate their mechanical properties and fracture patterns. Our results
indicate that Young’s modulus of MeG monolayers is about 511 GPa and in the
range of 583-664 GPa for the nanotubes investigated here. Moreover, MeG
monolayers exhibit a direct transition from elastic to complete fracture under
critical strain without a plastic regime, while MeG nanotubes show an extended
flat plastic regime before total fracture.

##  1 Introduction

Nanoscience is a rapidly growing area with great potential for developing new
materials [1, 2, 3]. The experimental realization of graphene [4] has been a
breakthrough. It has since been extensively studied in nanoelectronics due to
its lightweight nature, high mechanical strength, good transparency, and
excellent electrical and thermal conductivity [5, 6, 7]. This achievement has
inspired numerous experimental [8, 9] and theoretical [10, 11] investigations
on developing new 2D carbon allotropes. The primary objective is to create
all-carbon 2D materials sharing some of graphene’s exceptional properties,
which could overcome its limitations in flat electronics, particularly its
lack of an electronic bandgap [12].

Carbon allotropes such as graphite [13], carbon nanotubes [14], and graphene
[4] are mainly composed of sp2-hybridized carbon atoms forming hexagonal
rings, which provide excellent structural stability. In contrast, pentagonal
rings are less commonly found in stable all-carbon materials such as C60
molecules [15], primarily due to unfavorable bonding angles that may lead to
structural instability. However, 2D carbon allotropes with pentagonal rings
can exhibit semiconductor behavior [16], and the search for new stable carbon
structures containing pentagons is now attracting much interest [17].

A novel 2D carbon allotrope, similar to graphene, has recently been proposed
and named Me-graphene (MeG) [18]. It is composed of 5-6-8 rings of carbon
atoms with sp2 and sp3 hybridizations, resulting in a semiconductor structure.
Topologically, MeG can be formed by assembling C$-$(C3H2)4 molecules. It
exhibits a buckled layer but no negative frequencies in the phonon spectra,
suggesting its structural stability. MeG exhibits properties between graphene
and penta-graphene in terms of energy formation and non-null electronic
bandgap. It also presents a near-zero Poisson’s ratio behavior over a wide
range of strain values. It has an indirect band gap of 2.04 eV, which evolves
to a direct band gap of 2.62 eV under compressive strain. Furthermore, it is
predicted a high hole mobility of 1.60$\times 10^{5}$ cm2V-1s-1 at 300 K.
Since MeG presents electronic and structural features desirable in flat
nanoelectronics applications, it is important to investigate its mechanical
behavior, for both monolayer and nanotube forms, to improve our understanding
of these family of nanostructures.

In the present work, we have used full atomistic reactive (ReaxFF) classical
molecular dynamics (MD) simulations to investigate the mechanical properties
and fracture patterns of MeG monolayers (MeGM) and nanotubes (MeGNTs), see
Figure 1. Our simulations considered different tube diameters and temperature
regimes. Young’s modulus values of MeG monolayers were estimated to be about
511 GPa and can vary between 583-664 GPa for the nanotubes considered here.
Upon reaching a critical strain, MeGM undergoes complete fracture without
exhibiting a plastic deformation regime (brittle behavior). MeGNTs, on the
other hand, exhibit a flat plastic region between the elastic and the
completely fractured regimes.

##  2 Methodology

We carried out MD Simulations using the LAMMPS [19] code and using the
interatomic reactive force field (ReaxFF) [20]. ReaxFF, developed by van Duin,
Goddard III, and collaborators, connects quantum chemical force fields and
empirical bonding energy terms. It can be parameterized using available
experimental data and/or from density functional theory (DFT) simulations. It
goes beyond non-reactive force fields, allowing the description of forming and
breaking bonds, which is necessary to investigate non-linear structural
deformations and fracture dynamics. The ReaxFF parameter set utilized in this
study is described in references [21, 20].

Figure 1: Schematic representation of (a-b) Me-graphene monolayers (MeGMs) and
(c-e) Me-graphene nanotubes (MeGNTs) topologies. The MeGNTs studied here have
the same length (L) but distinct diameters (d). See text for discussions.

In Figure 1(a), we present the MeGM structure, where its lattice unit cell
(magnified in Figure 1(b)) is indicated by a white square and the rolling
direction used to generate the MeGNTS by a curvy arrow. The MeG carbon atoms
forming the different rings are indicated in the colored box (red, blue, and
yellow for fused pentagons, octagons, and hexagons). The ‘PaperChain’
representation implemented in the molecular visual dynamics (VMD) [22]
software was used to render the ring composition. This technique finds all
rings up to a user-defined maximum size by analyzing the molecular topology
and fitting a polyhedron to the selected atoms and the ring centroid [22].

In Figure 1(c), we present a side view of the MeGNTs, with their length
indicated by L. Figures 1(d) and 1(e) illustrate the front views of MeGNTs of
different diameters (d). The white rectangle in Figure 1(e) highlights MeGMs
and MeGNTs buckled morphologies. The MD snapshots were rendered using VMD
software, while the Sculptor [23] software suite was used to generate the
MeGNT structures.

The used MeGM structure models investigated here are composed of an 18x18
replication of the MeG unit cell ($5.74\times 5.74$ Å2 with 13 carbon atoms),
resulting in a membrane with dimensions of $103.3\times 103.3$ Å2 and
containing a total of 4212 carbon atoms. The MeGNT(3,0), MeGNT(12,0), and
MeGNT(24,0) structures were generated from $3\times 18$ (d${}_{(3,0)}=5.5$ Å),
$12\times 18$ (d${}_{(12,0)}=21.9$ Å), and $24\times 18$ (d${}_{(24,0)}=43.9$
Å) supercells, respectively. The length of all MeGNTs is the same and equals
103.3 Å. The total number of carbon atoms in MeGNT(3,0), MeGNT(12,0), and
MeGNT(24,0) are 702, 2808, and 5616, respectively.

MeGNTs are labelled using integer indices similar to those for graphene
nanotubes [24], i.e., MeGNT$(m,n)$. $m$ and $n$ are used in a linear
combination to determine the nanotube circumference. The MeG unit cell has
symmetric orthogonal lattice vectors so that $a$ and $b$ are equal. Therefore,
MeGNT$(m,0)$ and MeGNT$(0,m)$ nanotubes are identical. As mentioned above, we
have considered the MeGNT(3,0), MeGNT(12,0), and MeGNT(24,0) structures. The
diameter of each tube was calculated using d${}_{(m,0)}=ma/\pi$, where $a$ is
the lattice constant.

We used the velocity Verlet algorithm with a time step of 0.05 fs [25] to
integrate Newton’s equation of motion. Before carrying out the tensile loading
simulations on our model nanostructures, we thermalized them using a Nosé-
Hoover thermostat for 50 ps [26]. Afterwards, an NPT simulation with zero
pressure was performed for an additional 50 ps to eliminate any residual
lattice stress in the thermalized systems [27]. To obtain the elastic
properties and fracture patterns at 300 K, we stretched MeGM and MeGNTs using
a constant engineering tensile strain rate of $10^{-6}/$fs within an NVT
ensemble. The mechanical properties of MeGM and MeGNTs were studied by
applying stresses to these nanostructures along the x-direction (layers) and
z-direction (tubes), respectively, while the other directions are allowed to
relax freely. The stretching processes are maintained until the complete
structural failure (fracture).

##  3 Results

Initially, we analyzed the MeGM bond and angle distributions after lattice
thermalization at various temperatures, ranging from 10 K up to 400 K. In
Figure 2(a), we present the MeGM bond distribution profile, which exhibits two
prominent peaks centered at 1.41 Å and 1.46 Å at low temperatures. A third
peak appears at around 1.56 Å, representing the C$-$C bonds between fused
pentagonal rings (see Figure 1). These results are consistent with the
original data reported in reference [18] and indicate the MeGM structural
stability with the coexistence of sp2 and sp3-like bonds.

Figure 2: MeGM (a) bond and (b) angle distributions (in arbitrary units), as a
function of the temperature.

As expected, the average C$-$C bond-length values increase as the temperature
increases due to larger amplitude thermal vibrations of the carbon atoms. The
bond distribution profile changes significantly as the temperature increases.
The three peaks observed at 10 K are no longer distinguishable at temperatures
higher than 300 K, with only one peak at room temperature centered at
approximately 1.45 Å. It is worth noting that MeGM retains its structural
integrity even at temperatures as high as 3500 K, according to results
reported in reference [18]. No bond breaks or new structural arrangements were
observed, even under extreme conditions.

MeGM maintains its buckled morphology, as suggested by the results from Figure
2(b), where four well-defined angle distributions are present, with peaks
centered at around 92∘, 104∘, 119∘, and 132∘. The angle formed by the C$-$C
bonds among the fused pentagonal rings is represented by a small peak centered
at approximately 1.61 Å. While heating tends to spread the peaks, the initial
angle distribution profile remains, except for the peak centered at
approximately 161∘, which disappears for temperatures higher than 100 K (see
Figure 2(b)).

Figure 3 presents the stress-strain curves for MeGM and MeGNTS at 300 K, with
the stretching tensile loading applied along the x- and z-directions. As
expected, the system topology plays a significant role in defining the
mechanical properties. The stress-strain curves show that MeG monolayers
undergo a direct transition from elastic to a completely fractured regime
(zero stress regions) at a critical strain. In contrast, MeG nanotubes exhibit
a flat-like plastic stage before fracturing. The distinct curve profiles also
reveal that MeGNTs can stand more strain than MeGM before fracturing. Also,
for the MeNTs, the results suggest that the flat-like plastic regime increases
as the tube diameter increases. This can be explained by the fact that while
the monolayers do not have many ways to dissipate the accumulated stress, the
tubular topology allows many pathways to occur, thus preventing early and
brittle fractures.

Figure 3: Stress-strain curves for MeGM and MeGNTs as a function of the
applied strain.

In Table 1, we present a summary with some of the estimated elastic parameters
extracted from the stress-strain curves of Figure 3. The $Y_{M}$ denotes
Young’s modulus values. The $\sigma_{C}$ is the critical stress value.
$\varepsilon_{C}$ denotes the value at the moment of the first fracture. We
calculated the virial stress $\sigma_{C}$ along the stretching direction
according to the approach described in reference [28]. MeGNTs with smaller
diameters are stiffer and tend to present higher $Y_{M}$ and smaller critical
strains. In the cases studied here, the MeGNT(24,0) has a $Y_{M}$ value close
to the MeGM one. As the diameter of the tubes increases beyond MeGNT(24,0),
the decrease in Young’s modulus implies that the impact of curvature effects
becomes less significant.

As mentioned above, the MeGNTs stand higher strain values before breaking than
the MeGM, but can stand higher stress values as can be inferred from the
Figure 3 and the $\sigma_{C}$ ($\varepsilon_{C}$) values. The value is larger
than the corresponding graphene one (about 130 GPa (25%) [29]), but in the
same range as other 2D all-carbon allotropes theoretically proposed [30, 31,
32, 33]. The critical stress values for the MeGNTs are almost the same for the
different diameters considered here, while the critical strain is slightly
diameter dependent (see Table 1), thus indicating that one of the curvature
effects on the structure is making them more compliant.

System |  $Y_{M}$ (GPa) |  $\sigma_{C}$ (GPa) |  $\epsilon_{C}$ (%)
---|---|---|---
MeGM | 511$\pm$15 | 180 | 13
MeGNT(3,0) | 664$\pm$32 | 103 | 16
MeGNT(12,0) | 629$\pm$24 | 97 | 18
MeGNT(24,0) | 583$\pm$16 | 100 | 22
Table 1: Elastic properties for the MeGM and MeGNTs estimated from the stress-
strain curves shown in Figure 3.

Next, we examined the fracture patterns of MeGM obtained from the MD
simulations. In Figure 4, we present typical MD snapshots of the fracture
process for MeGM under x-directional stress at 300 K. The color scheme used in
this Figure represents the von Mises (VM) stress per-atom values [34], with
red and blue colors indicating high- and low-stress accumulation,
respectively. These VM values can better identify the point or regions where
the fracture starts and/or propagates. Further information about VM
calculations can be found in references [35, 36].

Figure 4: Representative MD snapshots for the fracture process of MeGM at
300K, stressed along the x-direction. Panels (a), (b-c), (d), (e-f), and (g)
show the MeGM bond arrangements when it is subjected 0.0%, 13.25%, 13.30%,
13.35%, and 13.50% of strain, respectively.

Figures 4(a) and 4(b-c) show the MeGM lattice at 0% and 13.25% of strain,
respectively. At a critical strain of 13.25%, the MeGM undergoes an abrupt
transition from an elastic to a completely fractured regime. Up to 13.25% of
strain, the deformation of the structure is entirely elastic, meaning that it
could return to its original configuration when released from that strain
without forming any defects. However, the first C-C bond breaking occurs at a
slightly higher strain, as shown in Figure 4(d). After that, several
additional bonds break, and linear carbon chain (LAC) formations, as shown in
Figure 4(e). The structure is wholly fractured by less than 0.2% of strain
above the critical one, justifying this transition to be considered abrupt.
The fracture occurs due to the breaking of two parallel bonds connecting a
six-membered atomic ring with two neighboring rings composed of eight atoms,
as illustrated in Figures 4(c) and 4(d). The crack then rapidly propagates
along the y-direction, opposite to the direction of tensile loading, as shown
in Figures 4(e-f). At the end of the stretching process, the MeGM fractures
almost into two symmetrical unstressed moieties, as seen in Figure 4(g), which
also shows the LAC formation.

Finally, we present some representative MD snapshots for the fracture process
of the MeGNTs studied here. Again, we use a color scheme representing the VM
stress per-atom values to help visualization of the fracture process. For
MeGNT(3,0), the critical strain is 15.80% (see Figures 5(a-e)). Figures 5(b-c)
show that the first bond breaking occurs for 13.30% of strain, with all rings
undergoing deformations along the longitudinal direction of the tube. The
breaking bonds align with the direction of the applied strain, as shown in
Figure 5(b). The nanotube fragmentation creates LACS, resulting in two
unstressed moieties, as shown in Figure 5(e). Similar fracture patterns were
observed for MeGNT(12,0) and MeGNT(24,0). Figures 5(g) and 5(j) show the
morphologies of these nanotubes in their critical stress regime. The first
bond breakings occur at 18.30% (Figure 5(h)) and 22.30% (Figure 5(k)) of
strain for MeGNT(12,0) and MeGNT(24,0), respectively. Subsequently, these
nanotubes separate into two parts at 31.00% (Figure 5(i)) and 26.00% (Figure
5(k)) of strain.

Figure 5: Representative MD snapshots for the fracture process of (a-e)
MeGNT(3,0), (g-i) MeGNT(12,0), and (j-l) MeGNT(24,0) at 300K, stressed along
their main axial direction.

##  4 Conclusion

In summary, we carried out extensive fully atomistic reactive (ReaxFF)
molecular dynamics (MD) simulations to investigate the mechanical properties
and fracture patterns of Me-graphene monolayer (MeGM) and ant their nanotubes
(MeGNTs) under different temperature regimes and considering tubes of
different diameters. Our results reveal that the MeGM Young’s modulus value is
approximately 511 GPa, while the corresponding ones for the MeGNTs are in the
range of 583-664 GPa. MeGM exhibits an abrupt transition from elastic to
completely fractured regimes (without a plastic regime), while MeGNTs display
flat-like plastic stages before total fracture. Additionally, the crack
propagation for all structures occurs opposite to the applied tensile loading.

As expected, the system topology plays a significant role in defining the
mechanical properties. The distinct profiles for the stress-strain curves
between MeGM and MeGNTs revealed that MeGNTs could stand higher stress values
than MeGM, although with lower stress values before breaking. These values are
higher than the corresponding graphene ones (about 130 GPa) but in the same
range as other 2D all-carbon allotropes theoretically proposed in the
literature.

## Acknowledgement

This work was financed in part by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior (CAPES) - Finance Code 001 and grant
88887.691997/2022-00, Conselho Nacional de Desenvolvimento Científico e
Tecnológico (CNPq), FAP-DF, FAPESP and FAPEPI. We thank the Center for
Computing in Engineering and Sciences at Unicamp for financial support through
the FAPESP/CEPID Grants #2013/08293-7 and #2018/11352-7. L.A.R.J acknowledges
the financial support from FAP-DF grants $00193-00000857/2021-14$,
$00193-00000853/2021-28$, $00193-00000811/2021-97$, and,
$00193.00001808/2022-71$ CNPq grants $302922/2021-0$ and $350176/2022-1$.
L.A.R.J. gratefully acknowledges the support from ABIN grant 08/2019 and
Fundação de Apoio à Pesquisa (FUNAPE), Edital 02/2022 - Formulário de
Inscrição N.4. L.R.S. acknowledges the National Institute of Science and
Technology of Complex Systems (Brazil). L.A.R.J. acknowledges Núcleo de
Computação de Alto Desempenho (NACAD) and for providing the computational
facilities. A.F.F. thanks the Brazilian Agency CNPq for Grant No.
303284/2021-8 and São Paulo Research Foundation (FAPESP) for Grant No.
#2020/02044-9. This work used resources of the Centro Nacional de
Processamento de Alto Desempenho em São Paulo (CENAPAD-SP). The authors
acknowledge the National Laboratory for Scientific Computing (LNCC/MCTI,
Brazil) for providing HPC resources of the SDumont supercomputer, contributing
to the research results reported within this paper. URL:
http://sdumont.lncc.br.

## References

  * [1] George M Whitesides.  Nanoscience, nanotechnology, and chemistry.  Small, 1(2):172–179, 2005.
  * [2] James E Hutchison.  Greener nanoscience: a proactive approach to advancing applications and reducing implications of nanotechnology.  ACS nano, 2(3):395–402, 2008.
  * [3] Tibor Braun, Andras Schubert, and Sandor Zsindely.  Nanoscience and nanotecnology on the balance.  Scientometrics, 38(2):321–325, 1997.
  * [4] Kostya S Novoselov, Andre K Geim, Sergei V Morozov, De-eng Jiang, Yanshui Zhang, Sergey V Dubonos, Irina V Grigorieva, and Alexandr A Firsov.  Electric field effect in atomically thin carbon films.  science, 306(5696):666–669, 2004.
  * [5] Andre Konstantin Geim.  Graphene: status and prospects.  science, 324(5934):1530–1534, 2009.
  * [6] Chun Hung Lui, Li Liu, Kin Fai Mak, George W Flynn, and Tony F Heinz.  Ultraflat graphene.  Nature, 462(7271):339–341, 2009.
  * [7] Xiangfan Xu, Luiz F. C. Pereira, Yu Wang, Jing Wu, Kaiwen Zhang, Xiangming Zhao, Sukang Bae, Cong Tinh Bui, Rongguo Xie, John T. L. Thong, Byung Hee Hong, Kian Ping Loh, Davide Donadio, Baowen Li, and Barbaros Özyilmaz.  Length-dependent thermal conductivity in suspended single-layer graphene.  Nature Communications, 5(1):3689, 2014.
  * [8] Bertil Sundqvist.  Carbon under pressure.  Physics Reports, 909:1–73, 2021.
  * [9] Hongjie Tang, Colin M Hessel, Jiangyan Wang, Nailiang Yang, Ranbo Yu, Huijun Zhao, and Dan Wang.  Two-dimensional carbon leading to new photoconversion processes.  Chemical Society Reviews, 43(13):4281–4299, 2014.
  * [10] Li-Chun Xu, Ru-Zhi Wang, Mao-Sheng Miao, Xiao-Lin Wei, Yuan-Ping Chen, Hui Yan, Woon-Ming Lau, Li-Min Liu, and Yan-Ming Ma.  Two dimensional dirac carbon allotropes from graphene.  Nanoscale, 6(2):1113–1118, 2014.
  * [11] Xinyu Fan, Jiao Li, and Gang Chen.  New carbon allotropes with metallic conducting properties: a first-principles prediction.  RSC advances, 7(28):17417–17426, 2017.
  * [12] Freddie Withers, Marc Dubois, and Alexander K Savchenko.  Electron properties of fluorinated single-layer graphene transistors.  Physical review B, 82(7):073403, 2010.
  * [13] DDL Chung.  Exfoliation of graphite.  Journal of materials science, 22(12):4190–4198, 1987.
  * [14] Sumio Iijima.  Helical microtubules of graphitic carbon.  nature, 354(6348):56–58, 1991.
  * [15] Harold W Kroto, James R Heath, Sean C O’Brien, Robert F Curl, and Richard E Smalley.  C60: Buckminsterfullerene.  nature, 318(6042):162–163, 1985.
  * [16] Shunhong Zhang, Jian Zhou, Qian Wang, Xiaoshuang Chen, Yoshiyuki Kawazoe, and Puru Jena.  Penta-graphene: A new carbon allotrope.  Proceedings of the National Academy of Sciences, 112(8):2372–2377, 2015.
  * [17] Susmita Jana, Arka Bandyopadhyay, Sujoy Datta, Debaprem Bhattacharya, and Debnarayan Jana.  Emerging properties of carbon based 2d material beyond graphene.  Journal of Physics: Condensed Matter, 34(5):053001, 2021.
  * [18] Zhiwen Zhuo, Xiaojun Wu, and Jinlong Yang.  Me-graphene: a graphene allotrope with near zero poisson’s ratio, sizeable band gap, and high carrier mobility.  Nanoscale, 12(37):19359–19366, 2020.
  * [19] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of computational physics, 117(1):1–19, 1995.
  * [20] Adri CT Van Duin, Siddharth Dasgupta, Francois Lorant, and William A Goddard.  Reaxff: a reactive force field for hydrocarbons.  The Journal of Physical Chemistry A, 105(41):9396–9409, 2001.
  * [21] Jonathan E Mueller, Adri CT Van Duin, and William A Goddard III.  Development and validation of reaxff reactive force field for hydrocarbon chemistry catalyzed by nickel.  The Journal of Physical Chemistry C, 114(11):4939–4949, 2010.
  * [22] William Humphrey, Andrew Dalke, and Klaus Schulten.  Vmd: visual molecular dynamics.  Journal of molecular graphics, 14(1):33–38, 1996.
  * [23] William Humphrey, Andrew Dalke, Klaus Schulten, et al.  Vmd: visual molecular dynamics.  Journal of molecular graphics, 14(1):33–38, 1996.
  * [24] Ge G Samsonidze, A R Saito, D A Jorio, E Pimenta, Filho Souza, F A Grüneis, D G Dresselhaus, MS Dresselhausa, et al.  The concept of cutting lines in carbon nanotube science.  Journal of nanoscience and nanotechnology, 3(6):431–458, 2003.
  * [25] Nicos S Martys and Raymond D Mountain.  Velocity verlet algorithm for dissipative-particle-dynamics-based models of suspensions.  Physical Review E, 59(3):3733, 1999.
  * [26] Denis J Evans and Brad Lee Holian.  The nose–hoover thermostat.  The Journal of chemical physics, 83(8):4069–4074, 1985.
  * [27] Hans C Andersen.  Molecular dynamics simulations at constant pressure and/or temperature.  The Journal of chemical physics, 72(4):2384–2393, 1980.
  * [28] J M De Sousa.  Nanostructures failures and fully atomistic molecular dynamics simulations.  2021\.
  * [29] Changgu Lee, Xiaoding Wei, Jeffrey W Kysar, and James Hone.  Measurement of the elastic properties and intrinsic strength of monolayer graphene.  science, 321(5887):385–388, 2008.
  * [30] ML Pereira Júnior, WF da Cunha, WF Giozza, RT de Sousa Junior, and LA Ribeiro Junior.  Irida-graphene: A new 2d carbon allotrope.  FlatChem, page 100469, 2023.
  * [31] ML Pereira Júnior, JM De Sousa, WHS Brandao, AL Aguiar, RA Bizao, LA Ribeiro Júnior, and DS Galvao.  On the elastic properties of single-walled phagraphene nanotubes.  Chemical Physics Letters, 756:137830, 2020.
  * [32] Chao Sui, Yushun Zhao, Zhisen Zhang, Jianying He, Zhiliang Zhang, Xiaodong He, Chao Wang, and Jianyang Wu.  Morphology-controlled tensile mechanical characteristics in graphene allotropes.  Acs Omega, 2(7):3977–3988, 2017.
  * [33] Hao Sun, Sankha Mukherjee, Matthew Daly, Ajay Krishnan, Manohar Harsha Karigerasi, and Chandra Veer Singh.  New insights into the structure-nonlinear mechanical property relations for graphene allotropes.  Carbon, 110:443–457, 2016.
  * [34] R Von Mises.  Mechanics of solid bodies in the plastically-deformable state.  Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1913.
  * [35] Molecular Massively Parallel Simulator.  Lammps.  Theater, 2012, 2012.
  * [36] WHS Brandão, AL Aguiar, and J M De Sousa.  Atomistic computational modeling of temperature effects in fracture toughness and degradation of penta-graphene monolayer.  Chemical Physics Letters, 778:138793, 2021.