(sec:optimization.typical)= # Geometry Optimizations, Surface Scans, Transition States, MECPs, Conical Intersections, IRC, NEB The usage of analytic gradients is necessary for efficient geometry optimization. In ORCA 5.0, the following methods provide analytic first derivatives - Hartree-Fock (HF) and DFT (including the RI, RIJK and RIJCOSX approximations) - MP2, RI-MP2 and DLPNO-MP2 - TD-DFT for excited states - CAS-SCF When the analytic gradients are not available, it is possible to evaluate the first derivatives numerically by finite displacements. This is available for all methods. The coordinate system chosen for geometry optimization affects the convergence rate, with redundant internal coordinates being usually the best choice. Some methods for locating transition states (TS) require second derivative matrices (Hessian), implemented analytically for HF, DFT and MP2. Additionally, several approaches to construct an initial approximate Hessian for TS optimization are available. A very useful feature for locating complicated TSs is the Nudged-Elastic Band method in combination with the TS finding algorithm (NEB-TS, ZOOM-NEB-TS). An essential feature for chemical processes involving excited states is the conical intersection optimizer. Another interesting feature are MECP (Minimum Energy Crossing Point) optimizations. For very large systems ORCA provides a very efficient L-BFGS optimizer, which makes use of the orca_md module. It can also be invoked via simple keywords described at the end of this section. (sec:optimization.geometry.typical)= ## Geometry Optimizations Optimizations are fairly easy as in the following example: ```{literalinclude} ../../orca_working_input/C05S02_118.inp :language: orca ``` An optimization with the RI method (the BP functional is recommend) would simply look like: ```{literalinclude} ../../orca_working_input/C05S02_119.inp :language: orca ``` An optimization of the first excited state of ethylene: ```{literalinclude} ../../orca_working_input/optimization-tddft.inp :language: orca ``` (sec:optimization.numgrad.typical)= ## Numerical Gradients If the analytic gradient is not available, the numerical gradient can simply be requested via: ```orca ! NumGrad ``` as in the following example: ```{literalinclude} ../../orca_working_input/C05S02_121.inp :language: orca ``` NOTE - Be aware that the numerical gradient is quite expensive. The time for one gradient calculation is equal to 6 $\times$ (number of atoms) $\times$ (time for one single point calculation). - The numerical gradient can be calculated in a multi-process run, using a maximum of three times the number of atoms (see section {ref}`sec:parallel.calling`). More details on various options, geometry convergence criteria and the like are found in section {ref}`sec:geomopt.detailed`. (sec:optimization.notes.typical)= ## Some Notes and Tricks :::{note} - `TightSCF` in the SCF part is set as default to avoid the buildup of too much numerical noise in the gradients. - Even if the optimization does not converge, the ORCA output may still end with \"\*\*\*\*ORCA TERMINATED NORMALLY\*\*\*\*\". Therefore do not rely on the presence of this line as an indicator of whether the geometry optimization is converged! Rather, one should instead rely on the fact that, an optimization job that terminates because the maximum number of iterations has been reached, will generate the following output message: ::: :::{Warning} The optimization did not converge but reached the maximum number of optimization cycles. Please check your results very carefully. ::: While a successfully converged job will generate the following message instead: ``` ***********************HURRAY******************** *** THE OPTIMIZATION HAS CONVERGED *** ************************************************* ``` :::{tip} - In rare cases the redundant internal coordinate optimization fails. In this case, you may try to use `COPT` (optimization in Cartesian coordinates). This will likely take many more steps to converge but should be stable. - For optimizations in Cartesian coordinates the initial guess Hessian is constructed in internal coordinates and thus these optimizations should converge only slightly slower than those in internal coordinates. Nevertheless, if you observe a slow convergence behaviour, it may be a good idea to compute a Hessian initially (perhaps at a lower level of theory) and use `InHess read` in order to improve convergence. - At the beginning of a TS optimization more information on the curvature of the PES is needed than a model Hessian can give. The best choice is analytic Hessian, available for HF, DFT and MP2. In other cases (e.g. CAS-SCF), the numerical evaluation is necessary. Nevertheless you do not need to calculate the full Hessian when starting such a calculation. With ORCA we have good experience with approximations to the exact Hessian. Here it is recommended to either directly combine the TS optimization with the results of a relaxed surface scan or to use the Hybrid Hessian as the initial Hessian, depending on the nature of the TS mode. Note that these approximate Hessians do never replace the exact Hessian at the end of the optimization, which is always needed to verify the minimum or first order saddle point nature of the obtained structure. ::: (sec:optimization.hessian.typical)= ## Initial Hessian for Minimization The convergence of a geometry optimization crucially depends on the quality of the initial Hessian. In the simplest case it is taken as a unit matrix (in redundant internal coordinates we use 0.5 for bonds, 0.2 for angles and 0.1 for dihedrals and improper torsions). However, simple model force-fields like the ones proposed by Schlegel, Lindh, Swart or Almlöf are available and lead to much better convergence. The different guess Hessians can be set via the `InHess` option which can be either `unit`, `Almloef`, `Lindh`, `Swart` or `Schlegel` in redundant internal coordinates. Since version 2.5.30, these model force-fields (built up in internal coordinates) can also be used in optimizations in Cartesian coordinates. For minimizations we recommend the `Almloef` Hessian, which is the default for minimizations. The `Lindh` and `Schlegel` Hessian yield a similar convergence behaviour. From version 4.1?, there is also the option for the `Swart` model Hessian, which is less parametrized and should improve for weakly interacting and/or unusual structures. Of course the best Hessian is the exact one. `Read` may be used to input an exact Hessian or one that has been calculated at a lower level of theory (or a "faster" level of theory). From version 2.5.30 on this option is also available in redundant internal coordinates. But we have to point out that the use of the exact Hessian as initial one is only of little help, since in these cases the convergence is usually only slightly faster, while at the same time much more time is spent in the calculation of the initial Hessian. To sum it up: we advise to use one of the simple model force-fields for minimizations. (sec:optimization.coords.typical)= ## Coordinate Systems for Optimizations The coordinate system for the optimization can be chosen by the `coordsys` variable that can be set to `cartesian` or `redundant`within the `%geom` block. The default is the redundant internal coordinate system. If the optimization with `redundant` fails, you can still try `cartesian`. If the optimization is then carried out in Cartesian displacement coordinates with a simple model force-field Hessian, the convergence will be only slightly slower. With a unit matrix initial Hessian very slow convergence will result. A compound job `two_step_opt.inp` that first computes a semi-empirical Hessian to start from is shown below: ```{literalinclude} ../../orca_working_input/C05S02_122.inp :language: orca ``` :::{Tip} - For transition metal complexes MNDO, AM1 or PM3 Hessians are not available. You can use ZINDO/1 or NDDO/1 Hessians instead. They are of lower quality than MNDO, AM1 or PM3 for organic molecules but they are still far better than the standard unit matrix choice. - If the quality of the initial semi-empirical Hessian is not sufficient you may use a "quick" RI-DFT job (e.g. `BP def2-sv(p) defgrid1`) - In semi-empirical geometry optimizations on larger molecules or in general when the molecules become larger the redundant internal space may become large and the relaxation step may take a significant fraction of the total computing time. ::: For condensed molecular systems and folded molecules (e.g. a U-shaped carbon chain) atoms can get very close in space, while they are distant in terms of number of bonds connecting them. As damping of optimization steps in internal coordinates might not work well for these cases, convergence can slow down. ORCA's automatic internal coordinate generation takes care of this problem by assigning bonds to atom pairs that are close in real space, but distant in terms of number of bonds connecting them. ```orca %geom AddExtraBonds true # switch on/off assigning bonds to atom pairs that are # connected by more than bonds and are less # than Ang. apart (default true) AddExtraBonds_MaxLength 10 # cutoff for number of bonds connecting the two # atoms (default 10) AddExtraBonds_MaxDist 5 # cutoff for distance between two atoms (default 5 Ang.) end ``` For solid systems modeled as embedded solids the automatically generated set of internal coordinates might become very large, rendering the computing time spent in the optimization routine unnecessarily large. Usually, in such calculations the cartesian positions of outer atoms, coreless ECPs and point charges are constrained during the optimization - thus most of their internal coordinates are not needed. By requesting: ```orca %geom ReduceRedInts true end ``` only the required needed internal coordinates (of the constrained atoms) are generated. OBS: If the step in `redundant` fails badly and only Cartesian constrains are set (or no constrains), ORCA will fallback to a `cartesian` step automatically. This can be turned off by setting CARTFALLBACK to FALSE. (sec:optimization.constrained.typical)= ## Constrained Optimizations You can perform constrained optimizations which can, at times, be extremely helpful. This works as shown in the following example: ```{literalinclude} ../../orca_working_input/C05S02_123.inp :language: orca ``` ``` Constraining bond distances : { B N1 N2 value C } Constraining bond angles : { A N1 N2 N1 value C } Constraining dihedral angles : { D N1 N2 N3 N4 value C } Constraining cartesian coordinates : { C N1 C } ``` :::{Note} - Like for normal optimizations you can use numerical gradients (see {ref}`sec:optimization.numgrad.typical`.) for constrained optimizations. In this case the numerical gradient will be evaluated only for non-constrained coordinates, saving a lot of computational effort, if a large part of the structure is constrained. - "value" in the constraint input is optional. If you do not give a value, the present value in the structure is constrained. For cartesian constraints you can't give a value, but always the initial position is constrained. - It is recommended to use a value not too far away from your initial structure. - It is possible to constrain whole sets of coordinates: ::: ``` all bond lengths where N1 is involved : { B N1 * C} all bond lengths : { B * * C} all bond angles where N2 is the central atom: { A * N2 * C } all bond angles : { A * * * C } all dihedral angles with central bond N2-N3 : { D * N2 N3 * C } all dihedral angles : { D * * * * C } ``` - For Cartesian constraints lists of atoms can be defined: ``` a list of atoms (10 to 17) with Cartesian constraints : { C 10:17 C} ``` - Coordinates along a single Cartesian direction can be frozen as described in section {ref}`sec:input.coords.special`. - If there are only a few coordinates that have to be optimized you can use the `invertConstraints` option: ``` %geom Constraints { B 0 1 C } end invertConstraints true # only the C-O distance is optimized # does not affect Cartesian coordinates end ``` - In some cases it is advantageous to optimize only the positions of the hydrogen atoms and let the remaining molecule skeleton fixed: ``` %geom optimizehydrogens true end ``` :::{Note} - In the special case of a fragment optimization (see next point) the `optimizehydrogens` keyword does not fix the heteroatoms, but ensures that all hydrogen positions are relaxed. - In Cartesian optimization, only Cartesian constraints are allowed. ::: (sec:optimization.fragment.typical)= ## Constrained Optimizations for Molecular Clusters (Fragment Optimization) If you want to study systems, which consist of several molecules (e.g. the active site of a protein) with constraints, then you can either use cartesian constraints (see above) or use ORCA's fragment constraint option. ORCA allows the user to define fragments in the system. For each fragment one can then choose separately whether it should be optimized or constrained. Furthermore, it is possible to choose fragment pairs whose distance and orientation with respect to each other should be constrained. Here, the user can either define the atoms which make up the connection between the fragments, or the program chooses the atom pair automatically via a closest distance criterion. ORCA then chooses the respective constrained coordinates automatically. An example for this procedure is shown below. (fig:12)= ```{figure} ../../images/612.* ``` The coordinates are taken from a crystal structure \[PDB-code 2FRJ\]. In our gas phase model we choose only a small part of the protein, which is important for its spectroscopic properties. Our selection consists of a heme-group (fragment 1), important residues around the reaction site (lysine (fragment 2) and histidine (fragment 3)), an important water molecule (fragment 4), the NO-ligand (fragment 5) and part of a histidine (fragment 6) coordinated to the heme-iron. In this constrained optimization we want to maintain the position of the heteroatoms of the heme group. Since the protein backbone is missing, we have to constrain the orientation of lysine and histidine (fragments 2 and 3) side chains to the heme group. All other fragments (the ones which are directly bound to the heme-iron and the water molecule) are fully optimized internally and with respect to the other fragments. Since the crystal structure does not reliably resolve the hydrogen positions, we relax also the hydrogen positions of the heme group. ```{literalinclude} ../../orca_working_input/optimization-fragments.inp :language: orca ``` :::{Note} - You have to connect the fragments in such a way that the whole system is connected. - You can divide a molecule into several fragments. - Since the initial Hessian for the optimization is based upon the internal coordinates: Connect the fragments in a way that their real interaction is reflected. - This option can be combined with the definition of constraints, scan coordinates and the `optimizeHydrogens` option (but: its meaning in this context is different to its meaning in a normal optimization run, relatively straightforward see section {ref}`sec:geomopt.detailed`). - Can be helpful in the location of complicated transition states (with relaxed surface scans). ::: (sec:optimization.wallpot.typical)= ## Adding Arbitrary Wall Potentials For some applications, it might be interesting to add arbitrary wall potentials during the geometry optimization. For example, if you want to optimize an intermolecular complex and need that both structures stick together, without one flying away during the optimization, or when using microsolvation. In ORCA you can add three kinds of arbitrary "wall potentials": an ellipsoid or spherical of the form $$\begin{aligned} V = \left(\frac{|\mathbf{R} - \mathbf{O}|}{ radius}\right)^{30}\end{aligned}$$ or a rectangular box potential with 6 walls of the form $$\begin{aligned} V = e^{5(\mathbf{R}-wall) }\end{aligned}$$ These can be given in two ways: by explicitly defining the origin of the potential and its limits, e.g: ``` %GEOM ELLIPSEPOT 0,0,0,5,3,4 # the last are the a,b and c radii ``` or: ``` %GEOM SPHEREPOT 0,0,0,5 # the last is the radius ``` or: ``` %GEOM BOXPOT 0,0,0,4,-4,3,-3,6,-6 # maxx, minx, maxy, miny, maxz and minz last ``` where the first three numbers are the center and the last is the radius for the sphere (or a,b and c for the ellipsoid) and the max and min x,y and z dimensions of the box. All numbers should be given in Ångström. In case a single number is given instead, the walls will be automatically centered around the centroid of the molecule and that number will be added to the minimum sphere or box that is necessary to contain the molecule. For example: ``` %GEOM SPHEREPOT 2 ``` or: ``` %GEOM BOXPOT 2 ``` will build a minimum wall centered on the centroid that encloses the molecule and add 2 Ångström on top of it. Still on the sphere case, a negative number like ``` %GEOM SPHEREPOT -2 ``` will make the total radius of the sphere to be Ångström. OBS: This will apply to regular geometry optimizations, as well as to the Global Optimizer (GOAT). (sec:optimization.relaxedsurface.typical)= ## Relaxed Surface Scans A final thing that comes in really handy are relaxed surface scans, i.e. you can scan through one coordinate while all others are relaxed. It works as shown in the following example: ```{literalinclude} ../../orca_working_input/C05S02_128.inp :language: orca ``` In the example above the value of the bond length between C and O will be changed in 12 equidistant steps from 1.35 down to 1.10 Ångströms and at each point a constrained geometry optimization will be carried out. :::{Note} - If you want to perform a geometry optimization at a series of values with non-equidistant steps you can give this series in square brackets, \[ \]. The general syntax is as follows: ``` B N1 N2 = initial-value, final-value, NSteps or: B N1 N2 [value1 value2 value3 ... valueN] ``` - In addition to bond lengths you can also scan bond angles and dihedral angles: ``` B N1 N2 = ... # bond length A N1 N2 N3 = ... # bond angle D N1 N2 N3 N4 = ... # dihedral angle ``` ::: :::{Tip} - As in constrained geometry optimizations it is possible to start the relaxed surface scan with a different scan parameter than the value present in your molecule. But keep in mind that this value should not be too far away from your initial structure. ::: A more challenging example is shown below. Here, the H-atom abstraction step from $\text{CH}_4$ to OH-radical is computed with a relaxed surface scan (*vide supra*). The job was run as follows: ```{literalinclude} ../../orca_working_input/C05S02_131.inp :language: orca ``` It is obvious that the reaction is exothermic and passes through an early transition state in which the hydrogen jumps from the carbon to the oxygen. The structure at the maximum of the curve is probably a very good guess for the true transition state that might be located by a transition state finder. You will probably find that such relaxed surface scans are incredibly useful but also time consuming. Even the simple job shown below required several hundred single point and gradient evaluations (convergence problems appear for the SCF close to the transition state and for the geometry once the reaction partners actually dissociate -- this is to be expected). Yet, when you search for a transition state or you want to get insight into the shapes of the potential energy surfaces involved in a reaction it might be a good idea to use this feature. One possibility to ease the burden somewhat is to perform the relaxed surface scan with a "fast" method and a smaller basis set and then do single point calculations on all optimized geometries with a larger basis set and/or higher level of theory. At least you can hope that this should give a reasonable approximation to the desired surface at the higher level of theory -- this is the case if the geometries at the lower level are reasonable. (fig:13)= ```{figure} ../../images/613.* Relaxed surface scan for the H-atom abstraction from CH$_4$ by OH-radical (B3LYP/SV(P)). ``` ### Multidimensional Scans After several requests from our users ORCA now allows up to three coordinates to be scanned within one calculation. ```{literalinclude} ../../orca_working_input/C05S02_129.inp :language: orca ``` :::{Note} - For finding transition state structures of more complicated reaction paths ORCA now offers its very efficient NEB-TS implementation (see section {ref}`sec:neb.detailed`). - 2-dimensional or even 3-dimensional relaxed surface scans can become very expensive - e.g. requesting 10 steps per scan, ORCA has to carry out 1000 constrained optimizations for a 3-D scan. - The results can depend on the direction of the individual scans and the ordering of the scans. ::: Simultaneous multidimensional scans, in which all scan coordinates are changed at the same time, can be requested via the following keyword (which brings the cost of a multidimensional relaxed surface scan down to the cost of a single relaxed surface scan): ```orca %geom Scan B 0 1 = 3, 1, 15 end Scan B 1 2 = 1, 3, 15 end Simul_Scan true end ``` (sec:optimization.multiplexyz.typical)= ## Multiple XYZ File Scans A different type of scan is implemented in ORCA in conjunction with relaxed surface scans. Such scans produce a series of structures that are typically calculated using some ground state method. Afterwards one may want to do additional or different calculations along the generated pathway such as excited state calculations or special property calculations. In this instance, the "multiple XYZ scan" feature is useful. If you request reading from a XYZ file via: ```orca * xyzfile Charge Multiplicity FileName ``` this file could contain a number of structures. The format of the file is: ```orca Number of atoms M Comment line AtomName1 X Y Z AtomName2 X Y Z ... AtomNameM X Y Z > Number of atoms M Comment line AtomName1 X Y Z ... ``` Thus, the structures are simply of the standard XYZ format, separated by a "\>" sign. After the last structure no "\>" should be given but a blank line instead. The program then automatically recognizes that a multiple XYZ scan run is to be performed. Thus, single point calculations are performed on each structure in sequence and the results are collected at the end of the run in the same kind of trajectory`.dat` files as produced from trajectory calculations. In order to aid in using this feature, the relaxed surface scans produce a file called `MyJob.allxyz` that is of the correct format to be re-read in a subsequent run. (sec:optimization.transitionstate.typical)= ## Transition States (sec:optimization.transitionstate.intro.typical)= ### Introduction to Transition State Searches If you provide a good estimate for the structure of the transition state (TS) structure, then you can find the respective transition state with the following keywords (in this example we take the structure with highest energy of the above relaxed surface scan): ```{literalinclude} ../../orca_working_input/C05S02_134.inp :language: orca ``` :::{Note} - You need a good guess of the TS structure. Relaxed surface scans can help in almost all cases (see also example above). - For TS optimization (in contrast to geometry optimization) an exact Hessian, a Hybrid Hessian or a modification of selected second derivatives is necessary. - Analytic Hessian evaluation is available for HF and SCF methods, including the RI and RIJCOSX approximations and canonical MP2. - Check the eigenmodes of the optimized structure for the eigenmode with a single imaginary frequency. You can also visualize this eigenmode with `orca_pltvib` (section {ref}`sec:properties.vib.animation.typical`) or any other visualization program that reads ORCA output files. - If the Hessian is calculated during the TS optimization, it is stored as basename.001.hess, if it is recalculated several times, then the subsequently calculated Hessians are stored as basename.002.hess, basename.003.hess, ... - If you are using the Hybrid Hessian, then you have to check carefully at the beginning of the TS optimization (after the first three to five cycles) whether the algorithm is following the correct mode (see TIP below). If this is not the case you can use the same Hybrid Hessian again via the inhess read keyword and try to target a different mode (via the `TS_Mode` keyword, see below). ::: In the example above the TS mode is of local nature. In such a case you can directly combine the relaxed surface scan with the TS optimization with the ```orca ! ScanTS ``` command, as used in the following example: ```{literalinclude} ../../orca_working_input/C05S02_136.inp :language: orca ``` :::{Note} - The algorithm performs the relaxed surface scan, aborts the Scan after the maximum is surmounted, chooses the optimized structure with highest energy, calculates the second derivative of the scanned coordinate and finally performs a TS optimization. - If you do not want the scan to be aborted after the highest point has been reached but be carried out up to the last point, then you have to type: ```orca %geom fullScan true # do not abort the scan with !ScanTS end ``` ::: As transition state finder we implemented the quasi-Newton like Hessian mode following algorithm.{cite}`schlegel1987ab,schlegel1995modern,schlegel1998encyclopedia,eckert1997comp,horn1991comp,baker1986comp,hess2008chem,harvey1998theor,li2002theor` This algorithm maximizes the energy with respect to one (usually the lowest) eigenmode and minimizes with respect to the remaining $3N-7$(6) eigenmodes of the Hessian. :::{Tip} - You can check at an early stage if the optimization will lead to the "correct" transition state. After the first optimization step you find the following output for the redundant internal coordinates: ``` --------------------------------------------------------------------------- Redundant Internal Coordinates (Angstroem and degrees) Definition Value dE/dq Step New-Value comp.(TS mode) ---------------------------------------------------------------------------- 1. B(H 1,C 0) 1.2857 0.013136 0.0286 1.3143 0.58 2. B(H 2,C 0) 1.1002 0.014201 -0.0220 1.0782 3. B(H 3,C 0) 1.1010 0.014753 -0.0230 1.0779 4. B(H 4,C 0) 1.1008 0.014842 -0.0229 1.0779 5. B(O 5,H 1) 1.2442 -0.015421 -0.0488 1.1954 0.80 6. B(H 6,O 5) 0.9803 0.025828 -0.0289 0.9514 7. A(H 1,C 0,H 2) 107.38 -0.001418 -0.88 106.49 8. A(H 1,C 0,H 4) 105.48 -0.002209 -0.46 105.02 9. A(H 1,C 0,H 3) 103.35 -0.003406 0.08 103.43 10. A(H 2,C 0,H 4) 113.30 0.001833 0.35 113.65 11. A(H 3,C 0,H 4) 113.38 0.002116 0.26 113.64 12. A(H 2,C 0,H 3) 112.95 0.001923 0.45 113.40 13. A(C 0,H 1,O 5) 169.26 -0.002089 4.30 173.56 14. A(H 1,O 5,H 6) 100.84 0.003097 -1.41 99.43 15. D(O 5,H 1,C 0,H 2) 17.02 0.000135 0.24 17.26 16. D(O 5,H 1,C 0,H 4) -104.09 -0.000100 0.52 -103.57 17. D(O 5,H 1,C 0,H 3) 136.64 0.000004 0.39 137.03 18. D(H 6,O 5,H 1,C 0) 10.52 0.000078 -0.72 9.79 ---------------------------------------------------------------------------- ``` ::: Every Hessian eigenmode can be represented by a linear combination of the redundant internal coordinates. In the last column of this list the internal coordinates, that represent a big part of the mode which is followed uphill, are labelled. The numbers reflect their magnitude in the TS eigenvector (fraction of this internal coordinate in the linear combination of the eigenvector of the TS mode). Thus at this point you can already check whether your TS optimization is following the right mode (which is the case in our example, since we are interested in the abstraction of H1 from C0 by O5. - If you want the algorithm to follow a different mode than the one with lowest eigenvalue, you can either choose the number of the mode: ```orca %geom TS_Mode {M 1} # {M 1} mode with second lowest eigenvalue end # (default: {M 0}, mode with lowest eigenvalue) end ``` or you can give an internal coordinate that should be strongly involved in this mode: ```orca %geom TS_Mode {B 1 5} # bond between atoms 1 and 5, end # you can also choose an angle: {A N1 N2 N1} # or a dihedral: {D N1 N2 N3 N4} end ``` :::{Tip} - If you look for a TS of a breaking bond the respective internal coordinate might not be included in the list of redundant internal coordinates due to the bond distance being slightly too large, leading to slow or even no convergence at all. In order to prevent that behavior a region of atoms that are active in the TS search can be defined, consisting of e.g. the two atoms of the breaking bond. During the automatic generation of the internal coordinates the bond radii of these atoms (and their neighbouring atoms) are increased, making it more probable that breaking or forming bonds in the TS are detected as bonds. ```orca %geom TS_Active_Atoms { 1 2 3 } # atoms that are involved in TS, e.g. for proton end # transfer the proton, its acceptor and its donor TS_Active_Atoms_Factor 1.5 # factor by which the cutoff for bonds is increased for # the above defined atoms. # (Default 1.5, i.e. increased by 50%) end ``` ::: (sec:optimization.transitionstate.hessian.typical)= ### Hessians for Transition State Calculations For transition state (TS) optimization a simple initial Hessian, which is used for minimization, is not sufficient. In a TS optimization we are looking for a first order saddle point, and thus for a point on the PES where the curvature is negative in the direction of the TS mode (the TS mode is also called transition state vector, the only eigenvector of the Hessian at the TS geometry with a negative eigenvalue). Starting from an initial guess structure the algorithm used in the ORCA TS optimization has to climb uphill with respect to the TS mode, which means that the starting structure has to be near the TS and the initial Hessian has to account for the negative curvature of the PES at that point. The simple force-field Hessians cannot account for this, since they only know harmonic potentials and thus positive curvature. The most straightforward option in this case would be (after having looked for a promising initial guess structure with the help of a relaxed surface scan) to calculate the exact Hessian before starting the TS optimization. With this Hessian (depending on the quality of the initial guess structure) we know the TS eigenvector with its negative eigenvalue and we have also calculated the exact force constants for all other eigenmodes (which should have positive force constants). For the HF, DFT methods and MP2, the analytic Hessian evaluation is available and is the best choice, for details see section Frequencies ({ref}`sec:frequencies.typical`). When only the gradients are available (most notably the CASSCF), the numerical calculation of the exact Hessian is very time consuming, and one could ask if it is really necessary to calculate the full exact Hessian since the only special thing (compared to the simple force-field Hessians) that we need is the TS mode with a negative eigenvalue. Here ORCA provides two different possibilities to speed up the Hessian calculation, depending on the nature of the TS mode: the Hybrid Hessian and the calculation of the Hessian value of an internal coordinate. For both possibilities the initial Hessian is based on a force-field Hessian and only parts of it are calculated exactly. If the TS mode is of very local nature, which would be the case when e.g. cleaving or forming a bond, then the exactly calculated part of the Hessian can be the second derivative of only one internal coordinate, the one which is supposed to make up the TS mode (the formed or cleaved bond). If the TS mode is more complicated and more delocalized, as e.g. in a concerted proton transfer reaction, then the hybrid Hessian, a Hessian matrix in which the numerical second derivatives are calculated only for those atoms, which are involved in the TS mode (for more details, see section {ref}`sec:geomopt.detailed`), should be sufficient. If you are dealing with more complicated cases where these two approaches do not succeed, then you still have the possibility to start the TS optimization with a full exact Hessian. Numerical Frequency calculations are quite expensive. You can first calculate the Hessian at a lower level of theory or with a smaller basis set and use this Hessian as input for a subsequent TS optimization: ```orca %geom inhess Read # this command comes with: InHessName "yourHessian.hess" # filename of Hessian input file end ``` Another possibility to save computational time is to calculate exact Hessian values only for those atoms which are crucial for the TS optimization and to use approximate Hessian values for the rest. This option is very useful for big systems, where only a small part of the molecule changes its geometry during the transition and hence the information of the full exact Hessian is not necessary. With this option the coupling of the selected atoms are calculated exactly and the remaining Hessian matrix is filled up with a model initial Hessian: ```orca %geom Calc_Hess true Hybrid_Hess {0 1 5 6} end # calculates a Hybrid Hessian with # exact calculation for atoms 0, 1, 5 and 6 end ``` For some molecules the PES near the TS can be very far from ideal for a Newton-Raphson step. In such a case ORCA can recalculate the Hessian after a number of steps: ```orca %geom Recalc_Hess 5 # calculate the Hessian at the beginning # and recalculate it after 5,10,15,... steps end ``` Another solution in that case is to switch on the trust radius update, which reduces the step size if the Newton-Raphson steps behave unexpected and ensures bigger step size if the PES seems to be quite quadratic: ```orca %geom Trust 0.3 # Trust <0 - use fixed trust radius (default: -0.3 au) # Trust >0 - use trust radius update, i.e. 0.3 means: # start with trust radius 0.3 and use trust radius update end ``` (sec:optimization.transitionstate.coords.typical)= ### Special Coordinates for Transition State Optimizations - If you look for a TS of a breaking bond the respective internal coordinate might not be included in the list of redundant internal coordinates (but this would accelerate the convergence). In such a case (and of course in others) you can add coordinates to or remove them from the set of autogenerated redundant internal coordinates (alternatively check the TS_Active_Atoms keyword): ```orca # add ( A ) or remove ( R ) internal coordinates %geom modify_internal { B 10 0 A } # add a bond between atoms 0 and 10 { A 8 9 10 R } # remove the angle defined # by atoms 8, 9 and 10 { D 7 8 9 10 R } # remove the dihedral angle defined end # by atoms 7, 8, 9 and 10 end ``` (sec:optimization.mecp.typical)= ## MECP Optimization There are reactions where the analysis of only one spin state of a system is not sufficient, but where the reactivity is determined by two or more different spin states (Two- or Multi-state reactivity). The analysis of such reactions reveals that the different PESs cross each other while moving from one stationary point to the other. In such a case you might want to use the ORCA optimizer to locate the point of lowest energy of the crossing surfaces (called the minimum energy crossing point, MECP). As an example for such an analysis we show the MECP optimization of the quartet and sextet state of \[FeO\]$^{+}$. ```{literalinclude} ../../orca_working_input/C05S02_145.inp :language: orca ``` - For further options for the MECP calculation, see section {ref}`sec:geomopt.mecp.detailed`. :::{tip} You can often use a minimum or TS structure of one of the two spin states as initial guess for your MECP-optimization. If this doesn't work, you might try a scan to get a better initial guess. ::: The results of the MECP optimization are given in the following output. The distance where both surfaces cross is at 1.994 Å. In this simple example there is only one degree of freedom and we can also locate the MECP via a parameter scan. The results of the scan are given in {numref}`fig:14` for comparison. Here we see that the crossing occurs at a Fe-O-distance of around 2 Å. For systems with more than two atoms a scan is not sufficient any more and you have to use the MECP optimization. ``` ***********************HURRAY******************** *** THE OPTIMIZATION HAS CONVERGED *** ************************************************* ------------------------------------------------------------------- Redundant Internal Coordinates --- Optimized Parameters --- (Angstroem and degrees) Definition OldVal dE/dq Step FinalVal ------------------------------------------------------------------- 1. B(O 1,Fe 0) 1.9939 -0.000001 0.0000 1.9939 ------------------------------------------------------------------- ******************************************************* *** FINAL ENERGY EVALUATION AT THE STATIONARY POINT *** *** (AFTER 8 CYCLES) *** ******************************************************* ------------------------------------- ---------------- Energy difference between both states -0.000000061 ------------------------------------- ---------------- ``` (fig:14)= ```{figure} ../../images/614.* Parameter scan for the quartet and sextet state of [FeO]$^+$ (B3LYP/SV(P)). ``` A more realistic example with more than one degree of freedom is the MECP optimization of a structure along the reaction path of the CH$_{3}$O $\leftrightarrow$ CH$_{2}$OH isomerization. ```{literalinclude} ../../orca_working_input/C05S02_146.inp :language: orca ``` :::{Note} - To verify that a stationary point in a MECP optimization is a minimum, you have to use an adapted frequency analysis, called by `SurfCrossNumFreq` (see section {ref}`sec:geomopt.mecp.detailed`). ::: (sec:optimization.conical.typical)= ## Conical Intersection Optimization OBS.: It is currently only available using TD-DFT, will be expanded in future versions. More details about the specific options on {ref}`sec:geomopt.ciopt.detailed`. A conical intersection (CI) is a complicated 3N-8 dimensional space, where two potential energy surfaces cross and the energy difference between these two states is zero. Inside this so-called "seam-space" minima and transition states can exist. Locating these minima is essential to understand photo-chemical processes, that are governed by non-adiabatic events, as e.g. photoisomerization, photostability - similar to locating transition states for chemical reactions. As an example for such an analysis we show the conical intersection optimization of the ground and first excited state of singlet ethylene. :::{Note} Even though locating the CI of a TD-DFT excited state and the reference state is supported, it is not the recommended way of finding the ground state-excited state CI, because such CIs are not described properly by TD-DFT (in particular, TD-DFT even predicts the wrong dimensionality for the intersection space). Instead, it is advised to use SF-TD-DFT for this purpose, e.g. use the $T_1$ state as the reference state, and calculate both the $S_0$ and $S_1$ states as excited states. (*vide infra*) ::: ```orca !B3LYP DEF2-SVP CI-OPT %TDDFT IROOT 1 END * xyz 0 1 C 0.595560237 -0.010483480 -0.000284187 C -0.831313750 0.167231832 0.001482505 H -1.381857976 0.227877089 0.963419721 H 1.265119434 0.874806815 0.006897459 H -1.382258208 0.243775568 -0.959090898 H 1.027489724 -1.032962768 -0.008829646 * ``` :::{Tip} You can often use a structure between the optimized structures of both states for your CI-optimization. If this doesn't work, you might try a scan to get a better initial guess. ::: The results of the CI-optimization are given in the following output. The energy difference between the ground and excited state is printed as E diff. (CI), being reasonabley close for a conical intersection. For a description of the calculation of the non-adiabatic couplings at this geometry, see section {ref}`sec:excitedstates.numnacme.typical`. ``` .--------------------. ----------------------|Geometry convergence|------------------------- Item value Tolerance Converged --------------------------------------------------------------------- Energy change 0.0000164283 0.0000050000 NO E diff. (CI) 0.0000025162 0.0001000000 YES RMS gradient 0.0000068173 0.0001000000 YES MAX gradient 0.0000136891 0.0003000000 YES RMS step 0.0000358228 0.0020000000 YES MAX step 0.0000821130 0.0040000000 YES ........................................................ Max(Bonds) 0.0000 Max(Angles) 0.00 Max(Dihed) 0.00 Max(Improp) 0.00 --------------------------------------------------------------------- Everything but the energy has converged. However, the energy appears to be close enough to convergence to make sure that the final evaluation at the new geometry represents the equilibrium energy. Convergence will therefore be signaled now ***********************HURRAY******************** *** THE OPTIMIZATION HAS CONVERGED *** ************************************************* --------------------------------------------------------------------------- Redundant Internal Coordinates --- Optimized Parameters --- (Angstroem and degrees) Definition OldVal dE/dq Step FinalVal ---------------------------------------------------------------------------- 1. B(C 1,C 0) 1.3254 -0.000005 -0.0000 1.3254 2. B(H 2,C 1) 1.1270 0.000004 -0.0000 1.1270 3. B(H 3,C 0) 1.1271 -0.000002 0.0000 1.1271 4. B(H 4,C 1) 1.1271 0.000000 -0.0000 1.1271 5. B(H 5,C 0) 1.1271 -0.000002 0.0000 1.1271 6. A(H 3,C 0,H 5) 106.00 0.000001 -0.00 106.00 7. A(C 1,C 0,H 5) 126.97 -0.000013 0.00 126.97 8. A(C 1,C 0,H 3) 127.03 0.000011 -0.00 127.03 9. A(C 0,C 1,H 4) 127.03 0.000013 -0.00 127.03 10. A(H 2,C 1,H 4) 106.01 0.000001 -0.00 106.01 11. A(C 0,C 1,H 2) 126.96 -0.000014 0.00 126.96 12. D(H 2,C 1,C 0,H 5) 73.60 0.000000 -0.00 73.59 13. D(H 4,C 1,C 0,H 3) 72.78 -0.000001 0.00 72.78 14. D(H 4,C 1,C 0,H 5) -106.81 0.000001 -0.00 -106.82 15. D(H 2,C 1,C 0,H 3) -106.82 -0.000002 0.00 -106.81 ---------------------------------------------------------------------------- ******************************************************* *** FINAL ENERGY EVALUATION AT THE STATIONARY POINT *** *** (AFTER 12 CYCLES) *** ******************************************************* ``` **CI minima between excited states** In an analogous way, the conical intersection minima between two excited states can be requested by selection both an IROOT and a JROOT, shown below. ```orca !B3LYP DEF2-SVP CI-OPT %TDDFT IROOT 2 JROOT 1 #IROOTMULT TRIPLET would search in the triplet PES #SF TRUE would search for the S0-S1 CI from a T1 reference, using SF-TD-DFT # (but remember to set the multiplicity as 3 instead of 1) END * xyz 0 1 C 0.595560237 -0.010483480 -0.000284187 C -0.831313750 0.167231832 0.001482505 H -1.381857976 0.227877089 0.963419721 H 1.265119434 0.874806815 0.006897459 H -1.382258208 0.243775568 -0.959090898 H 1.027489724 -1.032962768 -0.008829646 * ``` (sec:optimization.mechanochemistry.typical)= ## Constant External Force - Mechanochemistry Constant external force can be applied on the molecule within the EFEI formalism{cite}`marx:ChemRev2012` by pulling on the two defined atoms. To apply the external force, use the POTENTIALS in the geom block. The potential type is C for Constant force, indexes of two atoms (zero-based) and the value of force in nN. ```{literalinclude} ../../orca_working_input/C05S02_147.inp :language: orca ``` The results are seen in the output of the SCF procedure, where the total energy already contains the force term. ``` ---------------- TOTAL SCF ENERGY ---------------- Total Energy : -150.89704913 Eh -4106.11746 eV Components: Nuclear Repulsion : 36.90074715 Eh 1004.12038 eV External potential : -0.25613618 Eh -6.96982 eV Electronic Energy : -187.54166010 Eh -5103.26802 eV ``` (sec:optimization.irc.typical)= ## Intrinsic Reaction Coordinate The Intrinsic Reaction Coordinate (IRC) is a special form of a minimum energy path, connecting a transition state (TS) with its downhill-nearest intermediates. A method determining the IRC is thus useful to determine whether a transition state is directly connected to a given reactant and/or a product. ORCA features its own implementation of Morokuma and coworkers' popular method.{cite}`Morokuma:JCP1977` The IRC method can be simply invoked by adding the IRC keyword as in the following example. ```{literalinclude} ../../orca_working_input/C05S02_148.inp :language: orca ``` For more information and further options see section {ref}`sec:frequencies.irc.detailed`. :::{Note} - The same method and basis set as used for optimization and frequency calculation should be used for the IRC run. - The IRC keyword can be requested without, but also together with OptTS, ScanTS, NEB-TS, AnFreq and NumFreq keywords. - In its default settings the IRC code checks whether a Hessian was computed before the IRC run. If that is not the case, and if no Hessian is defined via the %irc block, a new Hessian is computed at the beginning of the IRC run. - A final trajectory (\_IRC_Full_trj.xyz) is generated which contains both directions, forward and backward, by starting from one endpoint and going to the other endpoint, visualizing the entire IRC. Forward (\_IRC_F\_trj.xyz and \_IRC_F.xyz) and backward (\_IRC_B\_trj.xyz and \_IRC_B.xyz) trajectories and xyz files contain the IRC and the last geometry of that respective run. ::: (sec:optimization.internalcoords.typical)= ## Printing Hessian in Internal Coordinates When a Hessian is available, it can be printed out in redundant internal coordinates as in the following example: ```{literalinclude} ../../orca_working_input/C05S02_151.inp :language: orca ``` :::{Note} - The Hessian in internal coordinates is (for the input ```printHess.inp```) stored in the file ```printHess_internal.hess```. - The corresponding lists of redundant internals is stored in ```printHess.opt```. - Although the ```!Opt``` keyword is necessary, an optimization is not carried out. ORCA exits after storing the Hessian in internal coordinates. ::: (sec:optimization.modelHess.typical)= ## Using model Hessian from previous calculations If you had a geometry optimization interrupted, or for some reason want to use the model Hessian updated from a previous calculation, you can do that by passing a ```basename.opt``` file, a ```basename.carthess``` file or the initial Hessian on a new calculation. ```orca %GEOM InHess READ InHessName "basename.carthess" END ``` (sec:optimization.lbfgs.geometry.typical)= ## Geometry Optimizations using the L-BFGS optimizer Optimizations using the L-BFGS optimizer are done in Cartesian coordinates. They can be invoked quite simple as in the following example: ```orca ! L-Opt ! MM %mm ORCAFFFILENAME "CHMH.ORCAFF.prms" end *pdbfile 0 1 CHMH.pdb ``` Using this optimizer systems with 100s of thousands of atoms can be optimized. Of course, the energy and gradient calculations should not become the bottleneck for such calculations, thus MM or QM/MM methods should be used for such large systems. The default maximum number of iterations is 200, and can be increased as follows: ```orca ! L-Opt %geom maxIter 500 # default 200 end *pdbfile 0 1 CHMH.pdb ``` Only the hydrogen positions can be optimized with the following command: ```orca ! L-OptH ``` But also other elements can be exclusively optimized with the following command: ```orca ! L-OptH %geom OptElement F # optimize fluorine only when L-OptH is invoked. # Does not work with the regular optimizer. end ``` When fragments are defined for the system, each fragment can be optimized differently (similar to the fragment optimization described above). The following options are available: FixFrags : Freeze the coordinates of all atoms of the specified fragments. RelaxHFrags : Relax the hydrogen atoms of the specified fragments. Default for all atoms if !L-OptH is defined. RelaxFrags : Relax all atoms of the specified fragments. Default for all atoms if !L-Opt is defined. RigidFrags : Treat each specified fragment as a rigid body, but relax the position and orientation of these rigid bodies. :::{Note} - The fragments have to be defined after the coordinate input. ::: A more complex example is depicted in the following: ```orca ! L-OptH %mm ORCAFFFILENAME "CHMH.ORCAFF.prms" end *pdbfile 0 1 CHMH.pdb %geom Frags 2 {8168:8614} end # First the fragments need to be defined 3 {8615:8699} end # Note that all other atoms belong to 4 {8700:8772} end # fragment 1 by default 5 {8773:8791} end # RelaxFrags {2} end # Fragment 2 is fully relaxed RigidFrags {3 4 5} end # Fragments 3, 4 and 5 are treated as rigid bodies each. end ``` (sec:optimization.neb.typical)= ## Nudged Elastic Band Method The Nudged Elastic Band (NEB) method is used to find a minimum energy path (MEP) connecting given reactant and product state minima on the energy surface. An initial path is generated and represented by a discrete set of configurations of the atoms, referred to as images of the system. The number of images is specified by the user and has to be large enough to obtain sufficient resolution of the path. The implementation in ORCA is described in detail in the article by Ásgeirsson et. al.{cite}`Asgeirsson2021` and in section {ref}`sec:neb.detailed` along with the input options. The most common use of the NEB method is to find the highest energy saddle point on the potential energy surface specifying the transition state for a given initial and final state. Rigorous convergence to a first order saddle point can be obtained with the climbing image NEB (CI-NEB), where the highest energy image is pushed uphill in energy along the tangent to the path while relaxing downhill in orthogonal directions. Another method for finding a first order saddle point is the NEB-TS which uses the CI-NEB method with a loose tolerance to begin with and then switches over to the OptTS method to converge on the saddle point. This combination can be a good choice for calculations of complex reactions where the ScanTS method fails or where 2D relaxed surface scans are necessary to find a good initial guess structure for the OptTS method. The zoomNEB variants are a good choice in case of very complex transition states with long tails. For more and detailed information on the various NEB variants implemented in ORCA please consult section {ref}`sec:neb.detailed`. In their simplest form NEB, NEB-CI and NEB-TS only require the reactant and product state configurations (one via the xyz block, and the other one via the keyword neb_end_xyzfile): ```orca !NEB-TS # or !NEB or !NEB-CI or !ZOOM-NEB-TS or !ZOOM-NEB-CI # or !Fast-NEB-TS (corresponds to IDPP-TS defined in the NEB-TS manuscript) # or !Loose-NEB-TS (corresponds to default NEB-TS in the NEB-TS manuscript) %neb neb_end_xyzfile "final.xyz" end ``` Below is an example of an NEB-TS run involving an intramolecular proton transfer within acetic acid. The simplest input is ```orca !XTB NEB-TS %neb neb_end_xyzfile "final.xyz" end *xyz 0 1 C 0.416168 0.038758 -0.014077 C 0.041816 0.011798 1.439610 O 1.524458 0.176600 -0.453888 O -0.654209 -0.127881 -0.803857 H -0.391037 -0.126036 -1.737478 H -0.913438 0.507022 1.585301 H -0.057787 -1.026455 1.750845 H 0.819515 0.485425 2.030252 * ``` Where the final.xyz structure contains the corresponding structure with the proton on the other oxygen. The initial path is reasonable and the CI calculation can be switched on after five NEB iterations. ``` Starting iterations: Optim. Iteration HEI E(HEI)-E(0) max(|Fp|) RMS(Fp) dS Switch-on CI threshold 0.020000 LBFGS 0 4 0.081017 0.073897 0.018915 3.2882 LBFGS 1 5 0.070244 0.056668 0.013913 3.2770 LBFGS 2 5 0.062934 0.038972 0.008763 3.3376 LBFGS 3 5 0.057358 0.032076 0.006535 3.3950 LBFGS 4 4 0.053260 0.019015 0.003599 3.4826 Image 4 will be converted to a climbing image in the next iteration (max(|Fp|) < 0.0200) Optim. Iteration CI E(CI)-E(0) max(|Fp|) RMS(Fp) dS max(|FCI|) RMS(FCI) Convergence thresholds 0.020000 0.010000 0.002000 0.001000 ``` The CI run converges after another couple of iterations: ``` *********************H U R R A Y********************* *** THE NEB OPTIMIZATION HAS CONVERGED *** ***************************************************** ``` Subsequently a summary of the MEP is printed: ``` --------------------------------------------------------------- PATH SUMMARY --------------------------------------------------------------- All forces in Eh/Bohr. Image Dist.(Ang.) E(Eh) dE(kcal/mol) max(|Fp|) RMS(Fp) 0 0.000 -14.45993 0.00 0.00011 0.00004 1 0.426 -14.44891 6.91 0.00092 0.00033 2 0.652 -14.42864 19.63 0.00084 0.00038 3 0.805 -14.41132 30.50 0.00075 0.00027 4 0.932 -14.40562 34.08 0.00057 0.00018 <= CI 5 1.044 -14.41047 31.03 0.00057 0.00024 6 1.153 -14.42200 23.80 0.00103 0.00034 7 1.280 -14.43666 14.60 0.00098 0.00037 8 1.476 -14.45106 5.56 0.00106 0.00033 9 1.869 -14.45988 0.03 0.00013 0.00006 ``` Additionally, detailed information on the highest energy image (or the CI) is printed: ``` --------------------------------------------------------------- INFORMATION ABOUT SADDLE POINT --------------------------------------------------------------- Climbing image .... 4 Energy .... -14.40561577 Eh Max. abs. force .... 9.5976e-04 Eh/Bohr ----------------------------------------- SADDLE POINT (ANGSTROEM) ----------------------------------------- C 0.040867 0.007347 -0.497635 C -0.075595 0.017879 0.979075 O 1.122340 0.126074 -1.145534 O -0.928470 -0.137946 -1.298318 H 0.165808 -0.021676 -2.055704 H -0.996979 0.514720 1.271668 H -0.116377 -1.013504 1.327873 H 0.788406 0.507105 1.418575 ----------------------------------------- FORCES (Eh/Bohr) ----------------------------------------- C -0.000646 -0.000111 0.000086 ... ----------------------------------------- UNIT TANGENT ----------------------------------------- C -0.246569 -0.031821 -0.019359 ... => Unit tangent is an approximation to the TS mode at the saddle point Next a TS optimization is performed on the CI from the NEB run. ``` Finally, a TS optimization is started, after which the MEP information is updated by including the TS structure: ``` --------------------------------------------------------------- PATH SUMMARY FOR NEB-TS --------------------------------------------------------------- All forces in Eh/Bohr. Global forces for TS. Image E(Eh) dE(kcal/mol) max(|Fp|) RMS(Fp) 0 -14.45993 0.00 0.00011 0.00004 1 -14.44891 6.91 0.00092 0.00033 2 -14.42864 19.63 0.00084 0.00038 3 -14.41132 30.50 0.00075 0.00027 4 -14.40562 34.08 0.00057 0.00018 <= CI TS -14.40562 34.08 0.00033 0.00013 <= TS 5 -14.41047 31.03 0.00057 0.00024 6 -14.42200 23.80 0.00103 0.00034 7 -14.43666 14.60 0.00098 0.00037 8 -14.45106 5.56 0.00106 0.00033 9 -14.45988 0.03 0.00013 0.00006 ``` Note that here both TS and CI are printed for comparison.