Materials Properties Prediction (MAPP):
Empowering the prediction of material properties solely based on chemical formulas

The datasets for both the melting temperature and heat of fusion were collected from a ten-volume compilation of thermodynamic constants of substances, the book “Thermodynamic properties of individual substances” glushko1988thermodynamic , which is openly available in database format russian_thermal_property_website .

The melting temperature dataset contains 9375 materials, among which 982 compounds exhibit high melting temperatures with melting points exceeding 2,000K. While the majority of the dataset is derived from experimental findings, a small portion originates from DFT-based data generated by our first-principle calculation tool, SLUSCHI Hong2012 ; hong_user_2016 . The dataset underlying this model predominantly consists of compounds with congruent melting temperatures. Hence, the melting temperature generated by this model is interpreted as the higher end of the solidus-liquidus temperature range. This interpretation similarly applies to compounds that decompose prior to melting. It is crucial to understand that this model does not predict solidus or liquidus temperatures.

Regarding the heat of fusion, the dataset comprises 774 data points. Upon inspecting the dataset, we identified an anomaly in the material $\mathrm{GdPd_{3}}$, leading to its manual removal from the dataset.

The bulk modulus and volume data were queried from the Materials Project database jain2013commentary , utilizing the Pymatgen package ong_python_2013 . To ensure data consistency, we then filtered out less stable structures, applying a threshold of 10 meV energy above the convex hull. Consequently, the bulk modulus dataset contains 4,236 entries, while the volume dataset consists of 49,213 entries. It is noteworthy that, despite materials initially querying spanning all dimensions (0D, 1D, 2D, and 3D materials), our finalized dataset is restricted to materials with 3D crystal structures, given that only these materials provide well-defined volume values. We note that the bulk modulus and volume data are based on DFT calculations at absolute zero.

In order to develop a universal model for predicting the superconducting critical temperature $T_{c}$ with a wide range of $T_{c}$ value, we utilized a comprehensive dataset from a previous study hamidieh_data-driven_2018 . This dataset encompasses a wide range of superconductors, including cuprate-based, iron-based, and Bardeen-Cooper-Schrieffer (BCS) theory superconductors, which originate from the Supercon material Database, maintained by the National Institute for Materials Science (NIMS) in Japan Supercon_dataset . This database is the most extensive and widely employed resource for data-driven research on superconductors, and it has been extensively utilized in previous studies stanev_machine_2018 ; hamidieh_data-driven_2018 ; zeng_atom_2019 ; konno_deep_2021 .

In our approach, the material’s chemical formula is depicted as a fully connected element graph $G_{e}$. Each element graph $G_{e}$ comprises nodes and edges, defined as $V$ and $E$ respectively. Figure 1 (a) shows the process of converting the chemical formula into an element graph. To exemplify, the material $\mathrm{Li_{7}Mn_{4}CoO_{12}}$ is visualized as an element graph with nodes corresponding to elements Lithium, Manganese, Cobalt, and Oxygen. Each node represents a specific element, connected to each neighboring node through a single edge. Each edge symbolizes the path for information exchange between paired elements. Specifically, the edges have no features and are treated equally, only indicating the connection of the information-exchanging route. In constructing the element graph, each element within the chemical formula is associated with a node feature vector. Each node feature vector has 14-dimensional features, composed of elemental properties, such as atomic mass, atomic number, melting temperature, boiling temperature, electronegativity, etc, as well as the composition of each element. The composition serves as an indicator of the relative importance of each element within the chemical formula of the material.

The node feature vector is denoted as $x_{v}$, with $v$ symbolizing a specific node within the node-set, $V$. To ensure the permutation invariance for all elements in the subsequent graph neural network section, each $x_{v}$ undergoes nonlinear transformation through an identical fully connected neural network, utilizing the same activation function. This transformation is expressed as: $h_{v}^{0}=\mathrm{ReLU}(W^{0}x_{v}+b^{0}),v\in V$. In this equation, $\mathrm{ReLU}$ denotes the rectified linear unit activation function. The term $h_{v}^{0}$ signifies the initial element embedding prior to the graph neural network layer, while the superscript indicates the specific layer within the graph neural network. Moreover, $W^{0}$ and $b^{0}$ denote the weight and bias parameters of the neural network.

After creating the element graph based on chemical formula, the graph neural network model is used to transform element embeddings and subsequently, the embedding of the whole element graph in order to capture more physical insights from the material and to perform the material property prediction better. Our task is to predict the material property using the element graph, as shown in Figure 1 (b), framing this as a graph-level regression task. Within the graph-level prediction task, the objective of the GNN is to learn a representation of the entire graph, the so-called material embedding $A$ in Figure 1 (b). This is achieved by iteratively updating the node embeddings based on the neighboring information, also termed as the message in literature gilmer2017neural and then aggregating the individual node representations to form the material representation, which is used to perform the regression. In the following paragraphs, the element embedding update process and the material embedding update process will be elaborated.

As shown in Figure 1 (b), element embedding of the target node is updated by incorporating messages from all the neighboring elements. This inference bias arises from the understanding that the information from the neighboring nodes is not only relevant to the target node but also enriches the target node’s information.

There are generally two phases in the element embedding update process, the message forming phase and the message element combining phase, as is shown in Figure 1 (c). The message is created based on both the embedding of the neighboring node and the current node and has the general form shown below:

$$m_{v}^{t-1}=\sum_{u\in N(v)}^{|V|-1}M_{t}\left(h_{v}^{t-1},h_{u}^{t-1}\right),$$

(1)

where $m_{v}^{t-1}$ is the message for target node $v$ at iteration $t-1$, $N(v)$ is the neighboring node set of $v$, $|V|$ is the number of the nodes in the entire graph. $h_{v}^{t-1}$ and $h_{u}^{t-1}$ denote the element embeddings of target node $v$ and neighboring node $u$ at iteration $t-1$. The message function $M_{t}$ can take multiple forms, for example, adding, averaging, concatenating, etc, which all serve as a way to preserve the information of the neighboring node. The target node’s overall message is updated by aggregating messages from all of the neighboring nodes, which preserves the properties of permutation invariance of the chemical formula. In this study, the specific form of message function is

$$m_{v}^{t-1}=\sum_{u\in N(v)}^{|V|-1}\Bigg{[}\mathrm{ReLU}\Bigg{(}W^{t-1}\mathrm{Add}\left(h_{v}^{t-1},h_{u}^{t-1}\right)+b^{t-1}\Bigg{)}\Bigg{]},$$

(2)

where $W^{t-1}$ and $b^{t-1}$ are the weight and bias of the neural network at the $t-1$ iteration. Figure 1 (c) illustrates the message-forming procedures. The messages from elements Co, Mn, and Li are first created by adding the element embeddings of the target node $h_{\mathrm{O}}$ and neighboring node $h_{u}$. Then all the individual messages are passed through an identical neural network layer. An overall message $m_{\mathrm{O}}$ is formed by aggregating all the individual messages.

After generating the message $m_{v}^{t-1}$, $h_{v}$ is updated by summing the message vector $m_{v}^{t-1}$ with the node embedding from previous layer $h_{v}^{t-1}$,

$$h_{v}^{t}=\mathrm{ReLU}\left(W^{t-1}\mathrm{Add}\left(h_{v}^{t-1},m_{v}^{t-1}\right)+b^{t-1}\right).$$

(3)

The material embedding $A^{t}$ for layer $t$ is constructed by aggregating all element embeddings, a process illustrated in Figure 1 (b),

$$A^{t}=\sum_{v\in V}^{|V|}h_{v}^{t}.$$

(4)

The element and material embedding update process will last for $T$ iterations, which is a hyperparameter to tune.

In our architecture, we record the material embeddings $A$ from every iteration, accumulating them in part to construct the final material embedding, $A^{\mathrm{final}}$. In general, $A^{\mathrm{final}}$ may be constituted in two ways. It could either be the material embedding $A^{T}$ derived from the final iteration within the graph neural network, or an aggregation of the material embedding $A$ created during the $T$-th iteration along with those from earlier iterations. This method helps in conserving information from preceding GNN layers, potentially mitigating the prevalent issue of over-smoothing often encountered in GNN applications.

After generating a material embedding through the GNN section, the material embedding is fed into the fully connected layers with ResNet architecture he_deep_2016 for further nonlinear transformation. The ResNet layer has the advantage over the plain dense layer in that it can mitigate the exploding and vanishing gradient issue, therefore, guaranteeing that the model can converge to the optimal solution easier. After $N$ number of ResNet transformations and a regression neural network layer, a final value of the predicted material property is given by the model.

To summarize the aforementioned model architecture and material prediction process, our GNN model’s standard workflow for predicting material properties is depicted in Figure 1 (d), demonstrating the comprehensive, end-to-end capability of the GNN. In this model, the chemical formula, as the sole input, is encoded as a series of elemental embeddings. Within the GNN layers, these elemental embeddings undergo $T$ iterations of updates, aggregating to form an ultimate material embedding. This embedding is subsequently employed for regression within the ResNet block.

To further increase the model accuracy, quantify prediction uncertainty, detect outliers, and perform comprehensive data analysis, we propose an ensemble model based on the bootstrap method. The model consists of 30 independent GNN models trained on different sampling from the original dataset, with a final model performance derived from aggregating the performance metrics of 30 individual models, such as coefficient of determination ($R^{2}$ score), root-mean-square error (RMSE), and mean absolute error (MAE). The training data for each individual model is randomly sampled using the bootstrap method from the original dataset. The testing set is composed of data not included in each sampling process.

The ensemble model can increase the robustness of the individual deep learning model. The diversity of the model can be ensured because all the models are trained on a different subset of the original dataset. This diversity can potentially improve the model’s accuracy since the errors made by different models may be canceled from each other when aggregated. Moreover, the overfitting tends to be averaged out as each individual model in the ensemble is trained on a different subset of the original dataset, which may help the model better generalize to the unseen data.

Additionally, the utilization of an ensemble model allows for the quantification of uncertainty and facilitates a comprehensive analysis of outliers within the dataset. By ensuring the inclusion of every data point in both training and testing sets at least once, we establish a systematic basis for thoroughly assessing the model’s predictive performance for each individual material. This approach guarantees that the training set closely represents the overall distribution of the original dataset, thereby providing a robust evaluation of the model’s performance under various scenarios. By analyzing the statistics of the predicted errors, specifically the mean, median, and maximum differences between predicted and actual values, we can identify materials that exhibit significant deviations in both the training and testing phases. Such anomalies in the data may necessitate a closer examination to determine the causes of such anomalies. These inconsistencies could stem from incorrect data entries, requiring their removal from the original dataset to enhance accuracy. Alternatively, they might represent unique statistical distributions that our current model fails to recognize, indicating a need for more sophisticated data handling or model adjustment to accommodate these exceptions.