Gaussian 原子电荷与布居分析实操

下载完整输入、实际数据和绘图脚本，解压到新目录：

bash run.sh

先从已有数据生成 analysis/、CSV和5张 figures/，不启动Gaussian。需要NumPy/Matplotlib，可使用集群ASE 3.29 Python环境，或通过PYTHON指定自己的解释器。

需要从头计算，先初始化个人Gaussian：

yskit gaussian-init

按Gaussian输入与CPU提交申请CPU资源，在分配的CPU作业里加载个人环境，再运行：

bash run.sh compute

顺序运行三份分子输入与formchk、九次Multiwfn电荷分析、三份cubegen静电势；首次还由Multiwfn调用个人Gaussian生成H/C/O三种ROHF原子参考。每步1线程、1GB。软件为个人授权Gaussian16 C.01 AVX厂商预编译CPU版本、内置NBO3.1，以及Multiwfn noGUI2026.10.18；没有公共Gaussian module或自编译器版本可套用。已有输出请另解压，完整work已有时可用 bash run.sh analyze。

公开的 native-fields.fchk.txt 是选定数值/基组/总密度字段，不能代替完整波函数供其他程序任意操作。需要分析自己的体系，请从自己的CHK生成完整FCHK。

1. 先选结构，再在同一密度上比较#

水使用固定几何，乙烯使用本组之前完成的基态优化结构，水二聚体固定O—O距离2.90Å。三者分别10、16、20电子，均为中性单重态。先优化自己的结构，再做这类分析；本组固定水二聚体用于演示电荷与成键，不是优化后的结合能结果。

水的完整输入：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) Pop=(Full,NBO,Hirshfeld,MK) SCF=Tight Int=UltraFine NoSymm

Gaussian population reference water

0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.758600000000 0.587100000000
H 0.000000000000 -0.758600000000 0.587100000000

Pop=Full打印更多布居信息；NBO调用本程序内置NBO，Hirshfeld输出Hirshfeld及CM5电荷，MK拟合分子静电势电荷。6-31+G(d)里的+为弥散函数，(d)为极化函数；不要把同名电荷方法在不同几何、泛函与基组上的结果当成同一条件。

2. 终端计算与提取#

在自己的新目录准备input.gjf，在分配的CPU作业内：

g16 < input.gjf > output.log
grep 'Normal termination' output.log
formchk result.chk result.fchk
sed -n '/Mulliken charges:/,/Sum of Mulliken charges/p' output.log
sed -n '/Hirshfeld charges, spin/,/Hirshfeld charges with hydrogens/p' output.log
sed -n '/Summary of Natural Population Analysis:/,/NATURAL BOND ORBITAL ANALYSIS:/p' output.log
sed -n '/ESP charges:/,/Sum of ESP charges/p' output.log

先确认正常结束，再核对原子编号、电荷之和与总电子数。本组Mulliken电荷还用真实AO重叠矩阵和总密度独立重建；不能把“有表格”直接当成正确读数。

3. 同一个氧原子为什么电荷不同#

水中O1的实际结果如下：

这些数值来自不同的电子密度划分或静电势拟合，不是六次计算对某个唯一实验原子电荷的六种测量。Mulliken对基组较敏感；NPA按自然原子轨道统计，Hirshfeld依赖自由原子参考密度，CM5是在相应电荷基础上的模型修正，MK关注拟合点上的静电势。

先确定用途：定性讨论电荷转移时，保持前后方法/基组/几何规则一致；用于力场点电荷时，还要选择与力场约定一致的拟合/约束流程。不能只挑电荷绝对值最大的方案，也不能把本例MK当作RESP拟合。

4. 电荷能否重建分子偶极#

Gaussian Hirshfeld同时打印原子偶极。本例水的z分量：电荷点贡献约0.346060a.u.，原子偶极之和约0.543292a.u.，加起来才得到总偶极0.889353a.u.。仅用ΣqARA会遗漏原子内部极化，不能据此宣布Hirshfeld输出错误。带电分子还要明确坐标原点。

完整逐原子表在 analysis/atomic-charges.csv；继续读NBO成键、静电势和Multiwfn参考密度。布居关键词背景见Gaussian作者手册存档，程序实际版本以自己的输出为准。