GPAW变分SIC的FD与PW操作

本篇补充变分PZ-SIC的两种表示：实空间有限差分FD与平面波PW。使用同一个很小的He原子，分别计算两个网格间距和两个截止能，展示输入、结果提取与密度绘图。

与水分子LCAO变分SIC相比，这里调用FDPWETDM，采用LDA和完整的Hartree、交换关联修正，两个缩放系数都是1。不能把两篇不同模型、泛函和缩放系数的绝对能量直接混在一起比较。

1 下载并计算#

下载FD与PW变分SIC算例包，保存为gpaw-sic-fd-pw.zip。从登录节点提交时，解压后使用sbatch job.slurm；已经进入分配好的CPU资源时，依次执行：

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

程序使用1 MPI进程、1线程，运行FD的0.20、0.15 Å，以及PW的400、600 eV四个条件。完整结果在results/VALIDATION.json，各条件目录中保留原生输出、自己生成的GPW文件和三维密度，图片在plots/。

2 He结构和边界条件#

模型只有一个He，放在8×8×8 Å盒子中心；He有2个电子，本例没有冻结芯层。FD使用非周期边界；PW使用周期晶胞和真空间隔近似孤立中性原子。这两种边界应明确区分。

from ase import Atoms

atoms = Atoms('He', positions=[[4, 4, 4]],
              cell=[8, 8, 8], pbc=False)  # FD
# PW条件使用相同晶胞与位置，pbc=True

本例固定原子位置，没有做几何优化。中心原子的解析力模小于0.02 eV/Å，作为这一对称小模型的受力检查。

3 FD关键输入#

PZ-SIC对轨道密度敏感，因此使用复数轨道，并由直接最小化方法优化：

from gpaw import GPAW, FD
from gpaw.directmin.etdm_fdpw import FDPWETDM

calc = GPAW(
    legacy_gpaw=True,
    mode=FD(force_complex_dtype=True), h=0.15,
    xc='LDA', nbands=2,
    occupations={'name': 'fixed-uniform'},
    eigensolver=FDPWETDM(
        localizationtype='PM_PZ', localizationseed=42,
        functional={'name': 'PZ-SIC', 'scaling_factor': (1.0, 1.0)},
        grad_tol_pz_localization=1e-4),
    mixer={'backend': 'no-mixing'}, symmetry='off',
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-7},
    maxiter=200, txt='FD.txt')
atoms.calc = calc
energy = atoms.get_potential_energy()
forces = atoms.get_forces()

两个缩放数分别控制Hartree和交换关联贡献，不是自旋上下通道。本例使用完整修正；如果改成0.5，就是另一种缩放SIC模型。

这里明确设置nbands=2，给固定均匀占据初始化留出一个空带。SIC修正来自占据轨道，附件不使用这个空带报告带隙或激发能。

4 切换到PW#

PW条件保留同一个直接最小化求解器，将表示改成下面的形式，并使用周期边界；不再传入FD的h参数：

from gpaw import PW

atoms.pbc = True
mode = PW(600, force_complex_dtype=True)
# GPAW参数中使用mode=mode，删除h=0.15；其他SIC设置保持一致

完整的两个PW条件都已包含在附件中。修改输入时，建议保留原目录，在新目录重新计算，避免把FD和PW的旧结果混用。

5 集群实际能量#

下表的原子总能包含PAW原子参考能；它不是反应能或结合能。

FD两档能量差0.013299 eV，PW两档差0.009566 eV。FD 0.15 Å与PW 600 eV的差为0.035226 eV。它们是在本例明确边界和参数下的实际对照，不能当作任意体系的精度保证；自己的工作应针对所需相对能和物性继续检查。

最终总能包含轨道优化后密度泛函部分与SIC部分。附件分别保存这两项，并核对它们之和；不要把普通LDA的旧总能与优化后的SIC修正项直接相加，冒充变分结果。

6 电子密度与数据导出#

左图来自实际FD 0.15 Å计算的全电子密度，右图分别使用FD、PW各自原生网格上的中心线，没有把不同网格未经处理地相减。完整三维密度保留在每个条件的fields.npz中，积分电子数均为2。

导出能量表：

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('he-sic-fd-pw.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'total_atomic_energy_eV',
                'SIC_correction_eV', 'electrons'])
    for r in v['results']:
        w.writerow([r['name'], r['total_atomic_energy_eV'],
                    r['SIC_orbital_correction_eV'], r['AE_ne']])
        print(r['name'], r['forces_eV_A'])
PYCODE

只重画附件原始结果，不重新计算时：

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

图片输出到reference-plots/。附件原始数据在reference/，自己计算的数据在results/；下载包不包含PAW数据集、GPW检查点或程序二进制。

7 环境与其他体系#

本例使用GPAW 26.7.0的legacy FD/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，没有修改公共安装。

换成多轨道体系时，还要检查局域化初猜、多个局域化种子和轨道优化是否足够。自相互作用修正的非自洽原子能量后处理、LCAO变分、FD/PW变分是不同操作路线，应按目标物性和所用表示选择。

方法与接口：GPAW PZ-SIC官方FD、PW和LCAO说明。