Comparison of predictive approaches to the dynamics of activated catalytic processes

Abstract

We compare two systematic approaches for constructing the kinetic transition network associated with a catalytic reaction, namely reactive global optimization (RGO) and discrete path sampling (DPS). We test convergence of pathways for selected steps of the dealloying processes occurring in Pt₂Mn slab models under oxidative conditions for a DFT-parametrized machine learning interaction potential (MLIP). We find close agreement between the approaches. In particular, both schemes resolve multistep transformations that appeared as single steps in a previous meta-dynamics (m-Dyn) treatment. RGO and DPS are therefore proposed as effective tools for the systematic exploration of reaction paths in catalysis.

---

1. Introduction

Defined within the Born–Oppenheimer approximation,¹ the potential energy surface (PES) of an atomically-described system plays a fundamental role in the physical description of chemically interesting systems. In particular, minima of the PES correspond to locally stable configurations, while saddle points with [(3N − 7, 1, 0) ≡ (+, −, 0)] signature in the internal Hessian² correspond to transition states (TS) that mediate the transition from one minimum to another via steepest-descent paths. Thus, the PES contains information about both thermodynamics (minima) and kinetics (TS) of chemical transformations.

To obtain predictive insights for such processes, information regarding minima and TS must be obtained for the kinetically relevant pathways, and a wide number of approaches are available in the literature, which we do not attempt to review here. A basic classification partitions these methods into single-ended and double-ended,^3,4 where the former start from a given minimum and identify neighboring structures (exploration), while the latter correspond to searches for paths between two known minima (connections, usually requiring intervening minima and the corresponding transition states). Examples of the former approach include metadynamics (m-Dyn),^5,6 the activation-relaxation technique,^7,8 our reactive global optimization (RGO) approach⁹ and the induced force method,¹⁰ while many double-ended schemes are based on the nudged (NEB)^11–14 or doubly-nudged¹⁵ (DNEB) elastic band chain of images. While some of these techniques are based on free energy surfaces (FES), and others on the PES, one should be able to establish a one-to-one correspondence between critical points, for example using the structures identified by any of these methods as an initial guess for a DFT minimization/transition state (TS) search with a consistent functional, basis set, etc.

In most chemical processes both thermodynamics and kinetics must be taken into consideration. For example, the global minimum for ethylene oxidation with O₂ is a mixture of CO₂ and H₂O. However, using suitable catalysts, the reaction stops at ethylene oxide or acetaldehyde (metastable states),^16,17 depending on the selectivity of the catalytic process. The different outcomes result from kinetic control of the reaction, indicating that exploration of the minima and TS of the PES should proceed together.

Among the possible approaches designed to properly and predictively account for kinetics, we focus here on the RGO algorithm, previously developed by some of us⁹ and described below. RGO starts from a given local minimum and climbs its way out of it with eigenvector-following (EF) searches, identifying the neighboring saddle points and the new minima on the opposite side of the TS, and selecting the kinetically accessible paths under the chosen conditions. The space of explored structures is expanded taking into account kinetic constraints, avoiding exploration of irrelevant directions when activation energies along the paths exceed selected thresholds. RGO is thus a single-ended, unbiased and kinetics-sensitive exploration approach, which is important in kinetics-controlled applications such as catalysis, where often, only the resting state is characterized and the higher-energy, low concentration reaction intermediates are unknown. We note that RGO is closely related to the activation-relaxation approach,^7,8 which also uses eigenvector-following for generating events and calculating barriers.

Although RGO has been applied in a number of cases both with first-principles methods^18,19 and empirical potentials,²⁰ a systematic assessment of its capabilities has not yet been attempted. Here we aim at providing such validation, by comparing RGO results against two well-established alternatives, m-Dyn^5,6 and discrete path sampling (DPS),^21,22 using a DFT-parametrized machine learning interaction potential (MLIP) to represent the PES. We now describe the rationale for these comparisons.

m-Dyn is a widely used single-ended method for enhancing sampling in molecular dynamics simulation using biases to reconstruct the free energy surface describing the given reaction mechanisms. The details of the m-Dyn technique and its application in various catalytic systems can be found in the original ref. 5 and 6.

On the other hand, DPS^21,22 corresponds to a double-ended approach that provides a well validated reference against which to measure the thoroughness of alternative approaches, thanks to the number of available schemes and systematic methods for expanding a kinetic transition network,^23–25 which are usually based on DNEB^15,26 searches to identify candidates for TS refinement with hybrid EF.^27–29 In this framework the connections for each TS are determined using approximate steepest-descent paths, and missing connections in multistep pathways are filled in by double-ended searches guided by the missing connection algorithm.³⁰ The resulting kinetic transition networks are refined by searches designed to shorten the overall path length between reactants and products, or to find alternative routes with lower barriers.^30,31 This approach has been used for clusters, biomolecules, condensed matter, and in surface science applications.^32,33 DPS has recently provided new insight into general optimisation problems with multiple minima, including machine learning loss functions,³⁴ quantum computing,³⁵ and self-consistent fields.³⁶ Moreover, while DPS is most often conducted using double-ended searches, it can also employ single-ended schemes. In fact, most of the RGO setup described below can be achieved within existing options of the PATHSAMPLE program, which is a driver for OPTIM. In this sense, RGO may be seen as a form of DPS based on a specific recipe of single-ended searches.

A point that will be relevant when discussing computational cost is that TS refinement within DPS is usually achieved via gradient-only methods,^27–29 which are much more practical for large systems. However, the dimension of the present application is small enough for EF with a full numerical Hessian to be feasible. This approach is especially convenient in the current RGO setup, where particular eigendirections are searched systematically. The corresponding eigenvectors can also be obtained within the gradient-only formulation used in hybrid EF,^27–29 if the Rayleigh–Ritz refinement is combined with successive orthogonalisation. However, there will be a progressive loss of precision, which makes diagonalisation of the full Hessian more attractive here. Furthermore, very recent advances,³⁷ developed after the completion of this work, should speed up the Hessian evaluation for MLIPs and dramatically reduce the computational cost of this approach.

As our test case, we revisit a previous study where some of us investigated the first steps of the dealloying process in a Pt₂Mn slab system under oxidative conditions.³⁸ Dealloying is the first discovered and ubiquitous catalyst restructuring process. However, despite its frequency and importance in catalysis, only a few atomistic studies are available in the literature, due to the intrinsic complexity of this process. Indeed, the atomistic mechanisms of the Pt–Mn system dealloying include the creation of (sub)surface vacancies, atom migration, element oxidation, …, providing representative examples of catalytic reaction mechanisms. A comparison and validation of different computational methods for this realistic example is thus a good test of the quality of both the different methods and of the comparison itself.

We find the present assessment of the RGO approach to provide favourable results. We show that, in the tested cases, the RGO exploration is putatively complete, in the sense that the network of neighboring minima of a starting minimum identified by the algorithm includes all of those obtained by the other methods, m-Dyn and DPS. Moreover, the lowest-energy path, as well as other potentially competitive alternative paths, identified by RGO precisely correspond to the ones characterised by DPS. Finally, both RGO and DPS are able to find atomistic transformation paths for the dealloying process that were not located by m-Dyn, and involve lower energy barriers leading to more stable configurations, thus both kinetically and thermodynamically favored. Hence the novelty of the present work: we demonstrate the congruence of two alternative complementary approaches to reactive dynamics (one basically single-ended and the other double-ended), and we test this agreement for an important dynamic phenomenon (restructuring) on an important process (catalyst de-alloying under oxidative conditions) and material (Pt-based catalyst models).

The report is structured as follows. In Section 2, we review the computational details. Then, in Section 3 we compare the results of RGO with those of m-Dyn and DPS, finally summarizing our conclusions in Section 4.

2. Computational details

RGO is an iterative algorithm that, starting from an initial minimum, builds the kinetic transition network^23–25 for the chemical system of interest. It requires as inputs a starting structure (min_I) and a way to calculate energies, forces and Hessians, either analytical or numerical. In its basic implementation considered here, it further requires at least two energy thresholds: one for the maximum acceptable activation energy above which a reaction is deemed too slow contribute significantly, and one for the maximum energy above the initial minimum below which a newly identified minimum is kept in the database. In other words, if either the activation energy is too high, or if the neighboring minimum is too high in energy, the corresponding branches of the graph are cut and not further investigated. In the present work we have not used algorithms to prune unfavorable links and accelerate the RGO search, so that the present simulations establish an upper bound to the RGO computational cost.

In the first step of the algorithm, the Hessian is calculated and diagonalized to yield its 3N − 6 non-0 eigenvectors. Then, we use an EF^27,29,39,40 algorithm to search for a TS along all the eigenvectors, both in the positive and in the negative directions (2 × (3N − 6) searches). If a TS is identified, the connected local minima are isolated by taking short steps parallel and antiparallel to the eigenvector corresponding to the negative eigenvalue and minimizing the energy to yield two structures (min_A and min_B). If either min_A = min_I or min_B = min_I, then the one found on the opposite side of the TS is added to the list of the nearest neighbours (NN) of the starting minimum, and the energy of the TS is saved. If the energy threshold requirements are satisfied, the minimum is checked for duplicates against all previously identified structures. If it is not already in the database, then the new minimum is saved in the list of the structures to be further analyzed. In case minimizing on either side of the TS doesn't yields the starting minimum, a second attempt is made with double-ended searches between min_I → min_A and min_I → min_B. If it is still impossible to confirm that a TS is directly connected to min_I, then the given normal mode is discarded. After all 2 × (3N − 6) directions have been analyzed, the algorithm selects a new structure from the list of the minima and iterates.

In the present work, we use a threshold of 0.9 eV for both TS and minima acceptance. This value is consistent with the reaction energy diagram predicted by m-Dyn for the Pt₂Mn oxidative de-alloying process we present in the Results section. More in general, an upper value of 0.9 eV is typical of many kinetic-controlled processes, such as the CH₃OH decomposition on Pd and Pt catalysts.⁴¹ Indeed, by assuming a typical value of about 10¹² s⁻¹ for the Arrhenius prefactor, a barrier of 0.9 eV corresponds to reactions with a relatively low turn-over-frequency of 1 s⁻¹ per catalytic site up to around 400 K. Finally, as will be clear below, this setup allows us to show that the algorithm is working properly, by correctly identifying all the barriers that are below the threshold and discarding the ones above. In applications, this threshold is a key hyperparameter controlling branch pruning and the size of the predicted final network of pathways: values too low can lead to unrealistically shallow minima as putative final products, whereas values too large explore unreasonable region of the PES that are not experimentally reachable.

In our implementation of the RGO, all minimisations, EFs, TS identifications and D-NEB searches use the OPTIM code.⁴² Energies and forces can be calculated either by OPTIM itself, if the selected forcefield model is included, or, as in our case, by the MD code LAMMPS. We exploit a purpose-developed interface between OPTIM and LAMMPS⁴³ to extend the applicability of the code to the large set of potentials available in LAMMPS, including recent MLIPs. We have recently made available on github examples of RGO input files and output for a test case.⁴⁴

As we will show in the results section, RGO is able to sample very thoroughly the PES thanks to its systematic single-ended search scheme based on second derivatives. For this reason, the computational cost of the algorithm is mainly related to the number of Hessian evaluations, which on its own it is related to the system size, the number of eigendirections explored, and the approach for calculating its entries (numerically, analytically, or automatically). In the present implementation, we stop the TS search if it requires more than 500 steps of EF, so that the theoretically worst situation would require the calculation of 500 steps × (3N − 6) × 2 directions Hessians. Eigenvalues are related to the number of steps needed to convergence, ranging from the 10–50 steps for the softest modes, to a few hundreds for the stiffer modes, indicating that an exploration of the lower barrier modes can be achieved significantly faster than a complete analysis, although this scheme can easily lead to misleading results in case of fast equilibria that do not lead to the actual products of interest (e.g. amide hydrolysis).

Since neither analytic nor automatic second derivatives are yet available for the BP-MLIP, we calculate the Hessians numerically. Also, for completeness, we search along all eigenvector directions (i.e., without using the pruning techniques explored in previous work). As such, we use the most computational expensive setup to show that for a system containing around 60 atoms, which is a typical size for slabs/unit cells used in QM models of catalytic processes on extended facets, the RGO approach is feasible in a couple of weeks of wall time on a 48-core CPU node, retrieving reaction paths missed by less systematic approaches. In the SI, we estimate that the same simulation could be 1–2 orders of magnitude faster using non-numerical second derivatives on a system of similar size. In this sense, recent advances³⁷ on fast Hessian calculations for MLIP potentials (not available when this work was completed) are very promising.

3. Results and discussion

As mentioned above, our test system is taken from one of our previous reports³⁸ in which we derived a Behler–Parrinello (BP) forcefield to study a complex oxidative dealloying of a Pt₂Mn catalyst surface by means of m-Dyn single-ended searches. Here, we use the most complex of these paths (the one from Fig. 5 of ref. 38) to compare RGO with m-Dyn. Also, since we now have a suggested path, we can also perform double-ended searches between previously determined local minima to compare the RGO results with the systematic DPS approach. Since our goal is to assess the performance of the RGO method and its reliability we adopt a consistent PES description by running all calculations (m-Dyn, RGO, DPS) on the same BP-MLIP PES. The MLIP parameters and a compatible LAMMPS input are given in the SI.

In Fig. 1 we show the reaction energy diagram and the structures identified in our previous work with the BP-MLIP energetics.

Fig. 1 Reaction steps leading to oxidation of the Pt₂Mn system starting from a configuration with two Mn in the surface layer and adsorbed oxygen adatoms, corresponding to Fig. 5 of ref. 38, with MLIP energetics. (top) Reaction energy diagram (RED), punctuated by local minima (IS = initial state, IntN = N-th intermediate, FS = final state) and transition states (TSN = N-th transition state). (bottom) Atomistic depiction of local minima along the oxidation path, from a side-view (top row) or top-view (bottom row). Color coding: O in red, Mn in pink, Pt in gray.

We focus on the Int1–Int2 step as a first test. Running an OPTIM connection search reveals that the path connecting these two configurations (directly linked in m-Dyn) is actually composed of three steps, as shown in the right half of Fig. 2. Such differences can occur in general, since m-Dyn forces the system out of its initial local minima, but does not guarantee that the predicted mechanism corresponds to a single-step process on the PES. Interestingly, running the RGO approach with Int1 as the initial state, we find and identify exactly the same minimum energy path as that predicted by DPS with double-ended searches, to end up in Int2. As a test, we also ran the test in the opposite direction, that is, starting RGO from Int2 towards Int1. In this case, the algorithm correctly identifies the first two reverse barriers, but does not find the third step back to Int1. However, this result was expected and is actually consistent with the kinetic strategy underlying RGO, since the third reverse barrier is higher than the 0.9 eV threshold, meaning that this branch must be successfully cut as it will not occur under the chosen kinetic conditions (irreversible oxidation).

Fig. 2 (top) RED and (bottom) atomistic depiction of local minima and saddle points predicted by RGO for the steps IS → Int1 (left half) and Int1 → Int2 (right half). Image scheme and color coding as in Fig. 1.

Analyzing the output of RGO starting from Int1, to support our claim of completeness, we should also be able to find a path leading from Int1 to IS. This path, as identified by the m-Dyn, contains a finite-size-model artifact: a plane slipping from a close-packed ABC to ABA layer structure. Starting from Int1, RGO is able to find IS after three iterations steps, along the path shown in the left half of Fig. 2. This is the same path found by double-ended DPS approach, with the last step being the plane slipping. Since RGO using the BP-MLIP potential is CPU intensive due to the numerical Hessian, we ran only one step of RGO starting from IS, to confirm that the first step of the IS → Int1 path is the same plane slip movement with the same activation energy, but in the opposite direction. We note that, in addition to Int1 → IS and Int1 → Int2, RGO identifies another 10 minima directly connected to Int1 after first iteration, and 162 minima after second one, 20 of which are up to 0.2 eV more stable than Int2, confirming that the m-Dyn search is limited to a subset of the actual reaction mechanisms, and that systematic explorations such as RGO or DPS are needed to predictively account for catalytic kinetics.

We now move our focus to the Int3–Int4 step. According to the OPTIM calculations, this pathway involves 12 TS to be crossed, as shown in Fig. 3. It thus represents a useful validation of RGO on a more complex and composite mechanism. To speed up the path search, we started RGO from both sides, from both Int3 and Int4 to “meet in the middle”, i.e. until we find one or more common minima to both searches. Also, after a few iterations, we selected the structures corresponding to the ones found by OPTIM along the path and reset the RGO algorithm from these structures, to avoid dedicating too much effort in investigating alternative paths that correspond to branchings not sampled by DPS (which certainly exist, given the complexity of the path and the systematic character of RGO). We found that the RGO algorithm was able to reconstruct the full path connecting Int3 to Int4, with the two ends of the search overlapping halfway through. This example shows that RGO can deal successfully even with long multi-step rearrangements well beyond the m-Dyn force-driven kinetics. To conclude our investigation, for all of the remaining steps, we limited ourselves to a simplified protocol in which we checked that, among all the structures identified after the first RGO iteration, there always were the ones that the OPTIM connection algorithm had identified as the first step towards both Int(n − 1) and Int(n + 1). We note that many of the intermediates identified by both RGO and DPS lie at lower energies than Int4 and are also kinetically favoured. For instance, Int3.3 lies at ≈−0.3/−0.4 eV with respect to both Int3 and Int4, and can be reached from Int3 with an overall energy barrier lower than ≈0.25 eV, which is lower than the energy barrier from Int3 to Int4 according to DFT, i.e., 0.36 eV. The correct dynamics is therefore Int3 → Int3.3 as predicted by both RGO and DPS, not Int3 → Int4 as predicted by m-Dyn.

Fig. 3 (top) RED and (bottom) atomistic depiction of local minima and transition states predicted by RGO for the step Int3 → Int4. Image scheme and color coding as in Fig. 1.

At the end of this discussion, we should comment on the selected energy thresholds. We note that all the energy barriers we have determined in the minimum-energy path are well below the chosen threshold value of 0.9 eV. This choice guarantees that we identified a network that is reachable significantly faster than the competing paths that we may have discarded because of the threshold. Indeed, we could have decreased the threshold value down to around 0.6 eV and still predict identical reaction paths.

Finally, we remark again that the comparison between FES based m-Dyn and PES based RGO or DPS is meaningful. Catalytic processes generally involve activated processes characterized by a complex PES, with well-defined local minima separated by energy barriers appreciably larger than the available thermal energy.⁴⁵ Under these conditions, there is a clear correspondence between the minima and transition states of the FES and PES, so that a comparison of approaches is useful.

4. Conclusions

The recent developments of increasingly more accurate MLIPs capable of predicting activated processes to useful chemical accuracy with respect to the underlying first-principles approaches⁴⁶ raise the need for similarly efficient methods that can characterise the kinetically relevant transformation paths under given conditions. A thorough understanding of kinetic transition networks would then allow us to fully and predictively capture kinetically controlled processes, such as the catalytic cycles, and would improve our understanding of MLIP reliability.^47,48 Here we take a paradigmatic example from our previous work, the m-Dyn^5,6 investigation of dealloying for a Pt₂Mn slab models under oxidative conditions on a BP-MLIP³⁸ PES, and we revisit it using systematic sampling algorithms that resolve all the underlying steps, namely, reactive global optimization (RGO)⁹ and discrete path sampling.^21,22 From this comparison, we show how these methods based on unbiased geometry optimisation can discover competing reaction paths and resolve the composite nature of transformations that appear as single-steps in m-Dyn runs. Alternative paths with lower energy barriers are characterised, leading to more stable configurations, which are both kinetically and thermodynamically favored. The multiple TS candidates are also unlikely to be resolved in single NEB or DNEB chain-of-states runs, but are easily determined by DPS with iterative double-ended searches and the missing connection algorithm.³⁰ The multistep paths have been cross-validated in the RGO and DPS runs, where we see a systematic coincidence of the kinetic networks for the given mechanistic steps produced by these two different but complementary sampling methods. Hence RGO and DPS are ready to be employed in the search of novel catalytic materials and processes that can solve the present global problems in the energy and environmental fields.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the major findings of this study are available in the supplementary information (SI) or from the corresponding authors upon reasonable request. Supplementary information is available, in which an estimate of the computational cost of RGO algorithm and LAMMPS input are provided. See DOI: https://doi.org/10.1039/d6cp01329e.

The BP-MLIP can be found in a Zenodo repository: https://zenodo.org/uploads/19473718. The software associated with RGO and path sampling approaches is freely available on github.

Acknowledgements

We thank the Italian CINECA supercomputing center for computational resources within the ISCRA program. We acknowledge financial support from ICSC, the Italian Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by the the European Union (Next Generation EU, grant number CN00000013) – PNRR, Mission 4 Component 2 Investment 1.4, and from the Italian Ministry of Environment and Energy Security POR H2 AdP MMES/ENEA project, funded by the European Union (Next Generation EU) – PNRR, Mission 2, Component 2, Investment 3.5 “Research and development on hydrogen”. Networking within the COST Action CA21101 “Confined molecular systems: from a new generation of materials to the stars” (COSY) supported by COST (European Cooperation in Science and Technology) and within the International Research Network IRN nanoalloys is also acknowledged.

Notes and references

1. M. Born, Ann. Phys., 1927, 84, 457–484.
2. J. N. Murrell and K. J. Laidler, Trans. Faraday Soc., 1968, 64, 371–377.
3. P. M. Zimmerman, J. Comput. Chem., 2015, 36, 601–611.
4. E. F. Koslover and D. J. Wales, J. Chem. Phys., 2007, 127, 134102.
5. A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
6. A. Barducci, M. Bonomi and M. Parrinello, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2011, 1, 826–843.
7. G. Barkema and N. Mousseau, Phys. Rev. Lett., 2000, 77, 4358–4361.
8. R. Malek and N. Mousseau, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, 7723–7730.
9. F. R. Negreiros, E. Aprà, G. Barcaro, L. Sementa, S. Vajda and A. Fortunelli, Nanoscale, 2012, 4, 1208–1219.
10. S. Maeda, T. Taketsugu and K. Morokuma, J. Comput. Chem., 2014, 35, 166–173.
11. G. Mills, H. Jónsson and G. K. Schenter, Surf. Sci., 1995, 324, 305–337.
12. H. Jónsson, G. Mills and K. W. Jacobsen, Classical and quantum dynamics in condensed phase simulations, World Scientific, Singapore, 1998, ch. 16, pp. 385–404.
13. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
14. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
15. S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2004, 120, 2082–2094.
16. A. Evnin, J. Rabo and P. Kasai, J. Catal., 1973, 30, 109–117.
17. T. Pu, H. Tian, M. E. Ford, S. Rangarajan and I. E. Wachs, ACS Catal., 2019, 9, 10727–10750.
18. F. R. Negreiros, A. Halder, C. R. Yin, A. Singh, G. Barcaro, L. Sementa, E. C. Tyo, M. J. Pellin, S. Bartling, K. H. Meiwes-Broer, S. Seifert, P. Sen, S. Nigam, C. Majumder, N. Fukui, H. Yasumatsu, S. Vajda and A. Fortunelli, Angew. Chem., Int. Ed., 2018, 57, 1209–1213.
19. C. R. Yin, F. R. Negreiros, G. Barcaro, A. Beniya, L. Sementa, E. C. Tyo, S. Bartling, K. H. Meiwes-Broer, S. Seifert, H. Hirata, N. Isomura, S. Nigam, C. Majumder, Y. Watanabe, S. Vajda and A. Fortunelli, J. Mater. Chem. A, 2017, 5, 4923–4931.
20. M. Asgari, F. R. Negreiros, L. Sementa, G. Barcaro, H. Behnejad and A. Fortunelli, J. Chem. Phys., 2014, 141, 041108.
21. D. J. Wales, Mol. Phys., 2002, 100, 3285–3306.
22. D. J. Wales, Mol. Phys., 2004, 102, 891–908.
23. F. Noé and S. Fischer, Curr. Opin. Struct. Biol., 2008, 18, 154–162.
24. D. Prada-Gracia, J. Gómez-Gardeñes, P. Echenique and F. Falo, PLoS Comput. Biol., 2009, 5, e1000415.
25. D. J. Wales, Curr. Opin. Struct. Biol., 2010, 20, 3–10.
26. D. Sheppard, R. Terrell and G. Henkelman, J. Chem. Phys., 2008, 128, 134106.
27. L. J. Munro and D. J. Wales, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 3969–3980.
28. G. Henkelman and H. Jónsson, J. Chem. Phys., 1999, 111, 7010–7022.
29. Y. Kumeda, L. J. Munro and D. J. Wales, Chem. Phys. Lett., 2001, 341, 185–194.
30. J. M. Carr, S. A. Trygubenko and D. J. Wales, J. Chem. Phys., 2005, 122, 234903.
31. B. Strodel, C. S. Whittleston and D. J. Wales, J. Am. Chem. Soc., 2007, 129, 16005–16014.
32. D. J. Wales, J. P. K. Doye, M. A. Miller, P. N. Mortenson and T. R. Walsh, Adv. Chem. Phys., 2000, 115, 1–111.
33. K. Röder, J. A. Joseph, B. E. Husic and D. J. Wales, Adv. Theory Simul., 2019, 2, 1800175.
34. D. J. Wales, Annu. Rev. Phys. Chem., 2018, 69, 401–425.
35. B. Choy and D. J. Wales, J. Chem. Theory Comput., 2023, 19, 1197–1206.
36. H. G. A. Burton and D. J. Wales, J. Chem. Theory Comput., 2021, 17, 151–169.
37. A. Burger, L. Thiede, N. Rønne, V. Bernales, N. Vijaykumar, T. Vegge, A. Bhowmik and A. Aspuru-Guzik, arXiv, 2025, preprint, arxiv:2509.21624 DOI:10.48550/arxiv.2509.21624.
38. T. Roongcharoen, X. Yang, S. Han, L. Sementa, T. Vegge, H. A. Hansen and A. Fortunelli, Faraday Discuss., 2023, 242, 174–192.
39. C. J. Cerjan and W. H. Miller, J. Chem. Phys., 1981, 75, 2800–2806.
40. J. Pancír̆, Collect. Czech. Chem. Commun., 1974, 40, 1112–1118.
41. T. Roongcharoen, G. Conter, L. Sementa, G. Melani and A. Fortunelli, J. Chem. Theory Comput., 2024, 20, 9580–9591.
42. OPTIM: A program for geometry optimisation and pathway calculations, https://www-wales.ch.cam.ac.uk/software.html.
43. Nicy, J. W. R. Morgan and D. J. Wales, J. Chem. Phys., 2024, 161, 054112.
44. G. Conter and A. Fortunelli, RGO-code, https://github.com/afloer/RGO-code, 2025.
45. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
46. T. Roongcharoen, G. Conter, L. Sementa, G. Melani and A. Fortunelli, J. Chem. Theory Comput., 2025, 21, 11164–11178.
47. G. Csányi, J. W. R. Morgan and D. J. Wales, J. Chem. Phys., 2023, 159, 104107.
48. V. Cărare, F. L. Thiemann, J. D. Morrow, D. J. Wales, E. O. Pyzer-Knapp and L. Dicks, npj Comput. Mater., 2026, 12, 9.

---