General optimization: http://www.sfu.ca/~ssurjano/optimization.html
For most electronic structure methods, the potential energy as a function of nuclear coordinates is not analytically available. To optimise the structure, the only feasible way is to keep trying different geometries until we are quite certain that the structure with the lowest energy has been found. Different geometry optimization methods have different ways to change geometries. Since the potential energy is a multi variable function, the most wanted optimization strategy aims at moving structure in the direction lowering the energy with respect to all coordinates. This direction is the opposite of the gradient vector, $\bold{g}$, which is defined as: $$ \bold{g} = \begin{bmatrix} \frac{\partial E}{\partial x _ 1} \\ \frac{\partial E}{\partial x _ 2} \\ \dots \\ \frac{\partial E}{\partial x _ n} \end{bmatrix} $$
Suppose the structure to be optimized has current position $\bold{x^k}$, which is a function has $n$ variables ($n=3N_{atoms}-6(5)$). The structure optimization method tries to place the structure in a new position $\bold{x^{k+1}}$ using the displacement $\Delta\bold{x}=\bold{x^{k+1}}-\bold{x^k}$. The energy at new position can be written in Taylor expansions truncating at second order: $$ E(\bold{x^{k+1}}) = E(\bold{x^k}) + (\frac{\partial E(\bold{x^k})}{\partial x^k_1} \Delta x_1 + \frac{\partial E(\bold{x^k}) }{\partial x^k_2} \Delta x_2 + \dots) + \frac{1}{2} (\frac{\partial^2 E(\bold{x^k})}{\partial {x^k_1}^2} \Delta x^2_1 + \frac{\partial^2 E(\bold{x^k})}{\partial {x^k_2}^2} \Delta x^2_2 + \dots) $$ In the matrix form, $E(\bold{x^{k+1}})$ is written as: $$ E(\bold{x^{k+1}}) = E(\bold{x^k}) + (\bold{g^k})^T \Delta\bold{x} + \frac{1}{2} \Delta\bold{x}^T \bold{H}^k \Delta\bold{x} $$ The optimization goal is to make the new structure a local minimum on PES of which the first derivative is zero: $\frac{\partial E(\bold{x^{k+1}})}{\partial {x^{k+1}}_i} = 0$. Applying this stationary condition to the Taylor expansions of $E(\bold{x^{k+1}})$: $$ \begin{aligned} \frac{\partial E(\bold{x^{k+1}})}{\partial {x^{k+1}}_1} = 0 + ( \frac{\partial E(\bold{x^k})}{\partial x^k_1} * 1 + 0 + \dots) + \frac{1}{2} (\frac{\partial^2 E(\bold{x^k})}{\partial {x^k_1}^2} * 2 \Delta x_1 + 0 + \dots) &= 0 \\ \frac{\partial E(\bold{x^{k+1}})}{\partial {x^{k+1}}_2} = 0 + (0 + \frac{\partial E(\bold{x^k})}{\partial x^k_2} * 1 + \dots) + \frac{1}{2} (0 + \frac{\partial^2 E(\bold{x^k})}{\partial {x^k_2}^2} * 2 \Delta x_2 + \dots) &= 0 \\ &\dots \\ \end{aligned} $$ To clean this up: $\frac{\partial E(\bold{x^k})}{\partial x^k_i} + \frac{\partial^2 E(\bold{x^k})}{\partial {x^k_i}^2} \Delta x_i = 0 \text{ (for i in 1, 2, \dots, n)}$. In the martrix form, it is written as: $$ \begin{aligned} \bold{g}^k + \bold{H}^k \Delta{\bold x} &= 0 \\ \Delta{\bold x} &= - {(\bold{H}^k)}^{-1} \bold{g}^k \end{aligned} $$
Quality of Newton-Raphson-based optimization:
在一个给定的压强下，结构优化可以以焓值（H = E + PV）最小化为目标。
很多DFT的程序支持这样功能，比如VASP（用PSTRESS来控制压强，单位为kBar，10kBar = 1GPa， 1GPa1Å$^3$=6.24x10$^{-3}$eV）。
- [2017 - PCCP - Pressure-driven phase transition mechanisms revealed by quantum chemistry: l-serine polymorphs]：在给定压强下，原子位置和晶格参数都可以优化（ISIF=3）
- GPAW从1.4.0之后开始支持Pulay stress，只能在plane wave下使用，GPAW(mode=PW(ecut, pulay_stress=...), ...)，单位是eV/Å^3。只有stress tensor里加了Pulay correction（加在diagonal term上），energy里没有加Pulay correction（pressurevolume）。只有调用了atoms.get_stress()，才会计算stress tensor