计算化学公社

标题: 吸附能应该以单胞的结果还是超胞上的为准？ [打印本页]

作者
Author: lwj 时间: 2016-11-8 09:21
标题: 吸附能应该以单胞的结果还是超胞上的为准？
我用MS中的DMol3做Co3O4（001）表面吸附CO2，先在一个晶胞表面吸附CO2，由于放置位置和显示模型的问题，认为没吸附上，又做了2*2的超胞模拟吸附，后来看了版主的文章（显示存在明显的距离），比较了分子自于表面的距离，知道已经吸附在表面了，计算了表面能，发现一个晶胞上吸附能为+0.748eV，而超胞上吸附能为-4eV。这是由于吸附浓度对吸附能的影响，一个晶胞上吸附一个分子过饱和了，所以吸附能为正，那么计算吸附能时，应该以多大的晶胞吸附一个分子计算结果为准呢？请老师指点一下，我应该以多大晶胞计算吸附能为准，还是继续增大晶胞到吸附能收敛，得到干净表面的吸附能。

作者
Author: lc00 时间: 2016-11-8 09:36
表面吸附的计算，选取的表面肯定不能太小，不然不能完全代表活性位点。选取太大的话，因为吸附本身就比较耗资源，计算比较困难。表面选取的大小，一般根据文献和自己的计算能力来定。

作者
Author: lc00 时间: 2016-11-8 09:42
根据我看的文献，就算使用不同大小的表面模型，小簇或者大的slab模型，吸附能也不会有这么大的差距。只能等大神来解释了。

作者
Author: lwj 时间: 2016-11-8 09:44

lc00 发表于 2016-11-8 09:36
表面吸附的计算，选取的表面肯定不能太小，不然不能完全代表活性位点。选取太大的话，因为吸附本身就比较耗 ...

谢谢，我看有的文献就直接用一个晶胞做吸附，所以最开始就用的一个

作者
Author: lwj 时间: 2016-11-8 09:56

lc00 发表于 2016-11-8 09:42
根据我看的文献，就算使用不同大小的表面模型，小簇或者大的slab模型，吸附能也不会有这么大的差距。只能等 ...

您这么说提醒了我，我再对能量用更大的K点试试，看是不是体系能量偏的太远。谢谢您了

作者
Author: 卡开发发 时间: 2016-11-8 11:05

lwj 发表于 2016-11-8 09:56
您这么说提醒了我，我再对能量用更大的K点试试，看是不是体系能量偏的太远。谢谢您了

Co3O4（001）吸附CO2的话，我估计单胞的晶格尺寸应该已经足够了。

吸附能算出来是正的话，建议你不妨检查计算参数问题。吸附物和载体的截断半径设置的值、smearing不合理之类的都可能导致吸附能算出来是正的。

作者
Author: lwj 时间: 2016-11-8 14:37

卡开发发发表于 2016-11-8 11:05
Co3O4（001）吸附CO2的话，我估计单胞的晶格尺寸应该已经足够了。

吸附能算出来是正的话，建议你不妨 ...

谢谢您，我检查了一下参数，吸附能为正的可能是因为按文献去真空层为13，后面再算一个单胞的25A真空层试试。谢谢。

作者
Author: lwj 时间: 2016-11-12 17:36

卡开发发发表于 2016-11-8 11:05
Co3O4（001）吸附CO2的话，我估计单胞的晶格尺寸应该已经足够了。

吸附能算出来是正的话，建议你不妨 ...

老师，您好，我用MS8.0 DMol3，请问模拟表面吸附带电阴离子如何设置电荷？是modify-charge还是设置磁矩的那栏的formal charge？还是计算里边设置？谢谢您。

作者
Author: 卡开发发 时间: 2016-11-12 20:24

lwj 发表于 2016-11-12 17:36
老师，您好，我用MS8.0 DMol3，请问模拟表面吸附带电阴离子如何设置电荷？是modify-charge还是设置磁矩的 ...

DMol3 Calculation->charge，不过外加电荷始终要进行修正，否则即便进行超胞外推吸附能也没法收敛。

作者
Author: lwj 时间: 2016-11-12 20:30

卡开发发发表于 2016-11-12 20:24
DMol3 Calculation->charge，不过外加电荷始终要进行修正，否则即便进行超胞外推吸附能也没法收敛。

谢谢您了。

作者
Author: lwj 时间: 2016-11-15 14:28

卡开发发发表于 2016-11-12 20:24
DMol3 Calculation->charge，不过外加电荷始终要进行修正，否则即便进行超胞外推吸附能也没法收敛。

老师，您好，我测试真空层厚度，增加真空层厚度，在厚度取30A时，同样垂直表面的初始结构，优化后由直立变成了斜着的（25A及以下优化后直立），吸附能比25A增加70%多，是不是之前的真空层不够，受另一层的作用才形成直立吸附？另外请问，我用DMol3，后续还有要计算TS，真空层取到多大之后就最好就不再加了？

作者
Author: 卡开发发 时间: 2016-11-15 17:11

lwj 发表于 2016-11-15 14:28
老师，您好，我测试真空层厚度，增加真空层厚度，在厚度取30A时，同样垂直表面的初始结构，优化后由直立 ...

你的分子尺寸多大？

作者
Author: lwj 时间: 2016-11-15 17:33

卡开发发发表于 2016-11-15 17:11
你的分子尺寸多大？

CO2，C-O键长1.17A，大概3.5A吧

作者
Author: 卡开发发 时间: 2016-11-15 18:44
本帖最后由卡开发发于 2016-11-15 18:45 编辑

lwj 发表于 2016-11-15 17:33
CO2，C-O键长1.17A，大概3.5A吧

按说差别不应该这么大，要么是两次的计算参数或者初始结构导致的，25A对于这种小分子其实很充分了，即便有偶极修正的问题，能量计算的误差也不会这么大。

作者
Author: lwj 时间: 2016-11-15 20:22

卡开发发发表于 2016-11-15 18:44
按说差别不应该这么大，要么是两次的计算参数或者初始结构导致的，25A对于这种小分子其实很充分了，即便 ...

我刚算了一个28A的，相同的添加分子方式，优化后接近躺着吸附，斜得比较厉害，吸附能与30A比较接近，但是更高，正在算35A的，现在比较纠结如何确定真空层厚度，原本以为增加之后能量能收敛，但实际能量呈增大的趋势，表面能有一个反常，已经检查过参数，还未发现问题，看来很多东西的在检查一下。

作者
Author: 卡开发发 时间: 2016-11-15 21:22

lwj 发表于 2016-11-15 20:22
我刚算了一个28A的，相同的添加分子方式，优化后接近躺着吸附，斜得比较厉害，吸附能与30A比较接近，但是 ...

加上偶极修正看一下，如果还不行把东西打包传来。

作者
Author: lwj 时间: 2016-11-15 22:44

卡开发发发表于 2016-11-15 21:22
加上偶极修正看一下，如果还不行把东西打包传来。

我一直都是用了偶极修正。。

作者
Author: lwj 时间: 2016-11-16 15:58

卡开发发发表于 2016-11-15 21:22
加上偶极修正看一下，如果还不行把东西打包传来。

算了35A的，结果仍然存在问题，我先整理一下数据，明天向老师汇报，汇报之后整理一下打包。

作者
Author: lwj 时间: 2016-11-17 14:18

lwj 发表于 2016-11-15 22:44
我一直都是用了偶极修正。。

我今天和老师讨论后，我想可能是因为真空层增加后吸附结构变化使结果没有对比性，而我要测试的就是真空层厚度对两层之间以及吸附的分子之间作用影响，能不能直接添加分子，不优化直接算能量，通过体系能量变化判断真空层是否足够，这样行吗？
PS:我将分子竖着是因为这个构型是分子与另一层距离最近的结构，如果这样能保证真空层足够，那么其他吸附构型也就满足了。

作者
Author: 卡开发发 时间: 2016-11-17 19:53
本帖最后由卡开发发于 2016-11-17 19:55 编辑

lwj 发表于 2016-11-17 14:18
我今天和老师讨论后，我想可能是因为真空层增加后吸附结构变化使结果没有对比性，而我要测试的就是真空层 ...

1、可行，实际上不必每个层高下都做结构优化，比较一下单点能就行，当然，更严格的话我们可以把gradient算出来（Calculate手工改成gradient）和能量一样进行层高的收敛性测试。

2、

优化后由直立变成了斜着的（25A及以下优化后直立），吸附能比25A增加70%多，是不是之前的真空层不够，受另一层的作用才形成直立吸附？

你发来的输出文件我看了，所以很早就建议你把输出文件发上来，你可以看到25A以下你采用的k-points是:

# Kpoint definition file (intervals/offset):                                                                                     <--
Kpoints                      file    9 9 1 0.0000 0.0000 0.0000                                                                <--
00125-.kpoints

到了28A以上的时候：

# Kpoint definition file (intervals/offset): <--
Kpoints file 5 5 1 0.0000 0.0000 0.0000 <--
00128.kpoints

我想这个差异应该是由k-points不同所导致的，所以你可以看到k-points同样应当严格测试，马虎不得。
3、还有，smearing我并不建议使用到0.003以上，你可以测试一下，检查一下

-kTS_e=

是否足够小，以及multipolar expansion使用oct和hex的能量差别是否在可接受的范围内，否则尽可能选择hex。

作者
Author: lwj 时间: 2016-11-18 08:04
本帖最后由 lwj 于 2016-11-18 08:18 编辑

卡开发发发表于 2016-11-17 19:53
1、可行，实际上不必每个层高下都做结构优化，比较一下单点能就行，当然，更严格的话我们可以把gradient ...

老师，非常感谢您的分析。不过第二个是我文件传错了，我把能量做了图，25A的能量变化趋势有问题，就用更高的K点算了一次，能量差为2×10^-7Ha，整理的时候的时候选错了，参数都是一致的，所以我觉得还是用吸附能做指标方法的问题，变量太多了。谢谢您的建议。另外gradient是直接在input文件改，还是计算出能量再自己手工算？

作者
Author: 卡开发发 时间: 2016-11-18 10:14

lwj 发表于 2016-11-18 08:04
老师，非常感谢您的分析。不过第二个是我文件传错了，我把能量做了图，25A的能量变化趋势有问题，就用更 ...

一般还是看相同结构的单点，简单而且影响的因素少。我没记错的话如果input中把calculate改成gradient的话应该能输出force。你再比较一下force的收敛性能，如果也是和能量这样差别这么大的话，只能认为DMol3中的dipole correction方法的对真空层收敛速度太慢。如果是这样的话，我建议再做一个更大的真空层看看差别，不过DMol3真空层增加对计算时间的增加微乎其微。

作者
Author: lwj 时间: 2016-11-18 11:37

卡开发发发表于 2016-11-18 10:14
一般还是看相同结构的单点，简单而且影响的因素少。我没记错的话如果input中把calculate改成gradient的话 ...

改成了gradient,结果里没有force.我按昨天的方案计算了能量,我现在打算拟合一个曲线,找到dE/dr<2*10-3Ha/A(离子间作用力),不过结果都是dE/dr在10-5Ha/A数量级,甚至30A后,平均的dE/dr小于10-5这一数量级,而且变化趋势也符合得很好.另外老师然我定一个收敛标准,所以我在想从那个点出发比较好,比如变化能量占的百分比还是能量变化小于某个值就行了.另外您建议用Hex代替Oct,可能因为初始结构没有用Hex优化,计算能量250圈能量只收敛到2*10-3左右.

作者
Author: 卡开发发 时间: 2016-11-18 12:13
本帖最后由卡开发发于 2016-11-18 12:21 编辑

lwj 发表于 2016-11-18 11:37
改成了gradient,结果里没有force.我按昨天的方案计算了能量,我现在打算拟合一个曲线,找到dE/dr

哦，记错了，不能够直接得到force。但是可以从这里换算出force

DFT-D corrected derivatives:
DERIVATIVES (au)

后三列xyz是导数，可以写个简单脚本或者程序处理，算出每个原子的force，然后取max force或者rms force作为比较好了。一般python或者bash几乎都可以很简单地解决，如果实在搞不定再联系。
我一般的经验，真空层应该能量不用收敛到这么严格，你可以关注一下force的变化。

作者
Author: lwj 时间: 2016-11-18 14:19

卡开发发发表于 2016-11-18 12:13
哦，记错了，不能够直接得到force。但是可以从这里换算出force

后三列xyz是导数，可以写个简单脚本或 ...

好的，谢谢您。

作者
Author: lwj 时间: 2016-11-19 09:25
本帖最后由 lwj 于 2016-11-19 09:27 编辑

卡开发发发表于 2016-11-18 12:13
哦，记错了，不能够直接得到force。但是可以从这里换算出force

后三列xyz是导数，可以写个简单脚本或 ...

不得不像您之前说的想吐槽一下DMol3，中间加CO2后，大于15A之后体系能量变化更慢，也就是说作用力更小，有点让人想不通，不知道是不是和电荷分布变化有关。过了15A之后变化趋势还反过来了

作者
Author: 卡开发发 时间: 2016-11-19 11:14

lwj 发表于 2016-11-19 09:25
不得不像您之前说的想吐槽一下DMol3，中间加CO2后，大于15A之后体系能量变化更慢，也就是说作用力更小， ...

实在不行你就直接把真空层取得足够大，也只能认为足够大的真空得到的结果可靠罢了。

作者
Author: lwj 时间: 2016-11-19 11:25

卡开发发发表于 2016-11-19 11:14
实在不行你就直接把真空层取得足够大，也只能认为足够大的真空得到的结果可靠罢了。

谢谢您的建议，其实在15A之后能量变化在0.001eV以内（15A和35A的差）。基本认为是收敛了，只不过我有强迫症，总想弄清楚。

作者
Author: 卡开发发 时间: 2016-11-19 12:32

lwj 发表于 2016-11-19 11:25
谢谢您的建议，其实在15A之后能量变化在0.001eV以内（15A和35A的差）。基本认为是收敛了，只不过我有强迫 ...

但是现在确实做出来吸附构型不同？我想还是检查一下force会比较保险。

作者
Author: lwj 时间: 2016-11-19 15:48

卡开发发发表于 2016-11-19 12:32
但是现在确实做出来吸附构型不同？我想还是检查一下force会比较保险。

好的，我试一下。

作者
Author: lwj 时间: 2016-11-19 16:05

卡开发发发表于 2016-11-18 12:13
哦，记错了，不能够直接得到force。但是可以从这里换算出force

后三列xyz是导数，可以写个简单脚本或 ...

老师，请问这个 DERIVATIVES (au)下面的x、y、z对应的是能量的导数还是力的导数（df）？

作者
Author: lwj 时间: 2016-11-19 17:04

卡开发发发表于 2016-11-19 12:32
但是现在确实做出来吸附构型不同？我想还是检查一下force会比较保险。

检查了一下结果，只有25A时力也收敛，其他的20，28，30，35，力都没收敛。看来25A是比较好的。

作者
Author: 卡开发发 时间: 2016-11-19 18:50

lwj 发表于 2016-11-19 17:04
检查了一下结果，只有25A时力也收敛，其他的20，28，30，35，力都没收敛。看来25A是比较好的。

不是这个意思，必然是z->∞的时候，force->定值，这样才叫做测试收敛。这个导数指的是每个原子能量对核坐标的导数，也就是力的分量。

作者
Author: lwj 时间: 2016-11-20 09:56
本帖最后由 lwj 于 2016-11-20 10:06 编辑

卡开发发发表于 2016-11-19 18:50
不是这个意思，必然是z->∞的时候，force->定值，这样才叫做测试收敛。这个导数指的是每个原子能量对核坐 ...

这个我算了两个，最大受力与方框里的很接近，其他的还没试，但是优化结果判断时，每个原子受力小于收敛条件的就25A时。

作者
Author: 卡开发发 时间: 2016-11-20 12:50

lwj 发表于 2016-11-20 09:56
这个我算了两个，最大受力与方框里的很接近，其他的还没试，但是优化结果判断时，每个原子受力小于收敛条 ...

你再好好理解我上面说的意思，包括k点和截断能之类的测试，都是建立在外推至无穷大的概念上。真空层也不例外，这是因为如果slab的上下表面没有电势差，我们认为随着两个层距离->∞的情况相互作用->0，才有这样的测试。

作者
Author: lwj 时间: 2016-11-20 13:17

卡开发发发表于 2016-11-20 12:50
你再好好理解我上面说的意思，包括k点和截断能之类的测试，都是建立在外推至无穷大的概念上。真空层也不 ...

哦，是不是意思是力收敛是结构稳定的要求，而较厚的真空层才是模拟表面的理论根本。

作者
Author: rainwatcher 时间: 2017-3-20 14:45
记得有文献提过，如果分子之间的作用较强，或是有明显的氢键作用，那么吸附密度应该较小，这样能够减小分子之间的作用对吸附能的影响。如果分子较小，特别是气体分子，使用单胞计算也是可以的。

作者
Author: lwj 时间: 2017-3-23 13:57

rainwatcher 发表于 2017-3-20 14:45
记得有文献提过，如果分子之间的作用较强，或是有明显的氢键作用，那么吸附密度应该较小，这样能够减小分子 ...

好的，谢谢。

作者
Author: 缠绕指 时间: 2017-4-2 16:01
不同的晶胞大小代表不同的覆盖率下的吸附能，这两者算的吸附能都有意义。

一般不注明吸附率时的吸附能说的都是覆盖率足够低（晶胞足够大）时的吸附能，也就是晶胞大到吸附能收敛时的吸附能。

不过你算的这个+0.748和-4，这个差距太大了，应该是你自己计算错误，而不是我上面说的覆盖率的影响。

作者
Author: lwj 时间: 2017-4-6 17:10

缠绕指发表于 2017-4-2 16:01
不同的晶胞大小代表不同的覆盖率下的吸附能，这两者算的吸附能都有意义。

一般不注明吸附率时的吸附能说 ...

谢谢您的回答，我为了后面的计算，计算体系稍大一些，现在结果都正常了。

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