GPAW杂化泛函与力的检查

本例用一个很小的H2分子，说明如何把普通PBE计算换成PBE0或HSE06，提取能量和力，再判断自己的计算参数是否够用。PBE0混入25%的完整精确交换；HSE06使用25%的短程精确交换，本例采用程序默认屏蔽参数0.11 Bohr⁻¹。这里跑的是各方法自己的自洽计算。

1 下载和运行#

下载杂化泛函实操包，保存为gpaw-hybrids.zip。从登录节点提交时，解压后使用sbatch job.slurm；如果已进入分配好的CPU资源，可依次执行：

python3 -m zipfile -e gpaw-hybrids.zip .
cd gpaw-hybrid-forces
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

程序使用1 MPI进程、1线程，依次计算16个小条件。完整数据在results/VALIDATION.json，本次计算图片在plots/。再次执行相同输入会读取完整结果；修改模型后，请保留旧目录，并在新目录中计算。

仅想重绘本篇原始数据时，在解压目录中先运行bash setup.sh，再运行bash plot.sh --reference；图片输出到reference-plots/，不会启动自洽计算。

2 结构和关键输入#

H—H距离为0.740 Å，两原子在6×6×6 Å盒中心，三个方向均不周期。它是固定几何的数值示例，没有做键长优化，也没有把这个盒大小作为任意分子的推荐值。

atoms = Atoms('H2',
              positions=[[3, 3, 3 - bond / 2],
                         [3, 3, 3 + bond / 2]],
              cell=[6, 6, 6], pbc=False)
calc = GPAW(mode=PW(cutoff),
            xc=method,                 # 'PBE'、'PBE0'或'HSE06'
            nbands=2, kpts=(1, 1, 1),
            convergence={'energy': 1e-8,
                         'density': 1e-7,
                         'eigenstates': 1e-10},
            maxiter=200, txt='native.txt')
atoms.calc = calc
energy = atoms.get_potential_energy()
forces = atoms.get_forces()

本集群GPAW 26.7.0的新接口已实际执行上述名称写法。其他接口或旧版构建可能需要显式指定PW交换后端，例如xc={'name': 'PBE0', 'backend': 'pw'}；不要把旧版FD杂化接口的能力直接套到本例PW接口上。本例已实际计算解析力。

输入中的能量、密度和本征态收敛标准分别控制不同误差，不能用一个energy值代替全部检查。默认本征态收敛针对占据态：本例第二条空带仅保留为原始诊断数据，不用它报告带隙或激发能。需要定量研究空态时，要另设convergence中的bands='all'，并检查空带数、盒大小和截断。

3 实际能量和力#

下表均为0.740 Å、800 eV时的实际结果。能量使用相同H的PAW参照，适合比较本例各方法；它们不是已经减去两个孤立原子能量的解离能。

两个H的受力方向相反，总力接近零。第二个原子的正向力表示在当前几何上有向更大z移动的趋势；不能仅凭一个点的力，就确定优化后键长。

左图只比较每种方法相对于其自身0.740 Å结果的能量变化，便于观察局部曲线。右图展示实际截断检查：

本例先按20 meV/原子的同一比较标准检查：400→600 eV的变化仍大，继续提高到800 eV后，600→800 eV的变化降到约2.5–2.8 meV/原子。这说明提高截断有实际必要；对自己的体系还应检查力、相对能以及所需物性，不能只看绝对能量变化。

4 用能量差分核对力#

在600 eV下，把键长分别设为0.738、0.740、0.742 Å，使用中心差分：

F₂,z ≈ −[E(0.742) − E(0.738)] / 0.004 Å。

这里两原子以盒中心对称移动，结合总力接近零，键长的能量导数正好对应第二个原子的z向力。

附件也计算了HSE06在同一600 eV条件下的实数与复数波函数表示，两者的能量差小于10⁻⁵ eV，力差小于10⁻⁴ eV/Å。复数表示由PW(600, force_complex_dtype=True)开启；本例只验证相同Gamma点分子问题的两种表示，没有把它当作额外的物理体系。

5 导出数据和重画#

计算完成后，可直接导出能量与力的CSV：

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('hybrid-energies-forces.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['method', 'cutoff_eV', 'bond_A', 'complex',
                'energy_eV', 'force_atom2_z_eV_A'])
    for r in v['conditions']:
        w.writerow([r['method'], r['cutoff_eV'], r['bond_A'],
                    r['complex'], r['energy_eV'],
                    r['forces_eV_A'][1][2]])
for method, s in v['summary'].items():
    print(method, 'force-error=', s['force_error_eV_A'],
          '600-to-800 eV/atom=',
          s['cutoff600to800_change_per_atom_eV'])
PYCODE

重新生成本次结果图和本篇归档图：

bash plot.sh
bash plot.sh --reference
ls plots/*.png reference-plots/*.png

用于新材料时，先改结构、周期性、自旋和k点，再逐项检查截断、真空、力以及目标物性。对晶体而言，Gamma点分子的计算成本与收敛规律不能直接迁移。

6 继续使用HSE03与B3LYP#

同一个H2模型也计算了HSE03和B3LYP。HSE03使用25%的短程精确交换，本例默认屏蔽参数为0.106 Bohr⁻¹；B3LYP使用20%的精确交换。补充包提供八个实际条件：每种方法在600 eV下的三个键长，以及800 eV下的中心键长。

下载HSE03与B3LYP补充包，保存到最初解压目录、与gpaw-hybrid-forces文件夹同级的位置，命名为gpaw-variants.zip。如果当前仍在主包目录，执行：

cd ..
python3 -m zipfile -e gpaw-variants.zip .
cd gpaw-hybrid-variants
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

HSE03使用新版PW接口；B3LYP采用先收敛PBE、保存波函数，再以legacy PW接口启动B3LYP的完整流程。run.py会自动处理这些步骤，在自己的results/中生成GPW，不需要用户手动准备检查点。B3LYP的关键顺序是：

# 已在相同结构、截断与k点下完成PBE
pbe.write('PBE-seed.gpw', mode='all')
hybrid = GPAW('PBE-seed.gpw', legacy_gpaw=True,
              xc={'name': 'B3LYP', 'backend': 'pw'},
              mixer=Mixer(beta=0.05, nmaxold=5, weight=50),
              convergence={'energy': 1e-8, 'density': 1e-7,
                           'eigenstates': 1e-10},
              maxiter=200, txt='B3LYP.txt')
atoms.calc = hybrid
energy = atoms.get_potential_energy()
forces = atoms.get_forces()

先读取旧波函数只是初始猜测，B3LYP仍会完成自己的自洽计算。补充包使用与主包相同的6 Å盒、0.740 Å中心键长及H的PAW数据。

不同泛函的总能不能仅按“越低越好”排序来选择方法；应根据要研究的物性和同一方法内部的能量差判断。图中的能量曲线各自减去本方法中心键长的能量。

补充包结果的终端读取方式：

source env.sh
"$YS_PY" - <<'PYCODE'
import json
from pathlib import Path
v = json.loads(Path('results/VALIDATION.json').read_text())
for r in v['conditions']:
    print(r['method'], r['cutoff_eV'], r['bond_A'],
          r['energy_eV'], r['forces_eV_A'][1][2])
for method, s in v['summary'].items():
    print(method, 'force-error=', s['force_error_eV_A'],
          'cutoff-change/eV-per-atom=',
          s['cutoff600to800_change_per_atom_eV'])
PYCODE
bash plot.sh --reference

程序为GPAW 26.7.0 CPU构建，GCC 13.3、HPC-X/OpenMPI 5、MKL 2023.2、Libxc 7.0.0；使用PAW数据26.8.0。接口与PW/FD差别见GPAW精确交换说明。纯EXX及非自洽接口可接着阅读精确交换与非自洽分析。