Geometry Optimization: Ver. 1.0

1 Newton method

The total energy $E$ of a system can be expanded by the Taylor series with respect to atomic coordinates $\{x_i\}$ around $E_0$ with $\{x_i^{(0)}\}$ as follows:

$E$ $=$ $E_0+{\displaystyle \sum _i^{3N}}{\left({\displaystyle \frac{\partial E}{\partial x_i}}\right)}_0(x_i-x_i^{(0)})+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j}^{3N}}{\left({\displaystyle \frac{{\partial }^{2}E}{\partial x_i\partial x_j}}\right)}_0(x_i-x_i^{(0)})(x_j-x_j^{(0)})+\mathrm{\cdots },$ (1)

where the derivatives ${()}_0$ mean the partial derivatives at $\{x_i^{(0)}\}$, and $N$ is the number of atoms. By differentiating Eq. (1) with respect to $x_k$, to the second order we have

${\displaystyle \frac{\partial E}{\partial x_k}}$ $=$ ${\left({\displaystyle \frac{\partial E}{\partial x_k}}\right)}_0+{\displaystyle \sum _i^{3N}}{\left({\displaystyle \frac{{\partial }^{2}E}{\partial x_k\partial x_i}}\right)}_0(x_i-x_i^{(0)}).$ (2)

In case the coordinates $\{x_i\}$ give a local minimum, assuming $\frac{\partial E}{\partial x_k}=0$, we have the following matrix equation:

(3)

The short notation is

$H\mathrm{\Delta }𝐱=-𝐠,$ (4)

where the matrix consisting of the second derivatives in the left-hand side is called Hessian $H$. Using Eq. (4), $\{x_i\}$ can be updated by

${𝐱}^{(n+1)}={𝐱}^{(n)}-{(H^{(n)})}^{-1}{𝐠}^{(n)}.$ (5)

This is the well known Newton method.

2 RMM-DIIS

In the OpenMX, ${𝐱}^{(n)}$ and ${𝐠}^{(n)}$ in Eq. (5) are replaced by ${\overline{𝐱}}^{(n)}$ and ${\overline{𝐠}}^{(𝐧)}$ given by the residual minimization method in the direct inversion of iterative subspace (RMM-DIIS) [1, 2] as follows:

${𝐱}^{(n+1)}={\overline{𝐱}}^{(n)}-\alpha {(H^{(n)})}^{-1}{\overline{𝐠}}^{(n)},$ (6)

where $\alpha$ is a tuning parameter for acceleration of the convergence, which can be small (large) for a large (small) ${\overline{𝐠}}^{(n)}$. $\overline{𝐠}$ in the RMM-DIIS can be found by a linear combination of previous upto p-th gradients $𝐠$ as

${\overline{𝐠}}^{(n)}={\displaystyle \sum _{m=n-(p-1)}^n}{a}_m{𝐠}^{(m)},$ (7)

where $a_m$ is found by minimizing $⟨{\overline{𝐠}}^{(n)}|{\overline{𝐠}}^{(n)}⟩$ with a constraint ${\sum }_{m=n-(p-1)}^n{a}_m=1$. According to Lagrange’s multiplier method, $F$ is defined by

$F$ $=$ $⟨{\overline{𝐠}}^{(n)}|{\overline{𝐠}}^{(n)}⟩-\lambda (1-{\displaystyle \sum _m^n}{a}_m),$ (8)

$=$ ${\displaystyle \sum _{m,m^{\prime }}}{a}_m{a}_{m^{\prime }}⟨{𝐠}^{(m)}|{𝐠}^{(m^{\prime })}⟩-\lambda (1-{\displaystyle \sum _m^n}{a}_m).$

Considering $\frac{\partial F}{\partial a_k}=0$ and $\frac{\partial F}{\partial \lambda }=0$, an optimum set of $\{a\}$ can be found by solving the following linear equation:

(9)

An optimum choice of ${\overline{𝐱}}^{(n)}$ may be obtained by the set of coefficients $\{a\}$ as

${\overline{𝐱}}^{(n)}={\displaystyle \sum _{m=n-(p-1)}^n}{a}_m{𝐱}^{(m)}.$ (10)

If the Hessian $H$ is approximated by the unity $I$, Eq. (6) becomes

${𝐱}^{(n+1)}={\overline{𝐱}}^{(n)}-\alpha {\overline{𝐠}}^{(n)}.$ (11)

This scheme in the Cartesian coordinate has been implemented as ’DIIS’ in OpenMX.

3 Broyden-Fletcher-Goldfarb-Shanno (BFGS) method

Define

$\mathrm{\Delta }{𝐠}^{(n)}$ $=$ ${𝐠}^{(n)}-{𝐠}^{(n-1)},$ (12)

$\mathrm{\Delta }{𝐱}^{(n)}$ $=$ ${𝐱}^{(n)}-{𝐱}^{(n-1)}.$ (13)

Then, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [3] gives the following rank-2 update formula for ${(H^{(n)})}^{-1}$:

${(H^{(n)})}^{-1}$ $=$ ${(H^{(n-1)})}^{-1}+{\displaystyle \frac{⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩+⟨\mathrm{\Delta }{𝐠}^{(n)}|{(H^{(n-1)})}^{-1}|\mathrm{\Delta }{𝐠}^{(n)}⟩}{{\left(⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩\right)}^2}}|\mathrm{\Delta }{𝐱}^{(n)}⟩⟨\mathrm{\Delta }{𝐱}^{(n)}|$ (14)

$-{\displaystyle \frac{{(H^{(n-1)})}^{-1}|\mathrm{\Delta }{𝐠}^{(n)}⟩⟨\mathrm{\Delta }{𝐱}^{(n)}|+|\mathrm{\Delta }{𝐱}^{(n)}⟩⟨\mathrm{\Delta }{𝐠}^{(n)}|{(H^{(n-1)})}^{-1}}{⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩}},$

where ${(H^{(0)})}^{-1}=I$. An optimization scheme using Eq. (6) and the BFGS update formula for the inverse of an approximate Hessian matrix in the Cartesian coordinate has been implemented as ’BFGS’ in OpenMX.

4 Rational function (RF) method

The BFGS update by Eq. (14) without any care gives an ill-conditioned approximate inverse of Hessian having negative eigenvalues in many cases. This leads to the optimization to saddle points rather than the optimization to a minimum. The rational function (RF) method [4] can avoid the situation in principle. Instead of Eq. (1), we may consider the following expression:

$E$ $=$ $E_0+{\displaystyle \sum _i^{3N}}{\left({\displaystyle \frac{\partial E}{\partial x_i}}\right)}_0(x_i-x_i^{(0)})+{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i,j}^{3N}}{\left({\displaystyle \frac{{\partial }^{2}E}{\partial x_i\partial x_j}}\right)}_0(x_i-x_i^{(0)})(x_j-x_j^{(0)})+{\displaystyle \frac{1}{2}}\lambda {\displaystyle \sum _i^{3N}}{(x_i-x_i^{(0)})}^2.$ (15)

Then, the equation corresponding to Eq. (4) becomes

$(H^{(n)}+\lambda I)\mathrm{\Delta }{𝐱}^{(n)}=-{𝐠}^{(n)}.$ (16)

Therefore, a large $\lambda$ assures that $(H^{(n)}+\lambda I)$ is positive definite. If ${\lambda }^{(n)}(=-\lambda )$ is given by

${\lambda }^{(n)}=⟨{𝐠}^{(n)}|{𝐬}^{(n)}⟩.$ (17)

With Eq. (17), Eq. (16) may be equivalent to

(18)

where the size of the matrix in the left-hand side is $(3N+1)\times (3N+1)$, and called the augmented Hessian. The lowest eigenvalue of the eigenvalue problem defined by Eq. (18) may give an optimum choice for $\lambda$, and the corresponding eigenvector, the last component is scaled to 1, gives an optimization step $\mathrm{\Delta }{𝐱}^{(n)}$. In Eq. (18), the approximate Hessian can be estimated by the following BFGS formula:

$H^{(n)}$ $=$ $H^{(n-1)}+{\displaystyle \frac{|\mathrm{\Delta }{𝐠}^{(n)}⟩⟨\mathrm{\Delta }{𝐠}^{(n)}|}{⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩}}-{\displaystyle \frac{H^{(n-1)}|\mathrm{\Delta }{𝐱}^{(n)}⟩⟨\mathrm{\Delta }{𝐱}^{(n)}|H^{(n-1)}}{⟨\mathrm{\Delta }{𝐱}^{(n)}|H^{(n-1)}|\mathrm{\Delta }{𝐱}^{(n)}⟩}},$ (19)

where $H^{(0)}=I$. An optimization scheme using Eq. (18) and the BFGS update formula Eq. (19) in the Cartesian coordinate has been implemented as ’RF’ in OpenMX.

5 Eigenvector following (EF) method

By diagonalizing the approximate Hessian given by Eq. (19), the ill-conditioned situation can be largely reduced [5]. The approximate Hessian is diagonalized as

$E^{(n)}={(V^{(n)})}^T{H}^{(n)}{V}^{(n)},$ (20)

where $E^{(n)}$ is a diagonal matrix of which diagonal parts are eigenvalues of $H^{(n)}$. If the eigenvalue of the approximate Hessian is smaller than a threshold (0.02 a.u. in OpenMX3.3), the eigenvalue is set to the threshold. The modification of eigenvalues gives a corrected matrix $E^{\prime (n)}$ instead of $E^{(n)}$. Then, we have the inverse of a corrected Hessian matrix $H^{\prime (n)}$ being a positive definite as

${(H^{\prime (n)})}^{-1}=V^{(n)}{(E^{\prime (n)})}^{-1}{(V^{(n)})}^T.$ (21)

A optimization scheme using the inverse Eq. (21) in Eq. (6) in the Cartesian coordinate has been implemented as ’EF’ in OpenMX.

In addition, there are two important prescriptions for the stable optimization: (i) If $⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩$ is positive in the update of Hessian by Eq. (19), it is assured that the updated Hessian is positive definite. Therefore, if $⟨\mathrm{\Delta }{𝐱}^{(n)}|\mathrm{\Delta }{𝐠}^{(n)}⟩$ is negative, the update should not be performed. (ii) The maximum step should be always monitored, so that an erratic movement of atomic position can be avoided.