Gaussian 密度梯度与拉普拉斯
真实水分子密度导数、任意点与有限差分实操；区分完整密度、价电子网格和原子受力。
阅读本篇目录
本篇使用鸢算科技集群实际 Gaussian 16 C.01 AVX CPU 版和水分子 B3LYP/6-31G(d)、UltraFine 单点结果。Gaussian 为厂商预编译版本；配套分析使用 ASE 3.29 环境中的 NumPy 和 Matplotlib。这里的电子密度导数不是原子所受的力，也不是分子动力学中的速度。
下载输入、真实数值网格、终端脚本和高清绘图程序，解压进入目录：
bash run.sh
先生成 analysis/ 与 figures/ 查看已有真实结果，无需Gaussian。要从自己的检查点计算，先 yskit gaussian-init，再在已分配CPU的作业中执行：
bash run.sh compute
流程依次是水分子单点、formchk、网格和向量运算、数值分析与三张图。每个输入1线程、1GB内存；不要在客户登录终端启动正式计算。已有输出时请另解压一份，脚本拒绝覆盖；提交方法见Gaussian输入与检查点。
1. 先得到自己的完整FCHK#
本例完整输入为：
%chk=result.chk
%NProcShared=1
%Mem=1GB
#p B3LYP/6-31G(d) SCF=Tight Int=UltraFine NoSymm
water-b3lyp-ultrafine
0 1
O 0.00000000 0.00000000 0.00000000
H 0.00000000 0.75700000 0.58700000
H 0.00000000 -0.75700000 0.58700000
在已分配的CPU作业中运行：
g16 < input.gjf > output.log
grep 'Normal termination' output.log
formchk result.chk water.fchk
本例固定几何，不是优化后的水结构；10电子、闭壳层，无ECP。二进制检查点和完整波函数只留在自己的目录，下载包只包含数值网格。
2. 密度梯度、梯度模与拉普拉斯#
cubegen 1 Gradient=SCF water.fchk gradient.cube -2 h
cubegen 1 NormGradient=SCF water.fchk normgradient.cube -2 h
cubegen 1 Laplacian=SCF water.fchk laplacian.cube -2 h
cubegen 1 FDensity=SCF water.fchk full-density.cube -2 h
cubegen 1 Density=SCF water.fchk valence-density.cube -2 h
| 输出 | 每个格点的值 | 量纲 |
|---|---|---|
| Gradient | 密度ρ、∂ρ/∂x、∂ρ/∂y、∂ρ/∂z，共4项 | ρ为e/Bohr³，导数为e/Bohr⁴ |
| NormGradient | √[(∂ρ/∂x)²+(∂ρ/∂y)²+(∂ρ/∂z)²] | e/Bohr⁴ |
| Laplacian | ∂²ρ/∂x²+∂²ρ/∂y²+∂²ρ/∂z² | e/Bohr⁵ |
先看cube头第三行的最后一项 NVal。本例 Gradient 是4，其余两项为1；不能用只支持一个数值的解析程序直接读取4分量网格。
**本次C.01实际导出的 Gradient=SCF 密度分量与 FDensity=SCF 的完整电子密度一致。**它不同于 Density=SCF 默认省去冻结芯层的价电子网格。不要因为参数里都有SCF，就默认这些分支具有同一冻结核选择；先与自己的完整/价电子密度作数值对照。
本例粗网格均为31×31×31，cube坐标与步长为Bohr，步长约0.333333 Bohr。完整密度体素积分约13.597412，价电子约8.170642。完整水分子的电子数仍是10；前者严重偏离来自粗网格采样芯层尖峰，不能把这个值当作电荷转移或新增电子数。定量积分要按密度网格教程分别收敛范围与步长。
白色箭头使用实际梯度的单位方向，显示密度增大的方向，箭头长度不代表梯度大小；不是电子运动轨迹或原子受力。图取最接近分子平面的原生格点，x约0.01548 Å；拉普拉斯色标截在±5，原始数值完整保留。正负拉普拉斯可以帮助描述密度局部变化，不能单凭一张切片确定全部化学键性质或临界点。
3. 在指定空间点计算#
cubegen -5 可读取任意点。本次C.01该输入与回显坐标为Å，输出导数仍采用原子单位。
cat > points.in <<'EOF'
0.45 0.30 0.35
-0.65 0.25 0.40
0.25 -0.35 1.20
0.20 1.10 0.60
EOF
cubegen 1 Gradient=SCF water.fchk point-gradient.txt -5 h < points.in
cubegen 1 Laplacian=SCF water.fchk point-laplacian.txt -5 h < points.in
保留 h 时输出含cube式头信息；头后每行依次是x、y、z和对应数值，Gradient 后有ρ与三个导数。本例四个离核点的实际结果：
| 点 | ρ / e·Bohr⁻³ | ∂ρ/∂x | ∂ρ/∂y | ∂ρ/∂z | ∇²ρ / e·Bohr⁻⁵ |
|---|---|---|---|---|---|
| 1 | 0.255917615 | -0.439653649 | -0.231511901 | -0.302823204 | 0.364361595 |
| 2 | 0.117302642 | 0.252968633 | -0.079410550 | -0.146482228 | 0.394730561 |
| 3 | 0.019729951 | -0.013361047 | -0.006402128 | -0.045096863 | 0.080730489 |
| 4 | 0.067752437 | -0.078241339 | -0.160354184 | -0.022438891 | 0.075619532 |
4. 如何知道这些导数与密度一致#
配套脚本在每个离核点，分别取中心及x/y/z方向的正负位移，从独立的 FDensity=SCF 数值构造中心差分，不直接把解析导数再抄一遍。步长h先0.02 Bohr，再减半为0.01 Bohr；生成输入点时乘0.529177210903转换成Å。
- 一阶导数：
[ρ(x+h)-ρ(x-h)]/(2h)。 - 拉普拉斯：三个方向的
[ρ(x+h)-2ρ(x)+ρ(x-h)]/h²相加。
第一点的实际梯度向量误差从0.000114183降至0.0000285431 e/Bohr⁴，拉普拉斯误差从0.000375600降至0.0000939352 e/Bohr⁵；四个点的误差都随步长减半下降约4倍，符合中心差分的二阶趋势。
bash run.sh
cat analysis/finite-difference.csv
这里展示离核点的数值一致性，不能代替近核空间网格收敛、QTAIM临界点搜索或完整成键分析。需要处理分量、叉乘或单位向量时继续看向量网格与分量选择。
工具参考：Gaussian编写的cubegen手册存档为较早版本说明，本文命令和数据以本次实际C.01行为为准。