关于TS校正和Grimme校正在算色散作用能量时的疑问

jiewei

请问一下各位老师和同学：
最近我用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校正吧？

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

卡开发发

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

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

jiewei

卡开发发 发表于 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

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

卡开发发

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

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

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

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

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

jiewei

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

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

卡开发发

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

jiewei

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

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

卡开发发

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

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

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

卡开发发

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），我们再进行讨论。

jiewei

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

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

卡开发发

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

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

jiewei

卡开发发 发表于 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。

jiewei

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

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

卡开发发

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方程的拟合曲线，偏离的都比较惨，所以一般都不用，这个建议仅做参考吧。

jiewei

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

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

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