Modelling the dielectric constants of crystals using machine learning

The relative permittivity of a crystal is a fundamental property that links microscopic chemical bonding to macroscopic electromagnetic response. Multiple models, including analytical, numerical and statistical descriptions, have been made to understand and predict dielectric behaviour. Analytical models are often limited to a particular type of compounds, whereas machine learning (ML) models often lack interpretability. Here, we combine supervised ML, density functional perturbation theory, and analysis based on game theory to predict and explain the physical trends in optical dielectric constants of crystals. Two ML models, support vector regression and deep neural networks, were trained on a dataset of 1,364 dielectric constants. Shapley additive explanations (SHAP) analysis of the ML models reveals that they recover correlations described by textbook Clausius-Mossotti and Penn models, which gives confidence in their ability to describe physical behavior, while providing superior predictive power.


I. INTRODUCTION
The dielectric response function is one of the fundamental properties of material, which can give an insight into optical and electric properties. Specifically, the electronic highfrequency component of the dielectric constant has been explored intensively in early studies and is still of interest among researchers today. 1,2 Numerous efforts have been made to model dielectric constants. Most notably, the Clausius-Mossotti and Penn models, and their variants, are widely adapted throughout the literature. [2][3][4]

A. Clausius-Mossotti Model
The Clausius-Mossotti (CM) equation expresses the dielectric constant ε as: Here, α i is polarizability of atomic species i, and v is the volume. In the case of molecular crystals, α i can be assigned to a constituent molecule, [5][6][7] and in ionic solids it is assigned to an ion. The inherently many-body nature of the dielectric constant is reduced to this simple relation through employment of two large assumptions.
a) Electronic mail: keith.butler@stfc.ac.uk b) Electronic mail: a.walsh@imperial.ac.uk Firstly, that an external electric field is screened by a dielectric medium before reaching an atom. In the CM model, cubic symmetry is assumed and this local field E loc (r) is expressed as where E ext (r) is the external electric field. This relation, the Lorentz relation, is frequently used for non-cubic cases and generally holds when the anisotropy is small. 8,9 Secondly, the atomic polarizability α i is assumed to be a constant additive quantity. In other words, α i is assumed to be unaffected by the surroundings. Due to its simple form, the CM model is still being used and in practice works for many materials, 2,10 despite the underlying assumptions. 8,[11][12][13] Numerous efforts have been made to improve the CM model. [14][15][16] In general, more effort was put into improving atomic polarizability α i rather then changing the functional form. 14-16 While Wilson and Curtis increased the accuracy of the model by considering electrostatic environment of an atom, 15 Jemmer et al. fitted parameters to Møller-Plessett calculations. 16 Nonetheless, the range of compounds and structures in which the model is valid was found to be limited. More recently, Shannon and Fischer have fitted atomic polarizabilities to sets of experimental and computational data, consisting of 650 values from 487 compounds. 2 This was one of the first studies to consider compounds throughout a wider variety of structures, and their values of atomic polarizability were more transferable. However, due to the limited dataset, some parameters were manually fitted using chemical intuition. Despite the success of the CM based models for many compounds, 2,15,16 they are far from being universal and poorly describe covalent systems.

B. Penn Model
The premise of the Penn model is the electronic band picture. When the overlap between electrons at different atomic sites is large, such a model is anticipated to be favorable. Treating the external electric field as a time-dependent perturbation, Penn derived the following relation: 1 where ω p is the plasma frequency and E g is the width of the band gap. The significance lies in the fact that the dielectric constant is described only by the plasma frequency and the band gap. Using the second order perturbation theory, polarizability α for an isolated atom can be written α ≈ ( ω p /E g ) 2 , where E g is the HOMO-LUMO gap. In this case, the screening of the external field can be ignored and Penn model is equivalent to the CM model.
As with the CM model, there are various assumptions for the Penn model. The first main assumption is to consider the electrons as a free electron gas. This allows the form of the equation to be greatly simplified, but is not valid when electrons are localized. 17,18 The second assumption is to consider only the bandgap and not the density of states around the valence and conduction band edges. This flat band approximation results in a large error when the valence or conduction band has a large dispersion. 4 Various efforts have been made to improve this model. Phillips added a correction for ionic bonds and showed that dielectric constants in zinc-blende and wurtzite were reproduced. 17,18 Furthermore, Van Vechten added a scaling factor taking into account that not all valence electrons contribute to the dielectric response. 19,20 These modifications did contribute to increasing the predictive power; however, a universal model was not achieved.

C. Data driven models
The CM model approaches the problem from a molecular picture, whereas the Penn model is based on the electronic band picture of a crystal. As such, the former has difficulty with delocalized electrons and the latter has difficulty with localized electrons. In both cases, the number of model parameters were increased in order to improve transferability. In all cases, this was done empirically and the physical justification was often limited.
Rapid progress in computer technology and first-principles modelling techniques has enabled the calculation of dielectric constants for large number of materials. [21][22][23] The modern calculation of dielectric constants are typically done using density functional perturbation theory (DFPT). 3,24-28 Therefore, many studies have involved analysis of datasets of dielectric constants aiming to understand trends or predict values for new materials. For example, Han and coworkers have performed high-throughput calculation of binary and ternary inorganic compounds, and compared the relation between the band gap and the dielectric constant. 3,24 Some studies employ machine learning (ML) methods where a statistical model is trained. 25,26 For example, Umeda et al. trained different ML models on a dataset of 3,382 compounds. 25 They obtain a good agreement with DFPT calculations; however, the reasons for this agreement are not discussed.
Overall, simple analytical models are limited by their generality, whereas ML models are difficult to interpret. By studying the difference between ML and analytical models, we may recover new underlying physics and deepen our understanding of the phenomenological behavior of the dielectric response.
In this study, we use two types of ML models and investigate the reasoning behind their predictions. Firstly, a deep neural network and support vector regression model are trained on the same dataset. Next, their results are compared with DFPT calculations. Lastly, we perform Game Theoretic analysis (Shapley additive explanations) to elucidate the characteristics of the trained ML models and to investigate the reason behind their predictions. We compare the result closely with the CM and Penn models and show that ML can indeed learn underlying physical trends, while having superior predictive power. Additionally the Shapley analysis allows us to identify reasons for poor model performance in particular cases, which is important when considering how relaible predictions are.

A. Machine learning
Supervised ML models were trained on a dataset of 1,364 dielectric constants. The dataset was prepared by combining two pre-existing datasets, 27,28 and averaging the data for the overlapping materials. Materials with a band gap less than 0.5 eV were removed because small gap materials require very dense sampling of the Brillouin zone, which is difficult to realize in high-throughput calculations. 33 The dielectric constant was calculated by taking a diagonal average of the electronic part of the dielectric matrix. Since this dataset only contained minimum features, it was augmented and processed using Materials Project, 29,30 pymatgen, 32 SMACT, 31 and scikit-learn. 34 The dataset was split into 8:2 the training:test sets. The same data was used when training the two different ML models.
The first ML model was support vector regression (SVR). Before training, we analyzed the available features and removed those unimportant or similar. This step was necessary, because having too many features increases the dimensions of the manifold that the model must learn in and thus lowers the performance. 35 Specifically, we calculated feature importance using the random forest ( Fig. S1) and removed low importance features: Space group (one hot encoded), atomic species (one hot encoded), number of elements, and number of sites inside the unit cell. Next, with the criterion of r 2 > 0.90, we removed maximum oxidation state and variation in Maledung energy, which had high correlation with variation in oxidation state and minimum in Madelung potential, respectively (Fig. S2). The refined features are presented in Table I.
Four machine learning models, the random forest, gradient boost regression, kernel ridge regression, and SVR, were trained and compared ( Fig. S4 and Table S1). The former two models are based on decision trees and latter two models are based on kernel methods. Since the random forest and gradient boost regression are ensemble methods, we surmise their performance to be better, however, they were overfitting and the kernel methods were able to generalise better. The kernel ridge regression had slightly lower performance compared to SVR, therefore SVR was employed for this study. The second ML model was a deep neural network (DNN). Specifically we have adopted the DNN architecture MEGNet developed by Chen, et al. 36 MEGNet can be trained only by using the crystal structure of materials. The exact structure of the network is presented in Fig. S3. It is not trivial to uniquely express a crystalline system in the form of a vector, however MEGNet overcomes this difficultly by representing bonding networks as graphs and using the set2set algorithm to consistently treat different sized graphs. 36 Information of atomic number and bond length are encoded into the graph representation. The model was trained for 400 epochs and the layer weights of the epoch epoch with the smallest Huber loss were employed.
There are two advantages and a disadvantage of using MEGNet over SVR. The first advantage is to avoid so called "feature engineering". Feature engineering is a process of feature selection, which was necessary for SVR. Although this can be a way of inputting our domain knowledge, 21 this relies heavily on experience and intuition and may obstruct systematic improvement of an ML model. The second advantage is that it is easy to perform transfer learning. As larger datasets are available for other materials properties, it is tempting to use that information in training. Since the upper layers in DNNs are known to learn general trends, it is possible to import the layer weights of the upper layers in advance and improve the accuracy. 36 Taking into account the Penn model, and the results of feature engineering for SVR (Figs. S2 and S1), we performed transfer learning from MEGNet model Chen et al. trained on band gaps. 36 The disadvantage of DNNs is that interpretability. Specially, because we used graph representation as the input, the relation between the features and dielectric constant were difficult to extract.

B. Density functional theory calculation
In order to validate our results, we have performed DFPT calculations on a total of 24 structures. These structures were the materials with largest error in the SVR and DNN models. The calculations were performed with projector-augmented wave scheme as implemented in VASP. [37][38][39] The structures were taken from Materials Project and were calculated using the PBE functional. 29,30,40 The reciprocal space were sampled so that the spacing between the k-point was about 2π × 0.03Å −1 . The energy cut off was set to at least 600 eV and the wavefunctions were optimized to a tolerance of 10 −6 eV.

C. Shapley additive explanation
For model interrogation, we perform Shapley additive explanations (SHAP) analysis. 41 The Shapley regression value is defined as: Here, S is a subset of the features (F ), i is a particular feature of interest, f is the ML model, and x S represent values of the input features in the set S. 41 φ i describes how much the model output changes when feature i is added to the model. Therefore, it can be used to quantify feature importance. As SVR cannot take features of different size, we used the median of the feature instead of removing them. Furthermore, if we fix data x S and model f , we can show which features are responsible for the prediction given subset x S . Using the additivity approximation suggested by Lundberg and Lee, 41 we calculated SHAP values of all the features and data in the dataset.

A. ML prediction
The dielectric constants predicted by SVR are plotted against the dataset values in Fig.  1(a). In general, materials with dielectric constants above 10 had larger error. DNN performed similarly to SVR as shown in Fig. 1(b). The main difference was suppression of error in the high dielectric constant range and slight increase of error in the low dielectric constant range in the predictions from the DNN.
Performance metrics of SVR and DNN are summarized in Table II. This performance is similar to the performance of ML models in previous studies. [42][43][44] Ideally, if there is no overfitting and the ML model is trained sufficiently, metrics for training and test data should be identical, however Table II imply tendency for both models to exhibit overfitting. For SVR, different sets of hyperparameters and loss functions were trialled and the result presented here represents the best achievable performance (Table S1).We built a series of SVR models on different sized subsets of the data (Fig. S6) and find that the difference in the loss function metric decreased as we increased the dataset size, but did not converge. We therefore attribute overfitting to the small size of our dataset. Compared to SVR, the difference in metric between the test and training data were smaller for DNN, suggesting that the latter model was able to generalize better, and exhibits overfitting to a lesser extent. For test data, the DNN demonstrated higher performance for all of the metrics than the SVR (Table II).
The 10 materials with largest prediction error are listed in Tables III and IV. Out of the 10 materials, seven of them were common amongst the two ML methods. It is worth noting that LiBC, Ga 2 Te 5 , LiAsS 2 and In 2 HgTe 4 have error large enough to be identifiable in Fig. 1. We also present the energy above hull (obtained from Materials Project) to show that these materials are not artefacts from computational screening and could possibly be stable. 29,30 Since the dielectric constant dataset we use for training was derived from a highthroughput workflow, 27,28 the precision is optimized to be convergent for most materials, 45 which results in some results being poorly converged. 33 Therefore, for comparison, we have performed higher precision calculations using DFPT on these compounds and the results are presented in Table III and IV.  TABLE III. List of 10 material with largest error in support vector regression (SVR) prediction. εSVR, εDataset and εDFPT(this work) is dielectric constant from SVR prediction, training dataset and our density functional perturbation theory (DFPT) calculation, respectively. ∆ε dataset is the difference between εSVR and εDataset. ∆εDFPT is the difference between εSVR and εDFPT. ∆E hull (meV) is energy above hull per atom. The larger error found when the large dielectric constant is greater is not surprising (Fig. 1). In addition to the numerical instability in DFPT calculations, 33 there are fewer materials with large dielectric constants in the dataset, therefore the ML models were not able to fully learn the trends in the large dielectric constant range.
Comparing SVR and DNN, the performance was similar, as shown in Fig. 1 and Table  II. This is surprising since the features used for training each model were different. It suggests that both models are sufficiently capturing physical trends in the training data. If the model is generalized well, we expect them to be able to detect anomalous data in the DFPT dataset. For eight (seven) materials out of 10 materials, SVR (DNN) was actually predicts values closer to our calculations, which were performed under higher precision. This is clearly the case for Ga 2 Te 5 and LiAsS 2 . For example in the case of Ga 2 Te 5 in Table  III, the SVR prediction was 10.48, whereas the dataset value was 17.19 and our calculation was 13.75. Since small band gap materials require fine Brillouin zone sampling, 33 the standardized sampling in the high-throughput calculation setup may be insufficient. To confirm this, we calculated the dielectric constant of LiAsS 2 with different Brillouin zone sampling density. Our converged value of dielectric constant was 8.57, while 11.46 was obtained using coarse sampling. 11.46 is closer to the dataset value of 12.39 (Tables III and S3). Although this is not a direct evidence, we speculate that materials which exhibited a large difference between our calculations and the dataset were especially sensitive to the Brillouin zone sampling and as a result, made the reported dielectric constants of these materials anomalous.
B. SHAP analysis SHAP analysis applies a game theoretic approach to calculate the importance of individual input features to a given model prediction. A positive (negative) SHAP value indicates that a given feature contributes to an increase (decrease) in the prediction with respect to the mean of the set. Fig. 2 shows calculated SHAP values for all features and data. The features are ordered by their importance. Note that the order of importance differs from Fig. S1, which was a random forest model trained on all features and data, whereas Fig. 2 is SVR built only on training data.
The high contribution of material density can be explained by both the CM and Penn expressions. If there are more electrons in a given volume, the dielectric response will become larger, and indeed SHAP analysis shows that dielectric constant monotonically increases with density.
The high contribution of the band gap is also no surprise from the Penn model (Eq. 3). Lower energy excitations result a larger dielectric constant. A large band gap gives a negative SHAP contribution. Interestingly, the magnitude of positive contributions from small band gaps has a longer tail in the distribution than the negative contribution from large gaps. This suggests that although a large gap decreases the dielectric response, this effect diminishes quickly so that continuously increasing band gap will not always decrease the dielectric constant.
To take a deeper look at these relationships, further analysis was performed in Fig. 3(a). The contribution from band gap decreases as the gap increases with an inverse power law relation. By coloring the individual points according to material density, the interplay of density and band gap can be observed. Specifically when the band gap is low, a low density increases the SHAP value, while when the gap is high a low density decreases the magnitude of the negative SHAP value. While it is not possible to derive rigorous analytical relationship from number of data points we have, Fig. 3(a) suggests a φ Eg ∼ n/E m g relation, where φ Eg is the SHAP value for the band gap, n is the material density, and m is an arbitrary constant. This is similar to the Penn relation (Eq. 3). Therefore, it is possible to interpret that ML is learning the Penn model while incorporating other feature relations as correction terms.
When the formation energy is high (low) the dielectric constant is raised (suppressed),  Table III. Starting with the dielectric constant predicted by a median of the features (ε = 5.58), each feature contributes to lowering or elevating the value of the dielectric constant in an additive manner. The final value of the dielectric constant ε pred is the output of the model. Fig. 2. This relation is opposite to that of the band gap and the negative correlation −0.65 between the two also agree with this trend (Fig S2). Within our ML model, two roles of formation energy can be suggested. Firstly, it is acting as a ionicity parameter. According to Pauling, the formation energy of binary compounds is ∆H ∝ (∆E Pauling ) 2 , where ∆E Pauling is the Pauling energy difference. 46 This trend was present in our dataset (Fig. S7); however, the variance was large which suggests contributions from other features. Secondly, that it is another feature for the band gap. In main group binary compounds, the relation ∆H ∝ (E 2 g /E f ) ln(E g /E f ) has been reported. 47 Quantitatively, this relation did not hold for our case (Fig. S7). Given these two relations and accounting for the fact that formation energy did not show obvious relation between the Pauling energy difference and the band gap, we suggest that it contains both information weakly.
The variation of oxidation states also exhibited high importance, Fig. 2. Since oxidation states have information about ionicity, we may expect it to lower the value of the dielectric constant when it is large. 19,20 Fig. 2 shows the opposite trend. This counter intuitive result can be explained by interpreting this feature as a correction the flat band approximation, described earlier for the Penn model. As the band gap is defined as a difference between the valence band maximum and the conduction band minimum, it does not have information about band dispersion. If ionic and covalent compounds with the same band gap exist, ionic compounds will have smaller effective band gap. When band gap is large (small) and contributes to making dielectric constant smaller (larger), lower variation in oxidation states amplifies (suppresses) this effect (Fig. S8). Theoretical studies also support the importance of this correction. 4 Other oxidation state and Madelung energy features did not show clear trends, which constitutes further evidence that the variation of oxidation state is acting as an ionicity parameter. This result would not have been available if we were only considering feature importance, and highlights the ability of SHAP to access deeper understanding of trained models.
Since the distribution of ionic species ratio features had broad SHAP value distributions and were difficult to interpret in Fig. 2, we continue their discussion based on Figs. 3(bd). SHAP values of O 2and Fratio are plotted against material density. In general, higher concentration of these ionic species tend to reduce the dielectric constant. This effect could be interpreted as originating from the strongly electronegative nature of O 2and F -, which thus stabilize valence electrons and suppress their dielectric response. In contrast, for the Li + cation, the SHAP value is positive, meaning it increases the dielectric constant. Furthermore, we can see the interplay of density and composition; when the density is higher, the effect of O 2and Fratio was enhanced, where O 2had slightly larger change. Due to the limited number of data points, further study is required to draw firmer conclusions.
Finally, we explore the limitations of the trained SVR model by examining LiBC and CdCN 2 in Table III. Since SHAP values are additive, 41 we are able to decompose them to each of the features. We begin with the dielectric constant predicted for the median of features (ε = 5.58) and add contribution of all the features sequentially (Fig. 4). Although only the 10 most important features are presented, other features had minor contribution ∼ 0.001, so we will only discuss features presented in Fig. III. For LiBC, which had largest error for both SVR and DNN, SVR predicted 3.54 where our calculated value was 11.96 (Table III). As shown in Fig. 4(a) the largest contribution was from the material density. Together with the fact that the material density was small (2.15 g/cm 2 ), the error is likely due to the incorrect contribution of material density. It is likely that in this range of material density, sufficient sampling was not achieved with our training dataset.
Most of the materials with large error feature a large contribution from material density. The exception was CdCN 2 . SVR predicted 6.55, where 4.45 was the result of our DFPT calculation (Table III). As shown in Fig. 4(b), the three largest contributions are N 3ratio, formation energy, and variation in oxidation state, which have reasonable values. We instead assign the error as arising from the large anisotropy, which is not taken into account. The full dielectric matrix is: The DNN performed better because the bond lengths are explicitly taken into account through the graph representation.

V. CONCLUSIONS
We showed that two different machine learning models, a support vector machine and a deep neural network, were able to predict the dielectric constants of crystals with reasonable precision. Comparison with our density functional perturbation theory calculations reveals that the machine learning models were able to detect erroneous results in the original dataset. We performed SHAP analysis of the support vector machine model, which illustrated that it is learning similar relations to the textbook Clausius-Mossotti and Penn expressions. Finally, we showed the limitation of our support vector model through detailed analysis of the predictions.
We suggest that as long as the dataset is sufficiently large to sample the crystal space of interest, machine learning models can be an effective approach not only for prediction of material properties but also to capture physical trends. The analysis approach used in this study is not restricted to dielectric response, and therefore has application potential for other relations, including properties that are not as intensively studied and where there is no existing analytical description.

VI. DATA AVAILABILITY STATEMENT
A set of electronic notebooks to reproduce the model training and analysis performed in this study are available from [DOI to be inserted].