计算化学公社

标题: 关于TS校正和Grimme校正在算色散作用能量时的疑问 [打印本页]

作者
Author: jiewei 时间: 2014-12-3 11:02
标题: 关于TS校正和Grimme校正在算色散作用能量时的疑问
请问一下各位老师和同学：
最近我用GGA/PBE-TS体系算了一个金属性类石墨烯体系和分子的范德华作用，算出的结合能为2.0eV左右，这个结合能是非常高的。
看文献：Chemical Reactivity of Electron-Doped and Hole-Doped Graphene 中提到G-TTF和G-F4TCNQ在M06和VDW-DF方法下计算的结合能大约为0.7-1.1.56 eV
文献Energy Level Alignment and Charge Carrier Mobility in Noncovalently Functionalized Graphene提到G和TTF的结合能大约为0.73 eV用的是vasp算的方法为GGA/PBE加了vasp里面的色散校正也就是第二代Grimme。

我的问题就是：TS校正（Alexandre Tkatchenko and Matthias Scheffler, Accurate Molecular Van Der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data, Phys. Rev. Lett., 2009, 102, 073005）
和Grimme校正比较是否会高估结合能？
另外SIESTA这个软件的范德瓦尔斯校正是用的Grimme校正吧？

非常感谢各位老师和同学。

作者
Author: 卡开发发 时间: 2014-12-4 12:39
SIETSA的DFTD采用的是Grimme校正即DFTD2（G06），私下用DMol3做了些测试（用了xfine的精度），测试做了如苯的堆叠和稀有气体相互作用以及晶格常数，普遍结果是高估弱作用，并且低估晶格常数（未估计应力），另外Langmuir 2010, 26(18), 14582–14593这篇文章就对s6因子进行了校正。从手册上看，SIESTA原则上也可根据测试结果也进行对MM.Grimme.S6的修正。

虽然使用的程序稍有不同，不知能否给您带来帮助。希望进一步讨论。

作者
Author: jiewei 时间: 2014-12-4 13:02

卡开发发发表于 2014-12-4 12:39
SIETSA的DFTD采用的是Grimme校正即DFTD2（G06），私下用DMol3做了些测试（用了xfine的精度），测试做了如苯 ...

感谢，昨天我搜索了一堆文献，看了看，发现SIETSA好像就是用的D2，今天专业人士来解答，非常感谢。
我用dmol的GGA/PBE-Grimme算一个体系的结合能大约是1.89eV，用castep算出来为0.715 eV。
用GGA/PBE-TS算的结合能在dmol中是-2.32，castep中是-1.037.
dmol中用的fine精度（k点我试过调比较高，影响不大），用的DNP和 all electron 。
现在看这个样子dmol里面确实就是高估弱作用了。
castep算出来的结合能还比较靠谱。用的是莫守恒赝势，680eV cutoff

我现在就是不知道这个产生这么大能量差异的原因是？
不知道楼上可知道么？谢谢。

作者
Author: 卡开发发 时间: 2014-12-4 13:52

jiewei 发表于 2014-12-4 13:02
感谢，昨天我搜索了一堆文献，看了看，发现SIETSA好像就是用的D2，今天专业人士来解答，非常感谢。
我用 ...

不能就这么说DMol3高估结合能，总觉得还是应该找一下原因，因为不知道实验值的情况不能说哪个程序做出来就不对。

DFTD按说是不依赖电子态的，只依赖离子结构，两个程序的DFTD参数应该是一样的。如果结构差异不是很大（也不该很大）的话这部分有可能是非DFTD部分本身导致的。建议可以找个化学键体系做个测试看看。

DMol3有时候Cutoff选取要相当小心，稳妥起见不建议<4A；另外多极展开最好是取到Hex（16极）。这些参数选取也并不是任意的，数值基组使用的时候还是需要小心，有时候做个测试比较稳妥。新手相当容易犯的毛病就是计算能差的时候，做差的体系采用的截断半径不一致导致作用能高估。

印象中自己以前用DMol3做的测试结果和文献中采用VASP、PWscf这些程序有小的出入，但普遍差异并不大。

作者
Author: jiewei 时间: 2014-12-4 21:18
本帖最后由 jiewei 于 2014-12-4 21:21 编辑

卡开发发发表于 2014-12-4 13:52
不能就这么说DMol3高估结合能，总觉得还是应该找一下原因，因为不知道实验值的情况不能说哪个程序做出来 ...

感谢你一直热心的讨论，我今天测试了一下， cutoff的值影响非常小，不大于0.2eV，8极和16极影响也非常小，也不大于0.2eV。
另外我测试了一下TS校正确实相对于grimme校正算的结合能大。我的体系差大约0.34eV。这个数值还是蛮大的，因为上两个都是不大于。再次感谢！

作者
Author: 卡开发发 时间: 2014-12-5 01:20
我做了个2x2的石墨烯体系，真空20A，不做优化，专门看了一下两个模块DFTD项，两个体系的能量差在5meV/atoms也不算是小数目了，DMol3是高一些，DFTD的方法用的是TS。这个问题倒是还真有些不可理解呢，有点伤脑筋。。因为看不到程序怎么写的，也是在搞不清到底是怎么处理的。

作者
Author: jiewei 时间: 2014-12-5 09:02

卡开发发发表于 2014-12-5 01:20
我做了个2x2的石墨烯体系，真空20A，不做优化，专门看了一下两个模块DFTD项，两个体系的能量差在5meV/atoms ...

嗯嗯，感谢你的测试，很有帮助。从我们的测试其实已经可以看出来dmol其实在计算能量的时候可能不是很准呢。

作者
Author: 卡开发发 时间: 2014-12-5 14:25

jiewei 发表于 2014-12-5 09:02
嗯嗯，感谢你的测试，很有帮助。从我们的测试其实已经可以看出来dmol其实在计算能量的时候可能不是很准呢 ...

这个暂时不能说到底是方法本身可能就是高估了，还是说哪个程序数值处理不当导致就是高估或是低估。另外，如果总作用能DMol3算的偏大也不表示仅仅只是DFTD两个程序上的差异，两者不进行DFTD的时候计算结果也有差别。这种影响太复杂，原因可能很难找。

不过经验上看TS肯定高估，至少同一个程序TS方法计算的结果就比G06大得多，建议试试看G06.如果必要的话只好做些测试对s6项做一些校正。

作者
Author: 卡开发发 时间: 2014-12-10 14:30

jiewei 发表于 2014-12-5 09:02
嗯嗯，感谢你的测试，很有帮助。从我们的测试其实已经可以看出来dmol其实在计算能量的时候可能不是很准呢 ...

你算的结果也许有问题，我用DMol3做出来的结果是-0.9056eV，按照Chemical Reactivity of Electron-Doped and Hole-Doped Graphene 这篇文献重复了一下，但计算参数没有取和文献一致的，我自己额外做了测试。

我采用的方法：

Quality=fine
Method=PBE-GD2

SCF
Basis Set=DNP (v3.5)
Cutoff Radii=4.5 Angstroms
Multi Polar Exp=16 (Hex)

KPOINTS=11 11 1
Occ=Fermi

k点经过测试，控制在收敛差<1meV/atoms。晶格常数处理的方式为固定真空层30A，测试晶格常数14.72~14.83 Angstroms及相应能量，将晶格常数-键能分别做二次方、三次方曲线，极小点均接近14.7925 Angstroms. 文献中的结果除了能量，分子与石墨烯的间距为3.5 Angstroms，我做出来为3.428 Angstroms. 我把部分结果上传上来，附件中包含着优化结构的原子坐标文件（PBC结构用cif，非PBC用pdb，为了防止xsd版本差异导致打不开）以及计算结果（outmol），我们再进行讨论。

作者
Author: jiewei 时间: 2014-12-10 17:02

卡开发发发表于 2014-12-10 14:30
你算的结果也许有问题，我用DMol3做出来的结果是-0.9056eV，按照Chemical Reactivity of Electron-Doped ...

好的，非常感谢，我最近在忙其他事情，等我忙完，看一下你的结果。以及发一下我的结果给你。预计两三天后，非常感谢！我的结果测试初步显示应该在真空层多加k点。

作者
Author: 卡开发发 时间: 2014-12-10 18:41

jiewei 发表于 2014-12-10 17:02
好的，非常感谢，我最近在忙其他事情，等我忙完，看一下你的结果。以及发一下我的结果给你。预计两三天 ...

k如果不充分甚至算出来的gap≠0，这个在我的测试的xlsx文件有一个数据就是这样的。

作者
Author: jiewei 时间: 2014-12-17 14:22

卡开发发发表于 2014-12-10 14:30
你算的结果也许有问题，我用DMol3做出来的结果是-0.9056eV，按照Chemical Reactivity of Electron-Doped ...

我测试的时候发现是这样的，661的和666的k点差别还是蛮大的，最后我选的666，你选的11*11*1和666的k点也差不多，这让我就有点费解，C方向上要不要加k点呢？ MS公司的人给我的答案是除了界面问题都要加。按照你的结果来看如果我a和b都充分了，似乎好像又可以不加了。
我的测试如下：

GGA/PBE-TS-DNP4 (Dmol3) -1.233
GGA/PBE-Grimme-DNP4 (Dmol3) -0.897
GGA/PBE-Grimme (Castep) -0.715
GGA/PBE-Grimme (VASP)a -0.727(J. Phys. Chem. Lett. 2013, 4, 2158㈢165)
M06-L/6-31G(D) (SIESTA)b -0.690(J. Phys. Chem. C 2013, 117, 38953902)
M06-L/6-311G(D) (SIESTA)b -0.871(J. Phys. Chem. C 2013, 117, 38953902)
从测试来看，数值基组和数值基组还是很匹配的，赝势和赝势还是很匹配的。 TS的校正比grimme约高0.3eV。

作者
Author: jiewei 时间: 2014-12-17 14:23

卡开发发发表于 2014-12-10 18:41
k如果不充分甚至算出来的gap≠0，这个在我的测试的xlsx文件有一个数据就是这样的。

另外，非常感谢开发卡卡，谢谢，论坛有卡发卡卡这样无私的成员必将发扬广大:)

作者
Author: 卡开发发 时间: 2014-12-17 15:38

jiewei 发表于 2014-12-17 14:22
我测试的时候发现是这样的，661的和666的k点差别还是蛮大的，最后我选的666，你选的11*11*1和666的k点也 ...

客气了，问题重在讨论，大家相互学习，营造好的学术氛围。

如果你做的是石墨烯，这样的话c轴不涉及周期性的问题，就是个界面的问题，k取0而不应再扩大，我想这有可能就是问题所在了。需要小心的是，做k点测试的时候最好还是看看outmol中的gap是不是=0，否则结果就无意义了。

v4.0基组我测试了一些金属格子的BM方程的拟合曲线，偏离的都比较惨，所以一般都不用，这个建议仅做参考吧。

作者
Author: jiewei 时间: 2014-12-19 10:47

卡开发发发表于 2014-12-17 15:38
客气了，问题重在讨论，大家相互学习，营造好的学术氛围。

如果你做的是石墨烯，这样的话c轴不涉及周 ...

哈哈，你和MS代理商的技术人员给的意见恰恰想反，他们建议是c方向上要加上k点。
DNP4的测试非常重要。非常感谢分享经验。

作者
Author: 卡开发发 时间: 2014-12-19 11:02

jiewei 发表于 2014-12-19 10:47
哈哈，你和MS代理商的技术人员给的意见恰恰想反，他们建议是c方向上要加上k点。
DNP4的测试非常重要。非 ...

石墨烯的话应该没啥争议才是啊，因为考虑的就是表面的计算，如果你算的是石墨的话就要考虑加上k了，难道你做的是石墨夹层？

我指的是4.4的基组。实际上DMol3还有3.8和4.0.0的基组，这个基组的性能我没做过测试，据说4.0.0主要是结合vwn-bp做溶剂化的计算的，究竟能有多好还得详细做测试，如果算Raman就不得不采用4.4的基组了。

DMol3的赝势比起SIESTA的赝势好很多，尤其是采用PBE泛函的时候。但是基组就不如SIESTA那样定义起来比较灵活了。

作者
Author: 沙场新兵001 时间: 2015-1-8 14:23

jiewei 发表于 2014-12-4 13:02
感谢，昨天我搜索了一堆文献，看了看，发现SIETSA好像就是用的D2，今天专业人士来解答，非常感谢。
我用 ...

其实我觉得不能只看几个不同软件计算出来的结果来相互比较，比如VASP和DMOL3相比较的话，所设置的参数的内容是不一样的，并且不同软件的算法近似是不一样的，dmol相对来说要粗糙一些，并且你在做A吸附到B上的时候，适合A的计算参数在计算B时可能就不会那么准确，对不同性质的材料来说不同的软件算出来的结果肯定是有差异的。

作者
Author: 卡开发发 时间: 2015-1-8 15:31

沙场新兵001 发表于 2015-1-8 14:23
其实我觉得不能只看几个不同软件计算出来的结果来相互比较，比如VASP和DMOL3相比较的话，所设置的参数的 ...

不同程序结果是会有些差异，但一般写程序只要保证的就是再现性，至少是半定量的，并且有人做了对比https://molmod.ugent.be/deltacodesdft 。虽然由于不同的计算、数值方法（非理论级别）会有些小的差异，比如LCAO和平面波、赝势或PAW等差异，但这些误差原则上不应该出现剧烈的差异。尤其是商业程序，这个问题原则上是要保证的。保证了半定量的差异的前提下，不同计算方法就仁者见仁智者见智了，对于k不太密集的体系LCAO方法显然非常有优势。

我不知道“DMol3是粗糙的”这个定论怎么来的（虽然几年前我曾这么认为过，因为那时候我也只会调调界面上那几个参数，我不知道这能不能算是偏见，我那时候就有些偏见），这个想法太过武断，算得快不一定精度差。A吸附到B上的时候，单独计算A和B两者本来选用的参数就有差异，无论VASP还是WIEN2K都不能避免，绝大多数情况下，A单独处理面临的结果必然是大的格子以及Gamma-k，显然很多时候难保障和A-B的体系一致。同时，大量的隐藏参数在Keywords中都有说明，到底怎么计算更精确不仅仅只是程序本身，还有操作者决定（如收敛标准、稀疏矩阵等）。

目前，DMol3是有些不足的地方，因为商业化不公开源代码导致一些功能比较受限也没办法修改，比如波函数文件是不开放的，从而没法做ELF（但grd开放，可以分析AIM）。

所以，如果您真的认为那个程序不准确建议从每一步的模拟过程和计算精度仔细验证之后再下定论也不迟，当然，我们清楚地了解计算中的近似就可以避开这些“不准”（如DMol3取得是外推能量的话，展宽就不能设定的过大）。

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