Two-dimensional RMSD projections for reaction path visualization and validation

This paper introduces a novel two-dimensional visualization framework that maps reaction trajectories onto a permutation-corrected RMSD plane with gradient-enhanced Gaussian Process energy interpolation, enabling more effective comparison and validation of optimization histories across different computational methods for complex chemical reactions.

The Gist

Imagine you are trying to describe a complex hiking trip from a valley (the Reactant) to a mountain peak (the Transition State) and down to a new valley (the Product).

In the world of computational chemistry, scientists use powerful computers to simulate these trips. However, the standard way they report the results is like giving you a one-dimensional graph that only shows your altitude (energy) against the number of steps you took.

The Problem with the Old Way:
If you only look at that graph, you don't know where you were. Did you take a straight, safe path? Did you wander off a cliff and climb back up? Did you take a scenic route or a dangerous shortcut?

• The "Step Count" Trap: If two hikers take different routes but take the same number of steps, the old graph looks identical. You can't tell if they found the same peak or if one of them got lost in a different valley.
• The "Altitude" Trap: It tells you how high you are, but not if you are on the right mountain or a different one entirely.

This makes it very hard to compare different hiking strategies (algorithms) or different maps (computer models) to see if they are actually finding the same path.

---

The New Solution: A 2D "Topographic Map"

The author, Rohit Goswami, proposes a new way to visualize these chemical reactions. Instead of a flat line, he creates a 2D topographic map (a landscape with hills and valleys) that shows the journey in a much clearer way.

Here is how the magic works, broken down into simple analogies:

1. The "Ruler" that Doesn't Care About Names (Permutation-Invariant RMSD)

Imagine you have a bag of LEGO bricks. You build a castle (Reactant) and then a spaceship (Product).
In the old way, if you swapped the names of two red bricks in your computer code, the computer might think the spaceship is completely different, even if it looks the same.

• The Fix: This new method uses a special "smart ruler" (called RMSD) that measures the distance between your current shape and the start/end shapes. It doesn't care which brick is which; it only cares about the overall shape. It's like saying, "How far is this pile of LEGOs from the Castle shape?" and "How far is it from the Spaceship shape?"

2. The Two-Axis Compass (The 2D Projection)

Instead of just tracking "steps taken," the new map uses two coordinates:

• Axis A: How far are you from the Start?
• Axis B: How far are you from the Finish?

By plotting every step of the journey on this 2D grid, you get a clear picture of the path.

• The "Turn": If the hiker takes a sharp turn, the line on the 2D map bends.
• The "Detour": If the hiker wanders off into a bush, the line moves away from the direct path.
• The "Shortcut": If they find a faster way, the line looks different.

3. Filling in the Blanks (The AI Painter)

The computer only calculates the energy at specific points (like taking photos every 10 meters). The spaces between the photos are empty.

• The Fix: The author uses a smart AI (a Gaussian Process) to "paint" the rest of the map. It looks at the photos and the forces (pushes/pulls) the atoms feel, then fills in the missing terrain with a smooth, colorful surface.
• The Safety Net: The AI also draws dashed lines around areas where it is guessing. If you see a dashed line, it means, "I'm not 100% sure what the terrain looks like here because we didn't take any photos nearby." This prevents scientists from trusting guesses too much.

4. Rotating the Map (Reaction Progress vs. Deviation)

Finally, the author rotates the map so it's easier to read:

• The "Progress" Line: This is the direct path from Start to Finish.
• The "Wiggle" Line: This shows how far off-course the hiker went.
  This helps scientists instantly see if a method is staying on the path or getting confused.

---

Why Does This Matter? (The Real-World Impact)

The paper tests this on three different chemical reactions:

1. A Simple Collision: Like two cars merging. The new map showed that two different computer models (one using AI, one using traditional math) found the exact same "mountain pass" (Transition State), even though they took slightly different routes to get there. The old 1D graph couldn't prove this.
2. A Complex Twist (Grignard Rearrangement): Imagine a dancer spinning and changing direction. The old graph just showed "energy went up, then down." The new 2D map showed a sharp 90-degree turn in the path, revealing that the computer was struggling to find the right move near the start.
3. A Bumpy Ride (Bicyclobutane): The new map showed a sharp "elbow" in the path right after the peak, confirming the AI model was accurate.

The Bottom Line

Think of the old method as a fitness tracker that only tells you your heart rate and total steps.
Think of this new method as a GPS map that shows you exactly where you walked, where you got lost, and whether you found the right destination.

In short: This paper gives chemists a better "map" to see if their computer simulations are actually finding the correct chemical paths, rather than just guessing based on a single number. It helps them trust their AI models and avoid wasting time on dead ends.

Technical Summary

1. Problem Statement

In computational chemistry, locating transition states and minimum energy paths (MEPs) is routine, typically using methods like the Nudged Elastic Band (NEB) or String methods. However, the standard analysis and visualization of these trajectories suffer from significant limitations:

• Dimensional Reduction: Standard plots rely on 1D "profile" graphs (Energy vs. Image Number or Energy vs. Cumulative Reaction Coordinate). These reduce high-dimensional configuration space ($3N$ dimensions) to a single scalar axis.
• Loss of Geometric Context: The 1D reaction coordinate is often history-dependent (based on the specific optimization path taken) and arbitrary. It obscures structural rearrangements, making it impossible to compare optimization histories from different algorithms (e.g., NEB vs. Frozen String) or different potentials beyond simple metrics like barrier height or total time.
• Validation Ambiguity: Distinguishing between numerical instability and physical relaxation into alternative basins is difficult. Validating whether a machine-learned potential (MLIP) captures the correct barrier topology often requires manual inspection of 3D structures, which is impractical for large datasets.

2. Methodology

The authors propose a framework to map high-dimensional optimization trajectories onto a 2D projection defined by permutation-invariant Root Mean Square Deviation (RMSD) distances to reactant and product configurations. The pipeline consists of four main stages:

A. Intrinsic Projection Coordinates $(r, p)$

Instead of using arbitrary reaction coordinates, the method defines coordinates based on the distance of any configuration $X$ from the reactant ($R$) and product ($P$) states:

• Metric: $d(X, X_{ref}) = \min_{Q, \Pi} \sqrt{\frac{1}{N} \|X - Q X_{ref} \Pi\|_F^2}$.
• Permutation Invariance: To handle varying atom indexing, the method uses the Iterative Rotations and Assignments (IRA) algorithm to find the optimal rotation $Q$ and permutation $\Pi$. This ensures the coordinates are invariant to rigid body motion and atom labeling.
• Coordinates: $r = d(X, R)$ and $p = d(X, P)$.

B. Reaction Progress Frame $(s, d)$

The raw $(r, p)$ plane contains unphysical regions. The authors apply a rigid rotation to decompose the plane into:

• $s$ (Reaction Progress): The coordinate along the path direction (tangent).
• $d$ (Orthogonal Deviation): The coordinate perpendicular to the path.
  This frame separates the progress of the reaction from deviations, analogous to path collective variables but without requiring pre-defined descriptors or smoothing parameters.

C. Synthetic Gradient Construction

Since the RMSD map involves non-differentiable minimization over permutations, direct projection of Cartesian forces is mathematically intractable. Instead, the method constructs synthetic gradients in the 2D space:

• The path tangent vector $(\tau_r, \tau_p)$ is computed using Savitzky-Golay smoothing of the $(r, p)$ coordinates.
• The parallel force component ($F_{\parallel}$) from the NEB calculation is projected onto this tangent to estimate the gradient of the energy surface in the 2D plane: $\nabla E \approx -F_{\parallel} \cdot \tau$.

D. Energy Surface Interpolation via Gaussian Processes (GP)

To visualize the continuous energy landscape $E(s, d)$ from sparse discrete data points:

• Model: A Gradient-Enhanced Gaussian Process is used.
• Kernel: The Inverse Multiquadric (IMQ) kernel, $k(x, x') = (c^2 + r^2)^{-1/2}$, is chosen for its strictly positive definiteness and polynomial decay (heavy tails), which helps capture long-range basin structures from sparse data.
• Inputs: The GP is trained on both energy values ($E$) and the synthetic gradients ($\nabla E$).
• Uncertainty Quantification: The GP posterior variance is calculated and overlaid as contours. Low variance indicates data-supported regions; high variance indicates extrapolation.
• Scalability: For large systems (e.g., crystalline slabs), a Nyström low-rank approximation is employed to reduce computational complexity from $O(N^3)$ to $O(NM^2)$.

3. Key Contributions

1. Coordinate-Free Visualization: A method to visualize reaction paths without requiring prior knowledge of collective variables (e.g., bond lengths, angles). It relies solely on Cartesian coordinates and endpoints.
2. Cross-Method Comparability: By using a universal metric (RMSD to endpoints), the framework allows rigorous comparison of trajectories generated by different algorithms or potentials, which was previously impossible with 1D profiles.
3. Geometric Validation: The 2D landscape reveals whether different methods converge to the same barrier region or distinct basins, even if the final saddle geometries differ slightly.
4. Reliability Indicators: The inclusion of posterior variance contours provides a built-in metric to distinguish between interpolated data and extrapolated artifacts.
5. Open-Source Implementation: The authors provide a command-line tool (rgpycrumbs) to generate these plots, making the method accessible for routine validation.

4. Results

The framework was validated on three distinct chemical systems:

• 1,3-Dipolar Cycloaddition (Ethylene + N2O):
  • Compared a Machine-Learned Potential (MLIP) path against a DFT (B3LYP) reference.
  • The 2D projection showed that despite geometric displacements, both saddles lie on comparable energy contours, validating the MLIP's barrier topology.
  • The posterior variance clearly delineated the reliable region near the path.
• Grignard Rearrangement:
  • Demonstrated the method's ability to capture complex, curved trajectories (90-degree turns in the $(s, d)$ plane).
  • Revealed optimization noise near the reactant basin that was invisible in 1D profiles.
• Conrotatory Bicyclobutane Ring Opening:
  • Showed a sharp geometric "elbow" post-transition state.
  • Confirmed that the MLIP path faithfully traversed the saddle region identified by the Free String Method (FSM).

In all cases, the 2D projection provided qualitative insights (e.g., "are we in the same basin?") that 1D scalar plots could not offer.

5. Significance and Limitations

Significance:

• Bridging the Gap: The method fills a critical gap between impractical 3D structural inspection and information-poor 1D scalar summaries.
• MLIP Validation: It offers a robust way to verify that machine-learned potentials reproduce the correct topology of reaction barriers, not just the final energy values.
• General Applicability: It is applicable to any double-ended path search (NEB, String) and can be extended to single-ended searches or molecular dynamics post-hoc.

Limitations:

• Lossy Projection: The map from $3N$ to 2D is many-to-one. Distinct configurations can share the same $(r, p)$ coordinates. The surface represents a conditional expectation, not the full Potential Energy Surface (PES).
• No Orthogonal Curvature: The method uses only tangential forces. It cannot recover or enforce the property that MEP images are minima in the orthogonal hyperplane.
• Permutation Complexity: For systems with large permutation groups (e.g., large molecules), the IRA algorithm's cost increases, though it remains tractable for the systems tested.
• Post-Hoc Only: The method visualizes existing calculations; it does not accelerate the NEB optimization itself.

Conclusion:
This work introduces a powerful, coordinate-free visualization tool that enhances the interpretability and validation of reaction path calculations. By transforming high-dimensional optimization histories into intuitive 2D landscapes with uncertainty quantification, it enables researchers to rigorously compare different computational methods and potentials, ensuring that the identified transition states are physically meaningful.