GPAW的vdW泛函与BEEF集合操作

本篇用一个很小的H₂结合能例子，练习BEEF-vdW自洽计算、ASE集合提取和其他vdW泛函的固定密度后处理。先分别计算H₂与自旋极化H原子，再按相同的泛函系数序号组合相对能，最后导出表格和图片。

H₂用于演示操作，不代表所有弱相互作用体系。得到一个集合标准差，也不等于已经检验真空、网格、结构或实验精度。

1 下载并运行#

下载BEEF-vdW算例包，保存为gpaw-beef-vdw.zip。登录节点解压后用sbatch job.slurm提交。已经进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-beef-vdw.zip .
cd gpaw-beef-vdw
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

脚本使用1个MPI进程、1个线程。它先完成两次自洽计算，再对已得到的密度计算2000组泛函集合；不需要为每一组集合重新做自洽计算。首次运行后，完整结果保存在results/。相同输入再次执行run.sh会复用已完成结果；修改输入时，请保留原目录，在新目录重新计算。

2 结构、网格和原子自旋#

H₂键长固定为0.75 Å，不做几何优化；非周期盒子为8×8×8.75 Å，目标网格间距0.20 Å。H原子使用相同盒子，置于中心，并设置一个未配对电子。

from ase import Atoms

h2 = Atoms('H2', positions=[[0, 0, 0], [0, 0, 0.75]])
h2.center(vacuum=4.0)
h = Atoms('H', cell=h2.cell, pbc=False, magmoms=[1])
h.center()

孤立H原子需要使用自旋极化态。不能为省事把它当成自旋配对原子，否则参考能发生变化，得到的结合能也会变化。

H₂关键计算参数如下；H原子的相应参数改为spinpol=True，并把fixmagmom设为True：

from gpaw import GPAW, FermiDirac, Mixer

calc = GPAW(
    legacy_gpaw=True, mode='fd', h=0.20,
    xc={'name': 'BEEF-vdW', 'backend': 'libvdwxc'},
    spinpol=False, nbands=2, symmetry='off',
    occupations=FermiDirac(0.01, fixmagmom=False),
    mixer=Mixer(0.25, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-7},
    maxiter=240, txt='native.txt')
h2.calc = calc
energy = h2.get_potential_energy()
calc.write('state.gpw', mode='all')

libvdwxc是这里明确选用的非局域关联后端。不同后端的核参数化可能存在差别，不能假定内置后端与libvdwxc结果逐位相同。

3 先看自洽结合能#

本例以正值表示结合：E_bind = 2 E(H) − E(H₂)。鸢算科技集群实际结果如下：

两项自洽能量使用相同H PAW参考，组合时原子参考能抵消。表中绝对能量属于本计算器的参考约定，研究其他反应时仍需保持各物种的设置一致。

4 用相同系数配对2000组集合#

自洽结束后，对H和H₂分别调用ASE的BEEFEnsemble，使用相同样本数和种子。完整程序还核对两者实际生成的2000×32系数矩阵相同。

from ase.dft.bee import BEEFEnsemble, readbee

ens = BEEFEnsemble(calc, verbose=False)
energies = ens.get_ensemble_energies(2000, seed=0)
ens.write('BEEF')  # 保存为BEEF.bee
energies_readback = readbee('BEEF')

本集群当前ASE接口返回的是集合总能量数组，可核对energies = 自洽能量 + ens.de。变量取名叫de不会把它变成能量差。readbee('BEEF')也返回总能量数组。

因此，保存两次计算的数组后，直接按同一序号组合：

binding_ensemble = 2 * energies_H - energies_H2
mean_binding = binding_ensemble.mean()
sigma_binding = binding_ensemble.std()  # ddof=0

不要在这个数组上再加一次自洽结合能。也不要把H与H₂的2000个样本随意打乱后配对。

这个标准差反映本例相对能在BEEF泛函集合中的分散程度。它不包括有限盒子、网格和几何误差，也不是任意体系的实验置信区间。H₂例子不涉及分子动力学采样，不能把它解释为MD轨迹统计误差。

5 六种vdW泛函的固定密度后处理#

对已经自洽的BEEF密度，还可以计算其他泛函与当前交换关联能的差：

delta = calc.get_xc_difference(
    {'name': 'vdW-DF2', 'backend': 'libvdwxc'})
energy_fixed_density = energy + delta

附件对H和H₂都做同样操作，再组合结合能：

这些结果使用的是BEEF-vdW密度上的非自洽能量。它们可以练习后处理、初步比较设置；需要目标泛函的自洽密度、力或结构优化时，应重新用目标泛函自洽计算。上表没有给出这些目标泛函的结构或受力结果。

6 检查密度和电子数#

两幅图使用各自原生网格上的实际全电子总密度。白色十字标出核位置，两图使用相同色标。H₂积分得到2个电子，H得到1个电子；自旋极化H的上、下自旋电子数之差为1。

读取总密度时，使用不指定spin的接口：

rho = calc.get_all_electron_density(gridrefinement=1)
rho_up = calc.get_all_electron_density(spin=0, gridrefinement=1)
rho_down = calc.get_all_electron_density(spin=1, gridrefinement=1)

对于当前legacy接口的自旋配对计算，spin=0只是半份总密度，spin=1给出另一半；计算器只有一个计算自旋通道，也不意味着把spin=0当成总密度。附件独立核对rho_up + rho_down = rho，避免在电子数积分时少算一半。

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())
e = np.load(r / 'binding-ensemble.npz')['binding_energies_eV']
np.savetxt('binding-ensemble.csv',
           np.column_stack([np.arange(len(e)), e]),
           delimiter=',', header='sample,binding_eV', comments='')
with open('vdw-binding.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['functional_on_BEEF_density', 'binding_eV'])
    w.writerows(v['summary']['nonselfconsistent_binding_energies'].items())
print('SCF binding:', v['summary']['SCF_binding_energy_eV'])
print('Paired mean and std:', e.mean(), e.std())
PYCODE

H/fields.npz和H2/fields.npz还保存完整密度、占据、32项贡献、2000×32系数和集合总能量。.bee文件保存原生ASE集合信息，方便独立读取。

只查看附件原始数据并重画图片时：

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

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

8 环境与下一步#

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

换成吸附能、反应能或层间结合能时，先统一各体系的泛函、PAW数据、网格、边界和自旋设置，再用同一系数序号组合每一组相对能。集合大小与相对能分布也应针对自己的问题检查。

方法参考：GPAW官方vdW与BEEF-vdW说明。本篇数组返回值同时按集群实际安装的ASE实现和.bee往返结果核对。