VASP 全电子密度与键临界点

用VASP算出电子密度后，想知道原子之间是否存在密度梯度路径、键临界点在哪里，可以交给Critic2继续分析。本例使用集群实际计算的H2全电子重构密度，完成从输入、提交、密度相加到临界点、单位与图片的流程。

1. 下载实例，先看图再重画#

下载H2全电子密度与QTAIM实操包，在新目录解压：

unzip vasp-H2-QTAIM-wave297-*.zip
cd vasp-H2-QTAIM-wave297
bash preview.sh
bash run.sh

preview/重画参考图；bash run.sh从真实全网格数值重建CUBE，实际运行三次Critic2，再提取和画图。work/保存新的CUBE与Critic2输出，results/data/为表格和全部选定数值，plots/为重画的三张图，随包figures/可直接查看。输出目录受保护，重复运行请使用新解压目录。

这一步使用已分配资源内1个CPU后处理，不重新运行VASP。数值包保存完整的120³、144³、192³选定密度，没有抽样或强制归一化；不含POTCAR、WAVECAR或原始完整VASP输出。生成的CUBE与本例用于分析的源CUBE逐字一致。

2. 准备结构与全电子密度#

本例几何取自已有气相H2优化，H–H约0.7501 Å，15 Å立方真空盒。后续密度使用VASP6.6.1 GPU/HDF5固定结构单点：PBE、PAW PBE64、550 eV、EDIFF=1E−6、Gamma单点，V10016G、1个MPI进程和1个OpenMP线程。GPU构建环境为NVHPC25.3、CUDA12.8及配套HPC-X；Critic2 1.4为GCC/GFortran13.3的CPU构建。VTST1040是Perl脚本，不需要编译；包内ASE3.29.0/NumPy/绘图使用pymatgen/2026.9.24环境。

自己的体系先优化结构，再做一致参数的静态密度。关键设置：

IBRION = -1
NSW = 0
LCHARG = .TRUE.
LAECHG = .TRUE.
NGXF = 192
NGYF = 192
NGZF = 192

AECCAR0为冻结芯密度，AECCAR2为重构的价电子密度；全电子参考为两者之和。H在这里没有冻结芯电子，因此全电子总数和价电子总数都是2。对于其他元素，不能将全电子积分人口直接减POTCAR的ZVAL；原子净电荷还要配套所用分区和电子计数定义。

mkdir H2-density
cp inputs/H2-request192-actual192/INCAR H2-density/
cp inputs/H2-request192-actual192/POSCAR H2-density/
cp inputs/H2-request192-actual192/KPOINTS H2-density/
cd H2-density
# 按H元素准备本人有权使用的POTCAR

本例192³输入读取前一条件的匹配WAVECAR。如果已有同结构、晶胞、基组和PAW的收敛波函数，将其复制到计算目录；首次独立计算没有该文件时，在INCAR中将ISTART改为0、ICHARG改为2，再提交。

保存下面的submit.sh后提交。核时模板使用V10016G分区；套餐用户通过yskit生成自己可用分区的模板，保留一致的输入参数：

#!/bin/bash
#SBATCH --partition=gpu-8v100-16g-hs
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --gres=gpu:v100-16g:1
#SBATCH --time=00:20:00
source /etc/profile
module purge
module load vasp/6.6.1-v100-hdf5
export OMP_NUM_THREADS=1
ulimit -s unlimited
mpirun --bind-to none -np 1 vasp_std > vasp.log

sbatch submit.sh
tail vasp.log
ls AECCAR0 AECCAR2 CHGCAR

不要只按INCAR中要求的值给网格命名。本例第一份输入要求96³，程序按照基组与FFT要求实际提高为120³；另两份实际为144³、192³。可检查原生输出：

grep 'dimension x,y,z' OUTCAR

后面的表格和图片一律使用真实的120³/144³/192³。

3. 用VTST相加，再让Critic2读取#

VASP完成后，在已分配资源中执行：

module load vtstscripts/1040 critic2/1.4
chgsum.pl AECCAR0 AECCAR2
mv CHGCAR_sum density.CHGCAR

把以下内容保存为qtaim.cri：

molecule density.CHGCAR
units bohr
load density.CHGCAR smoothrho
reference 1
auto
cpreport cps.json
cpreport bonds.xyz graph

critic2 qtaim.cri > qtaim.cro
tail qtaim.cro
ls cps.json bonds.xyz

这里先明确加载真实VASP密度，并设置reference 1，随后AUTO搜索这一密度的非核临界点。没有使用默认的原子密度叠加来代替DFT密度。SMOOTHRHO用于网格电子密度的平滑插值，本版本采用默认NENV=512。

这条VTST＋Critic2直接读取VASP文件的路线已在本例最细网格实际执行；资料包采用ASE把同一完整密度转成CUBE的路线，两者临界点密度的相对差约6.7E−8。两条路线都保留原始电子数，没有把密度调成自己希望的值。关于输入参见Critic2标量场与分子结构。

4. 读键临界点与梯度路径#

python - <<'PYCODE'
import json
a = json.load(open('cps.json'))
print('坐标单位:', a['units'])
for cp in a['critical_points']['nonequivalent_cps']:
    if not cp['is_nucleus'] and cp['signature'] == -1:
        print('BCP密度:', cp['field'])
        print('Laplacian:', cp['laplacian'])
        print('Hessian:', cp['hessian_eigenvalues'])
PYCODE

本例找到两个核和一个H–H键临界点。该点Hessian有两个负本征值、一个正本征值，即(3,−1)鞍点；沿键方向与垂直键方向的曲率符号不同。

左图路径点来自Critic2实际导出；右图为x=7.5 Å平面，保留原网格采样，显示区域取原子附近。颜色是电子密度而不是轨道相位；对数色标下非正值仅不显示，下载数据保持原值。

本输入中CPREPORT JSON的坐标为Bohr，密度为electron/bohr³、Hessian和Laplacian为electron/bohr⁵。bonds.xyz的路径坐标按XYZ格式为Å，不能再按JSON单位缩放一次。包内脚本分别处理这两种长度约定。

5. 检查总电子数，再看网格敏感性#

在资料包根目录提取后的表格：

cat results/data/H2-QTAIM.csv

最细网格相对2电子的误差约1.28E−4，满足本例0.001电子的积分检查。粗120³的误差约0.001626，表中保留原数，没有分摊误差或强制归一化。

BCP密度在这些网格中很接近，Laplacian更敏感。图中另外保留完整周期网格的独立傅里叶导数：它与SMOOTHRHO是不同的插值描述，不要求人为对齐。比较自己的不同结构时应保持相同密度定义、网格与插值设置；对二阶导数，需要继续做针对性的精度检查。

临界点和梯度路径描述密度的拓扑特征。定量键能需要另按一致方法计算能量差；成键趋势也可以结合COHP和原子电荷、实际ELF例子与NCI弱相互作用分别讨论。

完整192³输入#

其余两份输入保存在包内inputs/，网格要求分别为96与144；实际输出大小以程序为准。

SYSTEM=H2 PAW AE density 192
ENCUT=550
PREC=Accurate
EDIFF=1E-6
NELM=120
ALGO=Normal
ISMEAR=0
SIGMA=0.05
ISPIN=1
ISYM=0
LREAL=.FALSE.
LASPH=.TRUE.
LWAVE=.TRUE.
LCHARG=.TRUE.
LAECHG=.TRUE.
LH5=.FALSE.
NBANDS=8
NSIM=1
NGXF=192
NGYF=192
NGZF=192
NSW=0
IBRION=-1
ISTART=1
ICHARG=0

H2 isolated singlet                     
   1.00000000000000     
    15.0000000000000000    0.0000000000000000    0.0000000000000000
     0.0000000000000000   15.0000000000000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   15.0000000000000000
   H 
     2
Direct
  0.4999999999999780  0.4999999999999807  0.4749967658946567
  0.5000000000000220  0.5000000000000193  0.5250032341053433

Gamma molecule
0
Gamma
1 1 1
0 0 0