计算化学公社

标题: 求助：如何使用CP2K计算半导体功函数和fermi能级 [打印本页]

作者
Author: Ahey 时间: 2023-2-15 16:32
标题: 求助：如何使用CP2K计算半导体功函数和fermi能级
本帖最后由 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是不是也是合理的呢
谢谢各位老师！

作者
Author: jiangning198511 时间: 2023-2-15 17:25
半导体的费米能级就是最高占据轨道，你说的和其他文献不符合是指什么？

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

jiangning198511 发表于 2023-2-15 17:25
半导体的费米能级就是最高占据轨道，你说的和其他文献不符合是指什么？

老师好，首先感谢老师回答，我入门第一性原理时间不长，如果有什么概念性问题还请老师指正，我看到很多文献计算出的fermi能级都是在bandgap内的，并且我之前用QE计算的时候fermi能级也在带隙中（输入输出文件和dos文件如下，fermi能级为-3.7eV左右，禁带在-5~-2eV左右，确实在禁带内），同时之前学半导体物理时老师也讲过本征半导体fermi能级处在禁带内，所以感觉cp2k计算出的fermi能级很反直觉

作者
Author: jiangning198511 时间: 2023-2-15 17:54

Ahey 发表于 2023-2-15 17:28
老师好，首先感谢老师回答，我入门第一性原理时间不长，如果有什么概念性问题还请老师指正，我看到很多文 ...

这是Si的DOS结果，VASP计算的，费米能级在价带顶

作者
Author: Ahey 时间: 2023-2-15 18:32

jiangning198511 发表于 2023-2-15 17:54
这是Si的DOS结果，VASP计算的，费米能级在价带顶

再次感谢老师回复，也就是说可以认为HOMO能级就是Fermi能级么

作者
Author: sobereva 时间: 2023-2-16 02:23
半导体的费米能级不是HOMO（或者HOCO、价带顶）。这是常见的误区。只有导体在0K极限的情况下费米能级才是HOMO。
物理意义正确的确定方式是CP2K开smearing时用的确定方式：

(, 下载次数 Times of downloads: 8)

另外，最新版multiwfn也可以使用以CP2K开smearing时相同的方式严格按照费米能级的定义（参见wiki）计算费米能级，更为灵活，看最新版Multiwfn手册3.300.9节的详细解释。

一些文献计算非导体时还直接把HOMO能级当费米能级，这是明显错误的。

作者
Author: 卡开发发 时间: 2023-2-16 08:12
本帖最后由卡开发发于 2023-2-16 08:16 编辑

1、半导体物理书上本征半导体的Fermi能级确实是会位于价带顶和导带底之间（近1/2处），当中考虑了温度引起的电子-空穴激发并且使用了一些简化。计算当中为了方便只是简单将体系按照金属那样，对能级简单做展宽而已。

2、一般情况态密度还是讨论零温的，所以通常讨论不会去刻意加“零温”去修饰，此时希望展宽影响小一些，这样即便按价带顶去处理差异其实很小。

3、按上面讨论，若以50%占据作为轨道占据数求和的上限作为Fermi能级，Fermi能级又位于价带和导带之间，这样超出Fermi能级的轨道也有分数占据，逻辑上说仅对Fermi能级以下的轨道求和可能会不等于电子数。

作者
Author: jiangning198511 时间: 2023-2-16 09:28

卡开发发发表于 2023-2-16 08:12
1、半导体物理书上本征半导体的Fermi能级确实是会位于价带顶和导带底之间（近1/2处），当中考虑了温度引起 ...

有个问题，楼主用CP2K计算出的费米能级与HOMO一样，而用QE计算的费米能级处在能隙，其原因是两种软件的费米能级的定义不同吗？

作者
Author: 卡开发发 时间: 2023-2-16 11:09

jiangning198511 发表于 2023-2-16 09:28
有个问题，楼主用CP2K计算出的费米能级与HOMO一样，而用QE计算的费米能级处在能隙，其原因是两种软件的费 ...

QE那个我不知道啥原因，前面写了个脚本抓了一下确实如上所说，按道理大部分计算程序不会这样定义。有时间我去翻一下QE源码算了

。

作者
Author: jiangning198511 时间: 2023-2-16 11:30

卡开发发发表于 2023-2-16 11:09
QE那个我不知道啥原因，前面写了个脚本抓了一下确实如上所说，按道理大部分计算程序不会这样定义。有时间 ...

期待大佬的消息

作者
Author: Ahey 时间: 2023-2-16 11:58
本帖最后由 Ahey 于 2023-2-16 12:01 编辑

sobereva 发表于 2023-2-16 02:23
半导体的费米能级不是HOMO（或者HOCO、价带顶）。这是常见的误区。只有导体在0K极限的情况下费米能级才是HO ...

感谢老师的指导，请问学生是否可以这样理解：也就是说，
①我在上述的计算中，没有采用smearing，所以CP2K中输出的fermi能级是HOMO，并非真实的fermi能级；
②在CP2K中计算真实fermi能级需要开smearing才可以得到
③第一性原理软件中计算出的fermi能级事实上是零温下的fermi能级，温度变化fermi能级也会变化

作者
Author: Ahey 时间: 2023-2-16 12:04

卡开发发发表于 2023-2-16 08:12
1、半导体物理书上本征半导体的Fermi能级确实是会位于价带顶和导带底之间（近1/2处），当中考虑了温度引起 ...

学生理解了，看来之前还是对fermi能级有一知半解的地方，我再去好好学习，谢谢老师的回答

作者
Author: jiangning198511 时间: 2023-2-16 14:46

Ahey 发表于 2023-2-16 11:58
感谢老师的指导，请问学生是否可以这样理解：也就是说，
①我在上述的计算中，没有采用smearing，所以CP ...

1.你可以通过加smearing试试是不是对fermi能级有变化(估计温度很高才能看出差别)，但这个不影响对材料电子结果的影响
2.你这个体系算杂化是不是非常消耗内存，我这边试了一下直接是内存错误

作者
Author: Ahey 时间: 2023-2-16 19:13

jiangning198511 发表于 2023-2-16 14:46
1.你可以通过加smearing试试是不是对fermi能级有变化(估计温度很高才能看出差别)，但这个不影响对材料电 ...

谢谢老师，我今晚连夜就试试；确实很耗内存，我算的时候128G内存基本上全占了，还在输入文件中限制了每个核的内存，才没有报错

作者
Author: sobereva 时间: 2023-2-17 02:29

jiangning198511 发表于 2023-2-16 09:28
有个问题，楼主用CP2K计算出的费米能级与HOMO一样，而用QE计算的费米能级处在能隙，其原因是两种软件的费 ...

CP2K不同情况下给出的费米能级不一样，需根据理论知识判断什么时候给出的有意义。
开发者在CP2K的Google group说过
If you use diagonalisation and smearing it is the location of the chemical potential to get the correct number of electrons. For no smearing it is just given as the HOMO energy. At zero temperature it is not really defined in that case (anywhere between HOMO and LUMO would give same result).

实际中（并且经过我测试验证），仅对于不靠ADDED_MOS方式计算空轨道也不考虑k点的情况，用&PRINT输出轨道信息时，程序才简单、随意地把HOMO能级显示成费米能级。凭理论知识可知此时的这个值是没有实际意义的，因为原理上不可能在没有空轨道信息的前提下确定有意义的费米能级。费米能级是基于Fermi-Dirac分布在有限温度的前提下定义的（连0K也不能说，顶多说0K极限），因此势必要涉及空轨道的占据，若不计算出空轨道也自然就确定不了真正意义的费米能级。

作者
Author: sobereva 时间: 2023-2-17 02:36

Ahey 发表于 2023-2-16 11:58
感谢老师的指导，请问学生是否可以这样理解：也就是说，
①我在上述的计算中，没有采用smearing，所以CP ...

不同程序的规则不一样，同一个程序在不同情况下的输出也可能不一样，所以不能用“第一性原理软件中...”这么粗略地去说。只能说费米能级在物理上有确切的定义，见https://en.wikipedia.org/wiki/Fermi_level。
由于实际体系的态密度相对于费米能级不是对称分布的，所以费米能级是依赖于温度的。
CP2K里开smearing时可以根据费米能级严格的定义去给你输出当前你指定的smearing温度下对应的费米能级。通常在温度不是很高的情况下，费米能级对温度的依赖性不是很强。另外，如果基于CP2K给出的轨道能级信息用Multiwfn算费米能级的话，即便CP2K里计算时没用smearing也没关系。

作者
Author: 卡开发发 时间: 2023-2-17 05:14

Ahey 发表于 2023-2-16 11:58
感谢老师的指导，请问学生是否可以这样理解：也就是说，
①我在上述的计算中，没有采用smearing，所以CP ...

1、零温下基态情况把电子占据最高能级作为Fermi能级是没问题的（这是按照wiki那个链接的定义从化学势的角度看，电子气零温下能严格这么定义），此时你也可以理解为F-D在零温退化到阶跃函数（阶跃处值为1/2）。

2、对于CP2K，输出Fermi能级可能要通过展宽，但并不是说有温度的情况才能定义Fermi能级，其他有些程序中计算Fermi能级可以不使用展宽。展宽大多数情况不是为了得到温度的影响（而且也不是严格的），所以上面提到一般情况而言基态态密度的计算我们应该希望展宽对实际的计算尽可能小，才出现F-D之外形式的展宽或者不使用展宽的一些方案。

3、Fermi能级会随着展宽形式和参数而变化（F-D只是其中一种特例），不只是Fermi能级，电子结构也会受到些许影响。

作者
Author: 卡开发发 时间: 2023-2-17 05:39

jiangning198511 发表于 2023-2-16 14:46
1.你可以通过加smearing试试是不是对fermi能级有变化(估计温度很高才能看出差别)，但这个不影响对材料电 ...

1、电子结构会受到影响，但其实较低的电子温度也不会那么敏感，额外提及展宽形式也会有影响。
2、他的k点我看了下还算比较密的，再一个就是真空也会消耗不少计算量。我不确定算出来的晶胞就是这样还是啥原因，按道理g-C3N4存在更小的原胞。所以可能改进的方向可能是这样：
（1）可以使用更小的原胞（如果真的不重构的话）。
（2）使用截断库伦或者其他方法处理静电作用，一定程度减少真空层使用。
（3）k点和空轨道可以测试下是不是一定要这么大。

作者
Author: jiangning198511 时间: 2023-2-17 09:25

卡开发发发表于 2023-2-17 05:39
1、电子结构会受到影响，但其实较低的电子温度也不会那么敏感，额外提及展宽形式也会有影响。
2、他的k ...

1.请问这里所说的对电子结构的影响具体指什么？我所提及的对电子结构的不影响是说对于体系而言，费米能级的不同定义不影响其从能带角度判断其为导体，半导体或者绝缘体。比如LZ提到的例子，虽然两种软件给出的费米能级的位置不同(一个是价带顶，一个在能隙间)，但应该不影响其Gap的大小
2. 在什么情况下需要考虑电子温度的影响？，在常温下（300K）时，电子温度的影响应该非常小

作者
Author: 卡开发发 时间: 2023-2-17 09:53

jiangning198511 发表于 2023-2-17 09:25
1.请问这里所说的对电子结构的影响具体指什么？我所提及的对电子结构的不影响是说对于体系而言，费米能级 ...

1、如果你单纯平移能级，那确实无所谓。如果你指的是已经拿到一个自洽过程当中没使用展宽方法或者使用了很小的展宽的能级再去进行展宽画DOS，那也是没影响的。但是自洽计算阶段不同展宽大小会影响到体系的电子结构，甚至影响力的计算，过大的展宽值（也看具体展宽形式）对于窄带结构就可能被算成金属。

2、确实，一般300K时候展宽的影响是比较小的。

作者
Author: Ahey 时间: 2023-2-17 13:39
感谢各位老师们的回答，为了加深理解，我在昨晚和今天早上做了一些计算和测试，得出了一些结论，更新放在顶楼了，还请老师们指正

作者
Author: Ahey 时间: 2023-2-17 13:49

卡开发发发表于 2023-2-17 05:39
1、电子结构会受到影响，但其实较低的电子温度也不会那么敏感，额外提及展宽形式也会有影响。
2、他的k ...

感谢老师指导，我再去测试一下QE的计算，当时计算的时间确实很长，48核算了3天，也怪我当时没有经验，盲目地取非常密的k点，晶胞也取得不是原胞

作者
Author: Ahey 时间: 2023-2-17 13:50

sobereva 发表于 2023-2-17 02:36
不同程序的规则不一样，同一个程序在不同情况下的输出也可能不一样，所以不能用“第一性原理软件中...” ...

感谢老师指导，我这就去研读一下multiwfn的手册

作者
Author: 卡开发发 时间: 2023-2-17 17:37

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

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

作者
Author: Ahey 时间: 2023-2-17 20:32

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

好的，谢谢老师！

作者
Author: 汪杰 时间: 2023-2-28 18:57

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

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

作者
Author: sobereva 时间: 2023-3-1 16:13

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

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

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

作者
Author: ykhuang 时间: 2023-3-13 14:40
功函数print V_hartree.cube，然后做二维积分，工具应该是cubecruncher.x，近似真空区为真空能级=无限远，减去Efermi. 这部分内容在XPS原理

作者
Author: chen0201 时间: 2023-3-14 20:01
请问cp2k可以这样算功函数吗？正确吗？

作者
Author: sobereva 时间: 2023-3-23 07:43

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）

作者
Author: sobereva 时间: 2023-3-23 07:45

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

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

作者
Author: chen0201 时间: 2023-6-19 19:47
本帖最后由 chen0201 于 2023-6-19 19:54 编辑

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

作者
Author: chen0201 时间: 2023-6-20 15:12

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

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

作者
Author: sobereva 时间: 2023-6-20 18:55

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

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

作者
Author: sobereva 时间: 2023-6-20 18:56

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

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

作者
Author: chen0201 时间: 2023-6-20 19:40

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

明白了，谢谢老师回复！

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