VASP电子密度与NCI氢键图

已经得到一个吸附结构或分子二聚体，想看看弱相互作用主要位于哪里，可以计算电子密度，再画NCI图。本例用实际优化过的两个水分子，展示从结构、密度到图像的流程。

1. 先获得合理的结构#

本例水二聚体采用PBE-D3(BJ)、PAW64、550 eV、Gamma单点、EDIFF=1E-6，在约17.61×15.14×15.52 Å的真空盒中优化原子位置。优化后的最大原子力为0.0184 eV/Å，满足0.02 eV/Å要求。

实际O···O距离为2.885 Å，分子之间最近O···H距离为1.909 Å。先看几何是否合理，再分析密度；图像不能代替结构优化。

2. 对优化结构做静态电子密度计算#

把优化得到的CONTCAR作为新目录POSCAR，保持泛函、PAW、截断能和k点一致。附件已放入本例优化结构。静态步骤使用：

IBRION = -1
NSW = 0
LCHARG = .TRUE.
IVDW = 12

准备合法POTCAR后，进入密度计算目录提交：

sbatch run.sh

完成后得到CHGCAR。本例NCI使用其中的平滑价电子密度；不要把它当作包含全部冻结芯电子的密度做总电子数比较。色散修正用于本例几何和总能量，NCI图本身来自输出电子密度。

3. 用critic2生成NCI数据#

把自己完成的CHGCAR放入post目录。fragment1.xyz和fragment2.xyz分别定义两个水分子；换体系时须重写这两个片段和nci.cri中的观察区域，不能沿用本例坐标。

sbatch post.sh

输入先明确加载CHGCAR，再调用nciplot，并限制在两个片段之间。这样实际使用VASP密度，而不是在未加载场时自动采用原子密度叠加。

crystal CHGCAR
load CHGCAR
reference 1
nciplot
 oname water-dimer
 interpolation fourier
 fragment fragment1.xyz
 fragment fragment2.xyz
 nstep 80 80 80
 cutoffs 0.05 2.0
 cutplot 0.04 0.5
 nochk
end

附件还包含限定局部观察盒的cube一行。处理后得到散点数据water-dimer.dat、RDG和带符号密度的cube，以及water-dimer.vmd。可在VMD加载该状态文件查看三维等值面。critic2的NCI输入与输出说明

4. 如何看这张实算图#

左图横轴是密度Hessian第二本征值的符号乘以电子密度，纵轴是约化密度梯度RDG。右图截取穿过氢键附近的平面，只显示RDG≤0.5的区域，原子作平面投影；右侧水分子的两个H在投影中重叠。

本例低RDG区域位于分子之间、带符号密度为负的一侧，与该几何中的氢键吸引区相符。正值区域通常提示排斥特征，接近零的区域可对应较弱作用；颜色还取决于自己设定的色标和截断值，不应只凭蓝色面积比较两种材料的结合能。

本图采用240³的VASP密度网格，后处理局部网格80³，共输出31216个符合散点截断条件的点。这不是独立样本数，也不是氢键数量。密度网格从180³加密到240³时，同一局部后处理网格上RDG≤0.5的点从273个变为306个，公共区域的带符号密度最大差约7.86×10⁻⁵ e/bohr³。主要氢键区位置和负号保持，但边界对网格仍有敏感性，因此不把区域体积作为精确物理量。

5. 提取数据与避免单位错误#

module load ase/3.29.0
python extract_nci.py

提取脚本读取后处理输出，检查网格和数值，并给出结构距离及图像筛选条件。critic2的dens.cube存储的是100倍带符号密度，dat第一列直接给带符号密度；自行画图时必须区别处理，否则横轴会错100倍。

VASP使用6.6.1 P100构建（NVHPC25.3、CUDA12.8和配套MPI）；后处理使用critic2 1.4（GCC/GFortran13.3、OpenMP）。图中使用电子密度的导数，密度网格过粗会带来噪声，应先检查网格再讨论细节。NCI用于定性定位作用区域，结合能仍应按一致参照单独计算。

下载输入、结果与提取脚本（ZIP）

附件包含输入、片段、后处理与提取脚本、实际参考图和数据摘要，不含POTCAR或电子密度大文件。