计算化学公社

标题: 过渡态是自旋极化单重态时找TS和跑IRC流程分享（以坛友的实际问题为例） [打印本页]
作者
Author: Stardust0831    时间: 2025-3-26 02:25
标题: 过渡态是自旋极化单重态时找TS和跑IRC流程分享（以坛友的实际问题为例）
本帖最后由 Stardust0831 于 2025-3-26 02:38 编辑

昨晚，翻论坛时看见坛友@Solitude198 在帖子《给大家代算（不收费）》中发布的代算请求：
(, 下载次数 Times of downloads: 48)
我以此为例，分享一下过渡态是自旋极化单重态时找TS和跑IRC的流程。计算软件为 Gaussian 16:  ES64L-G16RevA.03。
本轮计算的完整输入输出文件链接（含chk文件）: https://pan.baidu.com/s/1pm2RibhDDQ7CwbJL6esX0g?pwd=ce8q 提取码: ce8q
本轮计算的gjf和out文件： (, 下载次数 Times of downloads: 88)

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

以下为计算流程：
首先运行ts_transform-guess.gjf产生稳定的波函数：
  1. %chk=ts_transform_guess.chk
  2. #p stable=opt 6-31g(d,p) em=gd3bj upbe1pbe guess=mix nosymm

  4. Generated by Multiwfn

  6. 0 1
  7. 坐标略
复制代码
然后运行ts_transform.gjf读取上一步产生的稳定的波函数搜索TS：
  1. %oldchk=ts_transform_guess.chk
  2. %chk=ts_transform.chk
  3. #P PBE1PBE/6-31G** em=gd3bj opt=(calcfc,ts,noeigen,nomicro,maxstep=5) freq guess=read Geom=AllCheck

  5. Generated by Multiwfn

复制代码
检查虚频无误后，运行IRC-forward_transform.gjf和IRC-Reverse_transform.gjf，读取TS任务的波函数和hessian矩阵分别跑两个方向的IRC（此时需要读取初猜波函数，所以两个方向必须分开跑）：
  1. %oldchk=ts_transform.chk
  2. %chk=IRC-forward_transform.chk
  3. #P PBE1PBE/6-31G**  em=gd3bj IRC(forward,rcfc,Maxpoints=200,LQA) guess=read Geom=AllCheck

  5. Generated by Multiwfn

复制代码
  1. %oldchk=ts_transform.chk
  2. %chk=IRC-Reverse_transform.chk
  3. #P PBE1PBE/6-31G** em=gd3bj IRC(Reverse,rcfc,Maxpoints=200,LQA) guess=read Geom=AllCheck

  5. Generated by Multiwfn

复制代码
最后取出两个方向的IRC的最后一个点，进行结构优化以得到反应物和生成物结构：
  1. %oldchk=IRC-Reverse_transform.chk
  2. %chk=opt-Reverse_transform.chk
  3. #P PBE1PBE/6-31G** em=gd3bj opt freq guess=read Geom=AllCheck

  5. Generated by Multiwfn
复制代码
  1. %oldchk=IRC-forward_transform.chk
  2. %chk=opt-forward_transform.chk
  3. #P PBE1PBE/6-31G** em=gd3bj opt freq guess=read Geom=AllCheck

  5. Generated by Multiwfn
复制代码

最后，可以在得到的结构的基础上用更高的计算级别算单点，结合shermo用RRHO模型计算反应物、生成物和过渡态的自由能，并基于自由能绘制能量台阶图即可。此处必须用RRHO模型而不能用准RRHO模型的原因是过渡态的虚频就是需要被投影掉的，此时准RRHO对虚频的处理就明显不合适了。

笔者能力有限，如有错漏，欢迎指正。




作者
Author: Solitude198    时间: 2025-3-26 13:25
首先运行ts_transform-guess.gjf产生稳定的波函数：


我觉得这是最关键一步，之前我在尝试的时候只想到用guess=mix，没有stable=opt，导致永远没能确定波函数稳定，运气再怎么好也不容易收敛到双自由基的情况，后面查看 S**2 也能发现自旋污染偏差到了1，也是符合预期了

264这篇基本可以解惑了我觉得 http://sobereva.com/264

btw, lz   nb
作者
Author: sobereva    时间: 2025-3-26 13:28
"此时准RRHO对虚频的处理就明显不合适" 这个描述有问题。如果这里的“准RRHO”是Shermo支持的比如Grimme、Minenkov、Truhlar的QRRHO，对于有虚频的情况是没有问题的，也是和RRHO一样直接忽略掉虚频对热力学量的贡献的。只不过对于过渡态，一般应该保持Shermo的settings.ini里的imagreal为默认的0，而如果设成比如70，则大小小于70 cm^-1的虚频会当成相应的实频处理。
作者
Author: Stardust0831    时间: 2025-3-26 13:44
sobereva 发表于 2025-3-26 13:28
"此时准RRHO对虚频的处理就明显不合适" 这个描述有问题。如果这里的“准RRHO”是Shermo支持的比如Grimme、M ...

学到了，感谢社长的指正！
作者
Author: Uus/pMeC6H4-/キ    时间: 2025-3-27 12:08
本帖最后由 Uus/pMeC6H4-/キ 于 2025-3-27 14:31 编辑

说起来，如果已知一步基元反应的反应物和产物，如何从一开始来判断按自旋极化/对称破缺单重态找过渡态的必要性呢？有的过渡态对称性高于反应物和产物的情况（如http://bbs.keinsci.com/thread-33890-1-1.html讲的经典环丁二烯异构化）或许可以预料到，以及有的按Woodward–Hoffmann规则属于concerted热禁阻的周环反应可以考虑走asynchronous双自由基机理（Stevens重排这种构型保持的[1,2]迁移似乎可以算一个例子，注意concerted和asynchronous都是描述单步的用词，有别于多步的stepwise），但其他情况呢？

考虑基态闭壳层分子A向闭壳层分子B的转化，如果单重态下能同时找到闭壳层过渡态TS和自旋极化的开壳层过渡态TS'（编辑：漏了一个前提，TS和TS'各自的结构下波函数均稳定且单一虚频振动方向均正确），且TS'跑IRC'全程能量低于TS跑的IRC，那确实可以下结论。但对一楼这个体系，从之前的一帖来看是找不到闭壳层的情况才开始考虑开壳层的，这确实也成为合适的理由但说服力感觉稍弱。
作者
Author: Stardust0831    时间: 2025-3-27 13:43
Uus/pMeC6H4-/キ 发表于 2025-3-27 12:08
说起来，如果已知一步基元反应的反应物和产物，如何从一开始来判断按自旋极化/对称破缺单重态找过渡态的必 ...

我实际计算中有两个细节。1、我并没有假设当前体系一定是自旋极化单重态。我实际计算中是只用了stable=opt，而没有guess=mix。这样第一次收敛到的波函数就是闭壳层了。只有当闭壳层波函数不稳定时才会尝试去收敛到对称破缺单重态。所以最终的稳定波函数既然是自旋极化单重态，那一定意味着当前结构为闭壳层时波函数不稳定。
2、我找到过渡态以后又判断了一次波函数稳定性。这样即使在最开始错把自旋极化单重态当成闭壳层算、也能在此时发现问题。
作者
Author: Uus/pMeC6H4-/キ    时间: 2025-3-27 15:09
Stardust0831 发表于 2025-3-27 13:43
我实际计算中有两个细节。1、我并没有假设当前体系一定是自旋极化单重态。我实际计算中是只用了stable=op ...

前帖有个地方没写清楚，刚才编辑了下。可能我理解有误，换个问法说吧：是不是因为都在单重态势能面上，所以连接A和B的（波函数稳定且单一虚频振动方向正确的）基元反应过渡态TS结构唯一，具体是闭壳层还是自旋极化开壳层也就确定了？上面说的“能同时找到”两种结构和自旋状态均相异的过渡态的情况能存在么？
作者
Author: RAL    时间: 2025-3-27 15:12
guess=read读取波函数时，如果是开壳层波函数，pbe1pbe前面不加u也会自动切换到UKS吗
作者
Author: Stardust0831    时间: 2025-3-27 16:37
RAL 发表于 2025-3-27 15:12
guess=read读取波函数时，如果是开壳层波函数，pbe1pbe前面不加u也会自动切换到UKS吗

chk里如果是UKS波函数，会自动当成UKS计算
作者
Author: TIANTTT    时间: 2025-4-4 20:28
本帖最后由 TIANTTT 于 2025-4-4 21:19 编辑

感谢老师指点！
作者
Author: Stardust0831    时间: 2025-4-4 20:37
本帖最后由 Stardust0831 于 2025-4-4 20:43 编辑
TIANTTT 发表于 2025-4-4 20:28
老师您好！我是计算化学新手，刷到您的帖子觉得您的讲解非常详细，冒昧前来求助。我目前在算一个55个原子的 ...

贴一下反应的方程式。我来看看。
IRC只跑两步就结束可能是因为势能面过于平缓，或者TS不够精确。
作者
Author: sobereva    时间: 2025-4-4 20:48
TIANTTT 发表于 2025-4-4 20:28
老师您好！我是计算化学新手，刷到您的帖子觉得您的讲解非常详细，冒昧前来求助。我目前在算一个55个原子的 ...

同一个问题不要在论坛里的两处同时问，否则会严重扰乱讨论秩序。http://bbs.keinsci.com/thread-52652-1-1.html里已经问过一次了
作者
Author: TIANTTT    时间: 2025-4-4 21:20
sobereva 发表于 2025-4-4 20:48
同一个问题不要在论坛里的两处同时问，否则会严重扰乱讨论秩序。http://bbs.keinsci.com/thread-52652-1- ...

感谢提醒！已做修改
作者
Author: xiachenche    时间: 2025-5-13 15:53
你好，我想请问一下”最后，可以在得到的结构的基础上用更高的计算级别算单点",  在用更高计算级别算过渡态的单点时，还需要确保过渡态是自旋极化单重态吗？
作者
Author: Stardust0831    时间: 2025-5-13 16:23
xiachenche 发表于 2025-5-13 15:53
你好，我想请问一下”最后，可以在得到的结构的基础上用更高的计算级别算单点",  在用更高计算级别算过渡态 ...

显然需要。
可以用stable=opt做波函数稳定性测试。
作者
Author: ionexchangeC    时间: 2025-5-14 02:39
最近正好在研究这类问题。
我的解决方法是，用一个Gaussian的external脚本，做到在结构优化/寻找过渡态的过程中，每一步都附加波函数稳定性检测。
(, 下载次数 Times of downloads: 6) (, 下载次数 Times of downloads: 5) (, 下载次数 Times of downloads: 5)
作者
Author: xiachenche    时间: 2025-5-14 14:29
Stardust0831 发表于 2025-5-13 16:23
显然需要。
可以用stable=opt做波函数稳定性测试。

好滴，感谢老师答复！
作者
Author: EddieSu    时间: 2025-8-25 19:24
老师您好，这里寻找过渡态时opt里写nomicro的目的是什么呢？
作者
Author: Stardust0831    时间: 2025-8-25 20:36
EddieSu 发表于 2025-8-25 19:24
老师您好，这里寻找过渡态时opt里写nomicro的目的是什么呢？

NoMicro forbids microiterations during ONIOM(MO:MM) optimizations.
我自己的工作流会用到Gaussian和xtb接口，需要用这个关键词来避免用到分子力场的优化器。
此处是转dft级别的时候忘删了，是个多余的关键词。
感谢指正。
作者
Author: 云边的小卖铺    时间: 2025-11-14 22:43
老师您好，您的帖子很详细。针对这个我还是有些疑惑，IRC两侧直接进行opt的问题，对于单侧原料端不是自旋极化单重态吗？所以可以直接进行对最低点进行几何优化和频率分析吗？
作者
Author: Stardust0831    时间: 2025-11-15 16:04
云边的小卖铺 发表于 2025-11-14 22:43
老师您好，您的帖子很详细。针对这个我还是有些疑惑，IRC两侧直接进行opt的问题，对于单侧原料端不是自旋极 ...

是否是自旋极化单重态以stable=opt的结果为准，不放心的话可以多stable=opt几次。
帖中直接读取了IRC的最后一个点的波函数当初猜，如果在IRC的时候维持住了自旋极化单重态，那此时opt的结果也是自旋极化单重态的。
作者
Author: 云边的小卖铺    时间: 2025-11-16 10:19
Stardust0831 发表于 2025-11-15 16:04
是否是自旋极化单重态以stable=opt的结果为准，不放心的话可以多stable=opt几次。
帖中直接读取了IRC的 ...

感谢您的回复，最近我尝试了一个自旋极化单重态的计算。找到TS，IRC趋势正确（分别跑两端，和您帖子过程一致），但在对原料优化过程中，有些问题。我第一次IRC走了80步，但结果很陡峭，又算了一次走了150步，仍然较陡峭，此时我读取IRC最后一点波函数初猜，进行opt，过程中优化的方向和我预期不符（两个独立片段中部分原子成键了），1.我猜测是否是这个opt过程中自旋极化单重态未维持住；2.有没有可能是TS没找对？（TS一个虚频，目前IRC右侧产物端opt已得到正确结果）3.此外我仍做了一个更长的IRC，走了250步，IRC结果趋势仍然没有缓和，较为陡峭，opt失败是否和这个有关系呢？（仍陡峭则最后一点与实际构型相差过大而导致难收敛）还请您指教。十分感谢！
作者
Author: Stardust0831    时间: 2025-11-16 12:30
云边的小卖铺 发表于 2025-11-16 10:19
感谢您的回复，最近我尝试了一个自旋极化单重态的计算。找到TS，IRC趋势正确（分别跑两端，和您帖子过程 ...

1.可以取出opt任务对应的chk文件，转成fchk等格式后用Multiwfn查看自旋密度来确定是不是仍然是自旋极化单重态。可以读chk中的波函数做stable=opt计算来确认波函数稳定性。
2.TS正确性需要结合具体反应和具体的计算输出文件讨论。
3.IRC的时候可以选择更小的步长（stepsize），或者删掉LQA关键词用默认的HPC（但用HPC可能遇到矫正步不收敛）。这俩方法都能得到更平滑的IRC。
作者
Author: 云边的小卖铺    时间: 2025-11-24 09:31
Stardust0831 发表于 2025-11-16 12:30
1.可以取出opt任务对应的chk文件，转成fchk等格式后用Multiwfn查看自旋密度来确定是不是仍然是自旋极化单 ...

好的，谢谢老师。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3