Hybrid localized graph kernel for machine learning energy-related properties of molecules and solids

Nowadays, the coupling of electronic structure and machine learning techniques serves as a powerful tool to predict chemical and physical properties of a broad range of systems. With the aim of improving the accuracy of predictions, a large number of representations for molecules and solids for machine learning applications has been developed. In this work we propose a novel descriptor based on the notion of molecular graph. While graphs are largely employed in classification problems in cheminformatics or bioinformatics, they are not often used in regression problem, especially of energy-related properties. Our method is based on a local decomposition of atomic environments and on the hybridization of two kernel functions: a graph kernel contribution that describes the chemical pattern and a Coulomb label contribution that 1encodes finer details of the local geometry. The accuracy of this new kernel method in energy predictions of molecular and condensed phase systems is demonstrated by considering the popular QM7 and BA10 datasets. These examples show that the hybrid localized graph kernel outperforms traditional approaches such as, for example, the smooth overlap of atomic positions (SOAP) and the Coulomb matrices.


Introduction
The past decade has seen an impressive growth in the development and application of machine learning techniques 1 in quantum chemistry and computational condensed matter physics [2][3][4] . These methods are of a great interest for theoreticians because they allow for the analysis, classification and prediction of various properties conventionally requiring a large amount of data generated via computationally demanding quantum mechanical calculations 5 . Indeed, machine learning techniques can be applied to a broad range of problems, including, among the others, potential energy surface fitting [6][7][8][9] , ab initio molecular dynamics [10][11][12] , prediction of various scalar properties [13][14][15] (e.g atomization energies, polarizability coefficients, highest occupied molecular orbital energies, electronic structure correlation energies, etc), and vectorial and tensorial quantities (e.g forces, polarizability tensors, etc) 16,17 .
The chemical compound space is characteristic by a huge dimensionality and complexity.
Datasets for molecules and solids proposed in the literature, which span only a small part of the chemical space, already contain impressively large numbers of compounds. The GDB-17 molecular dataset, for instance, contains 166 billions of molecules 18 . An example of a rich dataset of crystal structures is the set of elpasolites of Faber et al. 19 containing 2 millions crystals. This shows clearly that it is impossible to analyze or screen all these compositions and structures by demanding electronic structure calculations 2 . A solution to efficiently explore the chemical compound space (or at least significant parts of it) can be found in the development of new machine learning approaches to complement ab initio calculations. This field of research is known as quantum machine learning 20,21 (QML). The main idea of the QML approaches is to train a machine learning model on a subset of chemical structures (the smallest possible) for which simulations at a quantum mechanical level were done. The trained machine is then used to predict the target properties of the rest of the systems. This technique can be applied not only in the chemical composition space but also in the conformation space of a given system. For example, recently, a ∆-machine learning approach 22 based on thermodynamic perturbation theory has shown impressive results in the prediction of the adsorption enthalpy of small molecules adsorbed in chabazite at the random phase approximation (RPA) level of theory using only 10 training conformations 23 .
In general, the ultimate goal of QML would be to find a universal model which unifies both chemical and conformational spaces 21 .
One of the main challenges in the use of QML methods is the choice of proper descriptors -i.e. features -to represent the molecular and/or condensed matter systems. In this work we will focus on structural descriptors, namely approaches that assume the knowledge of the precise geometry of the systems under consideration beyond their chemical composition.
Currently, a large number of those descriptors have been proposed in the literature, with the symmetry functions 24 , the Coulomb matrices [25][26][27] , the smooth overlap of atomic positions (SOAP) 28,29 and the many-body tensor representation (MBTR) 30 being some of the most prominent examples. From a general point of view an efficient structural descriptor should satisfy certain prerequisites 3 . For instance, the descriptor should be unique (one-to-one correspondence between features and systems) and invariant with respect to overall rotations and rigid translations, as well as permutations of atoms of the same type.
In the present work we introduce a new descriptor based on the notion of graph that satisfies all the above-mentioned conditions. Graphs are already largely used in the domain of cheminformatics 31-34 and bioinformatics 35,36 to predict the activity or the toxicity of particular drugs or to recognize particular patterns of binding site of biomolecules (with the quantitative structure-activity relationships (QSAR) model 37 ). Nevertheless, in these cases, the use of graphs is restricted to classification problems and their applications in regres-sion problems are still rare, especially as far as energy-related properties are concerned [38][39][40][41] .
These quantities are particularly sensitive to inner changes of the molecular structure which cannot be described by simple graphs whose edges are represented by 0's and 1's for disconnected and connected nodes (namely atoms) 42 . Indeed, this simplified representation leads to isomorphic graphs even for rather different structures and this represents an issue for learning properties (e.g. the energy) which are very sensitive on the specific geometry of a system. To overcome this difficulty in the present work we used weighted graphs, whose values on the edges were obtained as superposition of atom-centered Gaussian functions.
In order to develop a machine approach applicable to both molecules and solids we considered graphs defined in localized environments around each atom of a given system. These localized graphs were used as input features of a machine learning approach based on the kernel ridge regression 43 (KRR). This method scales cubically with the size of the training set but typically requires fewer data point to be trained than other traditional approaches (e.g. neural networks). The role of the kernel in KRR is to provide a measure of similarity between different localized graphs belonging to different systems. A global "comparison" of two systems (molecules or solids) can then be obtained by summing localized contributions from each atomic environment 29 . The specific kernel used in this work combines two different parts: The first, based on a variant of the shortest-path kernel, is suited to describe structural motifs and chemical bonding between atoms; the second includes finer geometric details through a label enrichment of the nodes with Coulomb vectors. This approach can be viewed as an hybrid kernel function 44 and we will show that the introduction of both contributions provides a high level of accuracy. As a proof of principle this new methodology was applied to predict molecular atomization energies (QM7 dataset) and enthalpies of formations of solids (BA10 dataset).
The paper is organized as follows. In Sec. 2 the concept of labeled localized graph is introduced, which is represented trough maximum probability paths and enriched with a set of vectors that describe the Coulomb potential of the atomic environment. In the same section, the new kernel approach is also presented. In Sec. 3 the datasets QM7 and BA10, used to evaluate our method are reviewed. Finally, our results for the prediction of the atomization energies of molecular systems (QM7 dataset) and the formation enthalpies for the solid state systems (BA10 dataset) are presented in Sec. 4.

The Labeled Localized Graph Descriptor.
In this section we introduce the concept of labeled localized graph (LLG) for molecules and solids. By considering a local atomic environment X defined by a radius r cut , a LLG, de- -where N and M are associated with the number of atoms and the number of chemical bonds in the atomic environment, respectively 45,46 . The LLG is said "labeled " if a label l corresponding to the chemical symbol is attributed to each node Fig. 1).
A simple unlabeled graph is usually represented through an adjacency matrix A(G X ) ∈ R N,N , whose elements A ij are set to one if two different nodes i and j are linked by an edge or to zero otherwise 45,47 . However, the same graph would correspond to an infinite number of conformations through this definition: by changing the interatomic distances, as long as the edges are preserved (i.e. bonds are not broken), the corresponding graph will not change.
This issue is particularly problematic for the prediction of total energy related properties (e.g. atomization energy), which are strongly influenced by the precise atomic conformation. For this reason in this work we use weighted adjacency matrices 38 which include significantly more information on the structure of the atomic environment. Specifically, a Gaussian function is attributed to each atomic position r i in the atomic environment. Instead of considering only binary values for the off-diagonal elements of the adjacency matrix, weights w ij between each pair of vertices v i and v j are set to correspond the overlap of these Gaussian functions, such that: Hence the weights depend on covalent radii r cov i (coming from crystallographic data 48 ) of all elements and on the hyperparameter γ controlling the width of the Gaussian atomic distributions. Finally, a threshold is applied as follows: i.e., an edge between the two nodes v i and v j is defined when w ij is larger than epsilon (see Fig. 1). These weights can be interpreted intuitively as probabilities to "jump" from one atom to another atom.
Weigthed adjacency matrix : Labels of vertices : N, H 1 , H 2 , C Figure 1: Illustration of a labeled localized graph (LLG) for a nitrogen atomic environment defined by a radius r cut . The overlaps of the Gaussian functions associated to each atom in this environment determine the weights in the adjacency matrix of the graph.

Maximum Probability Paths and Label Enrichment.
From the definition of the weighted adjacency matrix A(G X ) in Eq. 2, we can determine a path π (v i , v j ) between each pair of vertices (v i , v j ) (we suppose here that each node is connected to at least one another node, as it is the case for all the systems studied in this work). This path is defined by a finite-length sequence of vertices 31,32,49 . In this paper we will consider a machine learning approach based on a kernel that measures the similarity between paths.
Specifically, we will focus on a variant of the shortest path (SP) kernel. While determining all paths of a graph is a NP-hard problem, computing the shortest path between pairs of vertices is a problem solvable in polynomial time (O(n 3 )) 49 .
Typically, the determination of the shortest path is based on the distance and on the triangular inequality (Floyd-Warshall algorithm) 50 . Rather then directly using an approach based on shortest paths, in the context of the present work we found more natural to introduce the analogous concept of maximum probability path (MPP). In particular, interpreting the elements of the weighted adjacency matrix as transition probabilities between two nodes, a path is identified with the greatest overall probability expressed as the product of the weights associated with all the edges that are involved in the path. Mathematically, the probability of a certain path can be written as: where the intermediate indexes in the product correspond to all the nodes visited along the path. In order to find the MPP between two vertices, we modified the original Floyd-Warshall algorithm 49 (see Alg. 1).
Through our variant of the Floyd-Warshall algorithm, the matrix A(G X ) is transformed into a new matrix P(G X ) that can be interpreted as a fully connected graph in which the

Algorithm 1 The Floyd-Warshall algorithm for MPPs
Require: Weighted adjacency matrix A(G X ) Ensure: MPPs matrix P(G X ) for k = 1 to N do end for end for end for matrix elements P ij correspond to maximized overall probabilities of transition from the node i to j. In our method, the MPPs are also labeled by a sequence of labels l(π(v i , v j )) = (l(v i ), . . . , l(v k ), . . . , l(v j )). As discussed below (see Sec. 2.3 A.III), this is necessary in order to include information on atomic species so that only the MPPs involving the same sequences of atoms are compared (see Fig. 2).
The concepts described up to this point (weighted graphs and MPPs) allow for a rather accurate description of geometric motifs and bonding between atoms. However, these approaches are not sufficient to capture all the geometric details in atomic environments. Specifically, modifications (e.g. rotations) which preserve interatomic bond distances in the local environment leave the matrices A(G X ) and P(G X ) unchanged. These finer structural details, which are not captured by the MPP approach, might have different relevance depending on the specific dataset but, nevertheless, contribute to the non-uniqueness of the connection between localized graph and property to predict. To overcome this difficulty, we introduced a label enrichment of the vertices 31,32 . Inspired by the Coulomb matrix descriptor introduced by Rupp et al. 25 , we defined for each node a Coulomb vector c (i) sorted according to the distances in increasing order whose components are defined as Weigthed adjacency matrix : Path matrix : The weighted adjacency matrix A(G X ) is transformed into a path matrix P(G X ). Each path between pairs of atoms is characterized by a sequence of labels l(π(v i , v j )) and a maximum probability path π(v i , v j ).
where Z i , Z j are the nuclear charges. The superscript i refers to the node which is labeled, while the subscript j is an index that indicates all the other atoms in the environment of i (see Fig. 3).
Because the atomic composition of environments can be different, to prevent any issue with dimensionality mismatch, the Coulomb vectors are padded with zeros to set their dimension to be equal and transferable. Finally, we obtain a new set of labels that takes into account the missing geometric information in the description of atomic environments. Figure 3: Illustration of the labels enrichment of nodes. In the atomic environment X to each node it is associated a sorted Coulomb vector c (i) . Example with the carbon atom :

The Hybrid Maximum Probability Path (HMPP) Kernel.
The main idea of kernel approaches consists in finding an accurate measure of the similarity between pairs of data points (atomic environments in our case). For example, in the SOAP kernel, the similarity is determined through an overlap calculation of the atomic density of two atomic environments 28,29 . When normalized, the kernel is close to one when two atomic environments are nearly identical and tends to zero if they are dissimilar.
In our approach, an atomic environment X is characterized by a set of MPPs be- . Being inspired by the iterative similarity for molecular graphs introduced by Rupp et al. 51 , where the use of a linear combination of two kernels applied to vertices and edges was proposed, the following hybrid kernel is introduced to compare the atomic environments (see also Fig. 4): This kernel is normalized via Tanimoto's normalization 45 : which is commonly used also in other graphs application [45][46][47]52 . This normalization can be seen as a Jaccard's distance -i.e. measuring how two sets intersect -between two discrete collections of paths and vertices. The kernel k(X , Y) is positive semidefinite 52 , as required for use in kernel-based machine learning algorithms.
To compare graph paths, we have chosen an approach similar to the labeled shortest path graph kernel 49 . In its original definition this kernel could take values different from 0 only when applied to pairs of (shortest) paths with the same initial and final labels 49 . In this work we found that this formulation does not reach a satisfactory level of accuracy. Indeed, even if the initial and final nodes of a pair of MPPs have the same "atomic" labels (namely they correspond to the same atomic species), the intermediate paths could be significantly different from a chemical point of view (i.e. involve different elements) and still be considered as similar according to this original definition of the kernel. For this reason in the present work we used an MPP kernel which is strictly 0 unless the two paths to be compared have exactly the same sequence of atomic labels l(π). This kernel can be expressed as follows: where the δ Label (l (π X ) , l (π Y )) is a delta kernel applied on the path labels such that: This kernel is equal to one if, and only if, the sequences of labels along π X and π Y are identical. The term k Path (π X , π Y ) is the Laplacian kernel 53 used to compare pairs of MPPs. Similarly, the Coulomb labels are compared through a sum of Laplacian kernel functions: where N and N denote the number of atoms in the atomic environments X and Y, respectively, and β 1 and β 2 are the hyperparameters that control the decay of the Laplacian kernel functions. They play a central role in the prediction quality and are estimated through a grid search using a validation dataset. The relative importance of the two kernels k MPP and k Coulomb in the description of a given system is controlled by the hyperparameter α in Eq. 5.

Benchmarks
To demonstrate the accuracy of our hybrid graph kernel approach we applied it to two regression problems of energy-related properties and compared its performance to that of the popular SOAP descriptor 28,29 and graph approximated energy (GRAPE) 38 . The machine learning models were based on in-house programs written in Python3 employing the   Using the five-fold partitioning originally proposed by Rupp et al. 25 and following the work of Ferré et al. 38 , we considered the first N molecules (with N = 100, 300, 500 and 1000 points) of the first fold as a training set. The second and third folds were combined to create a validation set for the grid search of the optimal hyperparameters, which, in the case of the HMPP kernel include the parameters r cut , γ, β 1 , β 2 and α. Finally, the two last folds were used as a test set to evaluate the accuracy of the model through the mean absolute error (MAE) and the root mean square error (RMSE) of predicted energies.

The BA10 dataset.
The BA10 dataset contains standard enthalpies of formation for a set of ten binary al-

Optimization of the hyperparameters
The determination of the hyperparameters was conducted on a hold-out validation set through a five-dimensional grid search for the parameters r cut , γ, β 1 , β 2 , and α. The cutoff radius r cut , that defines the size of each atomic environment, was evaluated in the interval between 1.6 and 2.4Å with an increment of 0.2Å for the QM7 dataset, and between 2.0 to 5.0Å with an increment of 1.0Å for the BA10 dataset. The values of 1.8Å and 3.0Å were identified as optimal for the QM7 and BA10 datasets, respectively. The scaling factor γ, that governs the width of the atomic Gaussian functions, was tested in the range from 0.1 to 0.9Å with a step of 0.2Å; the optimal values of 0.5Å and 0.5Å have been identified for the QM7 and BA10 datasets, respectively.
The two hyperparameters β 1 and β 2 that define the decay rate of the two exponential functions in the k MPP and k Coulomb kernels were evaluated on a decimal logarithmic grid defined on the interval between 10 −1 to 10 −4 . For both the datasets considered here, the optimal hyperparameters values were found to be 0.1 and 0.001 for β 1 and β 2 , respectively.
Compared to the other three parameters in the model, β 1 and β 2 seem to have a weaker system dependence and could be possibly transferred across different datasets. Although this observation should be confirmed by future investigations on several different datasets, the possibility of fixing β 1 and β 2 a priori could help to simplify the hyperparameter optimization procedure.
The α parameter, that describes the relative contribution of the kernels k MPP and k Coulomb , was tuned by considering values between 0 to 1 with a step of 0.1. For the QM7 dataset, a relatively small α (0.2) exhibited the best performance, while a significantly higher value (0.9) was identified for the dataset BA10. It is important to notice that the hybridization brings important improvements with respect to the sole use of the best performing between the k MPP and k Coulomb kernels (for the smallest training sets considered below the hybridization lowers the mean absolute error by 35 % and 17 % for QM7 and BA10, respectively).

Prediction of the atomization energies.
In this section we test the performance of the HMPP kernel in predicting the atomization energies of the molecules in the QM7 dataset. Table 1 and Figure 5 show the variation of the MAEs and RMSEs as a function of the training set size. For the sake of comparison we also show results for some other well established approaches, namely the SOAP kernel and GRAPE, the latter being a localized graph kernel method previously reported in the literature 38 and applied to regression problems. The results for the different methods have been generated using the same training sets and analogous procedures for the optimization of the hyperparameters. We observe that the HMPP kernel outperforms the other two methods regardless of the training set size. It is important to notice that in Ref. 37 the results reported for the GRAPE kernel slightly improved over the SOAP method (which is in contrast to our results from Table 1 and Figure 5). However, the authors of this work clearly stated that the tuning of the hyperparameters was limited and that, accordingly, it was not possible to conclude "that one method outperforms the other". The MAEs and RMSEs reported in Table 1 for GRAPE and SOAP are also sizeably smaller than those in Ref. 37, demonstrating the importance of a fine tuning of the model hyperparameters to fully establish the accuracy of a certain approach.  Table 1).
The GRAPE kernel uses a random walk kernel directly applied on the weighted adjacency matrix. Within this approach, each component of the adjacency matrix is defined through an overlap of atomic Gaussian functions of the same width. However, we can observe that this use of graphs does not reach a satisfactory level of accuracy, as shown by the high MAE and RMSE values in Table 1 and Figure 5. This is likely related to the lack of labels in the GRAPE approach. Tang and de Jong have proposed in 2019 a marginalized graph kernel to predict the atomization energies and have obtained impressive results with the introduction of labels 41 . However, the reported model was based on a global description of the molecular graphs and, differently from our present approach, the corresponding graph kernel cannot Table 1: Variation of mean absolute error (MAE) and the root mean square error (RMSE) (in kcal/mol) of atomization energies (QM7 dataset) obtained using the SOAP, GRAPE and HMPP kernels with the training set size (cf. Figure 5).    4.3 Prediction of the standard enthalpies of formation. Table 3: Variation of MAEs and RMSEs (in kcal/mol) with the training set size for enthalpies of formation of the BA10 dataset predicted using the SOAP and HMPP kernels (cf. Figure 6).  Table 3 and Figure 6 show the values of MAE and RMSE averaged over all the structures of the 10 binary alloys. The training sets, which are always excluded in the error evaluation, are obtained by randomly selecting 100, 300, 500, and 1000 configurations for each binary compound among the 1595 structures -i.e. without stratification of the data.
As for the dataset QM7 discussed in the previous section (see Sec Table 3).
is sizeably larger than that of SOAP (in this comparison it should also be kept into account that hyperparameters for the HMPP model are not reoptimized for each alloy). In Figure   7 we also report results from the previous work of Nyshadham et al. 60    1000 SOAP · · · GPR 0.10 · · · 60

Conclusions
In conclusion, we introduced in the framework of a local decomposition kernel a new similarity measurement based on molecular graphs. This kernel is composed by two parts: one that describes the local molecular pattern through a labeled graph, and a second one, that keeps into account some finer geometric information using Coulomb labels. These two kernels are hybridized through a hyperparameter α which controls their relative contributions.
The accuracy of this new kernel was tested on two datasets: The molecular QM7 dataset and the BA10 dataset, containing ten binary alloys. The first of these applications involves the prediction of the atomization energies of molecules. In this case our method outperforms previous descriptors proposed in the literature such as, for example, the smooth overlap of atomic positions and a previous approach based on unlabeled graph kernel introduced by Ferré et al. 38 . In the case of enthalpies of formations of solids (BA10 dataset), for small training set size -up to 300 datapoints -the accuracy of our new HMPP kernel is lower with respect to the SOAP approach. However, when the number of training configurations increases the learning ability of the SOAP descriptor saturates. The HMPP kernel does not suffer of this limitation and outperforms SOAP for 1000 and more training configurations. From these case studies the labeled graph kernel seems a particularly promising tool to improve the accuracy in machine learning regression problems in chemistry and materials science.

Data availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.