Geometry optimization: what happens in the algorithm?

Question

Geometry optimization using density functional theory (DFT) is done by moving the atoms of a molecule to get the most stable structure with the lowest possible ground state energy. My question is what are the equations involved in this procedure? Where are the electrons in this step?

Answer 1

GDIIS (Geometry Direct Inversion in Iterative Subspace)

This is a very popular method, implemented in almost all the major quantum chemistry software packages, which is quite different from the quasi-Newton and related approaches mentioned in Shoubhik's excellent answer.

DIIS was originally published in 1982 by Peter Pulay for accelerating SCF convergence, and has become the default method for SCF convergence (and also the convergence of many other things now, such as coupled cluster equations) in almost every electronic structure software I know, including (each link takes you to the software's DIIS page in the documentation, just search "DIIS") MOLPRO, Psi4, Q-Chem, CFOUR, MRCC, etc.

In 1984, Csaszar and Pulay published a paper titled "Geometry optimization by direct inversion in the iterative subspace", which has now known as "Geometry DIIS" or "GDIIS". The original paper is surprisingly simple, being only 4 pages (including all references!) of easy-to-follow steps, all of which I will summarize here, in the same notation as in the original paper:

Geometries $\mathbf{x}_i$ generated in the first $m$ cycles (each labeled by the index $i$) are used to generate an improved geometry for the $(m+1)^{\textrm{th}}$ cycle, in the following way:

$$\tag{1} \mathbf{x}_{m+1} = \sum_{i=1}^m c_i \mathbf{x}_i - \mathbf{H}^{-1}\left(\sum_{i=1}^m c_i\mathbf{g}_{i} \right), $$

where $\mathbf{g}_i$ is the gradient vector at iteration $i$ and $\mathbf{H}$ is the Hessian matrix or an approximation to it, and the $c_i$ coefficients are obtained by simple linear algebra:

$$ \begin{pmatrix} \langle \mathbf{e}_1 | \mathbf{e}_1\rangle & \cdots & \langle \mathbf{e}_1 | \mathbf{e}_m\rangle & 1\\ \vdots & \ddots & \vdots & \vdots\\ \langle \mathbf{e}_m | \mathbf{e}_1\rangle & \cdots & \langle \mathbf{e}_m | \mathbf{e}_m\rangle & 1\\ 1 & \cdots & 1 & 0\\ \end{pmatrix} \begin{pmatrix} c_1 \\ \vdots \\ c_m \\ -\lambda \end{pmatrix} = \begin{pmatrix} 0 \\ \vdots \\ 0 \\ -\lambda \end{pmatrix},\tag{2} $$

with $\lambda$ being the Lagrange multiplier and $\mathbf{e}_i$ being the error between $\mathbf{x}_i$ and the "true" final geometry $\mathbf{x}_f$, which in the case of a quadratic energy function can be obtained by:

$$ \mathbf{e}_i = -\mathbf{H}^{-1}\mathbf{g}_i.\tag{3} $$

That's it!

If you're interested in more detail I will also say:

• The reason why we are solving Eq. 2 is because we want $\sum c_i = 1$ and $\sum c_i \mathbf{e}_i$.
• You can get faster convergence if $\mathbf{H}$ is updated at each iteration, so we could replace it with $\mathbf{H}_i$ in all of the above.
• You don't need to include all previous $\mathbf{x}_i$. To save RAM and cost you can choose to use only the last few geometries, or perhaps only the geometries that are within some threshold of the most recent (ideally most accurate) geometry.
• Although the energy landscape is not exactly quadratic, every differential function with a minimum has a quadratic approximation near the minimum (e.g. the quadratic Taylor polynomial), so the approximation in Eq. 3 is quite accurate if the geoemtry at the $i^{\textrm{th}}$ iteration is close enough to the minimum (it doesn't even have to be very "close" to the minimum, just enough that energy landscape looks quadratic there).

Answer 2

To add to the other answer, the energy of a molecule (or any chemical species) can be expressed as a function of the coordinates of the atoms. Most softwares convert the 3N cartesian coordinates into internal coordinates (3N-6 for non-linear molecules, 3N-5 for linear), removing the rotational and translational degrees of freedom. Usually, the optimization is done in internal coordinates by most QM softwares. Internal coordinates mean bonds, angles, dihedrals etc.

If we have a set of coordinates $x_1, x_2, x_3, .... ,x_n$, then we can carry out the SCF procedure for that set of coordinates (i.e. Born-Oppenheimer approx.). The energy is then a function of the coordinates: $$E=f(x_1,x_2,x_3,...,x_n)$$ This is an optimization problem, i.e. we seek to minimize E by varying the variables $x_1, x_2,..., x_n$. There are a lot of algorithms that do this. Newton-Raphson (NR) has been mentioned. It is an efficient algorithm, however it needs accurate second derivatives (hessian). The second derivatives are time-consuming to calculate, so most QM codes use what are known are quasi-Newton methods.

These methods require only accurate gradient (first derivatives), and an approximate (guess) hessian. The gradient can be calculated directly from the SCF solution, so these quasi-Newton methods are much faster. Examples of quasi-Newton algorithms include BFGS, L-BFGS, Rational Function Optimization(RFO) etc.

The optimization procedure is iterative. In quasi-Newton methods, the hessian matrix is also updated after each step in optimization. Note that instead of guessing the hessian, we can also calculate it and that will lead to quicker convergence.

So in brief, the general steps of optimization are:

1. Start with a set of coordinates
2. Carry out SCF and get the energy, gradient (optionally hessian)
3. Feed the gradient and hessian into the optimizer algorithm, which will then tell how to vary the coordinates to go towards the minimum (this is called a step)
4. Check if the step is too large, if yes then scale it down (optional)
5. Take the step (i.e. change the coordinates)
6. Update the hessian
7. Go to step 1

The optimization procedure iterates through these steps until a certain criteria for convergence is reached. Usually this criteria would be that the gradient is lower than a certain threshold (close to zero), or the next step size is lower than a certain value etc. This type of optimization is usually done in QM softwares like Gaussian, GAMESS, Orca etc.

Additionally, there are modified versions of NR such as Eigenvector Following (EF) algorithm, which can use accurate hessian to actually move towards maxima (transition states).

---

In molecular dynamics, usually the system is so large that even storing the guess hessian is problematic. So, first-order optimizer algorithms are used e.g.- Conjugate gradient, Steepest descent etc. These optimizers do not require the hessian at all, and only work with the gradient. Therefore, they are slower to converge than second order methods. However, that's okay for MD as the energy/gradient calculation steps are much much faster than SCF calculation.

---

Edit: In your question, you have asked about the electrons. To make things clear, there are two separate things here—one is the optimization procedure and the other is the SCF. The optimization procedure needs only the coordinates, the gradient and the hessian and nothing else. It has nothing to do with electrons. The SCF part is where the energy, gradient (and hessian) is calculated keeping the coordinates fixed. The SCF equations (e.g. Hartree-Fock) will involve the electrons.

Answer 3

First, realize that optimization is a general thing that can be done for all types of problems.

In the geometry optimization of atoms and molecules, what we want to find is the configuration of the nuclei which is a minima on the potential energy surface. In other words, we want to find the positions, x which satisfy:

$\min(E[x])$

Secondly, the quantum mechanical energy of a molecule is dependent on the positions of the electrons and parameterically dependent on the position of the nuclei E[r,x] where r are the positions of the electrons. This is known as the Born-Oppenheimer approximation. So when we are dealing with geometry optimization we don't care where the electrons are because we can assume they are in their proper places for that geometry. In fact we can use molecular mechanics to give us the energy and there are no electrons in that model, only internuclear potentials.

There are many algorithms to find the minimum energy geometry. One of the most popular is a second order method called Newton-Raphson which constructs a local quadratic approximation to the potential energy surface:

$E(x) = E(x_{0}) + g_0^T\Delta x + 1/2\Delta x^TH_0\Delta x $

where $g_0$ is the gradient (dE/dx) at $x_0$, $H_0$ is the Hessian (d^2E/dx^2) at $x_0$ and $\Delta x = x-x_0$

So technically this is the equation which is "solved" when you are doing a geometry optimization. But there are many other ways of solving it as well, for example see my other answer here: https://mattermodeling.stackexchange.com/a/1069/52

For a good overview of geometry optimization I recommend Geometry optimization by H. B. Schlegel in Wiley Interdisciplinary Reviews: Computational Molecular Science (WIREs Comput Mol Sci, 2011, 790–809; doi 10.1002/wcms.34).