Gaussian 16 入门：输入文件的语法与算例

本讲导读

上一讲我们用差分法手写了一个求解一维 Schrödinger 方程的小程序。 从这一讲开始，我们换用真实科研中正在使用的工具——Gaussian 16。 它是量子化学领域使用最广的程序之一，而它的使用方式非常有代表性： 写一个纯文本输入文件，交上去，然后读一个纯文本输出文件。

这一讲要解决的就是这三个动作：

• 怎么写——一份 Gaussian 输入由哪几段组成，每一段的语法是什么；
• 关键词是什么意思——Opt、Freq、SCF、Guess 分别对应哪个近似、什么场景下要改；
• 怎么读——日志里哪几行才是关键，怎么判断作业是成功还是失败。

整讲围绕六个算例展开（H₂O₂ → C₂H₄ → CH₃F → NH₃ → 呋喃 → 乙炔）， 它们的技术主线只有一条：用最少的变量描述一个分子。

一、Gaussian 能做什么对应讲义 PPT 第 2–3 页 + 插页 1–5

1.1 能力清单

Gaussian 的官方能力页列了很长一串，但可以归成四类。值得注意的是： 这一整张清单背后只有一件事——在给定近似下求解薛定谔方程。

官网：https://gaussian.com 　·　能力总览页：gaussian.com/capabilities

1.2 内部结构：Link 与 Overlay

Gaussian 不是一个单体程序，而是由许多链接（Link）串起来的流水线。 每个 Link 负责一段工作，前一个的输出是后一个的输入。

1.3 软件的起源：1970 年

1970 年，约翰·波普尔（John A. Pople）及其在卡内基梅隆大学的研究团队推出 Gaussian 70，并通过量子化学程序交换组织（Quantum Chemistry Program Exchange，QCPE）向学术界提供。这一程序使其他研究者能够利用已有的计算工具开展分子电子结构研究，促进了量子化学计算方法的传播与应用。

Gaussian 70 的意义在于：它把量子化学理论、数值算法与计算机程序打包成一套可以直接使用的工具——研究者不必再各自重写一遍从头算程序，可以把精力留给分子本身。

1.4 为什么叫「Gaussian」：高斯型基函数

软件的名字来自它采用的高斯型基函数（Gaussian-type function）。在分子轨道计算中，分子轨道通常表示为有限个基函数的线性组合——高斯型基函数的典型形式是：

式中坐标以基函数中心为原点：N 是归一化常数，α 控制函数的空间分布（α 越大越紧凑），l、m、n 决定其角向特征。

与 Slater 型函数相比，高斯型函数的指数上是 r² 而不是 r，因此多中心积分可以解析求出，计算效率高得多；代价是单个高斯函数描述原子轨道的形状并不好——它在原子核附近的行为与真实轨道不符。

解决的办法是把若干高斯函数按固定系数组合成一个收缩基函数，既保留积分的便利，又改善对原子轨道的描述，从而在计算精度与计算成本之间取得平衡。这条技术路线为早期计算资源有限条件下的分子电子结构计算提供了重要基础。

1.5 版本演进与功能扩展

Gaussian 自发布以来持续迭代，主要版本以发布年份命名：

随着理论方法和数值算法的进步，Gaussian 的功能逐渐扩展到分子几何结构优化、振动频率与光谱计算、过渡态搜索、反应路径分析以及激发态性质研究等方面。电子相关方法和密度泛函理论的引入与发展，也进一步扩大了其应用范围。

1.6 密度泛函理论的加入

密度泛函理论（DFT）是现代电子结构计算的重要理论基础之一。沃尔特·科恩（Walter Kohn）及其合作者为该理论的建立作出了奠基性贡献。DFT 进入 Gaussian 之后，原本只能在小体系上做的计算得以在更大的分子上展开——这也是今天绝大多数 Gaussian 作业用 DFT 泛函（如 B3LYP）的原因。

1.7 模型化学与 1998 年诺贝尔化学奖

波普尔的贡献不仅在于软件开发，还在于推动「模型化学」（model chemistry）思想：将明确的理论近似与基组相结合，通过系统计算和实验比较检验其可靠性，进而用于预测分子的结构、性质及反应行为。

这套方法包含四个环节：① 选定近似——明确的理论方法 + 基组，写进 Route 段；② 系统计算——对同一系列分子做同样规格的计算；③ 与实验比较——检验这套「模型」到底可不可靠；④ 用于预测——预测分子的结构、性质与反应行为。这种研究思路使量子化学计算逐步形成了可检验、可重复的研究方法体系。

二、输入文件的五段结构对应讲义 PPT 第 4–16 页

2.1 总览与一个完整例子

一份 Gaussian 输入文件由若干「段落（section）」组成，段落之间用空行分隔， 顺序不能颠倒：

1. Link 0 段（以 % 开头）——指定暂存文件与资源，如 %chk、%mem、%nprocs；
2. Route 段（以 # 开头）——指定方法、基组、任务类型以及各种关键词；
3. Title 段——一行任意说明文字，必须有，但不参与计算；
4. 分子规格段——第一行是电荷与多重度，其后是坐标；
5. （可选）附加段——自定义基组、变量初值等，紧跟在坐标之后。

以水分子为例，把五段结构对应上：

%chk=water.chk          <- Link 0 段：以 % 开头
%rwf=water.rwf
%nprocs=1

#p HF/6-31g             <- Route 段：以 # 开头

Water                   <- Title 段：任意文字

0 1                     <- 电荷=0，多重度=1（单线态）
O                       <- 坐标段
H  1  R1
H  1  R1  2  a1

R1=1.04                 <- 变量段（内坐标写法才需要）
a1=104.5

怎么读这份文件：从上往下逐段问自己—— ① 结果要存到哪（%chk）？ ② 用什么方法、算什么（HF/6-31g，默认单点能）？ ③ 算的是谁（水分子）？ ④ 体系带什么电、什么自旋（0 1）？ ⑤ 原子在哪里（坐标或内坐标）？ 这五问能答全，输入文件就不会写错。

2.2 % 段与临时文件

这些文件的名字可以任意取，但后缀决定了它们的用途：

2.3 内存估算

内存不够是 Gaussian 失败最常见的原因。经验估算式为

2.4 Route 段与语法规则

Route 段的第一行以 # 开头，紧跟一个字母控制输出详细程度：

Route 段本身是自由格式的，规则如下：

实际写法建议：同一行的关键词用空格分隔；同一关键词的多个选项用逗号 包在括号里。这样最易读，也不容易触发 80 字符上限。

2.5 Title 段与电荷、多重度

2.6 分子规格段：三种坐标写法

H                       <- 第一个原子只需元素符号
O  1  R1                <- 与原子 1 成键，距离 R1
O  2  R2  1  A          <- 再指定与原子 1 的键角 A
H  3  R1  2  A  1  D    <- 再指定二面角 D

R1=0.9                  <- 变量在坐标之后集中定义
R2=1.4
A=105.0
D=120.0

三、常用关键词详解对应讲义 PPT 第 8–13 页

以下五个关键词覆盖了本课程绝大多数作业的写法。每个都按「选项 → 含义 → 什么时候改」 三层来读。

3.1 Guess：初始猜测

SCF 迭代需要一个初始的密度（或轨道），Guess 决定它从哪来：

3.2 SCF：收敛控制

3.3 Opt：几何优化

3.4 Freq：频率计算

3.5 热力学修正：从电子能到自由能

freq 作业在算完频率后会自动做热力学分析。总能量由四部分构成：

Zero-point correction=                           0.022502
Thermal correction to Energy=                    0.023452
Thermal correction to Enthalpy=                  0.023768
Thermal correction to Gibbs Free Energy=         0.017992
Sum of electronic and zero-point Energies=     -76.408954
Sum of electronic and thermal Free Energies=  -76.413464

四、算例 1：H₂O₂ 的几何优化对应讲义 PPT 第 17–23 页

4.1 两种坐标写法，同一个分子

过氧化氢 H₂O₂ 有 4 个原子，用内坐标只需要 6 个变量 （3 个键长 + 2 个键角 + 1 个二面角），比笛卡尔坐标的 12 个少一半。

完整的输入文件如下：

%chk=h2o2.chk
%rwf=h2o2.rwf
#p hf/6-31g opt

H2O2 energy calculation

0 1
H
O 1 0.9
O 2 1.4 1 105.0
H 3 0.9 2 105.0 1 120.0

4.2 逐项读日志

优化跑完会得到一个 .log 文件。以下四处是必须会看的。

（一）SCF 能收敛了吗——用 grep 抽出所有 SCF Done 行：

$ grep "SCF Done" h2o2.log

 SCF Done:  E(RHF) =  -150.710007216     A.U. after    9 cycles
 SCF Done:  E(RHF) =  -150.710007332     A.U. after    9 cycles
 SCF Done:  E(RHF) =  -150.710007337     A.U. after    8 cycles
 SCF Done:  E(RHF) =  -150.710007340     A.U. after    8 cycles

每一行对应优化过程中的一步几何。可以看到能量逐次下降并趋于收敛 （−150.710007216 → −150.710007340），迭代圈数稳定在 8–9 圈。 如果圈数在慢慢往上爬（9 → 15 → 30…），说明当前几何附近 SCF 越来越难收敛， 应当考虑加 SCF=QC。单位：A.U. 即 Hartree， 1 Hartree = 27.211 eV。

（二）优化前的初始参数：

                          !   Initial Parameters    !
                          ! (Angstroms and Degrees)  !
 ! Name  Definition              Value          Derivative Info.   !
 ! R1    R(1,2)                  0.9000         -DE/DX =   -0.0246   !
 ! R2    R(2,3)                  1.4000         -DE/DX =   -0.0312   !
 ! A1    A(1,2,3)              105.0000         -DE/DX =    0.0018   !
 ! D1    D(1,2,3,4)            120.0000         -DE/DX =    0.0000   !

Name 是内坐标变量名，Definition 说明它由哪几个原子定义 （R(1,2) = 原子 1 与 2 的键长），-DE/DX 是能量对该坐标的负梯度—— 也就是力。力为零的方向就是优化的方向。优化的目标很直白： 把所有 −DE/DX 都压到接近 0。

（三）收敛判据与最终结构：

 -- Stationary point found.

 ! R1    R(1,2)                  0.9661         -DE/DX =    0.0000   !
 ! R2    R(2,3)                  1.4028         -DE/DX =    0.0000   !
 ! A1    A(1,2,3)              100.1574         -DE/DX =   -0.0001   !
 ! D1    D(1,2,3,4)            121.2654         -DE/DX =    0.0000   !

「Stationary point found」是收敛的正式标志。 对照初始与优化后的数值可以学到物理：O—O 键从 1.400 Å 变到 1.4028 Å（几乎没变）， 而 O—H 键从 0.900 Å 明显拉长到 0.9661 Å，键角从 105° 收到 100.2°—— HF 方法能正确给出过氧键较长的特征。这也说明：初始值离真实值太远时， 优化的第一步往往很大，所以给一个合理的初始几何能省不少时间。

（四）Archive 段与终止信息：

 (Enter /export/home/gaussian/g09/l9999.exe)
 1\1\GINC-DIRAC2\FOpt\RHF\6-31G\H2O2\YLNIU\15-Oct-2018\0\\#p hf/6-31g opt\\...
 Normal termination of Gaussian 16 at ...

以 1\1\GINC-... 开头的那一长串是 Archive 段， 即机器可读的结果摘要——包含方法、基组、分子名、日期、路由段、 以及最终的优化坐标。实际用途：用脚本批量处理成百上千个 log 时， 只要 grep 这一行就能把结果全部抽出来，不必解析整个日志。

4.3 部分优化：只优化一部分坐标

有时我们想固定某些结构参数（例如冻结分子骨架、扫描某根键）。 做法是把变量分成两组：

五、算例 2–6：用对称性把输入写短对应讲义 PPT 第 24–32 页

这一节的五个分子，技术主线只有一条：把分子点群的对称性翻译成变量之间的约束。 点群越高，能共用的变量就越多，输入文件就越短，优化也越稳。

5.1 平面分子：乙烯 C₂H₄

# HF/STO-3G OPT

C2H4 opt

0 1
C
C   1  r1
H   1  r2  2  a1
H   1  r2  2  a1  3  180.0
H   2  r2  1  a1  3  0.0
H   2  r2  1  a1  4  0.0

r1=1.32
r2=1.09
a1=120.0

5.2 高对称分子：CH₃F（点群 C₃ᵥ）

先判断分子的点群，再用对称性把变量压到最少。 CH₃F 属于 $C_{3v}$：C—F 是一根轴，三个 C—H 完全等价。

# HF/STO-3G OPT

CH3F C3v opt

0 1
C
F   1  r1
H   1  r2  2  a1
H   1  r2  2  a1  3  b
H   1  r2  2  a1  3  -b

r1=1.38
r2=1.09
a1=110.6
b=120.0

5.3 同族分子：NH₃（点群 C₃ᵥ）

氨与 CH₃F 的点群相同（都是 $C_{3v}$），因此输入文件的结构完全一样， 只需换元素和初始数值：

# HF/STO-3G OPT

NH3 C3v opt

0 1
N
H   1  r1
H   1  r1  2  a1
H   1  r1  2  a1  3  b
H   1  r1  2  a1  3  -b

r1=1.01
a1=107.0
b=120.0

5.4 环状分子：呋喃（点群 C₂ᵥ）

呋喃（furan，C₄H₄O）是一个五元杂环芳香分子，属 $C_{2v}$ 点群—— 一条 $C_2$ 轴、两个镜面。

5.5 线性分子：乙炔（D_∞h）

线性分子（点群 $D_{\infty h}$）的特殊之处：键角恒为 180°， 二面角则完全无定义（绕轴旋转不改变任何东西）。

# HF/6-31G OPT

acetylene linear opt

0 1
C
C   1  r1
H   2  r2  1  180.0
H   1  r2  2  180.0

r1=1.20
r2=1.06

5.6 算例 2–6 的共同套路

5.7 读图：C_3v 的对称性信息与特征标表

查一次点群，拿到的是两样东西：这个分子有哪些对称元素（左图）， 以及它所属点群的特征标表（右图）。

六、作业与提交要求

作业

提交要求

本讲小结

附：PPT 页码对照表

本讲为讲义体系的第 2 讲，数字页码与源课件 PPT 页码一一对应（1–33）， 另加 6 张插页（PPT 第 3 页之后 5 张讲 Gaussian 发展简史、 第 27 页之后 1 张读点群信息的读图页），故放映文件共 39 页。 深链有两种写法：#p=N 按页码跳转（N = 1…33）， #s=M 按页序跳转（M = 1…39；本讲插页的页序为 4–8 与 33， 其后的页码相应后移）。