VASP与Hirshfeld原子电荷
用NaCl的三组真实密度网格，学习全电子Hirshfeld划分和电荷积分检查。
阅读本篇目录
Hirshfeld电荷按自由原子参考密度的权重分配电子密度。它与Bader分区的定义不同，因此不同方法得到的原子电荷不必相同。本例用NaCl展示原始Hirshfeld划分，并通过120³、240³、360³网格检查积分精度。
1. 准备静态密度计算#
本例沿用NaCl两原子原胞，PBE PAW64、ENCUT=550 eV、EDIFF=1E-6、7×7×7 Gamma网格，固定结构计算密度。输入包含：
LCHARG = .TRUE.
LAECHG = .TRUE.
NGXF = 360
NGYF = 360
NGZF = 360
LCHARG生成CHGCAR，LAECHG生成PAW全电子重构相关文件。这里用于划分的全电子密度是AECCAR0+AECCAR2，不能用仅含价电子的CHGCAR代替后仍减去原子核电荷Z。
三个目录分别准备与POSCAR顺序相符的合法POTCAR，再各自执行：
sbatch run.sh
参考120³是本例标准设置的实际输出网格，240³与360³显式指定细网格。核附近密度变化很快，普通能量已收敛不代表全电子密度的均匀网格积分也足够精确。不同材料不应照搬360³作为固定答案。
2. 用critic2分配全电子密度#
等待对应VASP计算完成后，进入postprocessing的对应网格目录，再执行sbatch run.sh。脚本连接该目录的CHGCAR、AECCAR0、AECCAR2，并调用critic2：
crystal CHGCAR
load AECCAR0
load AECCAR2
load as "$1+$2"
reference 3
integrable clear
integrable 3 f name AllElectron
hirshfeld
本例执行原始Hirshfeld方法，使用critic2的原子参考密度；没有使用迭代Hirshfeld-I，也没有执行DDEC。权重划分与核电荷扣除必须配套，方法区别见critic2原子积分说明。
critic2 1.4由GCC/GFortran13.3构建，在CPU上进行后处理；VASP6.6.1-p100由NVHPC25.3/CUDA12.8及配套MPI构建。全部算例经Slurm执行。附件后处理脚本使用现有GPU分区作为节点入口，但critic2这一步并没有GPU加速。
3. 提取原子电荷，先看总电子数#
module load ase/3.29.0
python extract_hirshfeld.py
脚本从三个hirshfeld.cro中读取每个原子的电子数N_A，计算q_A=Z_A−N_A。全电子NaCl应有11+17=28个电子；本例价电子数则为1+7=8，这两种总数不要混淆。
| 密度网格 | 全电子积分 | 28−积分 / e | Na电荷 / e | Cl电荷 / e |
|---|---|---|---|---|
| 120³ | 28.3786294 | -0.3786294 | +0.153293 | -0.531922 |
| 240³ | 28.0278355 | -0.0278355 | +0.202675 | -0.230510 |
| 360³ | 28.0061047 | -0.0061047 | +0.205573 | -0.211678 |
360³结果为Na +0.205573 e、Cl -0.211678 e。表格保留数值供复核，不能把所有小数位都视为材料预测精度。120³的总电子积分误差明显，直接读出其原子电荷会造成误判；比较240³与360³时还应结合原子电荷变化和总电荷误差判断。
240³到360³时，Cl电荷仍变化约0.018833 e，Na变化约0.002898 e；这些结果支持约±0.2 e量级的入门比较，不能声称各原子电荷已经达到0.001 e精度。
这里没有把电荷强行归一化到总和为零。强行平均分摊积分误差会掩盖粗网格问题，也不等价于重新完成密度划分。
4. 换成自己的材料#
先检查同一组输出的结构、晶胞和网格一致，再逐步加密并比较每个关注原子的电荷，而不只看程序是否正常结束。改变参考原子、赝势或划分方法时，要重新说明定义。
Hirshfeld电荷是密度划分量，不是整数氧化态；NaCl的形式氧化态±1不能用来要求本表也必须等于±1。研究趋势时尽量固定方法和数值设置，并与Bader及成键分析等不同定义的量区分。