GPAW水分子变分自相互作用修正

本篇演示水分子的变分PZ-SIC：先得到同一结构、同一基组的PBE结果，再优化修正泛函的轨道，最后比较能量、力和密度。它与原子非自洽SIC有区别：这里实际改变电子态；PAW冻结芯层仍保持冻结，没有额外做原子芯层NSC-SIC修正。

PZ-SIC依赖每条轨道的密度，轨道之间做酉变换也可能改变能量，因此局域化和直接变分优化是计算的一部分。不能简单地在普通PBE结果后加一个常数，就代替这条流程。

1 下载并按顺序运行#

下载水分子变分SIC算例包，保存为gpaw-water-sic.zip。从登录节点提交时，解压后使用sbatch job.slurm；已经进入分配好的CPU资源时，可以运行：

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

完整流程依次运行5个条件：PBE基准、局域化种子42与43的SIC、从种子42收敛态出发的两个位移结构。输出在results/，图片在plots/。两次位移各自重新收敛电子态，并非仅对一条拟合曲线求导。

2 结构、基组与复数轨道#

结构使用ASE的水分子，三个方向不周期，目标真空4 Å；adjust_cell按0.25 Å网格适配盒子。后面的位移计算保持这个晶胞不变：

from ase.build import molecule
from gpaw.utilities.adjust_cell import adjust_cell

atoms = molecule('H2O')
adjust_cell(atoms, 4.0, 0.25)
positions = atoms.positions.copy()
cell = atoms.cell.copy()

使用复数LCAO轨道、dzp基组、PBE以及固定均匀占据；H2O有10个物理电子，其中O 1s²冻结，显式价电子为8。普通PBE先按相同几何和数值参数收敛，保存PBE.gpw。

from gpaw import GPAW, LCAO

calc = GPAW(
    legacy_gpaw=True,
    mode=LCAO(force_complex_dtype=True), basis='dzp',
    h=0.25, xc='PBE', nbands='nao',
    occupations={'name': 'fixed-uniform'}, symmetry='off',
    convergence={'energy': 1e-6,
                 'density': 1e-5,
                 'eigenstates': 1e-6}, maxiter=200)

force_complex_dtype=True允许复数轨道；nbands='nao'使用与原子轨道基组相应的带数。本例没有用这些空态来报告激发能或带隙。

3 从PBE切换到变分PZ-SIC#

附件通过LCAOETDM优化轨道，采用PM_PZ局域化。两个0.5分别缩放Hartree和交换关联修正，不是两个自旋通道：

from gpaw.directmin.etdm_lcao import LCAOETDM

optimizer = LCAOETDM(
    localizationtype='PM_PZ', localizationseed=42,
    need_init_orbs=False, need_localization=True,
    functional={'name': 'PZ-SIC',
                'scaling_factor': (0.5, 0.5)},
    subspace_convergence=1e-5)
calc = GPAW('PBE.gpw', legacy_gpaw=True,
            eigensolver=optimizer,
            mixer={'backend': 'no-mixing'}, txt='SIC.txt')
atoms = calc.get_atoms()
calc.initialize(atoms)
atoms.calc = calc
energy = atoms.get_potential_energy()
forces = atoms.get_forces()

这里的calc.initialize(atoms)是原生初始化调用：读取检查点并更换求解器后，先建立新求解器，再请求能量或改变位置。尤其在位移重启时应保留这一行。附件没有修改GPAW安装文件。

4 集群实际能量与种子比较#

两次局域化得到的总能差为7.28×10⁻⁸ eV。它说明本例这两个起点落到了能量非常接近的解，没有证明任意初猜都能找到全局最低轨道解。

SIC总能变化不等于表中的轨道修正：轨道重新优化后，原密度泛函部分也会变化。应读取最终完整总能，不要把PBE基准与修正项直接相加当作变分结果，也不能通过跨泛函的绝对能量大小判断准确性。

5 用两个真实位移核对力#

附件固定晶胞和其他坐标，仅把第一个H的y坐标分别移动−0.002、+0.002 Å。两次都读取实际种子42的SIC检查点，保留已有局域化轨道，再进行电子态优化：

F(H1,y) ≈ −[E(y+0.002)−E(y−0.002)] / 0.004 Å。

本例按0.02 eV/Å比较标准检查两者一致性。这不是结构优化已经满足0.02 eV/Å停止条件；本例结构保持固定。需要优化结构时，可在合适的ASE优化器中调用同一SIC计算器，并另行检查基组、网格、盒子和轨道优化精度。

6 提取数据并画密度差#

左图是PBE的全电子密度，右图逐网格相减得到n(SIC)−n(PBE)。两份密度使用相同晶胞和网格，积分电子数均为10；因此右图描述同一分子的电子重分布，不是净增减了电子。

导出完整能量和受力表：

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('water-sic.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'total_energy_eV', 'SIC_correction_eV',
                'H1_y_force_eV_A', 'electrons'])
    for r in v['results']:
        w.writerow([r['name'], r['energy_eV'],
                    r['SIC_orbital_correction_eV'],
                    r['forces_eV_A'][1][1], r['all_electron_integral']])
for item in v['items'][-2:]:
    print(item)
PYCODE

只重画附件中的实际数据：

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

reference/包含全部5个条件的三维密度、晶胞、原子位置和能量受力表；自己计算的对应文件在results/。下载包不包含PAW数据集或GPW检查点。

7 环境与后续修改#

本例使用GPAW 26.7.0的legacy LCAO接口、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，任务使用1 MPI进程、1线程。

改为自己的分子时，先确定电荷、自旋和占据，再检查dzp基组、0.25 Å网格及4 Å真空是否够用。复杂体系宜比较多个局域化种子，结合能量和力检查电子态优化。这里演示的是缩放系数0.5的具体方法，不能把它写成未缩放的完整PZ-SIC。

方法与接口：GPAW原生PZ-SIC说明及FD、PW、LCAO示例。本例实际执行LCAO分支；其他表示需要各自的输入和数值检查。