初步了解第一性原理计算
我初步了解第一性原理计算是因为研一的时候，想要分析Pd表面吸附苯和苯酚的吸附能并作可视化的弱相互作用分析，当时的我只了解了ORCA，所以做的是Pd团簇模型，按照sob老师做Ag(111)吸附苯的案例尝试了一下，更换了很多泛函和基组，也尝试了很多解决收敛的方法，最终都以无法收敛告终。这让我认识到量子化学程序可能只适合部分过渡金属团簇模型的计算，真涉及周期性的计算模型，还是得使用第一性原理计算程序，但初出茅庐，第一性原理计算涉及的输入设置更为广泛，记得当时光是查阅相关名词和读懂概念就花了得2个多月，而真正让我将大部分碎片信息组装在一块的点是在我阅读了《DENSITY FUNCTIONAL THEORY：A Practical Introduction》，这本专著不像Levine主编的第七版《Quantum Chemistry》那么深入，机理讲得不多，对于我们了解专有名词和初步概念足矣，并且书中还阐述了孤立体系和周期性体系计算的差别及共性，前者主要用在量子化学计算中，采用空间局域化函数，计算量稍小，所以有很多高级的方法，而后者用于第一性原理计算中，采用空间周期性函数（又叫平面波函数），这种函数因为不存在边界，所以通常需要取截断来近似获取精确值，计算量大，不适合用于大体系，目前应用最多的就是广义梯度近似（GGA）泛函里的PBE泛函，精度不算高，但限于目前的计算机的算力，该泛函应用得最广。同时，书中也阐述了何为布里渊区（Brillouin Zone），也即第一性原理计算程序中需要设置的k点，这里只需要简单知道周期性模型三个矢量方向上的k点数目比例与模型里三个矢量方向的矢量长度成反比例关系即可，但一般以$15$ $\mathrm{\mathring{A}}$为最小的基数，即k点数量与对应矢量方向的矢量长度的乘积至少需要大于$15$ $\mathrm{\mathring{A}}$。最为常见的Slab模型在书中也有具体的讲解，并对k点数量的提升对体系能量变化的影响，对称性对k点数目的可约影响以及Slab为什么至少需要3层原子（涉及表面结构弛豫），做了带案例的描述，可以自行仔细阅读，了解更多知识。最后阅读完这本专著，就会明白，做第一性原理计算时，需要提供的参数至少包括平面波函数截断能、电子密度函数截断能、k点设置、模型大小、赝势类型等等，这里也提示我们在看别人论文的时候能有基本的判断，仔细观察别人计算细节部分有没有这些参数，不全的可以直接理解为作者藏着掖着，对其结果只能一看而不能全信。当然如果是我们自己在写论文时，自然也要具体的给出以上提到但不限于以上参数的描述，这对于其他读者是非常有用的计算信息，缺失任何一项就相当于我们做实验在实验步骤漏写一个关键步骤一样，其他人没法很好的重复这个结果，这是对自己和对ta人的不负责！！！
As we stated above, it is not appropriate to view methods based on periodic functions as ?right? and methods based on spatially localized functions as ?wrong? (or vice versa).【两类计算方式都没有对错，只有各自擅长和不擅长的地方】
《DENSITY FUNCTIONAL THEORY：A Practical Introduction》 p.27
《DENSITY FUNCTIONAL THEORY：A Practical Introduction》网上有影印版的，只需要在搜索引擎的搜索框搜索：
“DENSITY FUNCTIONAL THEORY” “A Practical Introduction” filetype:pdf
请务必使用谷歌或者必应搜索，百度完全搜不到，非常的拉跨，也不推荐使用百度搜索
目前常用的第一性原理计算程序：CP2K与Quantum Espresso
CP2K与Quantum Espresso（后边简称QE）都是免费开源的第一性原理计算程序，可以自行到其官网或者Github上获取程序。
CP2K: https://github.com/cp2k/cp2k
QE: https://github.com/QEF/q-e
先说说QE，QE和商业软件VASP和学术免费的CASTEP一样属于纯纯的平面波计算程序，优点就是精度高一些，缺点就是算得慢，对于56核的机子而言，大一点的体系100来个原子左右的就非常的吃力。CP2K采用的是高斯平面波混合的方式进行计算的，速度上往往会比QE快不少，而且能算很大的体系，250个原子左右的体系也没啥大问题，所以特别适合拿来跑动力学以及结构优化，其缺点呢就是平均精度不及QE这类平面波计算程序，但往往还可以接受。所以说这两个软件优势互补，利用好了可以很大程度的扩展自己想研究的体系，并节省计算时间。
输入文件
CP2K与QE的输入文件在形式上很类似，都是类似字典的结构还有Json文件结构，即键值对嵌套叠加。
1 | { |
CP2K的输入文件
接着来看CP2K的输入文件（这里用的是sob老师的Multiwfn生成的），这里用的是之前我建的一个CU(110)面吸附甲苯的结构：
1 | #Generated by Multiwfn |
QE的输入文件
下边是QE的输入文件，我稍微做了点修改（这里用叶老师的QEtoolkit生成的）：
1 | ! generated by QEtk 2022-04-29 20:58:25 UTC+8 |
可以直观感受一下输入文件对应的各个块体，我们修改相应的值和参数先找到大类再对应去修改小类即可。
CP2K的Manual: https://manual.cp2k.org/
QE的PW.x模块的Doc: https://www.quantum-espresso.org/Doc/INPUT_PW.html
收敛性测试
前面也提到了，第一性原理计算程序要获取准确的平面波函数是不可能的，一般需要进行截断以获取近似的值。截断能的选取对体系而言就非常重要了，一般来说对于首次计算的体系，必须进行对应的收敛性测试（省时间一般采用计算单点能的任务测试即可），当然也可以参考别的研究者的体系来选取初步的截断能设置，这是比较省事的做法，也很推荐大家使用这个方式。对于QE而言，可以参考Standard solid-state pseudopotentials (SSSP)的测试，一般来说很合理，对于存在多个元素的体系，一般可以选取最大值即可，但仍然建议做一下收敛限测试。
SSSP网站是类似一个测试集，评价了很多种赝势，并从中挑选了使用起来较为稳定的赝势，一般我们可以直接使用
说到赝势，这里必须提到，从原理上说波函数截断能（ecutoff）与电荷密度截断能（ecutrho）的关系是：
$$ecutrho = 4ecutoff$$
经过大家的测试发现，很多时候ecutrho的值往往要比4倍的ecutoff还要大，甚至是8倍10倍，尤其是对超软势（UPF）而言，所以大家在设置截断能时，若没有充分的参考，非常有必要进行收敛限测试，确保后续计算的偏差值不大才不会白白浪费机时，无论是使用个人服务器还是超算，这是对自己研究的负责！
还需要进行测试的有展宽设置、K点数量、磁序、偶极矩等，做对应的第一性原理计算需要考虑的东西比孤立体系的量子化学计算需要的更多。
Slab建模要求
- 建Slab模型一般而言原子层数至少需要3层，低于3层的对表面原子的结构弛豫影响非常大！真空层厚度至少$15$ $\mathrm{\mathring{A}}$如若是建立的吸附模型，还应该考虑吸附分子的结构，比如是长条的，那么Slab也应该是长条的，如果是圆的或方的，那么Slab也应该是方的。
- 吸附物质至少离吸附表面有$3$ $\mathrm{\mathring{A}}$的距离，这样让吸附物质通过弛豫和吸附作用使体系落到能量极小点内是非常关键的，如果距离太近，吸附物质与吸附表面间的范德华作用会大大增强，可能会偏离预期，也可能导致体系难以收敛。通常“成键”的距离在$1.6$ $\mathrm{\mathring{A}}$左右。
- 吸附表面至少需要大吸附物质一到两圈，以减少吸附物质间的相互作用影响。
k点的设置的说明
通常晶胞的3个方向的矢量长度至少是$15$ $\mathrm{\mathring{A}}$，如若不足，比如$|{\overrightarrow a}|=3$ $\mathrm{\mathring{A}}$，那么${\overrightarrow a}$方向的k点数量至少是5。
三个方向的k点数量与矢量的长度成反比关系。比如：
$$|{\overrightarrow a}|:|{\overrightarrow b}|:|{\overrightarrow c}|=5:5:10$$
那么对应的方向上的k点数量就应该是这样的倍数关系，且k点数量应该为整数：
$$k_{a}:k_{b}:k_{c}=n(3:3:1.5) ,n\in Z^+$$这里的a，b，c三个方向上的k点数量可以是3：3：2或者6：6：3。
建立的Slab模型在真空层方向一般不加k点数量，即设置为1，因为真空层一般足够大，而大多Slab模型如果建立的够大，一般使用$\Gamma$（gamma）点即可，即中心一个点。
展宽设置
展宽是在费米能级附近做展宽，所以在之前做收敛性测试时，一般展宽对体系的能量贡献越小越好，对QE而言即smearing contrib.的值应小于$1 {\times}10^{-4}$ $eV/atom$，金属或者导体的degauss值一般在$0.01$ $RY$左右，可以适当增大来使体系收敛，半导体设置尽可能小（$1 {\times}10^{-9}$ $RY$），而绝缘体应设置为0。