GPAW水二聚体与基组重叠误差修正

用局域基组计算分子间相互作用时，二聚体中的一个分子可以“借用”另一个分子的基函数，使二聚体相对单独分子显得结合得更强。这个效应称为基组重叠误差，BSSE。Counterpoise方法保留另一分子的基函数，同时去掉它的真实原子作用，构造可比较的片段参考能。

1 先跑完整例子#

下载水分子LCAO操作包，保存为gpaw-lcao-workflows.zip。普通作业可用sbatch job.slurm提交。下面的命令在已分配好的CPU资源中执行：

python3 -m zipfile -e gpaw-lcao-workflows.zip .
cd gpaw-lcao-workflows
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

结果在results/；本篇使用其中AB、A-own、B-own、A-ghost-B、B-ghost-A五组。另一个私有基组水分子和局域轨道步骤，见基组生成与局域轨道。

2 固定结构，建立五个目录#

本例使用ASE自带S22水二聚体的固定坐标，四周增加4 Å真空，不优化几何。A是原子0、1、2，B是3、4、5；每组均在同一个晶胞、同一个坐标系中计算，删除原子后不重新居中。

from ase.data.s22 import create_s22_system

AB = create_s22_system('Water_dimer')
AB.center(vacuum=4.0)
AB.pbc = False
A = AB[:3].copy()
B = AB[3:].copy()

ghost是虚原子中心：保留该元素的LCAO基函数，不保留真实原子的PAW核、电子数和原子势。不要用任意小核电荷或改元素符号来代替。

3 普通与ghost输入怎样写#

全部使用PBE、已安装DZP、0.22 Å网格、全部基组带、非周期Γ点和自旋配对。核心计算器如下：

from gpaw import GPAW, Mixer

calc = GPAW(
    legacy_gpaw=True, mode='lcao', basis='dzp',
    xc='PBE', h=0.22, nbands='nao', spinpol=False,
    occupations={'name': 'fixed-uniform'},
    mixer=Mixer(0.2, 5, 50),
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-7},
    maxiter=200, txt='native.txt')

计算A、保留B的基组时，仍传入全部6个位置，只把B对应的中心指定为ghost：

A_with_ghost_B = AB.copy()
A_with_ghost_B.calc = GPAW(
    legacy_gpaw=True, mode='lcao', basis='dzp',
    xc='PBE', h=0.22, nbands='nao', spinpol=False,
    occupations={'name': 'fixed-uniform'},
    setups={3: 'ghost', 4: 'ghost', 5: 'ghost'},
    convergence={'energy': 1e-6, 'density': 1e-5,
                 'eigenstates': 1e-7},
    txt='A-ghost-B.txt')

反过来计算B时，使用setups={0:'ghost',1:'ghost',2:'ghost'}。本例原子编号从0开始。自己的片段含更多原子时，应按实际编号修改，不能把这个三原子范围照搬到别的结构。

4 实际能量与公式#

采用同一个GPAW能量字段，单位eV：

未修正相互作用能为E(AB) − E(A-own) − E(B-own)，得到 −0.35666962 eV。

Counterpoise修正为E(AB) − E(A-ghost-B) − E(B-ghost-A)，得到 −0.20371010 eV。两者相差 +0.15295952 eV；按这个定义，负值表示相对于固定形状的两个分子具有吸引作用。

两个分子在输入中不是完全相同的孤立几何，因此A-own与B-own能量不必相同。这五个能量给出的是固定片段几何的相互作用能；没有另算孤立分子弛豫，所以不是已经包含形变能的完整结合能。

从结果直接提取：

python3 -c "import json; v=json.load(open('results/VALIDATION.json')); print([(r['name'],r['energy_eV'],r['valence_electrons'],r['nao']) for r in v['results'] if r['name']!='A-private']); print(v['summary'])"

5 密度和自己的计算#

附件为所有情况保存真实的赝密度、占据、基组重叠矩阵和本征值；仅有真实原子的情况还保存全电子密度网格。当前GPAW26.7.0的全电子网格导出接口不适用于含ghost的结构，含ghost时不要直接调用该接口；能量计算和counterpoise流程不依赖这一步。

source env.sh
python3 -c "import numpy as np; d=np.load('results/A-ghost-B/fields.npz'); print(d.files); print('valence electrons:',d['weighted_occupations'].sum()); print('real atom indices:',d['real_indices'])"
ls -lh plots/lcao-water-counterpoise.png

自己的计算仍需收敛基组、网格、真空层和相关能量差。尤其是弱相互作用体系，还要选择适当的交换关联与色散方法；counterpoise不能替代这些选择，也不是把PBE自动变成高精度氢键方法。

环境为GPAW26.7.0，GCC13.3/OpenMPI5.0.10rc2/MKL2023.2、Python3.11.16。原生ghost设置参考GPAW LCAO与基组说明。