GPAW金属占据、展宽与四面体积分

金属计算中，费米能级附近的占据会影响自洽收敛和布里渊区积分。本篇用一个很小的Al原胞，实际比较Fermi–Dirac、Methfessel–Paxton、Marzari–Vanderbilt冷展宽，以及普通、改进四面体积分。

先跑完整样例，查看能量和占据，再检查自己的k网格与展宽。这里的4³、6³网格是操作对照，不能当作Al已经收敛的生产设置。

1 下载与运行#

下载金属占据与展宽算例，保存为gpaw-metal-occupations.zip。登录节点解压后用sbatch job.slurm提交；已进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-metal-occupations.zip .
cd gpaw-metal-occupations
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

程序使用1个MPI进程、1个线程。六组占据设置分别在4×4×4和6×6×6 k网格上自洽，共12组计算。结果在results/，图片在plots/。verify.sh用于附件原样复跑；改变自己的体系、网格或方法后，应按新的问题检查物性，不能要求结果仍等于附件参考值。

2 Al原胞与基本输入#

采用fcc Al，常规立方晶格参数4.05 Å，实际计算使用只有1个原子的非正交原胞。PBE、PW 200 eV、8个带，自旋配对，不做结构优化。

from ase.build import bulk
from gpaw import GPAW, PW, Mixer

atoms = bulk('Al', 'fcc', a=4.05)
calc = GPAW(
    legacy_gpaw=True, mode=PW(200), xc='PBE',
    kpts=(6, 6, 6), nbands=8, spinpol=False,
    occupations={'name': 'fermi-dirac', 'width': 0.05},
    mixer=Mixer(0.2, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-7},
    maxiter=200, txt='Al.txt')
atoms.calc = calc
energy = atoms.get_potential_energy()
free_energy = atoms.get_potential_energy(force_consistent=True)
calc.write('state.gpw', mode='all')

附件同时保存默认能量字段和force_consistent=True的自由能字段。有限展宽下，两者定义不同；比较能量差时应统一所用字段，不能随意混用。

3 六组占据设置#

只替换occupations字典，其他参数保持一致：

展宽width的单位为eV。四面体两组没有人为添加width；它们利用规则k网格上的能带插值进行积分，不能把任意一串能带路径点直接当作体积积分网格。

展宽参数也不等于原子分子动力学的温度。MP和冷展宽首先是积分方法；其中的数值修正项不能直接解释成任意热力学温度下的真实电子熵。

4 集群实际能量#

下表使用同一个GPAW默认能量字段，单位eV/原胞；本原胞只有1个Al：

粗网格下，不同方法的能量差可以很大，普通四面体与改进四面体也不能混为一谈。表中两档数据没有证明任何一组已经达到所需精度；自己的研究应继续增加k网格，适当减小展宽，并针对所需能量差或响应量检查。

表中也没有包含孤立Al原子的参考计算，因此不要把这些绝对数值直接称为内聚能。

5 占据数怎样正确求和#

本例Al PAW有3个价电子，全电子密度积分为13。不要把价电子占据和与全电子密度积分混为一谈。

当前接口默认返回的占据已经包含IBZ k点权重和自旋简并因子：

weighted = calc.get_occupation_numbers(kpt=0)
raw = calc.get_occupation_numbers(kpt=0, raw=True)
weights = calc.get_k_point_weights()

对于本例自旋配对计算，每个IBZ点满足：

weighted[q,n] = 2 × weights[q] × raw[q,n]

所以，默认占据数组在全部IBZ点和带上直接求和，就是3个价电子。如果已经拿到默认数组，再乘一次k点权重就会重复加权。读取raw=True时，才需要自己加入权重与本例的简并因子2。

左图展示所用fcc非正交原胞的周期像，只用于结构说明，实际DFT仍是1原子原胞。右图使用6³计算得到的实际能带与原始占据，横坐标是各自能量相对于各自费米能级。

MP、冷展宽等广义占据可能超出0到1的概率区间。图中冷展宽也出现了大于1的值；不要手动裁剪数组，否则会破坏积分和电子数。它们不应被逐项解释成普通FD分布的概率。

6 独立计算占据时的参数#

附件还从实际能带重新计算占据、费米能级和数值熵项，核对它们与自洽结果一致。以下示意使用本例的当前单自旋简并前占据：

from gpaw.occupations import create_occ_calc

occupation = create_occ_calc({'name': 'fermi-dirac', 'width': 0.05})
raw, fermi_levels, entropy = occupation.calculate(
    nelectrons=1.5, eigenvalues=eigenvalues_eV, weights=IBZ_weights)

1.5来自本例3个价电子除以简并因子2。不要把这个数照搬到自旋极化或其他电子数体系。当前安装版本的实际调用以eV能量数组传入；脚本用GPAW返回的eV本征值重建，并独立核对费米能级与熵项单位。

四面体独立重建还需要倒格形状、规则网格大小及BZ→IBZ映射，附件给出完整调用。不能只有一个本征值数组就忽略这些几何信息。

本例保存的约7 eV费米能级依赖内部势能零点，它不是功函数。计算功函数还需要合适的表面模型及真空势参考。

7 导出能量与全部占据#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json, numpy as np
r = Path('results')
v = json.loads((r / 'VALIDATION.json').read_text())
with open('Al-smearing.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'energy_eV', 'free_energy_eV', 'Fermi_eV'])
    for row in v['results']:
        w.writerow([row['name'], row['energy_eV'],
                    row['force_consistent_free_energy_eV'], row['Fermi_eV']])
d = np.load(r / 'FD005-k6/fields.npz')
print('Default weighted sum:', d['weighted_occupations'].sum())
print('Raw + weights + degeneracy:',
      (2 * d['IBZ_weights'][:, None] * d['raw_occupations']).sum())
with open('Al-occupations.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['IBZ_point', 'band', 'weight', 'energy_eV', 'raw', 'weighted'])
    for q, weight in enumerate(d['IBZ_weights']):
        for n, energy in enumerate(d['eigenvalues_eV'][q]):
            w.writerow([q, n, weight, energy,
                        d['raw_occupations'][q, n], d['weighted_occupations'][q, n]])
PYCODE

两种电子数求和均应得到约3。fields.npz保留完整密度、本征值、IBZ点、权重、BZ映射、原始/带权占据以及独立重建结果。

只重画附件参考结果时：

bash setup.sh
bash plot.sh --reference
ls reference-plots

自己的数据和图片在results/、plots/，参考数据和图片在reference/、reference-plots/。附件不包含PAW数据集、GPW检查点或程序二进制。

8 环境与自己的体系#

使用GPAW 26.7.0 legacy PW、Python 3.11.16，CPU构建采用GCC 13.3、HPC-X Open MPI 5.0.10rc2、MKL 2023.2和Libxc 7.0.0。setup.sh加载gpaw/26.7.0-8173m。

处理自己的金属时，先选择合理的占据方法、带数和k网格，使费米面附近的能带充分覆盖，再检查网格与展宽。改变带数、结构或磁性状态后，也应重新检查电子数和所需能量定义。

方法参考：GPAW官方占据与展宽说明。