本文主要介绍多元醇酯（1-(癸酰氧基)-3-(壬酰氧基)丙-2-基十二烷酸酯）的构象搜索。该体系是一个柔性大分子体系，一定存在很多的构象，所以做构象搜索来分析其结构显得非常有必要。本实例采用sob老师开发的Molclus程序做构象搜索工作，先使用xTB做100 ps的动力学轨迹，先后经过GFN0-xTB和GFN2-xTB的结构优化筛选去除重复结构，再挑选能量差异在3 kcal/mol内的结构进行DFT方法的构象筛选工作，使用ORCA量子化学程序完成，最后使用Molclus统计了该多元醇酯在298.15 K下构象的的玻尔兹曼分布，并通过Laplace键级指数排序，根据键的类型挑选了最可能断键的位置并作了键能大小的分析。
多元醇酯1-(癸酰氧基)-3-(壬酰氧基)丙-2-基十二烷酸酯的介绍
1-(癸酰氧基)-3-(壬酰氧基)丙-2-基十二烷酸酯是一种合成润滑油用的基础原料之一，该物质容易在$300$ $^oC$高温下发生热解断键而失效。所以这里通过键能大小来初步分析一下该多元醇酯容易从哪哥地方进行断裂，进而判断其热解断键的初步机理。
不会绘画该分子结构可以看如何用Avogadro画多元醇酯结构以及如何用GaussView构建多元醇酯结构。
构象搜索流程
本示例采用的软件及版本如下：
- Molclus 1.9.9.9 # http://www.keinsci.com/research/molclus.html
- xTB 6.3.2 # https://github.com/grimme-lab/xtb/releases
- ORCA 5.0.2 # https://orcaforum.kofo.mpg.de/
- Multiwfn 3.8_dev MAR-2-2022 # http://sobereva.com/multiwfn/
构想搜索流程概览
本示例参考sob老师对瑞德西韦药物的构象搜索案例。构想搜索流程图如下所示，主要分为4个步骤： - 获取>100 ps的动力学轨迹
- 初筛构象
- 采用稍高精度方法筛选构象
- DFT方法筛选构象
不会调用程序的请收看如何在Linux操作系统上做计算。
使用xTB获取100 ps的动力学轨迹
使用xTB获取分子动力学轨迹是十分方便且便宜的，在按上述说法建立好模型以后，将其命名为01.xyz，分子动力学控制输入文件md.inp参数如下：
1 | $md |
解读一下就是，该分子会在$300$ $^oC$高温下发生热解断键，这里设置动力学温度为$400$ $K$，高一点的温度可以帮助分子在跑动力学时跨过一些较高势垒，以获取势能面较为全面的采样。动力学的时间为100 ps，按sob老师说法是轨迹至少100 ps，不然可能采不够具有统计意义的点做后续的构象搜索，有条件的话可以跑200 ps，一般来说足够了，这个示例采用的100 ps，也够。dump设置为50 fs，step设置1 fs，就是说轨迹会跑100 ps，但记录数据的话只有$100 (ps)\div 50 (fs)=2000$步的轨迹。
随后直接调用xTB，使用GFN0-xTB方法获取100 ps的分子动力学轨迹：
1 | xtb 01.xyz --omd --gfn 0 -I md.inp |
这里是使用12个核做xTB的分子动力学，整个过程大致一个半小时左右。运算结束后，当前目录下的xtb.trj就是我们要的100 ps的分子动力学轨迹。
使用GFN0-xTB做初步构象筛选
然后解压下载好的Molclus，在Molclus的目录下修改settings.ini，这里初步筛选使用GFN0-xTB，故修改其中的：
- iprog= 4 # 使用xTB程序
- ngeom= 0 # 处理所有的轨迹
- itask= 0 # 只做结构优化
- xtb_arg= “–gfn 0 –chrg 0 –uhf 0” # 设定xTB程序使用采用GFN0-xTB方法
- 其余的使用默认设置，具体的可以读每一项的注释说明，依次做自己额外需要的修改
再将100 ps的分子动力学轨迹文件复制到Molclus的目录下，并改名为traj.xyz，这个文件是molclus程序默认读取的轨迹文件。随后直接调用Molclus进行初步的构象筛选的批处理结构优化：
1 | ./molclus |
如果发现解压后无法执行，请添加执行权限：
1 | chmod +x ./molclus |
添加完执行权限后再运行molclus即可。优化速度比较快约1s一个结构，2000个结构就是约40 min跑完，当然你可以分成两段跑，充分利用服务器的资源，只需要对应修改ngeom参数为对应的0-1000以及1001-2000即可，这里省事就只接用一个任务跑了。
跑完初步筛选后，筛选后得到的信息molclus会帮我们存在isomers.xyz文件中，如果担心有问题，或者后续想查看文件，届时可以用cp命令备份一份（这里只演示一次，后续备份的操作一样）：
1 | cp isomers.xyz isomers.xyz.gfn0 |
随后运行Molclus文件夹里的isostat来对isomers.xyz里的各个优化结构做初步筛选：
1 | ./isostat |
如果发现解压后无法执行，请添加执行权限：
1 | chmod +x ./isostat |
添加完执行权限后再运行isostat即可。运行isostat，交互界面会提示我们：
1 | Number of CPU cores to be used: 24 |
即询问我们当前需要用几个核进行isomers.xyz文件的读取并处理，这里我用的是我24核的服务器，这里按下回车就会直接使用24个核做这个任务，一般来说不需要那么多，2-4个核足矣，届时大家在用的时候看看运行isostat时界面上出现的信息即可。这里我输入一个4再回车，然后来到下一个交互界面：
1 | Input the energy threshold for distinguishing different clusters |
这里是让我们填入判断不同结构间的能量阈值，直接回车的话默认是$0.25$ $kcal/mol$，默认的阈值比较小，获取的结构会多很多，这里按流程图里的设置，输入0.5再回车。随后来到下一个交互界面：
1 | Input the geometry threshold for distinguishing different clusters |
这里是让我们填入判断不同结构间的几何阈值，直接回车用的是默认的$0.1$ $\mathrm{\mathring{A}}$，也是比较小的，获取的结构会比较多，这里也是用示意图里的输入0.5再回车。之后就开始做筛选了。
以上的输入界面在后边不再重新演示，一律按示意图里的参数设置。
我这里摘取了我们需要关心的输出结果：
1 | ...... |
这里最后问我们做不做构象的在特定温度下的玻尔兹曼分布，这里不做，直接回车就退出了。
可以看到，2000个结构经过GFN0-xTB初步筛选后，得到了381个结构，筛选掉了很多，这里筛选得到的结构全部写入了Molclus目录下的cluster.xyz文件内。注意之前设置的两个阈值是当新结构与上一个结构之间的能量及结构差异有一项在阈值内就会被归为新的结构，所以会得到与能量最低点结构差异如此多的构象。但是GFN0-xTB方法太粗糙了，这里就需要换用更高精度一些但不贵的方法来做进一步的筛选，以减小最后用精确的DFT方法做筛选的工作量。
使用GFN2-xTB做更高精度的构象筛选
进一步筛选的方法采用GFN2-xTB，这时只需要将之前修改的settings.ini中的xtb_arg项即可：
- xtb_arg= “–gfn 2 –chrg 0 –uhf 0” # 设定xTB程序使用采用GFN2-xTB方法
随后将筛选后的结构重命名成traj.xyz：
1 | mv cluster.xyz traj.xyz |
然后运行molclus进行进一步的构象结构优化批处理：
1 | ./molclus |
因为采用的是GFN2-xTB，精度高了一些，所以耗时也上来了，优化一个结构2-8 s不等，总体跑完下来和初步筛选的一样，约40 min跑完。备份一份isomers.xyz为isomers.xyz.gfn2
之后同样在Molclus目录下执行isostat进行筛选，结果如下：
1 | ...... |
使用DFT方法进行构象筛选
进一步筛选后留下了237个构象，但这对DFT做结构筛选还是太多了，这里选取能量较低（通常来说比较稳定）的构象，截取$\Delta E < 3$ $kcal/mol$的结构做DFT的构象筛选，但如果说计算资源足够的话，可以保留$\Delta E < 5$ $kcal/mol$的结构，这样来说统计意义也更强一些，也能尽可能的避免构象的遗漏。但限于时间这里只用$\Delta E < 3$ $kcal/mol$的结构，即前25个结构。所以再对settings.ini做修改，这次改用ORCA程序做DFT方法的构象筛选：
- iprog= 3 # 使用ORCA程序
- ngeom= 25 # 处理前25个结构，即$\Delta E < 3$ $kcal/mol$的结构
- itask= 3 # 做结构优化+振动分析，即用吉布斯自由能来判定能量差距
- ibkout= 1 # 这里保存ORCA做完几何优化和振动分析的输出文件（.out），建议做DFT筛选都把输出文件保存一下，一是为了检查，二是可以保留更多有价值的信息
- orca_path= “/cal_soft/orca502/orca” # 这里一定是写入ORCA程序的绝对路径
- ibkgbw=1 # 这里也保存一下ORCA做完几何优化和振动分析的记录波函数的文件（.gbw），同样建议把这个波函数文件保存一下，也是比较有价值的文件，后期可以针对不同的构象利用Multiwfn进行详细的波函数分析
因为光采用示意图中的B3LYP(D3-BJ)泛函和def2-SVP基组做的振动分析与单点计算获取的吉布斯自由能偏差较大，所以还需额外利用高精度的双杂化泛函和高档的基组来获取高精度的电子能量，以获取校正后的吉布斯自由能，这样筛选的结果也更为可靠。按sob老师的案例，只要在Molclus目录下新建一个template_SP.inp，即可在结构优化与振动分析结束后接着计算高精度的单点能并获取校正后的吉布斯自由能。（如果是Gaussian用户好像是新建一个template_SP.gjf文件）
这里直接复制一份目录中已有的ORCA模板输入文件template.inp为template_SP.inp即可：
1 | cp template.inp template_SP.inp |
然后修改template.inp的内容为示意图中的参数所示：
1 | ! B3LYP D3 def2-SVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop |
修改template_SP.inp的内容为示意图中的参数所示：
1 | ! PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX tightSCF noautostart miniprint nopop |
然后将之前用GFN2-xTB筛选后的结构文件cluster.xyz重命名为traj.xyz。就可以调用molclus进行DFT方法的构象结构优化与振动分析的批处理了，这里因为采用DFT方法，可以保存一下molclus打印在终端上的信息：
1 | ./molclus | tee out.txt |
这样打印出来的信息就会自动保存到当前目录下的out.txt文件中，方便后期查看。这里我也截取了molclus打印的一些信息给大家看：
1 | ...... |
可以看到molclus按照预定的只对237个构象的前25个构象进行了DFT方法筛选，先做的结构优化+振动分析任务，随后会输出校正的吉布斯自由能，如果振动分析没有虚频就会输出计算正常结束，若有虚频出现，则会有Warning提示有几个虚频存在。随后会自动调用template_SP.inp计算高精度的电子能量，并与之前获取的校正的吉布斯自由能进行相加获取体系校正后的吉布斯自由能。
一般有虚频的结构，体系能量都会高不少，从计算的角度看不是极小点结构，也不会是我们需要找的构象，忽略或者删掉即可。
整个任务在56核的服务器上算完花了50个小时，内存占用最多的是双杂化泛函计算的部分，需要近90 GB，我的24核56GB内存的服务器还是做不了这么高精度的计算，降低一档采用$\omega$B97M-V泛函配上def2-TZVVP基组或许可以完成这个任务，但时间也相应会多近30 h。平均一个结构优化+振动分析40-60 min，高精度电子能量计算需要60 min，平均下来一个结构需要2 h的计算时间。由此可见DFT方法做筛选需要的时间较长，特别是对稍大的体系而言，这还是开了RI的ORCA，如若用Gaussian，可能时间还要多上几十个小时。所以做构象搜索的时候也要充分考虑服务器的算力和资源问题。
可以用grep命令筛选一下有虚频的任务：
1 | grep "Nimag = 1" isomers.xyz |
备份好isomers.xyz以后，将有虚频的任务结构删掉，最终剩22个结构，这时候再使用isostat按照流程图里的设置进行构象筛选：
1 | ...... |
最终剩下十六个结构，最后再在isostat的交互界面输入298.15，在室温$298.15$ $K$下进行构象的玻尔兹曼分布计算：
1 | ...... |
可以看到，大约有4个构象，最稳定的两个构象占比都很大。可以观察到能量最低的#1构象中，三个碳链朝同一方向拧转，能量第二低的#2构象则是三个碳链朝不同的方向伸展，以上两种构象占该化合物构象分布的绝大部分，而剩下少部分的#3与#4构象中均只有两个碳链朝同一方向拧转。
多元醇酯的断键分析
随后对最稳定的两个构象进行拉普拉斯键级计算，按sob老师博文所述（Multiwfn支持的分析化学键的方法一览），拉普拉斯键级指数可以估计同类型键的强度和极性大小。这里也应用一下拉普拉斯键级来分析该多元醇酯在那个地方键能相对较低，即预测容易断键的位置。这里使用之前高精度计算电子能量任务中得到的波函数文件进行分析。这里以#1号构象的波函数分析为例，直接调用Multiwfn软件读取波函数文件（比如是01.molden.input）并将结果输出到01_wfn.out中：
1 | Multiwfn 01.molden.input | tee 01_wfn.out |
在交互界面依次输入：
8 # 使用键级分析功能
8 # 计算拉普拉斯键级
q # 退出Multiwfn
我将分析的结果直接整理成了表格：
可以看到有三种键容易断裂，分别是C-O键、酯键、酯键上的$\alpha$-C与$\beta$-C的C-C键。该多元醇酯的断键标识也绘制了一番，如下图所示：
所以初步判断是三个类型的键会断裂，并且每个类型的键也必有三处断裂的位置对应该多元醇酯的三条尾巴，所以分别断裂这些位置的键并分别建模获取自选多重度为2，带电荷数为0的模型，共十八个模型，用于后续的键能计算，键能计算直接采用片段的电子能量总和减去整体化合物的电子能量即可，公式如下：
$$BE=E_{el,1}+E_{el,2}-E_{el,0}$$
把所有的键能计算任务的输入文件制作完成并上传到同一个文件夹内后，直接运行自动进行队列计算的脚本batch.sh。脚本会自动将目录下所有带.inp后缀的的文件进行数量统计，并调用orca进行计算。这里需要注意，ORCA 5版本默认使用TRAH进行收敛计算，但该体系在使用双杂化泛函PWPB95下使用4-zeta的def2-QZVVP基组时，20步内容易不收敛，届时可以在关键词行加入noTRAH来关闭TRAH的使用，保证计算流程的顺利进行。该脚本的的内容如下：
1 | #!/bin/bash |
最后汇总得到的#1构象断键分析总表如下：
可以看到除了C-C键外是10个碳链处的容易断裂外，其余的C-O键断裂均是12个碳链侧容易断裂，这与拉普拉斯键级指数大小对应较好。
总结
本示例首先利用了xTB程序进行分子动力学模拟，以对多元醇酯的势能面进行充分的采样。随后依次采用GFN0-xTB、GFN2-xTB、DFT方法在真空下对多元醇酯1-(癸酰氧基)-3-(壬酰氧基)丙-2-基十二烷酸酯进行了逐级构象筛选并在室温$298.15$ $K$下做了构象的玻尔兹曼分布，最终获得了4个较为稳定的构象。之后对最稳定的构象进行了拉普拉斯键级指数分析，以初步预测该多元醇酯在拥有的键类型中，容易断键的位置，最后介绍了如何使用自动队列计算脚本进行多个任务的计算，获取不同断键位置的键能。
通过这个示例也展示了柔性大分子有构象搜索的必要性。大家也可以尝试做些简单的构象搜索比如：
- 环己烷的两种构象（桥式和椅式）
- 苄基苯基醚的构象
这样可以深入的认识并掌握构象搜索的流程。
附
以下是我对苄基苯基醚在$\Delta E < 1$ $kcal/mol$，$\Delta G_{min} < 0.5$ $\mathrm{\mathring{A}}$的构象筛选，并在$298.15$ $K$下获取的构象玻尔兹曼分布结果（其中4号构象是存在一个虚频的，这里忽略就好）：
1 | ...... |
