费米软度的计算方法

linqiaosong

费米软度的计算方法

（请先看第七~九节）

原帖地址：Fermi-Softness-for-VASP：VASP计算费米软度的脚本

一、简介

费米软度是武汉大学黄冰、庄林等人在Angew上提出的一种新的金属表面描述符。其思想是通过费米-狄拉克分布的导数作为一个权重函数，来计算费米能级附近的有效态密度之和。一些文章上也利用费米软度来构建一些吸附能之类的经验关系。

费米软度的空间表示称之为局域费米软度，其表达式如下，f'是费米狄拉克分布的导数：

根据和@卡开发发 老师的讨论，可以把它写作为对k点和band的求和：

简缩费米软度则通过一些原子划分方法，把局域费米软度在原子划分区域中积分，以得到简缩到每一个原子上的费米软度。

全局费米软度则是直接把空间部分积分掉，得到一个数值，也称之为总费米软度。

二、安装方法

更详细的内容请参看：[Github] Fermi-Softness-for-VASP

1. 首先需要安装Anaconda3，安装方法略去。

2. 其次，需要安装bader。

3. 对于VASP用户，额外需要安装VASPKIT。

4. 最后，需要安装FermiSoftness包：

FermiSoftness的包可以采用pip直接安装：

1. pip install FermiSoftness

复制代码

或者通过Github下载源码安装：

1. git clone https://github.com/Linqiaosong/Fermi-Softness-for-VASP.git
   
2. cd Fermi-Softness-for-VASP
   
3. pip install -e .

复制代码

也可以手动下载源码压缩包安装：下载地址1  下载地址2

1. tar -zxvf FermiSoftness-1.2.0.tar.gz
   
2. cd FermiSoftness-1.2.0
   
3. pip install -e .

复制代码

安装过程中会自动安装ASE，如果ASE安装失败，可以通过以下命令手动安装ASE：

1. pip install --upgrade --user ase

复制代码

安装完成后，通过下列命令测试是否安装正常，如果未报错，则安装正常：

1. python -c "import ase; import FermiSoftness;"

复制代码

三、VASP用户的使用方法

对于VASP用户，请确认安装VASPKIT，并且VASPKIT "51-516"功能能够正常生成波函数cube文件。

首先计算流程和计算DOS一样，通过VASP完成一次大k点的非自洽计算，并在INCAR中设置LWAVE = .TRUE.，得到WAVECAR和vasprun.xml文件。

在非自洽计算的目录下，执行下列命令生成FermiSoftness的标准输入文件runfs.py：

1. python –c “from FermiSoftness import gen; gen(software=‘vasp’);”

复制代码

修改runfs.py内的参数：

1. #-------parameters----------
   
2. kbT=0.4                            # Electron temperature (eV): recommended 0.4 by B. Huang
   
3. 
   
4. dfdd_threshold=0.001               # Derivation of Fermi-Dirac distribution threshold: recommended 0.001 by B. Huang
   
5. 
   
6. intermediate_file_options=False     # Save intermediate files? False or True (default: False)
   
7. 
   
8. bader_dir='bader'                  # Path of bader, if bader is in your $PATH, you don't need to change it
   
9. 
   
10. vaspkit_dir='vaspkit'              # Path of vaspkit, if vaspkit is in your $PATH, you don't need to change it
    
11. 
    
12. band_gap={'VBM':[0.0],             # If band gap exists (You might need to confirm the occupation of VBM and CBM):
    
13.           'CBM':[0.0]}             # non-spin polarization: set as 'VBM':[E_VBM],'CBM':[E_CBM] (Do not minus E_fermi)
    
14.                                    # spin polarization: set as 'VBM':[E_VBM_UP,E_VBM_DW],'CBM':[E_CBM_UP,E_CBM_DW]
    
15.                                    # Otherwise: set as 'VBM':[0.0],'CBM':[0.0]
    
16. #----------------------------   

复制代码

随后执行下列命令，开始费米软度的计算，计算过程中会反复调用VASPKIT生成波函数的格点文件，通常需要花费比较长的时间：

1. python runfs.py

复制代码

完成后会在当前目录下生成LFS.cube和FSCAR文件。

四、Quantum-Espresso用户的使用方法（实验性功能）

对于QE用户，请确认编译了pp.x模块。

首先，通过QE完成一次大k点的nscf计算。

在nscf输入文件的目录下，执行下列命令生成FermiSoftness的标准输入文件runfs.py：

1. python –c “from FermiSoftness import gen; gen(software=‘qe’);”

复制代码

修改runfs.py内的参数：

1. #-------parameters----------
   
2. 
   
3. prefix='pwscf'
   
4. 
   
5. outdir='./tmp'
   
6. 
   
7. kbT=0.4                            # Electron temperature (eV): recommended 0.4 by B. Huang
   
8. 
   
9. dfdd_threshold=0.001               # Derivation of Fermi-Dirac distribution threshold: recommended 0.001 by B. Huang
   
10. 
    
11. intermediate_file_options=False     # Save intermediate files? False or True (default: False)
    
12. 
    
13. bader_dir='bader'                  # Path of bader, if bader is in your $PATH, you don't need to change it
    
14. 
    
15. pp_laucher='mpirun -np 4 pp.x'                 # Laucher of pp.x, e.g.: 'pp.x' or 'mpirun -np 4 pp.x'
    
16. 
    
17. band_gap={'VBM':[0.0],             # If band gap exists (You might need to confirm the occupation of VBM and CBM):
    
18.           'CBM':[0.0]}             # non-spin polarization: set as 'VBM':[E_VBM],'CBM':[E_CBM] (Do not minus E_fermi)
    
19.                                    # spin polarization: set as 'VBM':[E_VBM_UP,E_VBM_DW],'CBM':[E_CBM_UP,E_CBM_DW]
    
20.                                    # Otherwise: set as 'VBM':[0.0],'CBM':[0.0]
    
21. #----------------------------   

复制代码

注意：prefix和outdir需要设置的和nscf输入文件中一致，如果nscf中缺省了这些关键词，也应该填入默认值。如果outdir是当前目录，则应该输入'./'，此处不可以简写为' '。另外，pp_laucher应当设置为pp.x的启动命令，并确保这个命令是可用的。

随后执行下列命令，开始费米软度的计算，计算过程中会反复调用pp.x生成波函数的格点文件，通常需要花费比较长的时间：

1. python runfs.py

复制代码

完成后会在当前目录下生成LFS.cube和FSCAR文件。

五、CP2K用户的使用方法（实验性功能）

对于CP2K的用户，由于CP2K生成波函数只支持Gamma点，所以请尽可能使用大的超胞和不太寒酸的基组。

首先要完成一次单点能的计算，计算输入文件中应当包含MO_CUBES模块关键词，并设置NHOMO和NLUMO的值为-1：

1. &DFT
   
2. ......
   
3. ......
   
4.     &PRINT
   
5.       &MO_CUBES
   
6.         STRIDE 1
   
7.         NHOMO -1
   
8.         NLUMO -1
   
9.       &END MO_CUBES
   
10.     &END PRINT
    
11.   &END DFT

复制代码

完成单点能计算后会在当前目录下生成一系列以PROJECT名来命名的PROJECT-WFN_xxxxx_x-x_x.cube波函数格点文件。

在该目录下执行以下命令，生成FermiSoftness标准输入文件runfs.py：

1. python –c “from FermiSoftness import gen; gen(software=‘cp2k’);”

复制代码

修改runfs.py内的参数：

1. #-------parameters----------
   
2. 
   
3. filename='pt3y.out'
   
4. 
   
5. project_name='pt3y'
   
6. 
   
7. ispin=1                            # If you used UKS or LSD, set as 2; Otherwise, set as 1
   
8. 
   
9. kbT=0.4                            # Electron temperature (eV): recommended 0.4 by B. Huang
   
10. 
    
11. dfdd_threshold=0.001               # Derivation of Fermi-Dirac distribution threshold: recommended 0.001 by B. Huang
    
12. 
    
13. bader_dir='bader'                  # Path of bader, if bader is in your $PATH, you don't need to change it
    
14. 
    
15. band_gap={'VBM':[0.0],             # If band gap exists (You might need to confirm the occupation of VBM and CBM):
    
16.           'CBM':[0.0]}             # non-spin polarization: set as 'VBM':[E_VBM],'CBM':[E_CBM] (Do not minus E_fermi)
    
17.                                    # spin polarization: set as 'VBM':[E_VBM_UP,E_VBM_DW],'CBM':[E_CBM_UP,E_CBM_DW]
    
18.                                    # Otherwise: set as 'VBM':[0.0],'CBM':[0.0]
    
19. #----------------------------  

复制代码

注意：filename为单点能的输出文件，project_name应当与输入文件中PROJECT关键词后的值一致。

随后执行下列命令，开始费米软度的计算，这个过程一般比较快：

1. python runfs.py

复制代码

完成后会在当前目录下生成LFS.cube和FSCAR文件。

六、输出结果和后处理

FSCAR文件是Bader的ACF.dat格式的文件，里面的CHARGE对应每个原子在Bader划分下的简缩费米软度，最后的NUMBER OF ELECTRONS则是对所有原子求和后的总费米软度。勘误：目前版本的FSCAR的简缩费米软度的数值是错误的，详见第七节。

LFS.cube文件是局域费米软度的格点文件，可以使用VMD或者VESTA等可视化软件显示，也可以做成电荷密度等值面上的投影图。

具体的操作方法在此不再赘述。

七、一些问题

首先，FermiSoftness这个包里面我使用的是Bader划分方法做的原子划分，用于计算简缩费米软度，因为黄冰在博士论文中认为AIM划分方法会更好。勘误：目前的划分方法是错误的，FSCAR实际上得到的是按局域费米软度进行Bader划分，严格的做法应该是按全电子电荷密度进行Bader划分，因此正确的做法应该是去算一个全电子的电荷密度的格点文件（假设叫chgsum.cube），然后通过bader LFS.cube -ref chgsum.cube，这样得到的简缩费米软度数值才是正确的。

其次，费米软度通常只用来描述金属或者导体这种没有带隙的体系，对于有带隙的体系，根据和@卡开发发 老师的讨论，考虑过把费米狄拉克分布的导数劈成两半，一半放在价带顶，一半放在导带底，然后分别对价带和导带部分单独积分（或对k点和band求和），得到价带的费米软度（电子）和导带的费米软度（空穴）。不过暂时没有测试这么做是否合适。如果想要体验的话，可以修改输入文件的band_gap。

第三，QE和CP2K是新支持的功能，并没有经过严密的测试，所以我也不知道有没有bug。其实VASP的可能也有我没注意的潜在bug。对于数值上的结果，和一些文献报道的对应于吸附能的经验公式，不知道能不能够吻合。在Github上面的测试结果中，Pt3Y的QE和VASP计算出来的费米软度结果是不一样的，当然输入文件是我随手写的，也没有保证他俩严格一致，结果不一样也实属正常。

最后，手头的测试版本是VASPKIT 1.3, VASP 5.4.4, QE 6.4.1, CP2K 8.2，其它版本会不会由输出文件识别错误的问题暂不明确。

八、参考文献

[1] B. Huang, L. Xiao, J. Lu, L. Zhuang, Angew. Chem. Int. Ed. 2016, 55, 6239 –6243

[2] 黄冰. 固体表面反应性理论研究.武汉大学博士学位论文,2015.

[3] B. Wang, F. Zhang, Angew. Chem. Int. Ed. 2022, 61, e202111026.

九、声明

本人专业不是搞固体的，也不是搞计算的，所以写这个代码纯属业余无聊，没有任何经济利益，代码写的也很屎。所以代码和本文都不设版权，可以随意复制转载。

有bug或者建议之类的欢迎提出，反正也不一定改。由于课业紧张，本项目不定期维护。

如果论文中使用了本工具，也没必要引用什么东西，不过出于对原作者的尊重，请引用费米软度的原文，也就是参考文献第[1]条。

有任何关于费米软度原理和结果分析方面的问题，请找论文原作者讨论，你想找我讨论也不是不行，但是我不一定答得上来。

非常感谢@卡开发发 老师的指导，冬锅牛逼就是了。

也欢迎有想完善这个功能的大佬们进一步改进这个功能，反正公式也在上面有了，也不一定要用我的代码，反正我也不设置版权问题，都是学术交流和讨论嘛，大家都可以参与。

冰释之川

松哥太强了……

wypkdhd

这是真超级的技术质量。

cliuimr

请问一下，如果是半金属的话，那个VBM和CBM应该怎么设置呀？

linqiaosong

cliuimr 发表于 2022-8-30 11:50
请问一下，如果是半金属的话，那个VBM和CBM应该怎么设置呀？

费米软度对于半金属乃至非金属而言是否适用这个问题尚且是一个未知数，原文中应该只用来充当金属体系的描述符

chen0201

老师您好，我的费米软度已经开始运行，并且运行了一天多，最后只生成了.cube格式的文档 （和vaspkit生成的好像差不多），没有LFS.cube和FSCAR文件，请问这是为什么呢?哪里出了问题呢？

linqiaosong

chen0201 发表于 2022-11-11 20:31
老师您好，我的费米软度已经开始运行，并且运行了一天多，最后只生成了.cube格式的文档 （和vaspkit生成的 ...

目测是没有正常结束，查看一下vaspkit.log文件是否有报错

chen0201

linqiaosong 发表于 2022-11-11 21:22
目测是没有正常结束，查看一下vaspkit.log文件是否有报错

谢谢老师，我又重新测试了，确实生成了LFS.cube 和FSCAR,整个计算过程有三天左右，请问正常吗？计算过程中有警告（如图），这个影响大吗？

linqiaosong

chen0201 发表于 2022-11-16 11:52
谢谢老师，我又重新测试了，确实生成了LFS.cube 和FSCAR,整个计算过程有三天左右，请问正常吗？计算过程 ...

1. 慢是很正常的，因为这东西我没写并行，是串行的，而且反复调用vaspkit读写格点，IO也会花很长时间。
2. 这个警告不是费米软度的警告，你可以检查一下你的vasp计算本身有没有不合理的地方。
3. 此外，关于vasp算出来的FSCAR，请留意我的第七节。

首先，FermiSoftness这个包里面我使用的是Bader划分方法做的原子划分，用于计算简缩费米软度，因为黄冰在博士论文中认为AIM划分方法会更好。勘误：目前的划分方法是错误的，FSCAR实际上得到的是按局域费米软度进行Bader划分，严格的做法应该是按全电子电荷密度进行Bader划分，因此正确的做法应该是去算一个全电子的电荷密度的格点文件（假设叫chgsum.cube），然后通过bader LFS.cube -ref chgsum.cube，这样得到的简缩费米软度数值才是正确的。

4. 关于vasp算出来的电荷密度格点与vaspkit的波函数格点不一致，造成无法通过bader计算简缩费米软度的情况，我无法解决，你可以咨询一下vaspkit的开发者。

chen0201

linqiaosong 发表于 2022-11-16 21:36
1. 慢是很正常的，因为这东西我没写并行，是串行的，而且反复调用vaspkit读写格点，IO也会花很长时间。
...

谢谢老师，我有注意到您说的第七节，我目前只用了局域的费米软度，LFS.cube文件，另外您说的全电子的电荷密度的格点文件是算的CHGCAR吗？

linqiaosong

chen0201 发表于 2022-11-18 15:33
谢谢老师，我有注意到您说的第七节，我目前只用了局域的费米软度，LFS.cube文件，另外您说的全电子的电荷 ...

这个具体可以查看bader的官网

chen0201

linqiaosong 发表于 2022-11-18 17:27
这个具体可以查看bader的官网

好的，谢谢老师

huluobo123

报错 **** !!!! Running vaspkit error! check the vaspkit.log !!!! ****，我检查了vaspkit路径没有问题，有人知道报错原因吗？

linqiaosong

huluobo123 发表于 2023-10-19 17:25
报错 **** !!!! Running vaspkit error! check the vaspkit.log !!!! ****，我检查了vaspkit路径没有问题， ...

看一下生成的vaspkit.log文件，很多时候是生成波函数格点太密集，然后vaspkit写波函数cube的时候内存溢出了