小窍门：使用自旋极化单重态波函数作为NEB的初猜

Freeman · 发表于 Post on 2021-4-21 12:27:27

本帖最后由 Freeman 于 2021-4-22 11:23 编辑

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

假如我们现在有两个优化好的反应物和生成物坐标，reactant.xyz和product.xyz，其中reactant.xyz是自旋极化单重态。如果我们直接写neb1.inp

! neb ......
......
%neb
neb_end_xyzfile "product.xyz" #给出生成物坐标
end
* 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

! ......
......
%scf
HFTyp uhf #使用非限制性自洽场
STABPerform true #做波函数稳定性检测
STABRestartUHFifUnstable true #如果检测出波函数不稳定，那么再做一遍自洽场找稳定波函数（有极低可能最后仍然没找到，保险起见可以再创建一个inp文件以进一步检测波函数稳定性）
STABDTol 0.00001 #设置波函数稳定性的判定阈值，见下面“其他细节3”
STABRTol 0.00001
end
* xyz 0 1
......
*

复制代码

运行后获得一系列稳定的波函数文件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

......
>
2
image2
H 0.0 0.0 0.0
H 0.0 0.7 0.0
>
2
image3
H 0.0 0.0 0.0
H 0.0 0.8 0.0
>
......

复制代码

四、创建并运行最终的neb任务
最终的neb任务inp文件写为neb3.inp

! neb ......
%neb
neb_restart_xyzfile "neb1_MEP.allxyz" #给出allxyz文件作为初始中间点，这些中间点包括了反应物和生成物，但是下面第10行还要写一遍反应物坐标
neb_restart_gbwname "neb2" #给出波函数初猜文件，“neb2”表示波函数初猜文件为“neb2_imN.gbw”
end
%scf
hftyp uhf #仍然要指定使用非限制性自洽场
end
*xyz 0 1
[neb2_im0.inp对应的坐标]
*

复制代码

运行后得到结果。

其它细节：
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楼的回复。

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

喵星大佬 · 发表于 Post on 2021-4-21 13:29:18

orca不支持的是meta泛函的解析Hessian，像B3LYP，PBE0之类的普通杂化泛函是可以的

biogon · 发表于 Post on 2021-5-19 21:49:27

不过说实话能用高斯的opt=ts我肯定不会用这玩意

量化小菜鸡 · 发表于 Post on 2021-5-20 18:23:28

本帖最后由量化小菜鸡于 2021-5-20 18:25 编辑

我想问个比较白痴的问题，假如用DFT还是用hftype uhf就可以吗？不需要对应写UKS?
另外NEB优化的值建议用多少比较安全。

Freeman · 发表于 Post on 2021-5-22 12:35:20

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

就写uhf，不用写uks。没明白“NEB优化的值”指的是啥。

量化小菜鸡 · 发表于 Post on 2021-5-22 12:40:30

Freeman 发表于 2021-5-22 12:35
就写uhf，不用写uks。没明白“NEB优化的值”指的是啥。

谢谢，就是这一条
ORCA的neb和几何优化的收敛限都高得过头，所以要手动调低收敛限。见orca_manual_4_2_1.pdf第635页
想问问收敛限可以降到多少。

Freeman · 发表于 Post on 2021-5-24 20:14:08

量化小菜鸡发表于 2021-5-22 12:40
谢谢，就是这一条
ORCA的neb和几何优化的收敛限都高得过头，所以要手动调低收敛限。见orca_manual_4_2_1 ...

不好意思，这几天不经常上论坛，就没回复。
我的习惯是比默认值低一个数量级。

量化小菜鸡 · 发表于 Post on 2021-5-24 20:21:50

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

好的好的，谢谢。

量化小菜鸡 · 发表于 Post on 2024-9-22 23:38:28

你好老师，我最近用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

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[ORCA] 小窍门：使用自旋极化单重态波函数作为NEB的初猜

评分 Rate

评分 Rate