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). The collection of atoms might be a single molecule, an ion, a condensed phase, a transition state or even a collection of any of these. The computational model of chemical bonding might, for example, be quantum mechanics.

As an example, when optimizing the geometry of a water molecule, one aims to obtain the hydrogen-oxygen bond lengths and the hydrogen-oxygen-hydrogen bond angle which minimize the forces that would otherwise be pulling atoms together or pushing them apart.

The motivation for performing a geometry optimization is the physical significance of the obtained structure: optimized structures often correspond to a substance as it is found in nature and the geometry of such a structure can be used in a variety of experimental and theoretical investigations in the fields of chemical structure, thermodynamics, chemical kinetics, spectroscopy and others.

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. Instead of searching for global energy minimum, it might be desirable to optimize to a transition state, that is, a saddle point on the potential energy surface.^[1] Additionally, certain coordinates (such as a chemical bond length) might be fixed during the optimization.

Molecular geometry and mathematical interpretation

The geometry of a set of atoms can be described by a vector of the atoms' positions. 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.

Given a set of atoms and a vector, $r$ , describing the atoms' positions, one can introduce the concept of the energy as a function of the positions, $E (r)$ . 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, $\partial E /\partial r$ , is the zero vector and the second derivative matrix of the system, ${\begin{pmatrix}{\frac {\partial ^{2}E}{\partial r_{i}\,\partial r_{j}}}\end{pmatrix}}_{ij}$ , also known as the Hessian matrix, which describes the curvature of the PES at $r$ , has all positive eigenvalues (is positive definite).

A special case of a geometry optimization is a search for the geometry of a transition state; this is discussed below.

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:

calculate the force on each atom (that is, $-\partial E /\partial r$ )
if the force is less than some threshold, finish
otherwise, move the atoms by some computed step $∆ r$ that is predicted to reduce the force
repeat from the start

Practical aspects of optimization

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, $\partial E /\partial r$ and the second derivative matrix of the system, $\partial\partial E /\partial r i \partial r j$ , 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)$ , $\partial E /\partial r$ and $\partial\partial E /\partial r i \partial r j$ 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.

The choice of the coordinate system can be crucial for performing a successful optimization. Cartesian coordinates, for example, are redundant since a non-linear molecule with $N$ atoms has $3 N -6$ vibrational degrees of freedom whereas the set of Cartesian coordinates has $3 N$ dimensions. 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. However, in case that energy is expressed with standard force fields, computationally efficient methods have been developed ^[2] able to derive analytically the Hessian matrix in Cartesian coordinates while preserving a computational complexity of the same order to that of gradient computations. Internal coordinates tend to be less correlated but are more difficult to set-up and it can be difficult to describe some systems, such as ones with symmetry or large condensed phases.^[3] Many modern computational chemistry software packages contain automatic procedures for the automatic generation of reasonable coordinate systems for optimization.^[4]

Degree of freedom restriction

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. Sometimes these are referred to as being frozen degrees of freedom.

Figure 1 depicts a geometry optimization of the atoms in a carbon nanotube in the presence of an external electrostatic field. In this optimization, the atoms on the left have their positions frozen. Their interaction with the other atoms in the system are still calculated, but alteration the atoms' position during the optimization is prevented.

Figure 1: Electrostatic deflections of a carbon nanotube in an electric field.

Transition state optimization

Transition state structures can be determined by searching for saddle points on the PES of the chemical species of interest.^[5] A first-order saddle point is a position on the PES corresponding to a minimum in all directions except one; a second-order saddle point is a minimum in all directions except two, and so on. Defined mathematically, an nth order saddle point is characterized by the following: $\partial E /\partial r = 0$ and the Hessian matrix, $\partial\partial E /\partial r i \partial r j$ , has exactly n negative eigenvalues.

Algorithms to locate transition state geometries fall into two main categories: local methods and semi-global methods. 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. Some methods, such as the Dimer method (see below), fall into both categories.

Local searches

A so-called local optimization requires an initial guess of the transition state that is very close to the true transition state. Very close typically means that the initial guess must have a corresponding Hessian matrix with one negative eigenvalue, or, the negative eigenvalue corresponding to the reaction coordinate must be greater in magnitude than the other negative eigenvalues. 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.

Dimer method

The dimer method^[6] can be used to find possible transition states without knowledge of the final structure or to refine a good guess of a transition structure. The “dimer” is formed by two images very close to each other on the PES. The method works by moving the dimer uphill from the starting position whilst rotating the dimer to find the direction of lowest curvature (ultimately negative).

Activation Relaxation Technique (ART)

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.

Chain-of-state methods

Chain-of-state^[10] methods can be used to find the approximate geometry of the transition state based on the geometries of the reactant and product. The generated approximate geometry can then serve as a starting point for refinement via a local search, which was described above.

Chain-of-state methods use a series of vectors, that is points on the PES, connecting the reactant and product of the reaction of interest, $r reactant$ and $r product$ , thus discretizing the reaction pathway. Very commonly, these points are referred to as beads due to an analogy of a set of beads connected by strings or springs, which connect the reactant and products. The series of beads is often initially created by interpolating between $r reactant$ and $r product$ , for example, for a series of $N + 1$ beads, bead $i$ might be given by

$\mathbf {r} _{i}={\frac {i}{N}}\mathbf {r} _{\mathrm {product} }+\left(1-{\frac {i}{N}}\right)\mathbf {r} _{\mathrm {reactant} }$

where $i \in 0, 1, ..., N$ . Each of the beads $r i$ has an energy, $E (r i)$ , and forces, $-\partial E /\partial r i$ and these are treated with a constrained optimization process that seeks to get an as accurate as possible representation of the reaction pathway. For this to be achieved, spacing constraints must be applied so that each bead $r i$ does not simply get optimized to the reactant and product geometry.

Often this constraint is achieved by projecting out components of the force on each bead $r i$ , or alternatively the movement of each bead during optimization, that are tangential to the reaction path. For example, if for convenience, it is defined that $g i = \partial E /\partial r i$ , then the energy gradient at each bead minus the component of the energy gradient that is tangential to the reaction pathway is given by

$\mathbf {g} _{i}^{\perp }=\mathbf {g} _{i}-\mathbf {\tau } _{i}(\mathbf {\tau } _{i}\cdot \mathbf {g} _{i})=\left(I-\mathbf {\tau } _{i}\mathbf {\tau } _{i}^{T}\right)\mathbf {g} _{i}$

where $I$ is the identity matrix and $τ i$ is a unit vector representing the reaction path tangent at $r i$ . By projecting out components of the energy gradient or the optimization step that are parallel to the reaction path, an optimization algorithm significantly reduces the tendency of each of the beads to be optimized directly to a minimum.

Synchronous transit

The simplest chain-of-state method is the linear synchronous transit (LST) method. 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.

Nudged elastic band

In Nudged elastic band (NEB)^[11] method, the beads along the reaction pathway have simulated spring forces in addition to the chemical forces, $-\partial E /\partial r i$ , to cause the optimizer to maintain the spacing constraint. Specifically, the force $f i$ on each point i is given by

$\mathbf {f} _{i}=\mathbf {f} _{i}^{\parallel }-\mathbf {g} _{i}^{\perp }$

where

$\mathbf {f} _{i}^{\parallel }=k\left[\left(\left(\mathbf {r} _{i+1}-\mathbf {r} _{i}\right)-\left(\mathbf {r} _{i}-\mathbf {r} _{i-1}\right)\right)\cdot \tau _{i}\right]\tau _{i}$

is the spring force parallel to the pathway at each point $r i$ (k is a spring constant and $τ i$ , as before, is a unit vector representing the reaction path tangent at $r i$ ).

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. For systems with non-Euclidean (R^2) geometry, like magnetic systems, the method is modified to the geodesic nudged elastic band approach.^[14]^[15]

String method

The string method^[16]^[17]^[18] uses splines connecting the points, $r i$ , to measure and enforce distance constraints between the points and to calculate the tangent at each point. In each step of an optimization procedure, the points might be moved according to the force acting on them perpendicular to the path, and then, if the equidistance constraint between the points is no-longer satisfied, the points can be redistributed, using the spline representation of the path to generate new vectors with the required spacing.

Variations on the string method include the growing string method,^[19] in which the guess of the pathway is grown in from the end points (that is the reactant and products) as the optimization progresses.

Comparison with other techniques

Geometry optimization is fundamentally different from a molecular dynamics simulation. The latter simulates the motion of molecules with respect to time, subject to temperature, chemical forces, initial velocities, Brownian motion of a solvent, and so on, via the application of Newton's laws of Motion. This means that the trajectories of the atoms which get computed, have some physical meaning. 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.

References

^ "Input reference of CP2K version trunk, Section GEO_OPT, Keyword TYPE". CP2K. Retrieved 30 April 2015.
^ Chatzieleftheriou, S.; Adendorff, M. R.; Lagaros, N. D. (2016). "Generalized Potential Energy Finite Elements for Modeling Molecular Nanostructures". J. Chem. Inf. Model. 56 (10): 1963–1978. doi:10.1021/acs.jcim.6b00356. PMID 27653992.
^ Peng, C.; Ayala, P. Y.; Schlegel, H. B. (1996). "Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States". Journal of Computational Chemistry. 17 (1): 49–56. doi:10.1002/(sici)1096-987x(19960115)17:1<49::aid-jcc5>3.3.co;2-#.
^ "Home". gaussian.com.
^ Frank Jensen (1999). Introduction to Computational Chemistry. England: John Wiley and Sons Ltd.
^ Graeme Henkelman; Hannes Jónsson (1999). "A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives". J. Chem. Phys. 111 (15): 7010–7022. Bibcode:1999JChPh.111.7010H. doi:10.1063/1.480097.
^ G.T. Barkema; Normand Mousseau (1996). "Event-Based Relaxation of Continuous Disordered Systems". Phys. Rev. Lett. 77 (21): 4358–4361. arXiv:cond-mat/9607156. Bibcode:1996PhRvL..77.4358B. doi:10.1103/PhysRevLett.77.4358. PMID 10062518. S2CID 27932059.
^ Rachid Malek; Normand Mousseau (2011). "Optimized energy landscape exploration using the ab initio based activation-relaxation technique". Physical Review E. 135 (6): 7723–7728. arXiv:cond-mat/0006042. Bibcode:2000PhRvE..62.7723M. doi:10.1103/PhysRevE.62.7723. PMID 11138044. S2CID 119453527.
^ Eduardo Machado-Charry; Laurent Karim Béland; Damien Caliste; Luigi Genovese; Thierry Deutsch; Normand Mousseau; Pascal Pochet (2011). "Optimized energy landscape exploration using the ab initio based activation-relaxation technique". J. Chem. Phys. 62 (3): 034102–034112. Bibcode:2011JChPh.135c4102M. doi:10.1063/1.3609924. PMID 21786982.
^ Jensen, F. Introduction to Computational Chemistry; Wiley: 2 ed.; 2006
^ (a) G. Mills and H. Jónsson, Phys. Rev. Lett. 72, 1124 (1994) (b) Graeme Henkelman and Hannes Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113, 9978 - 9985 (2000)
^ "Nudged Elastic Band". UT Austin. Archived from the original on 2014-02-03.
^ Koistinen, Olli-Pekka; Dagbjartsdóttir, Freyja B.; Ásgeirsson, Vilhjálmur; Vehtari, Aki; Jónsson, Hannes (2017-10-21). "Nudged elastic band calculations accelerated with Gaussian process regression". The Journal of Chemical Physics. 147 (15): 152720. arXiv:1706.04606. Bibcode:2017JChPh.147o2720K. doi:10.1063/1.4986787. ISSN 0021-9606. PMID 29055305. S2CID 21822734.
^ P.F. Bessarab, V.M. Uzdin, H. Jónsson (2015). "Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation". Comput. Phys. Commun. 196: 335–347. arXiv:1502.05065. doi:10.1016/j.cpc.2015.07.001.{{cite journal}}: CS1 maint: multiple names: authors list (link)
^ Ivanov, A V; Dagbartsson, D; Tranchida, J; Uzdin, V M; Jónsson, H (2020-08-12). "Efficient optimization method for finding minimum energy paths of magnetic transitions". Journal of Physics: Condensed Matter. 32 (34): 345901. arXiv:2001.10372. Bibcode:2020JPCM...32H5901I. doi:10.1088/1361-648X/ab8b9c. ISSN 0953-8984. PMID 32316000. S2CID 210932577.
^ "Rare Events, Transition Pathways and Reaction Rates". and "The string method page".
^ Weinan E, Weiqing Ren, Eric Vanden-Eijnden (2002). "String method for the study of rare Events". Phys. Rev. B. 66 (5): 052301. arXiv:cond-mat/0205527. Bibcode:2002PhRvB..66e2301E. doi:10.1103/PhysRevB.66.052301. S2CID 119326534.{{cite journal}}: CS1 maint: multiple names: authors list (link)
^ Amit Samanta; Weinan E (2010). "Modified string method for finding minimum energy path". arXiv:1009.5612 [physics.comp-ph].
^ Baron Peters; Andreas Heyden; Alexis T. Bell; Arup Chakraborty (2004). "A growing string method for determining transition states: Comparison to the nudged elastic band and string methods". J. Chem. Phys. 120 (17): 7877–7886. Bibcode:2004JChPh.120.7877P. doi:10.1063/1.1691018. PMID 15267702.

External links

Numerical Recipes in Fortran 77

Additional references

Payne et al., "Iterative minimization techniques for ab initio total-energy calculations: Molecular dynamics and conjugate gradients", Reviews of Modern Physics 64 (4), pp. 1045–1097. (1992) (abstract)
Stich et al., "Conjugate gradient minimization of the energy functional: A new method for electronic structure calculation", Physical Review B 39 (8), pp. 4997–5004, (1989)
Chadi, "Energy-minimization approach to the atomic geometry of semiconductor surfaces", Physical Review Letters 41 (15), pp. 1062–1065 (1978)

[1] "Input reference of CP2K version trunk, Section GEO_OPT, Keyword TYPE". CP2K. Retrieved 30 April 2015.

[2] Chatzieleftheriou, S.; Adendorff, M. R.; Lagaros, N. D. (2016). "Generalized Potential Energy Finite Elements for Modeling Molecular Nanostructures". J. Chem. Inf. Model. 56 (10): 1963–1978. doi:10.1021/acs.jcim.6b00356. PMID 27653992.

[3] Peng, C.; Ayala, P. Y.; Schlegel, H. B. (1996). "Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States". Journal of Computational Chemistry. 17 (1): 49–56. doi:10.1002/(sici)1096-987x(19960115)17:1<49::aid-jcc5>3.3.co;2-#.

[4] "Home". gaussian.com.

[5] Frank Jensen (1999). Introduction to Computational Chemistry. England: John Wiley and Sons Ltd.

[6] Graeme Henkelman; Hannes Jónsson (1999). "A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives". J. Chem. Phys. 111 (15): 7010–7022. Bibcode:1999JChPh.111.7010H. doi:10.1063/1.480097.

[7] G.T. Barkema; Normand Mousseau (1996). "Event-Based Relaxation of Continuous Disordered Systems". Phys. Rev. Lett. 77 (21): 4358–4361. arXiv:cond-mat/9607156. Bibcode:1996PhRvL..77.4358B. doi:10.1103/PhysRevLett.77.4358. PMID 10062518. S2CID 27932059.

[8] Rachid Malek; Normand Mousseau (2011). "Optimized energy landscape exploration using the ab initio based activation-relaxation technique". Physical Review E. 135 (6): 7723–7728. arXiv:cond-mat/0006042. Bibcode:2000PhRvE..62.7723M. doi:10.1103/PhysRevE.62.7723. PMID 11138044. S2CID 119453527.

[9] Eduardo Machado-Charry; Laurent Karim Béland; Damien Caliste; Luigi Genovese; Thierry Deutsch; Normand Mousseau; Pascal Pochet (2011). "Optimized energy landscape exploration using the ab initio based activation-relaxation technique". J. Chem. Phys. 62 (3): 034102–034112. Bibcode:2011JChPh.135c4102M. doi:10.1063/1.3609924. PMID 21786982.

[10] Jensen, F. Introduction to Computational Chemistry; Wiley: 2 ed.; 2006

[11] (a) G. Mills and H. Jónsson, Phys. Rev. Lett. 72, 1124 (1994) (b) Graeme Henkelman and Hannes Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113, 9978 - 9985 (2000)

[12] "Nudged Elastic Band". UT Austin. Archived from the original on 2014-02-03.

[13] Koistinen, Olli-Pekka; Dagbjartsdóttir, Freyja B.; Ásgeirsson, Vilhjálmur; Vehtari, Aki; Jónsson, Hannes (2017-10-21). "Nudged elastic band calculations accelerated with Gaussian process regression". The Journal of Chemical Physics. 147 (15): 152720. arXiv:1706.04606. Bibcode:2017JChPh.147o2720K. doi:10.1063/1.4986787. ISSN 0021-9606. PMID 29055305. S2CID 21822734.

[14] P.F. Bessarab, V.M. Uzdin, H. Jónsson (2015). "Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation". Comput. Phys. Commun. 196: 335–347. arXiv:1502.05065. doi:10.1016/j.cpc.2015.07.001.{{cite journal}}: CS1 maint: multiple names: authors list (link)

[15] Ivanov, A V; Dagbartsson, D; Tranchida, J; Uzdin, V M; Jónsson, H (2020-08-12). "Efficient optimization method for finding minimum energy paths of magnetic transitions". Journal of Physics: Condensed Matter. 32 (34): 345901. arXiv:2001.10372. Bibcode:2020JPCM...32H5901I. doi:10.1088/1361-648X/ab8b9c. ISSN 0953-8984. PMID 32316000. S2CID 210932577.

[16] "Rare Events, Transition Pathways and Reaction Rates". and "The string method page".

[17] Weinan E, Weiqing Ren, Eric Vanden-Eijnden (2002). "String method for the study of rare Events". Phys. Rev. B. 66 (5): 052301. arXiv:cond-mat/0205527. Bibcode:2002PhRvB..66e2301E. doi:10.1103/PhysRevB.66.052301. S2CID 119326534.{{cite journal}}: CS1 maint: multiple names: authors list (link)

[18] Amit Samanta; Weinan E (2010). "Modified string method for finding minimum energy path". arXiv:1009.5612 [physics.comp-ph].

[19] Baron Peters; Andreas Heyden; Alexis T. Bell; Arup Chakraborty (2004). "A growing string method for determining transition states: Comparison to the nudged elastic band and string methods". J. Chem. Phys. 120 (17): 7877–7886. Bibcode:2004JChPh.120.7877P. doi:10.1063/1.1691018. PMID 15267702.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]