Tutorial: Saddlepoint ("TS") optimization via relaxed scan

- The very cheap (and not very accurate) semi-empirical method PM3 is used here (all these calculations can be run in seconds on any computer). Production calculations would of course be performed at the DFT level at least (use analytical Hessian in that case instead of the numerical one below). Using a cheap method like PM3 (more accurate alternatives are HF-3c or PBEh-3c) to crudely explore potential energy surfaces or try out and get familiar with surface scans and TS optimizations in ORCA, before switching to the more expensive method, is not a bad idea though. Just keep in mind that the potential energy surfaces can sometimes be completely different and it may not always save time.

Step 1: Relaxed Surface Scan

In order for the TS-optimization jobs to work at all, one needs to have previously found a geometry close to the TS and guessing it is usually too hard. One way is to find an internal coordinate that resembles the real reaction coordinate sufficiently well and run a relaxed surface scan.

More sophisticated methods such as the nudged elastic band method require no pre-defined choice of the reaction coordinate.

A relaxed surface scan involves constrained optimizations for different values of a reaction coordinate and could be done manually but is done more conveniently using the Scan keyword of the %geom block. For the S_N2 reaction above, the distance between F and C (or Cl and C) is is the obvious reaction coordinate. Starting with the F^- ion in 3.0 Å distance from the CH₃Cl molecule we run a relaxed scan until the distance is 1.2 Å.

Important: Many users forget in this step that ORCA counts atoms from 0, not from 1.

! PM3 Opt
%geom Scan
B 5 0 = 3.0, 1.2, 15 # Scanning the distance between atom 5 and atom 0 from distance 3.0 Å to 1.2 Å in 15 points.
end
end

*xyz -1 1
C -1.55384 -0.43723 0.00000
Cl -0.48384 -0.43723 0.00000
H -1.91051 -1.20940 0.64919
H -1.91051 -0.61336 -0.99331
H -1.91051 0.51106 0.34413
F -4.54939 -0.44897 0.08390
*

This job will trace a path (close to a minimum-energy path) from the reactants (F^- + CH₃Cl) to the products (Cl^- + CH₃F ). The ORCA output will give a final table giving the (relaxed) energy for each point of the scan. ORCA will also create xyz files that contain the optimized coordinates for each point of the scan. Note also that the Chemcraft program allows convenient visualization of relaxed surface scans, by reading in the ORCA outputfile directly.

**** RELAXED SURFACE SCAN DONE ***

SUMMARY OF THE CALCULATED SURFACE

----------------------------

RELAXED SURFACE SCAN RESULTS

----------------------------

Column 1: NONAME

The Calculated Surface using the 'Actual Energy'

3.00000000 -33.86955186

2.87142857 -33.87079549

2.74285714 -33.87167482

2.61428571 -33.87179810

2.48571429 -33.87087148

2.35714286 -33.86892508

2.22857143 -33.86651022

2.10000000 -33.86552921

1.97142857 -33.87141654

1.84285714 -33.88704580

1.71428571 -33.90914046

1.58571429 -33.93516450

1.45714286 -33.95871984

1.32857143 -33.96517287

1.20000000 -33.93045275

From Figure 1 we see that the highest point on the curve is at the coordinate with value 2.1 Å. This would make a good guess for an OptTS job although sometimes even finer scans are required (change 15 to a higher number in the %geom block).

Optimizing to a saddlepoint (often informally referred to as a "transition state") is often a bit tricky. Use of the OptTS keyword (eigenvalue following algorithm) directly, works only if the geometry is very close to the saddlepoint state already. A relaxed surface scan is one way to find an approximate minimum energy path connecting 2 minima and the saddle-point is then located by eigenvector-following method. This approach often works well for finding saddlepoints of simple bond making/breaking reactions or conformational changes, but sometimes more sophisticated minimum-energy-path methods are required, e.g. the nudged elastic band method.

In this tutorial, we demonstrate first how to find a saddlepoint using a relaxed surface scan and an eigenvector-following saddlepoint optimization.

In the end, we demonstrate how this can also be accomplished using the NEB method (new in ORCA 4.1). See Minimum energy path calculations and the NEB tutorial for more information on the NEB method and how to use it.

- As an example reaction, finding the saddlepoint of the S_N2 reaction CH₃Cl + F^- -> Cl^- + CH₃F is used in all input files below

Figure 1. Energy as a function of reaction coordinate. Scan ran from point 3.0 Å to 1.2 Å.

Step 2a: TS optimization from point close to TS

If you are close to a saddlepoin already, then starting a simple OptTS job will do the trick if you are lucky (uses an approximate Hessian). The following job (starting from the 2.1 Å point from the Scan job above) will converge to a TS and a subsequent NumFreq job confirms it. This job finds the TS in 14 optimization cycles but these jobs are not guaranteed to work and often you end up in a minimum instead of a TS.

! PM3 OptTS NumFreq

*xyz -1 1
C -1.99283303076949 -0.43906204115263 0.01309547131865
Cl 0.08654162182067 -0.43514846620583 -0.01623959206210
H -2.09531444966205 -1.26265127699215 0.70686814906647
H -2.12108831451492 -0.62708960496606 -1.04443065188875
H -2.10329090994441 0.57197871671963 0.38153156676382
F -4.09261491692978 -0.44315732740297 0.04308505680191
*

Step 2b: TS optimization with exact Hessian

What is usually necessary, even if you are close to saddlepoint, is to calculate the exact Hessian in the beginning at least. This job finds the TS in 3 optimization cycles but requires the calculation of the Hessian (usually quite expensive but is more likely to succeed than doing the calculation without the Hessian). A hybrid Hessian could also be used.

! PM3 OptTS NumFreq

%geom
Calc_Hess true # Calculate Hessian in the beginning
NumHess true # Request numerical Hessian (analytical not available)
Recalc_Hess 5 # Recalculate the Hessian every 5 steps
end

The NumFreq output confirms that a TS was indeed successfully found (1 imaginary frequency).

VIBRATIONAL FREQUENCIES

-----------------------

0: 0.00 cm**-1

1: 0.00 cm**-1

2: 0.00 cm**-1

3: 0.00 cm**-1

4: 0.00 cm**-1

5: 0.00 cm**-1

6: -503.24 cm**-1 ***imaginary mode***

7: 173.47 cm**-1

The imaginary vibrational frequency is easily visualized by using the orca_pltvib program. Either the ORCA output or the .hess file can be used.

orca_pltvib ch3cl-f-TSOpt-fromImage8-Hess.out 6

This would create an xyz trajectory file, ch3cl-f-TSOpt-fromImage8-Hess.out.v006.xyz that can be visualized by programs such as VMD, Chemcraft, Avogadro and others.

Opening the ORCA outputfile directly in Chemcraft or Avogadro also works for visualizing all frequencies.

New alternative: Finding the minimum energy path and/OR saddlepoint via NEB, NEB-CI or NEB-TS

A new feature in ORCA 4.1 is an efficient implementation of the nudged elastic band method (NEB), it's climbing-image variant (NEB-CI) and a new combination that combines NEB, CI-NEB and the OptTS procedure (NEB-TS). See NEB tutorial for more information. The inputs below show how to find the same saddlepoint using the NEB, NEB-CI and NEB-TS methods instead which would be the generally recommended way.

NEB job

If the NEB keyword is used, the elastic band is minimized until it converges to the minimum energy path.

Convergence tolerances for all images are by default: Max(|F|): 5.00E-4 Eh/Bohr and RMS(F): 2.00E-4 Eh/Bohr.

! PM3 NEB
%neb
NEB_End_XYZFile "prod.xyz"
Nimages 6
end
*xyzfile -1 1 react.xyz

This job converges the elastic band in 20 iterations (122 energy+gradient evaluations) which is remarkably good. Most systems will not behave like this. The highest energy image on this path should be reasonably close to the saddlepoint (here the energy is off by 0.6 kcal/mol) but is NOT the saddlepoint. This highest energy image could be used in a separate OptTS job to find the saddlepoint (like above).

NEB-CI job

If the NEB-CI keyword is used instead, the climbing-image variant of NEB is used. Here only an approximate minimum energy path is achieved (higher convergence thresholds on the images by default) but the calculation converts (after threshold max(|Fp|)=2.00E-2 Eh/Bohr is reached) one image into a "climbing" image that converges to the saddlepoint (with low convergence threshold). This is usually much more efficient than converging the full minimum energy path.

Convergence tolerances for regular images are by default: Max(|Fp|): 5.00E-2 Eh/Bohr and RMS(Fp): 2.50E-2 Eh/Bohr.

Convergence tolerances for climbing image ("saddle-point") are by default: Max(|F|): 5.00E-4 Eh/Bohr and RMS(F): 2.50E-4 Eh/Bohr.

Note that the thresholds for this job can also be lowered so that an accurate minimum energy path with an optimized saddlepoint is achieved. Usually only the saddlepoint is of main interest so this may not always be worth the effort.

! PM3 NEB-CI
%neb
NEB_End_XYZFile "prod.xyz"
Nimages 6
end
*xyzfile -1 1 react.xyz

This job converges the elastic-band+CI in 26 iterations (158 energy+gradient evaluations).

NEB-TS job

If the NEB-TS keyword is used instead, a CI-NEB job is carried out but with high convergence threshold on regular images, medium convergence threshold for CI. Once the NEB-CI part converges (with those higher thresholds), the climbing image with an approximate Hessian is fed into an OptTS job.

Convergence tolerances for regular images are by default: Max(|Fp|): 1.00E-1 Eh/Bohr and RMS(Fp): 5.00E-2 Eh/Bohr.

Convergence tolerances for climbing image are by default: Max(|F|): 1.00E-3 Eh/Bohr and RMS(F): 5.00E-4 Eh/Bohr.

Convergence tolerances for saddle-point (OptTS part) are by default: Max(|F|): 3.00E-4 Eh/Bohr and RMS(F): 1.00E-4 Eh/Bohr.

! PM3 NEB-TS
%neb
NEB_End_XYZFile "prod.xyz"
Nimages 6
end
*xyzfile -1 1 react.xyz

This job converges the elastic-band+CI in 21 iterations and one OptTS iteration (128 energy+gradient evaluations).