计算化学公社

标题: 在表面催化反应中用以寻找过渡态的方法有哪些？ [打印本页]

作者
Author: minagami 时间: 2015-12-12 22:56
标题: 在表面催化反应中用以寻找过渡态的方法有哪些？
本帖最后由 minagami 于 2015-12-13 18:23 编辑

前一阵子我一直用CASTEP重现一篇论文的结果，论文中用的是VASP中的NEB和CL-NEB方法寻找过渡态，而我用的是CASTEP中的LST/QST方法。经过多次尝试，我发现能否成功找到过渡态很大程度上取决于反应物和产物的优化后构型是否足够接近。我尝试了一下，在金属氧化物表面上，对于一个能用LST/QST方法寻找到过渡态的反应物和产物对(收敛标准为RMS=0.05)，如果把反应物表面的一个甲基旋转到别的角度再优化结构，重新与产物再次以同样的参数进行LST/QST寻找过渡态，结果就是RMS就一直在0.5左右震荡而不能收敛了。

但是反应物和产物的几何优化结构也很大程度取决于初始构型，如果人为地调节反应物和产物的初始构型使其足够接近，那么结果会不会与实际相差较远？因为反应物和产物的稳定构型不应该是（反应物和产物在保持结构的情况下）全势能面的最低点吗？何况对于多步反应，中间产物的构型就更难以进行调整了。因为一旦调整了之前的过渡态还得重新找一遍。而且对以RMS=0.05收敛标准寻找到的过渡态，用Dmol3进行频率计算后会发现还存在有不少虚频。

那么有没有更好的表面催化反应寻找过渡态方法呢？我曾拜读过Sob神的文章《过渡态、反应路径的计算方法及相关问题》http://sobereva.com/44 ，但是我发现实际在常用的能研究周期性体系的量子化学计算程序中，寻找过渡态的方法实际却是十分有限的。比如CASTEP的LST/QST，NEB（用于找IRC验证过渡态而不能用于寻找）；Dmol3的LST/QST，EF（好像要计算整个体系的Hessian矩阵？曾经使用过但是不知为何计算失败），NEB（同上）；和VASP的NEB，CL-NEB（刚装上VASP，打算用但没用过，不好评价）。

大家认为在这几种过渡态寻找方法中那种比较好呢？还有没有其它更有效的寻找表面催化反应寻找过渡态方法呢？欢迎大家进来踊跃讨论~

还有一个问题，能不能用CASTEP和Dmol3中的BFGS几何优化方法（不勾选use line search），根据过渡态初猜寻找更精确的过渡态初猜呢？或者对分别与反应物或产物较为接近的过渡态初猜结构进行优化，得到反应物和产物在晶体表面上的构型，再进行过渡态寻找（有时候对于找过渡态挺有效的，但我也很怀疑其合理性）。也欢迎大家前来讨论这种有点玄学的小技巧。

作者
Author: zyj19831206 时间: 2015-12-12 23:37
能否把你做的那个例子详细的发上来，我最近也在研究过渡态，也遇到很多问题。

作者
Author: minagami 时间: 2015-12-12 23:41

zyj19831206 发表于 2015-12-12 23:37
能否把你做的那个例子详细的发上来，我最近也在研究过渡态，也遇到很多问题。

Chin, Y.-H. (Cathy), Buda, C., Neurock, M., & Iglesia, E. (2013). Consequences of Metal–Oxide Interconversion for C–H Bond Activation during CH4 Reactions on Pd Catalysts. Journal of the American Chemical Society, 135(41), 15425–15442.
以上是原论文，你可以参照一下。我自己的计算过程就不太好发出来了。

作者
Author: minagami 时间: 2015-12-13 00:54

zyj19831206 发表于 2015-12-12 23:37
能否把你做的那个例子详细的发上来，我最近也在研究过渡态，也遇到很多问题。

话说兄弟你那创腾的Dmol3高级教程有电子版或者扫描版的话能分享一下不？最近我也打算研究研究Dmol3。。。

作者
Author: 卡开发发 时间: 2015-12-13 09:32
计算过程是否测试了KPOINTS和ENCUT，使得Total Energy和Force都能够收敛？如果Force的误差比设置的收敛大，显然不能够计算得到正确的结果。

DMol3搜索过渡态之后是会有超过一个的虚频，但原则上可以通过频率计算后进行EF来消除不合理的虚频。

只能说VASP的NEB、CI-NEB更被期刊认可一些，但原则上来说这些方法都能够得到可靠的结果。至于EF失败的原因还得看输出文件才能判断。

CASTEP和DMol3中的BFGS只能用于能量极小化的处理，程序默认的方式是先通过LST找到路径上的极大值，然后CG进行优化。另外，有些体系LST/QST可能解决不了问题，比如CH3Cl+Br-=CH3Br+Cl-的过程，就得利用一些技巧直接通过EF的方法做。原则上《过渡态、反应路径的计算方法及相关问题》http://sobereva.com/44 也应当适用于周期体系。

对于表面反应的处理，如果不涉及到DFT+U或是杂化泛函，尤其是真空层充分的情况，接近的精度下，DMol3比CASTEP或VASP来的快得多。

作者
Author: minagami 时间: 2015-12-13 13:06

卡开发发发表于 2015-12-13 09:32
计算过程是否测试了KPOINTS和ENCUT，使得Total Energy和Force都能够收敛？如果Force的误差比设置的收敛大， ...

计算用的基本上是论文中给出的参数值，另外KPOINTS，ENCUT，smear也用原胞对Total Energy进行过收敛测试了。分别用的是0.04(1/A)，396eV，和0.05，Force倒是没有测试过。

话说计算频率时只需计算表面分子的频率吗？我仔细看了输出文件计算失败原因迭代似乎是超过50次不收敛（之前没注意到，真是惭愧），不过在前50次迭代中三项参数也没有要收敛的趋势，而且一直显示类似于下面警告：
Warning: Magnitude of eigenvalue  22 too small.  Replaced by -0.000100
Warning: Magnitude of eigenvalue  23 too small.  Replaced by -0.000100
Warning: Magnitude of eigenvalue 1 too large.  Replaced by -25.000000
Warning: Magnitude of eigenvalue  99 too large.  Replaced by 25.000000
不知道这是什么原因造成的呢？看来即使增加迭代次数也难以收敛，那么有什么能帮助EF收敛的方法呢？

我的理解是，在NEB、CI-NEB中，反应路径的寻找是考虑到各个images的受力并不断进行优化的，相比而言QST中假设的二次反应路径假设就有点不科学。因为反应物和产物构象离得越远，势能面越复杂，真实反应路径和LST/QST中假设的反应路径就会相差得越大。这也是我认为LST/QST难以找到过渡态和不可靠的原因。

http://sobereva.com/44中介绍了拟牛顿法寻找过渡态的方法，那么同是拟牛顿法的BFGS为什么不能用来寻找过渡态呢？毕竟几何优化的收敛标准只是最大受力小于某个值，而与最后体系有没有虚频没有关系。我是觉得只要初猜和真实过渡态足够接近的话，几何优化的结果还是有可能能得到过渡态结构的。

虽说原则上《过渡态、反应路径的计算方法及相关问题》http://sobereva.com/44 也应当适用于周期体系。但是具体实践起来又应该用什么方法，或者什么软件包呢？求指教。

既然在接近的精度下DMol3比CASTEP或VASP计算得快的话，那么CASTEP或VASP与DMol3相比应该有什么不可取代的优点吧？可以把DMol3算出来的结果作为初猜结构再进行CASTEP中的计算从而减少计算时间吗？我最近刚开始学习DMol3还不太了解，不过用它来对PdO团簇进行计算后发现，虽然每步SCF过程快但是与CASTEP比起来SCF难以收敛。后来把参数改为限制自旋并且把smear提高到0.005Ha，才终于让SCF步收敛了，但是这样不就导致了后续的能量和力的计算结果不准确了吗？

作者
Author: 卡开发发 时间: 2015-12-13 13:28

minagami 发表于 2015-12-13 13:06
计算用的基本上是论文中给出的参数值，另外KPOINTS，ENCUT，smear也用原胞对Total Energy进行过收敛测试 ...

测试过的话一般在出问题就应该从反应物和产物的结构上找原因，频率当然是计算整个体系的频率。这个Warning一般不必理会。

LST/QST是有可能存在一些问题，要不然也不必额外再做过渡态的优化和确定了，你可以认为这个方法就是EF的初猜。所以CASTEP在这点上很没有优势。至于更多的方法为啥没用于其他程序中的原因我不怎么清楚。至于BFGS，Sob的文中已经写了“QN法泛指基于此原理的一类方法，常用的是BFGS(Broyden Fletcher Goldfarb Shanno)，此法对Hessian的修正保持其对称性和正定性，最适合几何优化，但显然不能用于找过渡态。”

CASTEP和VASP与DMol3相比采用平面波，对于高KPOINTS的3D固体优势比DMol3大得多，不仅相对准确而且速度反而快。但是对于2D的表面以及一些大体系，DMol3因为采用局域基组速度快得多，计算量几乎不依赖晶格的大小。除此之外，CASTEP和VASP可以使用杂化泛函和DFT+U，这对于窄带半导体能带的计算有很大的改善，而DMol3无法做到。当然，CASTEP使用杂化泛函只能使用模守恒，计算就更加缓慢，DFT+U还不能进行交换参数J的修正。

至于DMol3难以收敛，就应该再找原因，不妨一块研究一下。

作者
Author: minagami 时间: 2015-12-13 13:57

卡开发发发表于 2015-12-13 13:28
测试过的话一般在出问题就应该从反应物和产物的结构上找原因，频率当然是计算整个体系的频率。这个Warnin ...

反应物和产物的结构都是经过几何优化的。从反应物和产物的结构上找原因，是指要让反应物和产物的结构尽量接近过渡态吗？因为计算整个体系（大概有50来个原子）的频率十分耗时间，我还想着能不能只计算表面反应分子的频率呢。

BFGS不能用于找过渡态的“显然”只是对于Sob而言，我就觉得一点也不显然orz 专门查了一下BFGS的优化原理，我觉得会不会就是因为BFGS中的line search一步导致其无法用于寻找过渡态，如果不勾选use line search的话不知道结果如何。

高KPOINTS的3D固体是针对晶体较小的情形吗？我才知道原来对于2D的周期性体系DMol3的速度要快，原本以为DMol3只是在计算团簇上有优势的，受教了。

DMol3难以收敛我觉得会不会是我团簇的模型没建好，在团簇里面表面原子都是不饱和的，我觉得这会不会对SCF计算产生影响？我建立的团簇模型如下，主要是为了得到用于在CASTEP里进行计算的初猜结构。
(, 下载次数 Times of downloads: 46)

最近刚开始学习DMol3，可能会有不少理解得很肤浅的地方，还请前辈多多指教呀！

作者
Author: 卡开发发 时间: 2015-12-13 14:57

minagami 发表于 2015-12-13 13:57
反应物和产物的结构都是经过几何优化的。从反应物和产物的结构上找原因，是指要让反应物和产物的结构尽量 ...

也有可能反应物或是产物的结构不可靠，这种还得检查一下两者的频率等。BFGS不能用于优化过渡态目测应该不是line search的原因，但这个line search指的是哪种我不了解。

这种不饱和性可能确实会影响SCF收敛，文献上的簇确实取得大得多。CASTEP中按照表面计算？那就不如直接使用DMol3做表面体系的计算了。

谈不上指教，好多不懂得，望多讨论。

作者
Author: minagami 时间: 2015-12-13 15:19

卡开发发发表于 2015-12-13 14:57
也有可能反应物或是产物的结构不可靠，这种还得检查一下两者的频率等。BFGS不能用于优化过渡态目测应该不 ...

关于line search，我的理解是，BFGS中得到近似的Hessian矩阵后可确定步进方向和步长，line search就是在这个方向上对步长的一维优化的过程。

我待会试试再建一个小的表面体系用dmol3计算，但是因为有真空层和周期性的存在，不知道计算时间会不会大幅增加？毕竟这只是在测试一个初猜结构，可以的话我还是挺希望能尽量减少计算量的。一开始我还打算用分子力学方法来进行几何优化的，可惜MS里没找到有Pd2+的力场参数的力场。

我也感觉自己不懂的地方实在是太多了，希望以后也能多多讨论，共同进步

作者
Author: 卡开发发 时间: 2015-12-13 15:23
本帖最后由卡开发发于 2015-12-13 15:26 编辑

minagami 发表于 2015-12-13 15:19
关于line search，我的理解是，BFGS中得到近似的Hessian矩阵后可确定步进方向和步长，line search就是在 ...

line search的事情不太了解。DMol3的真空层几乎占不了多少计算量，仅仅只是Ewald求和的时候涉及一些，而平面波方法做FFT计算量会依赖晶格大小。CASTEP没办法进一步做EF，所以意思不大。
对于化学反应体系来说，分子力场方法几乎行不通，就算反应力场误差应该会很大。

还望多交流。

作者
Author: minagami 时间: 2015-12-13 15:39

卡开发发发表于 2015-12-13 15:23
line search的事情不太了解。DMol3的真空层几乎占不了多少计算量，仅仅只是Ewald求和的时候涉及一些，而 ...

寻找过渡态的话我也觉得用分子力场行不通，但如果是作为对产物和反应物结构的预优化呢？而且在研究吸附位建立反应物结构的时候，我觉得基于分子力场的蒙特卡罗方法还是挺有用的。不过前提在于得有一个好的分子力场。像在MS里我就没找到含有Pd2+参数的力场，只有Pd参数的力场。虽然用蒙特卡罗找了一下吸附位，但总觉得不可靠。

作者
Author: 卡开发发 时间: 2015-12-13 16:03

minagami 发表于 2015-12-13 15:39
寻找过渡态的话我也觉得用分子力场行不通，但如果是作为对产物和反应物结构的预优化呢？而且在研究吸附位 ...

如果吸附是物理吸附还凑合，化学吸附一般来说行不通，力场很难描述这种情况

作者
Author: minagami 时间: 2015-12-13 16:16

卡开发发发表于 2015-12-13 16:03
如果吸附是物理吸附还凑合，化学吸附一般来说行不通，力场很难描述这种情况

化学吸附实际上就是被吸附分子上的原子与吸附剂表面原子产生了化学键吧。如果分子力场有这两种原子之间的作用力参数，能描述这两种原子的成键规律的话，为什么就不能用来描述吸附过程呢？

作者
Author: 卡开发发 时间: 2015-12-13 16:18

minagami 发表于 2015-12-13 16:16
化学吸附实际上就是被吸附分子上的原子与吸附剂表面原子产生了化学键吧。如果分子力场有这两种原子之间的 ...

这样的力场不多吧，至少MS下面的力场应该都不是这么处理。

作者
Author: minagami 时间: 2015-12-13 16:29

卡开发发发表于 2015-12-13 16:18
这样的力场不多吧，至少MS下面的力场应该都不是这么处理。

我对MS的力场的了解仅限于按照MS中的教程“Using Materials Studio to edit a forcefield”改过一下力场，所以可能理解得有所偏差。在MS中力场中的参数里除了L-J势以外，不是还有一系列参数是专门用来描述原子之间的相互作用关系的吗？其中包括Bond stretch、Angle bend、Torsion和Inversion等，这些就是用来描述原子之间成键关系的吧。

作者
Author: 卡开发发 时间: 2015-12-13 17:05

minagami 发表于 2015-12-13 16:29
我对MS的力场的了解仅限于按照MS中的教程“Using Materials Studio to edit a forcefield”改过一下力场 ...

恩，这些必须手工把两个原子之间连上键才会这样描述。

作者
Author: minagami 时间: 2015-12-13 17:14

卡开发发发表于 2015-12-13 17:05
恩，这些必须手工把两个原子之间连上键才会这样描述。

话说在用分子力场进行结构优化的过程中不能实现自动判断两个原子之间是否成键吗?比如用Build-Bonds-Monitor bonding这个选项。

作者
Author: 卡开发发 时间: 2015-12-13 17:28

minagami 发表于 2015-12-13 17:14
话说在用分子力场进行结构优化的过程中不能实现自动判断两个原子之间是否成键吗?比如用Build-Bonds-Monit ...

应该不行吧，何况力场一开始分配之后每一步计算也不会重新分配力场。

作者
Author: minagami 时间: 2015-12-13 17:35

卡开发发发表于 2015-12-13 17:28
应该不行吧，何况力场一开始分配之后每一步计算也不会重新分配力场。

如果是这样的话MS中自带的分子力场方法的确是有严重的局限性呀。我觉得也许可以通过perl脚本实现每产生一个新构像后重新计算成键与否和分配力场，不过应该会大大减慢速度吧

作者
Author: 卡开发发 时间: 2015-12-13 17:40

minagami 发表于 2015-12-13 17:35
如果是这样的话MS中自带的分子力场方法的确是有严重的局限性呀。我觉得也许可以通过perl脚本实现每产生一 ...

速度慢，而且这种成键和非键作用各自的势函数接洽的地方也得处理，问题应该挺多的。

作者
Author: minagami 时间: 2015-12-13 17:51

卡开发发发表于 2015-12-13 17:40
速度慢，而且这种成键和非键作用各自的势函数接洽的地方也得处理，问题应该挺多的。

的确是呀。如果要处理这些的话就大大超出了我的能力范围了orz
不好意思占用了你这么多时间，感谢你的帮忙~

作者
Author: 卡开发发 时间: 2015-12-13 20:12

minagami 发表于 2015-12-13 17:51
的确是呀。如果要处理这些的话就大大超出了我的能力范围了orz
不好意思占用了你这么多时间，感谢你的帮 ...

没有，没解决啥问题，能够讨论就好。

作者
Author: sobereva 时间: 2015-12-13 22:53
BFGS保持Hessian的正定性，也就是说肯定会保持各个振动模式对应的频率都是正值，显然这种情况只能是搜索极小点。因为过渡态会有一个虚频。

模拟过程中自动判断成键重新分配力场并没有意义，因为程序里用的力场不是反应力场，而是要求模拟过程中拓扑关系是固定的，拓扑关系在模拟过程中不会发生变化，所以原子类型一直不会变，参数的分配也不会变。这些力场都没法描述键的形成与断裂，所以化学吸附肯定也没法描述。表面化学反应就不要试图考虑力场了。

作者
Author: minagami 时间: 2015-12-14 03:40

sobereva 发表于 2015-12-13 22:53
BFGS保持Hessian的正定性，也就是说肯定会保持各个振动模式对应的频率都是正值，显然这种情况只能是搜索极 ...

多谢你的回复和解答！我开始能理解BFGS为什么不能用于搜索过渡态了。

对表面化学反应本身想通过分子力场来解决看来也是不太可能了，不过毕竟分子在进行化学吸附之前得先进行物理吸附，化学吸附位点也应该就在所找到的物理吸附位点的附近吧？那么能否通过蒙特卡罗方法对众多可能的化学吸附位点进行一个预筛选呢？

作者
Author: sobereva 时间: 2015-12-14 03:42

minagami 发表于 2015-12-14 03:40
多谢你的回复和解答！我开始能理解BFGS为什么不能用于搜索过渡态了。

对表面化学反应本身想通过分子力 ...

是的，用力场研究物理吸附，以提供化学吸附量化研究的初始参考结构还是可以的。

作者
Author: zyj19831206 时间: 2015-12-16 22:03

minagami 发表于 2015-12-14 03:40
多谢你的回复和解答！我开始能理解BFGS为什么不能用于搜索过渡态了。

对表面化学反应本身想通过分子力 ...

兄弟，你举个例子出来看看啊，光说理论我们不懂啊

作者
Author: minagami 时间: 2015-12-16 23:42

zyj19831206 发表于 2015-12-16 22:03
兄弟，你举个例子出来看看啊，光说理论我们不懂啊

你指的例子是指哪方面的呢？如果是指蒙特卡罗方法找吸附点的话，官方帮助文件里有一个相应的教程“Determining the location of SO2 on the Ni(111) surface with Adsorption Locator”。

作者
Author: zyj19831206 时间: 2015-12-17 00:13

minagami 发表于 2015-12-16 23:42
你指的例子是指哪方面的呢？如果是指蒙特卡罗方法找吸附点的话，官方帮助文件里有一个相应的教程“Determ ...

你重复文献得到的结果，文献李省略了很多。

作者
Author: minagami 时间: 2015-12-17 00:41

zyj19831206 发表于 2015-12-17 00:13
你重复文献得到的结果，文献李省略了很多。

我主要是重复了其中的过渡态寻找工作，有什么操作上的问题的话我们可以一起讨论，不过具体过程的话就不方便发出来了。

作者
Author: zyj19831206 时间: 2015-12-17 00:47

minagami 发表于 2015-12-17 00:41
我主要是重复了其中的过渡态寻找工作，有什么操作上的问题的话我们可以一起讨论，不过具体过程的话就不方 ...

要发文章吗？那就不追问了。

欢迎光临计算化学公社 (http://ccc.keinsci.com/)