Abinit occopt：占据数选项完整解析

基于 Abinit 10.4.7 源码分析 核心模块：

• src/44_abitypes_defs/m_dtset.F90（占据数初始化 dtset_initocc_chkneu）
• src/57_iovars/m_invars1.F90（fband → nband 自动计算）
• src/57_iovars/m_chkinp.F90（occ 一致性检查）
• src/61_occeig/m_occ.F90（占据更新 newocc, getnel, occ_fd 等）
• src/79_seqpar_mpi/m_vtorho.F90（SCF 中的固定/可变占据分支）
• src/67_common/m_mkrho.F90（占据数 × 波函数 → 密度）
• abimkdocs/variables_abinit.py（变量定义）

引言：占据数在 DFT 中的角色

Kohn-Sham 密度的构造是

$$\rho(\mathbf{r}) = \sum_{\sigma}\sum_{n,\mathbf{k}} w_\mathbf{k},f_{n\mathbf{k}\sigma},|\psi_{n\mathbf{k}\sigma}(\mathbf{r})|^2$$

其中 $f_{n\mathbf{k}\sigma} \in [0, f_{\max}]$ 是占据数。它必须满足全电子数约束

$$\sum_{\sigma n \mathbf{k}} w_\mathbf{k}, f_{n\mathbf{k}\sigma} = N_{\text{elect}}$$

总能（变分）

$$E_{\text{tot}}[\rho, {f}] = \sum w_\mathbf{k},f_{n\mathbf{k}},\epsilon_{n\mathbf{k}} - E_{\text{dc}}[\rho] - T,S[{f}]$$

最后一项熵贡献只在有限温（金属型 smearing）时出现。

occopt 决定了 ${f_{n\mathbf{k}\sigma}}$ 如何确定：

1. 由用户手动给定（occopt = 0, 2）
2. 由代码按"半导体满填"自动给定（occopt = 1）
3. 由 Fermi-Dirac / Gaussian / Cold smearing 等公式根据当前本征值动态计算（occopt = 3..8）
4. 双 Fermi 能级，强制激发态（occopt = 9）

下文按"输入参数 → 物理原理 → 公式 → 代码实现 → 一步步流程"展开。

occopt 全部 10 种取值速查

判定标准：

• occopt < 3 → 固定占据（fixed_occ = .true.）：SCF 期间不调 newocc
• occopt >= 3 → 金属型（可变）：每次 SCF 步都调 newocc

关键代码（src/79_seqpar_mpi/m_vtorho.F90:527）：

fixed_occ = (dtset%occopt < 3 .or. electronpositron_calctype(electronpositron) == 1)

与占据相关的输入参数总览

下表收齐了所有与 occopt 直接/间接相关的输入变量（定义位置 abimkdocs/variables_abinit.py）：

nelect 的计算

代码位置：src/44_abitypes_defs/m_dtset.F90:1199-1208：

zval = 0
do iatom = 1, natom
   zval = zval + ziontypat(typat(iatom))
end do
nelect = zval - (cellcharge - nelectjell)

其中 ziontypat 来自每种伪势文件头部的 Zion（价电子数）。

maxocc 与每条带容量

maxocc = 2 / (nsppol * nspinor)

• nsppol=1, nspinor=1 → maxocc=2（每条带可装 2 个电子，自旋简并）
• nsppol=2, nspinor=1 → maxocc=1（自旋极化，每条带 1 个）
• nsppol=1, nspinor=2 → maxocc=1（非共线，自旋升降同等地位）
• nsppol=2, nspinor=2 → maxocc=0.5（理论上禁止，代码会报错）

占据数初始化：dtset_initocc_chkneu

代码位置：src/44_abitypes_defs/m_dtset.F90:1181-1444，是整个 SCF 启动前的占据数主入口。它做三件事：

1. 计算 nelect：核电荷数减去 cellcharge
2. 初始化 occ_orig（仅当 occopt = 1 或 3..9 时；occopt = 0, 2 时占据数由用户手动给）
3. 检查电荷平衡：Σ wtk × occ ≈ nelect，偏差超 tol8 报错

半导体填充算法（occopt = 1 或 3..9 的初始猜测）

代码 m_dtset.F90:1216-1289：

if (occopt == 1 .or. (occopt >= 3 .and. occopt <= 9)) then

   maxocc = 2.0 / (nsppol * nspinor)

   ! 全占满的能带数
   nocc = int((nelect - nh_qFD - 1e-8) / maxocc) + 1

   ! 最高占据带的"残量"（可能 < maxocc）
   occlast = nelect - nh_qFD - maxocc * (nocc - 1)

   ! 检查 nband 足够
   if (nocc <= nband(1) * nsppol .or. iscf == -2) then

      ! 情况 A: 无磁矩约束
      if (spinmagntarget < -99.98) then
         tmpocc(1 : nocc-1) = maxocc        ! 全满
         tmpocc(nocc)      = occlast        ! 最后一带部分填
         tmpocc(nocc+1 :)  = zero           ! 空带

      ! 情况 B: 固定磁矩 (nsppol=2)
      else
         do isppol = 1, nsppol
            sign_spin = 3 - 2*isppol         ! +1 spin-up, -1 spin-down
            nelect_spin = 0.5 * (nelect * maxocc + sign_spin * spinmagntarget)
            nocc = ceiling(nelect_spin / maxocc)
            occlast = nelect_spin - (nocc-1)*maxocc
            occ(1:nocc-1, isppol) = maxocc
            occ(nocc,    isppol) = occlast
            occ(nocc+1:, isppol) = 0
         end do
      end if
   end if
end if

示例：nelect=8 (Si), nsppol=1, nspinor=1：

• maxocc = 2
• nocc = int((8-0-1e-8)/2) + 1 = 4 → 4 条全满带
• occlast = 8 - 2*3 = 2 → 第 4 带也是满的
• occ_orig = [2, 2, 2, 2, 0, 0, 0, 0]（如果 nband=8）

有半填的例子：nelect=7 (Al, nsppol=1)：

• nocc = int(7/2-eps) + 1 = 4
• occlast = 7 - 2*3 = 1 → 第 4 带半填
• occ_orig = [2, 2, 2, 1, 0, 0, 0]

这种"非 0 非满"的半填会触发 occopt=1 报错（除非是 H 原子），告诉用户改用 occopt=7。检查代码 m_chkinp.F90:2780-2792。

occopt=9 的双 Fermi 分支

m_dtset.F90:1245-1268，先填价带（留 nh_qFD 个空穴）：

if (occopt == 9) then
   ! 填到 ivalence-1 满，第 ivalence 部分填到 nelect-nh_qFD
   tmpocc(nocc+1 : ivalence*nsppol) = zero    ! 价带顶部留空穴
   nocc2 = (ne_qFD - 1e-8) / maxocc + 1
   occlast2 = ne_qFD - maxocc * (nocc2 - 1)
   ! 然后在导带（ivalence+1..）填 ne_qFD 个电子
   tmpocc(ivalence*nsppol + 1 : ivalence*nsppol + nocc2 - 1) = maxocc
   tmpocc(ivalence*nsppol + nocc2)                          = occlast2
end if

得到例如：occ = [2, 2, 2, 1.5, 0, ..., 0, 0.5, 0, ...]（价带顶部 0.5 个空穴，导带底部 0.5 个电子）。

电荷平衡检查

代码 m_dtset.F90:1380-1442：

nelect_occ = 0
do isppol, ikpt, iband
   nelect_occ = nelect_occ + wtk(ikpt) * occ_orig(...)
end do

if (|nelect_occ - nelect| > tol11) → 警告
if (|nelect_occ - nelect| > tol8)  → 错误并停止

注意：

• occopt = 0, 2：用户给定占据。代码核对总数
• occopt = 1, 3..9：代码刚生成占据，这步检查只是个 sanity check
• occopt = 2 特殊：wtk 不自动归一化，必须用户保证 $\sum w_k = 1$（其他都自动归一化）

nband 的自动决定（与 fband 联动）

代码位置：src/57_iovars/m_invars1.F90:1933-2021。当用户没指定 nband（且 occopt 不是 0 或 2），由 fband 自动算出：

! fband 默认值
if (occopt == 1) fband = 0.125    ! 绝缘体只需少量空带
else            fband = 0.5       ! 金属需要更多空带做 smearing

! 自动 nband 公式
zelect = zval - cellcharge_min
mband_upper = nspinor * ( (ceiling(zelect-tol10) + 1) / 2          &
                          + ceiling(fband*natom - tol10) )         &
            + (nsppol - 1) * ceiling(half * sum_spinat)

nband(:) = mband_upper

含义：

1. 基础：填满 zelect 个电子所需的最小 band 数 ≈ ceiling(zelect/2)
2. 缓冲：再加 fband × natom 条空带（默认绝缘体 12.5%，金属 50%）
3. 自旋极化补丁：nsppol=2 时额外加 spinat 总和的一半

示例：Si (8 电子, natom=1)：

• occopt=1：mband_upper = 1 * (4 + ceiling(0.125*1)) = 5
• occopt=7：mband_upper = 1 * (4 + ceiling(0.5*1)) = 5

示例：Al (3 电子, natom=1)：

• occopt=7：mband_upper = 1 * (2 + 1) = 3

用户可以显式给 nband，覆盖 fband 默认，对金属强烈建议显式设大（例如 1.5 倍占据 band 数）。

occopt = 0：手动占据

输入要求

occopt 0
nband 8                          # 单标量
nkpt 10
nsppol 1
occ                              # nband*nsppol = 8 个数
  2.0 2.0 2.0 2.0  0.0 0.0 0.0 0.0
wtk                              # nkpt = 10 个数
  0.5  1.0 ... ...

所有 k 点共享同一组 occ，wtk 由用户保证 $\sum w_k = 1$。

代码处理

m_chkinp.F90:2743-2769 做：

if ((iscf > 0 .or. iscf == -1 .or. iscf == -3) .and. (occopt == 0 .or. occopt == 2)) then
   sumocc = 0
   do iband over all bands:
      if (occ < -tol8) → 报错"负占据"
      sumocc += occ
   end do
   if (sumocc < 1e-8) → 报错"未定义"
end if

进一步通过 dtset_initocc_chkneu 检查 $\sum w_k \cdot \text{occ} = N_{\text{elect}}$。

应用场景

• 强制磁矩高自旋态（如 Cr+5+ d¹ 配置）
• 验证波函数耦合（特殊几何）
• 与外部代码做对比测试

occopt = 1：半导体自动填充（默认值）

触发条件

需满足"整数填充"：(nelect - nh_qFD) 必须正好等于 maxocc * n 形式（n 整数），否则代码报错。

算法

dtset_initocc_chkneu 第 4.1 节算法的退化情形（spinmagntarget = -99.99）。

自旋极化时的强制 spinmagntarget

m_chkinp.F90:2780-2792：当 nsppol=2 且 occopt=1 时，必须显式给 spinmagntarget（除非是 H 原子）：

if (nsppol == 2 .and. occopt == 1 .and. abs(spinmagntarget+99.99) < tol8) then
   if (natom /= 1 .or. |znucl-1| > tol8) then
      ABI_ERROR("In nsppol=2 + occopt=1 case, must specify spinmagntarget. ...")
   end if
end if

原因：纯填充无法决定 spin-up / spin-down 各占多少；需要用户告诉系统期望的磁矩。

反铁磁的特殊处理

反铁磁体（如 NiO）用 nsppol=1, nspden=2：自旋密度分量在同一通道内交替分布，不需要 spinmagntarget。occopt=1 在这种情形下也能用。

occopt = 2：每 k 点自由占据

输入要求

occopt 2
nband     12 8 10 8 ...           # nkpt*nsppol 个值
occ       2.0 2.0 ...             # Σ_ikpt nband(ikpt) * nsppol 个数
wtk       0.05 0.10 ...           # NOT 自动归一化

每 k 点的 band 数可以不一样（其他 occopt 都要求 nband 是标量）。

代码分支

src/57_iovars/m_invars1.F90:2009-2017：

else if (occopt == 2) then
   ABI_MALLOC(reaalloc, (nkpt*nsppol))
   call intagm(reaalloc, nband, ..., 'nband', tnband, 'INT')
   if (tnband == 1) then
      do ikpt = 1, nkpt*nsppol
         if (nband(ikpt) > mband_upper) mband_upper = nband(ikpt)
      end do
   end if
end if

vtorho 在 k 点循环时有特殊分支（m_vtorho.F90:215, 267 等）：

if (occopt == 2) high_band_index = nband(ikpt + nkpt*(isppol-1))

应用场景

• 强迫某 k 点上的占据 (响应函数测试)
• 在特定 k 点用更多 band 提高精度
• 历史用法：早期版本算 SOC 的 workaround

不推荐普通用户使用，复杂且易错。

occopt = 3..9：可变占据（smearing）

这部分已在 abinit_smearing_report.md 中详细写过。这里仅给一个与本报告焦点相关的浓缩。

八种 smearing 核函数

设 $x = (\mu - \epsilon)/k_BT$，下表给出展宽 δ 函数 $\tilde{\delta}(x)$ 与对应占据数函数 $f(x) = \int_{-\infty}^x \tilde{\delta}(x')dx'$：

代码位置：src/61_occeig/m_occ.F90:1116-1168。所有公式被预先在 12001 点的查找表里离散化，运行时通过 splfit 三次样条插值得到 $f$、$s$、$df/dx$。

SCF 中的占据更新

每次 SCF 步在 vtorho 里：

1. 求解 $H\psi = \epsilon\psi$ 得新本征值 ${\epsilon_{n\mathbf{k}}}$
2. 调 newocc 用二分法（最多 120 次）找 Fermi 能级 $\mu$ 使 $\sum w_\mathbf{k} f_{n\mathbf{k}} = N_{\text{elect}}$
3. 更新 occ(:) 和 doccde(:)（占据对能量的导数）
4. 用新 occ 累加密度 $\rho(\mathbf{r})$

代码片段（m_vtorho.F90:1365-1369）：

call newocc(doccde, eigen, energies%entropy_ks, energies%e_fermie, energies%e_fermih, &
            ivalence, spinmagntarget, mband, nband, nelect, ne_qFD, nh_qFD, &
            nkpt, nspinor, nsppol, occ, occopt, prtvol, tphysel, tsmear, wtk, &
            prtstm=prtstm, stmbias=stmbias, extfpmd=extfpmd)

熵进入总能

occopt ∈ {3..8} 时（不含 9），总能多一项 $-TS$：

if (occopt >= 3 .and. occopt <= 8) then
   if (|tphysel| < tol10) then
      e_entropy = -tsmear * entropy
   else
      e_entropy = -tphysel * entropy
   end if
end if

这就是 Mermin 自由能 $F = E - TS$。代码 src/67_common/m_dft_energy.F90:1148-1156。

占据与密度的耦合（mkrho）

代码位置：src/67_common/m_mkrho.F90。每个 SCF 步收尾时由本子程序把占据数 × 波函数 → 密度：

do isppol, ikpt:
   weight = wtk(ikpt) / ucvol
   do iband = 1, nband_k:
      if (|occ(iband)| > tol8) then         ! 跳过零占据的 band
         weight_t(nband_occ) = occ(iband) * wtk(ikpt) / ucvol
         cwavef = cg(iband)
         nband_occ = nband_occ + 1
      end if
   end do

   ! 用 FFT 把 ψ(G) → ψ(r)，并累加 |ψ|² × weight 到 rhoaug
   call fourwf(option=1, rhoaug, cwavef, ..., weight_array_r=weight_t, ...)
end do

关键点：

• 零占据带不进入密度（节省 FFT）
• 每个 band 的权重 = occ × wtk / ucvol（这就是 BZ 积分的离散权重）
• PAW 还要额外加上"球内"贡献 rhoij，由 pawmkrhoij 处理

最后跨 MPI 进程求和（k 点并行）：

call xmpi_sum(rhoaug, mpi_comm_kpt, ierr)

占据与本征值能量

代码位置：src/61_occeig/m_occ.F90 调用 e_eigen（src/79_seqpar_mpi/m_vtorho.F90:727）：

e_eigenvalues = 0
do isppol, ikpt, iband:
   e_eigenvalues = e_eigenvalues + wtk(ikpt) * occ(iband) * eigen(iband)
end do

类似地，动能、非局域能、Fock 能等都用同样权重 wtk * occ 加和（m_vtorho.F90:1165-1171）：

energies%e_kinetic     += wtk(ikpt) * occ_k(iband) * ek_k(iband)
energies%e_nlpsp_vfock += wtk(ikpt) * occ_k(iband) * enlx_k(iband)
energies%e_fock        += 0.5 * wtk(ikpt) * occ_k(iband) * fock_eigen(iband)

occopt 对 SCF 流程的影响图

下图把 occopt 在三大阶段的分支条件用 Mermaid 流程图表达，便于直观对比"固定占据"（occopt < 3）与"可变占据"（occopt ≥ 3）两条主路径。

flowchart TD
    Start([输入文件 + 伪势]) --> A1

    subgraph PhaseInit["阶段一：初始化（istep = 0）"]
        direction TB
        A1["<b>chkinp</b><br/>检查 occ, wtk, fband, occopt 一致性"]
        A2["<b>invars1</b><br/>由 fband 自动决定 nband<br/>mband = ceiling(zelect/2) + ceiling(fband·natom)"]
        A3["<b>dtset_initocc_chkneu</b><br/>nelect = zval − cellcharge<br/>maxocc = 2 / (nsppol·nspinor)"]
        A4{"occopt 分支"}
        A5["验证用户给的 occ<br/>(总和需 = nelect)"]
        A6["自动生成 occ_orig<br/>nocc = int(nelect/maxocc) + 1<br/>occlast = nelect − maxocc·(nocc−1)"]
        A7["电荷检查：<br/>|Σ wtk·occ − nelect| < tol8"]
        A1 --> A2 --> A3 --> A4
        A4 -- "occopt = 0 或 2<br/>(固定占据，用户输入)" --> A5 --> A7
        A4 -- "occopt = 1 或 3..9<br/>(半导体填充作为初值)" --> A6 --> A7
    end

    A7 --> B1

    subgraph PhaseSCF["阶段二：SCF 主循环（istep = 1..nstep）— vtorho"]
        direction TB
        B1["fixed_occ = (occopt &lt; 3)"]
        B2["对每个 k 点：<br/>solve H ψ = ε ψ → eigen(k)<br/>(CG / LOBPCG / ChebFi)"]
        B3{"fixed_occ?"}
        B4["跳过占据更新<br/>occ 保持初始化时的值<br/>(occopt ≤ 2)"]
        B5["<b>call newocc</b><br/>二分法找 μ，重算 occ<br/>(occopt ≥ 3)"]
        B6["<b>get_entropy</b><br/>S = Σ wtk · s(x)"]
        B7["<b>e_eigen</b>：<br/>E_band = Σ wtk · occ · ε"]
        B8["<b>mkrho</b>：<br/>ρ(r) = Σ wtk · occ · |ψ(r)|²<br/>(零占据带跳过 FFT)"]
        B1 --> B2 --> B3
        B3 -- "true (固定)" --> B4
        B3 -- "false (可变)" --> B5 --> B6
        B4 --> B7
        B6 --> B7
        B7 --> B8
    end

    B8 --> C1

    subgraph PhaseEtot["阶段三：能量汇总 — etotfor"]
        direction TB
        C1{"3 ≤ occopt ≤ 8?"}
        C2["e_entropy = −tsmear × S<br/><i>Mermin 自由能 F = E − TS</i>"]
        C3["e_entropy = 0"]
        C4["etotal = e_band − e_dc + e_entropy + E_ewald + ..."]
        C1 -- "是" --> C2
        C1 -- "否 (occopt ≤ 2 或 = 9)" --> C3
        C2 --> C4
        C3 --> C4
    end

    C4 --> D1{"收敛？<br/>res2 &lt; tolvrs ?<br/>|Δetotal| &lt; toldfe ?"}
    D1 -- "否" --> Mix["Pulay/Anderson 混合<br/>(m_ab7_mixing)"]
    Mix --> B1
    D1 -- "是" --> End([输出 WFK, EIG,<br/>etotal, fermie, forces])

    classDef fixedOcc fill:#e3f2fd,stroke:#1976d2,stroke-width:1.5px,color:#0d47a1
    classDef varOcc   fill:#fce4ec,stroke:#c2185b,stroke-width:1.5px,color:#880e4f
    classDef phase    fill:#f3e5f5,stroke:#6a1b9a,stroke-width:2px
    classDef decision fill:#fff9c4,stroke:#f9a825,stroke-width:1.5px,color:#e65100
    class A5,B4,C3 fixedOcc
    class A6,B5,B6,C2 varOcc
    class A4,B3,C1,D1 decision

图例：

• 蓝色节点代表"固定占据"路径（occopt = 0, 1, 2）
• 粉色节点代表"可变占据"路径（occopt = 3..9，需在 SCF 中调 newocc）
• 黄色菱形是 occopt 控制的关键分支判断
• occopt = 9 落在阶段二的可变路径，但不进入 Mermin 自由能（C1 判断为"否"）

典型例子

例 1：Silicon 半导体（推荐 occopt=1）

nband 8           # 8 价电子，半导体只需 4 占据 + 几条空带
occopt 1
# 不需要 tsmear

dtset_initocc_chkneu 自动生成 occ = [2, 2, 2, 2, 0, 0, 0, 0]。SCF 不调 newocc，占据数永不改变。

例 2：Aluminum 金属（推荐 occopt=7）

nband 10          # 3 价电子 + 缓冲带（fband=0.5 默认）
occopt 7
tsmear 0.01       # Ha = 3160 K
ngkpt 16 16 16

初始 occ_orig = [2, 0.5, 0, 0, 0, 0, 0, 0, 0, 0]（来自 initocc_chkneu 第 4.1 节）。每次 SCF 步通过 newocc 重算：

istep 1 : occ = [2.0, 0.50, 0.001, 1e-8, 0, ...]
istep 2 : occ = [2.0, 0.48, 0.020, 1e-5, 0, ...]
...
istep 8 : occ = [1.99, 0.41, 0.45, 0.15, 0.001, ...]  # 真实 metallic 分布

fermie 在迭代中逐步收敛到正确值。

例 3：Fe 铁磁（occopt=7, nsppol=2, spinmagntarget=2.2）

nsppol 2
nspinor 1
occopt 7
tsmear 0.01
spinmagntarget 2.2     # μ_B per Fe atom

initocc_chkneu 第 4.1 节"情况 B"路径：

nelect_up   = (8 + 2.2)/2 = 5.1    ! 5.1 个 spin-up 电子
nelect_down = (8 - 2.2)/2 = 2.9    ! 2.9 个 spin-down 电子

注意 maxocc = 2/(2*1) = 1 了（每带每自旋通道最多 1 个）：

occ_spin_up   = [1, 1, 1, 1, 1, 0.1, 0, ...]
occ_spin_down = [1, 1, 0.9, 0, ...]

SCF 后 newocc 在分自旋通道上独立做二分（m_occ.F90:836-913 的固定磁矩分支），各自旋通道有不同的 Fermi 能级。

例 4：ΔSCF 激发态（occopt=9）

occopt 9
ivalence 4         # 价带最高指标
nqfd 1.0           # 强制 1 个电子激发
nband 10
tsmear 0.005

initocc_chkneu 第 4.2 节路径，初始 occ = [2, 2, 2, 1, 0, 1, 0, 0, 0, 0]（价带顶 1 空穴，导带底 1 电子）。SCF 中 newocc 同时维持 2 个 Fermi 能级 $\mu_e, \mu_h$（m_occ.F90:583-606, 718-737）。

熵不进入 Mermin 自由能（occopt=9 排除），因为它是约束态而非有限温。

例 5：STM 模拟（occopt=7, prtstm=1, stmbias=0.001）

m_occ.F90:786-829：在 Fermi 能级附近 [μ - stmbias, μ] 区间内的占据数被单独保留：

call getnel(... fermie_biased = fermie - stmbias ...)  ! 算偏置后的占据
occ(:) = occ_at_fermie(:) - occ_at_biased(:)            ! 只留差值

得到的 occ 仅在 STM 偏压窗口内非零，对应 STM 像的"差分密度"。

常见陷阱与诊断

"occopt=1 整数填充不通过"

症状：

ERROR initialization of occupation numbers ... charge balance ...

原因：nelect / maxocc 不是整数（半填的 band）。

对策：

• 改用 occopt 7 + tsmear 0.01（金属型）
• 或检查 cellcharge 是否设错

"nsppol=2 + occopt=1 但没给 spinmagntarget"

症状：

ERROR This is a calculation with spin-up and spin-down wavefunctions ...
must be specified, while the default value is observed.

对策：

• 加 spinmagntarget X.X（先做小型测试估算磁矩）
• 或改 occopt 7（让代码自己决定磁矩）

"金属用 occopt=1 收敛极慢/振荡"

症状：SCF 不收敛，能量上下抖动。

原因：Fermi 面附近的能带跨越导致占据数突变。

对策：必须用 occopt 7（或 3、4、6），加 tsmear 0.005..0.02。

"occopt=2 的 wtk 没归一化"

症状：nelect_occ ≠ nelect 警告。

原因：occopt=2 是唯一一个不自动归一化 wtk 的选项。

对策：

• 显式让 $\sum w_k = 1$
• 或改 occopt 1 让代码自动归一化

"Cold smearing (4, 5, 6) 能量跳变"

症状：几何优化中相邻几何能量出现台阶。

原因：cold smearing 的占据函数非单调，二分法可能落到不同 Fermi 解上。

对策：

• 改用单调的 occopt 7（Gaussian）
• 或 occopt 5（cold smearing 单调版）

"几何弛豫中 nband 不够"

症状：

ERROR There are not enough bands to get charge balance right

原因：体积膨胀后能带数下降，原本占据带超出 nband。

对策：增大 nband 留出余量；或用 fband 0.7 而非 0.5。

关键代码地图

附录：一次完整 SCF 中 occ 的逐步追踪

设输入：FCC Al（3 价电子），occopt=7，tsmear=0.01，ngkpt 8 8 8 → 10 IBZ k 点，nband=6。

Step 0: 初始化

1. invars1：fband=0.5，mband_upper = 1 * (ceiling(3) + ceiling(0.5)) = 1 * (3 + 1) = 4。用户给了 nband=6，覆盖之。
2. dtset_initocc_chkneu：
   • zval = 3, cellcharge = 0 → nelect = 3
   • maxocc = 2 / (1*1) = 2
   • nocc = int(3/2 - 1e-8) + 1 = 2
   • occlast = 3 - 2*1 = 1 → 第 2 带半填
   • occ_orig = [2, 1, 0, 0, 0, 0]（10 个 k 点都一样）
3. 电荷检查：Σ wtk × occ = 1 × (2 + 1 + 0 + 0 + 0 + 0) = 3 = nelect ✓

Step 1: 第一次 SCF

1. vtorho 设 fixed_occ = (7 < 3) = .false.
2. 求解 H ψ = ε ψ 得 eigen = [k=1: -0.42, -0.08, +0.15, ...; k=2: ...; ...]（10 k × 6 bands = 60 eigenvalues）
3. 调 newocc：
   • 边界：fermie_lo = min(eigen) - 0.06 = -0.48, fermie_hi = max(eigen) + 0.06
   • 二分到 fermie ≈ -0.05 Ha
   • 用 Gaussian smearing 表查得每个 (ε, k) 处的 occ
   • 结果约：occ_k1 = [2.0, 1.6, 0.05, 1e-7, ...]（k=1 处）；其他 k 略不同
4. mkrho：ρ(r) = Σ_{k,n} wtk(k) × occ(k,n) × |ψ_kn(r)|²
5. etotfor：E_tot = E_band - E_dc - tsmear × S = ... - 0.01 × 0.02 = ...

Step 2-7: 后续 SCF

每步 occ 都微调（因为 eigen 在变）。典型情况下：

istep    fermie       Σ entropy×wtk     occ at fermi level
  1     -0.0521      0.0184            ~0.5
  2     -0.0498      0.0156            ~0.5
  ...
  7     -0.0501      0.0152            ~0.49

occ 数组在每次 SCF 步都会重新计算，但因为 eigen 变化越来越小，occ 也越来越稳定。

Step 8: 收敛退出

res2 < tolvrs (1e-12) 连续 2 次 → quit
final fermie = -0.0501 Ha
final entropy = 0.0152
final etotal = ... - 0.0001 (= -tsmear × S)

报告生成于 2026-05-19，基于 Abinit 10.4.7。

References

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

Marzari, Nicola, and David Vanderbilt. "Maximally Localized Generalized Wannier Functions for Composite Energy Bands." Physical Review B, vol. 56, no. 20, 1997, pp. 12847–12865. link

Mermin, N. David. "Thermal Properties of the Inhomogeneous Electron Gas." Physical Review, vol. 137, no. 5A, 1965, pp. A1441–A1443. link

Methfessel, M., and A. T. Paxton. "High-Precision Sampling for Brillouin-Zone Integration in Metals." Physical Review B, vol. 40, no. 6, 1989, pp. 3616–3621. link

Paillard, Charles, et al. "Noncollinear Magnetism from First Principles." Physical Review B, vol. 100, no. 2, 2019, p. 024426. link

The Abinit Code. Version 10.4.7, 2024. link