Gaussian 后HF能量与梯度实操

下载23个条件的完整输入、原生数值摘录和绘图脚本，在新目录解压后先查看已有结果：

bash run.sh

得到 analysis/ 与 figures/。这一步不调用Gaussian，使用现有ASE 3.29 Python环境中的NumPy、Matplotlib分析和绘图。

需要重算时先准备自己的Gaussian环境：

yskit gaussian-init

按Gaussian输入与提交的CPU作业流程准备脚本，将作业执行命令改为 bash run.sh compute。已经获得CPU作业资源时，可在另一份新解压目录直接执行：

bash run.sh compute

程序依次运行23份原生输入、23次formchk，再从自己的 work/ 分析和画图。不要在客户登录终端启动正式计算。已有 work/analysis/figures 时请另解压一份；只有 work 而没有分析/图目录时可用 bash run.sh analyze 继续后处理。

软件为Gaussian 16 C.01 AVX厂商预编译CPU版本，按用户独立初始化，没有公共Gaussian module；本例没有自编译器/MPI配置可供套用。每份参考输入1线程、1GB，采用整节点模板时按实际分配与运行安排统一设置线程/内存。本例不作为GPU或Intel Xeon 8173M节点的性能数据。

1. 先固定几何、基组与电子状态#

本例使用中性单重态水分子，电荷/多重度为 0 1，10电子。两个O–H输入距离均0.96 Å，H–O–H角104.5°；HF、MP2、CCSD、CISD和CCSD(T)都在同一固定几何及6-31G(d)下比较，没有把几何优化差异混进方法比较。

2. MP2能量和梯度的完整输入#

inputs/mp2-force-0/input.gjf：

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

Water correlation/composite practical mp2-force-0

0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.759061990794 0.587728588833
H 0.000000000000 -0.759061990794 0.587728588833

Force在固定几何上计算梯度，不执行核坐标优化；SCF=Tight控制参考SCF，NoSymm方便保持本例坐标方向。%chk写入自己的检查点，%nprocshared与%mem分别指定线程和内存。

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

cd inputs/mp2-force-0
mkdir -p scratch
export GAUSS_SCRDIR="$PWD/scratch"
g16 < input.gjf > output.log
formchk result.chk result.fchk

这是对一个新输入目录的操作；整套23条件更方便的入口仍是包根目录的 bash run.sh compute。

本例MP2原生日志给出 NFC=1，冻结一个氧1s占据轨道，相关步骤涉及4个α与4个β占据轨道。总电子数仍是10，不是把氧芯电子从分子电荷中删掉。冻结核设置与全电子相关计算需要分别说明。

3. 读取正确的方法能量#

grep 'SCF Done' output.log
grep 'EUMP2' output.log
grep 'Total Energy' result.fchk

MP2日志里的 SCF Done 是参考HF能量；MP2总能量看 EUMP2 或相应FCHK能量字段。本例还将三个自旋配对相关能分项相加，核对 E_MP2 = E_HF + E2。

以下使用显式 Conver=10 的CCSD、CISD与CCSD(T)，其参数见相关能收敛实操：

这里的19个基函数适合快速学习方法和输入。不同近似的绝对能量不能直接当作准确度排行榜；研究反应能或物性时还应检查基组、参考波函数和目标能量差。

4. CCSD(T)单点的完整输入#

inputs/ccsdt-singlepoint/input.gjf：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p CCSD(T,Conver=10)/6-31G(d) SCF=Tight NoSymm

Water correlation/composite practical ccsdt-singlepoint

0 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.759061990794 0.587728588833
H 0.000000000000 -0.759061990794 0.587728588833

这次只做单点，最终能量包含微扰三激发修正。需要取数时看独立的 CCSD(T)=高精度行或FCHK的 Cluster Energy with triples/最终总能量字段，不要把参考HF、MP2或CCSD能量当成CCSD(T)。

FCHK里出现 Cartesian Gradient 并不自动证明它是当前最终方法的梯度；多步骤作业可能保留其他步骤字段。本篇只对明确请求 Force 的HF、MP2、CCSD、CISD展示实际梯度与差分，CCSD(T)保留为单点示例。正式优化时按所选方法支持的导数与资源需求选择作业类型。

5. 从能量检查一个原子的受力#

包内对H1的y坐标独立计算 y−0.0001 Å、原位置、y+0.0001 Å 三个条件：

F_y ≈ −[E(y+h)−E(y−h)]/(2h)。

Gaussian原生Cartesian Gradient单位为Hartree/Bohr，受力等于负梯度；换成eV/Å后，MP2本例为约0.33617210 eV/Å，能量差分约0.33617241 eV/Å。两者方向和单位一致。

在包根目录查看全部方法：

head analysis/energies.csv
head analysis/forces.csv

reference/*/output.log和 native-fields.fchk.txt是实际原生数值行/字段摘录，省略个人路径、程序正文和轨道系数；它们不是完整波函数文件。自己的重算结果在 work/，提供完整CHK/FCHK供后续使用。

继续看相关能与力的收敛以及G4/CBS复合方法。