ABINIT结构优化算法分析：ionmov=22 (L-BFGS)

算法概述

L-BFGS方法在晶体几何结构优化中扮演着关键角色。优化目标是最小化体系的总能量或势能面 $E(\mathbf{R})$，其中 $\mathbf{R}$ 表示原子的多维坐标向量（结合晶格的参数构成的向量）。根据泰勒级数展开，可以利用能量关于坐标的负梯度（即原子的广义受力） $\mathbf{F}_i = -\frac{\partial E}{\partial \mathbf{R}_i}$ 或 $\mathbf{g}_k = \nabla E(\mathbf{R}_k)$ 来指导弛豫。

数学原理推导

L-BFGS算法完整迭代格式

L-BFGS算法完整的迭代格式与主要流程如下：

1. 搜索方向计算： 在第 $k$ 步迭代中，利用近似的逆Hessian矩阵 $H_k \approx (\nabla^2 E)^{-1}$ 及当前梯度寻找搜索方向 $\mathbf{p}_k$： $$ \mathbf{p}_k = -H_k \mathbf{g}_k $$

2. 线搜索与步长更新： 从点 $\mathbf{R}_k$ 沿着方向 $\mathbf{p}_k$ 做线搜索（Line Search），找到满足Wolfe收敛准则的步长 $\alpha_k$（首次试探一般设为 $\alpha_k = 1$），使得势能存在充分的下降量，随后更新系统坐标位置： $$ \mathbf{R}_{k+1} = \mathbf{R}_k + \alpha_k \mathbf{p}_k $$

3. L-BFGS Hessian矩阵历史更新准则： 分别定义两步迭代之间的位置位移差和对应的梯度变化差： $$ \mathbf{s}_k = \mathbf{R}_{k+1} - \mathbf{R}_k $$ $$ \mathbf{y}_k = \mathbf{g}_{k+1} - \mathbf{g}_k $$

全秩显式BFGS方法利用如下公式对整张逆Hessian矩阵进行秩2修正迭代： $$ H_{k+1} = \left( I - \rho_k \mathbf{s}_k \mathbf{y}_k^T \right) H_k \left( I - \rho_k \mathbf{y}_k \mathbf{s}_k^T \right) + \rho_k \mathbf{s}_k \mathbf{s}_k^T $$ 其中归一化因子 $\rho_k = \frac{1}{\mathbf{y}_k^T \mathbf{s}_k}$。

为降低存储开销，在L-BFGS中，给定一条存储历史最大长度 $m$（由ABINIT自动规划），每次仅利用最近的 $m$ 组 $(\mathbf{s}_i, \mathbf{y}_i)$ 对来隐性迭代构建矩阵乘积 $H_k \mathbf{g}_k$的搜索结果，该算法被称为双循环算法（Two-Loop Recursion），极大缩减了 $O(N^2)$ 的Hessian矩阵空间。

关键Fortran代码实现

ABINIT中有关ionmov=22的核心实现主要位于两个文件：src/45_geomoptim/m_pred_bfgs.F90（主要做ABINIT体系坐标到求解变量的转换预测）和 src/45_geomoptim/m_lbfgs.F90（深度的基础线性求根及双循环标量化求解器）。

变量提取与准备 (m_pred_bfgs.F90::pred_lbfgs)

该子程序首先整合出被优化系统的完整维度 ndim，然后将广义原子坐标输入到工作向量 vin 中，广义受力传入 vout 中。同时会生成针对底层结构的对角预估计Hessian。

!Initialise the Hessian matrix using gmet
 if (itime==1)then

   ABI_MALLOC(diag,(ndim))
   do ii=1,3*ab_mover%natom
     ! 利用度规张量初始化Hessian对角线元素估计
     diag(ii) = gmet(MODULO(ii-1,3)+1,MODULO(ii-1,3)+1)
   end do
   ...
   ! 初始化底层L-BFGS规划，历史记录数设为5 (Nocedal标配)
   call lbfgs_init(ndim,5,diag)
   ABI_FREE(diag)
...

调用L-BFGS底层求解与物理量更新 (m_pred_bfgs.F90::pred_lbfgs)

真正的预测和坐标改变都源自于对 lbfgs_execute 的执行：

!##########################################################
!### 07. Compute the next values

 vin_prev(:) = vin
 vout_prev(:) = vout
 info = lbfgs_execute(vin,etotal,vout)

在该调用过程中，旧的坐标被传入并与总能量 etotal 一起估算推演，方法执行完毕后，新的原子坐标即保留在其覆写过的 vin 向量里。

核心逆矩阵近似与线搜索迭代 (m_lbfgs.F90::lbfgs)

在最底层的 L-BFGS 处理器中，可以清晰看到著名的双循环更新算法，用于重构正交搜索方向 $\mathbf{p}_k$。

subroutine lbfgs(N,M,X,F,G,DIAG,W,IFLAG,      &
                 GTOL,STPMIN,STPMAX,STP,ITER, & ...)
     ...
!
!  COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
!  "Updating quasi-Newton matrices with limited storage",
!  Mathematics of Computation, Vol.24, No.151, pp. 773-782.
     ...
   W(N+CP) = one / YS
   W(1:N)  = -G(1:N)

   CP = POINT
   do I= 1,BOUND
     CP = CP - 1
     if (CP ==  -1) CP = M - 1
     SQ = DOT_PRODUCT(W(ISPT+CP*N+1:ISPT+CP*N+N),W(1:N))
     INMC = N + M + CP + 1
     IYCN = IYPT + CP * N
     W(INMC)= W(N+CP+1) * SQ
     W(1:N) = W(1:N) - W(INMC) * W(IYCN+1:IYCN+N)
   enddo

   W(1:N) = DIAG(1:N) * W(1:N)

   do I=1,BOUND
     ...
     W(1:N) = W(1:N) + BETA * W(ISCN+1:ISCN+N)
     ...
   enddo

而在该方向搜寻结束后，再结合调用符合条件的 Wolfe / Armijo 准则函数用于确定最佳标量步幅 $\alpha_k$（STP）：

 call MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,MAXFEV,INFO,NFEV, & ...)

代码变量与数学公式符号映射表

这套基于 Fortran 开发的 L-BFGS 代码，其命名高度契合经典数学最优化方法。对应关系说明如下：