MECP简介(3)——实例分享

Outline

• Introduction - What is MECP?
• Algorithm - How to locate a MECP?
• Application - When do we need MECP?
• Theory - What can we do beyond MECP?

Demo

Ozone (O3)

From Wikipedia

O3的基态是单重态（事实上是一个单重态双自由基），三重态略高于单重态（约20 kcal/mol），键长键角都略有不同。

Preparation: Gaussian input

easyMECP提供了很方便的高斯接口，我们首先选取一个初始结构，构建高斯输入文件：

这里想要找的MECP是单重态和三重态的交叉点，所以自旋多重度的地方写为{1,3}，easyMECP会自动成两个高斯输入文件，自旋多重度分别为1和3，其它多重度的计算也是类似的。输入文件中大括号的部分会被easyMECP识别，不过笔者测试时发现只有分子坐标之前的大括号会被识别并写入各自自旋态的输入文件，坐标后面的部分不会被识别（不过通常坐标后的部分，例如基组赝势等，应当对不同自旋态都是一样的）。

输入文件中的force关键词是必须的，因为MECP算法需要知道当前结构两个的能量及梯度才能进行优化。在通常的高斯计算任务中，如果加入opt关键词，那么做完能量计算后高斯会自动计算梯度然后进行结构优化。如果不加入opt关键词，仅仅做单点计算，高斯默认是不计算梯度。在计算MECP时，我们需要利用高斯外部的优化算法，但是又想利用高斯的（单点）梯度计算，就需要指定force关键词了。guess=read不是必须的，对于起始结构，如果不加guess=read，easyMECP会首先让高斯计算两个单点，对于后续的计算自动补入guess=read节省计算时间。但是对于特殊情况，比如这里的O3，因为单重态最稳定的电子结构事实上是一个双自由基结构，需要事先用broken symmetry的计算算好，所以guess=read是必须的。对于MECP的计算笔者还是推荐事先做好单点计算，检查好各自自旋态的电子结构都是稳定的之后，再进行MECP的计算，因为过渡金属体系往往同一个多重度也能找到多个自旋态，如果盲目计算很可能得到的中间体、MECP都不是能量最低的。

Execution

执行easyMECP就十分简单了，只需指定好python、easyMECP脚本及输入文件的路径即可。如果高斯不能通过g09直接调用则还需指定高斯的路径。最后的--with_freq关键词指定easyMECP进行振动频率的计算，后面还会详细讲解。

Output

在shell里面跑上述命令，或者将上述命令写成一个bash脚本（下图的runmecp.sh）提交到服务器执行，即可得到下列输出文件：

其中下划线开头的及.conf文件都是easyMECP的配置文件，MECP.x是easyMECP用Fortran代码编译的执行文件。输出文件主要是ReportFile也就是easyMECP的log文件，每一步高斯的计算文件保留在了JOBS文件夹。

ReportFile主要内容如下：

输出文件的开头包含起始结构坐标，两个自旋态各自的能量，然后检查是否收敛并打印相应梯度。

优化若干步之后，MECP结构优化收敛，可以看到此时两个自旋态的能量十分接近，各个梯度的值也很小。

最终优化好的结构及能量如下：

可以看到MECP的能量（电子能）是必然高于两个自旋态的中间体能量的。然而在计算反应机理时，我们往往需要比较的是不同物种的自由能而非电子能，然而MECP并不是势能面上的驻点(stationary point)，正常的频率分析得到的能量没有物理意义，如何得到可比的MECP的自由能呢？

这里我们需要做一个额外的投影频率(projected frequency)的计算。因为MECP在两个自旋态势能面（3N-6维）的交叉面（3N-7维）上是最低点，只在一个维度上不是，所以我们可以在3N-7维的子空间内计算振动频率，而只忽略反应坐标(reaction coordinate)方向的频率。这一过程可以通过高斯关键词freq=projected来完成，不过如果加了--with_freq关键词easyMECP会自动进行投影频率的计算，从而得到MECP的自由能。与找过渡态一样，我们也可以以MECP为起点做IRC计算，确认MECP连接是我们想要的中间体。

除了计算自由能之外，频率计算还有一个重要的目的，就是检查优化得到的结构有没有多余的虚频。在这个例子中，我们可以看到easyMECP在提取自由能时同时提取了每个态的最小振动频率，其中一个(triplet)是有虚频的，这说明优化得到的MECP还不是最稳定的CP，算法之前没能得到最稳定的CP的原因是因为初始结构具有C2v的对称性，所以限制了算法的搜索范围。通过调节键长破坏对称性或者加nosymm关键词重新优化MECP的结构如下：

可以看到新的MECP与之前的相比，能量变化不大，仅降低了0.4 kcal/mol，但是结构差距很大，从对称的1.31 Å键长变成了明显不对称的1.26 & 1.39 Å键长。