Global atomic structure optimization through machine-learning-enabled barrier circumvention in extra dimensions

In the following, we explain the applied methodology in more detail. We start with the DFT calculations and follow this with the details of the machine-learning model, i.e. the fingerprint and the Gaussian process including optimization of hyperparameters. We then provide a description of the Bayesian search algorithm including the random searches with the GP potential. Finally, we discuss success curves and the identification of ground state structures.

Electronic structure calculation

All DFT calculations are performed using GPAW^82,83,84 and the Atomic Simulation Environment (ASE)^85,86. We apply the Perdew-Burke-Ernzerhof (PBE)⁸⁷ exchange-correlation functional, a planewave cutoff of 400 eV, and a Fermi temperature of 0.1eV. For clusters, only the Γ-point is used for k-point sampling while a k-point density of 6 Å was used for periodic systems in all periodic directions and a single k-point in the non-periodic directions. For the illustrative examples in Fig. 3, effective medium theory (EMT)^53,54 was used instead of DFT. All calculations are without spin-polarization except the ones in Fig. 16. We note that the approach presented in this paper is not dependent on any specific electronic structure method or exchange-correlation functional and that any method of calculating energies and forces could have been used instead.

Machine-learning model: fingerprint

The atomic structures are represented by a fingerprint ρ(x, Q) with x and Q being the full set of spatial and elemental coordinates respectively. ρ(x, Q) consists of a radial part, ρ^R(r; x, Q), and an angular part, ρ^α(θ; x, Q). The radial fingerprint is a function of the pair-wise interatomic distances, r_ij, of atoms i and j whereas the angular fingerprint is a function of the triplet-wise angles, θ_ijk, spanned between the distance vectors r_ij and r_ik from atom i to j and from atom i to k, respectively. Both r_ij and θ_ijk, given by Eqs. ((7) and (8)), are trivially extensible to more than three dimensions. The elemental and existence fractionalization of atoms is made possible by introducing the scalar values q ∈ [0, 1] for each atom in each term of the fingerprint sum as described in refs. ^43,44. ρ^R(r; x, Q) and ρ^α(θ; x, Q) are composed of several subfingerprints for each combination of two and three elements respectively concatenated together, the formula of which are given by Eqs. ((9), (10), (11)).

where \({R}_{c}^{R}\) and \({R}_{c}^{\alpha }\) are radial and angular cutoff radii, while δ_R = 0.4Å, δ_α = 0.4Å, and γ = 2 are constants. In general we use \({R}_{c}^{R}=5{r}_{co{v}_{max}}\) and \({R}_{c}^{\alpha }=3{r}_{co{v}_{max}}\), where \({r}_{co{v}_{max}}\) refers to the covalent radius of of the largest element in the system. For systems where the radius of the smallest element is 2/3 or less than that of the largest element \({R}_{c}^{\alpha }=2.5{r}_{co{v}_{max}}\) is used instead. The subscripts A, B, and C refer to elements with each radial and angular sub-fingerprint consisting of 200 and 100 entries for each elemental combination, respectively. A radial fingerprint containing two elements would thus have sub-fingerprints ρ_AA, ρ_AB, ρ_BA and ρ_BB (with the identity ρ_AB = ρ_BA) and thus a total length of 800 entries. A similar argument can be made for the angular part resulting in eight angular sub-fingerprints. In general, the radial and angular fingerprint will contain n² and n³ sub-fingerprints respectively where n is the number of elements in the system.

The fingerprint counts all pairs and triplets within the radial and angular cutoff radii of each atom. The formalism extends to periodic boundary conditions by counting all pairs and triplets of the atoms in the primary unit cell with the atoms in a set of adjacent copies of the unit cell for any of the three standard spatial dimensions. Any hyperspatial dimension is considered non-periodic just as one would consider the third dimension non-periodic in relation to two-dimensional materials.

We now show the “free-flow” property mentioned in Section “Fingerprint” We consider a situation where two atoms, say atoms 1 and 2, exchange chemical identity. We focus on the radial part of the fingerprint as written compactly in Eq. (5), but it holds more generally. We now write out explicitly the terms involving atoms 1 and 2:

We consider a situation where no atoms are moved in coordinate space, but where the chemical identity A is transferred from atom 1 to atom 2 by an amount Δq. We have the changes Δq_1,A = − Δq, Δq_1,B = Δq, Δq_2,A = Δq, and Δq_2,B = − Δq. We furthermore assume that the distance between atoms 1 and 2 is larger than the cutoff distance so that f_c(r₁₂) = 0. In that case, the first term in Eq. (12) vanishes, and the last (fourth) term is unchanged by the process. In the remaining two terms the values of q_i,A and q_i,B do not change, so we can write the change in the fingerprint as

We now see that if the environment of atoms 1 and 2 are identical, the last parenthesis vanishes, and the fingerprint is completely unchanged during the process. We also see, that if the environments are different, the change in the fingerprint is linear in Δq, which invites a smooth variation of the energy in the Gaussian process.

Machine-learning model: Gaussian process

Energy and forces μ = (E, − F) and their associated uncertainties Σ(x, Q) are predicted by a Gaussian process described by the following equations^88,89:

where μ_p(x, Q) is the prior mean, ρ(x, Q) is the fingerprint of the predicted structure, K and C = K + χ²I are the unregularized and regularized covariance matrices, respectively, with χ being a noise parameter, P is a vector of all fingerprints in the training data, y is the training energy and force targets, and μ_p(X) is the prior mean applied to all atomic structures in the training set. \(\tilde{K}(\rho [{\bf{x}},Q],\rho [{\bf{x}},Q])\) represents the covariance of the fingerprint with itself.

In this work, the kernel function in the covariance matrix has the form of the squared exponential function:

where ∣ρ₁ − ρ₂∣ is the Euclidean distance between two fingerprint vectors, l is the length scale, and σ² is the prefactor.

Machine-learning model: prior potential function

The prior is set to a constant, μ_c, plus a repulsive potential, U[x_ij(x, Q)], depending on the spatial and elemental coordinates as described by Eq. (17).

where q_i is the existence of atom i, \({r}_{co{v}_{e}}\) is the covalent radius of element e, r_min is the radius an atom will have at no existence and f is a scaling constant set to 0.8. r_min is set to the smallest covalent radius of any atom in the system. For the prior potential U[x_ij(x, Q)], we use a repulsive potential modified to go to zero at x_ij = 1 given by Eq. (17):

where \({\sigma }_{{p}_{rep}}\) is a strength constant set to 10 eV.

The associated forces, F(x, Q), element coordinate derivatives, dq_i,e, and stresses, S(x, Q), of Eq. 20 are given by:

where we in the equation for the stress made use of the virial theorem.

On top of any prior potential, extra potentials may be applied. Excessively large cell volumes were penalized by an extra potential:

where σ_V is a strength constant and V_high is a potential onset below which the potential is zero. We set σ_V to \(10eV/{V}_{0}^{2}\) and V_high = 3.5V₀ with \({V}_{0}={\sum }_{i}\frac{4}{3}\pi {r}_{cov,i}^{3}\).

Equation. (27) describes another extra potential punishing atoms being far into a non-periodic dimension of index d with coordinates x_d

where σ_NP is a strength constant set to 10 eV/Ų. \({{\bf{x}}}_{{d}_{high}}\) and \({{\bf{x}}}_{{d}_{low}}\) are system specific potential onset values between which the potential is zero. This potential was applied to the bulk systems in Fig. 5f and Fig. 14 with \({{\bf{x}}}_{{d}_{low}}\) and \({{\bf{x}}}_{{d}_{high}}\) set to 0 and \(3{r}_{co{v}_{max}}\) respectively, where \(3{r}_{co{v}_{max}}\) is the largest covalent radius in the systems.

Machine-learning model: force and stress predictions

According to Eq. (14), the predicted force on atom i is given by

where \({{\bf{F}}}_{i}^{(p)}\) is the prior force and the kernel function, k(ρ₁, ρ₂), is taken between two atomic structures with fingerprints ρ₁ and ρ₂. Here, atomic coordinates with indices i and j contribute in ρ₁ and ρ₂, respectively, and j runs over all atoms.

Similarly the element coordinate derivative for element e of atom i is given by:

where \(d{q}_{i,e}^{(p)}\) is the prior derivatives.

As the total energy is in the end described through a fingerprint, which has an explicit dependence on the interatomic vectors r_ij, the stress can be calculated using the virial theorem. The stress is given by

where ε denotes the strain, and S^(p) is the stress from the prior. The i-sum runs over the unit cell, while the j-sum runs over the surroundings within the interaction sphere defined by the cutoff of the fingerprint.

Machine-learning model: hyperparameter optimization

For each cycle in the global optimization algorithm (Fig. 2), the hyperparameters constituted by the length scale, l, the square root of the prefactor, σ, the noise, χ, and the prior mean constant, μ_c, are updated by maximizing the a posteriori probability p(l, σ, χ, μ_c∣y), given the training data y. The noise, prior mean constant, and the prefactor are set analytically as

with \({C}_{0}(P,P)={K}_{0}(P,P)+{\chi }_{r}^{2}I\), where y_eng,n is the energy of structure n in the database containing a total of N_DFT structures, Y is the total number of training targets, K₀(P, P) is the covariance matrix without the prefactor, and χ_r is a relative noise constant set to 0.001. The relative-noise is identical for energy and force contributions. The total number of training targets is equal to 1 energy and 3N_atoms forces for each structure in the training set, i.e., Y = N_DFT × (1 + 3N_atoms). As K₀ and hence C₀ depend on the length scale, the prefactor is always evaluated with respect to a given length scale, optimized by maximizing the log posterior \(\ln [p(l| y)]\):

where \(\ln [p(y| l)]\) is the log-likelihood and p(l) is a prior distribution for the length scale. The log-likelihood is expressed as

where we recognize the last term as the optimal prefactor at a given length scale from Eq. (34). The length scale is calculated in parallel by a nested grid search in the interval \([{\rm{median}}(\Delta {\rho }_{nn}),10\max (\Delta \rho )]\) in logarithmic space where Δρ marks the set of all euclidian distances between any two fingerprints in the training set, and Δρ_nn marks the set of nearest neighbor distances i.e. the shortest distances between a given fingerprint and all other fingerprints. This interval is chosen to seek a good compromise between accuracy and interpolatability between data points and new structures in the surrogate surface with the latter being of high importance when interpolating to fictive dimensions which can not be sampled by the model. To mitigate overfitting and secure interpolatability at low datasets a log-normal length scale prior distribution is applied in Eq. (36):

where μ_LN and σ_LN are the mean and the width in the logarithmic space, respectively. σ_LN is set to 2 and μ_LN is set from the equation: \({\rm{mode}}[p(l)]=\exp ({\mu }_{LN}-{\sigma }_{LN}^{2})=0.5[{\rm{mean}}(\Delta \rho )+\max (\Delta \rho )]\).

To make sure the radial and angular fingerprint had a reliable relative scaling across systems, the angular part was scaled by the following factor:

where ∣ρ^R∣_abs and ∣ρ^α∣_abs refer to the set of absolute differences between any two fingerprints in the training set for the radial and angular fingerprints respectively.

Bayesian search algorithm: overview

The overall structure of the Bayesian search algorithm shown in Fig. 2 has already been discussed in Section “Bayesian search algorithm”, but a number of details remain to be described. The following sections describe the generation of random structures, performing relaxations with the GP surrogate potential, selection of promising structures for database inclusion using an acquisition function, and discarding of undesired structures.

A potential issue with the Bayesian search method is that the surrogate PES could “degenerate” so that the global minimum never appears and all searches would lead to local minima, but not the global one. This behavior is counteracted by several means. Firstly, the use of an acquisition function instead of the bare energy will invite for “exploration” instead of only “exploitation” of previously investigated basins of the PES. Secondly, new suggested candidate structures obtained by relaxations on the surrogate PES are not selected if they are too close to already evaluated structures in the DFT database as described below. Thirdly, it might happen that at some stage in the optimization all relaxations in the surrogate surface lead to already known configurations (or gets discarded otherwise). In that case a new, truly random structure is created, directly evaluated with DFT (without relaxation on the surrogate PES), and included in the DFT database. So in principle, there is always a completely random element ensuring that the surrogate model will be improved in new regions of the configuration space.

Bayesian search algorithm: random structure generation

The following describes different ways of randomly placing atoms in a confined space. The atoms are afterwards repelled from one another. We found the potential Eq. (20) to be too strong and instead use a softer parabolic potential given by Eq. (40):

where the strength constant \({\sigma }_{{p}_{P}}\) is set to 10eV. The scaling constant f in Eq. (19) is set to 0.9 for structures entering the initial database and for generating random structures for surrogate relaxation.

Random cells are generated by generating a unit cube as represented by a 3 × 3 unit matrix and adding random numbers in the interval [− ξ_c, ξ_c] to all entries with ξ_c = 0.25 to secure an ensemble of cells with varying yet not extreme angles between the lattice vectors. The cell is next scaled to a volume in the range [1V_base, 3V_base] while maintaining the cell morphology, where V_base is a reference volume given by:

where V_D(r_cov,i) is the volume of atom i with covalent radius r_cov in D dimensions, Γ is the gamma distribution, and D_Hyper,k is the size of the non-periodic hyperspatial dimension k, set to \(3{r}_{co{v}_{max}}\) with \({r}_{co{v}_{max}}\) being the largest covalent radius of any atom.

This procedure is chosen to secure a similar span of initial atomic packing fractions in atomic systems of different dimensionality.

While working well for most compact materials this strategy is ill-suited for bulk systems with a lot of internal vacuum in which case one would have to come up with a larger guess for V_base and possibly a larger interval range.

The atoms are subsequently placed randomly inside the cell and the hyperspatial dimensions. The structure is relaxed by the repulsive potential of Eq. (40) thus potentially slightly expanding the cell. For the case of Fig. 5f and Fig. 14, Eq. (27) was also applied alongside Eq. (40).

For clusters, i.e. non-periodic systems, a cubic cell of length 25Å was set with a centrally centered cubic subvolume box of range [1V_box, 3V_box] with V_box = ∑_iV_D(r_cov,i) within which the atoms are placed and subsequently relaxed in the repulsive potential.

For dual-atom catalysts, initial structures were generated by creating a graphene layer, randomly substitute N_N carbon atoms with nitrogen and remove N_V carbon atoms. Adsorbate atoms were randomly placed above the substrate within 3Å. In surrogate relaxations, the graphene substrate remained intact, with nitrogen substitutions and vacancies generated via ICE and ghost methods, respectively.

Random elemental coordinates are generated by combinatorial use of the Dirichlet rescale algorithm^90,91 to satisfy the elemental constraints of Eq. (1).

Bayesian search algorithm: details of the surrogate relaxations

A key element in the procedure is the relaxation of randomly generated structures in the GP surrogate potential. Figure 4 illustrates such a relaxation process for a Cu₁₈Ni₅ cluster in the GP predicted potential energy surface. Due to the additional hyperspace and elemental coordinates, it is necessary to divide the relaxation process in four phases as we shall now discuss.

In the first phase, spatial and elemental atomic coordinates are updated simultaneously. As atoms embedded in the (3 + D_Hyper)-dimensional space will not spontaneously settle into the three dimensional space, all D_Hyper coordinates are punished by a potential, U_hs and its resulting force F_hs given by Eqs. (43) and (44)

where ∣x_hs,i∣ is the Euclidean norm of the vector of hyperspatial coordinates for atom i and ω(c) is a custom time-dependent strength factor. The relaxation is structured into \(n_{c}^{hs}\) cycles of index c each lasting \(n_{c_{sub}}^{hs}\) steps. In this paper, the strength factor was set to

with the parameters a and b tuned such that k(0) = 0.1 and \(k(n_{c}^{hs})\) = 1000 with \(n_{c}^{hs}\) = 100 to set aside 25 cycles for one order of magnitude, as what constituted a good magnitude and rate progression was observed to be system specific. Too slow progressions result in long run times whereas too fast disrupts hyperspatial relaxation. Likewise, insufficient final magnitude results in the failure of squeezing the atoms into three dimensions.

During this phase, the total existence of atoms in ghost-possessing elements/ICE-groups are restricted to the interval [q_low, 1] with 1 ≫ q_low > 0 since atoms with zero existence do not interact with other atoms at all. Hence, they become idle during the relaxation as argued in ref. ⁴⁴. Consequentially, the total elemental sum of any ghost-possessing element is temporarily set to \({N}_{e}+{N}_{{e}_{Ghost}}{q}_{low}\). The phase ends by projecting all atoms from (3 + D_Hyper) dimensions into three dimensions, which happens when ∣x_hs,i∣ < 0.01Å for all atoms.

In the second phase, the relaxation proceeds as in phase 1 but with all atoms embedded in three dimensions and the existence interval of ghost-possessing elements/ICE-groups kept at [q_low, 1] for \(n_{q_{low,1}}^{3D}\) steps.

In the third phase, the existence interval of atoms belonging to ghost-possessing elements/ICE-groups is changed from [q_low, 1] to [0, 1] and the total elemental sum of atoms belonging to ghost-possessing elements is changed from \({N}_{e}+{N}_{{e}_{Ghost}}{q}_{low}\) back to N_e by removing \({N}_{{e}_{Ghost}}{q}_{low}\) of elemental existence starting from the atoms with lowest existence and up. The spatial coordinates and elemental coordinates are then optimized with \(n_{q_{0,1}}^{3D}\) steps.

In the fourth phase, the atoms of any ICE-group are assigned to an element based on the highest atomic elemental coordinate subject to the elemental sum constraints and excess atoms of any ghost-possessing elements are deleted in order of lowest to highest atomic existence. The spatial coordinates are then relaxed with all elemental coordinates kept at unit identity for \(n_{ui}^{3D}\) steps.

In all steps, the lattice vectors of the unit cell may be optimized at the same time if desired. The relaxation steps terminate either when the total number of steps is reached or when the desired convergence criteria is met.

In this work, all relaxations were limited to a maximum of 700 steps with parameters as listed in Table 1. Figure 4 illustrates a surrogate relaxation of a Cu₁₈Ni₅ cluster extended to four spatial dimensions where Cu and Ni form an ICE-group possessing 11 ghost atoms. The atoms in Fig. 4a are seen to initially form a dense globule with seemingly overlapping atoms which, when comparing to Fig. 4b, is observed to be due to atoms being distant in the fourth dimension. As the relaxation progresses, the atoms are squeezed out of the fourth dimension, and the fractional elements are generally observed to converge to 0 or 1 except for a few atoms. The atomic energy of Fig. 4a can be divided into four segments: 1) initial decline due to relaxation in the four dimensions with low penalty constant, 2) a steady increase due to the increasing penalty constant, 3) a second decrease due to atoms being squeezed out of the fourth dimension hence eliminating the penalty due to Eq. (43), and 4) a final segment with no atomic penalty where the existence fractions are also allowed to go to 0. The jagged shape of the energy curve reflects the cycles and sub-steps of the hyperdimensional squeezing phase. Which atoms should exist or not is observed to be decided during the first few steps of the relaxation as seen from Fig. 4h.

Relaxations are generally performed using the SLSQP (Sequential Least Squares Programming), except in figures with only hyperspatial optimization without elemental coordinates in which cases the L-BFGS-B method (Limited-memory Broyden-Fletcher-Goldfarb-Shanno with Bounds) is used instead as it is in general more stable than the SLSQP method. Both methods are used as implemented in the scipy package⁹². The L-BFGS-B optimizer converged when all projected gradients were below 0.01 eV/Å, and SLSQP when the energy change between iterations was under 0.001 eV.

We use 40 parallel surrogate relaxations before we apply the acquisition function and perform a DFT calculation for the best candidate.

Bayesian search algorithm: acquisition function

Selection of the best candidate structure at the end of a cycle in the global optimization algorithm is determined by an acquisition function A(x) which in the present study is set to a lower confidence bound (LCB)

where κ is a constant set to 2 while E(x) and Σ(x) are the predicted energy and uncertainty of Eqs. (14) and (15), respectively. The dependency on Q is omitted as the acquisition function is only used on atoms with unit elemental identities.

Bayesian search algorithm: discarding structures

Some structures are discarded before being evaluated as candidate structures for the database. A bulk structure is discarded if the volume of the unit cell is outside the range of 0.5 to 5 times the sum of atomic volumes. A structure is discarded if any atom is closer than 0.5 times its covalent distance to another atom. Finally, a structure is discarded if the norm of the absolute difference in fingerprints between the relaxed structure and any structure in the database is smaller than 1, to avoid training on the same structure twice. This is a fast way to compare structures also for large datasets. In the case of Cu₁₂Ni₁₁, structures were discarded if any atom was over 1.25 times the sum of its and another atom’s covalent radii apart, indicating disconnection as such examples led to DFT convergence errors. In hyperspace runs, structures were discarded if one or more atoms failed to exit the hyperspatial dimensions.

To prevent memory issues, a surrogate relaxation is prematurely terminated and the structure discarded if any of the following happens: 1) An atom exits non-periodic unit cell boundaries, 2) a unit cell length is smaller than \(2{r}_{co{v}_{max}}\) or larger than 50 Å or 3) the unit cell volume falls below 0.3 of the total atomic volume.

Success curves

A success curve illustrates the cumulative fraction of optimization runs that have successfully identified the ground state structure as a function of the number of DFT calculations performed. Each success curve is based on 20 independent global optimization runs, each limited to 100 DFT calculations.

Success is declared in the success curves for clusters when a found energy is within a margin of 0.05 eV of the lowest found energy. This rather high value is chosen to ensure that different systems in the same DFT potential basin is in fact identified as being identical if the interatomic distances are a bit off. For Cu₁₂Ni₁₁ in Fig. 12 the energy margin is set to 0.01 eV to only find the lowest energy structure. The found structures are checked for structural agreement by eye inspection. Success for bulk materials is declared when a structure is found to be structurally equivalent to the lowest energy structure by the pymatgen package⁹³.

To estimate the uncertainty, we represent a success curve as n+m independent attempts to find the global minimum structure, where n and m indicate the number of successful and unsuccessful attempts respectively⁴⁴. By applying Bayes’ theorem with a uniform prior, the posterior probability of success p_s follows a Beta distribution B(p_s∣α = n + 1, β = m + 1). We take the mode of this distribution, given by mode(p_s) = n/(n + m), as the value of the success curve, and take the square root of the variance to express the uncertainty:

In the regions where the success curves are either zero (0% success) or one (100% success), the uncertainty is set to zero.