Efficient global structure optimization with a machine learned surrogate model

TL;DR

This paper introduces GOFEE, a global optimization method that combines machine-learned surrogate models with first-principles calculations to efficiently explore atomistic structures, significantly outperforming traditional algorithms.

Contribution

The paper presents a novel scheme integrating Gaussian Process-based surrogate modeling with active learning for efficient global structure optimization using first-principles energy calculations.

Findings

Outperforms traditional evolutionary algorithms by two orders of magnitude.

Successfully identifies surface reconstructions and initial oxidation stages.

Demonstrates efficiency in exploring complex atomistic potential energy landscapes.

Abstract

We propose a scheme for global optimization with first-principles energy expressions (GOFEE) of atomistic structure. While unfolding its search, the method actively learns a surrogate model of the potential energy landscape on which it performs a number of local relaxations (exploitation) and further structural searches (exploration). Assuming Gaussian Processes, an acquisition function is used to decide on which of the resulting structures is the more promising. Subsequently, a single point first-principles energy calculation is conducted for that structure. The method is demonstrated to outperform by two orders of magnitude a well established first-principles based evolutionary algorithm in finding surface reconstructions. Finally, GOFEE is utilized to identify initial stages of the edge oxidation and oxygen intercalation of graphene sheets on the Ir(111) surface.

Figures

Equations

$$E_{sur}({\mathbf{x}}_{\ast })$$

$${\sigma }_{sur}({\mathbf{x}}_{\ast }{)}^2$$

$$k(\mathbf{x},{\mathbf{x}}^{\prime })={\theta }_1\exp \left(-\frac{(\mathbf{x}-{\mathbf{x}}^{\prime }{)}^2}{2l_1^2}\right)+{\theta }_2\exp \left(-\frac{(\mathbf{x}-{\mathbf{x}}^{\prime }{)}^2}{2l_2^2}\right)$$

$$f(\mathbf{x})=E_{sur}(\mathbf{x})-\kappa \cdot {\sigma }_{sur}(\mathbf{x}),$$

Full text

In materials science and physical chemistry, the search for optimal structure is a recurring task, e.g. in describing crystalline defects, such as grain boundaries M van der Zande et al. (2013) and surface reconstructions Lazzeri and Selloni (2001); Merte et al. (2017), and in modeling heterogeneous systems such as binary compounds Flikkema and Bromley (2004); Ferrando et al. (2008) and supported nano-particles Zhai and Alexandrova (2018); Zandkarimi and Alexandrova (2019); Kolsbjerg et al. (2018). Depending on the search strategy and the complexity of a given problem, many thousands of energy and force evaluations may be required for the structural candidates in the course of the search. The results of these calculations constitute a set of structure-energy relation data points which represents a valuable resource that can direct the search. If the energy calculations are done at a first-principles (FP) level using density functional theory or quantum chemical methods, the computational bottle-neck lies in performing the individual energy/force evaluations and considerable speed-ups may be achieved by introducing machine-learning techniques that utilize this resource and provide tools to minimize the total amount of FP calculations.

As an example, a machine-learning technique may be utilized to estimate, based on the data set, the potential energy of each atom in a given structure at practically no computational cost Jacobsen et al. (2018). This information may then be used to guide further structure search, e.g. by shifting or replacing the least stable atoms, resulting in new structures that are more relevant to the search.

Greater is the potential for approaches where the entire search, including local relaxations, is performed on a machine learned surrogate energy landscape. As a first approach, the ML model can be trained in advance on a database of structures and FP properties. Here, regression models based on kernels, invariant polynomials and neural networks have all proven successful in a number of studies Behler and Parrinello (2007); Bartók et al. (2010); Shapeev (2016); Schütt et al. (2018) and are drastically changing the field of fitting force fields Zhai and Alexandrova (2016); Smith et al. (2017); Botu et al. (2017); Chmiela et al. (2017); Deringer and Csányi (2017).

Since the models are all interpolative and give reliable results only within their training domain, pre-fitted models have the drawback that they require the expensive generation of a large, diverse database of training data to be successfully applied to a structure search problem.

A more data efficient approach is to start from a small incomplete training database and then augment it on the fly only with the data deemed most relevant, at that time, for solving the task at hand e.g. structure optimization Zhai et al. (2015); Todorović et al. (2017); Yamashita et al. (2018); Deringer et al. (2018a, b); Tong et al. (2018); Bernstein et al. (2019); Gubaev et al. (2019); Van den Bossche (2019). This is the philosophy in the area of active learning Smith et al. (2018); Zhang et al. (2019). It was recently demonstrated in the context of an evolutionary algorithm (EA) structure search framework, where an artificial neural network was trained and used for local relaxation while the EA acted only on FP single-point energy evaluations Ouyang et al. (2015); Kolsbjerg et al. (2018). Active learning approaches have been extensively applied in molecular dynamics simulations Li et al. (2015); A. Peterson et al. (2017); Miwa and Ohno (2017); Podryabinkin and Shapeev (2017); Jinnouchi et al. (2019) with data efficient training databases as a byproduct. It has also been applied in local optimization problems such as local relaxation Garijo del Río et al. (2018); Denzel and Kästner (2018) and in minimum energy path determination with the Nudged Elastic Band (NEB) method Peterson (2016); Koistinen et al. (2017); Garrido Torres et al. (2019). Local optimization problems lend themselves particularly well to the construction of surrogate energy landscapes as a Cartesian coordinate representation of atoms may be adopted.

When a surrogate energy landscape is trained via active learning, the issue arises which next computationally expensive FP single point energy to evaluate. In this work, we present a strategy, GOFEE, that utilizes Bayesian statistics in the context of Gaussian Processes (GP) Rasmussen and Williams (2005). This framework allows for the estimation of the uncertainty in any prediction on the surrogate energy landscape, which provides the foundation for an acquisition function that guides the search. The method is first demonstrated to work up to two orders of magnitude faster for the identification of reconstructed surfaces of rutile $\text{SnO}_{2}(110)\text{-}(4\times 1)$ and anatase $\text{TiO}_{2}(001)\text{-}(1\times 4)$. Next, the method is used to tackle a set of otherwise prohibitively complex problems for graphene patches on Ir(111): the edge structure itself, the structure of the oxidized edges, and the pathways for oxygen intercalation at the graphene edges.

The degree of success achievable by any surrogate based search method is largely dependent on the quality of the surrogate model Deringer et al. (2018a, b). In this work, the GP regression method is adopted partly because of its tractable simplicity and partly because GPs are expected to behave well, as the number of training examples increase during the search, due to the adaptive predictive power inherent to non-parametric methods.

The task of GPs is to infer a distribution over functions, here $p(E_{sur}|X,\mathbf{E})$, that is consistent with a training set of observed atomic configurations $X=(\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N})^{T}$ and their corresponding energies $\mathbf{E}=(E_{1},E_{2},\dots,E_{N})^{T}$. To include the symmetries of the system, $\mathbf{x}_{i}$ is taken to be the feature vector for the i’th configuration rather than the Cartesian coordinates. We adopt the global fingerprint feature from Oganov and Valle Valle and Oganov (2010), however the method is expected to work equally well with other features. For GPs, the distribution is assumed to be normal, which enables estimation of not only the energy $E_{sur}(\mathbf{x})$ as the mean of the distribution, but also the predictive uncertainty $\sigma_{sur}(\mathbf{x})$. As will be discussed later, the predictive uncertainty is useful in a search context, as it allows the distinction between explored and unexplored regions of the search space.

A Gaussian process is specified by it’s prior mean function $\mu(\mathbf{x})$ and covariance function $k(\mathbf{x}_{i},\mathbf{x}_{j})$, which encodes prior assumptions about the target function. Given these, energy and uncertainty predictions for a new structure $\mathbf{x}_{*}$ is carried out using

[TABLE]

where $K=k(X,X)$ and $\mathbf{k}_{*}=k(X,\mathbf{x}_{*})$ and the target function is assumed noisy with uncertainty $\sigma_{n}^{2}=10^{-5}$\text{,}\mathrm{eV}$^{2}$, which acts as regularization. To include the repulsive atomic core generally present, the prior mean function is taken to be a conservatively chosen repulsive interatomic potential, specifically $\mu(\mathbf{x})\propto\sum_{ij}(0.7*r_{CD,ij}/r_{ij})^{12}$, where $r_{ij}$ and $r_{CD,ij}$ are the distance and covalent distance between atom $i$ and $j$. This is especially beneficial in a structure search context, where the fine details of the repulsive part of the potential is not crucial, unlike the near equilibrium part of the potential. The covariance function was chosen to be a sum of two Gaussian covariances

[TABLE]

one with a large characteristic length scale $l_{1}$ carrying most of the weight $\theta_{1}=(1-\beta)\theta_{0}$, and one with a smaller length scale $l_{2}$ and less weight $\theta_{1}=\beta\theta_{0}$, where $\beta=0.01$ and $\theta_{0}=\mathbf{y}^{T}(K+\sigma_{n}^{2}I)^{-1}/N$. In the $\text{SnO}_{2}$ and $\text{TiO}_{2}$ application we use $l_{1}=10$ and $l_{2}=2$ whereas for the graphene edge $l_{1}=10$ and $l_{2}=0.5$ were used. This enhances the models capacity to capture both the long scale features of the energy such as large energy funnels, as well as allowing for sufficient resolution of local minima.

Figure 1 shows, for a simple problem how the surrogate energy landscape improves as more data is added to the training set. The system considered is naphthalene, constrained, for illustrative purposes, to change only according to the two coordinates specified. The resulting 2d slice of the full energy landscape, contains four local minima including naphthalene itself. With only a few training examples near one local minimum, the model is able to predict the locations of the remaining minima to approximately coincide with those of the true energy landscape. In the search we will take advantage of this, and conduct most of the search, specifically all local relaxations, in the surrogate energy landscape, which is orders of magnitude faster than FP calculations. As illustrated in the figure, a structure relaxed with the model can then be evaluated with a single FP calculation and used to update the model.

Relying entirely on the surrogate model to guide the search has the drawback that the data collection process, vital to actively improving the model, is itself model dependent. This interplay has a tendency to cause under-exploration of the search space in turn leading to premature stagnation of the search. The minimum belonging to naphthalene in Fig. 1 is an example that the true depth of a minimum might be underestimated until appropriate data has been collected. To remedy this problem we bias data collection towards unexplored regions of the search space, using the predictive uncertainty $\sigma_{sur}(\mathbf{x})$ as a natural way to quantify this. Data collection is then performed according to an acquisition function $f(x)$, relying on both the predicted energy and uncertainty. There exist multiple choices for such an acquisition function Wang et al. (2017), expected improvement and probability of improvement are some, as well as the lower confidence bound

[TABLE]

used in this work due to its simplicity. Here $\kappa$ is a tunable parameter determining the emphasis on the predicted uncertainty and thus the degree of exploration in the search.

The surrogate model is central to the GOFEE search method, which is initialized with a small set of randomly generated structures, for which the FP energy is evaluated. They make up the first structures in a database used for training the surrogate model, and to which all subsequent FP evaluated structures are added.

A diagrammatic layout of GOFEE is given in Fig. S1. Training the surrogate model is the first step in each search iteration after which a diversified population of the best structures, currently found in the search, is used to generate a set of new candidate structures using simple rattle mutations as in Monte Carlo and evolutionary search strategies. To take full advantage of the computational inexpensiveness of the surrogate model, multiple new candidates are generated and relaxed in each search iteration instead of just one as shown for the example in Fig. 1. From all these relaxed candidates only the single most promising, as estimated by the acquisition function, is evaluated with FP. To accommodate some force information without training on forces directly, this structure is perturbed slightly in the direction of the force and a second FP calculation is performed on the resulting structure. A search iteration is concluded by adding these two structures to the training database.

As a first example using GOFEE we considered the $\text{SnO}_{2}(110)\text{-}(4\times 1)$ surface, for which the global minimum structure is known Merte et al. (2017) and is shown as an inset in Fig. 2a. The figure also shows the cumulative success curves for finding the global minimum with this method as well as with the well established EA originally used to find the structure Vilhelmsen and Hammer (2014). Noting the broken axis, the figure shows a two orders of magnitude decrease in the number of FP calculations required to reach, e.g., $80\%$ success. This is largely attributed to the fact that this method relies only on FP for single point calculations. To show the effect of the exploration promoting parameter $\kappa$, a specific search instance for each of $\kappa=1,2,4$ is shown in Fig. 2c. In the figure, principal component analysis (PCA) is used to project all structures visited in each of the three searches onto the same two principal components determined from a large set of structures. The chosen search instances showcase common behavior for the three values of $\kappa$. For $\kappa=1$ the search has a tendency to over-exploit local minima, and as a result get stuck in a local minimum before reaching the global minimum. For $\kappa=4$ the opposite is true and the search will superficially explore many local minima before starting to optimize the best of these. The $\kappa=2$ search represents the optimal compromise between the two, performing a necessary but sufficient amount of exploration before settling to optimizing the global minimum. In all three examples it is apparent that high energy structures are primarily sampled in the beginning of the search, when the surrogate model is still learning the rough features of the energy landscape. GOFEE was similarly applied to the $\text{TiO}_{2}(001)\text{-}(1\times 4)$ surface reconstruction problem, displaying the same degree of improvement as compared to the EA. The results are shown in Fig. S2.

To demonstrate the versatility of the GOFEE method we proceed to address the hitherto prohibitively complex problems of determining the edge structure of graphene patches on Ir(111) as well as that of the oxidized edge. The resulting structures are used to study the atomistic mechanisms involved in intercalation of oxygen in the system. The intercalation process has been intensively studied experimentally for Ir(111) Grånäs et al. (2012); Martínez-Galera et al. (2016); Larciprete et al. (2012) and involves dissociative adsorption and diffusion of oxygen as well as penetration of the graphene edge. Although experiments suggest Grånäs et al. (2012) that the limiting step for the intercalation process is this edge penetration, it is not well understood.

In our contribution to fill this gap, intercalation through the non oxidized edge was first considered. Figure 3a shows the most stable edge structures found, using our search method, when varying the number of carbon atoms present in the cell. Using the energy of carbon within the graphene patch as reference, the energies are compared in Fig. 3b, showing that the preferred structure is not the perfect edge (A), but instead the structures with one (A-1) and three (A-3) carbon atoms less on the edge, both of which feature pentagonal rings (see Fig. S3). This can be attributed to the fact that these structures avoid forming the unfavorable carbon iridium bond in the position of largest mismatch between the periodicities of the graphene edge and the iridium surface. Although this effect does cause small gaps in the graphene edge, it is not enough to allow for oxygen intercalation, with the calculated energy barriers being larger than $2\text{\,}\mathrm{eV}$ for all structures.

The structure of the oxidized graphene edge was also considered, as oxygen is naturally present during the intercalation process. Searches are performed with up to three oxygen atoms in the cell. The resulting structures and energies are depicted in Fig. 3c (B-D), Fig. S4 and Fig. 3d. They display a preference for oxidizing the edge with the oxidized region partially detaching from the surface when two or more oxygen atoms are present. This results in a significant gap in the edge likely of accommodating intercalation. Figure 3c (E-G) (Fig. S4) further shows the structures resulting from extending this trend up to six oxygen atoms. For the energies, atomic oxygen adsorbed on the iridium surface is used as the reference. Based on the energies, the size of the gaps are thermodynamically self limiting, as edge oxidation is only thermodynamically favored up to four oxygen atoms. Further oxidation requires breaking of increasingly strong C-Ir bonds.

To study whether these oxidized edge structures are likely to form and contribute to the intercalation process under typical experimental conditions, lowest energy paths were calculated using the climbing image elastic band (CI-EB) NEB type method Kolsbjerg et al. (2016). Figure 3e and f, (Fig. S5-6) respectively shows the energy profiles for inserting the first and third oxygen to the edge revealing the third oxygen to be the more expensive of the two with an energy barrier of $1.23\text{\,}\mathrm{eV}$. However oxygen intercalation experiments typically feature large oxygen coverages on the iridium surface, resulting in weaker bonding of the adsorbed oxygen as this coverage builds up. This effectively lowers the barriers, as the transition state for binding to and opening the graphene edge is expected to remain unchanged. The effect is depicted in Fig. 3f showing how the initial state energy, of the oxygen being inserted, is increased when sharing a single iridium atom with a neighboring adsorbed oxygen atom. Figure 3g (Fig. S7) shows the energy profile for the intercalation of an oxygen atom through the edge gap of structure (E), displaying a barrier of $0.93\text{\,}\mathrm{eV}$, which will also be lower at realistic oxygen coverages as discussed above.

In conclusion we have formulated a machine learning enhanced structure search method and used it to solve a previously prohibitively hard problem. This has provided insight into the atomistic structure of graphene island edges involving pentagonal rings and into the atomistic mechanisms of oxygen intercalation for graphene on Ir(111).

GOFEE search method

A sketch of the GOFEE search method is shown in Fig. S1 and depicts the key elements in the search.

$\text{TiO}_{2}$ benchmark and scaling

Improved structure search methods allow for the determination of increasingly complex structures within an appealing time frame. To show that this is indeed the case for the method presented in this work, the method is applied to the anatase $\text{TiO}_{2}(001)\text{-}(1\times 4)$ surface reconstruction and the complexity of the problem is increased by optimizing from one to three layers on top of a fixed layer as shown in the insets of Fig. S2. As for $\text{TiO}_{2}$ in the main text we compare to the well established evolutionary algorithm (EA) [L. B. Vilhelmsen and B. Hammer, J. Chem. Phys. 141, 044711 (2014)] compared to which it is orders of magnitude faster and handles the increased complexity better. The method is also compared to the same EA for which machine learning is used in the form of clustering to improve the candidate generation step as in [K. H. Sørensen, M. S. Jørgensen, A. Bruix, and B. Hammer, J. Chem. Phys. 148, 241734 (2018)]. This results is a significant improvement of the EA, but does not come close to the improvement achieved by GOFEE which uses machine learning to avoid the large number of first principles calculations traditionally spent on local relaxation.

Graphene edge structure

Larger plots as well as top and side views of the graphene edge structures presented in the main article is shown in Fig. S3 and Fig. S4. In addition the presented minimum energy profiles are shown in Fig. S5, Fig. S6 and Fig. S7 with snapshots of structures along the pathways.