GPAW剪刀算符与分区能级修正
用两个氢分子的实际计算演示整体能级平移、占据与空态分开修正、原子区域修正，以及自洽和固定密度后处理的区别。
阅读本篇目录
剪刀算符按用户给定的数值修正单粒子能级，例如把占据态下移、空态上移，或对不同原子区域施加不同平移。本篇用两个小氢分子，说明输入怎么写、修正后的谱怎样读取。
修正量是自己指定的模型参数，不是程序自动得到的GW准粒子修正。 先确认需要模拟的物理含义，再选择数值；不要把示例的±0.5 eV直接套用到实际材料。
1 下载与运行#
下载剪刀算符操作包，保存为gpaw-lcao-scissors.zip。普通作业可在解压目录执行sbatch job.slurm。已经分配好CPU资源时按顺序运行：
python3 -m zipfile -e gpaw-lcao-scissors.zip .
cd gpaw-lcao-scissors
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh
脚本用1个MPI进程、1个线程完成两组实际自洽：普通PBE，以及第一个H₂区域带剪刀修正的PBE；随后对每组密度做四种能级后处理。结果在results/，两张图片在plots/。
2 结构与自洽输入#
每个H₂键长0.74 Å，两个分子沿x相距4 Å，四周加3.5 Å真空。使用LCAO DZP、PBE、0.25 Å网格，4个价电子，自旋配对，不优化结构。
from ase import Atoms
from gpaw.new.ase_interface import GPAW
atoms = Atoms('H4', [[0, 0, 0], [0, 0, 0.74],
[4, 0, 0], [4, 0, 0.74]])
atoms.center(vacuum=3.5)
calc = GPAW(
mode='lcao', basis='dzp', nbands='nao', h=0.25,
xc='PBE', symmetry='off', spinpol=False,
occupations={'name': 'fixed-uniform'},
eigensolver={'name': 'scissors',
'shifts': [(-0.5, 0.5, 2)]},
convergence={'energy': 1e-6, 'density': 1e-5,
'eigenstates': 1e-7},
maxiter=200, txt='region-scissors.txt')
atoms.calc = calc
atoms.get_potential_energy()
calc.write('state.gpw', mode='all')
这里明确使用GPAW新后端。shifts中的一项是(占据态平移eV, 空态平移eV, 连续原子个数)：(-0.5,0.5,2)作用于前两个原子，也就是第一个H₂；后面未指定的原子不加修正。
如果接着写第二项，它从下一批原子开始。例如[(0.2,0.2,2),(-0.2,-0.2,2)]表示第一个分子整体上移、第二个整体下移。原子计数按当前结构顺序，不能当作任意原子编号列表。
区域修正可能破坏原子间的对称映射。本例关闭对称性，避免程序把应当不同的区域视为等价。
3 固定密度后处理#
对已经自洽的计算，可以只重新构造和对角化修正后的哈密顿量，不再做SCF：
from gpaw.lcao.scissors import non_self_consistent_scissors_shift
levels = non_self_consistent_scissors_shift(
[(-0.5, 0.5, 4)], calc.dft)
print(levels.shape) # 自旋数、不可约k点数、带数
print(levels[0, 0]) # eV
本例的四种输入为：
| 名称 | 输入 | 含义 |
|---|---|---|
| zero | [] |
当前密度的未修正哈密顿量 |
| uniform | [(0.2,0.2,4)] |
全部4原子的能级整体上移0.2 eV |
| global-gap | [(-0.5,0.5,4)] |
完整体系的占据态下移0.5、空态上移0.5 eV |
| region | [(-0.5,0.5,2)] |
只对前两个原子区域加修正 |
区域的局域轨道与分子能级之间可以杂化，所以分区修正不是简单地把某几条能带机械移动±0.5 eV。
4 实际谱线#
下面取普通PBE密度，列出最低四个带，单位eV，编号从0开始；带0、1占据，2、3为空态：
| 带 | 未修正 | 整体+0.2 | 整体开隙 | 第一个H₂区域 |
|---|---|---|---|---|
| 0 | −10.36089014 | −10.16089014 | −10.86089014 | −10.83645840 |
| 1 | −10.30945566 | −10.10945566 | −10.80945566 | −10.33384249 |
| 2 | 2.10495603 | 2.30495603 | 2.60495603 | 2.22581608 |
| 3 | 2.45891313 | 2.65891313 | 2.95891313 | 2.83803474 |
整体修正前HOMO–LUMO差为12.41441169 eV，整体开隙后为13.41441169 eV。它只是这组LCAO能级的差，不能直接解释为氢分子的光学吸收峰或实验电离能。
在带区域修正的自洽解上，原生自洽求解器的谱与同样区域参数的后处理谱一致；将后处理改成[]得到的是同一密度的未修正哈密顿量谱，不是“仍然保留原修正”。自洽与固定密度后处理也不能混称为同一个计算流程。
查看保存的谱：
python3 -c "import json; v=json.load(open('results/VALIDATION.json')); print([(r['name'],r['shifted']) for r in v['results']])"
5 占据读取和输出文件#
附件读取原始占据的副本，再在副本上乘自旋简并因子。GPAW26.7.0这一接口版本中，默认加权占据读取可能原地改变后续计算使用的数组；在还要调用剪刀后处理时，按下面方式读取：
raw = calc.get_occupation_numbers(raw=True).copy()
weighted = raw * calc.dft.ibzwfs.spin_degeneracy
本例只有Γ点，k点权重为1。多k点时还要按实际k点权重求和，不能把这里的单k点写法直接用于所有体系。
source env.sh
python3 -c "import numpy as np; d=np.load('results/plain/fields.npz'); print('raw occupation sum:',d['raw_occupations'].sum()); print('weighted electrons:',d['weighted_occupations'].sum())"
ls -lh plots/scissors-eigenvalues.png plots/scissors-structure-shifts.png
输出包含原始/修正本征值、实际占据、赝密度、结构与JSON结果。原样复跑可用verify.sh；换体系后需要自己检查带数、区域顺序、占据与所选修正量。
6 适用范围#
本篇演示自洽区域修正和非自洽能级后处理，没有用修正后的总能量比较结构稳定性，也没有验证相应的原子力。需要GW、激子光谱或受激电子结构时，应选择相应方法，而不是仅靠剪刀参数给出结论。
环境为GPAW26.7.0新后端、GCC13.3、OpenMPI5.0.10rc2、MKL2023.2、Libxc7.0、Python3.11.16。该功能在GPAW原生说明中仍标为开发中的功能，实际使用应按所装版本核对接口。