Two-dimensional RMSD projections for reaction path visualization and validation

Transition state or minimum energy path finding methods constitute a routine component of the computational chemistry toolkit. Standard analysis involves trajectories conventionally plotted in terms of the relative energy to the initial state against a cumulative displacement variable, or the image number. These dimensional reductions obscure structural rearrangements in high dimensions and are often history dependent. This precludes the ability to compare optimization histories of different methods beyond the number of calculations, time taken, and final saddle geometry. We present a method mapping trajectories onto a two-dimensional projection defined by a permutation corrected root mean square deviation from the reactant and product configurations. Energy is represented as an interpolated color-mapped surface constructed from all optimization steps using a gradient-enhanced Gaussian Process with the inverse multiquadric kernel, whose posterior variance contours delineate data-supported regions from extrapolated ones. A rotated coordinate frame decomposes the RMSD plane into reaction progress and orthogonal distance. We show the utility of the framework on a cycloaddition reaction, where a machine-learned potential saddle and density functional theory reference lie on comparable energy contours despite geometric displacements, along with the ratification of the visualization for more complex reactions, a Grignard rearrangement, and a conrotatory bicyclobutane ring opening.

💡 Research Summary

**
The paper introduces a novel two‑dimensional (2D) visualization framework for reaction‑path and transition‑state (TS) searches that overcomes the limitations of traditional one‑dimensional (1D) energy‑profile plots. Conventional analyses of nudged elastic band (NEB), climbing‑image, or string calculations typically display the relative energy versus image index or a cumulative Euclidean distance along the path. Such scalar reductions discard the high‑dimensional geometric information contained in the 3N Cartesian coordinates of each image, making it impossible to compare different algorithms or to detect subtle structural deviations that may indicate convergence problems or alternative reaction channels.

Key methodological steps

1. Permutation‑invariant RMSD coordinates – For each image X_i the authors compute two scalar distances: r_i = RMSD(X_i, R) and p_i = RMSD(X_i, P), where R and P are the reactant and product reference structures. To make these distances independent of atom ordering and overall rotation/translation, they employ the Iterative Rotations and Assignments (IRA) algorithm, which simultaneously finds the optimal permutation matrix Π and rotation matrix Q that minimize the Frobenius norm. This yields a unique, physically meaningful pair (r_i, p_i) for every geometry.

2. Rotation to reaction‑progress (s) and orthogonal deviation (d) – The raw (r, p) plane contains an unphysical region where both distances are simultaneously small. The authors therefore define a rigid rotation based on the line connecting the first and last images in the (r, p) plane. The rotated axes are:

   • s_i = (r_i – r_0)·ŝ_r + (p_i – p_0)·ŝ_p (progress along the reaction)
   • d_i = (r_i – r_0)·đ_r + (p_i – p_0)·đ_p (perpendicular deviation) where ŝ and đ are unit vectors parallel and normal to the line. The resulting (s, d) coordinates retain absolute distance information while providing an intuitive reaction coordinate (s) and a measure of how far a given geometry strays from the ideal path (d).

3. Synthetic gradients from NEB forces – NEB calculations provide the total forces F_i on each image. The component parallel to the path, F∥,i = (F_i·τ̂_i)τ̂_i, where τ̂_i is the unit tangent, is used to construct projected gradients in the (r, p) space:

   • ∇_r E ≈ –F∥·τ_r
   • ∇_p E ≈ –F∥·τ_p The authors first smooth the (r, p) trajectory with a Savitzky‑Golay filter to obtain stable tangents τ_r, τ_p.

4. Gradient‑enhanced Gaussian Process (GP) regression – The data set for the GP consists of the scalar energies E_i and the two projected gradients (∇_r E_i, ∇_p E_i). They adopt the inverse‑multiquadric (IMQ) kernel k(x, x′) = (c² + ‖x – x′‖²)^{‑½}, which is strictly positive‑definite for any dimension because its generating function is completely monotone. The IMQ kernel decays polynomially (∝ r^{‑1}) rather than exponentially, allowing it to capture long‑range basin structure from the sparse NEB samples. Derivative observations are incorporated analytically, tripling the information content per image without extra energy evaluations. Hyper‑parameters (c and noise levels) are learned on a small subset of final‑path images, reducing the computational cost from O(N³) to O(n³) with n ≈ 20.

5. Energy surface and uncertainty visualization – The GP posterior mean yields a continuous, color‑mapped energy surface E(s, d). The posterior variance is plotted as contour lines, clearly distinguishing regions where the model is well‑constrained by data from extrapolated zones. This provides an immediate visual diagnostic of convergence quality and of any “holes” in the sampled space.

Demonstrations

• Cycloaddition reaction – The authors compare a density‑functional theory (DFT) NEB path with a machine‑learned potential (MLP) path. Both lie on essentially the same energy contour, confirming that the MLP reproduces the correct barrier height, yet the RMSD projection reveals a small (~0.1 Å) geometric offset between the two saddles.
• Grignard rearrangement – The 2D surface displays multiple low‑energy “valleys” and intervening barriers that are invisible in a 1D plot, illustrating the method’s ability to expose competing pathways.
• Conrotatory bicyclobutane ring opening – The reaction progresses with a pronounced deviation in the d‑direction, and the GP variance highlights a region where the NEB images are sparse, suggesting where additional images would improve the description.

Comparison with other dimensionality‑reduction techniques

The authors argue that principal component analysis (PCA), t‑SNE, UMAP, and Sketchmap all require either a priori feature selection or a large number of samples (10⁴–10⁵) to learn a reliable mapping. In contrast, the RMSD‑based projection works directly on the raw Cartesian coordinates, needs only the ~10³ images generated by a typical NEB calculation, and preserves absolute distances to the reactant and product. Consequently, it enables quantitative cross‑method comparisons without any “circular” dependence on prior chemical intuition.

Limitations and future directions

Because the projection collapses 3N‑1 degrees of freedom onto two, information about orthogonal curvature of the true potential energy surface is lost; only the parallel force component can be projected without analytically differentiating the permutation‑invariant RMSD. The authors suggest extending the framework to multiple endpoints or to incorporate orthogonal force information via more sophisticated Jacobian approximations. They also envision applying the method to large catalytic surfaces, adsorbate networks, or reaction‑network explorations where conventional collective variables are unavailable.

Conclusion

The paper delivers a practical, mathematically rigorous, and chemically intuitive tool for visualizing and validating reaction‑path calculations. By mapping each NEB image onto a permutation‑invariant RMSD plane, rotating to a reaction‑progress frame, and interpolating energies with a gradient‑enhanced GP, the authors provide a color‑coded energy landscape together with uncertainty estimates. This approach surpasses traditional 1D energy profiles, facilitates direct comparison between different electronic‑structure methods or machine‑learned potentials, and offers a clear diagnostic of convergence and possible alternative pathways. It fills a methodological gap in the toolbox of computational chemists who need reliable, post‑hoc visual validation of transition‑state searches.