What are the methods that can be used to ensure that an optimized geometry is a local minimum?

Question

Usually it is suggested to do a phonon calculation. What are the other ways?

Answer 1

This depends on what you are studying. For molecular systems without periodicity, the simplest approach is to carry out a vibrational frequency analysis and confirm that there are no imaginary modes. It is considered standard to carry out a vibrational frequency analysis for all investigated structures, provided the number of investigated systems is not prohibitively large.

For extended solids, which it sounds like you're modeling, it is much less common to confirm the result of a structure relaxation is indeed a minimum, although it can't hurt and is always a good idea if feasible. A phonon calculation can indeed be used to confirm your location in the potential energy surface.

There are other options though for extended solids. If you are modeling the adsorption of a molecule onto a surface, you can likely assume that the surface atoms are near or at their minimum energy positions (if you have pre-relaxed the adsorbate-free structure). In this case, you could perform a vibrational frequency analysis of just the adsorbate atoms (or, even better, the adsorbate atoms and a few nearby surface atoms) to confirm there are no imaginary modes. ASE has an excellent vibrational analysis function that works with many of the most popular DFT packages, which I recommend if you're using a code like VASP. Beyond this, you can always try slightly modifying the atomic positions of the converged structure in several different ways, re-optimizing each structure to see if the same low-energy configuration is found. Of course, this does not guarantee you have arrived at a local minimum, but it can provide further support.

Answer 2

In order to check if a geometry is a local minimum, it is a necessary and sufficient condition that the Hessian is positive (semi)definite, i.e. that the lowest eigenvalue of the nuclear Hessian is non-negative.

Namely, expanding the energy $E({\bf R})$ around the reference point ${\bf R}_0$ you have the Taylor expansion $E({\bf R}) = E({\bf R}_0) - {\bf g} \cdot ({\bf R}-{\bf R}_0) + \frac 1 2 ({\bf R}-{\bf R}_0)^{\rm T} \cdot {\bf H} \cdot ({\bf R}-{\bf R}_0) + \mathcal{O}(|({\bf R}-{\bf R}_0)|^3)$ where ${\bf g} = -[\nabla E({\bf R})]_{{\bf R}={\bf R}_0}$ is the gradient and $H_{ij} = -\partial^2 E / \partial R_i \partial R_j $ is the nuclear Hessian, which has the size $3N_{\rm atoms} \times 3N_{\rm atoms}$.

At an extremum ${\bf g}={\bf 0}$; this is what the optimizer finds for you. However, in some cases you can find that the Hessian has negative eigenvalues $\lambda_i$ with corresponding eigenvectors $\boldsymbol{\rho}_i$. Now, if you set ${\bf R}={\bf R}_0+\epsilon \boldsymbol{\rho}_i$, you will find that $E(\epsilon)=E({\bf R})=E({\bf R}_0) + \frac 1 2 \epsilon^2 \lambda_i + \mathcal{O}(|\epsilon|^3)$. Since $\lambda_i<0$, this means you can find a lower energy by moving the atoms slightly along $\boldsymbol{\rho}_i$.

A frequency calculation means computing all eigenvalues of the Hessian, which becomes costly in large systems. (I'm not sure what is done in phonon calculations.) However, since you only need the lowest eigenvalue of the Hessian to check if you are at a local minimum, iterative diagonalizers such as the Davidson method can be used to focus only on the lowest roots of the nuclear Hessian.

Alternatively, there are also some methods that need only gradient information to estimate the lowest eigenvalue of the nuclear Hessian; a "dimer method" has been suggested e.g. by Henkelman and Jónsson in J. Chem. Phys. 111, 7010 (1999).