GPAW范围分离泛函与垂直电离能

本篇用固定键长的H₂和H₂⁺，练习GPAW的Yukawa范围分离泛函。操作顺序是：准备同一几何，分别自洽计算中性和阳离子，用总能差得到垂直电离能，再与中性体系的负HOMO能量比较。

附件包含两档LCY-PBE参数和一档CAMY-B3LYP。它们用于演示方法与输入，不构成完整的参数优化，也没有计算吸收谱或可靠的空带能量。

1 下载并计算#

下载范围分离泛函算例，保存为gpaw-range-separated.zip。在登录节点解压后，用sbatch job.slurm提交；已经进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-range-separated.zip .
cd gpaw-range-separated
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

程序使用1个MPI进程、1个线程。它完成六组范围分离泛函自洽计算；CAMY的两个电荷态各增加一次PBE初猜计算。结果保存在results/，图片在plots/。计算比普通PBE的小分子例子更慢，请等待作业正常结束。

2 结构和电荷态#

H—H键长固定为0.74 Å，使用12×12×12.74 Å非周期盒子，目标网格间距0.25 Å：

from ase import Atoms

atoms = Atoms('H2', positions=[[0, 0, 0], [0, 0, 0.74]])
atoms.center(vacuum=6.0)

中性H₂有2个电子，使用自旋配对；H₂⁺设置charge=1，保留1个电子并使用自旋极化。阳离子的初始磁矩可设为每个H各0.5。两者原子位置完全相同，因而总能差对应本几何的垂直电离能，不是两个分别优化结构之间的绝热电离能。

这里H没有冻结芯层。换成包含芯层的其他元素时，还要核对该范围分离方法所需的PAW数据条件，不能只替换元素名。

3 LCY-PBE输入#

下面是中性体系、默认γ=0.75 a₀⁻¹的实际设置：

from gpaw import GPAW, FermiDirac, Mixer
from gpaw.eigensolvers import RMMDIIS
from gpaw.poisson import PoissonSolver

calc = GPAW(
    legacy_gpaw=True, mode='fd', xc='LCY-PBE',
    h=0.25, nbands=2, charge=0, spinpol=False,
    symmetry='off', eigensolver=RMMDIIS(),
    occupations=FermiDirac(0.0, fixmagmom=False),
    poissonsolver=PoissonSolver(use_charge_center=True),
    mixer=Mixer(0.2, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-6},
    maxiter=240, txt='LCY-neutral.txt')
atoms.calc = calc
energy = atoms.get_potential_energy()
homo = calc.get_eigenvalues(spin=0)[0]
calc.write('state.gpw', mode='all')

本例只有一个占据空间轨道，因此spin=0的第0个本征值就是所用HOMO。对其他体系，不要直接照搬这个索引，应按电子数、自旋和实际占据识别HOMO。

切换到γ=0.50时，把泛函参数改为：

xc = 'LCY-PBE:omega=0.5'

接口参数叫omega，这里对应γ，单位是逆玻尔半径 a₀⁻¹，不是逆埃。LCY中的Y表示Yukawa形式；不能把它当成HSE中相同数值的屏蔽参数。

这里使用legacy有限差分实现、Γ点和非周期边界，并明确选择RMMDIIS。这组条件不要未经核对地推广到周期晶体或其他求解器。

4 CAMY-B3LYP的PBE初猜#

CAMY-B3LYP采用默认γ=0.34 a₀⁻¹。本例先用相同电荷、自旋、网格和盒子完成PBE，再从检查点切换泛函：

# shared是与当前电荷态相匹配的FD参数，见附件run.py
# 本例CAMY仅需要1个空间轨道：nbands=1
seed = GPAW(xc='PBE', txt='PBE-seed.txt', **shared)
atoms.calc = seed
atoms.get_potential_energy()
seed.write('PBE-seed.gpw', mode='all')

calc = GPAW(
    'PBE-seed.gpw', legacy_gpaw=True, xc='CAMY-B3LYP',
    eigensolver=RMMDIIS(),
    poissonsolver=PoissonSolver(use_charge_center=False),
    mixer=Mixer(0.05, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-6},
    maxiter=240, txt='CAMY.txt')
atoms = calc.get_atoms()
atoms.calc = calc
energy = atoms.get_potential_energy()

附件给出shared的完整定义，按中性与阳离子分别设置；阳离子使用spinpol=True和固定磁矩占据。这里保留PBE初猜、较小的密度混合系数，以及放在盒中心的中和高斯。收敛门槛与LCY例子相同。

CAMY的一个带足够容纳本例所需占据态；它不提供LUMO或带隙。LCY例子虽保留第二个带，本文同样不把空带当成已经计算可靠的激发态结果。

5 提取垂直电离能#

定义：

IP = E(H₂⁺; H₂几何) − E(H₂)

实际结果如下：

中性和阳离子使用相同原子参考能，总能差中该项抵消。不同泛函的绝对总能、IP和HOMO承担不同含义，不能把表中某一个较小的绝对总能视为该泛函更准确的依据。

两个量的残差IP + HOMO分别约为0.942873、1.763756和3.280240 eV。它们并未自动相等。只有两档LCY参数，不能证明最优γ已经找到；CAMY还改变了泛函本身，也不能把三行直接当成同一泛函的参数扫描。

如果研究目标确实是调参，需要在相同泛函、几何、边界和数值精度下增加γ采样，再检查残差及所需物性。本篇不把这些小模型数值当成实验精度保证。

6 电子密度和自旋#

前两图分别来自LCY、CAMY中性H₂的全电子总密度；第三图是CAMY阳离子的上、下自旋密度之差。中性积分为2个电子，阳离子为1个电子，自旋电子数差为1。图像由实际原生网格数组绘制，色标下限以下统一显示为最低颜色。

读取总密度时，使用get_all_electron_density(gridrefinement=1)；要分开读取物理自旋密度，则分别指定spin=0和spin=1。当前legacy接口的自旋配对计算中，单独spin=0只有半份总密度，不能把它直接用于总电子数积分。

7 导出IP与HOMO表格#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('rsf-ip-homo.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['setting', 'vertical_IP_eV', 'negative_HOMO_eV',
                'IP_plus_HOMO_eV', 'squared_residual_eV2'])
    for name, row in v['summary'].items():
        values = [name, row['vertical_IP_eV'], row['negative_HOMO_eV'],
                  row['IP_plus_HOMO_eV'], row['IP_HOMO_residual_squared_eV2']]
        w.writerow(values)
        print(*values)
PYCODE

各条件的fields.npz保留完整总密度、两份物理自旋密度、占据、晶胞和原子位置。原生输出在各条件的native.txt，CAMY还有PBE-seed.txt，方便查看两阶段收敛。

只重画附件参考结果时：

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

自己的结果在results/，原始参考结果在reference/。附件不包含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。setup.sh加载gpaw/26.7.0-8173m。

本篇只使用占据态和总能。需要范围分离泛函的非占据态或线性响应TDDFT时，还有Fock算符对非占据态的作用、PAW数据和响应空间等设置，不能仅用本文的HOMO计算替代。

方法与版本条件参考：GPAW官方范围分离泛函说明。