GPAW约束DFT与局部电荷

约束DFT可以在选定原子区域内指定局部电子数或自旋。本例固定H2结构，给第一个H设置0、+0.2或−0.2 e的局部电荷目标，观察约束态能量和拉格朗日乘子如何变化。

1 下载并计算#

下载约束DFT实操包，保存为gpaw-cdft.zip。新计算使用CPU资源；从登录节点提交可用sbatch job.slurm。在已分配的CPU资源中：

python3 -m zipfile -e gpaw-cdft.zip .
cd gpaw-constrained-dft
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

例程为1 MPI进程、4线程，结果在results/VALIDATION.json，实际权重与电子密度网格在results/*-FIELDS.npz。总体系始终有2个电子、保持中性；局部约束是在两个区域之间重新分配电子，不是改变整个体系的总电荷。

2 模型与关键输入#

H–H距离固定为1.5 Å，实际盒为8×8×10 ų，原子位置为(4,4,4.25)与(4,4,5.75) Å。采用PBE、FD间距0.25 Å，显式开启两个自旋通道，初始总磁矩为零；本例不做离子优化。

base = GPAW(
    legacy_gpaw=True, mode='fd', h=0.25, xc='PBE', spinpol=True,
    occupations=FermiDirac(0.01),
    convergence={'energy':1e-7,'density':1e-6,'eigenstates':1e-9},
    maxiter=200
)
cdf = CDFT(
    base, atoms, charge_regions=[[0]], charges=[0.2],
    charge_coefs=[1.0],
    mu={'H':0.5}, Rc={'H':20.0},
    minimizer_options={'gtol':1e-4,'maxiter':35}, tol=1e-4,
    method='BFGS', compute_forces=False
)
atoms.calc = cdf
energy = atoms.get_potential_energy()

原子索引从0开始，charge_regions=[[0]]选择第一个H。charges=[0.2]表示该区域少0.2个电子，因此目标区域电子数为0.8，而不是1.2。charge_coefs给出约束系数的初始值，程序还会优化它，不是最终固定施加的电压。

3 先明确“局部电荷”的定义#

这里使用原生高斯Hirshfeld型分区。高斯宽度mu和截断Rc均以Å输入；得到的区域电荷不是Bader电荷，也不是形式氧化数。

本例使用20 Å截断，是为了让高斯支撑覆盖这个有限盒中的全部网格，再逐点检查两个原子的权重和为1；它不是其他体系通用的最佳参数。较大体系应根据盒与分区定义选择范围，并检查被遗漏的电子密度。

实际对照中，3.1 Å截断遗漏约0.0014–0.0022个电子的尾部，导致正负局部电荷目标不是严格互补，受约束能量相差约6.865 meV。覆盖全盒后，该差异降为约9.3×10⁻⁸ eV。这个对照说明，不能只看优化器显示收敛，还要检查分区是否包含了要约束的电荷。

4 实际结果#

同样写“+0.2 e”，改变高斯宽度也会改变所定义的约束态。因此正式应用应说明分区方式和参数，不能把不同分区得到的能量当作完全相同状态来比较。

KS_constrained_energy_eV为该受约束密度对应的KS能量；CDFT_native_augmented_energy_eV还带有原生接口中的约束势项。比较状态能量时应先明确采用哪个定义，不能把两列混在同一能量差里。

5 提取电荷、约束系数与网格#

source env.sh
"$YS_PY" - <<'PYCODE'
import json, csv, numpy as np
from pathlib import Path
f=json.loads(Path('results/VALIDATION.json').read_text())
with open('h2-cdft-results.csv','w',newline='') as out:
    w=csv.writer(out)
    w.writerow(['target_charge_e','width_A','population_e','KS_energy_eV','lagrange_eV'])
    for r in f['results']:
        print(r['name'],'population=',r['region_population_e'],
              'target=',r['region_target_e'],'residual=',r['constraint_residual_e'])
        w.writerow([r['charge_target_e'],r['gaussian_width_A'],r['region_population_e'],
                    r['KS_constrained_energy_eV'],r['lagrange_eV']])
r=f['results'][1]
a=np.load('results/'+r['name']+'-FIELDS.npz')
print('actual weight grid:',a['weight'].shape)
np.save('H0-partition-weight.npy',a['weight'])
PYCODE

bash plot.sh --reference

后一条重绘归档图到reference-plots/。完整程序还把全网格权重、加上PAW增广修正后的区域电子数，以及约束能量关系独立核对，避免仅由显示的目标值推断约束已经满足。

本例演示固定结构的部分电荷约束，不计算电子转移速率、非绝热耦合或激发态寿命。自旋约束、分区电荷差和耦合计算需要相应模型与输入。程序为GPAW 26.7.0，GCC 13.3、HPC-X/OpenMPI 5、MKL 2023.2、Libxc 7.0.0，PAW数据26.8.0；原理与参数见GPAW约束DFT说明。