De novo exploration and self-guided learning of potential-energy surfaces

TL;DR

This paper introduces a largely automated method for building machine learning interatomic potentials that explore and fit potential-energy surfaces de novo, using a configuration-averaged kernel metric to select relevant structures efficiently.

Contribution

The authors develop a protocol that automates the creation of ML potentials from scratch, reducing manual intervention and enabling broader application in materials science.

Findings

Accurate potentials for diverse materials achieved with few DFT calculations.

Automated structure selection improves robustness and efficiency.

Method applicable to various chemical environments and coordination states.

Abstract

Interatomic potential models based on machine learning (ML) are rapidly developing as tools for materials simulations. However, because of their flexibility, they require large fitting databases that are normally created with substantial manual selection and tuning of reference configurations. Here, we show that ML potentials can be built in a largely automated fashion, exploring and fitting potential-energy surfaces from the beginning (de novo) within one and the same protocol. The key enabling step is the use of a configuration-averaged kernel metric that allows one to select the few most relevant structures at each step. The resulting potentials are accurate and robust for the wide range of configurations that occur during structure searching, despite only requiring a relatively small number of single-point DFT calculations on small unit cells. We apply the method to materials with diverse chemical nature and coordination environments, marking a milestone toward the more routine application of ML potentials in physics, chemistry, and materials science.

Tables

Table 1. Table 1: Hyperparameters for descriptors that we use in GAP fitting. For all descriptors: Gaussian width σ at subscript 𝜎 at \sigma_{\rm at} (squared-exponential kernel for 2- and 3-body; atomic density width for SOAP); number of sparse points N sp subscript 𝑁 sp N_{\mathrm{sp}} . For SOAP only: number of radial functions n max subscript 𝑛 max n_{\mathrm{max}} and angular momenta l max subscript 𝑙 max l_{\mathrm{max}} , and kernel exponent ζ 𝜁 \zeta . Cutoffs r cut subscript 𝑟 cut r_{\mathrm{cut}} are expressed in terms of the characteristic radius r 𝑟 r , listed for each material in the Interatomic potential fitting subsection.

						$r_{\mathrm{cut}}$ (Å)
	${\sigma }_{\text{at}}$ (Å)	$N_{\mathrm{sp}}$	$n_{\max }$	$l_{\max }$	$\zeta$	(covalent)	(metallic)
2-body	0.5	30				$9.0r$	$8.2r$
3-body	1.0	100				$2.925r$	$2.665r$
SOAP	0.75	2000	8	8	4	$4.5r$	$4.1r$

Table 2. Table S1: Additional information regarding the selected structures shown as insets in Fig.  5 : external pressure applied for this relaxation trajectory, p ext subscript 𝑝 ext p_{\rm ext} ; energy relative to the respective ground state, Δ ​ E Δ 𝐸 \Delta E ; maximum DFT-computed force component F i subscript 𝐹 𝑖 F_{i} in this structure.

		$p_{\mathrm{ext}}$	$\Delta E$	max $\left\{F_i\right\}$
		(GPa)	(eV / at.)	(eV / Å)
(i)	C (graphite)	0.111	+0.15	0.302
(ii)	C (buckled)	0.088	+0.23	1.678
(iii)	Si (dist. lon)	0.048	+0.32	0.742
(iv)	Si (unj)	0.000	+0.06	0.022
(v)	Si ($s{p}^3$ network)	0.124	+0.29	0.962
(vi)	Si (simple hex.)	0.208	+0.25	0.653
(vii)	Ti ($\omega$ phase)	0.185	+0.03	0.000
(viii)	Ti (hcp)	0.324	+0.05	0.238
(ix)	Si (RSS intermed.)	0.172	+0.59	0.857

Table 3. Table S2: Lattice parameters and atomic coordinates for silicon structure (v), as added to the GAP-RSS reference database in the final iteration (Fig.  5 ).

	$a=b=c=5.5869$ Å
	$\alpha =\beta =\gamma ={109.08}^{\circ }$
Si	-1.68605982	-0.41624606	3.95802416
Si	0.15435366	4.30188720	0.38981490
Si	4.01815977	-0.40190360	1.54106451
Si	-2.28421007	2.67006652	2.34218566
Si	3.26751299	2.67507753	-0.00486561
Si	1.46783885	-1.91088147	3.46795123
Si	1.26796258	1.77652573	3.00303591
Si	0.48155278	0.67469729	1.14050708

Table 4. Table S3: As Table S2 but for silicon structure (ix).

	$a=b=c=6.5595$ Å
	$\alpha ={84.76}^{\circ }$	$\beta ={130.49}^{\circ }$	$\gamma ={116.23}^{\circ }$
Si	-1.20744491	4.11837921	3.06869747
Si	1.15061078	2.70140772	0.21706451
Si	-2.74075733	-1.74587757	4.91801644
Si	-2.97936509	2.43953892	3.67216632
Si	-1.80907975	3.75555505	5.55352463
Si	-2.12808832	-0.09478641	3.03665790
Si	1.86693516	0.25171580	3.77895234
Si	-0.51394068	1.68239975	3.54412519
Si	3.12840615	4.18317367	0.57109224
Si	2.03701947	0.79695955	1.44756468

Equations

$$d_{\mathrm{AB}}=\sqrt{2-2k_{\mathrm{AB}}},$$

$${\rho }_i(\mathbf{r})=$$

$$=$$

$${\bar{c}}_{nlm}=\frac{1}{N}\sqrt{\frac{8{\pi }^2}{2l+1}}\sum _i{c}_{nlm}^{(i)},$$

$${\bar{p}}_{n{n}^{\prime }l}=\sum _m{\left({\bar{c}}_{nlm}\right)}^{\ast }{\bar{c}}_{n^{\prime }lm}.$$

$$k_{\mathrm{AB}}={\left(\sum _{n{n}^{\prime }l}{\bar{p}}_{n{n}^{\prime }l}^{(\mathrm{cell}\mathrm{A})}{\bar{p}}_{n{n}^{\prime }l}^{(\mathrm{cell}\mathrm{B})}\right)}^{\zeta },$$

$$P(p/p_0)=\frac{1}{\beta }\exp \left(-\frac{1}{\beta }p/p_0\right),$$

Full text

De novo exploration and self-guided learning of potential-energy surfaces

Noam Bernstein

Center for Materials Physics and Technology, U.S. Naval Research Laboratory, Washington, DC 20375, United States

  

Gábor Csányi

Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, United Kingdom

  

Volker L. Deringer

vld24@cam.ac.uk

Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, United Kingdom

Department of Chemistry, University of Cambridge, Cambridge CB2 1EW, United Kingdom

Interatomic potential models based on machine learning (ML) are rapidly developing as tools for materials simulations. However, because of their flexibility, they require large fitting databases that are normally created with substantial manual selection and tuning of reference configurations. Here, we show that ML potentials can be built in a largely automated fashion, exploring and fitting potential-energy surfaces from the beginning (de novo) within one and the same protocol. The key enabling step is the use of a configuration-averaged kernel metric that allows one to select the few most relevant structures at each step. The resulting potentials are accurate and robust for the wide range of configurations that occur during structure searching, despite only requiring a relatively small number of single-point DFT calculations on small unit cells. We apply the method to materials with diverse chemical nature and coordination environments, marking a milestone toward the more routine application of ML potentials in physics, chemistry, and materials science.

Introduction

Atomic-scale modeling has become a cornerstone of scientific research. Quantum-mechanical methods, most prominently based on density-functional theory (DFT), describe the atomistic structures and physical properties of materials with high confidence Lejaeghere et al. (2016); increasingly, they also make it possible to discover previously unknown crystal structures and synthesis targets Oganov et al. (2019). Still, quantum-mechanical materials simulations are severely limited by their high computational cost.

Machine learning (ML) has emerged as a promising approach to tackle this long-standing problem Behler and Parrinello (2007); Bartók et al. (2010); Li et al. (2015); Artrith and Urban (2016); Smith et al. (2017); Chmiela et al. (2017); Behler (2017); Huan et al. (2017); Chmiela et al. (2018); Zhang et al. (2018). ML-based interatomic potentials approximate the high-dimensional potential-energy surface (PES) by fitting to a reference database, which is usually computed at the DFT level. Once generated, ML potentials enable accurate simulations that are orders of magnitude faster than the reference method. They can solve challenging structural problems, as has been demonstrated for the atomic-scale deposition and growth of amorphous carbon films Caro et al. (2018a), for proton-transfer mechanisms Hellström et al. (2019) or dislocations in materials Fellinger et al. (2018); Maresca et al. (2018), involving thousands of atoms in the simulation. More recently, it was shown that ML potentials can be suitable tools for global structure searches targeting crystalline phases Deringer et al. (2017, 2018a, 2018b); Podryabinkin et al. (2019), clusters Ouyang et al. (2015); Tong et al. (2018); Kolsbjerg et al. (2018); Hajinazar et al. (2019), and nanostructures Eivari et al. (2017).

Assembling the reference databases to which ML potentials are fitted is currently mostly a manual and laborious process, guided by the physical problem under study. For example, hierarchical databases for transition metals have been built that start with simple unit cells and gradually add relevant defect models Szlachta et al. (2014); Dragoni et al. (2018); liquid and amorphous materials can be described by iteratively grown databases that contain relatively small-sized MD snapshots Eshet et al. (2010); Sosso et al. (2012); Deringer and Csányi (2017); Mocanu et al. (2018). A “general-purpose” Gaussian Approximation Potential (GAP) ML model for elemental silicon was recently developed Bartók et al. (2018) which can describe crystalline phases with meV-per-atom accuracy, treat defects, cracks, and surfaces Bartók et al. (2017), and generate amorphous silicon structures in excellent agreement with experiment Deringer et al. (2018c). Despite their success in achieving their stated goals, none of these potentials are expected to be even reasonable for crystal structures not included in their databases, say, phases that are stable at very high pressures.

In contrast, structure searching (that is, a global exploration of the PES) can be a suitable approach for finding structures to be included in the training databases in the first place Hajinazar et al. (2017); Deringer et al. (2018a, b); Podryabinkin et al. (2019). The principal idea to explore configuration space with preliminary ML potentials is well established: since the first high-dimensional ML potentials have been made, it was shown how they can be refined by exploring unknown structures Behler and Parrinello (2007); Behler et al. (2008); Sosso et al. (2012), and “on the fly” schemes were proposed to add required data while an MD simulation is being run Li et al. (2015); Jinnouchi et al. (2019a); Vandermause et al. (2019); Jinnouchi et al. (2019b). We have previously shown that the PES of boron can be iteratively sampled without prior knowledge of any crystal structures involved; we called the method “GAP-driven random structure searching” (GAP-RSS) Deringer et al. (2018a), reminiscent of the successful Ab Initio Random Structure Searching (AIRSS) approach Pickard and Needs (2006, 2011). Subsequently we demonstrated, by way of an example, that the crystal structure of black phosphorus can be discovered by GAP-RSS within a few iterations, and we identified several previously unknown hypothetical allotropes of phosphorus Deringer et al. (2018b).

In the context of ML potential fitting, so-called “active learning” schemes which detect extrapolation (indicating when the potential moves away from “known” configurations) are currently receiving much attention. A query-by-committee active-learning approach was suggested in 2012 by Artrith and Behler Artrith and Behler (2012). More recently, Shapeev and co-workers employed Moment Tensor Potentials Shapeev (2016) with active learningPodryabinkin and Shapeev (2017) to explore the PES and to fit ML potentials Podryabinkin et al. (2019); Gubaev et al. (2019), and E and co-workers described a generalized active-learning scheme for deep neural network potentials Zhang et al. (2018). So far, these studies mainly focused on specific intermetallic systems, namely, Al–Mg Zhang et al. (2018) and Cu–Pd, Co–Nb–V, and Al–Ni–Ti Gubaev et al. (2019), respectively. Furthermore, Podryabinkin et al. showed that their approach can identify various existing and hypothetical boron allotropes Podryabinkin et al. (2019). Finally, Jinnouchi et al. demonstrated how ab initio molecular dynamics (AIMD) simulations of specific systems can be sped up by active learning of the computed forces (in a modified GAP framework), using the predicted error of the Gaussian process to select new datapoints and to improve the speed of AIMD Jinnouchi et al. (2019a, b).

In this work, we present an efficient and unified approach for fitting ML potentials by GAP-RSS, exploring structural space from the beginning (de novo) by ML-driven searching and similarity-based distance metrics, all without any prior input of what structures are or are not relevant. We demonstrate the ability to cover a broad range of structures and chemistries, from graphite sheets to a densely packed transition metal. Our work provides conceptual insight into how computers can discover structural chemistry based on data and distance metrics alone, and it paves the way for a more routine application of ML potentials in materials discovery.

Results

A unified framework for exploring and fitting structural space

The overarching aim is to construct a ML potential with minimal effort: both in terms of computational resources and in terms of input required from the user. In regard to the former, we use only single-point DFT computations to generate the fitting database Deringer et al. (2018a). In regard to the latter, we define general heuristics wherever possible, such that neither the protocol nor its parameters need to be manually tuned for a specific system. The ML architecture is based on a hierarchical combination of two-, three-, and many-body descriptors Deringer and Csányi (2017). The two parameters that need to be set by the user are a “characteristic” distance and whether the material is primarily covalent or metallic. For the distance, we choose tabulated covalent (for C, B, and Si)Cordero et al. (2008) or metallic (for Ti) radii, depending on the nature of the system. These define the volume of the initial structures and the cutoffs for the ML descriptors (Methods section).

Our approach is based on an iterative cycle, as shown in the diagram in Fig. 1a. We generate ensembles of randomized structures as in the AIRSS framework Pickard and Needs (2006, 2011), a structure-searching approach that is widely used in physics, chemistry, and materials science Pickard and Needs (2008); Marqués et al. (2011); Stratford et al. (2017). In the first iteration, we generate 10,000 initial structures, from which we select the $N$ most diverse ones using the leverage-score CUR algorithm Mahoney and Drineas (2009). The distance between candidate structures, therein, is quantified by the Smooth Overlap of Atomic Positions (SOAP) descriptor Bartók et al. (2013), which has been widely used in GAP fitting Deringer and Csányi (2017); Bartók et al. (2018) and in structural analysis De et al. (2016); Mavračić et al. (2018); Caro et al. (2018b). While SOAP is normally used to discriminate between pairs of individual atoms and thus of local configurations, we here use a configuration-averaged SOAP descriptor that compares entire unit cells to one another (Methods section) Mavračić et al. (2018). We find this to be crucial for selecting the most representative structures, of which we can only evaluate a small number ($\ll 10,000$) with DFT. We also generate dimer configurations in vacuum at a wide range of bond lengths.

With the starting configurations in hand, we perform single-point DFT computations and fit an initial (coarse) potential to the resulting data; in subsequent iterations, we extend the database and thereby refine the potential Deringer et al. (2018a). In each iteration, we start from the same number of initial structures, and minimize their enthalpy using the GAP from the previous iteration. We then select the $N$ most relevant and diverse configurations from the full set of configurations seen throughout the minimization trajectories, for which we employ a combination of Boltzmann-probability biased flat histogram sampling (to focus on low-energy structures) and leverage-score CUR (to select the most diverse structures among those), as illustrated in Fig. 1b. These selected configurations are evaluated using single-point DFT calculations and added to the fitting database.

The iterative procedure runs until the results are satisfactory. Here we terminate our searches after 2,500 DFT data points have been collected, and our results show this to be sufficient to discover and describe all structures discussed in the present work. Other quality criteria, such as based on the distribution of energies in the database Deringer et al. (2018a), might be defined as well; the generality of our approach is not affected by this choice.

We demonstrate the method for boron, one of the most structurally complex elements Albert and Hillebrecht (2009). With the exception of a high-pressure $\alpha$-Ga type phase, all relevant boron allotropes contain B12 icosahedra as the defining structural unit Albert and Hillebrecht (2009). Boron has been the topic of structure searches with DFT Oganov et al. (2009); Wu et al. (2012); Mannix et al. (2015); Ahnert et al. (2017) and, more recently, with ML potentials for bulk allotropes Deringer et al. (2018a); Podryabinkin et al. (2019) and gas-phase clusters Tong et al. (2018). Our previous work showed how the PES for boron can be fitted in a ML framework Deringer et al. (2018a), leading to the first interatomic potential able to describe the different allotropes. However, at that time, we generated and fed back 250 cells per iteration (without further selection), and added the structure of $\alpha$-B12 manually at a later stage.Deringer et al. (2018a)

Our new protocol “discovers” the structure of $\alpha$-B12 in a self-guided way, as illustrated in Fig. 2. The increasingly accurate description of the B12 icosahedron is reflected in a gradually lowered energy error, falling below the 10 meV/atom threshold with less than 2,000 DFT evaluations, and below 4 meV/atom once the cycle is completed. This improvement is best understood by inspecting the respective lowest-energy structures that enter the database in a given iteration (Fig. 2). The lowest-energy structure at point A already contains several three-membered rings, but no B12 icosahedra yet. With one more iteration, there is a sharp drop in the GAP error (from 175 to 51 meV/at.), concomitant with the first appearance of a rather distorted $\alpha$-B12 structure (B). The final database has seen several instances of the correctly ordered structure (C).

Our experiment suggests that adding $N=100$ structures at each step slightly outperforms a similar cycle with $N=250$, although both settings lead to satisfactory results. In the remainder of this paper, we will perform all GAP-RSS searches with $N=100$ and up to a total of 2,500 single-point DFT evaluations.

Learning diverse crystal structures

Our method is not restricted to a particular chemical system. To demonstrate this, we now apply it to three prototypical materials side by side: carbon, silicon, and titanium, which all exhibit multiple crystal structures.

In carbon (Fig. 3a), both the layered structure of graphite and the tetrahedral network of diamond are correctly “learned” during our iterations. For graphite, the energy error reaches a plateau after only a few hundred DFT evaluations; for diamond, the initial error is very large, and after a dozen or so iterations we observe a rapid drop—concomitant with a drop in the error for the structurally very similar lonsdaleite (“hexagonal diamond”). The final prediction error is well below 1 meV/atom for the $sp^{3}$ bonded allotropes, and on the order of 4 meV/atom for graphite. We have previously shown that the forces in diamond show higher locality than those in graphite, making their description by a finite-ranged ML potential easier Deringer and Csányi (2017), given that sufficient training data are available. We also note that our method captures the difference between diamond and lonsdaleite very well: its value is 27 meV/atom with the final GAP-RSS version, and 28 meV/atom with DFT.

In silicon (Fig. 3b), the ground-state (diamond-type) structure is very quickly learned, more quickly so than diamond carbon, which we ascribe to the absence of a competing threefold-coordinated phase in the case of Si. We further test our evolving potentials on the high-pressure form, the $\beta$-tin type allotrope (space group $I4_{1}/amd$), which is easily discovered; the larger residual error for $\beta$-Sn-type than for diamond-type Si is consistent with previous studies using a manually tuned potential Bartók et al. (2018). We also test our method on a recently synthesized open-framework structure with 24 atoms in the unit cell (oS24) Kim et al. (2015), which consists of distorted tetrahedral building units that are linked in different ways, which the potential has not “seen”. Still, a good description is achieved after a few iterations.

In titanium (Fig. 3c), a hexagonal close packed (hcp) structure is observed at ambient conditions; however, the zero-Kelvin ground state has been under debate: depending on the DFT method, either hcp or the so-called $\omega$ phase is obtained as the minimum. Our method clearly reproduces the qualitative and quantitative difference between the two allotropes (22 meV/atom with the final GAP-RSS iteration versus 24 meV/atom with DFT) at the computational level we use.

Looking beyond the minimum structures, the DFT energy–volume curves are, by and large, well reproduced by GAP-RSS; see Fig. 3d–f. There is some deviation at large volumes for hcp and $\omega$-type Ti, but this is an acceptable issue as these regions of the PES are not as relevant, corresponding to negative external pressure. If one were interested in very accurate elastic properties, one would choose to include less dense structures by modifying the pressure parameters (Methods section, Eq. 6). Indeed, it was recently shown that a ML potential for Ti, fitted to a database of 2,700 structures built from the phases on which we test here ($\omega$, hcp, bcc) and other relevant structures can make an accurate prediction of energetic and elastic properties Takahashi et al. (2017).

Entire potential-energy landscapes

While the most relevant crystal structures for materials are usually well known and available from databases, we show that our chemically “agnostic” approach is more general. In Fig. 4, we show an energy–energy scatter plot for the last set of GAP-RSS minimizations, evaluated with DFT and with the preceding GAP version, and again across three different chemical systems. We survey both the low- and higher-energy regions of the PES—up to 1 eV per atom, which is very roughly the upper stability limit at which crystalline carbon phases may be expected to exist Aykol et al. (2018).

To analyze and understand the outcome of these searches in structural and chemical terms, we compute the distance between any two structures A and B as

[TABLE]

where $k_{\rm AB}$ is again the configuration-averaged SOAP kernel for A and B,i.e. a distance between two entire unit cells (Methods section). We then use a dimensionality reduction technique to draw a two-dimensional structural map reflecting these distances. Such SOAP-based maps have been used with success to analyze structural and chemical relationships in different materials datasets De et al. (2016); Engel et al. (2018); Caro et al. (2018b). Here, we use them to illustrate how different materials (including their allotropes as known from chemistry textbooks) are related in structural space.

To compare different materials with inherently different absolute bond lengths, we re-scale their unit cells such that the minimum bond length in each is $r_{0}=1.0$ Å, inspired by approaches for topological analyses of different structures Delgado-Friedrichs and O’Keeffe (2003). We then use SOAP ($r_{\rm cut}=2.5~{}r_{0}$, $\sigma_{\rm at}=0.1~{}r_{0}$) to determine the distance of all pairs of cells in this set, and use principal component analysis to represent this dataset in a 2D plane. Figure 5 shows the resulting plot, in which we have encoded the species by symbols and the average coordination number by color. (Coordination numbers are determined by counting nearest neighbors up to $1.2~{}r_{0}$.)

The results fall within four groups, moving from the left to the right through Fig. 5. The first group is given by graphite-like structures; they are three-fold coordinated and only carbon structures (circles) are found there. Roman numerals in Fig. 5 indicate examples, and in this first group we observe flat (i) and buckled (ii) graphite sheets. In the second group, we have four-fold coordinated (“diamond-like”) networks, made up by both carbon and silicon (recall that we are using a normalized bond length, so diamond carbon and diamond-type silicon will fall on the same position in the plot). The structures that are shown as insets are characteristic examples; from left to right, there is a distorted lonsdaleite-type structure (iii), the well-known unj framework (also referred to as the “chiral framework structure” in group-14 elements Pickard and Needs (2010); iv), and a more complex $sp^{3}$-bonded allotrope (v). While the axis values in our plot are arbitrary, they naturally reflect the structural evolution toward higher coordination numbers, and therefore we next observe a set of high-pressure silicon structures (squares), such as the simple-hexagonal one (vi), with an additional contribution from lower-coordinated titanium structures (circles). Finally, there is a set of densely packed structures, all clustered closely together; these are titanium structures including hcp (vii) and the $\omega$ type (viii). In the center of the plot, there is a structure that bears resemblance to none of the previously mentioned ones (ix), an energetically high-lying and strongly disordered intermediate from a relaxation trajectory that was added to the reference database, rather than a local minimum. This dissimilarity is reflected in relatively large distances from other entries in the SOAP-based similarity map.

Discussion

We have shown that automated protocols can be designed for generating structural databases and fitting potential-energy surfaces of materials in a self-guided way. This allows for the generation of ML-based interatomic potentials with minimal effort, both in terms of computational and user time, and it represents a step toward wide applicability of these techniques in computational materials science. Once a core (RSS-based) database has been constructed, it can be readily improved by adding defect, surface, and liquid/amorphous structural models, while at the same time being sufficiently robust to avoid unphysical behavior—even when taken to the more extreme regions of configuration space that are explored early on during RSS.

We targeted here the space of three-dimensional inorganic crystal structures, but conceptually similar approaches may be useful for nanoparticles Artrith and Kolpak (2014); Kolsbjerg et al. (2018). Finally, organic (molecular) materials are also beginning to be described very reliably with ML potentials Smith et al. (2017); Chmiela et al. (2018), and an interesting open question is how to use the structural diversity inherent in RSS in the context of organic solids Zilka et al. (2017).

Methods

Interatomic potential fitting

To fit interatomic potentials, we use the Gaussian Approximation Potential (GAP) ML framework Bartók et al. (2010) and the associated computer code, which is freely available for non-commercial research at http://www.libatoms.org. Compared to previous work, we have here developed suitable heuristics to automate and generalize the choice of fitting parameters where possible.

We use a linear combination of 2-, 3-, and many-body terms following Refs. 70 and 30, with defining parameters given in Table 1. The 2-body (“2b”) and 3b descriptors are scalar distances and symmetrized three-component vectors, respectively. For the many-body term, we use the Smooth Overlap of Atomic Positions (SOAP) kernel Bartók et al. (2013), which has been used to fit GAPs for diverse systems Szlachta et al. (2014); Deringer and Csányi (2017); Bartók et al. (2018); Mocanu et al. (2018). The overall energy scale of each descriptor’s contribution to the predicted energy (controlled by the parameter $\delta$) Bartók and Csányi (2015) is set automatically in our protocol. The 2b value is set from the variance of energies in the fitting database, the 3b value is set from the energy error between a 2b only fit and the fitting database, and the SOAP value is set from the energy error for a 2b+3b only fit.

The cutoffs for the three types of descriptors are expressed in terms of the characteristic radius $r$ (Table 1): that for 2b is longest range, while that for 3b is shortest (intended to capture only nearest neighbors), and the SOAP is intermediate in range. The resulting cutoff settings are listed in Table 1, the characteristic radii $r$ for the systems studied here being 0.84, 0.76, 1.11, and 1.47 Å for B, C, Si, and Ti, respectively.

The weights on the energies, forces, and stresses that are fit are set by diagonal noise terms in Gaussian process regression Bartók et al. (2010). We set these according to the reference energy of a given structure, to make the fit more accurate for relatively low-energy structures at each volume while providing flexibility for the higher-energy regions. The values are piecewise-linear functions in $\Delta E$, which is the per-atom reference energy difference relative to the same volume on the convex hull bounding the set of $(V,E)$ points from below (in energy). For the energy the error $\sigma_{E}$ is 1 meV/atom for $\Delta E\leq 0.1$ eV, 100 meV/atom for $\Delta E\geq 1$ eV, and linearly interpolated in between. For forces the corresponding $\sigma_{F}$ values are 31.6 and 316 meV/Å, and for virials the $\sigma_{V}$ values are 63.2 and 632 meV/atom.

Comparing structures

We also use SOAP, although with different parameters ($n_{\rm max}=l_{\rm max}=12$, $\sigma_{\rm at}=0.0875$ Å, $r_{\mathrm{cut}}=10.5$ Å), to compare the similarity of environments (as proposed in Ref. 54) in selecting from which data to train (in the CUR step). We obtain what we call a “configuration-averaged” SOAP by averaging over all atoms in the cell. In the SOAP framework Bartók et al. (2013), the neighbor density of a given atom $i$ is expanded using a local basis set of radial basis functions $g_{n}$ and spherical harmonics $Y_{lm}$,

[TABLE]

where $j$ runs over the neighbours of atom $i$ within the specified cutoff (including $i$ itself). To obtain a similarity measure between unit cells, rather than individual atoms, we then average the expansion coefficients over all atoms $a$ in the unit cell,

[TABLE]

and construct the rotationally invariant power spectrum for the entire unit cell Mavračić et al. (2018),

[TABLE]

Note that this is not equal to the average of the usual atomic SOAP power spectra used to describe the atomic neighbor environments. The final kernel to compare two cells, A and B, is then

[TABLE]

where $\zeta$ is a small integer number (here, $\zeta=4$).

Iterative generation of reference data

Randomized atomic positions are generated using the buildcell code of the AIRSS package version 0.9, available at https://www.mtg.msm.cam.ac.uk/Codes/AIRSS. The positions are repeated by 1–8 symmetry operations, and the cells contain 6–24 atoms. A minimum separation is also set, with a value of $1.8r$. The volumes per atom of the random cells are centered on $V_{0}=14.5\,r^{3}$ for covalent, and $V_{0}=5.5\,r^{3}$ for metallic systems. In the initial iteration, half of the structures are generated from the buildcell-default narrow range of volumes, and half from a wider range, $\pm 25$% from the heuristic value. In all later iterations, only the default narrow range is used. The wide volume range configurations are meant to simply span a wide range of structures Deringer et al. (2018a), and use only even numbers of atoms. The narrow volume range configurations are meant to be good initial conditions for RSS, and so for 80% (20%) of the seed structures, we choose even (odd) numbers of atoms, respectively. This is because for most known structures, the number of atoms in the conventional unit cell is even (eight for diamond and rocksalt, for example), although for some it is odd, including the $\omega$ phase Sikka et al. (1982). Biasing initial seeds toward distributions that occur in nature is a central idea within the AIRSS formalism Pickard and Needs (2011).

With the initial potential available, we then run structural optimizations by relaxing the candidate configurations with a preconditioned LBFGS algorithm Packwood et al. (2016) to minimize the enthalpy until residual forces fall below 0.01 eV/Å. As in Ref. 19, we employ a random external pressure $p$ with probability density

[TABLE]

here with $p_{0}=1$ GPa, and $\beta=0.2$. This protocol ensures that there is a small but finite external pressure, and also some smaller-volume structures are included in the fit Deringer et al. (2018a, b). We choose the same pressure range for all materials, for simplicity, although this value could be adjusted depending on the pressure region of interest.

The selection of configurations for DFT evaluation and fitting at each iteration involves a Boltzmann-biased flat histogram and leverage-score CUR, as illustrated in Fig. 1. To compute the selection probabilities for the flat histogram stage, the distribution of enthalpies (each computed using the pressure at which the corresponding RSS minimization was done) is approximated by the numpy num histogram function, with default parameters. The probability of selecting each configuration is inversely proportional to the density of the corresponding histogram bin, multiplied by a Boltzmann biasing factor. The biasing factor is exponential in the energy per atom relative to the lowest energy configuration, divided by a temperature of 0.3 eV for the first iteration, 0.2 eV for the second, and 0.1 eV for all remaining iterations. The leverage-score CUR selection is based on the singular-value decomposition of the square kernel matrix using the SOAP descriptors (with the dot-product kernel and exponentiation by $\zeta$, Eq. 5). Applying the same algorithm to the rectangular matrix of SOAP descriptor vectors was significantly less effective.

Computational details

Reference energies and forces were obtained using DFT, with projector augmented-waves (PAW) Blöchl (1994); Kresse and Joubert (1999) as implemented in the Vienna Ab Initio Simulation Package (VASP) Kresse and Furthmüller (1996). Valence electrons were described by plane-wave basis sets with cutoff energies of 500 (B), 800 (C), 400 (Si), and 285 eV (Ti), respectively. Exchange and correlation were treated using the PBEsol functional Perdew et al. (2008) for all materials except carbon, where the opt-B88-vdW functional Dion et al. (2004); Román-Pérez and Soler (2009); Klimeš et al. (2011) was chosen to properly account for the van der Waals interactions in graphitic structures. Benchmark data for energy–volume curves were obtained by scaling selected unit cells within given volume increments and optimizing while constraining the volume and symmetry of the cell.

Data access statement

Data supporting this publication will be made openly available via https://www.repository.cam.ac.uk.

Author contributions

N.B., G.C., and V.L.D. jointly designed the research, developed the approach, and analyzed the data. N.B. developed the computational framework and performed the computations with input from all authors. V.L.D. wrote the paper with input from all authors.

Acknowledgements.

We thank Profs. C. J. Pickard and D. M. Proserpio for ongoing valuable discussions. N.B. acknowledges support from the Office of Naval Research through the U. S. Naval Research Laboratory’s core basic research program. G.C. acknowledges EPSRC grant EP/P022596/1. V.L.D. acknowledges a Leverhulme Early Career Fellowship and support from the Isaac Newton Trust.