GPAW原子自相互作用能量修正

本篇用氦、铍两个小原子，演示“先做普通LDA，再计算非自洽自相互作用修正”的操作。自相互作用误差涉及近似泛函中电子对自身的非物理相互作用；PZ-SIC按轨道扣除相应的Hartree和交换关联贡献。

这里调用GPAW的NSCFSIC。它在已有LDA结果上给出能量修正，没有重新变分优化全部修正轨道。需要修正后的自洽密度、力或结构优化时，请看水分子变分SIC操作，不要把本篇能量后处理当作完整的变分SIC。

1 下载与计算#

下载原子非自洽SIC算例包，保存为gpaw-atomic-sic.zip。解压后，从登录节点用sbatch job.slurm提交；如果已进入分配好的CPU资源，可依次执行：

python3 -m zipfile -e gpaw-atomic-sic.zip .
cd gpaw-atomic-nscfsic
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

包内任务使用1个MPI进程、1线程，依次计算He和Be。results/He/、results/Be/保存本次计算，plots/保存两张图。完整结果在results/VALIDATION.json，重启用的LDA.gpw由自己的计算生成，下载包不包含PAW数据集或GPW文件。

同一输入再次运行会复用完整结果。修改元素、网格或泛函时，保留旧目录，在新目录计算。

2 原子结构与LDA输入#

He和Be各用一个原子，放在8×8×8 Å盒子中心，三个方向不周期；网格目标间距0.15 Å，采用自旋配对LDA。下面是附件使用的关键设置：

from ase import Atoms
from gpaw import GPAW

atoms = Atoms('He')                 # 另一个条件为'Be'
atoms.center(vacuum=4.0)
atoms.calc = GPAW(
    legacy_gpaw=True,
    mode='fd', h=0.15, xc='LDA', spinpol=False,
    convergence={'energy': 1e-6,
                 'density': 1e-5,
                 'eigenstates': 1e-8},
    maxiter=200, txt='LDA.txt')
relative_energy = atoms.get_potential_energy()
atomic_reference = atoms.calc.get_reference_energy()
atoms.calc.write('LDA.gpw', mode='all')

He有2个电子。Be实际有4个电子，本例PAW设置把1s²作为冻结芯层，显式求解的价电子数为2；所以get_number_of_electrons()对Be返回2是正常的，不能据此把它当作Be²⁺。

3 执行原生NSCFSIC#

普通LDA收敛后，核心调用是：

import numpy as np
from gpaw.utilities.sic import NSCFSIC

total_energy = NSCFSIC(
    atoms.calc, rng=np.random.default_rng(42)).calculate()
print(total_energy)                 # eV，包含原子参考能

附件固定随机种子42，使用原生默认的完整Hartree与交换关联修正。为保存程序内部的全精度芯层、价电子分解，附件仅在这次函数调用期间读取Python返回事件中的数值，并立即恢复原来的profile设置；没有修改GPAW源码或类。

GPAW常规能量输出和这里的原子总能采用不同的参考表达。附件明确按下式核对，避免把原子参考能漏掉或重复加上：

E(NSC-SIC总能) = E(LDA相对能) + E(PAW原子参考) + ΔE(SIC)。

4 鸢算科技集群实际结果#

这些原子总能包含各自的PAW原子参考，不能跨元素比较大小来判断谁“更稳定”，也不是已经扣除了其他物种能量的反应能。

GPAW随版本提供的原子示例采用取整后的参考值He −79.4 eV、Be −399.8 eV，比较容差0.1 eV。本例偏差分别为+0.02791、−0.04776 eV。它适合检查这条操作流程，不代表任意体系或物性都具有0.1 eV精度。

5 看电子密度与芯层数据#

左、中图来自实际LDA全电子密度，展示He和Be的空间分布；它们是修正前的LDA密度。右图是本次原生原子生成器产生的Be 1s轨道径向密度，用于芯层Hartree积分。电子密度积分分别得到2.000000和4.000000个电子，包含冻结芯层贡献。

导出能量表并查看电子数：

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('atomic-sic.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['element', 'LDA_total_eV', 'SIC_total_eV',
                'core_correction_eV', 'valence_correction_eV'])
    for r in v['results']:
        w.writerow([r['element'], r['LDA_total_energy_eV'],
                    r['NSC_SIC_total_energy_eV'],
                    r['core_correction_eV'], r['valence_correction_eV']])
        print(r['element'], 'all electrons =',
              r['all_electron_count_integral'])
PYCODE

只想查看附件原始数据、暂时不计算时，在解压目录执行：

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

图保存到reference-plots/，不会启动自洽计算。reference/He/density.npz、reference/Be/density.npz包含完整三维密度、晶胞与原子位置；reference/Be/core-radial.npz包含径向坐标、积分权重和芯层轨道密度。

6 环境与自己的体系#

本例使用GPAW 26.7.0、Python 3.11.16；集群CPU构建采用GCC 13.3、HPC-X Open MPI 5.0.10rc2、MKL 2023.2、Libxc 7.0.0。bash setup.sh加载gpaw/26.7.0-8173m，无需自行安装。

换成其他元素时，应先确认PAW的芯层与价电子划分、电子数、自旋状态和原子参考。再检查网格、真空及轨道局域化。本文没有计算NSC-SIC解析力，也没有使用这种能量修正进行几何优化。

方法来源：GPAW PZ-SIC说明；原子能量比较对应GPAW 26.7.0原生test_nscfsic.py中的He、Be示例及其引用的Perdew–Zunger原子参考。