第一原理计算简介和VASP入门

       借助边值问题的自洽求法和本征值问题的变分方法，可以有效处理单电子体系的定态薛定谔方程。在实际的材料计算中，我们需要采用第一原理方法处理多电子体系。以VASP软件为例，我们尝试用简单的语言展示密度泛函理论（DFT）的主要思路，希望对入门的同学有所帮助。首先声明相关内容图片来自网络，并参考了VASP官方网站的原理介绍部分：https://www.vasp.at/wiki/index.php/Lectures_and_presentations

       本期内容包含四个部分：（1）第一原理方法简介；（2） 密度泛函理论概述；（3） 赝势方法和展开函数；（4） VASP快速入门。

       首先，物理学中的第一原理方法是把已建立的物理定律作为出发点，而不作经验模型和拟合参数等假设。在材料模拟中，第一原理计算是不包括模型近似与实验数据拟合，通过求解薛定谔方程计算系统的电子结构。如下图所示，我们可以通过实验测量弹簧长度和回复力大小，拟合得到劲度系数；而在第一原理角度看来，不同应变对应的是原子位置的变化，我们求解薛定谔方程得到能量变化，就能得到体系的弹性系数。

       在第一原理计算中，描述系统的是量子力学基本原理，用到的基本参数是质子和电子的质量、电荷以及相关物理常数等。对于单电子的定态薛定谔方程，之前的自洽和变分方法都可以容易解决，而实际材料中，我们面对的是多电子系统，体系哈密顿量通常为：

       这里包含了离子、电子动能、离子-电子、电子-电子、离子-离子相互作用，非常复杂。狄拉克曾经说过，我们找到了绝大部分物理学和全部化学领域所需要的数学工具，但无法用含有几个摩尔数量级的偏微分方程组得到精确解。考虑到原子核质量远大于电子质量，我们可以认为原子核提供势场，其“位置”相当于输入参数，可以把原子核和电子的运动分离，这就是“绝热近似”。此外，忽略电子之间的关联和多体效应，在平均场的框架下，把基态波函数近似为单电子波函数之积，通过变分原理，我们可以得到Hartree的自洽场方程：

       多电子问题转化为单电子问题，数值求解成为可能。考虑到电子波函数的交换反对称，我们可以通过Slater行列式构造Fock态即可（相关内容参考量子力学教材）。在Hartree自洽场的基础上，看密度泛函理论就“亲切”不少。这里最主要的是Hohenberg-Kohn定理：多电子系统基态能量是电子数密度函数 ρ(r) 唯一泛函，在电子数不变时，正确的ρ(r)使能量泛函 E[ρ(r)] 达到极小。根据波函数和电荷密度的关系和变分，得到单电子方程（Kohn-Sham方程）：

       其中V_KS包含了外场V_ext（一般假定原子核不动，这里指电子感受到原子核的静电作用）、库伦作用V_coul（电子云之间的静电作用）、交换关联作用V_XC的贡献（这里波函数没有写出Fock态，交换反对称没有体现，和电子关联作用合并到V_XC中）。

        交换关联泛函V_XC的形式是关键的，无法通过基本规律推导得到。最初采用局域密度近似(LDA)，认为V_XC和局域电荷密度相关，通常会高估体系的结合能、低估体系的带隙；广义梯度近似(GGA)考虑电子密度梯度的修正，可以得到合理的结合能，GGA泛函可以运用拟合方法提高泛函精度，如BLYP，也可以通过物理限制得到较好的泛函形式，如PBE。但GGA计算的带隙仍偏小，可以使用杂化泛函HSE03,HSE06，得到与实验符合的带隙数值。对于具有局域的d、f 轨道，需要DFT+U来处理，此时参数U可以通过实验晶格常数等来确定合理范围。

       在实际的第一原理计算中，我们选择合理的展开函数，对k空间进行划分，考虑电子占据数等，这样才能准确有效的计算材料电子结构。在晶体对应的周期势场下，体系的波函数满足布洛赫定理，可以用平面波展开，而对应截断能是影响精度和速度的关键。

       体系的物理量（比如电荷密度等）计算是在k空间积分，需要对k空间进行划分得到离散化的k点和权重。而电子占据满足费米分布，在零温时对应阶跃函数，在结构优化时影响效率。

       一般来说，最外层电子对化学键贡献是主要的，芯电子通常不受影响。赝势方法可以把电子波函数分为外层电子和芯电子两部分，并使其正交，得到比全电子波函数平缓的赝波函数。赝波函数和全电子波函数对应的本征能级相同，在截断半径外，两类波函数重合。模守恒赝势要求，两类波函数在截断半径内对应电荷积分相同。超软赝势取消该约束，得到更平滑的赝波函数。PAW方法通过相对简单的操作，将赝波函数变换得到全电子波函数，逐渐成为主流。

       最后，我们看看VASP的输入文件（POSCAR、KPOINTS、POTCAR、INCAR）和以上理论基础的对应。POSCAR文件给出了体系的晶格矢量和原子位置，如果是团簇，可以选择一个比较大的晶格，使周期边界条件后的“镜像”分子相互作用可以忽略即可。下面是一个二元体系的POSCAR案例，中间的3x3矩阵是晶格矢量，Si C是对应的元素，下面的数字是原子个数，Direct对应晶格内坐标，再往下就是具体的坐标分量。POSCAR确定后，可以选择自动划分的MP方法，这里是4x4x4的k点划分。POSCAR是面心立方，k空间的倒格子是体心立方。

       根据POSCAR的元素，我们需要找到相应的POTCAR并组合。在PAW-PBE的POTCAR文件夹下，有各个元素对应的多种POTCAR文件。比如Na和Na_pv，前者只包含了一个电子，后者包含了7个电子，  VRHFIN =Na: s1p0，  VRHFIN =Na: p6s1，对应的截断能ENMAX分别是102 eV和260 eV。具体选择哪个文件更合理，速度更快，需要提前测试。INCAR文件包含了平面波展开的截断能设置、电子占据的设置、结构优化和电子自洽的收敛精度等。按要求准备好4个文件后，运行VASP，程序会根据POSCAR的原子位置，完成电子自洽，比如这里DAV从1到10，然后调整原子位置，再做电子自洽，最后达到用户设置结构优化收敛标准。

        在VASP的输入文件INCAR中，程序有很多默认的参数值，有需要的朋友可以查阅相关说明。总的来说，最关键的是POSCAR文件，确定初始结构合理后，剩下三个文件可以通过脚本自动配齐，比如INCAR里面选择一些通用设置（不是最优参数，但定性没有问题），可以实现简单的自动化计算流程。

        在做研究之前，我们应该系统测试POTCAR是否合理，如果计算氢气和氧气都有问题，那么接下来计算的水分子可能问题更大。第一原理计算有很好的可重复性，和实验中的晶格常数和内聚能都是符合的。参数测试需要符合数值方法原理，在后面的内容中再详细分析。