ABINIT结构优化`ionmov=2`算法：BFGS方法实现与分析

算法概述

从ABINIT最新代码架构来看，ionmov=2核心流程主要由 src/45_geomoptim/m_pred_bfgs.F90 模块中的 pred_bfgs 子程序主导位移预测，并通过 src/45_geomoptim/m_bfgs.F90 中的 hessupdt 来进行逆Hessian矩阵的实时更新。

数学原理推导

在结构优化过程中，核心目标是求取能够使得系统基态总能量 $E(\mathbf{R})$ 达到局部极小值的原子坐标向量集 $\mathbf{R}$。

定义全部原子所构成的位移向量的第 $k$ 步迭代为 $\mathbf{R}_k$。相应的，根据物理学定律，该系统在第 $k$ 步的原子力（Force）即相当于能量对该处坐标位置的负梯度：

$$ \mathbf{F}_k = -\frac{\partial E}{\partial \mathbf{R}_k} = -\nabla E(\mathbf{R}_k) $$

为方便数学表达及数值演化，定义体系无约束状态下受限的广义梯度向量 $\mathbf{G}_k = -\mathbf{F}_k = \nabla E(\mathbf{R}_k)$。

根据牛顿法的非线性优化迭代律：

$$ \mathbf{R}_{k+1} = \mathbf{R}_k - \mathcal{H}^{-1}_k \mathbf{G}_k $$

这里 $\mathcal{H}_k = \nabla^2 E(\mathbf{R}_k)$ 是真实的结构Hessian矩阵。然而在第一性原理物理计算中求解真实的Hessian极为缓慢。故在BFGS算法中，我们不去精确求取 $\mathcal{H}^{-1}_k$，而是构造一个正定对称的近似逆Hessian矩阵 $\mathbf{H}_k \approx \mathcal{H}^{-1}_k$。

迭代步长与坐标更新预测

每一次结构迭代步中，系统首先根据当前储存的近似逼近求出下一个预测的原子位置：

$$ \mathbf{R}_{k+1} = \mathbf{R}_k - \mathbf{H}_k \mathbf{G}_k = \mathbf{R}_k + \mathbf{H}_k \mathbf{F}_k $$

逆Hessian矩阵的BFGS修正更新法则

假定在第 $k$ 到迭代发生后，程序记录下了前后两步的信息。令第 $k$ 步产生的位移差和梯度变化量分别为：

$$ \Delta \mathbf{R} = \mathbf{R}_{k+1} - \mathbf{R}_k $$

$$ \Delta \mathbf{G} = \mathbf{G}_{k+1} - \mathbf{G}_k = -(\mathbf{F}_{k+1} - \mathbf{F}_k) $$

BFGS方法的核心在线性迭代地加上特定非对角约束矩阵修正 $\mathbf{H}$，其满足典型的割线方程条件：$\mathbf{H}_{k+1} \Delta \mathbf{G} = \Delta \mathbf{R}$。在ABINIT的代码实现中，针对此式给出了一种组合分解方案：

定义一些与能量地形曲率相关的标量归一化内积系数：

$$ D_1 = \frac{1}{\Delta \mathbf{G}^T \Delta \mathbf{R}} $$

$$ D_2 = \frac{1}{\Delta \mathbf{G}^T (\mathbf{H}_k \Delta \mathbf{G})} $$

并通过提取差异张量构造一个中间差分修正向量 $\mathbf{u}$：

$$ \mathbf{u} = D_1 \Delta \mathbf{R} - D_2 (\mathbf{H}_k \Delta \mathbf{G}) $$

经过化简重构，最终ABINIT的具体矩阵元素更新解析数学公式写为如下对称形式：

$$ \mathbf{H}_{k+1} = \mathbf{H}_k + D_1 (\Delta \mathbf{R})(\Delta \mathbf{R})^T - D_2 (\mathbf{H}_k \Delta \mathbf{G}) (\mathbf{H}_k \Delta \mathbf{G})^T + \frac{1}{D_2} \mathbf{u} \mathbf{u}^T $$

此过程反复循环。每次迭代更新后，新的 $\mathbf{H}_{k+1}$ 被平滑应用到下一次坐标迭代上。当系统最大的原子力（力图对应的梯度绝对值）低于预设定的判据阈值（在ABINIT中该参量名为 tolmxf）时，系统立刻结束迭代，判定物理结构已经成功处于收敛态。

关键Fortran代码实现与公式映射

经调取ABINIT最新的稳定版源代码档案进行检索配对（主要集中在 m_pred_bfgs.F90 和 m_bfgs.F90 中）相关实现提取如下。

变量名与数学符号对应映射表

为了帮助读懂源代码，以下罗列底层Fortran变量和上述算法推导的数学公式的符号对应关系：

原位Fortran核心实现源码片段展示

代码块1：坐标更新与期望预测

在主功能模块 pred_bfgs（位于 src/45_geomoptim/m_pred_bfgs.F90）中执行着核心的坐标空间预测迭代（映射公式 $\mathbf{R}_{k+1} = \mathbf{R}_k - \mathbf{H}_k \mathbf{G}_k$）：

! --- 源自 src/45_geomoptim/m_pred_bfgs.F90 的核心 BFGS 计算 ---
  if(ionmov==2 .or. itime==1)then

!  Previous cartesian coordinates (保存 R_k 到上一步变量空间)
   vin_prev(:)=vin(:)

!  New atomic cartesian coordinates are obtained from vin, hessin and vout
!  该行即对应方程：R_{k+1} = R_k - H_k * G_k，利用 matmul 作用 Hessian 到受力之上
   vin(:) = vin(:) - matmul(hessin(:,:), vout(:))

!  Previous atomic forces (保存当步受力/梯度为下步作为上一次记录)
   vout_prev(:)=vout(:)

代码块2：Hessian逆矩阵实时的重整及修正

若当前迭代步并未完全收敛，程序将调动 hessupdt （位于 src/45_geomoptim/m_bfgs.F90）根据刚才演化步骤拿到的增量数据反推构造下一阶段收敛更快的近似海森矩阵：

! --- 源自 src/45_geomoptim/m_bfgs.F90 的 hessupdt 更新函数 ---
subroutine hessupdt(hessin, iatfix, natom, ndim, vin, vin_prev, vout, vout_prev, nimage)
  ...
! Difference between new and previous vectors (求取变动增量 \Delta R 和 \Delta G)
  din(:) = vin(:) - vin_prev(:)       ! 相当于数学推导中的 \Delta \mathbf{R}
  dout(:)= vout(:) - vout_prev(:)     ! 相当于数学推导中的 \Delta \mathbf{G}
  ...
! Compute approximate inverse Hessian times delta fcart
! 求解 \mathbf{H}_k \Delta \mathbf{G} 项
  hdelta(:)=zero
  do ii=1,ndim
    hdelta(:)=hdelta(:) + hessin(:,ii)*dout(ii)
  end do

! Calculation of dot products for the denominators
! 分别映射 D_1 与 D_2 归一因子内面的求和点乘部分
  den1 = sum(dout(1:ndim) * din(1:ndim))
  den2 = sum(dout(1:ndim) * hdelta(1:ndim))

! Denominators are multiplicative (转化为倒数系数)
  den1=one/den1
  den3=one/den2

! Vectors which make a difference between the BROYDEN and the DAVIDON scheme.
! 构造中间调节增项向量 \mathbf{u} (代码内名为 bfgs)
  bfgs(:)=den1*din(:) - den3*hdelta(:)

! B.F.G.S. updating formula (根据数学公式进行 H_{k+1} 二维矩阵级别的加项更新)
  do ii=1,ndim
    do jj=1,ndim
      hessin(ii,jj) = hessin(ii,jj) + den1*din(ii)*din(jj) &
     &     - den3*hdelta(ii)*hdelta(jj) + den2*bfgs(ii)*bfgs(jj)
    end do
  end do
end subroutine hessupdt

补充注意事项与边界条件

• 原子的约束锁定处理 (Atom fixing constraints): 在实际代码执行中，为防止假性受力更新，程序定义了约束数组 iatfix(idir, iatom)==1。在BFGS对更新的判断期内，一旦涉及固定方向，受力变差会自动被置零(dout(idir)=0)与结构位移差清除（vin=vin_prev），屏蔽无效迭代提升速度。
• 混合策略预制: 前述的更新范式展现了典型的单步BFGS行为。为了增强健壮性和预防震荡，遇到特殊或复杂迭代地形时ABINIT也支持对 vin 应用更为复杂稳定的Pulay Multi-step混合修正补偿方案。

References

• ABINIT官方文档link