How to Model a Reaction

Why bother with computation?

If you've never used computational chemistry, you might wonder: why model a reaction when you could just run it in the lab?

Computation can give considerable insight by revealing mechanistic-level details that are difficult to probe experimentally. It can help answer questions like:

• What is the geometry of the transition state (TS)?
• Which interactions stabilize it?
• How do bonds form and break along the reaction path?
• Is the reaction stepwise or concerted?

Modeling can also be faster and higher-throughput than traditional lab work, particularly when studying late-stage intermediates. With Rowan, you can quickly tweak a reactant and recalculate properties of interest, predict enantioselectivity, or scale up to screen many reactions at once.

This post is a practical overview of how to use computational chemistry to model reactions. To learn more about the science behind modeling reactions and locating transition states, please see our post "Reactions from the Bottom Up."

Step 1: Prepare your endpoint structures

Before you can model a reaction, you need to obtain 3D structures of your reactants.

If you already have a SMILES string, Rowan can auto-generate 3D structures for downstream modeling. If you don't already have a SMILES string, it can easily be copied from ChemDraw or drawn from scratch using Rowan's 2D structure editor.

If your molecule has a high number of rotatable bonds, it can be helpful to run Rowan's "conformer search" workflow on your reactant and product structures to better understand their conformational landscapes.

To calculate the thermodynamic favorability ∆G₀ of a reaction (and not the kinetic barrier ∆G‡) all you need to do is calculate free energies for the reactants and products. You can skip straight to step 5!

Step 2: Find a guess transition state (TS)

Once you have your reactant and product structures handy, the next step is finding a guess transition-state structure. This is necessary because transition-state optimization algorithms work best when they start from a geometry that's already close to the desired transition state. Finding that geometry is often one of the hardest parts of modeling a reaction: a poor guess may fail to optimize or may converge to a transition state for a different process (like a methyl rotation).

There are two standard ways in Rowan to obtain a guess TS: a more manual way using scans and a more automatic way using guess-TS-finding algorithms.

Use a scan to find a guess TS

Scans can be used to calculate the energetic profile of moving along a reaction coordinate. For simple reactions, it's often easy to define a single bond coordinate that approximates the true reaction coordinate. In bimolecular association/dissociation reactions, it's often easiest to start from the unimolecular structure and scan towards the separated structures.

For more complex reactions, it can be helpful to run a coordinated scan that scans along multiple coordinates simultaneously. For example, when scanning to find the TS of this Diels–Alder reaction, two bond coordinates should be specified:

Once the scan is complete, identify the highest-energy point along the scanned path. If the original scan used coarse steps, rerun the region around this maximum with more, smaller steps to obtain a geometry closer to the TS. Then resubmit the highest-energy geometry from the refined scan as the starting point for a TS optimization.

For a worked example using scans to find TS structures, watch our video "Finding a Transition State."

Use a double-ended TS search to find a guess TS

Rowan's "double-ended transition-state (TS) search" workflow supports a number of automatic guess-TS-finding algorithms, including the freezing string method (FSM), the growing string method (GSM), and the nudged elastic band (NEB) method. These double-ended TS search methods start from both ends—the reactants and products—and find a feasible reaction path connecting them. The highest-energy geometry along that path is often a good guess for the TS.

For these algorithms to work properly, the atom indexing of the reactant and product system must match exactly. This is often best achieved by starting with your product structure, duplicating it, and manually separating the product into its component reactants before queueing the double-ended TS search. For help with atom reindexing, Rowan has a tool that can be used to reindex simple structures (on the double-ended TS search submit page); for complex cases, we recommend doing it manually or enlisting the help of your favorite LLM.

Another common failure mode is poor alignment between the reactant and product geometries. Visually inspect the resulting pathway: if it contains a large conformational change near either the reactant or product end, the endpoint geometries are probably poorly aligned. Reorient the structures so that corresponding atoms and unchanged parts of the molecules overlap as closely as possible, then rerun the search.

For more information about double-ended TS search algorithms, please see our post "Guessing Transition States."

For a worked example using the double-ended TS search, watch our video "Studying Azide–Alkyne Cycloaddition With the Freezing-String Method."

Step 3: Optimize your TS

Once you've obtained a guess TS structure, you'll need to optimize it and calculate its frequencies. Rowan performs both tasks through the "basic calculation" workflow: select the "optimize (TS)" and "frequency" tasks when configuring the workflow. The expected result of the TS optimization is one imaginary frequency (notated using a negative number) that corresponds to the reaction coordinate.

Look for one imaginary frequency

Frequencies in Rowan can be animated by clicking on the corresponding row. For the intended TS, the animated imaginary frequency should look like the reaction happening: the bonds you expect to form and break should move in the corresponding directions.

If multiple imaginary frequencies are found, you likely have an optimization problem; nudge the structure slightly along the problematic vibrational modes and re-optimize. It's common to find trivial transition states like methyl group rotations with the first few rounds of TS optimizations, so don't lose hope!

Step 4: Confirm you've found the right TS

Once you've optimized a TS structure, confirm that the imaginary frequency corresponds to the intended reaction coordinate and that the TS connects the expected reactants and products. The first check can be performed by visualizing the imaginary frequency, while Rowan's "intrinsic reaction coordinate (IRC)" workflow provides a more rigorous confirmation of the connected endpoints. If your molecule is conformationally complex, you should also consider whether lower-energy TS conformers exist.

Visualize the imaginary frequency

The simplest way to assess the identity of the TS is to visualize which atoms are moving in the imaginary frequency. Ideally, reactive atoms will dominate the imaginary mode, while ancillary bond rotations will be minor components at best. Some reactions (e.g. hydrocupration, E2-type eliminations, Alder ene reactions) have large and characteristic imaginary frequencies associated with light-atom motion.

Run an intrinsic reaction coordinate (IRC) calculation

An IRC calculation starts from an optimized TS and follows the path of steepest descent or "fall line" in both the forward and backward direction to find the reactant and product structures that the starting TS connects. If you've found a real TS structure, an IRC will yield your intended reactant and product structures.

For information about how IRC calculations work, read our science docs.

TS conformer searching

To run a TS conformer search, constrain all the atoms involved in the TS (whether using bond constraints, angle constraints, or by freezing atoms) when running your conformer search. Once you've obtained several conformers, you will need to optimize and confirm the validity of each independently.

Step 5: Calculate barrier heights

Once you've found a guess TS, optimized it, and confirmed that it actually connects your reactant and product structures, the final step is to optimize the IRC-derived endpoints to minima (this is an option in the workflow) and ensure that your level of theory is consistent across all your calculations.

For final calculations, it can be good to use a higher level of theory or a more diffuse basis set to reduce error. Transition-state theory predicts that the difference in barrier between the reactant structures and transition-state structure will determine the rate of reaction.

Include thermochemistry corrections

For accurate barrier-height energies, use the "total Gibbs free energy" value available in Rowan's thermochemistry tab. For more information about thermochemical corrections, read our thermochemistry documentation.

To properly account for translational entropy, it is necessary to optimize and calculate thermochemistry corrections for separate structures in separate calculations.

For worked example using Rowan to find barrier heights, watch our video "Computing Barrier Heights." Rowan maintains a simple web utility for turning Rowan calculations into potential-energy-surface graphs at https://labs.rowansci.com/graph.

Additional notes

Different classes of reactions

Some classes of reactions are much harder to model than others. In particular, take care with reactions that involve:

• Large changes in charge distribution. Charge-separated intermediates and transition states are often poorly described without an implicit- or explicit-solvent model.
• Strong solvent dependence. For example, you may struggle to find stable S_N2 or E2 transition-state structures without an implicit-solvent model.
• Large entropic contributions. Dissociative reactions such as S_N1 and E1 can be especially difficult because their barriers exist on the free-energy surface, not necessarily on the electronic-energy surface.

For more information about frequent challenges in reaction modeling, we recommend Ryu et al.'s "Pitfalls in Computational Modeling of Chemical Reactions and How To Avoid Them."

Choosing a level of theory

For many of the below steps, you'll need to choose a level of theory.

For qualitative work, we recommend using an OMol25-trained neural network potential (NNP) like OMol25's eSEN, UMA, or OrbMol. If you prefer physics-based methods or observe strange behaviors, we recommend using the GFN2-xTB semiempirical method for geometry finding and g-xTB for more accurate energies.

If you aim to publish your work, we recommend finding work in your field to serve as a reference, though we find that the r²SCAN-D4/vDZP level of theory is a good starting place for a lot of work (we suggest this in place of the oft-cited B3LYP/6-31G(d) combination). For computationally focused work, we recommend using a range-separated hybrid such as ωB97M-D3BJ for final publication values.

If your reaction is in solvent, be sure to use an implicit solvent model for energies. The inclusion of an implicit solvent model can sometimes slow down or confuse guess-TS-finding algorithms; omit if necessary. For low-cost methods like NNPs and semiempirical methods, we recommend the CPCM-X solvent model for correcting single-point energies.

Modeling multi-step reactions

This guide was written for single-step reactions; multi-step reactions are modeled as a series of single-step reactions, so TS structures for each step will need to be obtained and confirmed separately. For many applications, however, it's only necessary to model the rate-determining TS.