Abinit 结构优化（离子移动）：计算流程 + 代码实现详解

基于 Abinit 10.4.7 源码分析
内部 SCF 循环部分请参考姊妹报告：abinit_scf_convergence_report.md

总体架构：双层循环

Abinit 的结构优化是一个 双层嵌套循环：

外层（离子移动循环）：itime = 1 .. ntime
  └── 内层（SCF 自洽循环）：istep = 1 .. nstep

• 内层 SCF：固定原子位置，求电子基态 → 输出 总能、力、应力
• 外层离子移动：根据力/应力 → 预测新原子位置/晶格 → 再调内层 SCF
• 当 最大力 < tolmxf 或 达到 ntime 时退出外层

调用层级总览

abinit (主程序)
└── driver (95_drive/m_driver.F90)
    └── gstateimg → gstate (95_drive/m_gstate.F90)
        └── mover (95_drive/m_mover.F90)      ← 外层离子循环驱动
            ├── scfcv_run (94_scfcv/m_scfcv.F90)  ← 内层 SCF（详见 SCF 报告）
            ├── fconv / erlxconv               ← 收敛判定
            └── precpred_1geo (95_drive/m_precpred_1geo.F90)  ← 算法分派
                └── pred_* (45_geomoptim/m_pred_*.F90)  ← 各算法实现

外层离子移动主循环（mover）

文件：src/95_drive/m_mover.F90，第 485-964 行

关键变量

主循环结构

! m_mover.F90:485
do itime = 1, ntime
  do icycle = 1, ncycle

    ! 步骤 1：对称化原子坐标
    call symmetrize_xred(...)

    ! 步骤 2：调用 SCF 内层循环（核心电子计算）
    if (need_scfcv_cycle) then
      call scfcv_run(scfcv_args, electronpositron, itimes, &
                     rhog, rhor, rprimd, xred, xred_old, conv_retcode)
    end if
    ! 输出：fcart（笛卡尔力），strten（应力），etotal

    ! 步骤 3：将力/应力存入历史
    hist%fcart(:,:,hist%ihist) = ...
    hist%strten(:,hist%ihist)  = ...
    hist%etot(hist%ihist)      = etotal

    ! 步骤 4：收敛判定
    if (specs%isFconv) then
      if (dtset%tolmxf /= 0) then
        call fconv(hist%fcart, iatfix, iexit, itime, natom, ntime, &
                   optcell, strfact, strtarget, strten, rprim, tolmxf)  ! m_mover.F90:832
      else
        call erlxconv(hist, iexit, itime, itime_hist, ntime, tolmxde)   ! m_mover.F90:844
      end if
    end if
    if (itime == ntime .and. icycle == ncycle) iexit = 1   ! 强制退出

    ! 步骤 5：算法分派 + 预测下一步几何
    call precpred_1geo(ab_mover, ab_xfh, ...)   ! m_mover.F90:859

    ! 步骤 6：从历史结构中取出新 xred, acell, rprimd
    call hist2var(acell, hist, natom, rprimd, xred)

    ! 步骤 7：晶胞改变时更新晶格度规
    if (ab_mover%optcell /= 0) then
      call matr3inv(rprimd, gprimd)
      call initylmg(...)
    end if

    if (iexit /= 0) exit   ! 跳出 icycle
  end do
  if (iexit /= 0) exit     ! 跳出 itime
end do

收敛判定

主判据：tolmxf（最大力收敛）

文件：src/95_drive/m_mover.F90，子程序 fconv（第 1085-1207 行）

! m_mover.F90:1110
fmax = zero
do iatom = 1, natom
  do idir = 1, 3
    if (iatfix(idir, iatom) /= 1) then       ! 跳过固定原子
      fmax = max(fmax, abs(fcart(idir, iatom)))
    end if
  end do
end do

! 若做晶胞优化，把应力按 strfact 归一化后并入 fmax
if (optcell == 1) fmax = max(fmax, strfact * |diag(strten)|)
if (optcell == 2 .or. 3) fmax = max(fmax, strfact * |strten|)

if (fmax < tolmxf) iexit = 1     ! 收敛！

次判据：tolmxde（连续能量稳定）

文件：src/95_drive/m_mover.F90，子程序 erlxconv（第 1223-1272 行）

if (itime_hist >= 3) then
  ediff1 = etot(ihist)      - etot(ihist_prev)
  ediff2 = etot(ihist_prev) - etot(ihist_prev2)
  if (abs(ediff1) < tolmxde .and. abs(ediff2) < tolmxde) iexit = 1
end if

仅当 tolmxf == 0 时使用。

强制退出

if (itime == ntime .and. icycle == ncycle) iexit = 1   ! m_mover.F90:820

算法选择（ionmov）

文件：src/95_drive/m_precpred_1geo.F90，第 163-211 行

分派代码

! m_precpred_1geo.F90:163
select case (ab_mover%ionmov)
  case (1)      call pred_moldyn(...)        ! 分子动力学
  case (2,3)    call pred_bfgs(...)          ! BFGS（结构优化最常用）
  case (4,5)    call pred_simple(...)        ! 共轭梯度 / 简单松弛
  case (6,7)    call pred_verlet(...)        ! Verlet MD
  case (8)      call pred_nose(...)          ! Nose-Hoover MD
  case (9)      call pred_langevin(...)      ! Langevin MD
  case (10,11)  call pred_delocint(...)      ! 离域内坐标 BFGS/CG
  case (12)     call pred_isokinetic(...)    ! 等动能 MD
  case (13)     call pred_isothermal(...)    ! 等温/等焓 MD
  case (14)     call pred_srkna14(...)       ! 辛积分
  case (15)     call pred_fire(...)          ! FIRE 快速惯性松弛
  case (16)     call pred_langevin_pimd(...) ! Langevin PIMD
  case (20)     call pred_diisrelax(...)     ! DIIS 松弛
  case (21)     call pred_steepdesc(...)     ! 最陡下降
  case (22)     call pred_lbfgs(...)         ! L-BFGS（大体系推荐）
  case (24)     call pred_velverlet(...)     ! 速度 Verlet MD
  case (25)     call pred_hmc(...)           ! 杂化 Monte Carlo
  case (28)     call ipi_pred(...)           ! i-PI 外部协议
  case default  ABI_ERROR("Wrong value of ionmov")
end select

结构优化常用算法对照

BFGS 算法实现（ionmov=2,3 最常用）

文件位置

• 分派与封装：src/45_geomoptim/m_pred_bfgs.F90（pred_bfgs 第 89 行起）
• 数学核心：src/45_geomoptim/m_bfgs.F90（hessinit, hessupdt, brdene）
• L-BFGS 数学：src/45_geomoptim/m_lbfgs.F90

状态变量打包（vin / vout）

把所有自由度打包成一维向量送进 BFGS。

vin（位置向量，长度 ndim）：

! m_pred_bfgs.F90:158
ndim = 3*natom
if (optcell == 1)       ndim = ndim + 1   ! 仅各向同性体积
if (optcell == 2 .or. 3) ndim = ndim + 6  ! 6 个应变分量
if (optcell >= 4)        ndim = ndim + 3  ! 部分晶胞自由度

打包/拆包：src/45_geomoptim/m_xfpack.F90 的 xfpack_x2vin / xfpack_vin2x。

vout（梯度向量）= 力 + 应力（取负号），由 xfpack_f2vout 打包。

Hessian 初始化

文件：src/45_geomoptim/m_bfgs.F90，子程序 hessinit（第 68-158 行）

! 默认用倒空间度规 gmet 的对角块作为单位 Hessian
hessin(:,:) = zero
do iatom = 1, natom
  do idir1 = 1, 3
    do idir2 = 1, 3
      if (iatfix(idir1,iatom)==0 .and. iatfix(idir2,iatom)==0) then
        hessin(idir1+3*(iatom-1), idir2+3*(iatom-1)) = init_matrix(idir1, idir2)
      end if
    end do
  end do
end do

! 晶胞自由度的对角元
if (optcell /= 0) then
  diag = strprecon * 30.0_dp / ucvol
  if (optcell == 1) diag = diag / 3.0_dp
  ! 加到 ndim*ndim 矩阵的右下角
end if

固定原子的 Hessian 行列直接置零，等效于禁锢自由度。

Hessian 更新（核心 BFGS 公式）

文件：src/45_geomoptim/m_bfgs.F90，子程序 hessupdt（第 191-303 行）

! 位移 / 梯度差
din(:)  = vin(:)  - vin_prev(:)        ! Δx
dout(:) = vout(:) - vout_prev(:)       ! Δg

! 屏蔽固定原子
do iatom = 1, natom
  do idir = 1, 3
    if (iatfix(idir, iatom) == 1) dout(idir + 3*(iatom-1)) = 0
  end do
end do

! H·Δg
hdelta(:) = matmul(hessin, dout)

! 标量乘积
den1 = sum(dout * din)        ! Δg^T Δx
den2 = sum(dout * hdelta)     ! Δg^T H Δg

! BFGS 修正向量
bfgs(:) = (1/den1)*din(:) - (1/den2)*hdelta(:)

! ==== B.F.G.S. 反 Hessian 更新公式（m_bfgs.F90:295-301）====
do ii = 1, ndim
  do jj = 1, ndim
    hessin(ii,jj) = hessin(ii,jj) &
                  + (1/den1)*din(ii)*din(jj) &        ! Δx Δx^T 项
                  - (1/den2)*hdelta(ii)*hdelta(jj) &  ! H Δg Δg^T H 项
                  + den2*bfgs(ii)*bfgs(jj)            ! BFGS 矫正项
  end do
end do

对应的数学公式：

$$ H_{k+1} = H_k + \frac{\Delta x ,\Delta x^T}{\Delta g^T \Delta x} - \frac{H_k\Delta g, \Delta g^T H_k}{\Delta g^T H_k \Delta g} + \beta,\mathbf{u}\mathbf{u}^T $$

步长预测（产生新原子位置）

文件：src/45_geomoptim/m_pred_bfgs.F90

对 ionmov=2（仅用力的 BFGS）（第 348-423 行）：

vin_prev(:)  = vin(:)                          ! 保存上一步位置
vin(:)       = vin(:) - matmul(hessin, vout)   ! 牛顿步：x ← x − H⁻¹·g
vout_prev(:) = vout(:)

对 ionmov=3（带总能线搜索的 BFGS）（第 425-431 行）：

call brdene(etotal, etotal_prev, hessin, ndim, vin, vin_prev, vout, vout_prev)

brdene（在 m_bfgs.F90:336-402）调用 findmin 做一维抛物线插值，找出使总能最小的位置；然后用 BFGS 步前进。这对总能可信度较高的情形（如电子收敛严，但势能面平缓）特别有效。

L-BFGS 变体（ionmov=22）

• 在 m_pred_bfgs.F90:541-876 的 pred_lbfgs 中实现
• 用 m_lbfgs.F90 模块，只存储最近 M 步的 {Δx, Δg} 对
• 内存从 O(N²) 降到 O(MN)，适合 N > 50 的大体系
• 牺牲少量收敛步数换取大量内存节省

其它结构优化算法

FIRE（ionmov=15，m_pred_fire.F90）

思想：自适应分子动力学。把原子当带阻尼的粒子，根据 v·F 的符号动态调节时间步和阻尼。

核心规则：

优点：对大体系、长链、表面收敛比 BFGS 快很多，且不存储 Hessian。

最陡下降（ionmov=21，m_pred_steepdesc.F90）

x_new = x_old - step_size * gradient

最简单但最慢。一般只用作教学或粗优化。

DIIS 松弛（ionmov=20，m_pred_diisrelax.F90）

Direct Inversion in Iterative Subspace：保留最近 diismemory 步的 {x, g} 对，找系数 c_i 使残差线性组合最小：

$$ x_{\text{new}} = \sum_i c_i , x_i, \quad \text{满足} \sum c_i = 1,\ \min |\sum c_i g_i|^2 $$

接近收敛时加速效果显著，但远离极小点不如 BFGS 稳健。

离域内坐标（ionmov=10/11，m_pred_delocint.F90）

对分子体系，用键长、键角、二面角等内坐标做优化，能跳过笛卡尔坐标下的鞍点伪迹。底层仍是 BFGS（10）或 CG（11）。

共轭梯度（ionmov=4，m_pred_simple.F90）

经典 Polak-Ribière 公式，无 Hessian，仅靠两步力信息。鲁棒但收敛步数比 BFGS 多。

晶胞优化（optcell）

晶胞自由度被打包进同一个 vin/vout 向量，与原子坐标一起用 BFGS 求极小。

关键耦合（m_mover.F90:891-920）：晶胞改变后必须重新算度规：

if (ab_mover%optcell /= 0) then
  call matr3inv(rprimd, gprimd)   ! 倒格矢
  if (psps%useylm == 1) call initylmg(...)   ! 重算球谐基
end if

力与应力的数值平衡：strfact（默认 ≈ 0.05）把应力 ×strfact 后与力同量纲比较；过大或过小会破坏 BFGS 收敛。

历史结构（abihist）

文件：src/45_geomoptim/m_abihist.F90，第 94-131 行

type abihist
  integer :: ihist          ! 当前历史索引（环形缓冲）
  integer :: mxhist         ! 最大历史步数
  real(dp) :: acell(3, mxhist)
  real(dp) :: rprimd(3, 3, mxhist)
  real(dp) :: xred(3, natom, mxhist)
  real(dp) :: fcart(3, natom, mxhist)
  real(dp) :: strten(6, mxhist)
  real(dp) :: vel(3, natom, mxhist)
  real(dp) :: vel_cell(3, 3, mxhist)
  real(dp) :: etot(mxhist), ekin(mxhist), entropy(mxhist), time(mxhist)
end type

工作机制：

• hist2var：从历史取当前几何 → 喂给 SCF
• var2hist：SCF 算完力/应力 → 写回历史
• abihist_findIndex(hist, +1)：移到下一个时间槽（环形）
• ab_xfh_type：单独保存 BFGS 状态（Hessian 矩阵）以支持 重启（restartxf 输入）

外层离子循环流程图（Mermaid）

flowchart TD
    Start([gstate → mover<br/>读 ntime, ionmov, optcell, tolmxf]) --> InitMover["初始化<br/>abimover_ini / abihist_init<br/>必要时读取重启历史"]
    InitMover --> ITime{itime = 1 .. ntime}

    ITime --> ICycle{icycle = 1 .. ncycle}
    ICycle --> Sym["1. 对称化原子坐标<br/>symmetrize_xred"]
    Sym --> SCF[/"2. <b>调用 SCF 内层循环</b><br/>scfcv_run<br/>→ fcart, strten, etotal<br/>(详见 SCF 报告)"/]
    SCF --> Store["3. 力/应力/能量 → hist 历史"]

    Store --> ConvCheck{"4. 收敛判定<br/>tolmxf 不为零?"}
    ConvCheck -- "tolmxf > 0" --> Fconv["<b>fconv</b><br/>max|F| &lt; tolmxf ?"]
    ConvCheck -- "tolmxf = 0" --> Erlxconv["<b>erlxconv</b><br/>|&Delta;E| &lt; tolmxde 连续 2 次?"]

    Fconv --> Converged{iexit == 1 ?}
    Erlxconv --> Converged

    Converged -- 否 --> Precpred["5. <b>precpred_1geo</b><br/>(算法分派)"]
    Precpred --> Algo{ionmov 值?}

    Algo -- "2, 3" --> BFGS["pred_bfgs<br/>(BFGS / 带能量线搜)"]
    Algo -- 22 --> LBFGS["pred_lbfgs<br/>(L-BFGS 大体系)"]
    Algo -- 15 --> FIRE["pred_fire<br/>(FIRE 自适应)"]
    Algo -- 20 --> DIIS["pred_diisrelax<br/>(DIIS)"]
    Algo -- 21 --> SD["pred_steepdesc<br/>(最陡下降)"]
    Algo -- "4, 5" --> CG["pred_simple<br/>(共轭梯度)"]
    Algo -- "10, 11" --> Deloc["pred_delocint<br/>(内坐标)"]
    Algo -- "1, 6-9, 12-14, 24-25" --> MD["分子动力学预测器"]

    BFGS --> Predict[/"预测新 vin → xred, acell, rprimd<br/>写入 hist"/]
    LBFGS --> Predict
    FIRE --> Predict
    DIIS --> Predict
    SD --> Predict
    CG --> Predict
    Deloc --> Predict
    MD --> Predict

    Predict --> UpdateGeo["6. hist2var<br/>取出新几何 xred, rprimd"]
    UpdateGeo --> CellUpdate{"7. optcell &ne; 0?"}
    CellUpdate -- 是 --> Metric["更新度规 gprimd<br/>重算球谐基 Ylm"]
    CellUpdate -- 否 --> EndCycle
    Metric --> EndCycle

    EndCycle --> ICycleExit{iexit &ne; 0 ?}
    ICycleExit -- 否 --> ICycle
    ICycleExit -- 是 --> ITimeExit

    Converged -- 是 --> ITimeExit
    ITimeExit{itime == ntime<br/>or 收敛?} -- 否 --> ITime
    ITimeExit -- 是 --> Finalize["最终输出<br/>写 HIST / XML / netcdf"]
    Finalize --> EndNode([END])

    classDef coreStep fill:#fde68a,stroke:#b45309,stroke-width:2px,color:#000;
    classDef decision fill:#bfdbfe,stroke:#1d4ed8,color:#000;
    classDef terminal fill:#d1fae5,stroke:#047857,color:#000;
    classDef algo fill:#fce7f3,stroke:#9d174d,color:#000;
    class SCF,Sym,Store,Precpred,UpdateGeo,Metric,Predict,Finalize,InitMover coreStep;
    class ITime,ICycle,ConvCheck,Algo,Converged,CellUpdate,ICycleExit,ITimeExit decision;
    class Start,EndNode terminal;
    class BFGS,LBFGS,FIRE,DIIS,SD,CG,Deloc,MD,Fconv,Erlxconv algo;

BFGS 单步算法流程图（Mermaid）

flowchart TD
    A([进入 pred_bfgs<br/>(itime 步)]) --> B["xfpack_x2vin<br/>打包: xred,acell,rprim &rarr; vin<br/>(ndim = 3&middot;natom + 晶胞自由度)"]
    B --> C["xfpack_f2vout<br/>打包: fcart,strten &rarr; vout (梯度)"]
    C --> D{itime == 1 ?}
    D -- 是 --> E["<b>hessinit</b><br/>初始化 H = gmet 张量积<br/>+ 晶胞对角块"]
    D -- 否 --> F["<b>hessupdt</b><br/>用 (&Delta;x, &Delta;g) 更新 H<br/>BFGS 修正公式"]
    E --> G{ionmov 值?}
    F --> G

    G -- "ionmov = 2" --> H["Newton 步:<br/>vin_new = vin &minus; H &middot; vout"]
    G -- "ionmov = 3" --> I["<b>brdene</b><br/>能量抛物线插值 (findmin)<br/>+ Newton 步"]

    H --> J["xfpack_vin2x<br/>拆包 vin_new &rarr; 新 xred, acell, rprim"]
    I --> J
    J --> K["写入 hist<br/>(下一次外层迭代用)"]
    K --> Z([返回 mover])

    classDef step fill:#fde68a,stroke:#b45309,stroke-width:2px,color:#000;
    classDef decision fill:#bfdbfe,stroke:#1d4ed8,color:#000;
    classDef terminal fill:#d1fae5,stroke:#047857,color:#000;
    class B,C,E,F,H,I,J,K step;
    class D,G decision;
    class A,Z terminal;

关键文件汇总表

一句话总览

Abinit 的结构优化由外层 mover（src/95_drive/m_mover.F90）驱动 itime 循环：每步先调内层 scfcv_run 算电子基态拿到力/应力（详见 SCF 报告），再用 fconv 检查 max|F| < tolmxf；不收敛则由 precpred_1geo 按 ionmov 分派到具体预测器（BFGS=m_pred_bfgs.F90，FIRE=m_pred_fire.F90，DIIS=m_pred_diisrelax.F90 等）算出新坐标。BFGS 的核心数学在 m_bfgs.F90 中：hessupdt 用 Δx 和 Δg 按标准 BFGS 公式更新反 Hessian，pred_bfgs 用 vin ← vin − H·vout 一步牛顿更新位置；晶胞自由度通过 optcell 与 xfpack 打包进同一向量统一优化。

References

Bitzek, Erik, et al. "Structural Relaxation Made Simple." Physical Review Letters, vol. 97, no. 17, 2006, p. 170201. link

Broyden, C. G. "The Convergence of a Class of Double-Rank Minimization Algorithms." IMA Journal of Applied Mathematics, vol. 6, no. 1, 1970, pp. 76–90. link

Fletcher, R. "A New Approach to Variable Metric Algorithms." The Computer Journal, vol. 13, no. 3, 1970, pp. 317–322. link

Goldfarb, Donald. "A Family of Variable-Metric Methods Derived by Variational Means." Mathematics of Computation, vol. 24, no. 109, 1970, pp. 23–26. link

Gonze, Xavier, et al. "The Abinit Project: Impact, Environment and Recent Developments." Computer Physics Communications, vol. 248, 2020, p. 107042. link

Nocedal, Jorge. "Updating Quasi-Newton Matrices with Limited Storage." Mathematics of Computation, vol. 35, no. 151, 1980, pp. 773–782. link

Pulay, Peter. "Convergence Acceleration of Iterative Sequences. The Case of SCF Iteration." Chemical Physics Letters, vol. 73, no. 2, 1980, pp. 393–398. link

Shanno, D. F. "Conditioning of Quasi-Newton Methods for Function Minimization." Mathematics of Computation, vol. 24, no. 111, 1970, pp. 647–656. link

The Abinit Code. Version 10.4.7, 2024. link