Transition state search (Baker)

Simulate > Search Transition State (Baker)... refines a structure that is already close to a transition state into a first-order saddle point: a stationary point whose Hessian has exactly one negative eigenvalue. It runs pDynamo's eigenvector-following method of J. Baker (J. Comput. Chem. 7, 385–395, 1986), the same function used in Example 11 of the pDynamo book.

The method works in three steps:

• A numerical Hessian is computed once, at the starting structure.

• At every step the search goes uphill along one Hessian mode (the reaction coordinate) and downhill along all the others.

• The Hessian is not recomputed: after each step it is updated from the change in the gradient (unless Recompute Hessian every is used).

The method refines a transition state; it does not find a reaction path. The starting structure must already be near the saddle point. Good starting points are the energy maximum of a PES scan (Simulate > Reaction Coordinate Scans), the highest image of an NEB path (Simulate > Nudged Elastic Band) or of a CPR path (Simulate > Conjugate Peak Refinement...), or a hand-built guess with the forming and breaking bonds at intermediate lengths. The search runs as an ordinary job in the Process Manager, with the same log, trajectory and import handling as a geometry optimization.

The window (Figure 1)

• Coordinates: the starting structure, chosen from the coordinate sets of the active system.

• Log frequency (default 1): one line is written to the log every N steps, with the energy, RMS gradient and number of negative Hessian eigenvalues.

• Maximum iterations (default 250): if the search reaches this number, the job is reported as Finished (not converged).

• RMS Grad. Tolerance (default 0.001): the search stops when the RMS gradient falls below this value (kJ mol-1 Å-1).

• Follow mode (default 0): the Hessian mode followed uphill, counted from the lowest eigenvalue. At the first step the search climbs the chosen mode; after that it keeps following the same mode by its overlap with the previous step, even if the eigenvalue order changes. 0 is right when the starting structure is close to the transition state, where the lowest mode (with negative curvature) is the reaction coordinate. Use 1, 2, ... when the reaction coordinate is a higher mode, for example when starting from a minimum whose lowest modes are methyl rotations; run Simulate > Normal Modes on the starting structure and pick the mode that moves the atoms of the reaction.

• Hessian updating (default BOFILL): how the Hessian is updated between steps, BOFILL (recommended) or POWELL. BFGS and Murtagh–Sargent are not offered: in tests on cyclohexane, BFGS did not converge in 250 steps (it keeps the Hessian positive, which is wrong for a saddle point) and Murtagh–Sargent "converged" back to the chair minimum.

• Recompute Hessian every (default 0): recompute the exact Hessian every N steps; 0 means it is computed once and then only updated. Recomputing needs fewer steps and is more robust far from the transition state, but each recalculation is expensive: pDynamo has no analytic Hessian, so it costs 6 gradient evaluations per free atom (cheap with MM or semi-empirical QC, expensive with ab initio or DFT). For cyclohexane the search takes 123 steps with N = 0, 86 with N = 10, 80 with N = 5 and 27 with N = 1, all reaching the same transition state.

• Save Trajectory: writes the steps to a .ptGeo folder together with the log.

EasyHybrid runs its own copy of pDynamo's Baker optimizer, because the pDynamo3 version stops after 3 steps with an internal error when Follow mode > 0, and its "analytic Hessian" steps compute nothing for molecular systems. The algorithm is otherwise unchanged.

Worked example: cyclohexane chair → half-chair transition state

This is Example 11 of the pDynamo book, run from the EasyHybrid interface. The search starts from the chair conformation, a minimum, and climbs to the saddle point of the ring flip in a few seconds. The structure is $PDYNAMO3_HOME/examples/book/data/mol/cyclohexane_chair.mol (also in examples/tutorials/05_neb_cyclohexane of the EasyHybrid folder) and the OPLS parameters are in $PDYNAMO3_HOME/parameters/forceFields/opls/bookSmallExamples.

• Load the system: click New Project, set the force field to OPLS, select the bookSmallExamples parameter folder (the default, protein, has no cyclohexane types) and import cyclohexane_chair.mol. The energy of the chair is 82.56 kJ/mol (Simulate > Single Point Energy).

• Run the search: open Simulate > Search Transition State (Baker)..., keep the chair in Coordinates and all the defaults (250 iterations, tolerance 0.001, Follow mode 0, BOFILL), and click Run.

• Check the result: the search converges in 123 steps (final RMS gradient < 0.001), at about 177.8 kJ/mol, that is, +95.3 kJ/mol above the chair. The log ends with "Minimization converged.", which pDynamo writes for saddle-point searches too.

• Confirm the saddle point: run Simulate > Normal Modes on the new coordinate set. There must be exactly one negative (imaginary) frequency, about −311 cm-1 with these parameters (Figure 2 shows the same check for another cyclohexane model).

The barrier is a property of the small OPLS test parameters of the pDynamo book, not a reference value for the real ring flip.

When the search does not converge

If the search reaches Maximum iterations before the gradient falls below the tolerance, the Process Manager shows Finished (not converged) and the "Finished!" dialog and the log show a warning with the number of steps and the final RMS gradient. The structure from such a run is not a transition state (cyclohexane stopped after 50 steps with tolerance 0.1 has an RMS gradient of 3.1 and two imaginary frequencies). What to do:

• Restart from the last structure with more iterations. The search keeps no state between runs, so a restart recomputes the Hessian at the new point, which often helps.

• Start from a better guess, for example the maximum of a finer PES scan.

• Recompute the Hessian every few steps (1 to 10).

• Try POWELL instead of BOFILL.

If pDynamo itself reports an error, the job is marked Error! and no structure is imported.

After the transition state: the reaction path

From a converged transition state, Simulate > IRC / Reaction Path... follows the imaginary mode downhill in both directions, towards reactants and products. It checks again that the Hessian has exactly one negative eigenvalue, so a structure from a non-converged search is rejected there.

Practical notes

• Cost: each numerical Hessian takes 6N gradient calculations for N free atoms. For QC/MM systems, keep the number of free atoms small (fix distant atoms).

• Fixed atoms are honoured: the search runs only over the free coordinates.

• Restraints are part of the energy. Remove PES-scan restraints before searching, or the result is a saddle point of the restrained surface.

• Units: energies in kJ/mol and the gradient tolerance in kJ mol-1 Å-1, as everywhere in pDynamo.

Reference: J. Baker, "An algorithm for the location of transition states", J. Comput. Chem. 7, 385–395 (1986).