求助：如何使用CP2K计算半导体功函数和fermi能级

Ahey · 发表于 Post on 2023-2-15 16:32:27

本帖最后由 Ahey 于 2023-2-17 13:43 编辑

各位老师好，学生想用CP2K计算g-C3N4（真空层为15埃左右）的fermi能级和功函数，根据文献阅读，计算的思路如下：
①先通过CP2K计算，得到体系的fermi能级Ef
②产生静电势的cube文件，载入multiwfn后在Z方向取平均，得到真空静电势，作为真空能级E0
③用Ef和E0做差，得到功函数

主要是在计算fermi能级时遇到了问题，计算过程如下（输入、输出文件已上传）：

①用优化好的g-C3N4切面，然后载入multiwfn产生cp2k输入文件（PRINT_LEVEL设为high），选用HSE06泛函（用pbe收敛的波函数作为初猜），DZVP-GTH基组，采取对角化计算，仅考虑Gamma点，将体系扩大到4*4*1，并计算所有虚轨道（即ADDED_MOS -1，之前并没有设置这个，而是用的&PDOS关键词，在群里请教了sob老师，根据老师的回复，设置了计算空轨道），同时导出molden文件。

②进行计算，计算结束后在out文件中找“fermi”关键词，然后读取之；同时把molden文件载入multiwfn，在主功能0处读取HOMO-LUMO能级

然后我发现out文件中读取的fermi level严格等于multiwfn中的HOMO能级（下图中勾出），但是这个不符合其他文章中对g-C3N4的讨论和计算。

想要请问各位老师：
①上述计算功函数的流程是正确的么
②我用于计算fermi能级的输入文件是否存在问题，使得fermi level等于HOMO能级
③如何用CP2K计算出正确的fermi能级

还望各位老师不吝赐教，谢谢各位老师！

----------------------------------------------------------------------------------------------2023-2.17 12：46 更新---------------------------------------------------------------------------------------------------------------------

再次感谢各位老师赐教，昨天晚上以及今天早上学生又进行了一些测试，下面是测试内容、结果、和一些个人的理解，还请老师们指正
由于手头有一个较小的、优化好的黑磷晶胞，所以用黑磷做测试，黑磷也是一种半导体，文献中的带隙值为1.48eV

测试内容（excel表截图，图1中所示）：
所有输入文件均由multiwfn产生，smearing的设置采用multiwfn默认
①将优化好的黑磷slab模型扩胞后，分别进行HSE06(不smearing)、HSE06（smearing）、HSE06-ADMM、pbe（不smearing）、pbe（smearing）计算，除ADMM用的OT外，其他全用对角化
②每个计算任务，都计算所有空轨道（ADDED_MOS -1 ）、都导出静电势cub文件、导出molden文件，由于HSE06-ADMM产生的molden文件中不含轨道能量信息，因此加入&PDOS关键词，读取轨道能量
③基组：DZVP-GTH 截断 500 70
④真空能级：对于静电势cube文件，在multiwfn主功能13下，取Z方向积分平均，得到真空能级E0
⑤bandgap：对于molden文件中含有轨道能量信息的，在multiwfn中读取HOMO-LUMO轨道和gap，由于晶胞够大，姑且认为HOMO-LUMO gap就是bandgap，对于HSE06-ADMM，在输出的PDOS文件中读取gap（对于开了smearing的molden文件，有时候载入multiwfn后主功能0并不会直接输出gap，因此学生自行点击轨道查找，找到临近的两个轨道，占据数分别为1.999999和0.000001，把它们作为HOMO和LUMO，不知是否妥当，见图：2）
⑥fermi能级：不开smearing时，无论是在PDOS文件中、还是在out文件中，找到的fermi能级都是精确等于HOMO的；开smearing后，out中输出的fermi能级就在gap之中，详见excel表截图中
⑦功函数：用E0-Ef得到功函数，由于不开smearing的fermi能级等于HOMO，所以功函数显著高估

结论：对于我的体系
①HSE06不ADMM、HSE06-ADMM（PDOS中读取）、HSE06-smearing计算的gap基本上一样，并且比较准；pbe开和不开smearing得到的gap也基本上一样
②PBE显著低估带隙，但是真空静电势给的比较准，开smearing输出的fermi能级也接近于HSE06
③开smearing之后fermi能级在bandgap之间，是符合物理直觉的。
如有错误或者疏漏，还请老师们指正

还有几个问题想请教老师：
①根据老师所说，一般情况下fermi能级对温度不敏感，那么以后在计算fermi能级时，是否只需要打开smearing，电子温度为multiwfn默认就可以
②如上所示，开smearing和不开smearing的结果基本上差不多，那么以后结构优化等任务，是不是也可以开smearing进行计算呢；我之前看书的时候（DENSITY FUNCTIONAL THEORY：A Practical Introduction）中了解到smearing是专门为金属准备的，那么半导体，尤其是g-C3N4和黑磷这种带隙比较宽的，开smearing是不是也是合理的呢
谢谢各位老师！

chen0201 · 发表于 Post on 2023-6-20 19:40:12

sobereva 发表于 2023-6-20 18:56
当前版本Multiwfn没有晶胞是正交的要求
诸如你计算Z方向平面平均静电势，只要a、b晶轴方向与Z轴垂直即可 ...

明白了，谢谢老师回复！

sobereva · 发表于 Post on 2023-6-20 18:56:52

chen0201 发表于 2023-6-20 15:12
老师，请问一下，您之前在帖子说：当前的cube文件的三个平移矢量不是恰好是在x,y,z方向，即格点不是立方 ...

当前版本Multiwfn没有晶胞是正交的要求
诸如你计算Z方向平面平均静电势，只要a、b晶轴方向与Z轴垂直即可，不要求a、b之间正交

sobereva · 发表于 Post on 2023-6-20 18:55:13

chen0201 发表于 2023-6-19 19:47
请问一下，您说的产生静电势的cube文件，载入multiwfn后在Z方向取平均，得到真空静电势，是输入命令13，18 ...

是
和cube文件里用的单位一致。对于CP2K，是a.u.

chen0201 · 发表于 Post on 2023-6-20 15:12:50

sobereva 发表于 2023-2-17 02:29
CP2K不同情况下给出的费米能级不一样，需根据理论知识判断什么时候给出的有意义。
开发者在CP2K的Google ...

老师，请问一下，您之前在帖子说：当前的cube文件的三个平移矢量不是恰好是在x,y,z方向，即格点不是立方的，这种情况cube文件虽然能载入Multiwfn，但是数值不正常。而楼主的晶胞也不是立方的，为什么依然可以计算真空静电势呢？

chen0201 · 发表于 Post on 2023-6-19 19:47:11

本帖最后由 chen0201 于 2023-6-19 19:54 编辑

请问一下，您说的产生静电势的cube文件，载入multiwfn后在Z方向取平均，得到真空静电势，是输入命令13，18，Z吗？得到的如下图片，数值单位是多少呢？

sobereva · 发表于 Post on 2023-3-23 07:45:29

chen0201 发表于 2023-3-14 20:01
请问cp2k可以这样算功函数吗？正确吗？

用Multiwfn基于slab模型得到真空区的静电势，减去按我说的方式确定的费米能级，是标准的得到功函数的做法

sobereva · 发表于 Post on 2023-3-23 07:43:32

Eudaimonia 发表于 2023-3-9 19:44
请问一下你的真空能级是怎么得到的
我尝试了一下使用Multiwfn对cube文件进行处理，选择13号功能下的18，对 ...

对静电势格点数据应当在Multiwfn的此功能里选择绘制平面平均曲线，而不是绘制积分曲线。什么叫平面平均曲线下文有明确说明
使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线
http://sobereva.com/638（http://bbs.keinsci.com/thread-28225-1-1.html）

chen0201 · 发表于 Post on 2023-3-14 20:01:59

请问cp2k可以这样算功函数吗？正确吗？

ykhuang · 发表于 Post on 2023-3-13 14:40:09

功函数print V_hartree.cube，然后做二维积分，工具应该是cubecruncher.x，近似真空区为真空能级=无限远，减去Efermi. 这部分内容在XPS原理

Eudaimonia · 发表于 Post on 2023-3-9 19:44:51

请问一下你的真空能级是怎么得到的
我尝试了一下使用Multiwfn对cube文件进行处理，选择13号功能下的18，对Z方向的所有范围进行积分，之后点1得到积分曲线，但是纵轴上的数据和out文件中输入的能量对应不太上，应该如何取真空能级呢？

sobereva · 发表于 Post on 2023-3-1 16:13:29

汪杰发表于 2023-2-28 18:57
请教一下老师，物理学中有关于本征半导体的费米能级的推导，结果是费米能级应该位于1/2处，并且受到电子 ...

我前面回帖里已经清楚说了费米能级是怎么确定的。除非DOS相对于费米能级是对称分布的（事实上不是），否则费米能级不恰好在正中央。说在正中央的只是个近似而已

汪杰 · 发表于 Post on 2023-2-28 18:57:35

sobereva 发表于 2023-2-17 02:29
CP2K不同情况下给出的费米能级不一样，需根据理论知识判断什么时候给出的有意义。
开发者在CP2K的Google ...

请教一下老师，物理学中有关于本征半导体的费米能级的推导，结果是费米能级应该位于1/2处，并且受到电子与空穴的有效质量还有温度的影响（公式就不打了），但是对于本征半导体来说，电子与空穴都是成对产生的，也就是说，无论温度多少，费米能级都应该在禁带中央。
上面计算的材料，按照道理来说都是本征半导体，没有掺杂等因素，为什么会有费米能级偏离中央的情况？是因为软件的算法导致的吗？我对软件不太了解，但是从最基本的理论结果来看，软件给出的结果似乎偏离了理论。为什么偏离了这个理论，是我目前想不透的地方，希望老师能教一下

Ahey · 发表于 Post on 2023-2-17 20:32:44

卡开发发发表于 2023-2-17 17:37
原则上说展宽不宜太大，此时情况才不是很敏感，大了其实也会产生影响。展宽本身确实主要是针对金属的情况 ...

好的，谢谢老师！

卡开发发 · 发表于 Post on 2023-2-17 17:37:21

Ahey 发表于 2023-2-17 13:39
感谢各位老师们的回答，为了加深理解，我在昨晚和今天早上做了一些计算和测试，得出了一些结论，更新放在顶 ...

原则上说展宽不宜太大，此时情况才不是很敏感，大了其实也会产生影响。展宽本身确实主要是针对金属的情况，但不表示半导体不可以使用，尤其是要选择合理的展宽形式。按道理Gaussian型的展宽是万用一些的（专用于金属那些展宽反而不适合半导体使用），不过CP2K只支持F-D展宽，问题倒是也不大，尤其是gap宽的时候影响会更小一些。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[CP2K] 求助：如何使用CP2K计算半导体功函数和fermi能级