# Organic Structures Emerging From Bio-Inspired Graph-Rewriting Automata

Paul Cousin and Aude Maignan
Univ. Grenoble Alpes, CNRS, Grenoble INP, LJK
Grenoble, France

###### Abstract

Graph-Rewriting Automata (GRA) are an extension of Cellular Automata to a
dynamic structure using local graph-rewriting rules. This work introduces
linear algebra based tools that allow for a practical investigation of their
behavior in deeply extended time scales. A natural subset of GRA is explored
in different ways thereby demonstrating the benefits of this method. Some
elements of the subset were discovered to create chaotic patterns of growth
and others to generate organic-looking graph structures. These phenomena
suggest a strong relevance of GRA in the modeling natural complex systems. The
approach presented here can be easily adapted to a wide range of GRA beyond
the chosen subset.

Keywords — cellular automata, graph-rewriting automata, dynamical systems,
artificial life, complexity.

##  1 Introduction

The idea of extending the study of discrete dynamical systems like cellular
automata to systems with evolving typologies has already been approached in
several ways. DEM-Systems [1, 2, 3], Structurally Dynamic Cellular Automata
[4, 5, 6], Generative Network Automata [7, 8, 9] and Graph-Rewriting Automata
[10, 11, 12] are related concepts sharing this goal. These systems can produce
a range of new behaviors compared to their fixed topology counterparts,
however their implementation is much less straightforward.

In this paper, new linear algebra based tools will be presented, then applied
to study a natural subset of GRA. Links to the Mathematica source code, a more
recent GPU-accelerated Python implementation, and more information can be
found at paulcousin.github.io/graph-rewriting-automata.

A GRA consists of an initial graph $G_{0}$ defined at time step $t=0$ and a
rule to iteratively evolve it to any discrete time step $t>0$. $G_{t}$ is the
graph obtained at time $t$. The studied subset of GRA will be defined by the
following restrictions:

  1. 1.

3-regular, undirected and finite graphs,

  2. 2.

binary labeling of vertices $v$ with two possible states $s$:

    1. •

$s(v)=1$, called “alive” and colored purple  ,

    2. •

$s(v)=0$, called “dead” and colored orange  ,

  3. 3.

rules local to a vertex and its adjacent vertices, deterministic and applied
simultaneously on the entire graph,

  4. 4.

only one type of topology altering operation called division.

Divisions will be operated as follows:

In this model, vertices can, for example, be seen as cells that divide under
specific local conditions (internal and environmental). An initially simple
graph can thus grow into a large organism with a complex self-organized
structure.

|  $\rightarrow$ |  $\rightarrow$ |  $\rightarrow$ |  $\rightarrow$ | $\rightarrow\;\dots$ |
---|---|---|---|---|---|---

$G_{0}$ $G_{1}$ $G_{2}$ $G_{3}$ $G_{4}$

Figure 1: Example of evolution starting from a simple initial graph.

##  2 Graphs and Rules

A graph $G$ is fully described by an adjacency matrix $\mathcal{A}$ and a
state vector $\mathcal{S}$:
$G_{t}=\\{\,\mathcal{A}_{t}\,,\,\mathcal{S}_{t}\,\\}$

$G_{0}\equiv$ |  | $\mathcal{A}_{0}=\begin{pmatrix}0&1&1&1\\\ 1&0&1&1\\\ 1&1&0&1\\\ 1&1&1&0\end{pmatrix}$ |  | $\mathcal{S}_{0}=\begin{pmatrix}1\\\ 0\\\ 0\\\ 0\end{pmatrix}$
---|---|---|---|---

Figure 2: Example of the graph $G_{0}$ from Figure 1.

The neighborhood $N(v)$ of a vertex $v$ is the set of vertices which are
adjacent to it. The state $s$ of a vertex being either dead or alive and there
being from 0 to 3 alive vertices in its neighborhood, there are only 8
possible configurations $c$, which can be ordered in this way:

$c(v)=4\times s(v)+\sum_{i\in N(v)}s(i)$ (1)

For instance, if $v$ is a dead vertex surrounded by 3 dead vertices, $c(v)$
will be equal to 0 (see Figure 3 for all the other configurations).

Figure 3: Eight possible configurations.

The space of possible rules applicable in the subset of GRA defined in Section
1 is finite. Every rule must specify, for each configuration, whether the
vertex will be alive or dead at $t+1$ and whether it will have undergone a
division, leading to 4 possible final states. Thus, there are $4^{8}=65,536$
possible rules. Each rule can be described by two functions:

$\displaystyle
R:[[0,7]]\rightarrow\\{\;0\equiv\parbox{5.69054pt}{\includegraphics[width=7.39772pt]{graphics/rules/orange.pdf}}\;\,,\;\,1\equiv\parbox{5.69054pt}{\includegraphics[width=7.39772pt]{graphics/rules/purple.pdf}}\;\,\\}\qquad\qquad\qquad\qquad\quad\;$
(2) $\displaystyle R(c(v_{t}))=s(v_{t+1})$

$\displaystyle
R^{\prime}:[[0,7]]\rightarrow\\{\;0\equiv\parbox{17.07164pt}{\includegraphics[width=17.07164pt]{graphics/rules/non-
divised.pdf}}\;,\;1\equiv\parbox{17.07164pt}{\includegraphics[width=17.07164pt]{graphics/rules/divised.pdf}}\;\\}$
(3) $\displaystyle
R^{\prime}(c(v_{t}))=\begin{cases}\begin{array}[]{ll}1,&\text{if c leads to a
division at t+1}\\\ 0,&\text{otherwise}\end{array}\end{cases}$

Every rule can thus be labeled by a unique rule number $n$.

$n=\sum_{i=0}^{7}\left[2^{i}R(i)+2^{i+8}R^{\prime}(i)\right]$ (4)

This labeling system, inspired by the Wolfram code [13], is such that a rule
number in its binary form displays the behavior of the rule. Starting from the
right, the 8 first digits indicate the future state for each configuration as
they have been ordered previously. The 8 following digits show when a division
occurs.

Figure 4: Rule of the example in Figure 1 ($n=765=10,1111,1101_{2}$).

##  3 Implementation

Implementing GRA as described below must leverage a sparse array format to be
efficient.

The application of a rule to a graph will come in several steps which are
strictly equivalent to the rule being applied all at once. $o$ being the order
of the graph and @ being the operator applying a function to every element of
a vector, the first step consists of computing a configuration vector
$\mathcal{C}$.

$\mathcal{C}=\begin{pmatrix}c(v_{1})\\\ \vdots\\\
c(v_{o})\end{pmatrix}=4\times\mathcal{S}+\mathcal{A}\cdot\mathcal{S}$ (5)

$\mathcal{S}$ can then be updated as follows.

$\mathcal{S}=R\,\text{@}\,\mathcal{C}$ (6)

A division vector $\mathcal{D}$ is computed similarly.

$\mathcal{D}=R^{\prime}\text{@}\,\mathcal{C}$ (7)

Divisions can now be performed one by one as a combination of simple
operations on matrices. The first 1 in $\mathcal{D}$ signals the vertex to
divide. For the state vector, suffice to triple the line corresponding to this
vertex. For the adjacency matrix, both the line and the column must be
tripled. Ones then have to be spread across these lines and columns and the
intersection must be filled with a sub-matrix containing zeros on the diagonal
and ones everywhere else. Finally, the first 1 in the division vector is
turned into a 0 and then tripled. This process must be repeated until
$\mathcal{D}$ is a null vector.

$\begin{pmatrix}1\\\ {\color[rgb]{1,0,0}0}\\\ 0\\\
0\end{pmatrix}\rightarrow\begin{pmatrix}1\\\ {\color[rgb]{1,0,0}0}\\\
{\color[rgb]{1,0,0}0}\\\ {\color[rgb]{1,0,0}0}\\\ 0\\\ 0\end{pmatrix}$

(a) State vector

$\begin{pmatrix}0&{\color[rgb]{1,0,0}1}&1&1\\\
{\color[rgb]{1,0,0}1}&{\color[rgb]{.5,0,.5}0}&{\color[rgb]{1,0,0}1}&{\color[rgb]{1,0,0}1}\\\
1&{\color[rgb]{1,0,0}1}&0&1\\\
1&{\color[rgb]{1,0,0}1}&1&0\end{pmatrix}\rightarrow\begin{pmatrix}0&{\color[rgb]{1,0,0}1}&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,.5,0}0}&1&1\\\
{\color[rgb]{1,0,0}1}&{\color[rgb]{.5,0,.5}0}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,.5,0}0}\\\
{\color[rgb]{1,.5,0}0}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{.5,0,.5}0}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{1,0,0}1}&{\color[rgb]{1,.5,0}0}\\\
{\color[rgb]{1,.5,0}0}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{.5,0,.5}1}&{\color[rgb]{.5,0,.5}0}&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,0,0}1}\\\
1&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,0,0}1}&{\color[rgb]{1,.5,0}0}&0&1\\\
1&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,.5,0}0}&{\color[rgb]{1,0,0}1}&1&0\end{pmatrix}$

(b) Adjacency matrix

$\begin{pmatrix}0\\\ {\color[rgb]{1,0,0}1}\\\ 1\\\
1\end{pmatrix}\rightarrow\begin{pmatrix}0\\\ {\color[rgb]{1,0,0}0}\\\ 1\\\
1\end{pmatrix}\rightarrow\begin{pmatrix}0\\\ {\color[rgb]{1,0,0}0}\\\
{\color[rgb]{1,0,0}0}\\\ {\color[rgb]{1,0,0}0}\\\ 1\\\ 1\end{pmatrix}$

(c) Division vector

Figure 5: Division operated on the second vertex of the graph from Figure 2.

##  4 Exploration

An exhaustive exploration will be performed on the set of rules with exactly
one division. This set being color symmetric, if it is the case as well for
$G_{0}$, considering the half of it were the division is applied to a dead
cell will give the full picture. This makes a total of 1024 rules to be
explored.

$n\in\\{i+2^{j}\;|\;i\in[[0,255]],\;j\in[[8,11]]\\}$ (8)

To investigate the behavior of these rules, a simple graph $G_{0}$ depicted in
Figure 6 and containing all 8 configurations will serve as a starting point.
This ensures that $G_{1}\neq G_{0}$.

Figure 6: $G_{0}$: Initial graph with all 8 configurations.

All the results and plots in this paper were computed in Mathematica 13.0
[14], on a MacBook Pro 2020 with an Intel Core i5 Quad-Core 2.0GHz and 16GB of
memory. Graphs are displayed by the default graph plotting function in
Mathematica and can be better seen by enlarging this pdf.

##  5 Growth

One natural way to classify rules is by observing the pattern of growth they
produce. Table 1 shows a classification obtained using several methods, mainly
the least squares method.

Table 1: Number Of Rules By Growth Pattern.

Halted | Linear | Quadratic | Exponential | Unclassified
---|---|---|---|---
| strict | periodic | chaotic |  |  |
422 | 73 | 19 | 3 | 1 | 374 | 132

###  5.1 Halted growth

A total of 422 rules lead to a halt when reaching a cycle. These cycles have
periods of 1, 2, 3, 6, or 8.

Table 2: Period distribution in rules leading to cycles.

Period | 1 | 2 | 3 | 6 | 8
---|---|---|---|---|---
Number of rules | 310 | 102 | 2 | 6 | 2

###  5.2 Linear growth

73 rules produce a strictly linear growth, with the same number of vertices
gained at each time step. 19 rules create a repeating pattern of growth that
averages to linearity. The largest period spanned 42 time steps and was
produced by rule 2183. This comes from the two main structures of the graph: a
loop growing in 7 steps $\\{2,0,0,2,0,0,0\\}$ and a braid growing in 6 steps
$\\{4,0,0,4,2,0\\}$.

---
Figure 7: Rule 2183 at $t=600$ along with a close-up of the braid.

More surprising though are the 3 rules giving rise to a chaotic linear growth:
2222, 2238 and 2239. These rules produce ladder-like loops in which states
change according to an Elementary Cellular Automaton rule [15]. For the sake
of concision, only rule 2222 will be presented.

Figure 8: Rule 2222 at $t=200$.

This rule, although producing a rather simple graph structure, is interesting
by the rhythm of its growth. At any time step, it will either grow of 0, 2, 4
or 6 new vertices which averages to a linear growth. It can be noted that
there are no increase of 6 vertices at once for $t\in[[1217,\,20608]]$.

Figure 9: Rule 2222: linear growth.

Model fit: $o=72.8095+0.630684\times t$ Adjusted $R^{2}$: $0.999992$

What is remarkable here is the fact that this growth does not follow a
repeating pattern but, instead, seems to have chaotic properties. To
illuminate the peculiar aspect of this phenomenon, Figure 10 shows the
distribution of the lengths of time intervals without growth in a log-linear
plot.

Figure 10: Rule 2222: length distribution of time intervals without growth
computed for $0\leq t\leq 200,000$.

###  5.3 Quadratic growth

|  |
---|---|---
$t=200$ | $t=400$ | $t=600$
Figure 11: Rule 2182.

Rule 2182 was found to be unique in the set regarding its growth. It produces
a particular structure that grows almost quadratically, the exponent being
slightly above 2.

Figure 12: Rule 2182: quasi-quadratic growth.

Model fit: $o=0.0620341\times t^{2.01874}$ Adjusted $R^{2}$: $0.999999$

###  5.4 Exponential growth

Rule 256 produces a perfect exponential growth from $t=1$. Among the 373 other
exponential rules, some display simple fractal structures, others are more
complex.

|
---|---
Figure 13: Examples of fractal structures found in exponential rules.

###  5.5 Unclassified

132 rules have not allowed for a confident classification. Suspecting
interesting behaviors to be responsible for these less trivial growth
patterns, all of them were observed individually and a selection of rules
producing remarkable graph structures is presented in Section 6.

##  6 Graph structures

As Figure 11 already hinted at, some rules create organic-looking graph
structures. Figure 14 shows a selection of graphs with this property, labeled
with the number of the rule that produced them. These structures are
particularly evocative. Some even seem familiar, resembling macroscopic algae
or related photosynthetic eukaryotes.

|
---|---
549 | 618
|
1062 | 1111
|
2199 | 2236
Figure 14: Some organic-looking graph structures.

##  7 Conclusion

In this paper, a new approach to Graph-Rewriting Automata based on linear
algebra was presented. When leveraging a sparse array format, it can be used
to implement GRA much more efficiently than with previous methods. To
demonstrate the benefits of this approach, a natural subset of GRA was
explored only using a MacBook Pro 2020.

A focus was first put on the evolution of the number of vertices over time,
leading to a classification of rules by growth pattern. Several noteworthy
phenomena were discovered through this approach. It also led to a category of
unclassified growths in which interesting behaviors were suspected. The graphs
of this category were thus visually inspected on a individual basis and some
remarkable structures were discovered.

In future work, the ideas presented here could easily be adapted to cover a
variety of other types of GRA. Here are a few examples of possible
trajectories:

  1. •

using continuous valued or non-binary discrete states,

  2. •

working with d-regular or non-regular graphs,

  3. •

reaching an extended neighborhood with powers of $\mathcal{A}$,

  4. •

working with directed graphs using $\mathcal{A}\cdot\mathcal{S}$ and
$\mathcal{A}^{T}\cdot\mathcal{S}$,

  5. •

including other topology altering operations.

Other tools to analyse graphs directly from their adjacency matrix and state
vector could also be developed to gain new insights into the behaviors and
properties of GRA.

## Acknowledgment

This work was partially supported by the Université Grenoble Alpes, through an
Excellence Internship which took place between June 7th and July 15th 2022\.

## References

  * [1] R. Edwards and A. Maignan, “A class of discrete dynamical systems with properties of both cellular automata and l-systems,” _Natural Computing_ , vol. 19, no. 3, pp. 609–641, 2020.
  * [2] ——, “Dem-systems: a new type of adaptive system,” _Exploratory papers of automata_ , 2016.
  * [3] ——, “Complex self-reproducing systems,” in _ISCS 2013: Interdisciplinary Symposium on Complex Systems_. Springer, 2014, pp. 65–76.
  * [4] A. Ilachinski and P. Halpern, “Structurally dynamic cellular automata.” 2009.
  * [5] R. Alonso-Sanz and M. Martín, “A structurally dynamic cellular automaton with memory in the hexagonal tessellation,” in _International Conference on Cellular Automata_. Springer, 2006, pp. 30–40.
  * [6] R. Alonso-Sanz, “A structurally dynamic cellular automaton with memory,” _Chaos, Solitons & Fractals_, vol. 32, no. 4, pp. 1285–1295, 2007.
  * [7] H. Sayama, “Generative network automata: A generalized framework for modeling complex dynamical systems with autonomously varying topologies,” in _2007 IEEE Symposium on Artificial Life_. IEEE, 2007, pp. 214–221.
  * [8] H. Sayama and C. Laramee, “Generative network automata: A generalized framework for modeling adaptive network dynamics using graph rewritings,” in _Adaptive networks_. Springer, 2009, pp. 311–332.
  * [9] J. Schmidt and H. Sayama, “Designing and evaluating algorithms for automated discovery of adaptive network models based on generative network automata,” in _2013 IEEE Symposium on Artificial Life (ALife)_. IEEE, 2013, pp. 27–34.
  * [10] K. Tomita, S. Murata, and H. Kurokawa, “Self-description for construction and computation on graph-rewriting automata,” _Artificial Life_ , vol. 13, no. 4, pp. 383–396, 2007.
  * [11] ——, “Asynchronous graph-rewriting automata and simulation of synchronous execution,” in _European Conference on Artificial Life_. Springer, 2007, pp. 865–875.
  * [12] K. Tomita, H. Kurokawa, and S. Murata, “Graph-rewriting automata as a natural extension of cellular automata,” in _Adaptive Networks_. Springer, 2009, pp. 291–309.
  * [13] S. Wolfram, _A New Kind of Science_. Wolfram Media, 2002.
  * [14] W. R. Inc., “Mathematica, Version 13.0,” champaign, IL, 2021.
  * [15] E. W. Weisstein, “Elementary cellular automaton,” _mathworld.wolfram.com_ , 2002.


# Active learning based generative design for discovery of wide band gap
materials

Rui Xin, Edirisuriya M. D. Siriwardane, Yuqi Song, Yong Zhao, Steph-Yves
Louis, Alireza Nasiri
Department of Computer Science and Engineering
University of South Carolina
Columbia, SC 29201
&Jianjun Hu *
Department of Computer Science and Engineering
University of South Carolina
Columbia, SC 29201
<EMAIL_ADDRESS>

###### Abstract

Active learning has been increasingly applied to screening functional
materials from existing materials databases with desired properties. However,
the number of known materials deposited in the popular materials databases
such as ICSD and Materials Project is extremely limited and consists of just a
tiny portion of the vast chemical design space. Herein we present an active
generative inverse design method that combines active learning with a deep
variational autoencoder neural network and a generative adversarial deep
neural network model to discover new materials with a target property in the
whole chemical design space. The application of this method has allowed us to
discover new thermodynamically stable materials with high band gap (SrYF5) and
semiconductors with specified band gap ranges (SrClF3, CaClF5, YCl3, SrC2F3,
AlSCl, As2O3), all of which are verified by first principle DFT calculations.
Our experiments show that while active learning itself may sample chemically
infeasible candidates, these samples help to train effective screening models
for filtering out materials with desired properties from the hypothetical
materials created by the generative model. The experiments show the
effectiveness of our active generative inverse design approach.

_K_ eywords Active learning $\cdot$ generative design $\cdot$ deep neural
networks $\cdot$ band gap $\cdot$ inverse design $\cdot$ Bayesian optimization

##  1 Introduction

Computational discovery of new materials with desired functions has big
potential in a variety of industries. Currently there are two major categories
of approaches including 1) computational screening of existing materials in
materials databases [1] to find candidates with new properties and 2)
generation of new materials using rational design methods[2], crystal
structure prediction [3] or generative machine learning (ML) models [4, 5].
The computational screening approach, however, usually suffers from the
limited known materials in terms of both quantity and diversity. This is
especially true as most known materials in existing databases such as ICSD,
Materials Project[6], and OQMD[7] have only structural information annotated
with a few density functional theories (DFT) [8] calculated primitive
properties such as formation energies or band gaps. However, few of them are
annotated with more interesting macro properties such as thermal conductivity
(<2000) and ion conductivity (<200), which makes it challenging to train a
highly generalizable machine learning model for screening. The generation
approaches also face a variety of challenges, such as the difficulty to
generate diverse novel materials compositions and structures and the same
issue of lacking high-quality property data for generated hypothetical
materials. Essentially, both ML-assisted computational screening and
generation suffer from very limited training data with scarce annotated
property characterization data.

To improve the efficiency of materials screening, inverse design methods [9,
10, 11, 12, 13, 14, 15, 16] have been proposed to search materials for a given
target functionality, in which a sampling or optimization process is used to
select a subset of the whole design space for locating the optimal design
candidate. The sampling process in inverse design methods is usually guided by
a global optimization algorithm[13, 17] together with function performance
evaluations via atomic simulation. Global optimization methods used in inverse
design can be categorized into several types [18] including gradient-based
optimization methods, gradient-free optimization algorithms such as genetic
algorithms, particle swarm optimization, and Covariance Matrix Adaptation
Evolution Strategy (CMA-ES), neural network-based generative approaches, and
Bayesian optimization methods using surrogate evaluation models. Depending on
the discrete or continuous nature of the design targets, not all these
optimization methods can be applied directly. For example, the gradient
methods cannot be easily applied to search the inorganic materials composition
or crystal structure space as they tend to generate infeasible solutions
during optimization processes.

Due to the strong physicochemical and geometric constraints among the atoms in
periodic crystal materials, we find that gradient-based search operators or
crossover/random mutation-based genetic search operators might still be too
random in identifying the next feasible sampling point for expensive atomic
simulation to obtain its property value. Instead, this exploration step can be
guided by a machine learning model as the previously explored samples can be
used to learn the structure to properties relations. The resulting trained
surrogate model can then be used to predict the properties of unexplored
materials with a speed of 100 orders of magnitude faster than conventional DFT
calculation methods. This has led to the emergence of active learning methods
for materials design[8, 19, 20].

The basic idea of active learning is its capability to explore the design
space actively to search candidate with desired properties or identify
informative samples for building property prediction models for labeling so
that better prediction models can be trained with the existing and minimal
number of newly acquired/labeled samples[21]. An active learning model is
typically composed of the following components: an active learner with a
specific machine learning prediction model and a query strategy to recommend
new samples for labeling, and an oracle model to label the suggested samples
[22]. The oracle model can be experimental or computational validations such
as first principle DFT calculation or even machine learning models. In
general, active learning is very useful for tasks where data is not sufficient
or expensive to acquire, which is the case for materials design problems.
Since the active learner actively chooses the training data points, the number
of data points needed for training a model based on active learning is
typically lower or much lower than the number used by the passive learner.
Active learning algorithms can be categorized by the choices of their
components, especially the query strategies, which can be classified into
uncertainty-based, ensemble-based, and diversity-based approaches. The key
idea of an uncertainty-based method is to query those areas about which
predicted value the active sampler is least certain. Uncertainty sampling
strategies can be employed not only with probabilistic models, which help
sample the instance whose best labeling is the least confident, but also non-
probabilistic models, which can be modified to have probabilistic output. Such
as nearest-neighbor-based classifiers in which neighbors are voting on class
labels, and the ratio of the vote on a targeted class represents the posterior
probability[23]. Ensemble-based methods maintain a series of machine learning
models and help sampler decide valuable query points with which these models
most disagree. ACTIVE-DECORATE [24], which uses DECORATE committees to select
good training examples, outperforms traditional ensemble methods such as
Bagging and Boosting. DECORATE Algorithm is designed to use additional
artificially generated training data from data distribution and label them
different from current ensemble predictions in order to generate highly
diverse ensembles. Among all these methods, uncertainty-based active learning
has been most successfully applied to materials design.

Despite the success stories of a few active learning based materials design
cases, there exist several limitations to improve. For example, in [8], active
learning has been applied to find high band gap materials. Their model is
initially trained with only around 2% of the entire search space, and 20 new
samples suggested from the optimization schemes, are added at each
optimization cycle. Their experiments show that the active learning approach
is able to identify the materials with the highest property value after
searching only around 7.0% and 7.7% of the entire database. However as a
demonstration paper, this work does not identify any new materials and have
very limited design space compared to the real-world active inverse design
problem. In addition, current active learning for materials design studies
have been dominated by screening existing materials databases such as ICSD or
Materials Project [8], which limits their search space to a small portion of
the vast chemical design space which may unexpectedly exclude more promising
candidates[25].

We propose a method to combine active learning with generative machine
learning to build the active generative inverse design framework for new
materials discovery and apply it to the discovery of materials with a desired
bandgap. Compared to the existing active learning search of high band gap
materials, which is focused on exploring the known database, our study focuses
on an active learning-based exploration of the unknown design space for high
band gap materials. The framework includes four components, including a
materials composition generator, a variational autoencoder sampler, a band gap
predictor, a Bayesian optimization-based sampler, and a substitution-based
structure predictor. Our composition-based generative machine learning model
[5] can generate millions of new compositions that satisfy basic chemical
rules such as charge neutrality and balanced electronegativity and tend to
have low free energy. We then build a composition-based band gap prediction
model using a Random Forest regressor. A Gaussian process regressor-based
Bayesian optimization algorithm is used to train an ML-based sampling model,
which suggests the next sampling points in the active learning/search process
to maximize expected property improvement. To verify the predicted materials,
we use substitution-based structure predictors to get candidate structures
which can then be used to verify their structural stability and their band
gaps. Our active generative design (AGD) method has allowed us to find a high
band gap material SrYF5 with a band gap of 6.42 eV and seven semiconductor
materials with desired band gap values.

Our contributions include the following:

  * •

We propose a novel active generative inverse design approach for discovering
new materials in the whole design space rather than in a given database. We
find that samples suggested by active learning model may not be chemically
valid. But they are useful for building good property prediction models.

  * •

We show that active learning can be used to increase property prediction
model’s capability to screen out high band gap materials by sampling
informative data points for property evaluation, which is highly desirable for
screening materials with small datasets of property annotations.

  * •

Our evaluation experiments show that active learning based sampling can be
combined with generative models to speed up new materials discovery by
training a ML models from samples traversed by the active learning search
process for screening hypothetical material candidates created by the
generative models.

  * •

Our new pipeline has been able to discover eleven semiconductor materials
(Table2) and two high band-gap materials (SrYF5, RbSr2ClF4) with a DFT-
validated band gap of 6.42 eV and 5.64 eV.

##  2 Materials and Methods

###  2.1 Active Generative Design (AGD) Framework for Materials Discovery

The main idea of our AGD framework is to use the active learning approach to
learn a highly accurate property prediction model and combine it with a
generative model (e.g., an autoencoder model or a Generative Adversarial
Network (GAN model)). Our whole active generative design framework is shown in
Figure 1. It is composed of three modules: an active learning model for
training a high-performance property predictor using limited sample labeling,
and a candidate generation model, and a screening model for the discovery of
new materials. For the purpose of fast benchmark studies here, a machine
learning-based oracle model is trained to simulate the DFT calculations of
band gaps in this study.

The active learning module is the part where data augmentation via active
exploration of the design space is performed to acquire labels of the most
informative sampling points to improve the performance of the machine
learning-based Screening Model. The generative module uses VAE or GAN models
to generate candidate materials and feed them to the Screening model trained
by the enhanced dataset to filter out a number of promising materials
candidate formulas. Their corresponding candidate structures will subsequently
be generated by exploiting their similarity with existing structures using the
element-substitution approach. Formation Energy will also be predicted using
Global Attention Graph Neural Network (GATGNN) [26].

Within our active generative design framework, there are two ways to discover
hypothetical materials:

  * •

1) ALS: direct active learning search in the whole chemical space, which is
done by combining the active learning sampling model with the autoencoder
decoding model.

  * •

2) AL+RF/Roost with MATGAN: active learning-based screening of candidates
generated by the MATGAN generative model. Active learning is used to create
labeled samples for training ML screening models using RF or Roost, which are
then used to filter the samples generated by MATGAN.

Figure 1: Framework for active generative design for materials discovery.

Oracle Module

The Oracle Module is a machine learning model trained with material features
and band gap labels of the whole dataset in our experiments. This module aims
to simulate the DFT function and facilitate the prediction of band gap
property based on a vector of material features. In practice, this Oracle
model should be replaced with DFT calculations or experiments. The descriptors
for our Oracle model are calculated using Roost (Representation Learning from
Stoichiometry) [27]. Roost is a machine learning algorithm that learns the
stoichiometry-to-descriptor map directly from data, the key idea of which is
each material’s composition is represented as a dense weighted graph. Roost
uses Matscholar embedding [28] for embedding and CGCNN [29] for internal
representation. The output network used for the reference model is a Resnet
with five hidden layers and ReLU activation functions.

Active Learning Module

The function of the Active learning Module is to actively explore the whole
design space to identify the informative sampling points and obtain their
property annotations, which are fed to train a highly accurate property
prediction model. In previous active learning-based design studies such as
[8], the samples are directly represented as 145 common chemical features,
which are used to train the Gaussian Surrogate model. New sampling points are
not generated but selected from existing entries in a given database. However,
such physical feature representation of materials is not applicable for active
learning in the uncharted design space because the sampled feature
representations can not be easily mapped back to valid materials
formulas/structures. Actually, this reverse mapping procedure usually makes it
problematic for the Gaussian Surrogate model to generate meaningful samples
that can be decoded into chemical formulas. To address this representation
issue, we introduce the AutoEncoder deep neural network model as a latent
vector generator. It can be used to convert a chemical formula into a
128-dimension latent vector and vice versa. Thus, formulas are not only well
represented but also provide latent vector boundary for each dimension which
makes it easier for the optimization process. To get the property annotation
by the Oracle model, the recommended sampling points from the Surrogate model
will be first decoded into a formula, and then its Magpie descriptor features
are calculated, and the band gap property will be predicted by the Oracle
model.

Latent vectors as generated above, along with their property annotations, will
be fed as input to the Gaussian Surrogate model, which will give uncertainty-
based sampling point suggestions for Oracle Module to evaluate. This suggested
vector has to be converted back to a chemical formula in order to get material
features. However, cases might occur when some suggested latent vectors can
not be decoded correctly into formulas and will negatively affect the design
space exploration performance of the Active-Learning module. Afterward,
labeled new samples will be sent back to the Gaussian Surrogate model for a
new round of recommendations. It should be noted that different initial
datasets and the Gaussian Surrogate model constructed both have an influence
on the sampling points, which affect the final performance of the ML screening
model trained with these data points.

###  2.2 Autoencoder neural network for latent space sampling

Standard autoencoder is a popular representation learning and dimension
reduction method, which is composed of two parts: an encoder and a decoder.
Our autoencoder model is used to map materials compositions into a continuous
space, which can be used by the Bayesian optimizer to sample the design space.
The sampled points can then be decoded back into materials compositions. In
our model, the encoder has two convolutional layers followed by a fully
connected layer and a decoder consists of a fully connected layer and three
convolution layers. The autoencoder (AE) model is shown in Fig 2. The ReLU is
used as activation function in each neural layer except the last layer of the
decoder for which the Sigmoid function is used instead. We use the Adam
optimizer to train the model parameters and the learning rate is set as
$0.001$. We use the cross-entropy loss for the training. The batch size is set
as 128. The autoencoder is trained with 66324 inorganic materials selected
from the Materials Project database.

Figure 2: Autoencoder model.

The material compositions are represented using the one-hot encoding. Through
the simple statistical calculation of the materials in the MP dataset, 87
elements are found, and each element usually has less than 20 atoms in any
specific compound/formula. We then represent each material as a sparse matrix
with 0/1 cell values. Each column represents one of the 87 elements, while the
column vector is a one-hot encoding of the number of atoms of that specific
element. The whole encoding scheme is shown in Figure 3.

Figure 3: One hot encoding of material compositions

###  2.3 Active learning module for efficient sampling

The active learning module in our pipeline is used for efficient exploration
of the design space to identify informative samples and obtain their property
annotations, which allows generating more training samples to train an
accurate property prediction model. The sampled materials can also be used for
screening candidates from them with desirable properties. We use Bayesian
Optimization (BO) [30] as the active Learning sampling strategy to explore the
huge design space.

BO is an Efficient Global Optimization (EGO) algorithm aiming to find the
maximum value of a costly unknown function in as few iterations as possible.
BO builds Gaussian Process Regression, a surrogate model, to not only evaluate
the objective function but also quantify the uncertainty with each prediction.
Acquisition functions such as Expected Improvement (EI) [31] and Upper
Confidence Bound (UCB)[32] are conducted to find the next point to evaluate.

$\mathrm{EI}(\mathbf{x})=\left\\{\begin{array}[]{ll}\left(\mu(\mathbf{x})-f\left(\mathbf{x}^{+}\right)-\xi\right)\Phi(Z)+\sigma(\mathbf{x})\phi(Z)&\text{
if }\sigma(\mathbf{x})>0\\\ 0&\text{ if
}\sigma(\mathbf{x})=0\end{array}\right.$ (1)

Where $\mu(\mathbf{x}$)  and  $\sigma$($\mathbf{x}$)  are the mean and
standard deviation of GP posterior at point x.  $\Phi$ and  $\phi$ are the
Cumulative

Distribution Function (CDF) and Probability Density Function (PDF) of the
standard normal distribution, respectively whereas $Z$ is denoted in Equation
2. $f\left(\mathbf{x}^{+}\right)$ represents the best observed value whereas x
is the corresponding location.

The first term in the upper Equation 1 denotes exploitation, and the second
summation term is the exploration term. EI can adjust the level of balance
between exploitation and exploration through changing parameter $\xi$. The
higher $\xi$ is, the level of importance of posterior mean value improvement
will decrease, while the algorithm will be more inclined to explore unknown
areas.

$Z=\left\\{\begin{array}[]{ll}\frac{\mu(\mathbf{x})-f\left(\mathbf{x}^{+}\right)-\xi}{\sigma(\mathbf{x})}&\text{
if }\sigma(\mathbf{x})>0\\\ 0&\text{ if
}\sigma(\mathbf{x})=0\end{array}\right.$ (2)

Our framework adopted Acquisition function UCB in Equation 3 to probe
promising points. Intuitively, UCB explores the design space to locate the max
value by considering both expected performance and uncertainty across all
solutions. $\kappa$ factor here enables fine adjustments between exploration
and exploitation.

$\mathbf{x}_{t+1}=\underset{\mathbf{x}}{\arg\max}\left(\mu_{t}(\mathbf{x})+\kappa\sigma_{t}(\mathbf{x})\right)$
(3)

Fernando Nogueira’s package [33] is used as the sampler to provide suggestions
about which latent vectors will be generated and probed during the search
process. Depending on the experiments, a different number of initial sampling
points (prior) are used to build the Gaussian Process surrogate model. In
every iteration, the sampler will suggest 20 sampling points. Only those
candidates the can be decoded into a valid material composition will be
registered and used to update the Gaussian Process model. The sampler will
stop sampling when the setting budget is reached.

###  2.4 Template based structure prediction and DFT based validation

To validate the band gap property of promising candidates generated by our
framework, we perform a template-based structure prediction over the generated
compositions to create their structures for DFT calculation after charge
neutrality and electronegativity integrity checking.

We first measure the similarity between a candidate chemical formula and all
the formulas from the Materials Project(MP) database using Earth Mover’s
Distance (EMD) [34]. EMD computes the ratio of each of the elements and the
absolute distance between the elements on the modified Pettifor scale between
two compositions, which leads to greater alignment with chemical understanding
than the Euclidean distance. The MP formulas are then ranked by their EMD
distance, and the CIF file of the most similar formula is identified as the
structure template for element substitution to generate the structures for the
candidate compositions. The CIF files of these new structures will then be
used as input for a pre-trained Global Attention Graph Neural Network (GTGNN)
[26] to predict their formation energy. Only candidates with negative
formation energy will be kept. Finally, DFT calculations will be performed
over these filtered structures to calculate their band gaps and formation
energy. Phonon calculation is also applied to one of the discovered materials
to check its thermodynamic stability.

##  3 Experiments

###  3.1 Band gap Dataset

We download the Materials Project dataset from http://materialsproject.org and
remove duplicate composition entries to get 42,667 samples with band gap data.
We call this the MP dataset. The distribution of band gaps of the MP dataset
is shown in Figure 4.

Figure 4: Distribution of band gaps of the training set.

###  3.2 Experimental Design

To validate the performance of our active generative design framework,
1,180,000 potential material compositions are generated via the MATGAN network
[5]. Bayesian Optimization with the UCB utility function is used as active
learning to explore promising sampling points for machine learning model
training. $Kappa$ is set to 100 since the goal is to generate diverse samples.
The Roost band gap prediction model is trained based on the original and extra
generated sampling points so as to screen high band gap materials out of these
potential materials. In order to better evaluate the effectiveness of our AGD
framework, we designed four experiments according to the availability of the
different amounts of training data. Each model is evaluated by repeating the
experiment 10 times to report statistically significant results. In each
experiment, we split the whole dataset into 60% for training, 20% for
validation, and 20% for testing, respectively.

1) Exp1-oracle model: 42,667 materials from the Materials Project database are
used for training the Oracle model for band gap prediction, which is then used
to screen high potential band gap materials from MATGAN-generated candidates.
This is the ideal case where there are abundant annotated materials samples
for training an ML property prediction model for screening.

2) Exp2: We randomly select 1,000 labeled samples from the MP dataset to train
a baseline ML model (Exp2 BS). We then use the Gaussian Surrogate model for
active learning. After 1,000 sampling points are suggested by the Bayesian
optimizer, 741 samples can be decoded into valid compositions. We use these
newly-explored samples together with the 1000 initial samples to update our
training dataset and train an updated Screening Model (Exp2 AL). Both band gap
prediction models will be evaluated and used for screening MATGAN-generated
candidates.

3) Exp3: Only 300 samples are randomly selected from the MP dataset for the
baseline ML model (Exp3 BS) training and Gaussian Process initialization. We
then generate 1,000 new samples using the active learning model, which leads
to 853 valid compositions. In total, 1153 samples are used to train the
updated ML model (Exp3 AL). We compare the performance between the baseline
screening model (BS) and the auxiliary Active Learning trained screening model
(AL) as we do in Exp 2.

4) Exp4: Only 300 samples whose band gap values are below 3.0 eV are selected
for active learning. This control experiment aims at validating the
performance of our AGD framework when only small, poor quality dataset is
available.

To compare the performance of the model performances trained with a different
number of samples, as shown above, we pick the 8535 test samples from the Exp1
and remove the entries that appear in the Exp2 and Exp3 experiment datasets.
We then select 500 samples from the remaining samples as the test set for
evaluating all the models.

Additionally, a control experiment is conducted in which we adjust the
parameter kappa of the utility function UCB to inspect the active learning
behavior between exploration and exploitation.

###  3.3 Results

####  3.3.1 Prediction Model Capability with Active Learning

Table 1 shows model performance and the number of candidate materials with a
band gap greater than 6.0 eV found with or without active learning. Exp 1 has
the most available training data, which results in an average $R^{2}$
performance of 0.735 when evaluated on its 20% hold-out test set. The baseline
models trained with 1,000 and 300 samples in Exp 2 and Exp 3 have
corresponding $R^{2}$ scores as low as 0.361 and -0.243, showing the under-
performance due to lack of sufficient training samples. The negative $R^{2}$
score of the baseline model performance of Exp3 BS shows that 300 initial
samples are not sufficient to train a reliable model. In contrast, 1,000
initial samples can train a much better model. However, after we apply the
active learning process, the great performance gain is achieved as the model
metric $R^{2}$ scores are improved from 0.361 to 0.611 in Exp2 and from -0.243
to 0.757 in Exp3 when evaluated using the random hold-out test over the
dataset composed of initial samples and newly sampled points. However, the
improvements are a little bit over-estimated here as we can observe that the
performance improvement when evaluated with the external 500-sample test set,
the improvement due to active learning is reduced to 0.19 for Exp3 AL and is
only 0.005 for Exp2 AL. After close examination, we find that this is due to
the mismatch of the training set distribution and the external test set
distribution: most of the samples suggested by active learning during the
search process have high band-gap values greater than 2 eV while the majority
of the whole dataset from which the 500-sample test set are drawn from have
band gap values between 0 and 2 eV (See Figure 4). This shows the samples
generated during the active-learning search are highly biased toward
candidates with higher band gaps.

During active learning of Exp2 and Exp3, we find that 853 decodable new
materials out of 1,000 Bayesian Optimization iterations are found and
evaluated in Exp3 AL, while 741 new materials are found in Exp2 AL. The reason
that Exp2 AL finds fewer new materials given the same budget is that more
initial samples will lead to a better Gaussian Process model with less
uncertainty so that the optimizer tends to suggest smaller steps on the latent
vector. Such latent vectors with small changes tend to result in decoding into
the same materials.

Table 1: Average Model performance and screening results Exp No. | #samples |  $R^{2}$(hold out) |  $R^{2}$ (test set) | MAE | RMSE | #candidates | Recovery Rate
---|---|---|---|---|---|---|---
Exp1 | 42,677 | 0.735$\pm$0.021 | 0.850$\pm$0.023 | 0.505$\pm$0.016 | 0.826$\pm$0.034 | 102 | 100%
Exp2 BS | 1,000 | 0.361$\pm$0.022 | 0.409$\pm$0.016 | 0.882$\pm$0.024 | 1.208$\pm$0.021 | 123 | 51.96%
Exp2 AL | 1,741 | 0.611$\pm$0.009 | 0.414$\pm$0.013 | 0.539$\pm$0.008 | 0.869$\pm$0.010 | 171 | 78.43%
Exp3 BS | 300 | -0.243$\pm$0.082 | 0.227$\pm$0.042 | 1.202$\pm$0.057 | 1.633$\pm$0.054 | 33 | 2.94%
Exp3 AL | 1,153 | 0.757$\pm$0.015 | 0.417$\pm$0.005 | 0.309$\pm$0.019 | 0.592$\pm$0.019 | 81 | 57.84%

In Exp4 AL, we perform active learning based on randomly sampled 300 materials
whose band gap is below 3.0 eV. Figure 5(a) shows our framework can easily
find high band gap potential materials (>6.0 eV) within 30 iterations.
However, these found materials may not be valid if chemical rules are applied.
As a result, most of these generated labeled samples will function as data
augmentation measurements for screening model training.

####  3.3.2 Exploration versus exploitation in active learning: control
experiments

We conduct several control experiments to find the relationship between
exploration and exploitation during the active learning process. To find the
suitable trade-off between them, we conduct experiments using the Upper
Boundary Confidence (UCB) acquisition function, for which the Kappa parameter
can be used to control the exploration rate. Larger Kappa leads to more
exploration, while smaller one leads to more exploitation.

The relationship between exploration and exploitation is shown in Figure 5(a).
The blue line represents pure exploitation, which flattens after about 150
iterations since nearly all candidates with high expected values have been
utilized without exploration. As Kappa gradually increases, which means UCB
explores more, the algorithm will more likely jump out of local minima. As a
result, the number of identified high band gap structures almost always
increases when we increase Kappa from 0 to 10. However, the number of found
structures decreases when Kappa increases from 10 to 100 as active learning
continues to lean towards exploration. The grey line (Kappa=100) and pink line
(Kappa=50) become flatten almost from the beginning since uncertainty has a
rather large enhancement as shown in formula 3. As a result, active learning
arrives at its trade-off and maximizes the discovery of promising structures
when we set Kappa to 10.

####  3.3.3 Finding new materials with targeted band gaps using the active
generative design framework

To understand the effect of using active learning, we apply the trained ML
models to screen the MATGAN-generated hypothetical materials and inspect their
true band gaps calculated using the Oracle model .

The last two columns of Table 1 shows the number of candidate compositions
with predicted band gaps greater than 6 eV and the accuracy measured by the
overlap percentage of candidates with respect to the candidates screened out
by the Oracle model (Exp 1). The recovery rate column of Table 1 also shows
models trained with active learning (AL) augmented data achieve a much higher
recovery rate than the baseline model (BS). Specifically, with 1,000 search
budget, Exp2 AL and Exp3 AL improved the recovery rates of the baseline models
by 26.47% and 54.90% respectively, indicating that active learning can improve
the screening capability of high band gap materials with the ML models trained
with suggested samples and initial samples although they $R^{2}$ score as
evaluated by the external dataset is improved just a little.

(a) Band gap trace of the sampled materials during the search process

(b) Exploration vs Exploitation (UCB utility function).

Figure 5: Active Inverse Design Process

From the screening sets, we have found a high band gap material SrYF5 in Exp1,
Exp2 AL, and Exp3 AL whose predicted band gap is 6.71, 7.02, 7.11eV,
respectively. We use the element substitution to get its structure
(Figure6(a))and calculate its band gap using DFT calculation, which shows a
band gap value of 6.42eV. However, SrYF5 does not appear in the screening
results of models from Exp2 BS and Exp3 BS, which are the baseline models
constructed without conducting active learning. This demonstrates the benefit
of active learning in new materials discovery.

Additionally, our active learning process has found a new material RbSr2ClF4
(structure shown in Figure6(b)) directly during the active learning process.
Its DFT validated band gap value is 5.64eV. This material is a minor variation
of the material RbSr2Cl5 (band gap: 5.20eV), a member of the initial dataset
for active learning training. Altogether, the results show that our active
generative design framework is capable of working well when only small
training dataset is available by efficiently exploring the design space and
improving ML model performance.

(a) Structure of SrYF5

(b) Structure of RbSr2ClF4

Figure 6: Structures of the discovered new materials.

####  3.3.4 Finding new materials with MATGAN+screening approach

The main idea of active generative design is to use active learning to train a
good screening property prediction model, starting with limited training
samples. When there is a sufficient number of training samples, one can simply
build an ML model as the filter to screen the hypothetical samples generated
by GAN models. When there is limited labeled training samples, it is better to
use active learning to improve the ML model using as few expensive (e.g., DFT)
sample annotation as possible. To prove the utility of active learning for
such screening-based new materials discovery, we applied it to the high band
gap materials screening problem starting with only 300/1000 labeled samples.

We first use the MATGAN, a generative adversarial network (GAN) for efficient
generation of new hypothetical inorganic materials [5]. Trained with materials
from the ICSD database, the MATGAN model generates 2.99 million inorganic
formulas, most of which are charge-neutral and electronegativity-balanced. To
simplify the screening process, we trim the generated formula set by filtering
out those formulas with the total atom numbers more than 20 or the number of
elements greater than five, or the number of atoms for each element more than
7. After checking chemical validity, the total number of candidates for
screening is 1.18 million, and we call these samples the MATGAN dataset. We
then use the whole MP dataset to train a Random Forest model to screen this
dataset for potential semiconductors with a band gap between 2.0 eV and 4.0
eV. Similarly, we use the initial samples in Exp1, and Exp2 plus the active
learning explored samples to train two Roost ML models and use them to screen
the MATGAN dataset for band gaps greater than 6.0 eV. From this final
candidate list, we select 18 for DFT formation energy and band gap validation.

##### DFT based Calculations

The density functional theory (DFT) based first principle calculations were
carried out based on the Vienna ab initio simulation package (VASP) [35, 36,
37, 38]. The projected augmented wave (PAW) pseudopotentials were considered
for the electron-ion interactions with 400 eV plane-wave cutoff energy [39,
40]. The exchange-correlation functional was treated by using the generalized
gradient approximation (GGA) based on the Perdew-Burke-Ernzerhof (PBE) method
[41, 42]. The energy convergence criterion was set as 10-5 eV. The atomic
positions were optimized until the forces become less than 10-2 eV/Å. The
Brillouin zone integration for the unit cells was performed using the
$\Gamma$-centered Monkhorst-Pack $k$-meshes. The Formation energies (in
eV/atom) of several materials were calculated based on Eq. 4, where
$E[\mathrm{Material}]$ is the total energy per unit formula of corresponding
material, $E[\textrm{A}_{i}]$ indicates the energy of $i^{\mathrm{th}}$
element of the material, $x_{i}$ represents the number of Ai atoms in a unit
formula, and N is the total number of atoms in a unit
formula($N=\sum_{i}x_{i}$). Those values are reported in Table 2. CaClF5
(-2.09 eV/atom), SrYF5 (-4.10 eV/atom), and SrClF3 (-2.55 eV/atom) show very
low formation energies, implying our model is capable of finding new
semiconductor materials, which are thermodynamically highly stable relative to
the parent compounds of their elements.

$E_{\mathrm{form}}=\frac{1}{N}(E[\mathrm{Material}]-\sum_{i}x_{i}E[\textrm{A}_{i}])$
(4)

Table 2: DFT verified formation energies and the band gaps of several semiconductors Material |  $E_{\mathrm{form}}$(eV/atom) |  $E_{\mathrm{gap}}$ (eV) | Method
---|---|---|---
SrYF5 | -4.10 | 6.42 | AL+Roost with MATGAN
RbSr2ClF4 | -3.35 | 5.64 | Direct Active Learning Search
SiO2F | -1.69 | 3.99 | Direct Active Learning Search
PO2F | -2.10 | 3.84 | Direct Active Learning Search
SH2O3 | -0.10 | 3.74 | Direct Active Learning Search
SiF3 | -4.27 | 3.79 | Direct Active Learning Search
NaF5C | -1.78 | 4.80 | AL+RF with MATGAN
CaClF5 | -2.09 | 3.46 | AL+RF with MATGAN
YCl3 | -1.94 | 2.18 | AL+RF with MATGAN
SrClF3 | -2.55 | 3.18 | AL+RF with MATGAN
SrCl2F3 | -1.91 | 2.27 | AL+RF with MATGAN
AlSeCl | -0.58 | 3.70 | AL+RF with MATGAN
As2SO3 | -1.16 | 2.72 | AL+RF with MATGAN

Moreover, we find that SrClF3 is a stable semiconductor with P1 space-group
symmetry (Figure7). The unit cell lengths of the material are $a=4.44$ Å,
$b=4.75$ Å and $c=8.30$ Å, while the unit cell angles are $\alpha=102.3$,
$\beta=90.1$ and $\gamma=90.1$. The positive frequencies in the phonon
calculations (see Figure 7 (b)) show that the material is dynamically stable
at 0K temperature. The energy above hull, which is 0.046 eV/atom, is obtained
against the competing phases SrF2 and ClF using Pymatgen code [43].

We also carry out density functional perturbation theory (DFPT) calculations
as implemented in VASP code for studying the mechanical properties of SrClF3.
Due to the P1 space-group symmetry of this material, there are 21 independent
elastic constants. The calculated $C_{11}=64.19$ GPa and $C_{12}=62.87$ GPa
elastic constants are almost equal because of $a\approx b$, and $\alpha=\beta$
lattice parameter relationships. Bulk modulus (41.26 GPa), Shear modulus
(15.61 GPa), and Young’s modulus (41.587) were calculated based on Hill
approach [44]. When an infinitesimal strain ($\epsilon$) is applied on a
material, the total energy of that material can be written as
$E=E_{0}+\frac{1}{2}V_{0}\sum_{i,j=1}^{6}C_{i,j}\epsilon_{i}\epsilon_{j}+O(\epsilon^{3})$.
Here, C represents the matrix of second-order elastic constants. The Born
elastic stability conditions show that the matrix should be definite positive,
all eigenvalues of $C$ should be positive, and all the principal components
should be positive. Even though most crystal symmetries have simple
relationships between the elastic constants to show whether the materials are
mechanically stable, triclinic systems have very complicated equations because
of having 21 independent elastic constants [45]. Thus, we used VASPKIT [46]
code to analyze the elastic properties, and it confirms that SrClF3 is very
likely to be a mechanically stable semiconductor.

(a) The structure of SrClF3

(b) Phonon dispersion curves for SrClF3

Figure 7: Properties of the discovered new materials.

##  4 Discussion

Figure 8: Comparison of the sample distributions of the base line training set
and the AL enhanced dataset using t-sne visualization. It is found that the
samples traversed by the active learning algorithm are much more diverse and
span the majority of the design space.

Our active learning experiments have shown that it can find new materials with
high band gap values as validated by DFT, demonstrating that active learning-
based generative design has the capability to discover new material with
target properties using guided/directed search. Figure 8 shows how active
learning explores the inorganic material space. Active-learning probed samples
( circles) are obviously more widely distributed than the original training
set (+). However, there are several aspects to improve for our active
generative design framework.

First, currently, we are using Bayesian Optimization as the Efficient Global
Optimization (EGO) strategy to search the 128-dimension design space
containing the entire inorganic materials. Yet, many latent representations
sampled during active learning search cannot be decoded into valid chemical
formulas, which has caused a significant waste of the computational budget and
decreases the efficiency of the Bayesian Optimization sampler. Additional
physicochemical constraints need to be incorporated into the sampling strategy
so as to achieve more efficient sampling. Second, we encountered the out-of-
distribution prediction issue of standard machine learning: our random
forest/Roost models cannot predict band gap values that are larger than the
maximum band gaps in the training set [47, 48]. Surrogate ML models with
higher generalization capability are needed to future improve the prediction
accuracy of our screening ML models.

Another key observation of our active learning-based screening experiments is
that the machine learning models such as the graph neural network model Roost
model do not perform too well on the whole dataset since the dataset is highly
biased (too many data points around 0 eV). However, the ML screening model
trained by the augmented samples is capable of identifying high/wide band gap
materials during the screening process. The reason is that the suggested
samples by the active learning tend to focus on high-bandgap materials.

We also tried to directly use active learning to search the whole design space
for discovering high band gap materials, which indeed has found several new
materials such as RbSr2ClF4 and four others. Compared to previous studies
which focus on screening known materials databases, this approach has the
advantage of search in the complete chemical space. However, we find that many
candidates generated via this approach by the active learning model do not
satisfy basic physical requirements such as neutrality and balanced electric
negativity. However, these suggested samples are useful for building an ML-
based screening model as they have effectively explored the design space. The
resulting model can then be applied to the generated samples by MATGAN. So for
semiconductor materials design, we use the MATGAN+screening approach, which is
more suitable to discover semiconductors since the intuition of optimization
is to find the maximum or minimum of a black-box function.

##  5 Conclusion

We propose an active learning-based generative design approach for the
discovery of new inorganic materials and apply it to the discovery of high
band gap materials and semiconductor materials. Our active learning-based
screening is unique due to its search in the whole design space rather than
known databases as done in the previous studies. Our active generative design
(AGD) framework can perform an efficient guided search in the design space to
identify informative samples for building screening ML models with higher
accuracy compared to randomly sampled points. The ML screening model is much
more computationally efficient compared to the conventional DFT-based
screening model. We have successfully applied our AGD approach to the inverse
design of materials with high band gap value and semiconductors and found
several potential new materials, including SrYF5, RbSr2ClF4, and seven other
semiconductors as validated using DFT based principles. Our approach is thus
expected to be applicable to finding materials with other properties such as
thermal conductivity or superionic conductors.

##  6 Contribution

Conceptualization, J.H.; methodology, R.X. and J.H.; software, R.X.;
validation, R.X., J.H.; investigation, R.X, J.H., D.S., Y.S., Y.Z., A.N., S.L.
; resources, J.H.; writing–original draft preparation, R.X., J.H.;
writing–review and editing, R.X. and J.H; visualization, R.X. and Y.S.;
supervision, J.H.; funding acquisition, J.H.

##  7 Acknowledgement

The research reported in this work was supported in part by NSF under the
grant and 1940099 and 1905775 and by NSF SC EPSCoR Program under award number
(NSF Award OIA-1655740 and GEAR-CRP 19-GC02). The views, perspectives, and
content do not necessarily represent the official views of the National
Science Foundation.

## References

  * [1] George S Fanourgakis, Konstantinos Gkagkas, Emmanuel Tylianakis, and George E Froudakis.  A universal machine learning algorithm for large-scale screening of materials.  Journal of the American Chemical Society, 142(8):3814–3822, 2020\.
  * [2] Seung Hwae Heo, Seungki Jo, Hyo Seok Kim, Garam Choi, Jae Yong Song, Jun-Yun Kang, No-Jin Park, Hyeong Woo Ban, Fredrick Kim, Hyewon Jeong, et al.  Composition change-driven texturing and doping in solution-processed snse thermoelectric thin films.  Nature communications, 10(1):1–10, 2019.
  * [3] Artem R Oganov, Chris J Pickard, Qiang Zhu, and Richard J Needs.  Structure prediction drives materials discovery.  Nature Reviews Materials, 4(5):331–348, 2019.
  * [4] Juhwan Noh, Geun Ho Gu, Sungwon Kim, and Yousung Jung.  Machine-enabled inverse design of inorganic solid materials: promises and challenges.  Chemical Science, 11(19):4871–4881, 2020.
  * [5] Yabo Dan, Yong Zhao, Xiang Li, Shaobo Li, Ming Hu, and Jianjun Hu.  Generative adversarial networks (gan) based efficient sampling of chemical composition space for inverse design of inorganic materials.  npj Computational Materials, 6(1):1–7, 2020.
  * [6] Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, et al.  Commentary: The materials project: A materials genome approach to accelerating materials innovation.  Apl Materials, 1(1):011002, 2013.
  * [7] Scott Kirklin, James E Saal, Bryce Meredig, Alex Thompson, Jeff W Doak, Muratahan Aykol, Stephan Rühl, and Chris Wolverton.  The open quantum materials database (oqmd): assessing the accuracy of dft formation energies.  npj Computational Materials, 1(1):1–15, 2015.
  * [8] Kyoungmin Min and Eunseog Cho.  Accelerated discovery of novel inorganic materials with desired properties using active learning.  The Journal of Physical Chemistry C, 2020.
  * [9] Alex Zunger.  Inverse design in search of materials with target functionalities.  Nature Reviews Chemistry, 2(4):1–16, 2018.
  * [10] Piyush M Tagade, Shashishekar P Adiga, Shanthi Pandian, Min Sik Park, Krishnan S Hariharan, and Subramanya Mayya Kolake.  Attribute driven inverse materials design using deep learning bayesian framework.  npj Computational Materials, 5(1):1–14, 2019.
  * [11] Chun-Teh Chen and Grace X Gu.  Generative deep neural networks for inverse materials design using backpropagation and active learning.  Advanced Science, 7(5):1902607, 2020.
  * [12] Juhwan Noh, Jaehoon Kim, Helge S Stein, Benjamin Sanchez-Lengeling, John M Gregoire, Alan Aspuru-Guzik, and Yousung Jung.  Inverse design of solid-state materials via a continuous representation.  Matter, 1(5):1370–1384, 2019.
  * [13] Baekjun Kim, Sangwon Lee, and Jihan Kim.  Inverse design of porous materials using artificial neural networks.  Science advances, 6(1):eaax9324, 2020.
  * [14] Kyungdoc Kim, Seokho Kang, Jiho Yoo, Youngchun Kwon, Youngmin Nam, Dongseon Lee, Inkoo Kim, Youn-Suk Choi, Yongsik Jung, Sangmo Kim, et al.  Deep-learning-based inverse design model for intelligent discovery of organic molecules.  npj Computational Materials, 4(1):1–7, 2018.
  * [15] Mingkun Chen, Jiaqi Jiang, and Jonathan A Fan.  Design space reparameterization enforces hard geometric constraints in inverse-designed nanophotonic devices.  ACS Photonics.
  * [16] Jiaqi Jiang and Jonathan A Fan.  Multiobjective and categorical global optimization of photonic structures based on resnet generative neural networks.  Nanophotonics, 1(ahead-of-print), 2020.
  * [17] Yunwei Zhang, Hui Wang, Yanchao Wang, Lijun Zhang, and Yanming Ma.  Computer-assisted inverse design of inorganic electrides.  Physical Review X, 7(1):011017, 2017.
  * [18] Mahmoud MR Elsawy, Stéphane Lanteri, Régis Duvigneau, Jonathan A Fan, and Patrice Genevet.  Numerical optimization methods for metasurfaces.  Laser & Photonics Reviews, 14(10):1900445, 2020.
  * [19] Ruihao Yuan, Zhen Liu, Prasanna V Balachandran, Deqing Xue, Yumei Zhou, Xiangdong Ding, Jun Sun, Dezhen Xue, and Turab Lookman.  Accelerated discovery of large electrostrains in batio3-based piezoelectrics using active learning.  Advanced materials, 30(7):1702884, 2018.
  * [20] Turab Lookman, Prasanna V Balachandran, Dezhen Xue, and Ruihao Yuan.  Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design.  npj Computational Materials, 5(1):1–17, 2019.
  * [21] Burr Settles.  Active learning literature survey.  Technical report, University of Wisconsin-Madison Department of Computer Sciences, 2009.
  * [22] Christopher Schröder and Andreas Niekler.  A survey of active learning for text classification using deep neural networks.  arXiv preprint arXiv:2008.07267, 2020.
  * [23] Michael Lindenbaum, Shaul Markovitch, and Dmitry Rusakov.  Selective sampling for nearest neighbor classifiers.  Machine learning, 54(2):125–152, 2004.
  * [24] Prem Melville and Raymond J Mooney.  Diverse ensembles for active learning.  In Proceedings of the twenty-first international conference on Machine learning, page 74, 2004.
  * [25] Zahed Allahyari and Artem R Oganov.  Coevolutionary search for optimal materials in the space of all possible compounds.  npj Computational Materials, 6(1):1–10, 2020.
  * [26] Steph-Yves Louis, Yong Zhao, Alireza Nasiri, Xiran Wang, Yuqi Song, Fei Liu, and Jianjun Hu.  Graph convolutional neural networks with global attention for improved materials property prediction.  Physical Chemistry Chemical Physics, 22(32):18141–18148, 2020.
  * [27] Rhys EA Goodall and Alpha A Lee.  Predicting materials properties without crystal structure: Deep representation learning from stoichiometry.  arXiv preprint arXiv:1910.00617, 2019.
  * [28] Vahe Tshitoyan, John Dagdelen, Leigh Weston, Alexander Dunn, Ziqin Rong, Olga Kononova, Kristin A Persson, Gerbrand Ceder, and Anubhav Jain.  Unsupervised word embeddings capture latent knowledge from materials science literature.  Nature, 571(7763):95–98, 2019.
  * [29] Tian Xie and Jeffrey C Grossman.  Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties.  Physical review letters, 120(14):145301, 2018.
  * [30] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  Advances in neural information processing systems, 25:2951–2959, 2012.
  * [31] Jonas Mockus, Vytautas Tiesis, and Antanas Zilinskas.  The application of bayesian methods for seeking the extremum.  Towards global optimization, 2(117-129):2, 1978.
  * [32] Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger.  Gaussian process optimization in the bandit setting: No regret and experimental design.  arXiv preprint arXiv:0912.3995, 2009.
  * [33] F Nogueira.  Bayesian optimization: Open source constrained global optimization tool for python.  URL https://github. com/fmfn/BayesianOptimization, 2014.
  * [34] Cameron Hargreaves, Matthew Dyer, Michael Gaultois, Vitaliy Kurlin, and Matthew J Rosseinsky.  The earth mover’s distance as a metric for the space of inorganic compositions.  2020\.
  * [35] G. Kresse and J. Hafner.  ab initio.  Phys. Rev. B, 47:558–561, Jan 1993.
  * [36] G. Kresse and J. Hafner.  ab initio.  Phys. Rev. B, 49:14251–14269, May 1994.
  * [37] J. Furthmüller G. Kresse.  Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set.  Comput. Mater. Sci., 6:15–50, jul 1996.
  * [38] G. Kresse and J. Furthmüller.  Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.  Phys. Rev. B, 54:11169–11186, Oct 1996.
  * [39] P. E. Blöchl.  Projector augmented-wave method.  Phys. Rev. B, 50:17953–17979, Dec 1994.
  * [40] G. Kresse and D. Joubert.  From ultrasoft pseudopotentials to the projector augmented-wave method.  Phys. Rev. B, 59:1758–1775, Jan 1999.
  * [41] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple.  Phys. Rev. Lett., 77:3865–3868, Oct 1996.
  * [42] John P. Perdew, Kieron Burke, and Matthias Ernzerhof.  Generalized gradient approximation made simple [phys. rev. lett. 77, 3865 (1996)].  Phys. Rev. Lett., 78:1396–1396, Feb 1997.
  * [43] Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geoffroy Hautier, Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent L. Chevrier, Kristin A. Persson, and Gerbrand Ceder.  Python materials genomics (pymatgen): A robust, open-source python library for materials analysis.  Computational Materials Science, 68:314–319, 2013.
  * [44] R Hill.  The elastic behaviour of a crystalline aggregate.  Proceedings of the Physical Society. Section A, 65(5):349–354, may 1952.
  * [45] Félix Mouhat and François-Xavier Coudert.  Necessary and sufficient elastic stability conditions in various crystal systems.  Physical Review B, 90, 09 2014.
  * [46] V. Wang, N. Xu, J. Liu, Gang Tang, and W. Geng.  Vaspkit: A user-friendly interface facilitating high-throughput computing and analysis using vasp code.  2019\.
  * [47] Zheng Xiong, Yuxin Cui, Zhonghao Liu, Yong Zhao, Ming Hu, and Jianjun Hu.  Evaluating explorative prediction power of machine learning algorithms for materials discovery using k-fold forward cross-validation.  Computational Materials Science, 171:109203, 2020.
  * [48] Christian Loftis, Kunpeng Yuan, Yong Zhao, Ming Hu, and Jianjun Hu.  Lattice thermal conductivity prediction using symbolic regression and machine learning.  The Journal of Physical Chemistry A.