Abinit K点（布里渊区采样）：参数 + 公式 + 代码实现详解

基于 Abinit 10.4.7 源码分析（核心：src/56_recipspace/m_kpts.F90、shared/common/src/29_kpoints/m_symkpt.F90、src/57_iovars/m_inkpts.F90）

引言：K点的物理含义

固体物理中，由 Bloch 定理，单粒子态 $\psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n\mathbf{k}}(\mathbf{r})$ 用波矢 $\mathbf{k}\in\text{BZ}$ 标记。任何物理量（密度、能量、力）都要在第一布里渊区上做积分：

$$n(\mathbf{r}) = \sum_n \frac{V}{(2\pi)^3} \int_{\text{BZ}} f_{n\mathbf{k}} |\psi_{n\mathbf{k}}(\mathbf{r})|^2 d^3\mathbf{k}$$

DFT 求解只能在离散的 k 点上做。Abinit 把这个积分变成加权和：

$$\int_{\text{BZ}} F(\mathbf{k}),d^3\mathbf{k} ;\approx; \frac{(2\pi)^3}{V}\sum_{i=1}^{N_k} w_i,F(\mathbf{k}_i),\qquad \sum_i w_i = 1$$

K点采样涉及三件事：

1. 生成均匀网格（Monkhorst–Pack 方案）
2. 用对称操作折叠到不可约布里渊区 IBZ（symkpt）
3. 分配权重 $w_i$（每个 IBZ 点代表多少个 BZ 等效点）

特殊情形：能带结构计算需要沿高对称线采样（kptopt<0 路径模式）。

输入参数总览

下表按使用频率列出所有 k 点相关变量（定义位置：abimkdocs/variables_abinit.py）。

主要场景对照：

• 半导体/绝缘体 SCF：kptopt 1, ngkpt 4 4 4, nshiftk 1, shiftk 0.5 0.5 0.5
• 金属 SCF：用更密的 ngkpt（如 12 12 12），加 tsmear（见 smearing 报告）
• 能带结构：先 SCF（prtden 1），再 iscf -2, kptopt -N, kptbounds ..., ndivsm 20
• DFPT 响应函数：kptopt 2（仅时间反演）或 kptopt 3（无对称）

Monkhorst–Pack 公式与代码实现

数学定义

对每个方向 $\alpha\in{1,2,3}$，取 $N_\alpha = \tt{ngkpt}(\alpha)$ 个点：

$$\mathbf{k}{ijk} = \sum{\alpha=1}^3 \frac{i_\alpha - 1 + s_\alpha}{N_\alpha},\mathbf{b}\alpha,\qquad i\alpha \in {1,\ldots,N_\alpha}$$

其中 $\mathbf{b}\alpha$ 是倒格矢，$s\alpha = \tt{shiftk}(\alpha)$ 是平移。

更一般的形式用 kptrlatt：k 点是晶格倒格矢

$$\mathbf{B}{\text{kpt}} = (B{\text{rec}})\cdot(\tt{kptrlatt})^{-T}$$

的子格。一个简单的对应关系：当 nshiftk=1 时，kptrlatt 是对角矩阵 diag(ngkpt)。

代码：smpbz（Sampling of Brillouin Zone）

位置：src/56_recipspace/m_kpts.F90:2088-2615。

主要分支按 brav 区分：

• brav=1：通用（任意 kptrlatt，包含简单晶格）
• brav=2：面心立方专用
• brav=3：体心立方专用
• brav=4：六方专用

通用分支的核心循环（m_kpts.F90:2274-2354）：

nn = 1
do kk = boundmin(3), boundmax(3)
   coord(3) = kk
   do jj = boundmin(2), boundmax(2)
      coord(2) = jj
      do ii = boundmin(1), boundmax(1)
         coord(1) = ii
         do ikshft = 1, nshiftk
            ! 试探 k 点在 k 晶格上的整数坐标
            k1(1) = ii + shiftk(1, ikshft)
            k1(2) = jj + shiftk(2, ikshft)
            k1(3) = kk + shiftk(3, ikshft)
            ! 变换到倒格矢约化坐标
            k2(:) = k1(1)*klatt(:,1) + k1(2)*klatt(:,2) + k1(3)*klatt(:,3)
            ! 只保留 [0,1) 内的点
            if (k2(1) < -tol10 .or. k2(1) > one-tol10) cycle
            if (k2(2) < -tol10 .or. k2(2) > one-tol10) cycle
            if (k2(3) < -tol10 .or. k2(3) > one-tol10) cycle
            ! 折叠到 ]-1/2, 1/2]
            call wrap2_pmhalf(k2(1), k1(1), shift)
            call wrap2_pmhalf(k2(2), k1(2), shift)
            call wrap2_pmhalf(k2(3), k1(3), shift)
            spkpt(:, nn) = k1(:)
            nn = nn + 1
         end do
      end do
   end do
end do
nkpt = nn - 1

整个 BZ 上的 k 点总数应满足：

$$N_{\text{BZ}} = \det(\tt{kptrlatt}) \times \tt{nshiftk}$$

这是通过 nkptlatt = det(kptrlatt)（m_kpts.F90:2179-2184，按矩阵展开公式）然后乘 nshiftk 来核对的。

自动选优：testkgrid

如果输入既没给 ngkpt 也没给 kptrlatt，Abinit 会枚举大量候选 kptrlatt + shiftk 组合，按以下准则选优：

• 首要标准：所谓 k 点格的实空间长度 kptrlen（k 格的倒格 = 一个超晶格，其最短实空间向量长度即 kptrlen）。比 kptrlen（默认 30 Bohr）大的网格才被接受。
• 次要标准：使用 kptopt=1 全对称约化后 IBZ 点数最少的组合获胜。

代码位置：m_kpts.F90:2654（testkgrid 子程序）。这一步耗时较多（枚举可达上万组合），所以建议 prtkpt=1 单独跑一次得到推荐网格，记下来再用。

对称约化：BZ → IBZ

数学原理

固体的点群 ${S_i}$（含可能的时间反演 $T$）对 k 空间有作用：

$$\mathbf{k}' = S_i,\mathbf{k} \pmod{\mathbf{G}}$$

如果 $\mathbf{k}_a$ 与 $\mathbf{k}_b$ 由某个对称操作连接（可能差一个倒格矢 $\mathbf{G}$），则它们是等价的，只需计算其中一个，权重相加。

经过 BZ→IBZ 约化后：

$$\sum_{\mathbf{k}\in\text{BZ}} F(\mathbf{k}) = \sum_{\mathbf{k}\in\text{IBZ}} w_\mathbf{k},F(\mathbf{k}),\qquad w_\mathbf{k} = \frac{|\text{star}(\mathbf{k})|}{|G_0|}$$

其中 $|\text{star}(\mathbf{k})|$ 是 k 点的 star 大小，$|G_0|$ 是阶（用于归一化使权重总和=1）。

代码：symkpt

位置：shared/common/src/29_kpoints/m_symkpt.F90:93-404。

kptopt 的语义控制对称使用范围（见 m_kpts.F90:1401-1424）：

代码片段：

if (kptopt == 1 .or. kptopt == 4) then
   nsym_used = 0
   do isym = 1, nsym
      if (symafm(isym) == 1 .or. kptopt == 4) nsym_used = nsym_used + 1
   end do
   ABI_MALLOC(symrec, (3,3,nsym_used))
   ! 把实空间对称矩阵 symrel 取逆转置得到倒空间矩阵 symrec
   do isym = 1, nsym
      if (symafm(isym) == 1 .or. kptopt == 4) then
         call mati3inv(symrel(:,:,isym), symrec(:,:,nsym_used))
      end if
   end do
else if (kptopt == 2) then
   nsym_used = 1
   ! 只放一个单位矩阵；时间反演由 timrev=1 单独开启
   symrec(1:3, 1:3, 1) = identity
end if

核心约化算法（m_symkpt.F90:269-340）：

! 按 k 向量长度排序，方便快速跳出
do ikpt = 1, nkbz - 1
   ind_ikpt = list(ikpt)               ! 排序后的索引
   if (wtk_folded(ind_ikpt) < tol16) cycle ! 已被合并
   do ikpt2 = ikpt + 1, nkbz
      if (length2(ikpt2) - length2(ikpt) > tol8) exit  ! 长度不同，必不等价
      ind_ikpt2 = list(ikpt2)
      if (wtk_folded(ind_ikpt2) < tol16) cycle
      ! 尝试用每个对称操作把 ikpt 映到 ikpt2
      do isym = 1, nsym
         do itim = 1, (1 - 2*timrev), -2     ! itim = +1 然后 -1（若 timrev=1）
            if (isym /= identi .or. itim /= 1) then
               ksym(:) = itim * matmul(symrec(:,:,isym), kbz(:,ind_ikpt))
               ! 检查 ksym - kbz(:,ind_ikpt2) 是否为整向量（差一个倒格矢）
               if (sum(|reduce|) < tol8) then
                  ! 找到对称等价，把 ikpt2 的权重并到 ikpt
                  wtk_folded(ind_ikpt)  = wtk_folded(ind_ikpt) + wtk_folded(ind_ikpt2)
                  wtk_folded(ind_ikpt2) = zero
                  ! 记录映射：S*k1 + G = k2
                  bz2ibz_smap(1, ind_ikpt2) = ind_ikpt
                  bz2ibz_smap(2, ind_ikpt2) = isym
                  bz2ibz_smap(3:5, ind_ikpt2) = nint(-ksym(:) + kbz(:,ind_ikpt2) + tol12)
                  bz2ibz_smap(6, ind_ikpt2) = (itim == -1) ? 1 : 0
                  exit
               end if
            end if
         end do
      end do
   end do
end do

! 取出权重非零的点作为 IBZ
nkibz = 0
do ikpt = 1, nkbz
   if (wtk_folded(ikpt) > tol8) then
      nkibz = nkibz + 1
      ibz2bz(nkibz) = ikpt
   end if
end do

输出：

• nkibz：IBZ 内 k 点数（最终的 nkpt）
• ibz2bz(1:nkibz)：IBZ 索引 → BZ 索引
• wtk_folded(1:nkbz)：聚合后的权重（IBZ 点上非零，其余为零）
• bz2ibz_smap(6, 1:nkbz)：BZ 上每点对应的 IBZ 点 + 用了哪个对称操作 + 偏移倒格矢 + 时间反演标志

权重的最终归一化

getkgrid 完成后，所有权重之和应为 1：

$$\sum_{i\in\text{IBZ}} w_i = 1$$

代码会自动除以 $N_{\text{BZ}} = \det(\tt{kptrlatt}) \cdot \tt{nshiftk}$ 来归一化。

能带路径（kptopt < 0）

输入解析

kptopt = -N 表示路径有 N 段、N+1 个端点。代码位置：src/57_iovars/m_inkpts.F90:285-373。

ABI_MALLOC(kptbounds, (3, nsegment+1))    ! 端点
ABI_MALLOC(ndivk, (nsegment))             ! 每段分点数

! 路径生成：
jkpt = 1
kpt(:, 1) = kptbounds(:, 1)
do ii = 1, nsegment
   dkpt = ndivk(ii)
   do ikpt = 1, dkpt
      fraction = dble(ikpt) / dble(dkpt)
      kpt(:, ikpt+jkpt) = fraction * kptbounds(:, ii+1) &
                       + (one-fraction) * kptbounds(:, ii)
   end do
   jkpt = jkpt + dkpt
end do

ndivsm 自动均匀分段

如果用户给的是 ndivsm（推荐），由 mknormpath（m_kpts.F90:3318-3398）根据每段在倒空间中的物理长度自动决定每段的 ndivk：

do ii = 1, nbounds-1
   dd(:) = bounds(:,ii+1) - bounds(:,ii)
   ! 倒空间度规：lng = sqrt(d^T * gmet * d)
   lng(ii) = sqrt( dd(1)*gmet(1,1)*dd(1) + dd(2)*gmet(2,2)*dd(2) + dd(3)*gmet(3,3)*dd(3) &
                 + 2*(dd(1)*gmet(1,2)*dd(2) + dd(1)*gmet(1,3)*dd(3) + dd(2)*gmet(2,3)*dd(3)) )
end do
fct = minval(lng) / ndiv_small
ndiv(:) = nint(lng(:) / fct)
npt_tot = sum(ndiv) + 1

这样最短段恰好有 ndivsm 个分点，其他段按比例放大，可视化时每两个 k 点在图上等距。

kptopt<0 路径模式下权重不需要（NSCF 计算，不积分）；getkgrid 跳过 symkpt。

getkgrid：高层封装

位置：m_kpts.F90:1257-1842。这是 SCF 设置阶段调用的真正入口（由 src/57_iovars/m_inkpts.F90:468 调起）。

它的工作流：

1. 预处理 nshiftk：尝试用素因子分解减少 shiftk 数量（m_kpts.F90:1432-1690，约 250 行的优化算法）。如果两个 shift 之间差一个 k 格矢量，可以合并成更大的 kptrlatt + 更少的 shift。
2. 调 smpbz 生成完整 BZ 上的所有 k 点（spkpt(3, nkptlatt*nshiftk)）和均匀权重 $1/N$。
3. kptopt=1,2,4：调 symkpt 做对称约化，得到 IBZ k 点和权重。
4. kptopt=3：直接保留 BZ 全部点。
5. 检查对称性破缺（chksymbreak=1 时）：所有对称操作下网格应自映射，否则报错。
6. 输出：kpt(3, nkpt) 是 IBZ 点，wtk(nkpt) 是归一化权重。

istwfk：时间反演压缩存储

当 k 点的所有分量都是 0 或 1/2 时（"半整数点"），波函数满足额外的对称性：

$$u_{\mathbf{G}0/2}(\mathbf{G}) = u{\mathbf{G}_0/2}(-\mathbf{G} - \mathbf{G}_0)^*$$

只需存一半的 G 向量，内存和 FFT 时间减半。

代码：src/44_abitools/m_cgtools.F90:919-950：

integer pure function set_istwfk(kpoint) result(istwfk)
   real(dp), intent(in) :: kpoint(3)
   integer :: bit0, ii, bit(3)
   bit0 = 1
   do ii = 1, 3
      if (DABS(kpoint(ii)) < tol10) then
         bit(ii) = 0
      else if (DABS(kpoint(ii) - half) < tol10) then
         bit(ii) = 1
      else
         bit0 = 0
      end if
   end do
   if (bit0 == 0) then
      istwfk = 1   ! 普通 k 点，不能压缩
   else
      istwfk = 2 + bit(1) + 4*bit(2) + 2*bit(3)
   end if
end function

返回值含义：

• istwfk = 1：通用 k 点，不压缩
• istwfk = 2：$\Gamma$ 点 (0, 0, 0)
• istwfk = 3..9：八个 BZ 边角点（如 (1/2, 0, 0)、(1/2, 1/2, 1/2) 等）

用户通常不设 istwfk，让 Abinit 在 m_inkpts.F90:508-520 自动判断。响应函数（RF）、自旋极化（nspinor=2）、SOC 等情形强制为 1。

K点在 SCF 中的并行化

进程分布：distrb2

位置：src/51_manage_mpi/m_mpinfo.F90:2299-2574。

把每个 (k 点, 自旋, 能带) 三元组分配到一个 MPI 进程：

mpi_enreg%proc_distrb(ikpt, iband, isppol) = 该任务所属进程号

策略：

• paralbd=0：按 k 点平均分配，要求 nkpt*nsppol 能被 nproc_spkpt 整除（否则警告"不高效"）。
• paralbd=1：再把每个 k 点的能带切到多个进程上。

之后 mpi_enreg%my_kpttab(ikpt) 给出本地 k 点索引（仅属于本进程的那些）。

SCF k 点循环（vtorho）

位置：src/79_seqpar_mpi/m_vtorho.F90:826-1300（标注为"BIG FAT k POINT LOOP"）。骨架：

do isppol = 1, nsppol
   ! 装载该自旋通道的有效势 vlocal
   call gs_hamk%load_spin(isppol, vlocal=vlocal, ...)

   ikpt_loc = 0
   do while (ikpt_loc < nkpt1)
      ikpt_loc = ikpt_loc + 1
      ikpt = ikpt_loc
      my_ikpt = mpi_enreg%my_kpttab(ikpt)

      ! 跳过非本进程的 k 点
      if (proc_distrb_cycle(mpi_enreg%proc_distrb, ikpt, 1, nband_k, isppol, me_distrb)) cycle

      nband_k = dtset%nband(ikpt + (isppol-1)*dtset%nkpt)
      istwf_k = dtset%istwfk(ikpt)
      npw_k   = npwarr(ikpt)

      ! 调对角化器（CG / RMM-DIIS / LOBPCG / ChebFi），算出 nband_k 个本征值
      call vtowfk(...)

      ! 累加每个量到全局：能量、密度
      do iband = 1, nband_k
         energies%e_kinetic     = energies%e_kinetic     + wtk(ikpt)*occ_k(iband)*ek_k(iband)
         energies%e_eigenvalues = energies%e_eigenvalues + wtk(ikpt)*occ_k(iband)*eig_k(iband)
         energies%e_nlpsp_vfock = energies%e_nlpsp_vfock + wtk(ikpt)*occ_k(iband)*enlx_k(iband)
      end do
      ! rho(r) += wtk(ikpt) * occ_k * |psi_k(r)|^2  (在 fft 网格上累加)
   end do

   ! MPI 规约：把各进程的 rho、E 加在一起
   call xmpi_sum(rhoaug, spaceComm_distrb, ierr)
end do

wtk(ikpt) 在这里就是 BZ 积分的离散权重。所有物理量都通过这个权重做 BZ 平均。

占据数计算（衔接 smearing）

每次 SCF 迭代后：

1. 收集所有 $\epsilon_{n\mathbf{k}}$；
2. 调 newocc（见 smearing 报告）二分法求 Fermi 能级和占据数 $f_{n\mathbf{k}}$；
3. 用 $w_\mathbf{k} \cdot f_{n\mathbf{k}}$ 加权累积下一轮 $\rho(\mathbf{r})$。

Abinit 中 k 点的完整生命周期

下面按时间顺序，描述一次 SCF 计算从输入到输出，k 点的每个阶段。

阶段 1：输入解析（m_inkpts.F90）

1. 读 kptopt、ngkpt、shiftk、kptrlatt 等
2. 若 kptopt < 0：直接构造路径 kpt(:,:)（m_inkpts.F90:285-380）
3. 若 kptopt > 0 且未给 ngkpt/kptrlatt：调 testkgrid 自动搜索最优网格（m_inkpts.F90:456）
4. 调 getkgrid 生成 IBZ k 点 + 权重（m_inkpts.F90:468）
5. 对每个 k 点调 set_istwfk 决定存储模式（m_inkpts.F90:511）

阶段 2：MPI 初始化（m_mpinfo.F90）

1. 调 distrb2 把 k 点分配到 MPI 进程：mpi_enreg%proc_distrb
2. 建立本地映射 mpi_enreg%my_kpttab(ikpt)

阶段 3：每个 k 点的内存预分配

• 在 G 空间球内填入 G 向量：m_kg.F90:kpgio，每个 k 点的 npw_k 和 kg_k(3, npw_k)
• 计算非局域投影：nonlocal/、fock/ 等

阶段 4：SCF 主循环（m_dft_energy/m_vtorho）

对每次迭代 istep = 1..nstep：

do isppol = 1, nsppol               # 自旋外循环
   load V_xc, V_H, V_ext            # 当前势函数
   do ikpt = 1, nkpt                # k 点循环
      if (该进程不负责 ikpt) cycle
      H_k = T_k + V_loc + V_nl(k)   # 构造 k 处 Hamiltonian
      solve H_k |psi_nk> = e_nk |psi_nk>      # 对角化（npw_k * nband_k 维）
      do iband = 1, nband_k
         rho(r) += wtk(ikpt) * occ_k(iband) * |psi_nk(r)|^2
         E_band += wtk(ikpt) * occ_k(iband) * e_nk
      end do
   end do
end do
xmpi_sum(rho, E_band, ...)           # 跨进程规约

# 算 V_xc[rho]、E_total = E_band - E_dc + ... - TS
# 调 newocc 更新 occ 和 fermie
if (rho 已收敛) exit

阶段 5：输出（m_iowf、m_clnup1）

• 写 WFK 文件：每个 k 点的 (kpt, wtk, npw_k, kg_k, eigen_k, occ_k, cg_k)
• 写 EIG 文件：本征值数据
• 写 DOS（若 prtdos=1）：在所有 IBZ 点上用 smearing 函数做加权求和

典型例子

FCC Si（半导体），8×8×8 网格

输入：

kptopt 1
ngkpt 8 8 8
nshiftk 4
shiftk 0.5 0.5 0.5
       0.5 0.0 0.0
       0.0 0.5 0.0
       0.0 0.0 0.5

执行：

• BZ 全部点数：det(kptrlatt) * nshiftk = 8³ * 4 = 2048
• Si 的 Fd-3m 对称群有 48 个操作 + 时间反演 = 96 个
• symkpt 把 2048 → 通常约 60 个 IBZ 点（具体数依赖网格规整性）
• 每个 IBZ 点权重 $w_i \in {1/2048, 2/2048, ...}$，总和=1

BCC Fe（金属磁性），16×16×16

kptopt 1
ngkpt 16 16 16
nshiftk 1
shiftk 0.5 0.5 0.5
nsppol 2
occopt 7
tsmear 0.001

执行：

• BZ 点数：4096
• O_h 群 + 时间反演 = 96
• IBZ ≈ 126 点（具体由 Abinit 计算）
• 每个 IBZ 点要算 2 套自旋（共 252 个独立对角化任务）
• 用 MPI 16 进程：每进程平均 16 个 k 点，建议设置 paral_kgb 1 让 nproc=nproc_kpt × nproc_band × nproc_fft

GaAs 能带结构

ndtset 2
# dataset 1: SCF
kptopt1 1
ngkpt1 4 4 4
shiftk1 0.5 0.5 0.5
prtden1 1
# dataset 2: NSCF along L-Γ-X-W-K-Γ
iscf2 -2
getden2 -1
kptopt2 -5
kptbounds2 0.5  0.5  0.5  # L
           0.0  0.0  0.0  # Gamma
           0.0  0.5  0.5  # X
           0.25 0.5  0.75 # W
           0.375 0.375 0.75 # K
           0.0  0.0  0.0  # Gamma
ndivsm2 20
tolwfr2 1d-10

mknormpath 算出路径总长，按最短段 20 分点比例分配，生成约 100 个 k 点的等距路径。

关键代码地图

参数选择建议

• 绝缘体/半导体：ngkpt 4 4 4～8 8 8，nshiftk 1，shiftk 0.5 0.5 0.5。
• 金属：从 12 12 12 起步，根据收敛性增加到 24 24 24 甚至更密。
• 小元胞超导/磁性体：16 16 16～32 32 32 不罕见，配 smearing 收敛。
• 大元胞（>50 原子）：从 2 2 2 或 Gamma-only 开始，必要时加密。
• 二维材料：ngkpt N N 1（真空方向 1 点），用 slabwsrad 切断。
• 能带结构 ndivsm：粗略图 10、出版图 20–40。
• shiftk 的选择：
  • 默认 0.5 0.5 0.5 在很多体系收敛最快，但破缺四方/六方对称，会被 Abinit 拒绝。
  • 兼容性最好但收敛慢：0 0 0（含 Γ）。
  • FCC 优选：4 个 shift (½ ½ ½) (½ 0 0) (0 ½ 0) (0 0 ½)。
  • 六方推荐：0 0 ½。
• chksymbreak：建议保留默认 1，否则错误的网格会静默给出错的结果。

附录：一次完整 SCF 中 k 点处理逐步追踪

设输入：FCC Si, ngkpt 4 4 4, nshiftk 4（默认 FCC 4-shift）, kptopt 1。

Step 1: 输入解析

m_inkpts.F90 读到 kptopt=1, ngkpt=4 4 4, 自动构造 kptrlatt = diag(4,4,4)。

Step 2: 进入 getkgrid

• nshiftk_reduction：发现 4 个 shifts 可以合并成 kptrlatt = ((-4 4 4), (4 -4 4), (4 4 -4)) + nshiftk=1（FCC k 格的紧致表示）。
• 调 smpbz，扫 BZ 上 $\det(\tt{kptrlatt}) = 128$ 个 k 点。

Step 3: symkpt 约化

• Si 有 48 个空间对称 + 时间反演 = 96
• 排序 128 个 k 点的长度（m_symkpt.F90:178-204）
• 双重循环找等价对：每个对称操作作用后看能否落到其他 k 点上
• 最终 nkibz = 10（FCC Si 的标准结果），权重比例如 (1, 8, 6, 24, 8, 6, 12, 8, 24, 3) / 128

Step 4: 设置 istwfk

• 对每个 IBZ k 点，调 set_istwfk
• Γ 点得到 istwfk = 2，X 点得到 istwfk = 3..4（看具体哪个分量），其余多为 1
• 这影响 G 球内 npw_k 的实际存储数

Step 5: G 球预分配

• 对每个 k，从 0.5*|k+G|² < ecut 选出 npw_k 个 G 向量
• 调 kpgio 建立 kg_k(3, npw_k)

Step 6: MPI 分发（假设 5 个 MPI 进程）

• distrb2：10 IBZ kpt / 5 proc = 每进程 2 个 k 点
• 例：proc 0 拿 (k=1,2)，proc 1 拿 (k=3,4)，…

Step 7: SCF 迭代

每步 SCF：

1. 每进程并行求解自己负责的 2 个 k 点的 Hamiltonian 矩阵： $H_\mathbf{k} \psi_{n\mathbf{k}} = \epsilon_{n\mathbf{k}} \psi_{n\mathbf{k}}$（约 200 plane waves × 20 bands）
2. 局部累积 rho_proc(r) += wtk(k) * occ * |psi|²
3. xmpi_sum：把 5 个进程的 rho 加起来
4. 算 V_new(r) = V_H[rho] + V_xc[rho] + V_ext
5. 调 newocc 二分找 Fermi 能级（只用 IBZ 的 10 个 k 点，因为权重已经包含了 BZ 折叠）
6. 检查 tolvrs / toldfe / ... 是否达标

Step 8: 输出

• 写 WFK：保存 kpt(3, 10), wtk(10), eigen(nband, 10), occ(nband, 10), cg(2, npw_k*nband, 10)
• 写 EIG：纯本征值文本
• 若 prtdos=1：在 10 个 IBZ 点上用展宽函数做加权和得到 DOS

报告生成于 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

Monkhorst, Hendrik J., and James D. Pack. "Special Points for Brillouin-Zone Integrations." Physical Review B, vol. 13, no. 12, 1976, pp. 5188–5192. link

The Abinit Code. Version 10.4.7, 2024. link