Gaussian 带电分子的PCM与SMD实操

下载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. 先明确要比较什么#

本例使用OH⁻与H₃O⁺，分别在气相、IEFPCM水和SMD水中做B3LYP/6-31+G(d)优化与频率，再在实际气相优化几何上补算相同连续介质单点。OH⁻的电荷/多重度是 -1 1，H₃O⁺是 1 1，两者都是10电子。

先比较“固定同一气相核几何时换模型”，再比较“在同一溶剂模型中继续优化核几何”。这两步不能混成同一个量。对阴离子，基组中的弥散函数尤其需要注意；本例的 + 保留弥散函数，但小基组并非最终收敛结果。

2. 气相、PCM和SMD完整输入#

OH⁻气相：

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

Solvent and electric response oh-anion-gas-optfreq

-1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.000000000000 0.980000000000

OH⁻的IEFPCM水：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) Opt Freq SCF=Tight Int=UltraFine NoSymm SCRF=(IEFPCM,Solvent=Water)

Solvent and electric response oh-anion-pcm-water-optfreq

-1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.000000000000 0.980000000000

OH⁻的SMD水：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) Opt Freq SCF=Tight Int=UltraFine NoSymm SCRF=(SMD,Solvent=Water)

Solvent and electric response oh-anion-smd-water-optfreq

-1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.000000000000 0.980000000000

H₃O⁺气相：

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

Solvent and electric response hydronium-gas-optfreq

1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.960000000000 0.000000000000 0.200000000000
H -0.480000000000 0.831384000000 0.200000000000
H -0.480000000000 -0.831384000000 0.200000000000

H₃O⁺的IEFPCM水：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) Opt Freq SCF=Tight Int=UltraFine NoSymm SCRF=(IEFPCM,Solvent=Water)

Solvent and electric response hydronium-pcm-water-optfreq

1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.960000000000 0.000000000000 0.200000000000
H -0.480000000000 0.831384000000 0.200000000000
H -0.480000000000 -0.831384000000 0.200000000000

H₃O⁺的SMD水：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) Opt Freq SCF=Tight Int=UltraFine NoSymm SCRF=(SMD,Solvent=Water)

Solvent and electric response hydronium-smd-water-optfreq

1 1
O 0.000000000000 0.000000000000 0.000000000000
H 0.960000000000 0.000000000000 0.200000000000
H -0.480000000000 0.831384000000 0.200000000000
H -0.480000000000 -0.831384000000 0.200000000000

如需单独学习一个条件，复制对应输入到自己的新工作目录，在分配的CPU作业中运行：

g16 < input.gjf > output.log
grep 'Normal termination' output.log
grep 'Maximum Force' output.log
grep 'Frequencies --' output.log
formchk result.chk result.fchk

先确认整个作业正常结束、优化收敛，再检查频率。OH为双原子分子应有1个振动模式，H₃O⁺非线性四原子应有6个；本例相应模式均为正。

图中只为显示将O移到原点，不改变键长。连续介质没有显式水分子；需要研究水合壳层或具体氢键结构时，应另建显式溶剂模型。

3. 在实际气相几何上补溶剂单点#

以OH⁻的SMD为例，下列坐标来自上面已经完成的气相优化。这里不再写Opt/Freq：

%nprocshared=1
%mem=1GB
%chk=result.chk
#p B3LYP/6-31+G(d) SCF=Tight Int=UltraFine NoSymm SCRF=(SMD,Solvent=Water)

Derived solvent/charge reference oh-anion-smd-water-fixed

-1 1
O 0.000000000000 0.000000000000 0.002567304609
H 0.000000000000 0.000000000000 0.977432695632

流程：气相优化/频率→导出自己的优化坐标→在该坐标上算PCM/SMD单点→同模型优化/频率→按相同定义取能量。更换分子或重新优化后，必须更新对应单点坐标；不能继续使用本例坐标。

固定核几何不等于非平衡溶剂响应。本篇固定几何单点仍是基态平衡连续介质；激发态非平衡PCM是另一种响应过程。

4. 实际能量如何取数#

定义Egas为气相优化能量，Efixed为气相核几何上的溶剂模型单点，Erelaxed为溶剂中重新优化能量：

固定几何模型变化=Efixed−Egas；溶剂中核几何弛豫=Erelaxed−Efixed；二者相加为 Erelaxed−Egas。

cat analysis/continuum.csv
grep 'SCF Done' output.log
grep 'PCM non-electrostatic energy' output.log
grep 'Nuclear repulsion after' output.log
grep 'SMD-CDS' output.log

本例SMD日志把非静电项加入了核排斥能，随后SCF Done与FCHK Total Energy已经包含这项。不要再把SMD-CDS加到总能量上，否则会重复计算。附带脚本按实际核坐标重建裸核库仑能，再与非静电项、原生核排斥能逐项对照；不是只根据表头猜取数。

5. 这些结果能说明什么#

表中是本模型电子/模型总能差，未包含完整溶液化学势、标准态、完整构象与水合统计。带电物种的绝对溶剂化自由能还涉及参照约定，本表不能直接当实验单离子溶剂化自由能，更不能直接换成pKa或电极电势。

带电分子的偶极还依赖坐标原点；比较前必须约定原点。连续介质和固定核响应的区别继续见外场实操；脱电子的垂直/绝热定义见电子脱附。