Coarse-grained crystal graph neural networks for reticular materials
  design
Reticular materials, including metal-organic frameworks and covalent organic
frameworks, combine relative ease of synthesis and an impressive range of
applications in various fields, from gas storage to biomedicine. Diverse
properties arise from the variation of building units$\unicode{x2013}$metal
centers and organic linkers$\unicode{x2013}$in almost infinite chemical space.
Such variation substantially complicates experimental design and promotes the
use of computational methods. In particular, the most successful artificial
intelligence algorithms for predicting properties of reticular materials are
atomic-level graph neural networks, which optionally incorporate domain
knowledge. Nonetheless, the data-driven inverse design involving these models
suffers from incorporation of irrelevant and redundant features such as full
atomistic graph and network topology. In this study, we propose a new way of
representing materials, aiming to overcome the limitations of existing methods;
the message passing is performed on a coarse-grained crystal graph that
comprises molecular building units. To highlight the merits of our approach, we
assessed predictive performance and energy efficiency of neural networks built
on different materials representations, including composition-based and
crystal-structure-aware models. Coarse-grained crystal graph neural networks
showed decent accuracy at low computational costs, making them a valuable
alternative to omnipresent atomic-level algorithms. Moreover, the presented
models can be successfully integrated into an inverse materials design pipeline
as estimators of the objective function. Overall, the coarse-grained crystal
graph framework is aimed at challenging the prevailing atom-centric perspective
on reticular materials design.
Vadim Korolev, Artem Mitrofanov
# Coarse-grained crystal graph neural networks for reticular materials design

###### Abstract

Reticular materials, including metal-organic frameworks and covalent organic frameworks, combine relative ease of synthesis and impressive range of applications in various fields, from gas storage to biomedicine. Diverse properties arise from the variation of building units, metal centers and organic linkers, in an almost infinite chemical space. Such a variability substantially complicates experimental design and promotes the use of computational methods. In particular, the most successful artificial intelligence algorithms for predicting properties of reticular materials are atomic-level graph neural networks with optional domain knowledge. Nonetheless, the data-driven inverse design utilizing such models suffers from incorporating irrelevant and redundant features such as full atomistic graph and network topology. In this study, we propose a new way of representing materials, aiming to overcome the limitations of existing methods; the message passing is performed on the coarse-grained crystal graph that comprises molecular building units. We assess the predictive performance and energy efficiency of neural networks built on different materials representations, including composition-based and crystal-structure-aware models, to highlight the merits of our approach. Coarse-grained crystal graph neural networks show decent accuracy at low computational costs, making them a valuable alternative to omnipresent atomic-level algorithms. Moreover, the presented models can be successfully integrated into the inverse materials design pipeline as estimators of the objective function. Overall, the coarse-grained crystal graph framework aims to challenge the prevailing atomic-centric perspective on reticular materials design.

_Keywords_ reticular design metal-organic frameworks covalent organic frameworks coarse graining graph neural networks

## 1 Introduction

Artificial intelligence (AI) is becoming the next game changer in materials science[1, 2, 3]. Nowadays, supervised learning algorithms represent a cutting-edge tool for resolving complex structure-property relationships that determine material functionality. AI methods have been utilized to examine a plethora of physicochemical parameters, including thermodynamic[4, 5, 6], electronic[7, 8, 9], mechanical[10, 11, 12], adsorption[13, 14, 15], and catalytic[16, 17, 18] properties. Deep learning models[19] provide state-of-the-art predictive performance by extracting crucial features from input data, e.g., chemical composition and crystal structure. Most representation schemes related to deep learning models in materials informatics require researchers to perceive materials from an atomistic viewpoint. In particular, crystal graphs (where nodes and edges correspond to atoms and chemical bonds, respectively) enhance AI toolset of materials scientists with a diverse set of graph neural networks[20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Structure-agnostic neural networks[30, 31, 32]
handle chemical composition through assigning vector representations (so-called embeddings) to the chemical elements present in the material.

Atomic-level models are highly effective in representing materials if atoms are considered the most natural choice for the elementary structural unit. In the inorganic domain, chemical subspaces can be formed through element substitution in the prototype crystal structure; the resulting material sets have been processed with predictive models in a high-throughput manner[33, 34, 35]. In contrast, materials representation at a scale of molecular building units is preferable for crystalline extended structures governed by reticular chemistry principles[36, 37]. The most prominent examples of such materials are metal-organic frameworks[38] (MOFs) and covalent organic frameworks[39] (COFs). Specifically, MOFs are formed by linking organic molecules and metal-containing entities through coordination bonds, whereas COFs are made by stitching organic molecules through covalent bonds. Some of synthesized compounds show exceptional adsorption[40, 41] and catalytic[42, 43] properties, making them promising candidates for energy-related applications. Moreover, the modular structure of reticular materials provides a great opportunity for further tuning of relevant properties. There have been recent efforts to develop specialized featurization schemes for reticular design[44, 45, 46, 47, 48]. Global descriptors (e.g., topology, volumetric attributes, and energy grids) incorporated into neural network architecture improve predictive performance and leaves room for interpretability analysis. On the other hand, the aforementioned attributes constrain the scenarios where the corresponding models can be applied. The most intriguing data-driven strategy for developing functional materials--inverse design (i.e., "from property to structure")--requires the synchrony between discriminative and generative models. In particular, the input data modalities (i.e., materials representation _in toto_) used in the former models should be validly reproduced by the latter ones. Given this context, the incorporation of energy-grid embeddings in predictive models (e.g., MOFTransformer[47]) presents a challenge for inverse materials design owing to the limited reconstruction ability of existing algorithms. Out of one million structures, only a few zeolite shapes created by generative adversarial network (ZeoGAN[49]) successfully passed all cleanup operations. Other predictive models[48] for reticular materials directly consider framework topology, which is unknown _a priori_ for specified building blocks, e.g., organic linkers and metal-containing unit in MOFs. Consequently, the unpredictable synthetic accessibility of frameworks (in the form of the linkers-nodes-topology triad) created by generative models has become a cornerstone for practical data-driven design of reticular materials. Recent findings have confirmed concerns that most relevant studies have overlooked: only 136 frameworks out of one thousand hypothetical structures with the highest hydrogen working capacity have been identified as highly synthesizable[50].

By viewing experimental synthesis and characterization of _in silico_ generated frameworks as the ultimate goal of AI-assisted materials discovery, we can formulate the main challenge in the field as follows: there exists a fundamental disparity in how synthetic chemists and AI-practitioners perceive reticular materials. As a result, atomic-level and topology-aware models leverage materials modalities that are not relevant or accessible in reticular design. To address this issue, we introduce the coarse-grained crystal graph framework. The performance of neural networks that utilize sets of molecular building units as input is analyzed in terms of accuracy, energy efficiency, and transferability. In our comparative analysis, we also take into account widely used architectures, which are the dominant paradigm in materials representation learning; composition-based and crystal-structure-aware models are considered. In addition to material screening applications, the feasibility of incorporating coarse-grained crystal graph neural networks into the inverse reticular design pipeline is explored as well.

## 2 Results

### Coarse-grained crystal graph framework

The modern landscape of predictive models in materials informatics[51] is mostly shaped by neural networks built on crystal graph \((\mathcal{V},\mathcal{E})\), which is a set of vertices (atoms) \(v\in\mathcal{V}\) and edges (bonds) \((u\to v)\in\mathcal{E}\). Despite an impressive diversity of neural network architectures, most of them have a common foundation in the message-passing paradigm[52]. The vector representation \(h_{v}\) of a node \(v\) is generated by propagating messages \(m_{u\to v}\) from source nodes \(u\) to a destination node \(v\); the messages from all the nodes forming a neighborhood \(\mathcal{N}(v)\) of node \(v\) are taken into account. The message-passing procedure may be performed multiple times, so a node representation \(h_{v}^{t+1}\) at step \(t+1\) is computed as follows:

\[m_{u\to v}^{t+1}=\phi(h_{u}^{t},h_{v}^{t},e_{u\to v}^{t}) \tag{1}\]

\[m_{v}^{t+1}=\rho(\big{\{}m_{u\to v}^{t+1},\forall u\in\mathcal{N}(v)\big{\}}) \tag{2}\]

\[h_{v}^{t+1}=\psi(h_{v}^{t},m_{v}^{t+1}) \tag{3}\]
where \(\phi\), \(\rho\), and \(\psi\) are (learnable) message, reduce, and update functions, respectively, and \(e_{u\to v}^{t}\) is a vector representation of edge \(u\to v\) at step \(t\); the calculation of node features can optionally incorporate edge features.

Our approach involves applying the message-passing method to learn representations of materials, which consist of molecular building units rather than distinct atoms. The proposed framework heavily relies on subgraph neural networks[53, 54] and includes steps inherent to their construction. First, the strategy for sampling/selecting subgraphs needs to be explicitly stated. Following chemical intuition, we decompose MOF structures (Figure 1a) into their inorganic ("metallic") and organic components (Figure 1b); a similar approach known as the "standard simplification" algorithm[55] has been previously employed to differentiate MOF structures. Trivial subgraphs of metal atoms (\(k\in\mathcal{K}\)) (Figure S1) are comprised of single nodes in the crystal graph; intermetallic bonds are ignored. Another set of subgraphs (\(l\in\mathcal{L}\)) is formed by the connected components of nonmetal atoms. Domain heuristics enable the development of alternative partition schemes that deconstruct MOFs into organic linkers and inorganic substructures, commonly known as second-building units (SBUs). Our preliminary experiments show that those methods are not universally applicable and cannot extract molecular building units for a considerable fraction of synthesized MOFs. Thus, 94.5% and 71.8% of entities from the Computation-Ready, Experimental Metal-Organic Framework (CoRE MOF) database[56, 57] are correctly processed by the MOFid[58] and moffragmentor[59] packages, respectively. On the other hand, the implemented scheme achieves identification of metal centers and organic linkers in 99.0% of structures. The decomposition rates of other databases are similar (Figure S2). In the case of COFs, we apply a similar methodology, with the only variation being the composition of inorganic nodes in terms of chemical elements; here boron and silicon are also taken into account (Figure S1). Consequently, the current stage of development does not support the adequate processing of a diverse family of metal-free COF structures within the proposed framework.

The next step is to define subgraph representations. Crystal graph neural networks appear to be the most relevant reference in this context; during training, learnable atomic embeddings are commonly generated in the input layers of a model. In opposition, we rely on predefined vectors to characterize both the inorganic \(\mathcal{K}\) and organic \(\mathcal{L}\) subsets. Specifically, the mol2vec model[60] is implemented to represent molecular fragments, and mat2vec embeddings[61] serve as features for single-node metal subgraphs (Figure 1c).

The considered sets of entities (\(\mathcal{K}\) and \(\mathcal{L}\)) form a bipartite graph (\(\mathcal{K},\mathcal{L},\mathcal{E}\)) by design; each edge in the graph connects species of different types, i.e., organic linker to metal center and _vice versa_. Unfortunately, reticular design does not provide prior knowledge of the specific bonds that will form in an experimental crystal structure. In order to eliminate arbitrariness in selecting initial connectivity, we postulate every potential edge, forming a complete bipartite graph (\(\mathcal{K},\mathcal{L},\tilde{\mathcal{E}}\)) called the coarse-grained crystal graph. The term refers to the coarse-graining modeling approach, which involves the simplification of complex atomic systems, e.g., reducing the number of degrees of freedom by representing groups of atoms as a single pseudo-atom. Similar techniques have been applied to quantify MOF diversity[62, 63] (in an unsupervised manner) and establish structure-property relationships in a specific class of hybrid materials[64] (zeolitic imidazolate frameworks). To the best of our knowledge, this study is the first to use the message-passing paradigm to model reticular materials in the coarse-grained regime.

Coarse-grained crystal graphs can be easily integrated into graph neural networks designed for heterogeneous graphs. To handle structure-property relationships in reticular materials, we implement a model with simple architecture (Figure 1d), hereafter designated as Coarse-Grained Crystal Graph Neural Network (CG\({}^{2}\)-NN). Predefined embeddings initiate the sequential update of node representations through the message-passing procedure, occurring in three interaction blocks; each block includes a convolutional layer, layer normalization[65], and nonlinearity, Exponential Linear Unit[66] (ELU). Owing to the universal nature of the proposed structure representation, a wide range of convolutional operations can be integrated into the interaction block. In this study, we perform experiments using four message-passing techniques: Graph Convolutional Network layer[67] (CG\({}^{2}\)-GCN), SAmple and aggreGatE layer[68] (CG\({}^{2}\)-SAGE), Graph ATtention layer[69] (CG\({}^{2}\)-GAT), and Graph ATtention layer with a universal approximator attention function[70] (CG\({}^{2}\)-GATv2). It should be also emphasized that the aggregation phase of message passing is obviously influenced by the specific topology of coarse-grained crystal graph (Equation 2). Since each linker has a neighborhood comprised of all metal centers and each metal center has a neighborhood comprised of all linkers, this phase looks as follows:

\[m_{k}^{t+1}=\rho(\{m_{l\to k}^{t+1},\forall l\in\mathcal{L}\}) \tag{4}\]

\[m_{l}^{t+1}=\rho(\{m_{k\to l}^{t+1},\forall k\in\mathcal{K}\}) \tag{5}\]
The calculation of graph-level representation is based on the following readout function:

\[h_{g}=\frac{1}{|\mathcal{K}|}\sum_{k\in\mathcal{K}}h_{k}+\frac{1}{|\mathcal{L}|} \sum_{l\in\mathcal{L}}h_{l} \tag{6}\]

Finally, the output of CG\({}^{2}\)-NN is obtained by passing the graph-level embedding \(h_{g}\) through a dense layer.

Figure 1: Coarse-grained crystal graph framework. Overview of data processing and scalability. (a) Original crystal structures of notable reticular materials, metal–organic framework HKUST-1 and covalent organic framework COF-5. Atoms are colored based on their chemical element. (b) The crystal graphs of HKUST-1 and COF-5 are colored based on their substructure type. The metal and nonmetal atoms are colored by violet and yellow, respectively. (c) Scheme of the coarse-grained crystal graph construction. (d) Schematic diagram of the coarse-grained crystal graph neural network architecture implemented in this study. (e) Number of nodes in the coarse-grained crystal graph as a function of number of atoms in the corresponding unit cell. The top panel contains the corresponding distribution of number of atoms.
In terms of message-passing scalability, the coarse-grained crystal graph framework is expected to outperform the full crystal graph. The Quantum MOF (QMOF) database[71, 72] contains structures with hundreds of atoms in the primitive cell (Figure 1e). It should be noted that the initial set of candidate materials was formed considering the limitations of high-throughput density functional theory (DFT) calculations. The structures in other experimental MOF subsets, such as the Crystal Structure Database (CSD) MOF Collection[73] and CoRE MOF database, can contain as many as ten thousand atoms (Figure S3). In contrast, the coarse-grained crystal graph have a maximum of 9, 17, or 35 nodes in the three databases mentioned above; five or fewer vertices are common in most coarse-grained crystal graphs (98.1%, 91.1%, and 92.8%). Probably the most illustrative quantity--the average ratio of nodes in the crystal graph to nodes in the coarse-grained crystal graph--equals 35.3, 125.6, and 101.8 in the QMOF database, CSD MOF Collection, and CoRE MOF database, respectively. However, the impressive scalability is only partially attributed to coarse-graining molecular subgraphs. Another factor is the occurrence of multiple identical subgraphs in the primitive cell; the coarse-grained crystal graph is intentionally free of duplicates, which means that a specific ratio of building units in the original reticular structure is disregarded. As a result, CG\({}^{2}\)-NNs are likely to have lower predictive performance compared to crystal graph neural networks. The following sections primarily address the efficiency-vs-accuracy dilemma of the proposed computational framework.

### Predictive performance of CG\({}^{2}\)-NNs

To evaluate the predictive performance of CG\({}^{2}\)-NNs, we consider a diverse set of practically important MOF properties, including band gap, thermal decomposition temperature, heat capacity, and Henry coefficients of eight gases (N\({}_{2}\), O\({}_{2}\), Kr, Xe, CH\({}_{4}\), CO\({}_{2}\), H\({}_{2}\)O, H\({}_{2}\)S); the examined tasks are all in a regression setting. In addition to the proposed reticular-specific architecture, several other neural networks that are composition-based and (crystal-)structure-aware are benchmarked as well. The algorithm designated as Representation Learning from Stoichiometry[31] (Roost) and model under the framework of Wyckoff Representation regression[74] (Wren) form the first group. It is interesting to note that Wren incorporates Wyckoff representations in addition to stoichiometry, introducing the concept of coordinate-free coarse graining into materials discovery. We categorize this approach as composition-based because the model does not directly incorporate the full crystal graph. Another group of methods includes four crystal graph neural networks: Crystal Graph Convolutional Neural Network[20] (CGCNN), MatErials Graph Network[21] (MEGNet), Global ATtention-based Graph Neural Network with differentiable group normalization and residual connection[28] (DeeperGATGNN), and Atomistic Line Graph Neural Network[27] (ALIGNN). In addition to "general purpose" structure-aware architectures, a few recent studies have introduced neural networks that integrate domain knowledge related to reticular chemistry; MOFNet[46], MOFTransformer[47], and MOFormer[48] are worth mentioning. MOF-related models are not part of our benchmark study; the original studies provide an overall picture of accuracy by comparing with crystal graph neural networks, e.g., CGCNN.

Despite acknowledging the criticism of using the coefficient of determination (\(R^{2}\)) as a primary predictive performance measure[75], we still find it useful for identifying general trends (Figure 2a). To ensure completeness, we have provided the corresponding mean absolute error (MAE) and root mean squared error (RMSE) values for all tasks (Figures S4 and S5). Based on the metric values, we can conclude that CG\({}^{2}\)-NNs are surprisingly effective in predicting DFT band gap. A coefficient of determination for the best models (CG\({}^{2}\)-SAGE) is 0.86, 0.85, and 0.76 at three different levels of theory (see details in the Methods section). In comparison, the best crystal graph neural networks (ALIGNN) have a coefficient of determination of 0.90, 0.92, and 0.84 for the same tasks. The accuracy of composition-based models is substantially lower, with a coefficient of determination of 0.72, 0.75, and 0.54, respectively. The prediction of thermal decomposition temperature presents a challenge for all implemented models. In particular, CG\({}^{2}\)-SAGE and ALIGNN reach a coefficient of determination of 0.38 and 0.44, respectively. We can speculate that a low accuracy is mainly associated with a relatively high uncertainty in determining the target value from the thermogravimetric analysis (TGA) data. A typical rounding step of 25 \({}^{\circ}\)C (applied to decomposition temperature values[76]) is comparable to MAE of predictive models: 50 \({}^{\circ}\)C for CG\({}^{2}\)-SAGE and 47 \({}^{\circ}\)C for ALIGNN. Therefore, the semi-quantitative agreement between experimental and predicted values can hardly be enhanced without expanding a set of target values and improving the resolution limit of TGA data. Nevertheless, CG\({}^{2}\)-NNs and structure-aware models have similar coefficients of determination, 0.32-0.38 vs. 0.30-0.44. The next endpoint, thermal capacity, poses a different challenge for predictive models. Unexpected model rankings and ineffective training are caused by a shortage of data points in the dataset (214). Composition-based Roost and structure-aware ALIGNN have comparable performance, with an average coefficient of determination of 0.68 for four temperature values. On the other hand, CGCNN and MEGNet, both incorporating crystal graph data, demonstrate nearly zero predictive performance. CG\({}^{2}\)-NNs exhibit moderate accuracy; CG\({}^{2}\)-GCN achieves the highest average coefficient of determination of 0.52 among these models. Finally, Henry coefficients are considered in the
analysis (Figure 2a); the following metric values were calculated by averaging eight endpoint-wise quantities. ALIGNN with a coefficient of determination of 0.67 significantly outperforms all the other models; the second best MEGNet has a coefficient of determination of 0.54. CG\({}^{2}\)-SAGE achieves a coefficient of determination of 0.48, surpassing structure-aware CGCNN (0.46). Both composition-based models, Roost and Wren, show limited predictive ability with a coefficient of determination of 0.29 and 0.30, respectively. The difference in predictive performance between ALIGNN and other models may be due to the significance of the specific geometry of adsorption sites in MOF; the explicit inclusion of three-body interaction terms in ALIGNN enables precise handling of atomic environments. To sum up, CG\({}^{2}\)-NNs demonstrate accuracy comparable to simplistic structure-aware models, e.g., CGCNN, and generally outperform composition-based algorithms.

Figure 2: Predictive performance of three model classes. (a) Coefficients of determination (\(R^{2}\)) of composition-based, coarse-grained crystal graph, and crystal-structure-aware neural networks are colored by green, blue, and violet, respectively. (b–d) Mean absolute error (MAE) in the band gap prediction as a function of training dataset size.
Further, we analyze the scalability of the proposed computational framework with respect to the size of the training dataset (Figure 2b-d), focusing on the PBE band gap; the models from each of the mentioned above class (composition-based, coarse-grained, and structure-aware) are examined. We approximate how MAE depends on training-dataset size \(N\) using a linear fit in a logarithmic scale for three chosen models: Wren, CG\({}^{2}\)-SAGE, and ALIGNN. As it follows from the fitted coefficients \(a\) and \(b\) in equation \(\lg MAE=a+b\lg N\) (Figure 2b-d), CG\({}^{2}\)-SAGE and ALIGNN show similar behavior, whereas the increase in training-dataset size leads to a much smaller boost in performance for Wren. At the same time, most structure-aware models are outperformed by Wren and another composition-based model (Roost) in the small-data regime (200 points in the training dataset); the observed tendency for heat capacity prediction is reproduced here.

The lack of data for the endpoint of interest can be addressed by leveraging advanced techniques, including transfer learning[77, 78, 79] and self-supervised learning[80, 81]. To explore the potential of applying CG\({}^{2}\)-NNs in conjunction with one of such algorithms, we consider the task of band gap prediction across different classes of reticular materials: the models pretrained on the MOF band gap are used to evaluate the same property in COFs (Figure 3a). The models trained from scratch, i.e., without pretraining, serve as a baseline. CG\({}^{2}\)-GAT and ALIGNN show the lowest MAE of 0.31 eV; most of the models perform at a similar level. In contrast, two structure-aware models--CGCNN and MEGNet--fail to quantitatively reproduce the target value, as observed in the thermal capacity prediction of MOFs (Figure 2a). From the other hand, the transfer-learning technique results in the largest increase in accuracy (i.e., decrease in MAE) for these models: 0.16 eV and 0.34 eV, respectively. The lowest MAE of 0.26 eV is achieved by the fine-tuned ALIGNN model; DeeperGATGNN with MAE of 0.29 eV is the second-best model owing to the moderate performance improvement (0.05 eV) associated with fine-tuning. Pretraining with MOF data has minimal effect on the accuracy of all CG\({}^{2}\)-NNs. Minor benefits from transfer learning may stem from incorporating pretrained embeddings into the proposed neural network architecture; at the concept level, this approach is essentially equivalent to freezing the weights of the input layer. The mol2vec (mat2vec) embeddings have been learned by exploring the vast chemical

Figure 3: Cross-domain (from metal–organic frameworks to covalent organic frameworks) transferability of coarse-grained crystal graph neural networks. (a) Mean absolute error (MAE) in the band gap prediction. MAE values for models pretrained on MOF data are shown in color, while gray represents models exclusively trained on COF data. (b) Scatter plot of calculated and predicted values of band gap. (c) Scatter plot of calculated and predicted values of band gap; linear scaling to minimize MAE is applied. (d) The two-dimensional projection of linker chemical space is produced within the Uniform Manifold Approximation and Projection (UMAP) algorithm from mol2vec embeddings; the linkers from the Quantum MOF database and from the subset of CURATED COFs are colored by gray and blue, respectively.
space, which extends beyond the organic linkers (metal centers) in reticular materials. In other words, the generalizability of low-level representations of molecular building units allows for the description of both MOFs and COFs from scratch. The band gap of COFs can be semi-quantitatively reproduced by CG\({}^{2}\)-GAT trained on MOF data, without additional fine-tuning; only a linear scaling of the model outputs is needed (Figure 3b,c). From the other perspective, a common chemical space is formed by the two-dimensional projections of organic linkers in MOFs and COFs, generated within the Uniform Manifold Approximation and Projection[82] (UMAP) approach from mol2vec embeddings (Figure 3d).

Figure 4: Interplay between predictive performance and energy efficiency of composition-based, coarse-grained crystal graph, and crystal-structure-aware models. The following endpoints are considered: (a) PBE band gap, (b) HLE17 band gap, (c) HSE06 band gap, (d) thermal decomposition temperature, (e) heat capacity at 250 K, (f) heat capacity at 300 K, (g) heat capacity at 350 K, (h) heat capacity at 400 K, (i) N\({}_{2}\) Henry coefficient, (j) O\({}_{2}\) Henry coefficient, (k) Kr Henry coefficient, (l) Xe Henry coefficient, (m) CH\({}_{4}\) Henry coefficient, (m) CO\({}_{2}\) Henry coefficient, (o) H\({}_{2}\)O Henry coefficient, and (p) H\({}_{2}\)S Henry coefficient. For the illustrative purposes, the areas of the MAE vs. GHG emissions space accessible by different models are painted in two colors: yellow-colored regions correspond to the set of composition-based and structure-aware models, whereas red-colored fields reflect how the Pareto front is reshaped by introducing CG\({}^{2}\)-NNs.
### Carbon footprint of training CG\({}^{2}\)-NNs

The main focus of the materials informatics community has been on enhancing the predictive performance of developed neural network architectures[83], paying little attention to other critical aspects such as computational efficiency, explainability, transferability, and scalability. As it was shown in our recent study[84], an excessive focus on model accuracy has led to an exponential growth in trainable parameters and greenhouse gas (GHG) emissions from model training. Taking into account a decent predictive ability and high scalability (in terms of input data size) of CG\({}^{2}\)-NNs, we also estimate the carbon footprint of our models in the hope that the presented framework can tackle the accuracy-efficiency dilemma in materials property prediction. In Figure 4, the interplay between GHG emissions (measured in kilograms of carbon dioxide equivalents, kgCO\({}_{2}\)eq) and predictive performance (expressed in terms of MAE) is analyzed for MOF-related models from the previous section. It is worth noting that the carbon footprint of the model lifecycle is solely determined by the electrical energy consumption of the hardware in use. Accordingly, GHG emissions can provide a description of energy efficiency, in addition to quantifying environmental impacts. In contrast to the issue of ranking models by the accuracy objective, here we deal with a set of so-called non-dominated solutions, i.e., models that provide a tradeoff between target quantities. Thus, ALIGNN is positioned on the low-error section of the Pareto front in all cases, except for thermal-capacity prediction. The composition-based Roost model is the only exception, as it shows high predictive performance in the small-data regime. Roost and CG\({}^{2}\)-SAGE dominate the low-emission region of the Pareto front, confirming the impressive computational efficiency of the presented framework. As we can see, CG\({}^{2}\)-NNs offer nearly state-of-the-art performance at a fraction of the computational cost of previous algorithms. For instance, the PBE band gap MAE of CG\({}^{2}\)-SAGE can be reduced by 20% (22%) with MEGNet (ALIGNN), causing GHG emissions to increase 13 (89) times; the similar estimates apply to other endpoints as well. Therefore, the coarse-grained crystal graph concept offers a strong alternative to current dominant methods for representing reticular materials, if predictive performance and energy efficiency are equally important.

### Inverse materials design accompanied with CG\({}^{2}\)-NNs

CG\({}^{2}\)-NNs are readily accessible as discriminative algorithms for _in silico_ high-throughput screening of reticular materials. However, the presented computational framework can be integrated into inverse design pipelines as well. As a proof-of-concept study, we explore maximizing hydrogen storage capacity in MOFs; the optimization problem is given as follows:

\[x^{*}=\operatorname*{argmax}_{x\in\mathcal{X}}f(x) \tag{7}\]

where \(\mathcal{X}\) is a design space represented in the coarse-grained crystal graph framework, \(f(x)\) is an objective function (hydrogen storage capacity) approximated by CG\({}^{2}\)-NN, and \(x^{*}\) is an optimal set of metal centers \(\mathcal{K}\) and organic linkers \(\mathcal{L}\) encoded as a complete bipartite graph (\(\mathcal{K},\mathcal{L},\tilde{\mathcal{E}}\)). Using IRMOF-20[85] as a prototype, we apply an iterative evolutionary procedure to modify the linker (thieno[3,2-b]thiophene-2,5-dicarboxylic acid) and maximize the target property; the metal center (zinc) and another linker (oxygen atom) are maintained unchanged. In other words, we move from a global optimization problem to a local maximization one; our intention is to confine the design space to the region surrounding the original structure. At each optimization step (Figure 5a), SELF-referencIng Embedded Strings[86] (SELFIES) representation of molecular building block undergoes one of three operations: a single character addition, deletion or replacement. SELFIES strings are 100% chemically robust by design, satisfying valence-bond rules under any mutation, but other specific requirements should also be fulfilled by a potential linker. In particular, the "linker-likeness" and synthetic-accessibility filters are implemented into the pipeline (see details in the Methods section). The linkers that meet the above criteria are passed through CG\({}^{2}\)-SAGE as a part of the coarse-grained crystal graph for predicting hydrogen storage capacity. Further, we evaluate the uncertainty in predictions using the deep-ensemble method[87]. The variance of model outputs calculated for a set of neural networks may serve as a reliable measure of epistemic uncertainty, which indicates the model's limitations in reproducing structure-property relationships outside its domain of applicability. As we see in Figure 5b, narrower prediction intervals in terms of percentiles (defined by means of quantile regression[88]) are associated with a lower standard deviation of model ensemble outputs. Accordingly, in-domain MAE decreases with lowering the upper bound of standard deviation (Figure 5c); the structure is considered inside the domain of applicability if the uncertainty measure is below the predefined threshold. In the context of optimization problems, the ranking capability of model appears to be a key indicator; the implied threshold value (2.5 g/l) enables reaching a Pearson correlation coefficient of 0.86 (Figure 5c). Finally, the hydrogen storage capacity of the generated linker must exceed the lowest one in the previous population.
Additionally, the top-performing linkers undergo filtering based on their 3D structure. Rigid molecules with antiparallel binding sites, i.e., carboxylic groups, are considered as potential precursors. The set consists of molecules that possess both the original framework (linker **4**) and a new one (linkers **1**, **2** and **3**); they are all located in the vicinity of the initial compound (thieno[3,2-b]thiophene-2,5-dicarboxylic acid) in the chemical space (Figure 5d). The parent MOF structure (IRMOF-20) has a nearly record-breaking performance (experimental volumetric capacity of 51 g/[89]), so the generated linkers only slightly enhance the predicted target property value (up to 3.8 g/l). Expanding the design space of the optimization problem, such as by loosening the uncertainty-quantification criterion, may potentially improve the quality of found solutions; simultaneous optimization of metal center and organic linker achieves the same goal. Nonetheless, the molecules that have been identified show potential as building units for MOFs with a high hydrogen working capacity. The presented inverse-design pipeline provides outputs that can be readily utilized to predict MOF synthesis parameters[90, 91].

## 3 Discussion

As it was shown in the Results section, CG\({}^{2}\)-NNs are able to compete in the predictive performance with models incorporating atomic connectivity information. This suggests that the underlying topology of reticular material can be learned during model training. We intentionally consider datasets containing experimentally resolved MOFs and COFs; the relationship between a set of molecular building blocks and their self-assembled

Figure 5: Inverse reticular design driven by coarse-grained crystal graph neural networks. (a) Schematic diagram of the linker optimization workflow. SELFIES is an abbreviation of SELF-referencIng Embedded Strings; UQ is an abbreviation of uncertainty quantification. (b) Predictive error as a function of model ensemble standard deviation. The colored areas correspond to the predictive intervals estimated by quantile regression: within one standard deviation (green), within two standard deviation (orange), and over two standard deviations (red). (c) In-domain mean absolute error (MAE) and Pearson correlation coefficient as a function of model ensemble standard deviation. The threshold value of standard deviation (2.5 g/l) is indicated by vertical line. (d) The initial linker (taken from IRMOF-20) and four generated molecules with the highest hydrogen storage capacity. The two-dimensional projection of linker chemical space is produced within the Uniform Manifold Approximation and Projection (UMAP) algorithm from mol2vec embeddings; the linkers from MOFs used for training objective function predictors are shown in gray.
structure is usually straightforward. On the contrary, widely used databases of hypothetical structures expand the reticular materials genome by providing multiple structures of various topologies for specific organic linker and SBU. Obviously, CG\({}^{2}\)-NNs will fail to distinguish such structures, although most of _in silico_ generated MOFs thermodynamically inaccessible and, as a result, are not of practical interest. The synthesized MOFs, which exhibit polymorphism[92, 93], will cause the same issue. Incorporating global state attributes, e.g., temperature and pressure, into CG\({}^{2}\)-NN architecture can potentially overcome the challenge of polymorphic phase transitions.

By eliminating topology from explicit consideration, we address the key shortcoming of current approaches to generating reticular materials. The entire design space is defined by the linkers-nodes-topology triad, but the optimization task is constrained by reticular chemistry principles: topology is determined by molecular building units and is not an independent variable. The issue has been mostly ignored in previous studies[44, 94, 95], whereas the energy ranking of polymorphs demands considerable computing resources[50, 96]. We anticipate that the coarse-grained crystal graph framework will shift community focus from navigating inaccessible regions of the reticular materials genome to improving methods for predicting synthesis conditions.

## 4 Conclusion

The introduction of the coarse-grained crystal graph framework aims to make data-driven design of reticular materials more accessible for synthetic chemists. The efficiency of the presented approach is achieved by integrating relevant domain knowledge into the neural network architecture. Specifically, pretrained embeddings are used to represent organic linkers and metal centers, demonstrating the impact of previous community efforts on AI tools for materials design. We hope that this research will not be the final link in the chain, but rather bolster the further enhancement of predictive models.

## 5 Methods

### Datasets

The QMOF database[71, 72] served as a source for MOF crystal structures and their corresponding band gap values obtained at three levels of theory: generalized gradient approximation (GGA), meta-GGA, and screened-exchange hybrid GGA; the density functionals in use were Perdew-Burke-Ernzerhof[97] (PBE), High Local Exchange 2017[98] (HLE17), and Heyd-Scuseria-Ernzerhof[99] (HSE06), respectively. The datasets consisted of 20237 (PBE endpoint), 10664 (HLE17), and 10718 (HSE06) entities. Another collection of neural networks was trained using a subset of structures (3038 points) from the CoRE MOF 2019 database and their corresponding decomposition temperatures extracted from the TGA data by Nandy et al.[100] For the heat capacity prediction, we utilized the values obtained within the harmonic approximation at the GGA level of theory[101]; 214 MOFs from the CoRE MOF 2019 and QMOF database were taken into account. The Henry coefficients of eight gears--N\({}_{2}\), O\({}_{2}\), Kr, Xe, CH\({}_{4}\), CO\({}_{2}\), H\({}_{2}\)O, H\({}_{2}\)S--computed using the grand canonical Monte Carlo (GCMC) simulations and the corresponding structures from the QMOF database (1431, 1552, 1297, 1205, 1268, 1538, 1482, and 1352 compounds, respectively) were taken from the dataset presented by Jablonka et al.[59] The PBE band gap values[102] of 61 compounds from the CURATED COFs database[103] containing boron or silicon atoms were used to estimate the transferability of neural networks of interest. Usable volumetric hydrogen capacity of MOFs under temperature-pressure swing conditions (77 K/100 bar and 160 K/5 bar) was implemented as an objective in the optimization task; 4146 structures from the CoRE MOF 2019 database and GCMC values were taken from the study[104] to train an ensemble of CG\({}^{2}\)-SAGE models.

### Coarse-grained crystal graph processing

The procedure outlined below was applied to construct the coarse-grained crystal graph and assign features to its nodes. Metals (in MOFs and COFs) and metalloids (in COFs) were completely removed from the initial crystal structure (the full list of atoms types is provided in Figure S1). The removed atoms are considered as metal centers, each of them was featurized using 200-dimensional mat2vec embeddings[61]. The crystal structure containing only nonmetal atoms was processed within the OpenBabel routines[105] to produce Simplified Molecular Input Line Entry System[106, 107] (SMILES) strings from connected components in the reduced crystal graph. SMILES strings were then converted into 300-dimensional mol2vec embeddings[60]. Metallic and organic units were built into the complete bipartite graph, where each node of one type is
connected to a node of another type. In this form, the coarse-grained crystal graph was passed through a graph neural network.

### Model training

CG\({}^{2}\)-NNs were built with PyTorch[108] and Deep Graph Library[109] (DGL). The choice of deep learning framework has a substantial effect on energy efficiency of model training[110]; therefore, PyTorch-based implementations were utilized for all other models as well. In the predictive performance (transfer learning) study, we used the Adam optimizer[111] with a learning rate of 10\({}^{-3}\) (10\({}^{-4}\)) and a batch size of 64 (16); the maximum number of epochs and early stopping criterion were set to 500 and 50, respectively. The five-fold cross-validation was performed for model evaluation; one-eighth of training data were used for early stopping. To evaluate an objective in the optimization task (hydrogen storage capacity), the ensemble of 20 CG\({}^{2}\)-SAGE models was trained by means of the Adam optimizer with a learning rate of 10\({}^{-3}\) and a batch size of 64. The holdout cross-validation technique with a training-validation-test ratio of 80:10:10 was applied. The mean and standard deviation of model ensemble outputs were used to estimate the target property and the corresponding uncertainty[87]. The Eco2AI library[112] was utilized to estimate GHG emissions of model training; the emission intensity coefficient was set 240.56 kgCO\({}_{2}\)e/MWh (Moscow). All experiments in the study were performed on the workstation equipped with two Intel(r) Xeon(r) CPUs E5-2695 v4 @ 2.10GHz, 144 GB RAM, and NVIDIA GeForce RTX 3090 Ti.

### Inverse reticular design

The linker optimization algorithm was heavily inspired by the previously developed methods[113, 114, 115]. Three symbolic operations (addition, deletion and replacement) were applied to modify SELFIES strings[86]; each optimization step included from one to twenty randomly selected mutations and several filtering procedures described below. Initial population was produced by generating one hundred thousand mutants from the SELFIES string, corresponding to the parent structure (thieno[3,2-b]thiophene-2,5-dicarboxylic acid); one hundred strings with the highest hydrogen storage capacity value were taken for further processing. Then, one thousand mutations steps were carried out. At each iteration, the population was updated by including a newly generated molecule if its objective was higher than the lowest value in the current generation; the number of molecules in the population was maintained at one hundred. Three filters were implemented to exclude irrelevant molecules. First, we assumed exactly two binding sites (carboxylic groups) in molecule by applying the "linker-likeness" filter. Second, the molecules with a low synthetic accessibility (SAscore[116] higher than 5.0) were also ignored. Third, 20 CG\({}^{2}\)-SAGE models were used to access uncertainty in predictions of target property by means of deep-ensemble learning[87]. If the standard deviation of model ensemble outputs exceeded 2.5 g/l, the corresponding molecule was discarded as well. One hundred launches of the optimization procedure were done in total.

The SELFIES strings included in the last generations were transferred into 3D atomic coordinates. Specifically, the Experimental-Torsion basic Knowledge Distance Geometry[117] (ETKDG) approach implemented in the RDKit library was employed to generate 30 conformers for each of the 30 molecules with the highest hydrogen storage capacity. Then, the local geometry optimization was performed using the ANI-2x neural-network potential[118] and Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm with the convergence criterion of 10-3 eV/A for forces. Two quantities were calculated for each ensemble of conformers: 1) the mean angle between two vectors, orienting from the carbon atom to the midpoint between oxygen atoms in carboxyl groups, 2) standard deviation of pairwise distance between carbon atoms in carboxyl groups. A molecule was included into the final set if the former value was higher than 150deg and the latter one was lower than 0.5 A; the linker optimization was completed on four structures.

## 6 Acknowledgements

V.K. was supported by the Fellowship from Non-commercial Foundation for the Advancement of Science and Education INTELLECT.

## 7 Data and code availability

The source code and data accompanying this study will be published soon.