Gaussian 偶极与极化率张量实操

下载25份完整输入、原生数值与绘图脚本，解压到新目录后：

bash run.sh

默认读取已有数据，生成 analysis/ 和5张 figures/，不启动Gaussian。NumPy、Matplotlib来自集群ASE 3.29 Python环境；也可通过PYTHON指定自己的解释器。

需要重算时，先初始化个人Gaussian环境：

yskit gaussian-init

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

bash run.sh compute

程序依次完成25份Gaussian输入、逐个运行formchk，再从自己的输出分析与画图。本例是Gaussian16 C.01 AVX厂商预编译CPU版本，个人独立初始化，没有公共Gaussian module或自编译器/MPI版本可以套用。每份输入1线程、1GB，按实际分配资源调整；不是GPU或Intel Xeon 8173M性能测试。

已有输出时另解压一份，避免覆盖。已有完整work而尚无analysis/figures时可用 bash run.sh analyze。公开 native-fields.fchk.txt 只是数值摘录，不是完整波函数；自己的CHK须先用formchk转换。

1. Polar先算什么#

极化率描述小外场下偶极的响应。对于孤立中性分子，先在固定核几何上算电子极化率；这不同于有场优化后的总响应，也不等于周期材料的介电常数。

本例零场完整输入：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) SCF=VeryTight Int=UltraFine NoSymm Polar

Solvent and electric response water-zero-polar

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

Polar请求解析响应，SCF=VeryTight有助于与很小场的有限差分进行比较。这里固定的是给定水结构，不把其电子响应称为弛豫离子介电响应。

2. 正负场与零场一起准备#

X方向一对完整输入：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) SCF=VeryTight Int=UltraFine NoSymm Field=X+10

Solvent and electric response water-field-xp10

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

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) SCF=VeryTight Int=UltraFine NoSymm Field=X-10

Solvent and electric response water-field-xm10

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

对应Y/Z、±20输入均在下载包。不要只算一个正场就把偶极除以场强：永久偶极需要消去，且要保持核坐标完全相同。

按外场符号约定，q为Gaussian算符系数，本例F=−q：

αᵢⱼ≈−[μᵢ(qⱼ=+h)−μᵢ(qⱼ=−h)]/(2h)。

沿X施场时同时取μX、μY、μZ，得到张量的一整列；Y、Z同理。只取三个平行分量只能得到对角项，不能声称已经核对完整张量。

3. 终端提取与绘图#

在分配的CPU作业中：

g16 < input.gjf > output.log
formchk result.chk result.fchk
grep -A2 'Polarizability ' result.fchk
grep -A1 'Dipole Moment' result.fchk

本例FCHK的6个解析分量按下三角顺序 xx,xy,yy,xz,yz,zz恢复对称3×3张量；不要按3行3列直接读取原始6个数。每个分量单位为原子单位。

bash run.sh
cat analysis/RESULTS.json

4. 解析结果与两种步长对照#

本例平均电子极化率 Tr(α)/3=6.80140959 a.u.。非对角项在本坐标下接近零，仍从实际三个方向响应中保留。

脚本另外用 −[E(+h)+E(−h)−2E(0)]/h²核对相应对角极化率。这种二阶能量差比单点能量更敏感，步长太小可能被SCF/打印误差淹没，太大则混入更高阶响应。两套步长只是本例的局部数值比较，不是任意体系已收敛的保证。

如将极化率写成ų，指极化率体积 αSI/(4πε₀)；本例原子单位数乘以a₀³即可，不能把它当作盒子体积或直接推出固体介电常数。带电物种绝对偶极的原点依赖也不能省略说明。

流程：统一核坐标→零场Polar→三个方向正负小场→逐方向取三个偶极分量→中心差分恢复张量→能量二阶差分与不同步长对照→再选择适合自己体系的精度。