计算化学公社

标题: 小窍门：使用自旋极化单重态波函数作为NEB的初猜 [打印本页]
作者
Author: Freeman    时间: 2021-4-21 12:27
标题: 小窍门：使用自旋极化单重态波函数作为NEB的初猜
本帖最后由 Freeman 于 2021-4-22 11:23 编辑

ORCA支持的NEB法可以用来找能量最低反应路径、过渡态。一般来说，要用NEB的话，我们需要输入一个反应物坐标和一个生成物坐标，然后ORCA就会自动产生一系列中间点（image）及其对应的波函数初猜。但是如果这些中间点里有比较特殊的电子结构（例如是自旋极化单重态），NEB自动产生的波函数初猜经过SCF后可能会收敛到不稳定波函数上去，反应路径就不对了。我最近在这个地方踩过一些坑，如今终于折腾出来了，于是就写了这篇小窍门，希望能帮到遇到我这种情况的小伙伴儿。

假如我们现在有两个优化好的反应物和生成物坐标，reactant.xyz和product.xyz，其中reactant.xyz是自旋极化单重态。如果我们直接写neb1.inp
  1. ! neb ......
  2. ......
  3. %neb
  4.     neb_end_xyzfile "product.xyz" #给出生成物坐标
  5. end
  6. * xyzfile 0 1 reactant.xyz #给出反应物坐标
复制代码
就会导致（1）提示reactant不处在势能面极小点以及（2）反应路径完全错误，因为ORCA这里并不会自己找到正确的稳定波函数。我们要采取如下策略。

一、创建一系列中间点
ORCA的NEB可以自动创建中间点。我们先运行上面的neb1.inp，然后等到当前文件夹下出现neb1_initial_path_trj.xyz后杀死任务。这个neb1_initial_path_trj.xyz保存了ORCA/NEB通过差值产生的中间点。

二、获取中间点的稳定波函数
把每一个中间点保存为一个带波函数稳定性追踪的inp文件neb2_imN.inp（其中N为中间点的序号，反应物对应0，依次排号至生成物）。注意：ORCA/NEB产生的第一个和最后一个坐标分别为反应物和生成物，但该坐标和neb1.inp里的product.xyz与reactant.xyz可能采取不同的坐标系，最好统一以新产生的坐标为准。在这些inp的文件里实现波函数稳定性追踪的写法是neb2_imN.inp
  1. ! ......
  2. ......
  3. %scf
  4. HFTyp uhf #使用非限制性自洽场
  5. STABPerform true #做波函数稳定性检测
  6. STABRestartUHFifUnstable true #如果检测出波函数不稳定，那么再做一遍自洽场找稳定波函数（有极低可能最后仍然没找到，保险起见可以再创建一个inp文件以进一步检测波函数稳定性）
  7. STABDTol 0.00001 #设置波函数稳定性的判定阈值，见下面“其他细节3”
  8. STABRTol 0.00001
  9. end
  10. * xyz 0 1
  11. ......
  12. *
复制代码
运行后获得一系列稳定的波函数文件neb2_imN.gbw。

三、创建allxyz文件
与neb1_initial_path_trj.xyz一样，ORCA/NEB的allxyz文件保存了所有中间点坐标，但是前者是初始坐标，经过优化（不论是否收敛）后的坐标才会保存进allxyz。allxyz是自定义中间点/自定义初猜波函数NEB所必须的，因此我们需要把neb1_initial_path_trj.xyz里的信息复制到neb1_MEP.allxyz，后者将用于下一步的neb初猜。allxyz与xyz文件格式差不多，唯一的区别是两个中间点之间由>连接，如下：neb1_MEP.allxyz
  1. ......
  2. >
  3. 2
  4. image2
  5. H 0.0 0.0 0.0
  6. H 0.0 0.7 0.0
  7. >
  8. 2
  9. image3
  10. H 0.0 0.0 0.0
  11. H 0.0 0.8 0.0
  12. >
  13. ......
复制代码


四、创建并运行最终的neb任务
最终的neb任务inp文件写为neb3.inp
  1. ! neb ......
  2. %neb
  3. neb_restart_xyzfile "neb1_MEP.allxyz" #给出allxyz文件作为初始中间点，这些中间点包括了反应物和生成物，但是下面第10行还要写一遍反应物坐标
  4. neb_restart_gbwname "neb2" #给出波函数初猜文件，“neb2”表示波函数初猜文件为“neb2_imN.gbw”
  5. end
  6. %scf
  7. hftyp uhf #仍然要指定使用非限制性自洽场
  8. end
  9. *xyz 0 1
  10. [neb2_im0.inp对应的坐标]
  11. *
复制代码
运行后得到结果。


其它细节：
1、neb分为普通的neb、neb-ci和neb-ts，分别代表只做neb（只找最小反应路径）、做完neb再做climbing_image和做完climbing_image后再当作过渡态优化一遍。可想而知，neb-ts是最稳妥的方法。另外还有zoom-neb，意思是在最高的几个中间点之间再差值，做neb。以上流程适用于neb、neb-ci和neb-ts，但不保证适用于zoom-neb，因为后者会多出来一些中间点，而要给这些多出来的中间点找稳定波函数，操作就比较繁琐了。
2、ORCA的neb和几何优化的收敛限都高得过头，所以要手动调低收敛限。见orca_manual_4_2_1.pdf第635页。
3、ORCA的波函数稳定性检测判定阈值也是太窄，导致经常把不稳定波函数当成稳定波函数了，需要手动拓宽。见http://bbs.keinsci.com/thread-16217-1-1.html我在15楼的回复。
4、neb-ts可以和freq一起用，意思是找到过渡态后再算一下频率。注意ORCA不支持meta-gga杂化泛函的解析频率（@喵星大佬 指正），所以对于M06-2X这类泛函，要写作numfreq，意为使用数值频率计算。若要观看振动模式，可以使用社长的OfakeG工具，见http://bbs.keinsci.com/thread-13952-1-1.html。
5、以上流程中都没有写并行核数。用户要根据中间点的个数选定并行核数，见http://bbs.keinsci.com/thread-20961-1-1.html我在12楼的回复。

欢迎讨论。如有纰漏，欢迎批评指正。


作者
Author: 喵星大佬    时间: 2021-4-21 13:29
orca不支持的是meta泛函的解析Hessian，像B3LYP，PBE0之类的普通杂化泛函是可以的
作者
Author: biogon    时间: 2021-5-19 21:49
不过说实话能用高斯的opt=ts我肯定不会用这玩意
作者
Author: 量化小菜鸡    时间: 2021-5-20 18:23
本帖最后由 量化小菜鸡 于 2021-5-20 18:25 编辑

我想问个比较白痴的问题，假如用DFT还是用hftype uhf就可以吗？不需要对应写UKS?
另外NEB优化的值建议用多少比较安全。
作者
Author: Freeman    时间: 2021-5-22 12:35
量化小菜鸡 发表于 2021-5-20 18:23
我想问个比较白痴的问题，假如用DFT还是用hftype uhf就可以吗？不需要对应写UKS?
另外NEB优化的值建议用多 ...

就写uhf，不用写uks。没明白“NEB优化的值”指的是啥。
作者
Author: 量化小菜鸡    时间: 2021-5-22 12:40
Freeman 发表于 2021-5-22 12:35
就写uhf，不用写uks。没明白“NEB优化的值”指的是啥。

谢谢，就是这一条
ORCA的neb和几何优化的收敛限都高得过头，所以要手动调低收敛限。见orca_manual_4_2_1.pdf第635页
想问问收敛限可以降到多少。
作者
Author: Freeman    时间: 2021-5-24 20:14
量化小菜鸡 发表于 2021-5-22 12:40
谢谢，就是这一条
ORCA的neb和几何优化的收敛限都高得过头，所以要手动调低收敛限。见orca_manual_4_2_1 ...

不好意思，这几天不经常上论坛，就没回复。
我的习惯是比默认值低一个数量级。
作者
Author: 量化小菜鸡    时间: 2021-5-24 20:21
Freeman 发表于 2021-5-24 20:14
不好意思，这几天不经常上论坛，就没回复。
我的习惯是比默认值低一个数量级。

好的好的，谢谢。
作者
Author: 量化小菜鸡    时间: 2024-9-22 23:38
你好老师，我最近用ORCA6跑老是出现这个报错，不知道您有见过吗？

Starting iterations:

Optim.  Iteration  HEI  E(HEI)-E(0)  max(|Fp|)   RMS(Fp)    dS
Switch-on CI threshold               0.020000
[file orca_numfreq/numfreq_utils.cpp, line 1523]: Error (ORCA_NUMCALC): Child  terminated with kill-11
作者
Author: CCLI    时间: 2025-5-8 00:16
老师好，
reactant.xyz是自旋极化单重态，仅仅给出坐标吗？
作者
Author: Freeman    时间: 2025-5-9 17:46
CCLI 发表于 2025-5-8 00:16
老师好，
reactant.xyz是自旋极化单重态，仅仅给出坐标吗？

第一步并不涉及SCF计算，只是想获得一个初猜的NEB路径，所以只要给个坐标即可。生成初猜路径之后，第一步的任务就可以停了。
作者
Author: CCLI    时间: 2025-5-10 12:36
老师我这个体系产物为对称破缺态，产物为单重态，out文件显示错误，是不是要按照你的步骤来
作者
Author: Freeman    时间: 2025-5-10 14:31
CCLI 发表于 2025-5-10 12:36
老师我这个体系产物为对称破缺态，产物为单重态，out文件显示错误，是不是要按照你的步骤来

信息给太少了，看不出来为什么出错。
作者
Author: lushiyue    时间: 2025-5-13 11:33
只生成初始路径可以直接写 ! neb-idpp ，这样不用杀任务，想要一大批初始路径的时候就很方便
作者
Author: CCLI    时间: 2025-5-15 19:10
量化小菜鸡 发表于 2024-9-22 23:38
你好老师，我最近用ORCA6跑老是出现这个报错，不知道您有见过吗？

Starting iterations:

你好，我遇到了相同的报错，请问有什么好的方法解决？
作者
Author: CCLI    时间: 2025-5-30 10:00
你好老师，ORCA6是不是不支持这个方法，还是我这个体系没有过渡态
作者
Author: Freeman    时间: 2025-6-11 16:28
CCLI 发表于 2025-5-30 10:00
你好老师，ORCA6是不是不支持这个方法，还是我这个体系没有过渡态

这么说，可能这个方法确实不适用于ORCA6吧。我还没有用过ORCA6，所以不知道是什么情况。
作者
Author: CCLI    时间: 2025-6-12 14:34
Freeman 发表于 2025-6-11 16:28
这么说，可能这个方法确实不适用于ORCA6吧。我还没有用过ORCA6，所以不知道是什么情况。

好的，谢谢老师



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