计算化学公社

标题: 请教：CP2K做过渡态、电子结构等静态计算时的可靠性 [打印本页]

作者
Author: djjj148 时间: 2021-2-18 17:59
标题: 请教：CP2K做过渡态、电子结构等静态计算时的可靠性
本帖最后由 djjj148 于 2021-2-18 17:59 编辑

听说CP2K运行速度快，但是擅长做AIMD的它做静态计算的精度没有VASP和QE准。于是测试了自己的研究体系：同一个初猜结构下，CP2K在做结构优化速度十分的快，而且收敛性非常好(监控每一个离子步的能量变化发现几乎都是逐渐降低)，反观VASP在结构优化时，最大受力总是会不断振荡才缓慢收敛。70多原子体系(晶格参数：16.7308 12.7134 25)的结构优化任务，CP2K(OT算法111K点)可以在1-3小时做完，而VASP(K点：1 2 1)需要1-2天，就算改为111 K点搭配vasp_gam版本也大概需要一天多。结构优化尚且如此，如果做过渡态计算或更大体系的话用时应该会更加悬殊。
然而，对于优化收敛的结构，分别做bader电荷分析发现，后处理VASP的输出文件生成的ACF.dat文件中比较正常，而CP2K的ACF.dat竟然出现了这样的情况：大部分C原子的CHARGE都是接近赝势中的价电子，没有差很多，然而小部分的C原子的CHARGE却是0，这似乎已经到了离谱的状态。相应的CP2K的inp文件已上传，请教一下大家：
1. 是否是因为inp中的设置的精度不够才导致如此“离谱”的结果？
2. 我在ACS.cata这样的期刊上看到的纯模拟的文章都是用VASP做的静态计算(如过渡态)，CP2K做AIMD或者全部用VASP。然而，由于自己能用的计算硬件不够多和硬，还是眼馋CP2K的速度，不知道在恰当设置inp后，质量较好的期刊能否接受"CP2K做结构优化和过渡态计算得到的能量"？
3. 如果用CP2K优化极小点或过渡态结构后，用VASP算单点能量和进行下一步电子结构分析，应该可以提高结果的可靠性吧？

作者
Author: sobereva 时间: 2021-2-18 19:47
1 不是精度的事。没有具体结构，也不好去判断。用CP2K没必要用bader电荷（AIM电荷）说事，本身CP2K自己就能算诸如Hirshfled-I、REPEAT这样通常比很垃圾的AIM电荷有意义得多的原子电荷，还可以结合Multiwfn算CM5电荷。而且CP2K用medium的print级别直接就输出Mulliken电荷，也能定性说明问题。

2 完全可以。TZV2PX-MOLOPT-GTH结合足够高的截断能已经能达到很高的精度。

3 没这个必要。CP2K结合恰当的设置产生的波函数做电子结构分析已经很好了。

顺带推荐个我之前审的一篇算沸石的文章，Materials Today Communications 26 (2021) 102028（https://doi.org/10.1016/j.mtcomm.2021.102028），全都是CP2K算的，还用Multiwfn做了IGM分析直观考察了沸石对芳香分子的吸附。

作者
Author: djjj148 时间: 2021-2-18 21:45

sobereva 发表于 2021-2-18 19:47
1 不是精度的事。没有具体结构，也不好去判断。用CP2K没必要用bader电荷（AIM电荷）说事，本身CP2K自己就能 ...

谢谢sob老师指点。还想请教一下，因为TZV2PX-MOLOPT-GTH是常用的BASIS_MOLOPT中精度最高的基组，我刚才测试了一下，其他参数不变，用TZV2PX-MOLOPT-GTH的耗时：
  Step    Update method    Time Convergence       Total energy Change
  ------------------------------------------------------------------------------
   1 OT CG    0.12E+00 9.1    0.01652627    -457.3860222786 -4.57E+02
   2 OT LS    0.17E+00 8.5                   -470.3961998629
   3 OT LS    0.22E+00 8.5                   -472.8613114815
   4 OT CG    0.18E+00 15.7    0.01742506    -473.7860931582 -1.64E+01

用DZVP-MOLOPT-SR-GTH：
   1 OT CG    0.12E+00 2.5    0.02632822    -459.3078393090 -4.59E+02
   2 OT LS    0.17E+00 1.7                   -471.9480597834
   3 OT LS    0.16E+00 1.6                   -472.8442457990
   4 OT CG    0.13E+00 3.1    0.03062221    -472.8461111027 -1.35E+01
耗时大概变成了原来3.6-5倍，幅度还是不小的。我想请教一下，像过渡态这样的优化任务用DZVP或TZVP-MOLOPT-GTH得到结构，然后用TZV2PX-MOLOPT-GTH和更高的截断能算单点和分析波函数可以吗？

作者
Author: 卡开发发 时间: 2021-2-18 21:50
本帖最后由卡开发发于 2021-2-18 22:00 编辑

1、CP2K的GPW是结合赝势的，所以赝波函数只有价层的且是smooth的，有可能不能正确识别NCP，也许得考虑缀加芯密度甚至是额外价层的赝波函数变换到全电子价层密度才行。实在不行也可以试试GAPW的方式。
3、并没有办法保证CP2K算得的结构能量差（反应能和能垒）能够和VASP一致，如果出现这种不一致可能会比较难解释（比如DOI: 10.1021/acs.jpcc.9b03554）。

如果只讨论精度，CP2K其实做了deltatest（https://www.cp2k.org/_media/events:2019_ghent:automation.pdf；其他程序的deltatest测试可以看https://molmod.ugent.be/deltacodesdft），可靠性方面我不做评价，可自行比对。

作者
Author: sobereva 时间: 2021-2-19 07:22

djjj148 发表于 2021-2-18 21:45
谢谢sob老师指点。还想请教一下，因为TZV2PX-MOLOPT-GTH是常用的BASIS_MOLOPT中精度最高的基组，我刚才测 ...

可以

另外，对于大多数波函数分析的目的，比如RDG、IRI、ELF、AIM分析等，其实TZVP-MOLOPT-GTH就够了，明显没有能量对基组敏感性那么高。

PS：对于在Multiwfn中做波函数分析，我建议额外用pob-TZVP基组做个GAPW单点计算得到molden文件。一方面pob-TZVP是基于片段收缩的，在Multiwfn里分析耗时低得多得多（比DZVP-MOLOPT-SR-GTH等价的GTF数能少好几倍），另一方面这样产生的波函数直接描述了内核电子，便于涉及到内核电子结构的分析。

作者
Author: djjj148 时间: 2021-2-19 09:47

sobereva 发表于 2021-2-19 07:22
可以

另外，对于大多数波函数分析的目的，比如RDG、IRI、ELF、AIM分析等，其实TZVP-MOLOPT-GTH就够了 ...

谢谢sob老师！请问现在Multiwfn是可以分析CP2K产生的基于周期性结构的波函数了吗？我看了一下好像还没有相应的博文。还是目前按照分析Gaussian和ORCA的fchk、molden文件一样直接载入Multiwfn中分析呢？如果是这样，晶格边界的电子性质是不是不可用呢？

作者
Author: djjj148 时间: 2021-2-19 09:48

卡开发发发表于 2021-2-18 21:50
1、CP2K的GPW是结合赝势的，所以赝波函数只有价层的且是smooth的，有可能不能正确识别NCP，也许得考虑缀加 ...

谢谢卡卡老师~！

作者
Author: 卡开发发 时间: 2021-2-19 10:09

djjj148 发表于 2021-2-19 09:47
谢谢sob老师！请问现在Multiwfn是可以分析CP2K产生的基于周期性结构的波函数了吗？我看了一下好像还没有 ...

最新的Multiwfn手册2.9.2节有提到。

作者
Author: gog 时间: 2021-2-19 12:20
可以用 QZVP开头的基组和适合的截断能获得cp2k能达到的最高精度

作者
Author: k64_cc 时间: 2021-2-19 12:53
CP2K的能量精度很可以的。我以前测过几个有机反应体系，在调好截断能的情况下计算两个构象间能量差。同个泛函、相同水平基组下，CP2K和Gaussian的dE差距甚小。

MOLOPT基组的使用细节我记得前几天有帖子。

作者
Author: 卡开发发 时间: 2021-2-19 13:35

k64_cc 发表于 2021-2-19 12:53
CP2K的能量精度很可以的。我以前测过几个有机反应体系，在调好截断能的情况下计算两个构象间能量差。同个泛 ...

如果只是CHO体系，从delta test的结果反映来看误差也确实不大，但到了含过渡金属或者其他元素的情况就不是很有保障了，比如前面提到的CO@ZnO的情况。

作者
Author: k64_cc 时间: 2021-2-20 00:51
本帖最后由 k64_cc 于 2021-2-20 00:53 编辑

卡开发发发表于 2021-2-19 13:35
如果只是CHO体系，从delta test的结果反映来看误差也确实不大，但到了含过渡金属或者其他元素的情况就不 ...

啊这个，我看他输入文件里也没有过渡金属，最多就一个B，就觉得应该还行吧……

作者
Author: sobereva 时间: 2021-2-20 15:23

djjj148 发表于 2021-2-19 09:47
谢谢sob老师！请问现在Multiwfn是可以分析CP2K产生的基于周期性结构的波函数了吗？我看了一下好像还没有 ...

最近会写一系列帖子。前阵子太忙没空写

作者
Author: sobereva 时间: 2021-2-20 15:34
顺带发一个最近看到的Chem. Sci., 2017, 8, 3554–3565的补充材料的对比，可见CP2K可以和量子化学程序吻合得相当好，即便是用GTH赝势做GPW也是如此

(, 下载次数 Times of downloads: 78)

作者
Author: djjj148 时间: 2021-2-20 23:46

sobereva 发表于 2021-2-20 15:23
最近会写一系列帖子。前阵子太忙没空写

给卢老师打call，Multiwfn支持对CP2K的分析让这把瑞士军刀更为锋利了！

作者
Author: 卡开发发 时间: 2021-2-21 04:21
上面的deltatest我补充解释一下：使用相同的赝势情况下，abinit基于平面波可以通过调整波函数的动能截断来逼近基组极限，也就是说abinit可以认为是GTH赝势下的基组极限结果，随着zeta增加其实也并未趋于基组极限。S22测试中只包含了一部分主族元素，但过渡金属可能会出现问题，比如（DOI: 10.1021/acs.jctc.8b00957的SI）
(, 下载次数 Times of downloads: 63)
CP2K做过渡金属体系的文章其实也并不是很多，具体原因不明，可能还有其他奇怪的例子吧。其实用什么程序无所谓，但建议各个程序使用前先对精度大致做好调研，然后对自己体系最好做好充分的测试。

作者
Author: 邱美佳 时间: 2021-3-6 14:48

sobereva 发表于 2021-2-19 07:22
可以

另外，对于大多数波函数分析的目的，比如RDG、IRI、ELF、AIM分析等，其实TZVP-MOLOPT-GTH就够了 ...

请问sob老师cp2k如果是算DOS那要选什么级别的基组结果比较可靠呢？

作者
Author: sobereva 时间: 2021-3-7 05:34

邱美佳发表于 2021-3-6 14:48
请问sob老师cp2k如果是算DOS那要选什么级别的基组结果比较可靠呢？

通常DZVP-MOLOPT-SR-GTH就可以，想更好点就用TZVP的

作者
Author: 卡开发发 时间: 2021-3-9 09:59

邱美佳发表于 2021-3-6 14:48
请问sob老师cp2k如果是算DOS那要选什么级别的基组结果比较可靠呢？

建议上TZ以上的（或者小心测试），DZ可能会高估gap，参考DOI: 10.1063/1.4736998。

作者
Author: sobereva 时间: 2021-3-16 08:16
CP2K的MOLOPT系列比同等zeta档次的其它基组精度往往高一个档，如根据原文测试，DZVP-MOLOPT-GTH算能量问题上就已经略胜于TZVP-MOLOPT-GTH，因此有较高概率DZVP-MOLOPT-SR-GTH就已经够用。不放心的话可以和pob-TZVP或TZV2P-GTH对比一下。

作者
Author: 丁越 时间: 2021-10-19 14:14

卡开发发发表于 2021-2-18 21:50
1、CP2K的GPW是结合赝势的，所以赝波函数只有价层的且是smooth的，有可能不能正确识别NCP，也许得考虑缀加 ...

请问一下卡开发发老师，我用CP2K优化结构，用VASP计算能量，这样的方式审稿人能接受么？主要是CP2K算的很快，想利用这个优势。
其次，我算了一个含Cu、Fe的体系吸附NO、CO等，但是计算下来的结合能和VASP利用CP2K优化完的结构算结合能相差太离谱了。。肯定是VASP的结果更靠谱一点

作者
Author: Aridea 时间: 2021-10-19 18:45

丁越发表于 2021-10-19 14:14
请问一下卡开发发老师，我用CP2K优化结构，用VASP计算能量，这样的方式审稿人能接受么？主要是CP2K算的很 ...

我最近也算了一套cp2k的弱相互作用，结合能就比较离谱，1-6eV的都有

作者
Author: 卡开发发 时间: 2021-10-23 18:29

丁越发表于 2021-10-19 14:14
请问一下卡开发发老师，我用CP2K优化结构，用VASP计算能量，这样的方式审稿人能接受么？主要是CP2K算的很 ...

建议你可以找个稍微小一些的体系测试一下吸附能之类的数据。目前我也说不好是赝势、基组还是算法的问题，影响的因素很多，例如前面看到赝势大小芯也会显著影响精度。

作者
Author: 卡开发发 时间: 2021-10-23 18:42

Aridea 发表于 2021-10-19 18:45
我最近也算了一套cp2k的弱相互作用，结合能就比较离谱，1-6eV的都有

因素比较多，我暂时不好说。或者可以把数据输入拿出来讨论。

作者
Author: Aridea 时间: 2021-10-23 21:00

卡开发发发表于 2021-10-23 18:42
因素比较多，我暂时不好说。或者可以把数据输入拿出来讨论。

大大，我感觉有点像反应了，必经是带电的体系，电子转移啥的也是有可能的嘛

作者
Author: 卡开发发 时间: 2021-10-23 21:23

Aridea 发表于 2021-10-23 21:00
大大，我感觉有点像反应了，必经是带电的体系，电子转移啥的也是有可能的嘛

带电怎么个带电法？你直接体系加1电荷？周期性边界的计算这个要非常小心。

作者
Author: Aridea 时间: 2021-10-23 21:59

卡开发发发表于 2021-10-23 21:23
带电怎么个带电法？你直接体系加1电荷？周期性边界的计算这个要非常小心。

嗯嗯，大大我是用cp2k计算的2Dxy slab加的MT求解器，应该可以吧

作者
Author: 卡开发发 时间: 2021-10-24 13:21

Aridea 发表于 2021-10-23 21:59
嗯嗯，大大我是用cp2k计算的2Dxy slab加的MT求解器，应该可以吧

可以。

作者
Author: Aridea 时间: 2021-10-24 16:49

卡开发发发表于 2021-10-24 13:21
可以。

嘿嘿，当初我发帖子大大还有社长指点过了，感谢感谢~

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