The non-serendipitous discovery of materials with targeted properties is the ultimate goal of materials research, but to date, materials design lacks the incorporation of all available knowledge to plan the synthesis of the next material. This work presents a framework for learning a continuous representation of materials and building a model for new discovery using latent space representation. The ability of autoencoders to generate experimental materials is demonstrated with vanadium oxides via rediscovery of experimentally known structures when the model was trained without them. Approximately 20,000 hypothetical materials are generated, leading to several completely new metastable V_xO_y materials that may be synthesizable. Comparison with genetic algorithms suggests computational efficiency of generative models that can explore chemical compositional space effectively by learning the distributions of known materials for crystal structure prediction. These results are an important step toward machine-learned inverse design of inorganic functional materials using generative models

Aspuru-Guzik, Alán

Gregoire, John M.

Jung, Yousung

Kim, Jaehoon

Noh, Juhwan

Sanchez-Lengeling, Benjamin

Stein, Helge S.

Caltech Authors - Main

Matter, Volume 1Supplemental InformationInverse Design of Solid-State Materials
via a Continuous Representation
Juhwan Noh, Jaehoon Kim, Helge S. Stein, Benjamin Sanchez-Lengeling, John M.
Gregoire, Alan Aspuru-Guzik, and Yousung Jung
1 
 
S1. Machine learning model details 
 
Figure S1. The proposed hierarchical two-step Image-based Materials Generator (iMatGen); (a) image 
compression for unit cell and (b) basis, and (c) materials generator (MG) using a variational autoencoder 
(VAE) model. 
 
The structure of the 3D convolutional autoencoder was inspired from the 3D GAN introduced 
by Wu and coauthors1, and details are described in S1.1 and S1.2 where all machine learning 
models were developed using Tensorflow2, an open-source machine learning framework using 
Python. 
 
S1.1 Cell convolutional autoencoder (CELL) 
For the CELL (Figure S1a), the encoder consists of four 3D convolutional layers with kernel 
sizes {4,4,4,4}, strides {2S,2S,2S,1V} and number of channels {64,64,64,20}. For stride, S 
and V mean padding option SAME and VALID used in Tensorflow, respectively. We used a 
leaky-ReLU activation function with a parameter of 0.2 between convolutional layers and 
hyperbolic tangent function at the end layer. The input is a 32×32×32 tensor, and the output (or 
2 
 
the latent vector, 𝑣𝑐𝑒𝑙𝑙) is a 20-dimensional vector the value of which is in [-1, 1].  
The decoder has a mirrored structure to that of encoder where the kernel sizes are {4,4,4,4}, 
the strides are {1V,2S,2S,2S} and the number of channels are {64,64,64,1}. We used the leaky-
ReLU activation function with a parameter of 0.2 between convolutional layers and sigmoid 
function at the end. The decoder takes the input a 20-dimensional vector (𝑣𝑐𝑒𝑙𝑙), and the output 
a 32×32×32 tensor the value of which is in [0, 1].  
For training of CELL, we adopted mean-squared-error loss function (or L2-loss function). 
 
S1.2 Basis convolutional autoencoder (BASIS) 
For the BASIS (Figure S1b), the encoder consists of five 3D convolutional layers with kernel 
sizes {4,4,4,4,4}, strides {2S,2S,2S,2S,1V} and number of channels {64,64,64,64,200}. We 
used a leaky-ReLU activation function with parameter of 0.2 between convolutional layers, 
and hyperbolic tangent function at the end. The input is a 64×64×64 tensor, and the output (or 
the latent vector, 𝑣𝑏𝑎𝑠𝑖𝑠) is a 500-dimensional vector the value of which is in [-1, 1]. 
The decoder has a mirrored structure to that of encoder where kernel sizes {4,4,4,4,4}, strides 
{1V,2S,2S,2S,2S} and number of channels {64,64,64,64,1}. We used a leaky-ReLU activation 
function with the amount of leakage is 0.2 between convolutional layers, and sigmoid function 
at the end. The decoder takes the input a 200-dimensional vector (𝑣𝑏𝑎𝑠𝑖𝑠), and the output is a 
64×64×64 tensor the value of which is in [0, 1]. 
For training of BASIS, we adopted mean-squared-error loss function (or L2-loss function). 
 
S1.3 Materials Generator (MG) – Variational AutoEncoder (VAE) 
The architecture of VAE (Figure S1c) is following. The encoder consists of four 2D 
convolutional layers with kernel sizes {4,4,4,1}, strides {2S,2S,2S,1S} and number of channels 
{100,100,100,50}, and we used leaky-ReLU activation function with parameter 0.2 between 
convolutional layers. The input is a 1×200×6 tensor, and the output is a 1250-dimensional 
vector. The output of the convolutional layers are transformed to a 500-dimensional vector (𝑧) 
by fully connected layers. Between the last convolutional layer and fully connected layers, a 
leaky ReLU with parameter 0.2 was used.  
The decoder takes a 500-dimensional vector (𝑧 ) as input which is transformed to a 1250-
dimensional vector by a fully connected layer with leaky ReLU with parameter 0.2. Next, four 
2D convolutional layers were followed with kernel sizes {1,4,4,4}, strides {1S,2S,2S,2S} and 
number of channels {100,100,100,6}. A leaky ReLU with parameter of 0.2 was used between 
convolutional layers, and hyperbolic tangent activation function at the end. 
For the formation energy classification task from the latent space, we used a single hidden layer 
neural network where the input is a 500-dimensional vector (𝑧) and the number of node in 
hidden layer is 500. A Hyperbolic tangent function was used between layers and sigmoid 
function was used at the end layer. 
For training of VAE, we adopted following loss function, ℒ𝑉𝐴𝐸: 
3 
 
 
ℒ𝑉𝐴𝐸 = ℒ𝑟𝑒𝑐𝑜𝑛. + 𝛼ℒ𝐾𝐿 + 𝛽ℒ𝑐𝑙𝑎𝑠𝑠., 𝑧~𝛿 ∗ 𝑁(0, 𝐼) 
where ℒ𝑟𝑒𝑐𝑜𝑛. is the reconstruction loss term (L2-loss) and ℒ𝐾𝐿 is the KL divergence term to 
regularize the latent space (𝑧 ) to the normal distribution. Furthermore, ℒ𝑐𝑙𝑎𝑠𝑠. represents 
classification loss where we used the cross entropy binary classification loss with the weighting 
parameter 𝛽 to control learning speed, similar to 𝛼 for the ℒ𝐾𝐿. 
 
S1.4 Training notes 
For all training procedures, we randomly sampled 90% of the total data as a training set. Before 
training iMatGen with the VO dataset, we trained iMatGen using the MP subset consisting of 
28,242 crystal structures taken from the MP database (up to 5 atom types in the unit cell) 
satisfying two constraints as mentioned in the main text to obtain pre-trained model parameters 
for IC and MG (without formation energy classification task). Training of iMatGen with VO 
dataset was done by using the latter pre-trained model parameters, and as mentioned in the 
main text, formation energy classification task was also additionally included. All the used 
hyper-parameters including batch size, learning rate and optimizer are listed in Table S1. 
 
Table S1. Lists of the all hyper-parameters used in this work 
 
 
S1.5 Latent space of the trained MG 
Since the normal distribution is assumed as prior distribution of the latent space, we checked 
that the trained latent space is shaped appropriately. As shown in Figure S2, the mean and 
standard deviation of the latent space taken from the trained MG are well matched to the normal 
distribution. 
 
4 
 
 
Figure S2. Mean and standard deviation distribution of the 500-latent vector elements after training the 
Materials Generator (MG) 
 
 
Figure S3. Effect of including stability classification task from the latent space of the Materials 
Generator (MG) 
 
 
Figure S4. Statistical analysis of the trained latent space with respect to the stability and structural 
diversity 
  
5 
 
S2. DFT calculation details 
S2.1 DFT calculations for VO dataset  
We used the PBE functional3 and the PAW4-PBE pseudopotentials for the V and O atoms as 
implemented in the program package VASP5. Structures (both atomic positions and cell 
parameters) were fully relaxed using the conjugate gradient descent method with convergence 
criteria of 1.0e-5 for energy and 0.05 eV/Å for force. For computational efficiency, we used 
relatively sparse reciprocal space meshes (0.5 Å−1 as the grid spacing) and the cut-off energy 
for the plane-wave expansions of 500 eV since the main goal of constructing the VO dataset 
was collecting diverse atomic configurations. The formation energy (or formation enthalpy, 𝐸𝑓) 
was calculated using EVxOy − (xEV − yEO)/(x + y) (in eV/atom). 
 
S2.2 DFT calculations for the generated materials 
For all generated structures (~20,000 materials), we performed spin polarized GGA+U 
calculations, with the same U parameter for V used in the Materials Projects database6. We 
relaxed both atomic positions and cell parameters using conjugate gradient descent method 
with convergence criteria of 1.0e-5 for energy and 0.05 eV/Å for force with 500 eV cut-off 
energy. To compare the phase stability among the generated structures for all generated 
materials we first used sparse reciprocal lattice grid with a grid spacing of 0.5 Å−1 (See Figure 
4 and Figure 6). Then, for a smaller set of materials that satisfy 𝐸hull ≤ 0.1 eV/atom for 
Slerp and 𝐸hull ≤ 0.2 eV/atom for Random), we refined the formation energy calculations 
using a denser reciprocal lattice grid with grid spacing of 0.25 Å−1 (Figure 7). 
 
S2.3 Phase Diagram using GGA/GGA+U-mixed approach 
As mentioned in the main text, we followed the approach proposed by Anubhav Jain et al.7 
which used GGA/GGA+U-mixed approach for the formation energy calculation especially for 
transition metal oxides. In the latter approach, the energy of all GGA+U calculations is 
corrected by equation S1. 
 
𝐸𝑉,𝑂 compound
GGA+U renorm. = 𝐸𝑉,𝑂 compound
GGA+U − 𝑛𝑉∆𝐸𝑉  (S1) 
 
In equation (S1), 𝐸𝑉,𝑂 compound
GGA+U  is the GGA+U energy for V, O compound and 𝑛𝑉 is the 
number of 𝑉 atoms in the compound. Here, our target is to obtain a correction term ∆𝐸𝑉 for 
𝑉. In addition, we included an O2 energy correction term (equation S2) taken from the MP 
database for the 𝐸𝑉,𝑂 compound
GGA+U  energy. 
 
𝐸(𝑂2)𝑐𝑜𝑟𝑟. = −1.4046𝑒𝑉/𝑎𝑡𝑜𝑚    (S2) 
6 
 
 
After a simple rearrangement using equation (S1), the correction energy for 𝑉 atom is then 
obtained using equation (S3). 
 
∆𝐸𝑉 = (𝐸𝑓
ref. − Calc. 𝐸𝑓) 𝑟𝑉⁄    (S3) 
 
In equation (S3), 𝐸𝑓
ref. is the formation energy for the 5-reference structures (V2O5, V3O7, VO2, 
V3O5 and V2O3) taken from the MP database all of which are at the ground state for the 
corresponding composition, Calc. 𝐸𝑓 is the formation energy calculated from this work, and 
𝑟𝑉 is the fraction of 𝑉 atom in a compound. Therefore, ∆𝐸𝑉 can be interpreted as the slope 
of the (𝐸𝑓
ref. − Calc. 𝐸𝑓) vs. 𝑟𝑉 plot as shown in Figure S5a. We calculated ∆𝐸𝑉 for two 
cases of reciprocal lattice grid space (0.5 Å−1 and 0.25 Å−1). The value of ∆𝐸𝑉 for each case 
is shown in Figure S5a, and the adjusted formation energy values are described in Figure S5b. 
 
 
Figure S5. Results of GGA/GGA+U mixing scheme applied in this work: (a) Deriving correction 
energy for V (∆𝐸𝑉) using the GGA/GGA+U-mixed scheme for the 5-VxOy materials at the convex hull, 
and (b) results after applying the latter correction energy terms for the 5-VxOy materials at the convex 
hull (in eV/atom). 
 
7 
 
As shown in the Figure S5b, we find that all the formation energies after correction are well 
matched to the reference values. For further validation, we next applied the same correction 
schemes to calculate phase stability (formation energy and energy above the convex hull) of 
all 112 vanadium oxide structures existing in MP database. We find a good agreement between 
the present results and those from the MP database as shown in Figure S6. 
 
 
Figure S6. Comparison of the calculated phase stability values for 112-VxOy materials existing in MP 
database: (a) formation energy (𝐸𝑓) and (b) energy above the convex hull (𝐸ℎ𝑢𝑙𝑙). 
 
  
8 
 
S3. Materials Generation 
S3.1 Latent space sampling 
In addition to Gaussian random sampling of 10,000 vectors from the latent space, we 
additionally sampled the latent vector by interpolating between the known VxOy structures in 
MP database with the assumption that the interpolated points between two existing structures 
can generate more reasonable structures. Here, we applied the spherical linear interpolation 
(slerp, (equation S4)): 
 
𝑧𝑚 =
sin(t∙Ω)
sin Ω
𝑧𝑠𝑡𝑎𝑟𝑡 +
sin(1−t)∙Ω
sin Ω
𝑧𝑒𝑛𝑑, 𝑡 = 𝑚/𝑁 (S4) 
 
where 𝑧𝑚 is the interpolated vector, Ω is the angle between two latent vectors (𝑧𝑠𝑡𝑎𝑟𝑡 and 
𝑧𝑒𝑛𝑑). 𝑡 is the interpolation grid defined as 𝑚/𝑁, where the maximum number of divisions 
(𝑁) is 20 and 𝑚 is positive integer in [1,19] (𝑚 = 0 is corresponding to 𝑧𝑠𝑡𝑎𝑟𝑡 and 𝑚 = 20 
is corresponding to 𝑧𝑒𝑛𝑑). 
 
S3.2 Inverse transformation of the input representation  
 
 
Figure S7. Description of unit cell reconstruction from input representation: Herein, we show 
reconstruction of ab-plane as an example. 
 
9 
 
 
Figure S8. Dissimilarity value for the reconstructed materials with respect to the reference materials as 
shown in Fig 4 of main text. 
 
 
Figure S9. Error histograms for the reconstruction scheme applied to the basis images. 
 
 
Figure S10. Error histograms for the inverse transform scheme applied to the cell images 
 
S3.3 Post-process for the generated images 
In post-processing of the generated structures, we first removed the unfavorable position 
10 
 
overlap within a single image under periodic boundary conditions (PBC). In other words, we 
imposed the PBC of each generated image in a post-hoc manner by merging the V atoms within 
the covalent radius (1.6 Å) as a single V atom at the unit cell boundary, and likewise for the 
oxygen images using the 1.3 Å criterion. We used the same criteria to handle two atoms of the 
same atom type in close proximity inside the unit cell. For structures in which V and O are in 
close proximity, we used 1.2 Å criterion to remove oxygen atom in a way to meet the valid 
oxidation state of V. 
 
S3.4 Materials generation from the iMatGen – numerical statistics 
Figure S11 shows the distribution of all generated VxOy compositions within valid 𝑂𝑆𝑉 range 
as shown in Figure 5 of the main text, and 4-different sections are also described in the inset of 
Figure S11.  
 
Figure S11. All numerical results for the generated materials within the valid 𝑂𝑆𝑉 range. Each section 
is depicted in the inset coordination figure. 
  
11 
 
References 
1 Wu, J., Zhang, C., Xue, T., Freeman, B., and Tenenbaum, J. (2016). Learning a probabilistic latent 
space of object shapes via 3d generative-adversarial modeling. Paper presented at: Advances in Neural 
Information Processing Systems. 
2 Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., 
and Isard, M. (2016). Tensorflow: a system for large-scale machine learning. Paper presented at: OSDI. 
3 Perdew, J.P., Burke, K., and Ernzerhof, M. (1996). Generalized gradient approximation made simple. 
Physical review letters 77, 3865. 
4 Blöchl, P.E. (1994). Projector augmented-wave method. Physical Review B 50, 17953. 
5 Kresse, G., and Furthmüller, J. (1996). Software VASP, vienna (1999). Phys Rev B 54, 169. 
6 Jain, A., Ong, S.P., Hautier, G., Chen, W., Richards, W.D., Dacek, S., Cholia, S., Gunter, D., Skinner, 
D., and Ceder, G. (2013). Commentary: The Materials Project: A materials genome approach to accelerating 
materials innovation. Apl Materials 1, 011002. 
7 Jain, A., Hautier, G., Ong, S.P., Moore, C.J., Fischer, C.C., Persson, K.A., and Ceder, G. (2011). 
Formation enthalpies by mixing GGA and GGA+Ucalculations. Physical Review B 84. 
 


Inverse Design of Solid-State Materials via a Continuous Representation

https://authors.library.caltech.edu/99012/2/1-s2.0-S2590238519301754-mmc1.pdf

Inverse Design of Solid-State Materials via a Continuous Representation

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main