CP2K 使用与 MOF 材料计算

@tkaray CC-BY-NC(4.0) 2022.11.07 / 2024.10.04

CP2K 作为免费, 开源, 快速的第一性原理计算软件已经得到了广泛的应用, 但目前中文资料较少, 本文是一些个人使用 CP2K 处理 MOFs 材料问题的简略总结, 有新资料将会持续更新.

另外, 这里有一些关于杂化泛函相关的问题仅存在于 CP2K 2022 版本中, 所以可能和最新版有微小差别, 推荐使用最新版本的 CP2K, 能够减少很多麻烦.

官方链接:

cp2k website

cp2k GitHub Repository

1 一些值得关心的问题

Q1. 为什么要使用 CP2K?

MOFs 材料由金属离子和有机配体组成, 晶胞尺寸一般都超过 1000 Å^3, 原子数一般大于 200. 因此, MOFs 材料的计算对效率的要求极高, 如果单个晶胞有 1000 以上的原子, 基本也就告别轻轻松松使用 VASP 了. 还好, CP2K 效率极高, 使用 64 核心节点对 500 原子的 MOFs 结构优化, 一般也可以在 2-3 天完成, 1000 原子以上也不是很离谱 (可能要一星期以上, 想想 300 原子的 Gaussian 杂化泛函优化时间).

正因为 CP2K 的高效率, 越来越多的 MOFs 材料计算会选择 CP2K 进行从头算分子动力学 (AIMD) 的软件. 因此, 学习 CP2K 对于精准描述 MOFs 材料的性质大有帮助.

Q2. 在哪里能够学习 CP2K?

CP2K 的官网 reference manual 写的比较"鸡肋", 只对懂得如何做相关计算的人有用, 对于新手帮助不大. 但是官网的一些 workshop 值得参考, 只要在官网右上角搜索就能查到相关信息.

由于 CP2K 是开源免费的软件, 因此没有国内的销售代理, 也就基本没有人组织相关培训, 那么第三方资料就很重要. 计算化学公社 的第一性原理板块有大量的讨论内容, 卢天老师也有基于 CP2K 的第一性培训班. 除此之外, 有一些 CP2K 的中文资料也能够参考一下, 比如 https://wiki.cheng-group.net.

2 CP2K 的安装

目前这类资料已经有很多, 在 bilibili 也能搜到一些录屏教程. 个人推荐使用 singularity 版本, 一次编译在哪里都能用 (如果用别人的 sif 镜像连编译都省下了). 需要注意的是超算中的 CP2K 经常缺少功能, 无法跑通全部的测试. 即使在超算中使用, 也可以通过 singularity 方式调用自己编译好的版本.

我自己的记录可以见 CP2K 在 Docker / Singularity 容器中的安装 (v2022.1)

更新: 现在可以直接用 docker pull 的办法拉取了, 流程见链接

macOS 也可以使用 brew 链接

如果使用的是虚拟化方案, 参考以下运行脚本:

export SIF=/home/tkaray/cp2k-2023.2
export cpucores=36
singularity run -B ${PWD}:/mnt -H /mnt ${SIF} \
mpiexec -np ${cpucores} cp2k.psmp $1 |tee ${2:-cp2k.out}

建议用沙盒模式运行虚拟环境, 方便调整 OMP_NUM_THREADS.

3 MOFs 材料的计算准备

MOFs 材料需要特别注意的是 cif 文件的处理和自旋多重度的选择.

3.1 cif 文件的处理

晶体结构数据一般都会作为结构优化的初猜, 并基于优化后的结果进行后续计算. MOFs 材料中经常会解析出游离溶剂或离子, 骨架有无序结构, 这些问题都需要处理.

对于游离的溶剂分子, 如果对于目标的材料性质不重要且基本和骨架没作用, 可以去除. 但骨架若带有电荷, 游离离子则不能去除. 去除过程可以在 Olex 2, Diamond, MS 等软件中进行. 如果空腔中的阴阳离子不能解出, 也可以考虑根据合成过程在骨架中添加假的阴阳离子 (比如氯离子, 二甲胺阳离子等) 放置在和骨架不发生作用的位置, 但请注意添加原子不应改变材料的性质.

2024.10.04 更新: 对于 CP2K 来说, 直接设置骨架带电荷可能比设置离子更可靠, 但可能有些性质计算会存在问题, 需要进一步测试.

CP2K 的计算中不允许分数占有率的存在, 无序结构必须删掉部分原子保留单套结构. 删除过程也可以在前文所述软件中进行, 重新保存 cif 即可. 如果需要考察掺杂性质, 则可能需要构建超胞.

3.2 金属的处理

如前文所述, MOFs 材料均包含金属, 很多都是过渡或稀土金属, 因此需要修改自旋多重度设定. 如果材料中没有单电子, 那么自旋多重度一般为 1, 如果生成工程文件 (见下一部分) 后发现自旋多重度为 2, 那么结构肯定有问题. 如果材料中有单电子, 那么需要根据单个晶胞中单电子数目 S 设置自旋多重度为 2S+1. 比如, 材料中有 4 个 Mn(II), 那么自旋多重度就应该是 41. 对于具有高低自旋可能性的材料, 需要考虑将自旋多重度分别设定并进行能量对比, 具体情况可能更为复杂. 当金属带有单电子时, 需要在基组设置部分添加 MAGNETIZATION 项指定单电子数目.

另外, 如果需要进行 UKS 三重态结构优化, 也可以选择将自旋多重度设置为 3 并优化, 此时能够得到 T1 态的晶体结构. 晶体的 TDDFPT 和团簇的 TDDFT 计算有很大区别, 这点将在后续进一步说明.

4 输入文件的生成

4.1 快速生成可运行的 CP2K 输入文件

在详细了解输入文件内容前, 有必要先生成一个能够正常运行的文件. CP2K 的工程文件比较"反人类", 但 Multiwfn 能够有效解决该问题. 卢天老师的文章 使用Multiwfn非常便利地创建CP2K程序的输入文件 详细描述了 Multiwfn 的 CP2K 文件生成功能, 非常好用.

4.2 输入文件结构与报错检查

刚开始用 CP2K 时会觉得学习曲线非常陡峭, 但是用了一段时间后会发现该软件的输入文件非常清晰, 比起 VASP 较为平面化, 基于数字设置项的 INCAR 形式, 当需要使用功能较多时, 以多层级, 关键词的设置方式会更加清晰可靠.

CP2K 的设置文件一般只有一个, 其中同时包括了晶胞信息和工程内容. 常用部分与功能如下:

• GLOBAL: 全局设置
• FORCE_EVAL: 包含基于 DFT 能量计算的主要设置
  • DFT: 包含 DFT 控制项
  • DFT-PRINT: 包含能带结构, PDOS 等项目
  • PROPERTIES: 包含原子电荷, 线性响应, TDDFT 等性质
• MOTION: 包含结构优化, MD, MC 等内容

不清楚性质项目时, 可以在 CP2K 官网搜索 tutorial. 当程序报错退出时, 可以根据报错信息搜索源代码, 也可以在 Google 上搜索 CP2K Google Group 的讨论.

4.3 基组的选择

CP2K 的推荐基组是针对该软件专门优化的 MOLOPT 基组 (配合 GTH 赝势), 一般都可以直接使用. 当使用 Multiwfn 时, MOFs 材料的金属或重无机元素基组设置需要进行一定修改, 包括更换基组名称和改变引用项目. 方法很简单, 只需要在 CP2K 的源代码, data 文件夹中搜索对应的元素, 并在工程文件中写入该基组名和对应文件名即可.

例如, 如果我们想使用 Zn 的基组 TZVP-MOLOPT-SR-GTH, 那么将基组部分 Zn 所使用的基组名改为 TZVP-MOLOPT-SR-GTH 后, 再在 DFT 部分中增加一行:

BASIS_SET_FILE_NAME BASIS_MOLOPT_UCL

如果使用超算节点遇到找不到 BASIS_SET 文件的情况, 一般是 CP2K_DATA_DIR 环境变量设置存在问题, 这时候需要自己准备一套文件 (把源代码拷贝上去), 然后在

module load cp2k
# 或对于自己配置的cp2k
export PATH=$PATH:/path/to/cp2k/bin

下一行补上源代码 data 文件夹的位置:

export CP2K_DATA_DIR=/path/to/cp2k/data

由于MOFs 材料的晶胞一般较大, 为了保证精度与速度的平衡, 我们一般会选择 2-zeta 基组, 即 DZV/DZVP 开头的基组. 常见的包括 DZVP-GTH-PBE, DZVP-MOLOPT-SR-GTH和DZVP-MOLOPT-PBE-GTH 等, 这三种基组都是 2-zeta 基组, 精度依次提升, 计算量也逐渐增长. 如果想要高精度就选择最后一类. 当然, 如果想要更好的结果, 可以将重金属原子的基组提升到 3-zeta (即 TZVP 开头的基组).

需要注意的是, 配合 GTH 赝势使用的 GPW 基组不能与与全电子 GAPW 基组同时使用.

4.4 输入文件的调整

晶胞尺寸, 基组, 色散校正, 截断能等是 MOFs 材料需要首先关注的设置项.

对于 PBE 纯泛函优化, 超过 500 原子就比较慢了. CP2K 的 OT 方法相比 Diagonalization 更为迅速, 但仅支持 Gamma 点优化. 因此, 为了在速度和精确度中取得平衡, 当 MOFs 材料的各方向晶胞尺寸大于等于 20 Å (此时晶胞尺寸约 8000 Å^3) 时可以选择单 Gamma 点优化, 而某个方向小于 15 Å 时应当考虑选择扩胞后使用 OT 算法, 或使用 Diagonalization 方法并在该方向增加 K 点数量, 哪种方法能收敛 (而且够快) 就用哪种.

对于截断能, 一般需要进行截断能测试, 官网的资料写的已经很清楚 How to Converge the CUTOFF and REL_CUTOFF. 对于结构优化, 截断能设的太高意义不大, 速度变慢效果也没什么变化. 但是能量计算可能需要设定较高的截断能 (不知道选择 2000 Ry. 截断能的文章算了多久…).

色散校正选择 DFT-D3(BJ) 即可, 对于有机配体较多的 MOFs 材料结构优化以及弱相互作用处理都有一定增益效果.

4.5 截断能的设定

截断能需要进行收敛性测试, 但偷懒一般也就省下来了. 测试方法见 How to Converge the CUTOFF and REL_CUTOFF

5 能量计算与结构优化

5.1 能量计算

CP2K 的全局计算模式有多种设置, 其中最基本的是能量计算. CP2K 的 PDOS, 能带, TDDFT 等计算一般都在能量计算模式下进行. 设置方法为 RUN_TYPE energy.

&GLOBAL
  PROJECT ch2o_pbe0_rks_s_tddfpt_admm2
  RUN_TYPE energy
  PRINT_LEVEL low
&END GLOBAL

此外, 在 GLOBAL 中, 还可以设置程序输出的信息量. 许多内容需要在 PRINT_LEVEL MEDIUM 层级才能看到, 此时需要进行调整.

5.2 结构优化

CP2K 的结构优化包括 GEO_OPT 和 CELL_OPT 两种, 分别对应优化结构与优化晶胞与结构. 只要 SCF 能收敛, 一般也能进行正常的结构优化. MOFs 材料簇模型的结构优化经常要加入一定限制, 但在CP2K中做的周期性结构优化不加入限制性条件也能获得较好的结果. 一般来说 MOF 的单胞是可靠的, 所以可以用 GEO_OPT

CP2K 支持的限制条件包括原子的位置, 晶系, 空间群, 晶胞角度项, 一个或多个晶胞尺寸. 除此之外, 外加电场, 压力, 温度 (感觉好像对 DFT 结构优化不生效, 没具体测试过) 也能进行设定. 对于 CELL_OPT, 所有的设置项都在 CP2K_INPUT / MOTION / CELL_OPT 内, 可以在官网 Manual 查看具体设置方法.

如果需要优化压力条件下的结构, 可以选择输入一个分量, 也可以分别指定压力张量的九个数值. 在目前版本的 Multiwfn 可以直接在生成过程中指定. 然而必须要注意的是, CP2K 优化加压结构每次结果都可能不一样, 不是很可靠, 只能做极为有限的参考. 当然如果已经有了确定单胞, 还是值得试试看的.

6 PDOS 与能带计算

6.1 PDOS

PDOS 计算设置项在 FORCE_EVAL/DFT/PRINT/PDOS .

&PDOS
  NLUMO 50
&END PDOS

PDOS 的要点包括:

• 默认模式下能够输出不同元素的 PDOS, 加入 COMPONENTS 关键词能够输出不同轨道的分量.
• 为了输出空轨道的分布, 需要设置 NLUMO 关键词为所需空轨道数目. 设置为 -1 即输出全部能量, 但此时数据量会非常大.
• 当需要按照不同组分输出时, 可以设置 LDOS Subsection, 通过设置不同原子的列表进行 PDOS 的计算. 可以考虑使用 GaussView 获取原子列表, 它的原子选择功能较为方便.

PDOS 计算结束后可以使用 cp2k-output-tools 进行转换. 链接:

• cp2k-output-tools
• cp2k-output-tools PyPI

以下是批量转换脚本:

#!/bin/bash -l

# conda activate cp2k-parser # 使用 conda+pip 安装了cp2k-output-tools

icc=0

nfile=`ls ./*.pdos|wc -l`

for inf in *.pdos

do

((icc++))

echo Converting ${inf} ... \($icc of $nfile\)

cp2k_pdos --sigma=0.005 --total-sum ${inf} > ${inf//pdos/csv}

echo ${inf} has finished

echo

done

该软件包包括四个部分:

• cp2kparse … parse CP2K output (for restart & input files look at the cp2k-input-tools project) and allow easy selection of common values.
• xyz_restart_parser … when restarts occur during an MD you may end up with duplicated frames in the trajectory, this tool filters them (and can easily handle huge files)
• cp2k_bs2csv … convert a CP2K band structure file to multiple (one-per-set) CSV files for easier plotting. There is also an API available if you need to import bandstructure data into your application.
• cp2k_pdos … apply a convolution with Gaussians on a regular grid on the CP2K PDOS output and generate a CSV file for further processing or plotting. The same grid is used for all input files with the min/max of the grid automatically determined, but no summation of the different projections is done.

cp2k_pdos 的具体设置项, 较小的 Sigma 对应较窄的展宽.

$ cp2k_pdos -h
usage: cp2k_pdos [-h] [--sigma SIGMA] [--de DE] [--scale SCALE]
                 [--total-sum] [--no-header] [--output OUTPUT]
                 <PDOS-file> [<PDOS-file> ...]

Convert the discrete CP2K PDOS points to a smoothed curve using convoluted
gaussians. If you have separate alpha/beta-spin pdos files (spin-
unrestricted calculation), pass both files as arguments to get one common
grid for both of them.

positional arguments:
  <PDOS-file>           the PDOS file generated by CP2K, specify two files
                        (alpha/beta) in a spin-unrestricted case

optional arguments:
  -h, --help            show this help message and exit
  --sigma SIGMA, -s SIGMA
                        sigma for the gaussian distribution (default: 0.02)
  --de DE, -d DE        integration step size (default: 0.001)
  --scale SCALE, -c SCALE
                        scale the density by this factor (default: 1)
  --total-sum           calculate and print the total sum for each orbital
                        (default: no)
  --no-header           do not print a header (default: print header)
  --output OUTPUT, -o OUTPUT
                        write output to specified file (default: write to
                        standard output)
$ cp2k_pdos -s 0.01 -o output.csv k1-1.pdos

6.2 能带

CP2K 的能带计算非常迅速, 但目前 CP2K 仅支持纯泛函下的能带计算, 如果需要杂化泛函能带计算要选择 VASP 或 QE . 不过对于大的 MOFs 体系使用其他软件计算杂化泛函能带的计算量比较大, 计算的文章也较少.

注: 2024版本的CP2K已经支持 k 点计算了, 但是 MOF 体系还是算不动

设置项比较简单, 首先使用 vaspkit 或 seek-path 基于优化后的结构生成 CP2K 文件, 然后修改基组, 截断能即可.

(2024.10.04补充) 在logzzz – cp2k利用纯泛函进行能带计算的简单方法中, 有一段调用vaspkit和multiwfn用于自动生成能带计算的脚本, 可以参考一下. 具体使用方法见以上链接:

#!/bin/bash

###############################
###zhaicg,[email protected]
##############################

echo " "
echo "利用vaspkit,multiwfn生成band计算文件"
read -p "输入文件(需要加.inp):" flag
echo " "
read -p "输入维度1)1D 2)2D 3)3D:" table
echo " "

cp $flag $flag.bk
Multiwfn $flag << EOF >/dev/null
100
2
27
POSCAR
0
q
EOF

(echo 306; echo $table) |vaspkit >/dev/null

sed -n '/&PRINT/,/&END PRINT/p' KPATH.cp2k >cpkp
sed 's/<NAME_HERE>.bs/BAND.bs/g' cpkp -i
sed '/BAND.bs/a\\ADDED_MOS 20' cpkp -i
sed '4s/^/            /' cpkp -i
sed  '5d' -i cpkp

sed '$d' -i $flag
sed '$d' -i $flag
cat cpkp >> $flag
echo "  &END DFT" >>$flag
echo "&END FORCE_EVAL">>$flag
rm cpkp KPATH.cp2k  POSCAR  SYMMETRY

6.3 轨道等值面图

轨道等值面图一般都需要 virtual MOs 信息以绘制 LUCO, 所以必须使用 Diagonalization 方法进行计算, 并使用 ADDED_MOS 添加 virtual MOs. 随后使用 Multiwfn 输出 cub 文件, 处理方法见 详谈使用CP2K产生给Multiwfn用的molden格式的波函数文件 然后使用喜欢的绘图软件作图. TDDFT 计算的 out 文件配合任意一个 molden 文件也能绘制空穴-电子分析图.

6.4 数据后处理

CP2K 官方有一个 PDOS/能带后处理的 py 包. 也可以调用 multiwfn 中的 cp 功能进行能带处理. 推荐后者.

基于 cp2k-output-tools

cp2k_bs2csv 的输入较为简单, 一般直接按照以下输入运行即可:

$ cp2k_bs2csv output.bs

但是 cp2k_bs2csv 输出的数据与 vaspkit 不同, 坐标为 k 点位置, 因此直接画图是等宽的. 为了根据倒易空间中实际的 k-path 长度进行绘图, 我们需要知道每一条 k-path 的长度, 这可以基于倒易空间晶胞来计算.

较为简单的方法是基于 seek-path, seek-path 在生成路径时同时会输出倒易格子基矢信息与高对称点坐标, 基于这两者即可计算出不同路径的长度 (如下图)

![[Pasted image 20221128163951.png|350]]
因此, \(\Gamma – X\) 的距离为 0.157, 根据不同的路径手动设置 X 值位置即可.

基于 Multiwfn

首先用 Multiwfn 打开 .bs 文件, 然后输入 cp, 回车. 对于能带来说, 可以先用功能 2 转换能带文件, 然后使用子功能 -3 输出曲线数据, 然后画图.

7 杂化泛函计算

使用 Multiwfn 设定使用杂化泛函十分方便, 需要注意的是基组需要基于 data 文件进行一定调整. 常用的杂化泛函可以参考 Multiwfn 给出的选项, 也可以参考 VASP – List of hybrid functionals

参考资料:

• CP2K Tutorial – DFT and beyond – Hybrid
• CP2K:杂化泛函 (关于 hfx_eps_warning, ADMM 基组选择, Cutoff Radiis Warning)

CP2K 进行杂化泛函计算的优势在于其可以使用 Auxiliary Density Matrix Methods (ADMM) 降低内存使用. 杂化泛函工程的决速步与内存占比最大的计算在于 four center electron replusion integrals(ERI), 如果体系较大内存不足, 就会面临 ERI on-the-fly 的情况, 严重降低计算效率. 因此, ADMM 能够有效提升杂化泛函的计算速度 (大约只比 GGA 高2个数量级).

目前测试得到的计算速率, 对于 400 原子 5000 基函数的体系, 截断能设置为 600 Ry – 55 Ry, 基于收敛的 PBE 波函数作为初猜, 平台配置双路 8375c, 2 进程-32线程并行的计算, PBE0-ADMM 第一步需要 4 分钟, 总计大概 10 分钟能够进行一次能量计算. 内存消耗量为 150 GB. 事实上这样的体系也应该能够进行基于杂化泛函的结构优化, 时间预估在 2-3 天左右.

当进行杂化泛函-ADMM 计算时, 建议使用 MOLOPT 基组. CP2K 的杂化泛函计算支持 OT 和对角化, 但不支持 K 点, 因此不能算能带.

8 TDDFT 计算

CP2K 的 TDDFT 能够计算单线态, 三线态 (以及开壳层的 N 线态)的能量, 跃迁偶极矩与 NTO 分析, 现在也支持激发态结构优化.

但是, TDDFT 计算不支持 K 点, 因此想要得到不受质疑的结果, 需要对晶体结构进行扩胞. 然而扩胞会导致激发态的数目加倍, 原本只有 1 个的激发态会因为一个方向加倍变为 2 个, 导致为了捕捉更高能量的激发态也需要将所计算的激发态数目加倍. 好在 TDDFT 对于内存的消耗量只会因为激发态数量的增加而线性提升, 计算耗时也和激发态数目线性相关. 对于基于 GGA 的激发态计算, CP2K 的内存消耗量比较小, 但是对于杂化泛函的激发态计算则较大. 因此, 5000 基函数的体系如果需要计算 200 个激发态, 每个进程需要至少 256 GB 的内存, 少于 ORCA 计算所需的内存量级, 但远多于 Gaussian 所需的最低内存.

TDDFT 的计算中有许多限制, 可以参见源代码 qs_tddfpt2_methods.F/SUBROUTINE tddfpt_input 部分中的警告. 目前已经发现的包括不支持 K 点, 杂化泛函-ADMM 限制等. 对于杂化泛函-ADMM, 要求以下设置:

&AUXILIARY_DENSITY_MATRIX_METHOD
      METHOD BASIS_PROJECTION
      ADMM_PURIFICATION_METHOD NONE
    &END AUXILIARY_DENSITY_MATRIX_METHOD

对于杂化泛函的 TDDFPT 也要求必须输入 XC Section:

   &XC                            ! If choosing kernel FULL, the underlying functional can be
     &XC_FUNCTIONAL PBE0             ! specified by adding an XC section
     &END XC_FUNCTIONAL              ! The functional can be chosen independently from the chosen
   &END XC                        ! GS functional except when choosing ADMM

成功的TDDFPT-PBE0示例 (不包含数据文件):

2024.10.04 补充: 这里的XC部分少了交换能的25% HF 成分, 但当时的版本(2022)如果写上这部分在计算时会有无法继续的错误. 在后续版本重新尝试时该错误已不存在, 且不再必须在此处写 XC 部分了. 但如果运行过程中报错, 可以考虑检查一下此部分的内容. 现在推荐使用 Multiwfn 直接生成 TDDFT 文件, 最为简便.

  &PROPERTIES
    &TDDFPT #TDDFT calculation
      NSTATES  400
      MAX_ITER 100
      RKS_TRIPLETS .F. #If calculating triplet rather than singlet excited states
      CONVERGENCE [eV] 1E-7 #Convergence criterion of all excitation energies
      MIN_AMPLITUDE 0.01 #The smallest excitation amplitude to print
      RESTART .T. #If restarting TDDFT calculation. If true, WFN_RESTART_FILE_NAME should be set to previous .tdwfn file
      WFN_RESTART_FILE_NAME 1-360-RESTART.tdwfn
      &XC
        &XC_FUNCTIONAL
          &PBE
            SCALE_X 0.75 #75% GGA exchange
            SCALE_C 1.0 #100% GGA correlation
          &END PBE
        &END XC_FUNCTIONAL
        &XC_GRID
          XC_DERIV SPLINE2_SMOOTH #The method used to compute the derivatives
        &END XC_GRID
      &END XC
      &PRINT
        &MOS_MOLDEN
          NDIGITS 8
        &END MOS_MOLDEN
      &END PRINT
    &END TDDFPT
  &END PROPERTIES

使用 CP2K 计算MOF发光性质的限制

CP2K 有把所有的发光过程都当做 CT 态的倾向, 对于分析机制来说不够可靠.
CP2K 能够计算 TDDFT 过程中的旋轨耦合常数, 但对于 MOF 来说很不靠谱 (有磷光性质的材料 SOC 却是0)

9 吸附结合能计算

结合能计算比较简单, 但能不能得到想要的结果有点看运气 (不知为何之前算过的结合能不是极小就是负值, 很奇怪). 首先进行复合物结构优化, 然后在此基础上分别计算二者的能量. 再用复合物能量减去单体能量之和即可.

10 红外光谱计算

红外光谱计算与振动模式模拟参考以下文章: 使CP2K计算的振动模式可以被GaussView观看的程序：MfakeG

补充

2024版本后支持带k点的杂化泛函计算, 但是需要注意的是能带计算时, 必须设置 DFT-KPOINTS 选项, 否则程序会卡住.