使用柔性扫描与半经验DFT辅助搜索有机反应过渡态经验分享

Stardust0831

本帖为我之前的《过渡态是自旋极化单重态时找TS和跑IRC流程分享（以坛友的实际问题为例）》延伸与补充。
在看论坛的时候，发现坛友@LBH0532 在帖子《过渡态优化时六元糖环意外开环》中尝试重复的体系很有代表性，适合当作例子来介绍过渡态搜索中的柔性扫描与半经验DFT技巧。故单开一个帖子，用于分享这些计算经验。该机理较普通且已被发表过，原帖也是公开问的，故不再考虑保密问题。
本轮计算的gjf和out文件： TS1-LBH0532.zip (8.86 MB, 下载次数 Times of downloads: 4)

2025-8-1 02:52 上传 Uploaded

点击下载
Click to download

如果大家有类似的计算需求，也可以在本贴留言，我会在精力允许的情况下尽量帮忙解决大家计算中遇到的问题。（不收费，给点eV就行）（需要提前确认好自己有对应软件的版权）
如果想系统学习高斯搜索过渡态的技巧，欢迎大家报名：北京科音中级量子化学培训班。sob老师在这个培训班中有非常细致的过渡态搜索技巧的讲解。

待计算的反应是类似文献https://link.springer.com/article/10.1007/s00214-020-02693-x 中的葡萄糖衍生物叠氮炔click反应，涉及到了一根碳氮键的形成。

对于这个体系，结构有一定复杂性，不好确定opt=TS找过渡态的初猜结构，此时借助柔性扫描帮助判断合适的初猜结构往往是有帮助的。由于柔性扫描过程本质上相当于一大批限制性优化，计算成本相当高昂，因此本例使用了社长开发的gau_xtb接口，实现Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析，通过比DFT计算便宜2-3个数量级的计算成本，即可完成柔性扫描到TS和IRC的初期摸索。

1.结构整理
初始结构来自手动建模，可能存在诸多不合理，故调整到合适的键级，使用gview内置的clean快速、粗糙的调整一次结构。并手动微调结构至合理。对于opt=TS，在手动微调的时候，可以将要断裂/形成的键调整至成键键长的1.5倍左右，本例是柔性扫描，保持这个C-N的距离较长即可。


2.柔性扫描
初始结构中，N和碳键距离较远，考虑由远及近扫描。首先找到希望通过柔性扫描实现靠近的两个原子。直接用测量工具选择他们，在gview窗口的左下可以看到他们是3号原子N和4号原子C：


随后，这样准备输入文件：

1. %chk=TS1-LBH0532-scan.chk
   
2. %nprocshared=1
   
3. #p opt=(modredundant,nomicro) external='./xtb.sh'
   
4. 
   
5. ****
   
6. 
   
7. 2 2
   
8. N                  0.50378900    0.12533200   -0.45515100
   
9. N                  0.40712500    1.25044600   -0.04678000
   
10. N                  0.09065500    2.26645400    0.29409200
    
11. C                 -2.74812400    1.86153300   -0.50654000
    
12. （中间坐标略）
    
13. H                  6.94015500    0.41090700    0.07463000
    
14. 
    
15. 4 3 S 17 -0.1

复制代码

需要注意的是，使用gau_xtb接口的时候必须使用nomicro关键词，否则Gaussian在优化的时候会试图调用分子力学的optimizer去搞，达不到我们的目的。写%nprocshared=1是为了实现串行计算，避免在Gaussian通过external方式调用外部脚本期间会造成很高无意义的CPU资源消耗。最后一行的“4 3 S 17 -0.1”意思就是对4号和3号原子的距离每次都缩小0.1A，一共缩小17次。
调用xTB计算很便宜，1分钟就拿到了结果。

可以发现曲线有一个很明显的极大点，这便是在N-C键形成过程中能量最高的位置，也是后续计算中TS的初猜。（在距离更近的地方能量再次上升，这是由于距离过近，原子间斥力导致的。）（不需要管此时的成键显示，原因见《谈谈原子间是否成键的判断问题》，gview给的成键不影响此时的计算、也没有意义）
3.半经验档次过渡态搜索
取上一步找到的能量极大点结构，保存成gjf文件，重写关键词：

1. %chk=TS1-LBH0532.chk
   
2. %nprocshared=1
   
3. #p opt=(calcfc,ts,noeigen,nomicro,maxstep=5) external='./xtb.sh'

复制代码

此时，使用maxstep=5是为了避免优化时的结构变化过于剧烈，导致结构不再接近过渡态。柔扫的极大点通常都挺接近TS的，maxstep=5不会给收敛带来太多计算成本负担。
计算很快收敛，再从%chk中读取结构跑freq（使用脚本时，TS和freq需要分开），发现虚频已经符合预期：

1. %oldchk=TS1-LBH0532.chk
   
2. %chk=TS1-LBH0532-freq.chk
   
3. %nprocshared=1
   
4. #p freq geom=allcheck external='./xtb.sh'
   

复制代码

4.DFT档次过渡态搜索

1. %oldchk=TS1-LBH0532-freq.chk
   
2. %chk=TS1-LBH0532-dft.chk
   
3. #p opt=(rcfc,ts,noeigen,maxstep=5) PBE1PBE def2svp freq em=gd3bj geom=allcheck
   

复制代码

此时，%oldchk写成半经验档次振动分析的chk文件，我们从其中用geom=allcheck读取结构、电荷和自旋多重度的设定。用rcfc是为了读取半经验的hessian结果，精度通常够用，且避免了算DFT级别hessian的昂贵耗时。过渡态搜索通常不需要大基组，此时选用关键词def2svp来使用def2-SVP基组。由于含有Cu元素，我倾向于直接用PBE0-D3(BJ)，这个泛函和B3LYP-D3(BJ)都是公认适合结构优化的泛函，优化精度轮流拿第一的程度，但比起B3LYP-D3(BJ)，有证据表明PBE0-D3(BJ)更适合这种含有过渡金属的体系。我也曾在2群请教过社长，有幸得到社长的答疑：

计算花了些时间，但我们半经验档次找到的TS的特征被很好的保留下来了，顺利收敛后，检查发现虚频也是正确的：


5.DFT档次IRC
由于机时问题，我没有真的跑这个，但不妨说一下关键词：

1. %oldchk=TS1-LBH0532-dft.chk
   
2. %chk=TS1-LBH0532-IRC.chk
   
3. #p irc=(rcfc,maxpoints=50,lqa) PBE1PBE def2svp em=gd3bj geom=allcheck

复制代码

rcfc和geom=allcheck是读取上一步的结构、hessian等参数，由于找TS的时候已经算过freq，我们直接读取、避免重算hessian。LQA会比默认HPC方法快一些，且不会出现校正步不收敛问题导致IRC断掉，代价是IRC曲线可能会更不光滑一些。此处只是为了确认反应物和生成物的正确性，也就无需讲究曲线的好看。算完IRC后，取出两个方向的IRC的最后一个点，进行opt+freq以得到反应物和生成物结构，freq用于使用shermo算热力学量。


总结
通过使用社长的gau_xtb接口，使用目前最好的半经验之一的GFN2-xTB来算有机反应机理可以在计算成本大大降低的情况下得到还不错的结构和能量结果，帮助我们快速弄清势能面。柔性扫描是这种成键/断键很重要的反应的过渡态中很重要的技巧，能辅助找到合适的过渡态结构。使用本文的思路，对于每个过渡态，只需要计算一次DFT档次的hessian，这对于一些很大的体系而言可以明显节约耗时。
希望本文能对你有帮助，笔者能力有限，如有错漏，欢迎指正。

loovfnd

有收获，有理有据，收藏了。