Gaussian ONIOM边界与受力实操

下载完整输入、选定原生数据与绘图脚本，解压到新目录：

unzip gaussian-oniom-practical*.zip -d oniom-example
cd oniom-example
bash run.sh

这会用随附的实际计算结果生成 analysis/energies.csv、analysis/finite-difference.csv 和 figures/ 下的5张图。脚本需要NumPy和Matplotlib，集群默认使用ASE的Python环境；其他电脑可用 PYTHON=python3 bash run.sh 指定自己的环境。

本组采用个人有权使用的 Gaussian16 C.01 AVX CPU程序，单例1核心、1GB；资料包不包含Gaussian程序。完成自己的Gaussian环境初始化、申请CPU资源后检查入口：

command -v g16
command -v formchk

需要全部重算时，在新解压目录执行 bash run.sh compute，顺序运行23个很小的计算；已有 work/ 同名目录时会拒绝覆盖。该命令不负责申请队列资源，应在已分配的CPU任务中执行。

1. link H的位置由哪两个真实原子决定#

本例link H位于O3→C2连线上：Rcap=RO3+g(RC2-RO3)。从Gaussian实际生成的模型坐标得到 g≈0.686808847。这个数值来自本例的类型和键，不是通用固定常数。

InputFiles 选项会在日志中打印各个分量的输入，便于核对link H和模型的连接。模型有3个活动原子：link H、O3、H4；打印诊断里其他完整体系中心可能作为非活动占位出现，不能据打印行数误判量子区的原子数。

bash compute.sh ethanol-hf-uff-force
grep -A 15 'Point  2 -- high level on model system' work/ethanol-hf-uff-force/output.log

2. 显式修改缩放因子#

下面实际计算采用两个相同的显式因子0.65，分别用于模型低层和模型高层；资料包还有0.72的对照。

%nprocshared=1
%mem=1GB
%chk=result.chk
#p ONIOM(HF/STO-3G:UFF) Force Geom=Connectivity NoSymm SCF=Tight

Small layered-model link-scale-65

0 1 0 1 0 1
C-C_3-0.00000000 0 1.171592164149 -0.403050689766 -0.000000000006 L
C-C_3-0.00000000 0 -0.006726532503 0.556916685590 -0.000000000006 L H-H_-0.0 3 0.650000000 0.650000000
O-O_3-0.00000000 0 -1.201076743568 -0.221755587510 -0.000000000001 H
H-H_-0.00000000 0 -1.953062098208 0.389286604240 0.000000000014 H
H-H_-0.00000000 0 0.039662126550 1.209262019623 0.887716250593 L
H-H_-0.00000000 0 0.039662126550 1.209262019623 -0.887716250593 L
H-H_-0.00000000 0 2.118986840378 0.147543167958 0.000000000001 L
H-H_-0.00000000 0 1.140320061155 -1.044590109792 0.886956786043 L
H-H_-0.00000000 0 1.140320061155 -1.044590109792 -0.886956786043 L

1 2 1.0 7 1.0 8 1.0 9 1.0
2 3 1.0 5 1.0 6 1.0
3 4 1.0
4
5
6
7
8
9

bash compute.sh link-scale-65
bash compute.sh link-scale-72
grep 'ONIOM: extrapolated energy' work/link-scale-65/output.log

改变link H距离会改变模型能量和边界受力。这里的0.65、0.72用于说明敏感性，不是通过试出最低能量来替代合理边界的选择。高低两层的link atom定义应保持一致。

3. 边界受力不是只把模型力贴到高层原子#

模型link H的力需要按坐标关系传给O3和C2：模型高层减模型低层的cap力中，g部分传给C2、1-g部分传给O3，随后加上完整体系低层受力。资料包用三份独立计算重建完整受力，最大差约 6×10⁻¹¹ Hartree/Bohr。

还对C2的x、H4的y和H7的x分别做了正负位移，每个方向采用0.001Å和0.0005Å两个步长，共12个有限差分输入。采用 F≈-ΔE/(2h)，将Å换算为Bohr后比较。

4. 提取数据和画图#

head -n 7 analysis/finite-difference.csv
ls figures/03-link-forces-finite-difference.png
bash run.sh

默认 bash run.sh 用随附真实数据提取CSV和画图；全部重算采用 bash run.sh compute。只算一个位移输入时可用 bash compute.sh fd-2xp500，在自己的新目录内与负位移一起读能量，不把正位移或不同模型的能量单独当作力。

link atom和缩放顺序见Gaussian ONIOM手册。