GPAW的DFT+U与投影占据

DFT+U在选定的局域轨道子空间加入修正。这里用四原子NiO模型比较U=0与U=6 eV，让读者看到输入怎么写、带隙与磁矩怎么提取，以及“投影占据”具体保存在哪。

1 先运行完整例程#

下载NiO DFT+U实操包，保存为gpaw-u.zip。新计算使用CPU资源；从登录节点提交可用sbatch job.slurm。在已分配的CPU资源中：

python3 -m zipfile -e gpaw-u.zip .
cd gpaw-hubbard-u
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

例程采用1 MPI进程、4线程，依次计算两种U设置。results/VALIDATION.json保存带隙、全部采样k点本征值、原子磁矩，以及U=6 eV的原生投影矩阵。图片保存在plots/。

2 结构、磁态与输入#

实际模型可用下面的ASE输入构造：

a = 4.19
b = a / np.sqrt(2)
atoms = Atoms(
    'Ni2O2', cell=(b,b,a), pbc=True,
    positions=[(0,0,0), (b/2,b/2,a/2),
               (0,0,a/2), (b/2,b/2,0)],
    magmoms=[2,-2,0,0]
)
calc = GPAW(
    legacy_gpaw=True, mode=PW(340), kpts=(2,2,2),
    xc='oldPBE', setups={'Ni': '10:d,6.0'},
    occupations=FermiDirac(0.05),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-8},
    maxiter=220, txt='NiO-U6.txt'
)

10选择Ni的10价电子PAW数据，:d,6.0为d子空间指定6 eV的Hubbard参数；本例修正形式是U/2×Tr(N−N²)，对应有效U形式，不能将它与另行指定U、J的其他输入直接混用。U=0对照使用setups={'Ni':'10'}。

legacy_gpaw=True选择兼容计算后端；oldPBE是本例的PBE实现标识，程序仍是26.7.0。两个Ni的初始磁矩相反。结构固定，不进行离子优化；这里的四原子反铁磁排布不代表已经比较了NiO的全部磁相。

3 结果怎样解读#

U=6 eV时，两个Ni的局域磁矩约为+1.785和−1.785 μB，总磁矩接近零。局域磁矩取决于投影区域，不能只把总磁矩为零理解成每个原子都没有磁性。

这个小例用相同几何、340 eV截断与2×2×2网格比较参数影响，不把4.799 eV当作经过充分收敛的材料预测值。自己的体系要检查U的依据、磁态、k点、截断和投影定义。

4 导出带隙与磁矩#

source env.sh
"$YS_PY" - <<'PYCODE'
import json, csv
from pathlib import Path
f = json.loads(Path('results/VALIDATION.json').read_text())
with open('nio-u-results.csv', 'w', newline='') as out:
    w = csv.writer(out)
    w.writerow(['U_eV','gap_eV','energy_eV','Ni0_muB','Ni1_muB'])
    for r in f['results']:
        print('U=',r['U_eV'],'gap=',r['gap_eV'],
              'moments=',r['local_moments_muB'])
        w.writerow([r['U_eV'],r['gap_eV'],r['energy_eV'],
                    *r['local_moments_muB'][:2]])
PYCODE

5 读取原生投影占据#

source env.sh
"$YS_PY" - <<'PYCODE'
import json, numpy as np
from pathlib import Path
f = json.loads(Path('results/VALIDATION.json').read_text())
r = next(r for r in f['results'] if r['U_eV']==6)
for site in r['native_Hubbard_sites']:
    n = np.array(site['projected_occupations'])
    print('Ni atom:', site['atom'], 'shape:', n.shape)
    print('spin traces:', np.trace(n,axis1=1,axis2=2))
    print('native on-site U energy/eV:',site['native_Hubbard_energy_eV'])
    for spin in range(2):
        np.savetxt(f"Ni{site['atom']}-spin{spin+1}-occupation.dat",n[spin])
PYCODE

每个自旋通道是5×5的d投影矩阵。它的迹是该投影定义下的占据，不是化学式中的形式氧化数。配套程序还独立计算Tr(N−N²)，并用密度矩阵有限差分核对原生Hubbard矩阵导数，帮助检查输入确实进入了所选子空间。

只看归档图可运行bash plot.sh --reference。附件给出完整Python输入与数值，不分发PAW数据。环境为GPAW 26.7.0、GCC 13.3、HPC-X/OpenMPI 5、MKL 2023.2、Libxc 7.0.0，PAW数据26.8.0。模型与接口可参照GPAW官方Hubbard小例源码；磁性入门另见自旋与SOC。