Gaussian Counterpoise与幽灵基组实操

下载完整输入、原生数值摘录与高清图脚本，在新目录解压后先复画已有数据：

bash run.sh

得到 analysis/ 与5张 figures/，此步骤不启动Gaussian，使用集群ASE 3.29 Python环境中的NumPy、Matplotlib。

需要重算时，先在个人账号准备Gaussian：

yskit gaussian-init

按Gaussian输入与CPU提交申请CPU作业，将执行命令设置为 bash run.sh compute。已取得CPU资源并加载自己的Gaussian环境后，也可在另一份新解压目录执行：

bash run.sh compute

程序顺序完成12个HF/BSSE输入与3套完整G4任务，运行各自的formchk，再从自己的 work/ 分析和画图。水的G4数据复用此前完整实算；本包保留来源和参考输入，不重复计算这份水。客户正式任务应在分配到的CPU资源内运行。

已有 work/analysis/figures 时另解压一份；只已有完整 work/ 而无分析和图片时，可用 bash run.sh analyze继续后处理。Gaussian为16 C.01 AVX厂商预编译CPU程序，按用户独立初始化，没有公共Gaussian module或自编译器/MPI版本需要套用。参考输入每份1线程、1GB，按实际CPU资源调整，不作为GPU或Intel Xeon 8173M节点性能数据。

1. 先区分相互作用能与结合自由能#

本例在固定水二聚体几何上计算HF/6-31G(d)电子相互作用能，分别使用O–O=2.90 Å和5.00 Å。片段始终保持复合物里的几何，不单独优化片段，也没有振动/熵修正；因此结果不是结合自由能，也不是加上形变能后的完全弛豫结合能。

有限基组中，一个片段可以“借用”另一个片段的基函数，造成基组重叠误差（BSSE）。Counterpoise用幽灵中心保留基函数、去掉另一片段的核电荷和电子，从而统一片段和复合物使用的基组空间。

2. 最直接的原生Counterpoise输入#

inputs/dimer-2p9-counterpoise/input.gjf完整内容：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p HF/6-31G(d) SCF=Tight NoSymm Counterpoise=2

Reaction energy practical dimer-2p9-counterpoise

0 1 0 1 0 1
O(Fragment=1) 0.000000000000 0.000000000000 0.000000000000
H(Fragment=1) 0.000000000000 0.000000000000 0.960000000000
H(Fragment=1) 0.000000000000 0.929400000000 -0.240400000000
O(Fragment=2) 0.000000000000 0.000000000000 2.900000000000
H(Fragment=2) 0.000000000000 0.759060000000 3.487730000000
H(Fragment=2) 0.000000000000 -0.759060000000 3.487730000000

Counterpoise=2声明两个片段，Fragment=1/2标记归属。电荷行 0 1 0 1 0 1依次是总体系、片段1、片段2的电荷和多重度；三个对象在这里均为中性单重态。换体系时要按自己的电子数和自旋修改，不能只改原子名字。

在已分配CPU资源、加载个人Gaussian环境后，可单独运行这份输入：

cd inputs/dimer-2p9-counterpoise
mkdir -p scratch
export GAUSS_SCRDIR="$PWD/scratch"
g16 < input.gjf > output.log
formchk result.chk result.fchk
grep -E 'Counterpoise corrected|BSSE energy|sum of fragments' output.log

整套15输入更方便的入口是包根目录 bash run.sh compute。

3. 独立计算每个组成能量#

下载包除了原生Counterpoise，还在每个距离独立运行复合物AB、孤立片段A/B、带另一片段幽灵基函数的A/B，共5个单点。以“A是真实原子，B为幽灵中心”的完整输入为例：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p HF/6-31G(d) SCF=Tight NoSymm

Reaction energy practical dimer-2p9-ghostA

0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.000000000000 0.960000000000
H 0.000000000000 0.929400000000 -0.240400000000
O-Bq 0.000000000000 0.000000000000 2.900000000000
H-Bq 0.000000000000 0.759060000000 3.487730000000
H-Bq 0.000000000000 -0.759060000000 3.487730000000

O-Bq、H-Bq保留对应原子的基函数，但该中心核电荷为0。此作业只有真实水分子的10电子，使用整个二聚体的38个基函数；孤立水是19个基函数。不要把幽灵中心当成带氧/氢核的额外真实原子。

孤立A的完整输入：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p HF/6-31G(d) SCF=Tight NoSymm

Reaction energy practical dimer-2p9-isolatedA

0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.000000000000 0.960000000000
H 0.000000000000 0.929400000000 -0.240400000000

另一份ghostB保留B的真实核、将A改为幽灵中心，完整文件在 inputs/dimer-2p9-ghostB/input.gjf。复合物和5 Å对照也都提供完整输入。

4. 按同一几何做能量差#

记 E_A^A为孤立A基组下的片段能量，E_A^AB为带B幽灵基函数的A能量：

ΔE_raw = E_AB^AB − E_A^A − E_B^B

ΔE_CP = E_AB^AB − E_A^AB − E_B^AB

BSSE = E_A^A + E_B^B − E_A^AB − E_B^AB

所以 ΔE_CP = ΔE_raw + BSSE。本例BSSE为正，修正后相互作用的绝对值减小。

2.90 Å的实际组成能量，单位Hartree：

5. 两个距离的实际结果#

这些是同一组固定方向、有限基组几何的结果，不是水二聚体全局最优构型或完整势能曲线。5 Å的相互作用与BSSE都更小，但不等于严格无相互作用。

6. 读对原生字段#

本例日志 Counterpoise corrected energy是修正后的复合物总能量，还需减去 sum of fragments才能得到相互作用能。BSSE energy本身也不是结合能。

原生Counterpoise包含5次SCF，最后一条 SCF Done可能属于片段，不能直接当复合物能量。本例最终FCHK的 Total Energy则等于Counterpoise修正总能量，而其 SCF Energy仍是最后一个片段步骤；字段要结合实际作业类型解释，不能机械地选最后一条能量。

grep -E 'Total Energy|SCF Energy' result.fchk
head analysis/counterpoise.csv

最后一条在包根目录分析完成后执行。流程为：定义片段与固定几何→准备Counterpoise输入→提交CPU作业→读取三行摘要→按相同参照求差→需要时增加优化、频率与热修正。反应量的系数和参照处理见反应能与自由能台阶。