Post

Gaussian入门：几何优化、振动分析和能量计算

Posted Updated
By Bane Dysta
27 min read
Gaussian入门：几何优化、振动分析和能量计算

Gaussian是最流行的量子化学计算软件，功能涵盖了绝大部分常规任务。且Gaussian的输入文件绝对是所有主流量子化程序中最简洁的，这一点笔者在学过ORCA、Dalton、NMChem等软件后深有体会。因此，Gaussian是入门量子化学应当首先掌握的软件。而几何优化又是所有任务中最基础的任务，绝大多数量子化学计算都要基于几何优化后的结构来进行。总之，这是一个入门级的几何优化教程贴。

0. 软件安装

Gaussian安装包各位八仙过海各显神通，笔者这里就不提供了。

Windows下安装Gaussian，只要跟着安装程序走就行。需要注意的问题有：

  • 由于笔者是做荧光的，相信来看笔者blog的也是对计算激发态有需求的，那么强烈建议安装Gaussian 16以及配套的GaussView 6，而非Gaussian 09和GaussView 5。因为G09不支持TDDFT的二阶解析梯度，做激发态优化时若要进行频率分析，则必须使用数值梯度，纯属自找麻烦。
  • 如果坚持要使用Gaussian 09，一定要找64位Gaussian的安装包，不要32位的。
  • 要先安装Gaussian，再安装GaussView，并把GaussView安装在Gaussian的文件夹内。
  • 不建议安装在C盘。
  • 安装路径一定要避免中文或空格，且尽量避免除了”-“和”_“以外的特殊字符。
  • 不要乱改G16创建的主文件夹的名称，想自定义可以建立一个上级文件夹，把Gaussian安装在里面。
  • 安装完毕后，打开G16W检查有没有scratch文件夹，若没有可以手动建立。

Linux版Gaussian是预编译版本，只需要使用tar -xvf name.tar将压缩包解压到一处，然后设置环境变量：

1
2
3

echo 'export g16root=/path/to/g16' >> ~/.bashrc
echo 'export GAUSS_SCRDIR=/path/to/g16/scratch' >> ~/.bashrc
echo 'source /path/to/g16/bsd/g16.profile' >> ~/.bashrc

注意将/path/to/g16替换为你的实际解压路径。最后在g16目录下输入chmod 750 -R *加上可执行权限，Gaussian就安装完毕了。

1. 搭建结构

我们平时看的lewis式或键线式都是2D结构，不包含3D信息。因此需要首先将其转换成相对合理的3D坐标，才能拿给程序去计算。这里笔者以香豆素为例，介绍一下转换方法。

1.1 使用Chem Office套件转换

首先，在ChemDraw上绘制香豆素的结构。能看到这的人应该不用我教ChemDraw了，如果你不会，请你关闭此网站并认真完成今天的家庭作业。在ChemDraw上画好2D结构后，Ctrl+C复制，打开Chem3D，然后Ctrl+V粘贴，就得到了一个3D结构。如果你用的不是老爷机，可以使用Ctrl+M进行MM2初步优化，得到一个相对合理一些的结构。完成后，使用Ctrl+Shift+S另存为Gaussian Input文件（在后缀名下拉菜单里找到gjf并保存）。

最后以记事本格式打开产生的gjf文件，把#后面自动产生的sp test删掉，这样结构就准备好了。使用Chem3D转换结构的坏处是Chem3D提供的初始结构可能相当不合理，需要依照化学直觉检查有没有什么离谱的地方。

注意：如果你使用了MM2优化，你还需要额外查找文件内以”Lp”字样开头的行并将其全部删掉。

1.2 使用GaussView搭建结构

这个笔者自己也用不太明白，暂时略过，以后看心情补。GaussView的好处是可以自己调整立体结构，比Chem3D灵活得多，且很适合用点群工具搭建对称性体系，但是笔者已经太习惯ChemDraw延申碳碳键的方式了，换GaussView很不适应。

1.3 使用OpenBabel转换smiles字符串

有些网站会提供smiles字符串，可以使用OpenBabel将其转换为3D结构，适合写批量转换脚本，有一定命令行基础的读者可以尝试。去OpenBabel官网下载并安装好OpenBabel后，在命令行中运行：

1

obabel -:"O=C1OC2=CC=CC=C2C=C1" -O Coumarin.gjf --gen3d

即可在当前目录生成Coumarin.gjf文件。注意这里生成的的gjf文件没有title行，也没有计算参数行，得手动补上。这种方式适合批量获取化合物的3D初始结构。

1.4 从数据库中下载

很多已知的分子可以从一些化学数据库，如PubChem上找到已经使用力场预优化过的较为合理的3D结构。在PubChem中搜索Coumairn，点进第一个链接，在右侧点击Download，会弹出一个选项卡。在3D Conformer下，可以看到PubChem提供SDF、JSON、XML、ASNT格式的几何文件。笔者习惯下载SDF文件，然后用OpenBabel命令：

1

obabel name.sdf -O name.gjf

将其转换为gjf格式。这里生成的的gjf文件也没有title行和计算参数行。

Chemdraw绘制出来的结构可以使用Ctrl+Alt+C复制为smiles字符串，可以粘贴进PubChem来搜索。

2. 输入文件准备

刚刚拿到的输入文件只有几何信息，没有计算相关信息，还需要告诉程序计算相关的设置才能开始计算。以香豆素为例，进行结构优化的完整输入文件应该是类似这样的：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

%mem=16gb                           //分配内存
%nproc=6                            //分配CPU核心数
%chk=D:\test\Coumarin.chk           //chk文件保存路径
# opt freq wB97XD/def2SVP           //关键词行

title                               //title行

 0 1                                //电荷与自旋多重度
O          -1.05750         1.12680         0.00040
O          -3.31910         0.76440        -0.00070
C           0.44750        -0.78620        -0.00010
...
H           3.45870         1.67230         0.00000
H          -2.81290        -1.79660         0.00060

                                    //几何坐标输入结束后，留空两行

2.1 电荷与自旋多重度

虽然几何结构已经有了，但是程序自动给出的电荷和自旋不一定是对的，需要手动检查。

电荷就是体系的总电荷数，这个不必多说。

自旋多重度是α电子数-β电子数+1，这个解释起来稍微麻烦一些，可以姑且认为常见有机体系自旋均为1。但是如果涉及自由基、金属、卡宾、氧气等情况，就另当别论了。可以参考这个帖子里的讨论。

2.2 基组的选择

基组是描述电子在空间中分布的基函数集，用来近似表示分子体系中电子的波函数。笔者常用的基组系列有：

  • pople系列基组：如6-31g、6-311g等
  • Ahlrichs系列基组：如def-TZVP、def2-SVP等
  • Dunning’s相关一致性基组：如cc-pVTZ、cc-pVQZ等

其中，DFT计算常用前两个系列的基组，而Dunning’s相关一致性基组是给post-HF计算设计的，并不适合DFT计算。

如果想要完美描述电子的分布，所需要的代价太大，因此实际计上常常使用尺寸小一些的基组进行计算。由于价电子是最活跃的电子，参与了大多数的化学键和分子间相互作用，应当用更精确的方式来描述价电子的行为。因此在实际计算中，内层常用单组基函数来描述，而价层则常用多组基函数来描述，这称为基组的分裂价ζ。如果用两组基函数描述价层，则称为2-ζ基组，如果用三组基函数描述价层，则称为3-ζ基组，依此类推。分裂价越高，基组尺寸越大，计算代价也越大。但如果选择的基组分裂价太低，则计算精确性会受到严重影响。因此，适当选用基组档次是很重要的。

除了分裂价，还有两个地方值得注意：

  • 极化函数。添加极化函数表示给原子添加一个更高角动量的轨道，允许电子在不同方向更灵活地分布。通常重原子和氢原子的极化是可以分开加的。例如，6-31g给重原子加一个d轨道，就是6-31g(d)，再给氢加一个p轨道，就是6-31g(d,p)
  • 弥散函数。添加弥散函数表示给原子添加一个衰减指数较小的基函数，来描述远离原子核的电子分布。通常重原子和氢原子的弥散也是可以分开加的。例如6-31+g*表示给重原子加一层弥散函数，6-31++g*表示再给氢原子加一层弥散函数。

具体基组档次怎么选择、弥散函数和极化函数做什么任务的时候要加、怎么加的问题很复杂。简单来说：

  • 几何优化和频率分析用2-ζ基组就行，单独计算能量可以用3-ζ或4-ζ基组；
  • 重原子的极化一定要加，而氢原子极化在计算直接与氢相关的问题才必须加；
  • 重原子弥散在算阴离子能量时等特殊情况需要加，而氢原子弥散几乎不用考虑。

更详细的说明可以参考谈谈量子化学中基组的选择。这里对初学者给出常用基组尺寸的大致排序和简单的建议：

1
2

STO-3G < 3-21G < 6-31G(d) < def2-SVP ≈ 6-31G(d,p) << 6-311G(d,p) < def-TZVP < def2-TZVP << def2-QZVP
          2-ζ                                           3-ζ                  几何优化封顶      4-ζ
  • 除非只是算着玩玩，否则但凡你想得到任何有实际价值的结果，绝对不要考虑STO-3G。
  • 前四周期体系的几何优化和振动分析用6-31g*，如果体系中有第五周期以上的元素，需要改用def2-SVP。
  • 如果是含大量第四周期过渡金属的体系，建议使用SDD赝势基组，参考详解Gaussian中混合基组、自定义基组和赝势基组的输入第6节。
  • 基于优化好的几何结构计算单点能时用def2-TZVP，计算阴离子能量等时使用6-311+g**或ma-TZVP(需要自己引用基组定义)。

基组在Gaussian里的关键词不一定与基组名相同，需要查阅Gaussian的基组手册获取正确的关键词。有的基组名和关键词差别很大，如def-TZVP的关键词是TZVP，不看手册根本猜不到的。此处附常用基组名和关键词、同义词对照表：

Basis SetKeywordsynonymBasis SetKeyword
6-31g(d)6-31g(d)6-31g*def2-SVPdef2SVP
6-31g(d,p)6-31g(d,p)6-31g**def-TZVPTZVP
6-311+g(d,p)6-311+g(d,p)6-311+g**def2-TZVPdef2TZVP

2.3 泛函的选择

密度泛函是基于Kohn-Shan方程用于计算体系电子相关能量的函数表达式，核心思想是通过电子密度来确定体系的能量，而不是直接处理多电子波函数。密度泛函是一种映射关系，将电子密度映射到体系的能量上，因此不同的泛函会影响到计算结果的精确度。

由于DFT理论对电子库伦相关与交换作用描述并不充分，可以引入Hartree-Fock相关项与MP2交换项来改进计算结果。引入Hartree-Fock相关项的泛函称为杂化泛函，同时引入Hartree-Fock相关项与MP2交换项的泛函称为双杂化泛函。

泛函的选择也很复杂，同样可以参考sobevera老师的博文。这里也给出大致的建议：

  • 纯泛函：一般任务不建议使用，会高估电子离域性，且低估激发能。对于ORCA用户，可以尝试基于纯泛函B97引入多项校正的组合方法B97-3c，比较便宜。
  • 杂化泛函：基态非共轭体系几何优化和振动分析推荐B3LYP，大共轭体系推荐ωB97X-D、cam-B3LYP；LE激发推荐PBE0，CT激发推荐ωB97X-D、cam-B3LYP；计算主族体系能量推荐M06-2X；过渡金属推荐TPSSh；计算核磁推荐revTPSS；对于ORCA用户，还推荐使用ωB97M-V计算高精度能量。
  • 双杂化泛函：仅建议高精度计算能量时使用。Gaussian对B2PLYP-D3的支持比较好；DSD-PBEP86-D3(BJ)可以用来计算激发态；对于ORCA用户，PWPB95-D3(BJ)是比较稳健的选择。

泛函在Gaussian里的关键词同样不一定与泛函名相同，需要查阅Gaussian的泛函手册获取正确的关键词。有的泛函名和关键词差别很大，如PBE0的关键词是PBE1PBE。更糟糕的是你直接写PBE0甚至不会像基组一样写错了会报错，而是会调用精度很差的PBE0-DH双杂化泛函，要特别注意。此处附Gaussian常用杂化泛函和关键词对照表：

FunctionalKeywordFunctionalKeyword
B3LYPB3LYPcam-B3LYPcam-B3LYP
PBE0PBE1PBEM06-2XM062X
ωB97X-DwB97XDTPSShTPSSh

2.4 计算关键词

在gjf文件中，关键词行是控制量子化学计算内容的核心部分，包含了计算任务、计算方法、计算参数等所有计算相关的选项。书写规则：

  • 位于核数、内存、检查点文件设置之后，位于title行之前
  • 以”#”开头，关键词之间用空格分隔，且不区分大小写、不区分顺序
  • 在不会引起歧义的前提下，可以进行省略。如opt=ModRedundant可简写为opt=modr，因为opt的参数中以”modr”开头的只有ModRedundant。
  • 由于结构优化和频率分析的计算级别需要严格一致，多数量化程序都支持结构优化和频率分析合并成一个任务。但除这两个可以一起做以外，其他的计算任务关键词不能写在一起。

首先应当明确的是计算任务。结构优化的关键词为opt，振动分析的关键词是freq，能量的关键词是sp或留空。如前所述，几何优化和频率分析一般是一起做的，可以写opt freq。而如果你要计算能量，习惯上是不写任务关键词的。

然后是计算水平，即泛函/基组，如前所述，找到对应关键词后添加进关键词行即可。香豆素是共轭体系，笔者选择使用ωB97X-D和def2-SVP进行计算，则关键词为wB97XD/def2SVP

定下任务，选好泛函和基组之后，关键词行剩下的就是些边边角角的琐碎选择了。

  • 色散校正：由于泛函对色散行为描述很差而引入的一个校正项，计算弱相互作用体系时需要添加。常用的D3色散校正号称计算器手算都能算出来，耗时可以忽略不计，所以如果不确定体系有没有弱相互作用，可以总是添加D3校正，尤其是对于B3LYP。通常使用的D3校正关键词为em=gd3bj。

    注意不同泛函对色散校正的支持情况不同，B3LYP支持BJ阻尼形式的D3校正em=gd3BJ，而M06-2X只支持0阻尼形式的D3校正em=gd3，ωB97X-D则自带D2校正，不用额外添加。

  • 溶剂模型：使用scrf=(solvent=water)来添加隐式溶剂，默认模型为IEFPCM。在计算能量时，可以使用scrf=(SMD,solvent=water)`切换为SMD溶剂模型，以更好地描述非极性部分。water可以更改为Gaussian支持的任何溶剂，通常关键词为该化合物的系统命名。Gaussian的溶剂手册列出了Gaussian支持的所有溶剂与对应的关键词。

2.5 计算资源分配

如果是个人笔记本或台式机来运行Gaussian，笔者建议核数按照(计算机物理核心数-2)来给即可，而内存最好每个核分配2-3gb，最少每个核要1gb，不然可能会拖慢计算速度。

为了提供一些参考，假设配置为12核，内存24GB，计算级别选择B3LYP/6-31g*，这样算小型体系是算的动的，而中型体系就有些吃力了，较大体系基本没戏。笔者平时计算惯用配置是32核+96gb，用比较好的级别如wb97XD/TZVP计算中型体系也没什么问题，而较大体系换小级别也能应付一下。

一般chk文件就保存在当前目录即可。注意linux版Gaussian输入文件的chk行可以不写绝对路径。

3. 提交计算

提交计算是很简单的。windows提交计算可以直接把gjf拖进Gaussian程序，然后点开始；linux版Gaussian需要命令行启动，命令为：

1

g16 coumarin.gjf > coumarin.log

当计算结束时，log或者out文件末尾会打印类似于这样的输出：

1

Normal termination of Gaussian 16 at Fri Nov  1 18:40:19 2024.

此时就可以用GaussView打开log文件或者out文件查看结果了。多数简单的信息可以在右键菜单的Results选项卡中找到对应的输出，但也有少数信息需要阅读输出文件才能得到，这里就不再多讲。

能量计算的输出可以在Results - Summary处查看，在Overview选项卡列出的表格中，有一行以E(计算级别)开头，后面的数值即为计算出的绝对能量。绝对能量值不是化学上感兴趣的量，应当横向与其他化合物进行对比。

对于做了opt+freq的任务，要判断几何结构是否真正收敛到了势能面极小点，可以在右键菜单中找到Results - Vibration，检查最上面有没有频率值为负数的虚频。如果没有虚频，说明当前几何结构确实是能量最低结构。而如果存在虚频，则说明当前结构处于势能面上的鞍点，没有收敛到真正的极小点。这种情况下，还需要将log文件Ctrl+S保存为gjf文件，更改计算设置后重新进行一次几何计算。推荐尝试opt=calcfc。

报错是很正常的。笔者初学时，花了一个星期才写出来第一个能从头到尾不报错跑完整的计算任务。对于计算时遇到的报错问题，可以在这个帖子里查找解决办法，或按照帖子所述提取关键信息发到论坛或者群里求助别人。

Quantum Chemistry, Gaussian
Optimization
This post is licensed under CC BY 4.0 by the author.