Automating Transition-State Search
strings and bands; interpolating between states; searching for transitions; confirmation; Vicena integration; recent blogs
We’re releasing a completely rewritten reaction toolkit today, coinciding with our presentation about it at the Gordon Research Conference on “Molecular Discovery and Simulation in the AI Era” in Barcelona this week. In brief:
We’ve rewritten everything from scratch to be more robust and performant, using only NumPy and SciPy as dependencies.
We’ve added the minimum-energy-path-generation methods GSM and CI-DNEB to our existing FSM.
We’ve also added geodesic and IDPP interpolation methods.
Our IRC code now implements the HPC/EulerPC method and is significantly faster.
Both of these codes have been updated to support both periodic and non-periodic systems.
Here’s a bit more on the science behind this, what we’ve changed, and the impact that this has for practical tasks in Rowan. (Also see the end of this newsletter for details on our new Vicena integration and two new blog posts on protein–ligand complexes.)
Transition-State Searching
Locating a transition state (TS) for a reaction is a difficult problem. Even when the precise mechanism is known, properly aligning the forming and breaking bonds at the correct distances takes considerable practice. TS optimizers have small convergence radii, and small errors in placement can lead to the optimization failing.
Worse, small changes in the substituents can shift the transition state earlier or later, causing a good guess to fail when an R-group is swapped. This difficulty is compounded by the cost of the quantum-chemical methods commonly used for optimizations. TS optimization starts and ends with the calculation of a Hessian (matrix of mixed second derivatives); a semi-numerical Hessian on a 100-atom system requires 300 gradient calls. Thus, computational chemists typically spend significant amounts of time setting up a guess geometry; having to repeat the calculation multiple times is costly and time-consuming.
Even with the advent of low-cost methods like NNPs, good guess geometries are still hard to come by. Simply interpolating between the reactants and products will shove atoms through each other on their way to their final destination, producing nonsense. Smart interpolation methods such as IDPP and geodesic can steer the reactions around obvious clashes and avoid breaking/forming too many bonds, while interpolation in redundant internal coordinates can better express the large-scale motions that are happening (e.g. backbone rotations). However, even the best interpolation methods still produce highly strained species. Information from the true potential-energy surface is needed to help the path avoid high-energy intermediates. Methods of doing so are often divided into two classes: string methods and band methods.
Band methods stretch an interpolation across the surface, sample images along it, and then optimize these images subject to spring forces from adjacent images (preventing images from collapsing into the minima). Nudged elastic band (NEB) is the dominant method; it adds spring forces between images and projects out spring forces orthogonal to the tangent at the image that would otherwise pull atoms off the path. Images are relaxed under these constraints, approximating the true path. Many variants exist, including doubly-nudged NEB (DNEB) and climbing-image (CI-NEB), which accelerate the determination of the path.
In contrast, string-based methods iteratively grow strings from each endpoint towards the opposite endpoint. They take a step along the interpolation, place a node, and optimize the nodes orthogonal to the existing string. This process is repeated until the string ends meet. The growing-string method (GSM) loosely optimizes the whole of each string every time new nodes are added, and performs a tighter post-connection sweep. The freezing-string method (FSM) eschews post-placement optimization, only loosely optimizing new nodes before adding them to the string. FSM produces paths that do not follow the minimum-energy path and are often significantly more strained, but it is often an order of magnitude cheaper than GSM, and doesn’t risk nodes sliding back into the endpoint region.
We’ve implemented FSM and GSM with BFGS Hessian updates on each node and copy the Hessian from the previous node when adding new nodes, accelerating optimization. We included both climbing image and double nudging in our NEB implementation, and we use the FIRE optimizer to accelerate band convergence due to the non-conservative band forces. This rewrite also allows us to treat periodic and non-periodic systems on equal footing and tightly integrate the code with our existing optimizer and NNP infrastructure, allowing us to realize significant improvements in speed and robustness. For a longer write-up on how double-ended TS search methods work, see our blog from last year.
Intrinsic Reaction Coordinates (IRC)
Once a putative transition state has been found, it is important to use an IRC to verify that it truly connects the reactants and products. It is not uncommon for a proposed transition state to instead connect two minima that are nothing like the reactants or products, or connect one conformer to another. IRCs are particularly important for complex or highly congested reactions, where multiple branching paths could be possible, but we hope that they can become standard practice for all published transition states.
IRCs follow the mass-weighted minimum-energy path (the same path as a perfectly damped molecular-dynamics trajectory). Our IRC code implements the Hessian-based predictor–corrector (HPC/EulerPC) algorithm of Hratchian & Schlegel in mass-weighted Cartesian coordinates:
Compute the Hessian at the TS.
Displace along the imaginary mode.
Find the point on the hypersphere that is the lowest in energy:
Build a distance-weighted interpolant (DWI) by taking the energy and gradient at a predictor step.
Use the DWI surface to analytically determine the true minimum on the hypersphere.
Determine the length along the path via the DWI.
Repeat until max steps or a minimum is reached.
Repeat in the opposite direction.
This method is effectively 4th order, improving on our previous IRC workflow which used the 2nd-order Gonzalez–Schlegel (GS2) IRC from geomeTRIC. The use of Cartesian coordinates as opposed to TRIC avoids the difficulties from back-transforming internal coordinates, allowing it to quickly solve for the next step and to work well for periodic systems.
We’ve also added a few additional features to our workflow to simplify the job of a chemist. The IRC workflow provides initial TS optimization, endpoint optimization, and Hessian refinement instead of full Hessian recompute when the Bofill-update Hessian becomes ill-conditioned. We’ve also added support for periodic systems and support a wide variety of NNPs, semiempirical methods, and even DFT.
Benchmarks
The Baker TS set is commonly used for benchmarking TS-optimization and TS-search methods. It comprises 25 small-molecule transition states generated with HF/3-21G. We re-optimized the TSs with UMA-S-1.2-OMol, generated endpoints with IRCs, and benchmarked how various search and interpolation combos performed.
Our geodesic interpolation method (implemented following Zhu, Thompson, and Martínez) outperforms all other interpolations on the Baker set. It penalizes bond rearrangement, favoring interpolations that generate the fewest bond changes (technically, it finds the geodesic on a Riemannian manifold with a Morse-scale interatomic-distance metric). The systems in this benchmark set are too small to see the benefits of our redundant internal coordinates and often fail the iterative back transformation due to the tightly curved nature of the paths (it has a comparative advantage for large systems with dihedral-driven motion, e.g. biphenyl rotations). FSM is an order of magnitude faster than other methods, with only marginally reduced accuracy. It can produce highly strained species (see figure below), but much of the strain is orthogonal to the true reaction coordinate, and thus the subsequent partitioned rational-function TS optimization quickly optimizes it out.
These methods also work well for larger systems, or when significant rearrangement occurs. The oxidative addition of methyl bromide to a gold complex splits the methyl group and bromine and places them on opposite ends of the complex. While FSM produces a highly strained initial guess (>100 kcal/mol above the TS), the subsequent optimization is still able to locate the TS.
Reactions in periodic systems often involve complex motions that cannot simply be decomposed into bond stretches or bends. Along with the large expense of periodic DFT gradients, TS search is often difficult to perform. However, geodesic interpolation does an excellent job of steering atoms, while avoiding the back-transformation difficulties inherent in redundant internal coordinates. For an interstitial carbon migration in BCC Fe, FSM and NEB trace similar paths and find the same transition state, with FSM completing significantly faster.
We’ve added numerous graphical updates to our double-ended TS-search and IRC workflows, including the ability to compare workflows, stream the results in near real time without the save call blocking the main compute thread, and included the post-search TS optimization and frequencies as tabs in the GUI. All of these search, interpolation, and IRC methods are now available to our users in both the GUI and via our rowan-python API.
We hope to continue improving our reaction modeling tools and build new ones for modeling reactivity at scale. If you are interested in reaction modeling at scale or find any reactions where these TS-search methods consistently fail, please reach out. We are always trying to improve our tools, and establishing a compendium of difficult systems is helpful for continuing to make them better.
Vicena Integration
We’re excited to share that Rowan now integrates with Vicena, an agentic scientific workspace that (in their words) brings research tools, evidence, and results into one continuing project, so researchers can focus on the scientific question rather than managing each tool separately. Rowan’s suite of simulation tools integrates naturally with Vicena’s search and reasoning capabilities; you can read more on specific use cases on Vicena’s website.
Here’s what Vicena CEO Mehdi Farzanepour has to say on the integration:
The value of bringing Rowan and Vicena together is that researchers get the depth and reproducibility of computational simulations alongside the breadth and flexibility of a unified agentic research. This helps close the research loop by connecting ideas, literature, data, calculations, and results, then using what is learned to refine the question and decide what to do next.
As part of the integration, Mehdi is offering Rowan users a 50% discount for 4 months on subscriptions to Vicena. If you’re interested, you can sign up for Vicena with code LETSDOMORESCIENCE (valid through September 30th).
Blogs
We’ve also recently published two blogs that might be of interest to our readers in drug discovery:
“What to Do With a Pose,” an introduction to how to analyze protein–ligand poses in Rowan. We’ve noticed a lot of early Rowan users get stuck at this stage so we decided to compile all our advice and relevant features into one place.
Some case studies benchmarking Rowan’s hydration-site-analysis functionality, demonstrating that our implementation recovers known hydration sites from the literature and showing off compatibility with larger molecular glues. This hydration-site analysis feature is now live for all Rowan subscribers in the pose-analysis MD workflow!
Until next time, happy computing!