烯丙基正离子旋转异构反应的计算化学实验设计
王亚妮，张学鹏
陕西师范大学化学化工学院，西安 710119
随着计算机技术的不断突破和量子化学理论的不断完善，理论计算在材料、催化、合成、生物医药等不同领域都得到了广泛应用[1]。计算化学已受到广大科研工作者的高度重视，使学生学习并掌握一定计算化学实验技能也愈发必要。
目前，国内很多高校都已将计算化学引入本科或低年级研究生教学[2-6]。自新冠肺炎疫情爆发以来，为响应国家抗击疫情的号召，全国各大高校按照教育部提出的“停课不停教、停课不停学”的要求，积极开展了多种形式的线上教学。在化学专业课程的线上教学过程中，一些知识点比较抽象，学生较难快速地理解和掌握，同时教师也较难及时了解学生的学习进度。特别是在化学实验课的线上教学过程中，学生无法通过线上课堂进行实验操作，教师无法对仪器设备、药品进行实物讲解，使得学生对实验现象的了解以及实验整体的把握都不好。相对地，计算化学可以较好地适应线上教学，实验安全性高，教学成本低，并且可以帮助学生更好地了解和掌握微观抽象的化学知识点。在实验过程中，教师与学生只需通过计算机，就可以共同完成实验，教师可以根据计算结果，实时掌握学生的实验进展，对存在的问题进行第一时间的沟通与处理。并且计算化学可以将抽象的反应势能面、过渡态、分子轨道、原子电荷、静电势等概念进行图形化处理，让学生可以直观、深入地掌握这些知识[7]。计算化学实验的开设不仅可以帮助学生掌握计算化学的基本概念与实验技能，还可以与无机化学、物理化学、有机化学等专业课程的理论知识结合，帮助学生加深对一些复杂的化学概念的理解[8-12]。
本文以烯丙基正离子的旋转异构反应为例，设计了一个完整的计算化学实验，该实验设计可以较好地进行线下与线上教学。烯丙基正离子的分子结构简单，计算耗时较少，可以在规定课程时间内完成教学。学生通过该实验设计能够进一步加深旋转异构反应机理的掌握，包含反应物结构的确定、反应势能面的构建、活化能的计算。我们希望能够通过此计算化学实验课程设计，帮助学生掌握基本的计算化学概念与实验技能，从而对其后面的课程学习以及科研有所帮助。
1)了解计算化学基本概念，了解结构优化、单点能计算、过渡态寻找、布局分析等计算化学常用研究手段。
2)了解和掌握烯丙基正离子旋转异构反应机理。
随着计算机性能的不断提升和计算化学理论与程序的不断完善，计算化学已经成为化学、材料、环境等领域不可或缺的研究手段之一。目前，常用的电子结构理论基于量子力学的基本原理和方法，通过求解薛定谔方程得到分子的能量和相关性质。薛定谔方程可以写为：Ψ = EΨ，式中为哈密顿(Hamiltonian)算符，Ψ为波函数，E为能量。分子体系的算符由以下几项组成：
基于不同的处理手段，现今常用的电子结构理论有从头算方法(ab initio)、半经验方法以及密度泛函理论(Density Functional Theory，DFT)。从头算方法在求解薛定谔方程时对哈密顿算符不做过多的简化，不引入经验参数来拟合双电子积分[13]。目前，从头算方法包括基于Hartree-Fock方程的HF方法[14]，以及在HF方法基础上引入电子相关作用而发展起来的后HF (Post-HF)，主要有多级微扰方法(Many-body Perturbation method，MPn)[15]、耦合簇方法(Coupled Cluster method，CC)[16]、组态相关作用法(Configuration Interaction method，CI)[17-19]等计算方法。HF方法引入了波恩-奥本海默绝热近似、非相对论近似和单电子近似，考虑了自旋电子之间的交换作用，但忽略了电子间的库伦相关作用，对于弱相互作用(比如氢键、卤键、π-相互作用等)比较明显的体系描述不好。半经验方法在求解薛定谔方程时通过省略多中心积分或者直接使用拟合的经验参数来近似求解自洽场分子轨道方程，极大提高了计算效率。PM6是目前使用比较广泛的半经验方法[20]，其支持过渡金属(支持上限到铋元素)，在大分子量体系计算、结构初猜与筛选等方面有较广泛的应用。密度泛函理论使用电子密度函数取代多电子波函数来描述相互作用，极大简化了薛定谔方程的求解。密度泛函理论[21-23]依然引入了绝热近似、非相对论近似以及单电子近似，与HF方法相比，DFT在单电子近似的基础上进一步考虑了交换能与关联能，对电子相关作用描述得更好一些，但依然不能较好地描述弱相互作用。DFT方法的计算速度比HF方法快或者相当，但计算精度大大增加。本文中使用的B3LYP泛函[24]是最常使用的DFT泛函之一，经常用于体系的结构优化以及能量计算。
量化计算中用基组来构建体系的波函数，常用的基组由多个高斯型函数(GTO)组成。一般基组越大，计算的精度也越高。例如6-311++G(d,p)，其表示用6个GTO来构建一个内层轨道，外层价轨道劈裂为三组，分别用3个、1个和1个GTO来构建，“d”表示给非氢原子添加d轨道，“p”表示给氢原子添加p轨道，第一个和第二个“+”分别表示给非氢原子和氢原子添加弥散函数。通过给轨道添加d或p角动量可以改变轨道形状，从而易与其他轨道成键。另一方面，通过给轨道添加弥散函数，可以使轨道占据更大的空间，从而有利于弱相互作用的描述。
Gaussian是目前使用最广泛的计算化学软件之一，它功能强大，简单易学，支持多种计算方法。GaussView是辅助Gaussian程序使用的图形界面软件，主要用于分子建模和计算结果分析。考虑到Gaussian和GaussView为收费软件，本实验也可以通过Chem3D、HyperChem等免费软件进行分子建模，使用开源软件ORCA[25]、GAMESS-US[26]、Dalton[27]等完成量化计算。
烯丙基正离子常在SN1亲核取代反应中产生，是有机反应重要的活泼中间体，其可以发生氧化、偶联、亲电取代等多类化学反应。著名的Friedel-Crafts烷基化反应[28]是典型的芳烃烯丙基化的方法。Suji-Trost烯丙基取代反应[29]可以合成出多种用于反应的过渡金属和亲核试剂。深入了解烯丙基正离子的反应活性和化学性质有利于帮助我们更好地掌握相关有机反应的机理，进一步设计合成更多天然化合物。目前，烯丙基正离子旋转异构反应的机理一般认为通过以下三个步骤进行：① 烯丙基正离子中间碳原子上的氢原子迁移到α位碳原子上，形成甲基；
② 发生甲基与烯烃之间C―C键旋转；
③ 氢原子从α位碳原子迁移回中间碳原子上。
根据分子结构的大小，以及计算耗时的考虑，本实验主要使用DFT的B3LYP泛函和6-31G(d,p)基组对反应物结构进行优化和频率计算，并在此基础上使用更高精度的B3LYP/6-311+G(d,p)方法计算反应物的热力学能。另外，我们还将B3LYP的结果与HF、PM6计算的结果进行了对比，从而让学生更加直观地了解计算化学方法的差异。需要注意的是，本实验设计的教学内容较多，常规4小时的实验时间可能无法完成全部的实验内容，指导教师可以基于学情酌情选择适当的教学内容。
台式机或者笔记本电脑(安装Windows XP，Vista 8, 10或11)，Gaussian 09[30]软件，GaussView 5.0或6.0软件，Origin软件，Chemdraw软件。
4.1 构建分子初始结构
使用GaussView软件构建反应物(reactant，RC)、反应过渡态(transition state，TS1、TS2和TS3)、反应中间体(intermediate，IM1和IM2)以及产物(product，PC)的初始分子结构，其结构如补充材料中图S1所示。具体操作如下：点击Create New Molecule Group按钮建立工作区，然后点击Element Fragment按钮在弹出的元素周期表中选择需要的元素以及构型，在工作区构建各反应分子。烯烃基团也可以通过点击R-Group Fragment按钮在弹出的窗口中直接选取烯烃基团。分子建模提示：由于烯丙基正离子的α位碳原子失去了一对电子，带有部分正电荷，其p轨道和C＝C双键的π键之间具有p-π共轭效应，所以反应物RC应呈平面结构。过渡态TS1涉及氢迁移，具有即将断裂和即将形成的C―H键，所以应有一个氢原子位于中间碳原子和α位碳原子的上方。氢迁移完成后，产生一个中间体IM1。IM1通过甲基旋转过渡态TS2产生中间体IM2。IM2再发生氢迁移，得到过渡态TS3。最后氢原子重新回到中间碳原子上后，得到仍是平面结构的产物PC。
4.2 编辑输入文件
使用GaussView完成分子建模后，点击Calculate → Gaussian Calculation Setup进行计算命令设置，如图1(a)。
图1 (a) GaussView软件的计算命令设置界面；
(b) Gaussian 09软件执行量化计算界面
(1) Job Type部分：其中Energy表示能量计算，Optimization表示结构优化计算，Frequency表示分子频率计算，Opt + Freq表示分子结构优化后继续计算分子的频率，IRC表示计算过渡态分子的内禀反应坐标，Scan表示对分子的键长、键角、二面角等参数做势能面扫描，Stability和NMR分别表示稳定性计算和核磁谱图计算。本实验中依次选择Opt + Freq和Energy完成结构优化和热力学能的计算。本界面其他选项采用默认设置。
(2) Method部分：用来选择实验所用的计算方法。当Job Type选择Opt + Freq时，依次选择DFT →B3LYP泛函和6-31G(d,p)基组，设置结构优化的命令；
当Job Type选择Energy时，依次选择DFT →B3LYP泛函和6-311+G(d,p)基组，设置能量计算的命令。方法对比时可以选择Semi-empirical → PM6或Hartree-Fock设置半经验或从头算计算。本界面Charge统一设为1，表示整个体系带一个正电荷，Spin设为Singlet，表示体系为单线态。
(3) Title部分：可以自由设置计算任务的标题。
(4) Link 0部分：计算控制设置。内存(Memory)和CPU核数(Shared Processors)可以采用默认设置，或者自定义(内存最大可设为2 GB，CPU核数建议为电脑核数的一半)。Chkpoint File指定chk文件的命名和保存路径(chk文件在本实验中可以用来生成分子轨道、静电势图)。其他选项不设置。
(5) General、Guess部分本实验不设置。
(6) Pop部分：设置原子电荷的计算命令，也可以在Additional Keyword部分中直接输入原子电荷计算的关键词，比如Pop = NPA，表示计算NPA电荷。
(7) Preview部分：查看编辑完成的输入文件的内容，RC结构优化的输入文件可参考表1。为了简化学生建模与输入文件编辑过程，本实验使用的所有输入文件都提供在补充材料中(Supporting Information)。
表1 RC分子结构优化计算输入文件的编辑和解释
输入文件设置好以后点击左下角的Submit按钮运行Gaussian 09进行量化计算。或者也可以先打开Gaussian 09，选择File，点击Open，打开编辑完成的输入文件，检查各命令行的关键词无误，点击Run命令，开始量化计算，如图1(b)。
4.3 寻找过渡态和IRC计算
反应起始物RC、中间体IM1和IM2，以及产物PC可以通过4.2小节中计算步骤的描述直接进行结构优化计算得到它们稳定的分子结构。对于反应过渡态TS1、TS2和TS3，本实验使用TS法寻找其稳定的结构[31]，在编辑过渡态计算的输入文件时，Job Type部分选择Opt + Freq，设定Optimize to a为TS(Berny)，Force Constants为Calculate at First Point。得到优化后的结构后，点击Result → Vibrations查看其是否有且只有一个虚频(过渡态结构有且只有一个数值为负的振动频率，即虚频(imaginary frequency))。所有分子优化结束后，按4.2小节所描述的步骤对它们进行能量计算。
点击File → Save将优化后的过渡态结构保存为新的输入文件，在Job Type中选择IRC计算，设定Force Constants为Calculate Once，点击Compute more points，设定N = 200，完成IRC计算的输入文件编辑。最后可以得到由过渡态结构出发，连接反应物和产物两个极小点的势能面曲线，进一步确认该过渡态结构是否正确。需要注意的是，如果IRC计算报错，可以一开始在Additional Keyword部分输入IRC = LQA。IRC计算比较复杂，可以作为教师演示实验或者给有兴趣的同学额外计算。
4.4 原子电荷计算
本实验使用Mulliken、NPA、Merz-Kollman (MK)原子电荷计算方法[32-37]，计算过渡态TS1分子的原子电荷。原子电荷计算可以通过Calculate → Gaussian Calculation Setup → Pop设置，也可以直接在Additional Keyword里添加Pop = NPA或Pop = MK关键词。Mulliken电荷由Gaussian 09软件自动计算。上述三种不同的原子电荷计算方法在下面的章节会有简单介绍。
5.1 分子的结构优化
图2为各反应分子优化后的结构示意图。由此图可知，此旋转异构反应的RC分子的C1―C2和C1―C3键长大小接近，分子呈平面构型，证实了此分子具有p-π共轭效应。在TS1分子的结构中，C1―C2键长减小更接近C＝C，C1―H6键长增加，逐渐生成C3―H6键，键角A (C1―H6―C3)为56.6°，表示H6原子开始向C3迁移。IM1分子的结构中，C1―C2键长为1.27 Å (1 Å = 0.1 nm)，为C＝C双键；
C1―C3键长为1.41 Å，为C―C单键，此时H6原子迁移到C3原子上，H6与H5同侧，C3原子由一开始的sp2杂化转换为sp3杂化。甲基旋转过渡态TS2中发生C1―C3单键旋转，二面角D (H5―C2―C3―H6)由IM1中的-0.2°旋转为TS2中的90.2°。IM2相比IM1，发生了甲基旋转，H6与H5对侧。H6原子由甲基迁移回到C1原子上的TS3与TS1类似。烯丙基正离子旋转异构反应的PC回到与RC类似的结构，但相比于RC，PC中H6与H4同侧，H6的位置发生了改变。
图2 B3LYP/6-31G(d,p)方法下优化的分子结构以及原子序号、关键键长(Å)和二面角数据
表2为B3LYP/6-31G(d,p)、B3LYP/6-31G(d)、HF/6-31G(d)和PM6方法下，优化TS1分子的结构参数和耗时数据。由表2可知：PM6方法优化的TS1分子中将要断裂的C1―H6键长明显比HF和B3LYP方法大，HF和B3LYP方法差异较小。耗时方面，PM6、HF、B3LYP方法耗时依次增加。在B3LYP泛函下，6-31G(d)和6-31G(d,p)基组得到的结构参数差别较小，但后者的耗时较大。该计算方法对比内容可以作为教师演示实验或者给有兴趣的同学额外计算。
表2 B3LYP/6-31G(d,p)、B3LYP/6-31G(d)、HF/6-31G(d)和PM6方法优化TS1分子的结构参数和耗时数据＊
5.2 热力学数据和反应能垒图
在GaussView软件中，点击Result → Summary → Thermo查看各分子的热力学数据。其中Electronic Energy (EE)是分子的电子能，Thermal Correction to Energy/Enthalpy/Free Energy是分子的内能/焓/吉布斯自由能的校正因子，EE + Thermal Energy/Enthalpy/Free Energy Correction是分子的内能/焓/吉布斯自由能的大小。高精度的焓/吉布斯自由能=高精度单点能+结构优化的焓/吉布斯自由能校正因子。在此旋转异构反应中，各分子的单点能、焓值、吉布斯自由能值以及相对于反应物的能量差值如表3所示。
由表3中的各分子热力学数据和图3的旋转异构反应吉布斯自由能图可知，RC和PC分子的吉布斯自由能大小一致，互为旋转异构体，而且TS1和TS3分子以及IM1和IM2分子的吉布斯自由能大小也分别一致，显示了前后两次氢迁移过程的对称性。甲基旋转过渡态的反应能垒几乎可忽略不计，此旋转异构反应的决速步为氢迁移过程，其吉布斯自由能垒108.8 kJ·mol-1(26.0 kcal·mol-1)，表明此旋转异构反应在室温下较难发生。同时，由表3可知，该烯丙基正离子旋转异构反应各步骤的焓变与自由能变的程度接近，说明反应过程的熵变对反应影响不大。
表3 B3LYP/6-311+G(d,p)//B3LYP/6-31G(d,p)方法下各反应分子的热力学数据＊
图3 B3LYP/6-311+G(d,p)//B3LYP/6-31G(d,p)方法下烯丙基正离子旋转异构反应的吉布斯自由能图
5.3 过渡态和IRC分析
使用GaussView软件查看TS1、TS2和TS3分子的输出文件，点击Result → Vibrations查看频率计算结果，发现这三个分子都只有一个虚频，分别为-835.44 cm-1、-143.62 cm-1和-835.38 cm-1，点击Start Animation可以看出TS1和TS3的虚频对应的分子振动模式都是氢迁移过程，TS2的虚频对应的分子振动模式是甲基旋转过程。点击Result → IRC/Path分析IRC计算结果并导出IRC数据，通过Origin软件作图得到图4所示的以过渡态分子结构为坐标零点，连接反应物和产物的IRC曲线。IRC计算可以进一步确认过渡态TS1、TS2和TS3分子结构的正确性。
图4 B3LYP/6-31G(d,p)方法下TS1、TS2和TS3分子的内禀反应坐标(IRC)
5.4 原子电荷
本实验选取了三种原子电荷计算方法。其中Mulliken方法是Gaussian系统默认的计算原子电荷的方法，其将成键电子平均分配给成键原子；
NPA方法按照体系的波函数划分电子云空间分布，得到原子电荷，常用于分子体系的电荷计算；
MK方法是典型的拟合静电势电荷的计算方法，常用于周期性体系的电荷计算。通过GaussView软件查看原子电荷计算的输出文件，点击Result → Charge Distribution得到在不同原子电荷计算方法下TS1分子上各原子的电荷值。由表4数据知，不同原子电荷计算方法对于TS1分子中不同碳原子和氢原子上电荷值大小的分配存在较大差异。因此，在进行原子电荷计算时，要根据体系的性质和原子电荷计算方法的特点，选择合适的电荷计算方法。
表4 B3LYP/6-31G(d,p)方法下TS1分子在不同电荷计算方法下的原子电荷值＊
5.5 HOMO、LUMO轨道和静电势图
打开Gaussian 09软件，点击Utilities → formchk，将在B3LYP/6-311+G(d,p)下得到的chk文件转换为fchk文件。使用GaussView软件打开fchk文件，点击Tools → Mos查看分子的轨道信息。选中目标分子轨道后，轨道会被标黄，设置Isovalue值等于0.07，点击Visualize → Update，生成相应的分子轨道图并保存，即得到图5所示的各分子HOMO、LUMO轨道图。RC分子的HOMO轨道是一个大共轭π键，表明其3个C原子间存在p-π共轭效应。随着此旋转异构反应的进行，TS1、TS3分子和IM1、IM2分子的HOMO轨道大π键在消失，表明分子的p-π共轭效应也在消失。
图5 B3LYP/6-311+G(d,p)方法下各分子的HOMO、LUMO轨道
当打开fchk文件时，依次点击Result → Surfaces/Contours → Cube Actions → New Cube，在Tpye选项中选择Total Density，生成分子的Cube文件，调整Density值为0.11，点击Surface Actions → New Mapped Surface，在Type选项中选择ESP，得到图6所示的各分子静电势图。由图6可知，RC分子的正电荷被两侧的碳原子平分，这两个碳原子具有相同的反应活性。随着旋转异构反应的进行，TS1分子的两侧碳原子上的正电荷变小，IM1分子上正电荷已被转移到中间碳原子上，之后随着反应的进行正电荷又重新回到两侧碳原子上。
图6 B3LYP/6-311+G(d,p)方法下各分子的静电势图
HOMO、LUMO轨道和静电势图的计算操作较为复杂，可以作为教师演示实验或者给有兴趣的同学额外计算。
本实验面向高年级本科生或低年级研究生开设，若为线上教学需要在课前指导学生下载并安装所需软件，计划使用4课时完成实验教学。具体安排如下：1课时介绍计算化学基本概念和烯丙基正离子旋转异构反应机理；
0.5课时完成分子建模，编辑输入文件；
2.5课时完成各部分量化计算。课后安排学生独立撰写实验报告，完成思考题。
(1)根据以上计算出的反应吉布斯自由能垒大小，判断烯丙基正离子的旋转异构反应在常温常压下能否发生？
(2)根据本实验得到的热力学值和反应势能面图，结合物理化学中的相关知识，思考反应热力学与动力学的区别。
此烯丙基正离子旋转异构反应的计算化学实验课程，包括了计算化学实验的基本概念与实验步骤。整个实验教学内容比较基础，方便没有量化计算基础的学生学习。通过本实验的学习，可以帮助学生认识计算化学，激发学生对计算化学的学习兴趣，还可以更好的理解课本中一些较为抽象的概念。
笔者在三年间多次将此实验应用于大三本科生创新实验、大四本科生开展毕业设计前以及低年级研究生入门计算化学的教学实验中，均取得较好的效果。每个学生单独进行实验，并撰写实验报告，通过较短的时间就可以让学生了解并掌握一定的计算化学实验基本操作。另外，学生只需准备个人电脑就可以完成此实验，十分适合疫情期间的线上教学以及线下课堂的教学。
补充材料：可通过链接http://www.dxhx.pku.edu.cn免费下载。
相关推荐：