Geometry optimizations

Geometry optimizations are usually performed at the DFT level as DFT has a very favorable cost-accuracy ratio and because analytical gradients are available for most functionals. It's typically a considerable effort to get much better geometries than those obtained from DFT calculations (i.e. via electron correlation WFT methods) and in practice only possible for small systems. The PBE0 functional with a dispersion correction and a triple-zeta basis is probably the most accurate all-round method for geometry optimizations of molecules from the whole periodic table (note that open-shell 3d transition metal systems might be an exception to this).

• RI-MP2 optimizations (MP2 or SCS-MP2) is another option for main-group elements (MP2 is a very robust method for organic molecules) while for transition-metal chemistry, DFT optimizations are generally preferred.

• Sometimes NDDO-type semiempirical methods can be useful for pre-optimizations of closed-shell organic molecules (typically not for transition metal compounds) but these methods are not as robust as DFT methods.

• Alternative methods for cheap (but often accurate) geometry optimizations (or at least pre-optimizations) come from the Grimme group:

  • the HF-3c method which is a minimal-basis HF calculation, corrected for dispersion, BSSE and basis set incompleteness. It's even cheaper than a DFT calculation and more robust than a standard AM1/PM3 semiempirical calculation (but not as fast). See more here.

  • DFT-3c methods such as: r2SCAN-3c, B97-3c and PBEh-3c which

  • The GFN-xTB tightbinding methods (available in ORCA 4.2). See GFN1-xTB paper , GFN2-xTB paper and a study on large TM complexes with GFN-xTB methods.

• ORCA, by default uses a Quasi-Newton optimiser using the BFGS update (Powell and Bofill also available) and the optimization is carried out in redundant internal coordinates. The initial Hessian is by default an approximate model Hessian by Almlöf (other possibilities are model Hessians by Lindh and Schlegel, a diagonal Hessian or a previously calculated Hessian). TightSCF convergence criteria is enforced by default in ORCA for geometry optimizations (reduces numerical noise in the gradient). The default convergence criteria (NormalOpt) are: TolE=5e-6, TolRMSG=1e-4, TolMaxG=3e-4, TolRMSD=2e-3, TolMaxD=4e-3. Note that in ORCA 4.0, if a geometry optimization does not converge within the Max. number of cycles, ORCA will not go on to do a frequency calculation.

General recommendations:

• Choose a good starting geometry. If you are building your starting coordinates from scratch it makes sense to look them over carefully to check for reasonable bond lengths etc. This will speed up your optimization. Programs like Chemcraft can sometimes be used to "clean up" a structure so that the bond lengths/angles are more reasonable. Doing a preliminary optimization with a cheap semi-empirical method such as GFN-xTB can also work (runs in seconds to minutes for small to medium molecules) or perhaps HF-3c/B97-3c.

• Redundant internal coordinates are automatically used (regardless of coordinate input) for efficient geometry optimizations. In rare cases, the internal coordinates run into trouble and ORCA either stops or the optimization takes a very long time. In such cases it is possible to switch to optimization in Cartesian coordinates instead (COPT) which is slower but will usually converge in the end. Should that also lead to trouble, other options can be tried (ZOPT, GDIIS-COPT, GDIIS-ZOPT) although calculating the full exact Hessian to help the optimizer along, is probably more likely to succeed (see below).

• Use GGA DFT functionals if they are accurate enough (depends on your system), with the RI-J approximation (default) as that is often the fastest useful optimization one can do. Use of the RI-J approximation leads to minimal geometrical errors. Often the slightly higher accuracy from hybrid functionals is not worth the effort.

• Double-zeta basis sets are often good enough for a decent optimization of organic molecules. However, for transition metal systems, a triple-zeta basis set on the metal is usually a minimum requirement and one should certainly aim for triple-zeta quality geometries in general for all elements (for really bulky organic ligands, this may not be necessary). Geometry optimizations with a quadruple-zeta basis set is almost never required (usually practically identical to triple-zeta geometries).

• Hybrid DFT optimizations can be sped up a lot by use of the RIJCOSX or RIJK approximation. RIJCOSX/RIJK errors for geometries tend to be minimal.

• Use a dispersion correction, like Grimme’s DFT-D3(BJ) method which includes dispersion in DFT calculations free of charge (D3BJ keyword). Grimme's DFT-D3(BJ) method is very robust, is crucial for reliably describing noncovalent interaction dominated potential energy surfaces and dispersion becomes more and more important as the molecule grows larger in size.

• The XC grid in DFT calculations is an important aspect of DFT optimizations. For systems with heavy atoms, the default grid settings are no longer applicable and can result in considerable geometrical errors (See Numerical accuracy). A too small grid can sometimes even slow down geometry optimizations or even prevent one from ever reaching the minimum. In my experience (RB) using a larger than default grid often speeds up the geometry optimization as fewer optimization steps are needed so it is often a good idea.

• Rarely needed, but occasionally one wants a tightly optimized geometry (gradient even closer to zero than the defaults). Use the TIGHTOPT keyword for that: set TolE=1e-6, TolRMSG=3e-5, TolMaxG=1e-4, TolRMSD=6e-4, TolMaxD=1e-3. You can also set the thresholds manually, see manual). Increasing the SCF convergence and DFT grids might then also be a good idea if the intent is to reach maximum accuracy.

Optimization using the BP86 functional, a double-zeta basis, the RI-J approximation and DFT-D3 dispersion correction

This is pretty much the fastest way of running a geometry optimization that gives fairly reliable geometries for a good portion of the periodic table. Increasing the basis set to def2-TZVP (after converging the def2-SVP calculation) is then an obvious next step if a more accurate geometry is required, and/or possibly switching to a hybrid functional (depends on the system). If the system contains a metal, it is wise to put a triple-zeta basis set on the metal straight away ( e.g. %basis newgto Fe "def2-TZVP" end end ).

! RI BP86 def2-SVP def2/J D3BJ TIGHTSCF Opt

Nowadays one might also consider the use of the r2SCAN functional, the r2SCAN-3c composite method etc. as alternative choices for a quick optimization.

Using the exact Hessian in optimizations

For tricky potential energy surfaces it may be desirable to calculate the exact Hessian in the beginning. The Hessian information is used to make a better guess for the next optimization step. For really tricky potential energy surfaces (a very flat PES for example) one might even want to calculate the exact Hessian regularly (in every step for a really tricky PES). Calculation of the exact Hessian may be quite expensive and so it will depend on the molecular system whether requesting this is optimal or not. Depending on the availability of the analytical Hessian for a method, the analytical or numerical one can be chosen.

For methods with an analytical Hessian available (Hartree-Fock method, GGA and hybrid-GGA DFT functionals):

! BP86 def2-SVP def2/J Opt
%geom
Calc_Hess true # Request an exact Hessian (here analytical) in the first optimization step
Recalc_Hess 10 # Recalculate the exact Hessian every 10 steps.
end

For methods with only a numerical Hessian available (meta-GGA DFT functionals, semiempirical methods, ab initio correlation methods etc.) the NumHess option needs to be turned on:

! TPSS def2-SVP def2/J Opt

%geom
Calc_Hess true # Request an exact Hessian in the first optimization step
NumHess true # Request the numerical exact Hessian (only one available for this method).
Recalc_Hess 10 # Recalculate the exact Hessian every 10 steps.
end

Reading in an old Hessian

If you have a previously calculated Hessian (from a CalcHess job or a previous frequency calculation) you can read it in to your geometry optimization job to speed up the geometry optimization (should always be better than the default Almlöf model Hessian). Even an Hessian calculated at a different level of theory is a good option. Reading in an old Hessian like this is especially useful for OptTS jobs (see later) where a reliable Hessian is necessary for locating saddlepoints.

! TPSS def2-SVP def2/J Opt
%geom
inhess Read
InHessName "test.hess"
end

Using hybrid Hessians

Hybrid Hessians are an alternative to calculating the exact full Hessian. By selecting only a few atoms of the systems for which to calculate the Hessian one can get a much cheaper, yet often useful Hessian to aid geometry optimizations or saddle point searches (see below).

! TPSS def2-SVP def2/J Opt

%geom
Calc_Hess true
Hybrid_Hess [0 1 2 3] # Requests a Hessian calculation for only the atoms 0, 1, 2 and 3.
end
end

Using hybrid Hessians does not always work and one should e.g check during the OptTS job whether the correct eigenmode is being followed (see Redundant Internal Coordinates table after 1 opt cycle). If an incorrect mode is being followed, then one can stop the job, read in the previous hybrid Hessian again (see above) and select a different mode using the TS_Mode keyword

! TPSS def2-SVP def2/J Opt
%geom
TS_Mode {M 1} # Choosing mode number 1 instead of mode 0 (default) to follow.
end
end

Optimization to a saddlepoint

If one already has a decent guess for the saddlepoint ("transition state"), e.g. from a relaxed surface scan or NEB calculation (see below), then the OptTS keyword (eigenvalue-following TS finding algorithm) can be used. Using exact Hessians are often useful here. One can also read in previously calculated Hessians (see above) or hybrid Hessians (see above).

See Tutorial: Transition state Optimization for a more detailed discussion on finding saddlepoints.

See the Tutorial: NEB calculations for information on how to find saddlepoints using the NEB method instead.

See the IRC tutorial on confirming saddlepoints by calculating the intrinsic reaction coordinate.

! BP86 def2-SVP def2/J OptTS

%geom
Calc_Hess true # Request an exact analytical Hessian in the first optimization step.
Recalc_Hess 5 # Recalculate the exact Hessian every 5 steps.
end

By default, the OptTS algorithm will follow the mode with the most negative eigenvalue. If your geometry has been chosen wisely (e.g. the highest energy structure from a relaxed scan) then typically the saddlepoint mode will be the one with the most negative eigenvalue. In some cases, this is not the case (check by visualizing the modes) and the mode of interest is instead the second one. In this case it is possible to tell OptTS to follow that one instead.

Here choosing to follow the second lowest eigenmode (1) instead of the first eigenmode (0).

%geom
TS_Mode {M 1}
end
end

Constrained optimizations

Keep in mind when defining constraints that ORCA counts atoms from no. 0 (not 1).
Applying a distance constraint between atoms 0 and 1. Distance will be fixed at current value.

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{B 0 1 C} # B for Bond, C for Constraint
end
end

Applying two distance constraints (between atoms 0 and 1 and between atoms 0 and 2):

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{B 0 1 C}
{B 0 2 C}
end
end

Applying a distance constraint with a specific value between atoms 0 and 1 (here constraining the bond between atoms 0 and 1 to be fixed at 1.25 Å):

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{B 0 1 1.25 C}
end
end

Applying an angle constraint between atoms 0,1 and 2. Angle will be fixed at current value.

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{A 0 1 2 C } #A for Angle
end
end

Applying a dihedral angle constraint specified by the atoms 0, 1, 2 and 3. Dihedral angle will be fixed at current value.

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{D 0 1 2 3 C } # D for Dihedral angle
end
end

Applying a Cartesian constraint:

! BP86 def2-SVP def2/J Opt
%geom
Constraints
{C 5 C} # Constraining atom no. 5 in space.
end
end

Constrain everything but hydrogen positions:

! BP86 def2-SVP def2/J Opt
%geom
optimizehydrogens true
end

Relaxed and Unrelaxed Surface Scans

Surface scans are used for many purposes: PES exploration, transition state searches, finding new minima etc. In a surface scan (relaxed or unrelaxed), an internal coordinate (bond, angle or dihedral) is held fixed (applying a constraint like above) and then gradually changed (scan step) until it reaches a certain value. In a relaxed scan, a geometry optimisation is performed in each scan step, relaxing the geometry (all other coordinates) but keeping the constraint in place. In an unrelaxed scan, there is no optimisation performed in a scan step and hence no relaxation (much cheaper).

Keep in mind when defining constraints that ORCA counts atoms from no. 0 (not 1).

Relaxed scan

Perform a relaxed surface scan, scanning the bond length between atoms 0 and 1 from 1.0 Å to 3.0 Å in 12 steps.

! BP def2-SVP def2/J Opt

%geom Scan
B 0 1 = 1.0, 3.0, 12
end
end
*xyz 0 1
H 0.0 0.0 0.0 # This is atom zero
F 0.0 0.0 1.0 # This is atom one
*

Scanning an angle (here scanning the H₂O angle ( from 104.4° to 3° in 12 steps):

! BP def2-SVP def2/J Opt

%geom Scan
A 0 1 2 = 104.4, 3.0, 12
end
end
*xyz 0 1
O 0.000000 -0.005910 0.000000
H 0.765973 0.587955 0.000000
H -0.765973 0.587955 0.000000
*

Scanning a dihedral (here scanning the dihedral of ethane from 0° to 90° in 10 steps):

! BP def2-SVP def2/J Opt
%geom Scan
D 3 0 4 6 = 0, 90, 10
end
end
*xyz 0 1
C 0.277784124 0.000000000 -4.653971743
H 0.277784124 1.019690000 -5.053227743
H 1.160862124 -0.509845000 -5.053227743
H -0.605293876 -0.509845000 -5.053227743
C 0.277784124 0.000000000 -3.123983743
H -0.377660283 0.781128227 -2.724727743
H -0.070970657 -0.958195073 -2.724727743
H 1.281983313 0.177066846 -2.724727743
*

Unrelaxed scan

There are a few ways of performing an unrelaxed scan. The simplest is via the %paras block that allows you to define internal coordinates that can be varied. Below the parameter R is defined initially to be 1.05 Å and is then gradually changed (in 8 steps) to become 1.13 Å. As can be seen in the coordinate block below, R corresponds here to the Z coordinate for the second nitrogen atom in Cartesian space. Varying this parameter will thus change the bond length.

! BP def2-SVP def2/J
%paras
R= 1.05,1.13,8
end
* xyz 0 1
N 0 0 0
N 0 0 {R}
*

For more complicated unrelaxed scans, internal coordinates might have to be defined. Internal coordinates can be created in multiple ways (e.g. by ORCA, Molden, Chemcraft etc.). Here, the %paras block in combination with internal coordinates is used for twisting the double bond of ethylene (dihedral is changed from 0° to 180° in 18 steps).

! BP def2-SVP def2/J

%paras
Alpha=0,180,18
end

* int 0 1
C 0 0 0 0 0 0
C 1 0 0 1.3385 0 0
H 1 2 0 1.07 120 0
H 1 2 3 1.07 120 180
H 2 1 3 1.07 120 {Alpha}
H 2 1 3 1.07 120 {Alpha+180}
*