计算化学公社

标题: 用CP2K做过渡金属001表面HER台阶图时遇到的数值问题 [打印本页]

作者
Author: Eudaimonia 时间: 2023-2-15 10:26
标题: 用CP2K做过渡金属001表面HER台阶图时遇到的数值问题
得到的结果和文献报道的结果有非常大的不同，经过人为粗暴的数值校正后，看起来比较类似；是否计算设置相关的精度需要进行进一步的提升？
体系是Ag Au Cu Pb Pd Pt Ru Ti V W金属的001面，均使用221的超晶胞，Z轴方向原子层数＞5以保证上表面接近实际表面，真空层厚度设置为15埃，k点选取为331，进行D3-BJ校正，打开smearing，其他设置均为Multiwfn 3.8 dev默认参数

进行H原子吸附计算的时候，体系自旋多重度为2，打开UKS，其他设置与001 slab一致；

对不同位点的H吸附计算完毕后选取能量最低的模型，固定与H原子距离＞3埃的所有原子，与几何优化参数计算相同，但是由于设置k点会导致无法计算部分数据导致频率计算不能正常完成，在freq中使用Gamma点进行计算

计算频率后输出文件使用Shermo 2.3.5提取ΔG，参数为默认参数，未经任何更改(将输入文件用Shermo处理的时候会提示周期性体系中imode需要设置为1，测试了一下是否设置为1会带来约0.05eV能量的误差，对本文主要问题差至少1个数量级，故忽略)

H吸附能定义为Gads = Eele + ΔG Note：cp2k的频率计算如果固定了原子，输出结果中会检测不到电子能量，所以直接使用几何优化得到的 Total energy作为电子能量

之后吸脱附能垒，根据Gads = -neU进行换算，得到热力学上的起始电位U，以此电位为y轴，对Gads作图得到火山图，注意此处Gads和U是通过公式计算出来的，所以围成的曲线必为y = ±ax，对于n使用1或者2，只会影响a的大小，但是对于金属的相对位置没有影响

结果如下作图所示：可以看到位于火山图最顶端的是Au，Ag，与文献报道并经过实验验证的Pt最好完全不一致(文章可以找一下Nørskov的相关文献)；然后人为将所有吸附能+0.7eV(实际上是将Pt加到0eV附近)后，图片如右图所示，此时大部分金属位置与传统火山图比较类似（Note: 传统火山图一般计算111面）

之后使用相同参数计算了Pt111面的HER台阶图，能垒高达1.5eV，相差更远

在之前的帖子 [CP2K] cp2k 算吸附能和vasp相比负了大约0.8eV，还能被审稿人接受么？卡开发发老师提到”赝势和基组可能会导致对全电子结果再现性问题并不应该忽视“，并给出了参考文献10.1021/acs.jpcc.9b03554，大概浏览了一下，说使用不同基组会带来0.4eV左右的误差，CP2K计算HER火山图的误差来源是否是基组呢？

相关输入文件见附件

作者
Author: jiangning198511 时间: 2023-2-15 11:28
本帖最后由 jiangning198511 于 2023-2-15 11:42 编辑

1. 几何优化用的K点是331，频率计算采用的是伽马点计算，这个是否会带来误差
2. 这些金属的吸附能数据是否有文献报道过？
2. 文献里面Y轴用的是电流，你用的是电位，这两个可以直接对比吗？

作者
Author: kubingxia23 时间: 2023-2-15 15:40
为什么是 “固定与H原子距离＞3埃的所有原子” 呢，看输入文件似乎是还放开了表面的两个金属原子，可能这里导致你的能量校正有点问题

全部选择(001)面有什么理由吗，建议先做调整使Pt(111)的计算结果和文献相符

作者
Author: jiangning198511 时间: 2023-2-15 16:06

kubingxia23 发表于 2023-2-15 15:40
为什么是 “固定与H原子距离＞3埃的所有原子” 呢，看输入文件似乎是还放开了表面的两个金属原子，可能这里 ...

一般好像是固定slab，只放开吸附分子

作者
Author: Eudaimonia 时间: 2023-2-15 16:33

jiangning198511 发表于 2023-2-15 11:28
1. 几何优化用的K点是331，频率计算采用的是伽马点计算，这个是否会带来误差
2. 这些金属的吸附能数据是否 ...

1.这个不是很确定频率和k点之间的关系了，个人感觉很容易带来误差

2.有很多文献报道过，但是一时间可能不太好找

3.纵轴的获得方式有好几种，a.实验值 b.经由计算结果结合动力学进行模拟 c.使用能斯特方程进行校正，此处使用的是方法c，只考察趋势的话应该是可以对比的，可以参考本文的相关结论 10.1039/c1cp22271f

作者
Author: Eudaimonia 时间: 2023-2-15 16:38

kubingxia23 发表于 2023-2-15 15:40
为什么是 “固定与H原子距离＞3埃的所有原子” 呢，看输入文件似乎是还放开了表面的两个金属原子，可能这里 ...

此处是根据该帖子来设置的，主要是cp2k计算频率如果全放开计算量太大实在承受不了所以才如此进行的

作者
Author: jiangning198511 时间: 2023-2-15 16:50

Eudaimonia 发表于 2023-2-15 16:33
1.这个不是很确定频率和k点之间的关系了，个人感觉很容易带来误差

2.有很多文献报道过，但是一时间可 ...

1.这个不是很确定频率和k点之间的关系了，个人感觉很容易带来误差
因为频率计算不能使用K点，你可以都用超胞计算以保持计算参数稳定性，用一个体系分析对吸附能的影响

2.有很多文献报道过，但是一时间可能不太好找
如果是发表文章的话，还是要找一下，不然会有审稿人质疑计算合理性

3.纵轴的获得方式有好几种，a.实验值 b.经由计算结果结合动力学进行模拟 c.使用能斯特方程进行校正，此处使用的是方法c，只考察趋势的话应该是可以对比的，可以参考本文的相关结论 10.1039/c1cp22271f
HER的文章很多直接给出过电位的对比图来分析不同金属的催化活性

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