Energy minimization

In the field of computational chemistry, energy minimization (also called energy optimization, geometry minimization, or geometry optimization) is the process of finding an arrangement in space of a collection of atoms where, according to some computational model of chemical bonding, the net inter-atomic force on each atom is acceptably close to zero and the position on the potential energy surface (PES) is a stationary point (described later).

Typically, but not always, the process seeks to find the geometry of a particular arrangement of the atoms that represents a local or global energy minimum.

[1] Additionally, certain coordinates (such as a chemical bond length) might be fixed during the optimization.

This could be the set of the Cartesian coordinates of the atoms or, when considering molecules, might be so called internal coordinates formed from a set of bond lengths, bond angles and dihedral angles.

Geometry optimization is then a mathematical optimization problem, in which it is desired to find the value of r for which E(r) is at a local minimum, that is, the derivative of the energy with respect to the position of the atoms, ∂E/∂r, is the zero vector and the second derivative matrix of the system,

The computational model that provides an approximate E(r) could be based on quantum mechanics (using either density functional theory or semi-empirical methods), force fields, or a combination of those in case of QM/MM.

Using this computational model and an initial guess (or ansatz) of the correct geometry, an iterative optimization procedure is followed, for example: As described above, some method such as quantum mechanics can be used to calculate the energy, E(r) , the gradient of the PES, that is, the derivative of the energy with respect to the position of the atoms, ∂E/∂r and the second derivative matrix of the system, ∂∂E/∂ri∂rj, also known as the Hessian matrix, which describes the curvature of the PES at r. An optimization algorithm can use some or all of E(r) , ∂E/∂r and ∂∂E/∂ri∂rj to try to minimize the forces and this could in theory be any method such as gradient descent, conjugate gradient or Newton's method, but in practice, algorithms which use knowledge of the PES curvature, that is the Hessian matrix, are found to be superior.

For most systems of practical interest, however, it may be prohibitively expensive to compute the second derivative matrix, and it is estimated from successive values of the gradient, as is typical in a Quasi-Newton optimization.

Additionally, Cartesian coordinates are highly correlated, that is, the Hessian matrix has many non-diagonal terms that are not close to zero.

This can lead to numerical problems in the optimization, because, for example, it is difficult to obtain a good approximation to the Hessian matrix and calculating it precisely is too computationally expensive.

[4] Some degrees of freedom can be eliminated from an optimization, for example, positions of atoms or bond lengths and angles can be given fixed values.

Figure 1 depicts a geometry optimization of the atoms in a carbon nanotube in the presence of an external electrostatic field.

Transition state structures can be determined by searching for saddle points on the PES of the chemical species of interest.

Defined mathematically, an nth order saddle point is characterized by the following: ∂E/∂r = 0 and the Hessian matrix, ∂∂E/∂ri∂rj, has exactly n negative eigenvalues.

Local methods are suitable when the starting point for the optimization is very close to the true transition state (very close will be defined shortly) and semi-global methods find application when it is sought to locate the transition state with very little a priori knowledge of its geometry.

Further, the eigenvector with the most negative eigenvalue must correspond to the reaction coordinate, that is, it must represent the geometric transformation relating to the process whose transition state is sought.

Given the above pre-requisites, a local optimization algorithm can then move "uphill" along the eigenvector with the most negative eigenvalue and "downhill" along all other degrees of freedom, using something similar to a quasi-Newton method.

The Activation Relaxation Technique (ART)[7][8][9] is also an open-ended method to find new transition states or to refine known saddle points on the PES.

The method follows the direction of lowest negative curvature (computed using the Lanczos algorithm) on the PES to reach the saddle point, relaxing in the perpendicular hyperplane between each "jump" (activation) in this direction.

The generated approximate geometry can then serve as a starting point for refinement via a local search, which was described above.

For this to be achieved, spacing constraints must be applied so that each bead ri does not simply get optimized to the reactant and product geometry.

where I is the identity matrix and τi is a unit vector representing the reaction path tangent at ri.

It operates by taking interpolated points between the reactant and product geometries and choosing the one with the highest energy for subsequent refinement via a local search.

The quadratic synchronous transit (QST) method extends LST by allowing a parabolic reaction path, with optimization of the highest energy point orthogonally to the parabola.

In Nudged elastic band (NEB)[11] method, the beads along the reaction pathway have simulated spring forces in addition to the chemical forces, -∂E/∂ri, to cause the optimizer to maintain the spacing constraint.

In a traditional implementation, the point with the highest energy is used for subsequent refinement in a local search.

There are many variations on the NEB method,[12] such including the climbing image NEB, in which the point with the highest energy is pushed upwards during the optimization procedure so as to (hopefully) give a geometry which is even closer to that of the transition state.

There have also been extensions[13] to include Gaussian process regression for reducing the number of evaluations.

Geometry optimization, by contrast, does not produce a "trajectory" with any physical meaning – it is concerned with minimization of the forces acting on each atom in a collection of atoms, and the pathway via which it achieves this lacks meaning.

Different optimization algorithms could give the same result for the minimum energy structure, but arrive at it via a different pathway.