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

Ahey

各位老师好，学生想用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是不是也是合理的呢
谢谢各位老师！

jiangning198511

半导体的费米能级就是最高占据轨道，你说的和其他文献不符合是指什么？

Ahey

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

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

jiangning198511

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

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

Ahey

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

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

sobereva

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



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

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

卡开发发

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

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

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

jiangning198511

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

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

卡开发发

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

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

jiangning198511

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

期待大佬的消息

Ahey

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

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

Ahey

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

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

jiangning198511

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

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

Ahey

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

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

sobereva

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极限），因此势必要涉及空轨道的占据，若不计算出空轨道也自然就确定不了真正意义的费米能级。