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

Stardust0831

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

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

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

zxn1318

老师您好，您找过渡态的思路是不是先做柔性扫描，找到极值，以此为初始结构，做过渡态优化呢，但我看到您的过渡态结构优化有一个条件没有收敛，这样也可以么，谢谢老师

Stardust0831

zxn1318 发表于 2025-8-19 15:46
老师您好，您找过渡态的思路是不是先做柔性扫描，找到极值，以此为初始结构，做过渡态优化呢，但我看到您的 ...

可以，键长过近/过远的情况是不被关心的，感兴趣的键长的部分收敛就可以。扫描的目的是为了拿到更精确的TS对应的键长+弛豫结构。

如果提前知道TS对应的键长大约是多少，也可以约束那个键的键长为猜的TS键长，先跑限制性优化，然后解除冻结跑TS。

当前结构不太好猜TS对应的键长是多长，所以需要柔性扫描扫一下断键过程，扫出第一个极大点就行。

luoxuping

你使用#p opt=(modredundant,nomicro) external='./xtb.sh'  做柔性扫描的时候，输出的.log文件Scan of Total Energy的能量在.log文件哪个位置的？

Hugo_314cat

或许可以尝试拓展一下，让它能够去做过渡态的构象搜索？
以我的经验来看，xtb2是完全可以胜任这种任务的，也成功帮我做了大量的过渡态构象搜索工作。如果能用简单脚本就可以实现的话就更好了。

luoxuping

Hugo_314cat 发表于 2025-8-27 13:38
或许可以尝试拓展一下，让它能够去做过渡态的构象搜索？
以我的经验来看，xtb2是完全可以胜任这种任务的， ...

这个怎么做过渡态搜索，你是怎么做构象搜索的

Daniel_Arndt

luoxuping 发表于 2025-8-27 16:57
这个怎么做过渡态搜索，你是怎么做构象搜索的

实际操作中是给研究者关心的若干个键长（往往就是成键、断键对应的键）设置constraining potential，就是在运行crest的目录增加一个constrain.inp文件，此文件的格式可以参考 https://xtb-docs.readthedocs.io/ ... training-potentials 和运行“man 7 xcontrol”后的输出（得先加载GFN-xTB相应的module）中相关内容。然后在crest的命令中添加一项“--cinp constrain.inp”（参考 https://crest-lab.github.io/cres ... onstraining-options ）。

kikiw

非常感谢分享，困扰我几天的过渡态终于有指导方法了

nanokey

老师您好，我现在想做一个柔性扫描看单电子转移过程中C-F断裂过程，期望是看到F-远去，但是我使用柔性扫描得不到预想的结果，想请您看看有什么问题

%chk=/home/PTFE6-scan.chk
#p opt=modredundant def2svp m062x scrf(SMD,solvent=THF)

Title Card Required

-1 2
C                 -0.64131100    0.16679800   -0.29131000
C                 -1.89641800   -0.60515100   -0.09766500
F                 -0.60613900    0.85888900   -1.43669300
F                 -1.77180700   -1.42890300    0.95869600
F                 -2.15100900   -1.38673400   -1.19158400
F                 -0.79474941    1.19469451    0.53316133
C                  0.65264900   -0.43138900    0.12432800
C                 -3.19275600    0.22999100    0.11668300
F                  0.79022500   -1.72532300   -0.32034800
F                  0.77550800   -0.50895500    1.46576000
F                 -4.25329500   -0.54525500   -0.13765500
F                 -3.23658400    1.26062700   -0.71491800
F                 -3.30826900    0.65202000    1.36045600
C                  1.87855600    0.34082900   -0.43720200
F                  1.73346900    1.65706100   -0.28502500
F                  2.00109700    0.06584100   -1.74441800
C                  3.22409600   -0.04304700    0.24052400
F                  3.38901500   -1.35739900    0.27151400
F                  3.28426700    0.42956700    1.47440300
F                  4.21695300    0.49092200   -0.46152500

1 6 S 20 0.05

Stardust0831

nanokey 发表于 2026-1-14 10:04
老师您好，我现在想做一个柔性扫描看单电子转移过程中C-F断裂过程，期望是看到F-远去，但是我使用柔性扫描 ...

尝试通过广义化内坐标(GIC)进行柔性扫描，把F拉离这个分子的几何中心。详见《详谈使用Gaussian做势能面扫描》
http://bbs.keinsci.com/forum.php ... 0&fromuid=58653
(出处: 计算化学公社)
扫描前确保电荷和自旋多重度正确，我个人不推荐M06-2X+SMD用于柔性扫描和过渡态搜索等任务理由是格点精度要求高、计算耗时高、易出虚频。

nanokey

Stardust0831 发表于 2026-1-14 15:46
尝试通过广义化内坐标(GIC)进行柔性扫描，把F拉离这个分子的几何中心。详见《详谈使用Gaussian做势能面扫 ...

好的，感谢老师！