费米软度的计算方法 - 计算化学公社

一、简介

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

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

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

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

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

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

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

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

二、安装方法

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

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

2. 其次，需要安装bader。

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

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

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

pip install FermiSoftness

复制代码

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

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

复制代码

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

tar -zxvf FermiSoftness-1.2.0.tar.gz
cd FermiSoftness-1.2.0
pip install -e .

复制代码

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

pip install --upgrade --user ase

复制代码

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

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：

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

复制代码

修改runfs.py内的参数：

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

复制代码

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

python runfs.py

复制代码

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

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

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

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

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

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

复制代码

修改runfs.py内的参数：

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

复制代码

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

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

python runfs.py

复制代码

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

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

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

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

&DFT
......
......
&PRINT
&MO_CUBES
STRIDE 1
NHOMO -1
NLUMO -1
&END MO_CUBES
&END PRINT
&END DFT

复制代码

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

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

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

复制代码

修改runfs.py内的参数：

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

复制代码

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

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

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]条。

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

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

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