Global analysis of energy landscapes for materials modeling: A test case for C₆₀

The motivation to develop empirical interatomic and intermolecular potentials to provide more efficient computer simulations of molecules and materials has a long history. The earliest work generally considered pairwise isotropic representations, such as the Lennard–Jones¹ and Morse potentials,² with functional forms inspired by quantum mechanics. The theory of interatomic and intermolecular forces^3,4 developed in parallel with efforts to provide increasingly accurate descriptions of systems ranging from van der Waals complexes to condensed phases. In turn, these potentials were employed in numerical simulations designed to predict and understand key thermodynamic and dynamical properties, and they continue to underpin major research efforts in chemical physics, biophysics, and materials science.

Recent developments have exploited the ability to fit flexible functional forms to large and very diverse structural databases using machine learning methodologies.^5–8 This work has made a particular impact in materials modeling, where new developments arrive at a remarkable pace.⁹ The much lower computational cost of fitted analytical potentials, compared to density functional theory, has facilitated studies of larger and more complex systems, including reactive processes.¹⁰ 

Our aim in the present work is to present a detailed account of how global exploration, using computational tools from the energy landscape approach, can reveal weaknesses in a fitted potential and provide a systematic and efficient means to remove them. Such problems might not be immediately evident in conventional simulations because unphysical parts of the potential can lie behind relatively large barriers. An enhanced sampling methodology is then required. The usual wisdom for empirical potentials is that they are only accurate for the few key observables they have been fitted to, but at the same time, they should not be unphysical anywhere. Machine learned potentials that are not based on physically motivated functional forms can be very accurate for a wide range of configurations that lie not too far from the fitting domain. However, there is no general theory to predict what “too far” might be. The issue of extrapolation is likely to be even more important, and building in the correct limits for the potential by including appropriate configurations in the fitting database is likely to extend the domain of useful applicability.

We focus here on the global analysis of a particular set of Gaussian Approximation Potentials (GAPs)^6,11 for carbon and their ability to reproduce the energy landscape for C₆₀. The fundamental importance of carbon, including the wide variety of finite structures that elemental carbon can adopt, along with a complex phase diagram, makes this system a key target. In a previous report, we explored the C₆₀ landscape beyond the fullerene domain to predict pathways between buckybowls and buckminsterfullerene.¹² We harvested kinetic transition networks^13–15 using two efficient tight-binding potentials, namely GFN2-xTB¹⁶ and SCC-DFTB.¹⁷ Interfaces to these potentials were added to the Cambridge software for landscape exploration, specifically GMIN¹⁸ for basin-hopping global optimization,^19–21, OPTIM²² for locating transition states and characterizing the associated pathways, and PATHSAMPLE,²³ which is a driver for OPTIM and also features a wide variety of methods for calculating emergent thermodynamic and kinetic properties from a stationary point database. In the present work, new interfaces were constructed for the QUIP package,²⁴ which provides access to the Gaussian approximation potentials from GMIN and OPTIM. These programs are all available for download under the GNU General Public License. The first GAP model that was published for carbon focused on the most challenging amorphous phase²⁵ and was later significantly extended to describe many other phases and nanostructures,^26,27 fitting over 6000 configurations. It is this latter model that is the focus of our study here. The corresponding kernel is a smooth, analytic function of the atomic positions, including the analytic cutoff function, so that gradients and second derivatives are also analytic.

We found that the published carbon GAP potential^26,27 reproduced the fullerene landscape leading to buckminsterfullerene quite well. However, further exploration using discrete path sampling^28,29 revealed two additional families of C₁@C₅₉ and C₂@C₅₈ structures with lower energies than buckminsterfullerene. Although the original potential was fitted to a database of configurations that was iteratively constructed using random structure search³⁰ and high temperature molecular dynamics driven by the potential itself, clearly there are regions of configuration space that are poorly described but are separated from the usual structures by high barriers and are, thus, not normally explored. Seven refined GAP potentials were then considered, and the key pathways were recomputed for each of them. As detailed below, we find that a fit including just the lowest energy C₂@C₅₈ minimum and the corresponding minimum after relaxation with density functional theory yields new potentials where both the C₁@C₅₉ and C₂@C₅₈ families are shifted to much higher and more realistic energies.

The computational energy landscape approach involves geometry optimization to locate local minima and the transition states and pathways that connect them; these methods have been described in detail and reviewed elsewhere.^15,31–33 Post-processing involves the standard methodologies of statistical mechanics and unimolecular rate theory.^34,35 A brief summary is provided here.

We employ the geometrical definition of a transition state as a stationary point (vanishing gradient) that has one negative Hessian eigenvalue.³⁶ For local minimization in GMIN and OPTIM, we use a customized L-BFGS³⁷ (limited memory Broyden,³⁸ Fletcher,³⁹ Goldfarb,⁴⁰ and Shanno⁴¹) algorithm.⁴² The same method is employed to characterize pathways and minimize a chain of images in the doubly-nudged⁴³ elastic band^44,45 (DNEB) approach and the Rayleigh–Ritz and tangent space minimizations in hybrid eigenvector-following.^46–49 DNEB is used for double-ended connection attempts between selected pairs of minima, with images corresponding to local maxima taken as transition state candidates for accurate refinement using hybrid eigenvector-following. In multistep pathways where a complete connection requires more than one cycle, we use the missing connection algorithm⁵⁰ to choose new pairs of minima.

In the first phase of this study, we reoptimized all the transition states obtained with the two tight-binding representations in previous work¹² using hybrid eigenvector-following.^46–49 These transition states and the corresponding connected minima were used to construct an initial stationary point database, which was then expanded using PATHSAMPLE to manage OPTIM jobs running double-ended searches. We also attempted to connect minima to the principally connected component of the network, containing buckminsterfullerene, if they were initially disconnected. These schemes correspond to the discrete path sampling^28,29 approach for generating a kinetic transition network via geometry optimization.

The global optimization searches conducted for the original landscapes all employed basin-hopping.^19–21 The lowest minimum for the tight-binding potentials in the absence of periodic boundary conditions was buckminsterfullerene, as expected.

As the baseline model, we use the Gaussian approximation potential for carbon introduced in Ref. 26 and later corrected to use consistent k-point sampling.²⁷ The unique identifier of this potential is GAP_2021_4_21_0_21_17_24_635. The detailed construction of GAP models has been described in depth in Ref. 8, and the one we use is a sum of a 1/r⁶ long range pair potential, which represents dispersion, along with two-body, three-body, and many-body short range terms, the latter utilizing the SOAP descriptor.¹¹ The potential has been fitted to density functional theory (DFT) data computed using the Vienna ab initio simulation package (VASP). The exchange-correlation functional was optB88-vdW.⁵¹ 

We explore several enhancements to the above baseline potential, and here we describe some details of these modifications. The original training dataset contained carbon dimers down to an interatomic distance of 1.12 Å, which was sufficient to yield a potential that had a core repulsion strong enough so that atoms did not unphysically collide during the high temperature molecular dynamics that was used in the original publication. However, the global optimization tools used in the present work require the core repulsion to extend to smaller interatomic distances. Therefore, we calculated the dimer curve of carbon with the reference DFT method down to a distance of 0.18 Å and included a spline fit to this data in the pair potential that forms part of the GAP model. The models that include this modification have “repcore” in their labels. A different kind of enhancement is the utilization of improved radial basis functions, in particular the “turbo soap” construction introduced in Ref. 52, and we identify such models with the label “turbo.” This additional short-range repulsion has little effect on the physically relevant pathways since they do not sample such distances. It basically serves to ensure that unphysical distances have appropriate high energies and are avoided in the “repcore” potentials.

The clusters were placed in a periodically repeated cubic cell with a side length of 25 Å, which is large enough so that the periodic images do not interact. The initial network for the first GAP potential considered⁵⁶ contained 54 772 minima and 46 777 transition states after relaxation from the previous tight-binding results¹² and subsequent additional sampling. Two additional families of structures emerged during subsequent runs to explore the landscape more thoroughly, both containing minima below buckminsterfullerene.

The C₆₀ landscape for the low-lying fullerenes^57–59 reproduces structures and energetics quite well for local minima in the stacks⁶⁰ defined by successive pyracylene rearrangements.⁶¹ However, the global minimum is a C₂@C₅₈ structure, and there are also numerous low energy minima with C₁@C₅₉ morphology, where one surface atom adopts a bridged, low coordination configuration. There are many low energy C₂@C₅₈ minima with alternative orientations of the interstitial C₂ moiety and different C₅₈ shells. There are also various single shell minima with heptagonal and square faces that should lie at higher energy. The global organization of the landscape is illustrated in Fig. 1. The magnification of the low-lying region shown in Fig. 2 highlights the lowest fullerene structures and selected low energy minima that do not correspond to a C₆₀ fullerene shell. These minima are artifacts of the original GAP potential, and we discuss our strategy to deal with them in Sec. IV. The energetics are summarized in Table I. DFT results consistently place the C₂@C₅₈ minimum about 10 eV above buckminsterfullerene, whereas the original GAP model produces a C₂@C₅₈ global minimum about 2 eV lower in energy.

To analyze the GAP results in more detail, we examined three key pathways to buckminsterfullerene. These pathways were extracted from the network using Dijkstra’s shortest path algorithm⁶² with edge weights based on branching probabilities calculated using transition state theory rate constants between connected minima^34,35 obtained for harmonic vibrational densities of states. This approach finds the path with the largest contribution to the overall rate constant when intervening minima are placed in a steady state.^28,63

The first pathway considered was the selected bowl structure used in analogous calculations for GFN2-xTB and SCC-DFTB in previous work.¹² The energetic profile is shown in Fig. 3 and is comparable to the results obtained for the tight binding potentials. The other two pathways we examined link the lowest C₂@C₅₈ and C₁@C₅₉ minima to buckminsterfullerene. The fastest path from the C₂@C₅₈ global minimum involves three steps, where the cage first rearranges to make a heptagonal ring, and C₂ then reorients and inserts in the cage (Fig. 4). The pathway from the lowest C₁@C₅₉ minimum to buckminsterfullerene has five steps (Fig. 5). The endohedral carbon first inserts into the cage to give the second lowest C₆₀ fullerene, which rearranges to the lowest energy icosahedral fullerene in the last step via the pyracylene C₂ local rotation mechanism.⁶¹ 

Having identified artifacts in the original potential, our objective is to improve the fitting with feedback from surveys of the landscape until the minima in question are removed or shifted to appropriate energies. In particular, we chose to augment the training set with the lowest energy C₂@C₅₈ minimum and the corresponding minimum obtained from it by relaxation using density functional theory. This relaxation also employed the OPTIM program, with energies and gradients supplied by CASTEP⁶⁴ for the PBE⁶⁵ functional. The two additional C₂@C₅₈ structures are identified in Table I. Altogether, seven new GAP models were considered, as follows:

1. soap_2023_jan (original parameters, mpi parallelization)

2. soap_turbo (original training set, improved radial basis functions)

3. soap_with_extra_configs (two new configurations in the training set)

4. soap_turbo_with_extra_configs (two new configurations in the training set and an improved radial basis)

5. soap_turbo_repcore (original training set, improved radial basis functions, enhanced repulsive core)

6. soap_with_extra_configs_repcore (two new configurations in the training set, enhanced repulsive core)

7. soap_turbo_with_extra_configs_repcore (two new configurations in the training set and an improved radial basis, enhanced repulsive core)

The first of this set gives very similar results to the original potential, as expected, and was included to validate the GAP code that has significantly evolved since the original publication of the baseline model. The other models involve different combinations of fits that include the two C₂@C₅₈ structures obtained with the original potential, an enhanced repulsive core with dimer data extended down to 0.18 Å (see Sec. II), and improved radial basis functions. In each case, we first relaxed the transition states and minima for each of the three paths to buckminsterfullerene considered in Sec. III and illustrated in Figs. 3–5. However, potentials (2) and (4) with improved radial basis functions exhibit unphysical behavior for most of the stationary points in the three-step C₂@C₅₈ path (buckminsterfullerene itself is not a problem) and for the two highest energy minima in the bowl to buckminsterfullerene pathway. These cases were not investigated further but instead led us to develop potentials (5), (6), and (7), which include an additional repulsive core term.

Some of the resulting landscapes did not have complete connections between the chosen end minima, and these missing pathways were completed using PATHSAMPLE to organize appropriate double-ended searches with OPTIM using the missing connection algorithm.⁵⁰ For the bowl to buckminsterfullerene path, all the low energy steps are essentially the same for potentials (1) and (3) (Fig. 6). The overall paths for these cases are shorter than for the original potential, but it is not surprising to find different routes for pathways with over 50 steps. The path for potential (6), where an enhanced dimer core repulsion was added, is similar to the result for the original potential. However, the pathways for potentials (5) and (7), with the extended radial basis set, are significantly longer.

In the C₂@C₅₈ to buckminsterfullerene path, the original and soap_2023_jan paths are very similar. For potential (3), trained on the two extra configurations, the pathway changes significantly, and the C₂@C₅₈ minimum is now much higher in energy (Fig. 7). The C₂@C₅₈ minimum is also higher in energy for all three potentials with the extended repulsive core, even for (5), where the additional configurations are not included in the fit. As for the bowl to buckminsterfullerene paths, potentials (5) and (7) with the extended radial basis correspond to much longer paths.

The results for the C₁@C₅₉ to buckminsterfullerene path are similar: the C₁@C₅₉ decoy minimum is now a high energy structure for new potentials (3) and (4) trained on the decoy structures (Fig. 8). In this case, the new parameterizations (2) and (4) with improved radial basis sets did not hit any discontinuities. Both of these potentials, including the one with only the original training data, shift the C₁@C₅₉ minimum to high energy. The profile for the first step of the path appears very similar for all eight parameterizations, with a systematic shift in the energy origin when the improved radial basis is used. All three potentials with the extended repulsive core move the C₁@C₅₉ minimum to a more realistic high energy. Once again, the paths for potentials (5) and (7) with the extended radial basis are significantly longer in terms of the number of steps and the integrated path length.

The number of steps in the paths to buckminsterfullerene from the selected bowl, C₂@C₅₈, and C₁@C₅₉ minima are shown in Table II for all the potentials. The results for the bowl can be compared with the landscapes for GFN2-xTB and SCC-DFTB from our previous report.¹² The longer paths for potentials (5) and (7) suggest a systematic effect. For pairwise potentials, catastrophe theory suggests that minima and transition states will merge and disappear as the range of the pair interaction increases, producing path segments with increasingly asymmetric barriers.⁶⁶ This prediction has been investigated for clusters and bulk systems in previous work.^67,68 In the present case, however, potentials (5) and (7) with the “turbo soap” radial basis have the same range as the ordinary SOAP representation, and so further analysis is required, which is beyond the scope of the current report.

To provide one further check on the global landscape for one of the new potentials, we relaxed all the transition states from the original landscape for potential (3), where the three key pathways to buckminsterfullerene now have appropriate energetic signatures. We then conducted additional discrete path sampling^28,29 runs. The PATHSAMPLE program was then used to harvest further OPTIM connection attempts for low-lying minima separated from the global minimum by relatively high barriers. The lowest minimum located for potential (3) is now buckminsterfullerene, and the second-lowest minimum, which is the correct stack 2 structure⁶⁰ with two edge-sharing pentagons, lies 1.35 eV higher. A disconnectivity graph^53,54 for the new landscape is shown in Fig. 9, with a magnification including only the lowest 1000 minima (plus the two minima that correlate with the lowest C₂@C₅₈ and C₁@C₅₉ structures for the original potential) in Fig. 10. We therefore, conclude that this potential provides a much better representation of carbon, as judged specifically for C₆₀.

As discussed in the Introduction, it is important to fit an interatomic potential using data that reports on all the configurations that are sampled in calculating the properties of interest. It is striking how small modifications of the potential fitting procedure can result in large changes in the energy landscape in regions far from the training data. In our particular case, using the “turbo soap” radial basis eliminated the erroneous global minimum, and when using the original SOAP radial basis, adding just two new configurations achieved the same objective. The use of an enhanced repulsive core potential helps carry out global exploration but also changes the details of the shortest paths between selected configurations, highlighting once again the sensitivity of the landscape paths to subtle details of the potential. Our results show how energy landscape tools can provide efficient validation and key additional data for the refinement of fitted potentials. Such refinement may be particularly important if the functional form of the interaction is not based on theoretical principles and limiting forms for the underlying interatomic or intermolecular interactions, in contrast to the case of typical empirical potentials.³ 

We have shown how feedback between landscape exploration and training can improve the potential by removing or shifting unphysical low energy states. A detailed analysis of structure, energetics, and pathways in this approach is a much more stringent test than an analysis of properties that depend on ensemble averages. Such quantities may not reveal problems with the potential manifested in calculations that depend on different regions of configuration space. For example, tunneling splittings depend on the height and shape of barriers between minima. Ideally, we seek a high resolution test of the potential throughout the range that may be sampled in the computation of all the emergent properties of interest. The field of ML force fields is currently evolving rapidly. Since the potential that we used here was published, several new potentials for carbon have appeared, using both GAP⁵⁶ and more advanced mathematical techniques, such as the atomic cluster expansion⁶⁹ and MACE,⁷⁰ an equivariant message passing network.⁷¹ We will report on global landscape exploration for these representations in future work.

The test set for validation of a fitted potential following training should really contain the global minimum. However, given the dimensionality of the problem, it is unlikely that low-lying unphysical states will be included without global optimization or enhanced sampling of some kind to locate them. A detailed analysis of the global energy landscape provides a powerful tool for identifying problems with a fitted potential and for guiding refitting to fix unphysical artifacts.

We thank Miguel Caro for some DFT calculations that helped to extend the GAP model.

The authors have no conflicts to disclose.

Gábor Csányi: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). John W. R. Morgan: Methodology (equal); Software (equal); Writing – review & editing (equal). David J. Wales: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

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