3.3 结构优化

什么是结构优化（relaxation）？

一句话概括就是：寻找当前结构对应势能面所在位置附近的局域最小值。结构优化的目的是为了在给定限制（压强、是否固定晶格等条件）下得到合适的晶格常数和原子位置。

什么是势能面？

当一个体系的元素以及原子个数确定时，想象一个由原子坐标（对于复杂体系还有磁矩等参数）和体系能量构成的高维空间，体系能量完全取决于原子坐标。所以不同的原子坐标（以及磁矩等）就会对应不一样的能量，如果遍历全部的坐标可能，就得到了能量-坐标的一个高维曲面，即势能面（PES，Potential Energy Surface），下图：

如果原子的坐标不够好，使得该结构落在势能面的坡上（而非谷底），那么该结构必然会有自发向低能谷底的趋势（通常情况下自然总是喜欢能量低的，因为稳定）。对应在物理上就是该结构的原子间的受力没有平衡掉，原子为了让体系趋向低能的势能面谷底，就会在力的作用下略微移动；或者该结构的晶格过大（过小），体系在当前压强下要么收缩要么膨胀，直到它落在势能面谷底（局域最小值）。

原则上，我们可以对一个落在势能面坡上的结构手动修改每一个原子的位置或是手动对晶格放缩，不断重复直到它落到势能面的局域最小值（此时原子间受力为零或者微弱不计，体系能量达到极小值）。但是没有必要，VASP有成熟的结构优化算法，由参数IBRION控制（请自行查阅）：

https://www.vasp.at/wiki/index.php/IBRION

结构优化算法流程：

这其实是QE的结构优化流程，懒得重新画，就直接拿来用了.VASP的区别就在于右侧的算法用的不是BFGS，而是前面IBRION提到的那些。实际上就是“离子步”和scf嵌套在一起的for循环：每移动一次离子，做一次scf计算得到每个原子的受力和体系的能量（这一过程称为离子步）。再借助算法得到下一个离子位置，进入下一次scf计算，直到达到收敛阈值。

3.3.1 结构优化参数介绍及实例

（1）首先你需要一个POSCAR，内容如下：

注意，这里的晶格常数是人为随意取的一个较为合理的值（为了演示用）。

（2）结构优化采用的INCAR内容如下：

相比于前面的静态自洽计算，结构优化多了IBRION，ISIF，NSW，POTIM，EDIFFG，PSTRESS等参数。

首先需要强调，这里对于Si来说，不需要考虑磁性、强关联（+U）、SOC、VDW等效应，也不采用高阶泛函（如PBEsol、SCAN等），因此INCAR中暂时不写出这些参数。其次我们对这些参数做一些说明：

PREC=Accurate，涉及到力的计算时（结构优化/声子频率等）取Accurate，其他计算可以注释掉。在这里我们注释掉了以减少计算量，可以自行比较注释前后的结果差异。

控制结构优化的主要参数：

NSW，最大离子步，当300个离子步还没有达到收敛标准就自动停止。

POTIM，控制结构优化过程中的原子移动距离，单位Å，默认结构优化时取0.5

IBRION，控制结构优化算法，在结构优化时可以取1 2 3三种，一般取2

ISIF，控制结构优化过程是否变晶格/胞。常用的是ISIF=2和3；2对应固定晶格优化原子位置（晶格长度不变），3是同时优化晶格与原子位置（晶格长度也会变化）。

EDIFFG= -0.01，单位是eV/Å，当每个原子在xyz方向的受力都小于该值时，优化停止。大部分情况下-0.01已经足够（对于一些金属表面吸附分子的问题，由于原子数很多，可以适当放宽精度，EDIFFG用-0.05）；在算声子谱前的结构优化建议这个值设置为-0.001。

（EDIFFG可以取正值，那个对应的是前后两次优化的能量，这种情况用得不多）

请自行查阅官方手册的参数介绍：

https://www.vasp.at/wiki/index.php/Category:INCAR_tag

最后两行多了NCORE和KPAR参数，这两个参数是加速计算的参数，还有NPAR、KPAR参数。一般设置这两个参数比默认情况下，计算速度会更快。

https://mp.weixin.qq.com/s/uiWuvILE1WeiB5SC3_Frvg

（3）KPOINTS采用vaspkit-102生成，密度取0.025

（4）POTCAR用vaspki-103生成。

然后确保这四个文件以及任务提交脚本在同一个路径下，通过脚本提交计算任务。

当任务顺利结束时，你应该要在result文件一行最后看到如下内容：

或者在OUTCAR最后看到如下信息：

接下来我们讲解一下结构优化计算过程

首先是result文件，会看到左边栏目有一系列的1F，2F等；每一个F代表一次循环，这里叫离子步，每一次离子步中的电子步都收敛了。这里一共有4F，即进行了四个离子步最终结构优化完成。

现在我们查看relax过程中每一次离子步的体系的压强。

用如下命令grep ' in kB ' OUTCAR 屏幕会输出OUTCAR中的in kB关键词（自行在屏幕输入命令内容，直接复制可能有格式问题）

这里输出了四行，对应四个离子步的压强，从左到右依次是Pxx，Pyy，Pzz，Pxy，Pyz，Pzx，可以自行翻看OUTCAR文件找到这些内容相应的位置。第一个离子步输出的体系压强是418kB，即42GPa（10kBar=1GPa），意味着当前晶格常数太小，体系会膨胀（有一个正压强）。最后压强变成了0（预设PSTRESS=0kBar），晶格常数到了5.47Å。

同时我们还关心每个离子步中每个原子的受力情况，

自行翻看OUTCAR，找到这部分内容：

右边八行就是八个原子的受力，分别是xyz三个方向。这里每个原子受力都是0，是因为原子位置非常规整，没有偏离理想格点。也可以使用awk命令抓取力的部分：awk '/POSITION/,/total drift:/' OUTCAR （自行在屏幕输入命令内容，直接复制可能有格式问题）

结构优化的目的是为了得到合适的晶格常数和原子位置。该信息在CONTCAR文件中。我们查看CONTCAR文件得到relax后的结构的晶格常数，可以看到晶格常数变成了5.47Å。CONTCAR文件格式和POSCAR一样，因此可以直接把relax成功后输出的CONTCAR文件复制为POSCAR文件（新建文件夹操作，不要覆盖原来的relax中的POSCAR），用于后续其它计算。

XDATCAR文件里面有每一次离子步的结构信息，可以自行打开查看relax过程中的晶格常数变化。OVITO软件可以直接查看原子位置和轨迹，请自行网上下载（https://www.ovito.org/），并自学或向同门请教如何使用基本功能。

3.3.2 结构优化结果依赖于初始构型

时刻记住一句话，计算结果很大程度上取决于计算的初始条件。我们一开始给定的POSCAR中，八个Si原子位置都是规整的金刚石格点位置，如果有一些原子偏离了的话，会怎么样？relax能帮我们优化回完整格点吗？

现在我们构造一个这样的POSCAR

注意，我们修改了红框里原子坐标，使得这两个原子偏离了规整的金刚石格点。

其它INCAR，KPOINTS和POTCAR不变，提交任务，做relax计算。

首先我们会发现此时计算速度明显变慢。grep LOOP OUTCAR可以查看每一次电子步的时间，LOOP+是离子步的时间。这是因为结构对称性下降了，使得要计算的不可约k点也变多了。可以比较此时的IBZKPT文件与前面计算的IBZKPT文件中不可约k点数。

其次，当relax成功结束时（你需要查看result或者OUTCAR中相应的信息，确保relax成功结束），我们会看到CONTCAR似乎不是“金刚石”结构，好像是离金刚石结构差一点点，如下

运行vaspkit-608命令查看CONTCAR的对称性，会看到此时显示的是P1空间群，为什么不是金刚石的Fd-3m呢？因为此时我们给定的初始结构就不是Fd-3m了（有两个原子偏离了一点点），可以用vaspkit-601查看POSCAR的对称性。但可以想象，这个结构在构型空间中距离完整的金刚石结构应该“很近”，此时我们需要在寻找对称性时取一些宽松的标准。只需要vim ~/.vaspkit，修改文件中的对称性参数，如下，默认1E-5，我们改成1E-2。

然后保存并退出，再运行一次vaspkit-608找CONTCAR的对称性，就可以看到此时空间群是Fd-3m了。注意，608命令只是找对称性，并不会生成施加对称性后的新结构。这需要用602或者603功能实现：新建一个文件夹，把刚刚的CONTCAR复制进去并命名为POSCAR，再进入这个文件夹，运行vaspkit-603寻找当前结构的晶胞；此时生成一个CONVCELL.vasp文件，就是vaspkit读取CONTCAR的信息，并施加宽松（刚刚我们设置的1E-2）的对称性标准之后，重新找到的具有对称性的晶胞。运行vaspkit-602可以得到施加对称性之后获得的标准原胞。用完该功能后，记得把~/.vaspkit里面的参数改回去。

有一些体系的结构优化计算，如某些结构在不同压力下的焓，需要知道优化后的CONTCAR对称性是否还和POSCAR保持一致，可以通过vaspkit的608功能快速查看CONTCAR的对称性。VASP在计算时会默认开启对称性，这个对称性并不是使得结构在优化过程中空间群保持不变，由ISYM参数控制，具体含义见官网。一般的结构优化计算，ISYM取默认即可（注释）；有些特殊的计算可能需要设置ISYM=0或者-1，完全关闭对称性（但是计算速度也明显变慢），这属于进阶内容，不在本教程范畴。

3.3.3 限制下的结构优化

前面我们涉及到的结构优化主要是“所有原子坐标和三条轴全放开”或者“原子移动，轴固定”的结构优化。有些特殊的结构优化，可能需要固定住某些原子，或者固定住某些轴。这时候就需要在一定限制条件下做relax计算。

（1）固定POSCAR中某些原子的结构优化？

有些时候，我们需要研究一个气体分子在某个材料表面的扩散、吸附或表面重构问题时，此时的POSCAR结构原子很多，事实上只有表面附近的原子会移动（前提是你已经单独对材料表面、衬底优化过了）。因此为了减小计算量我们可以对远离表面的原子固定，只优化表面附近几层的原子位置。此时POSCAR写法如下：

注意看红框重点，首先多了一行Selective Dynamics内容；其次，每一个原子坐标的右边多了F和T，F表示固定。第一个F对应该原子的坐标的第一个值。这里这个POSCAR表示：在relax过程中，固定住第一个(0.25,0.25,0.25)原子，剩下七个原子的坐标在relax过程中可以移动。可用手动设置T和F，或者借助vaspkit-402功能实现。

如果固定的原子不合理，可能会导致有些原子受力始终没法低于EDIFFG，此时需要适当放开更多的原子继续优化。

（2）固定c轴的结构优化？（二维材料）

当我们研究二维材料时，通常在第三个维度Z方向加上15-20Å的真空层，来避免Z方向上周期性所导致的层与层之间的相互影响（这带来的后果就是在相同ENCUT下，有了真空层就需要更多的平面波数量）。在优化过程中我们只需要对该胞做固定c轴（放开ab轴）的结构优化，VASP目前本身不支持这一点，需要人为另外编译。请自行搜索“VASP固定c轴的结构优化”。

3.3.4 结构优化计算总结

（1）结构优化本质就是一系列依次进行的静态自洽计算，所以它相比于静态自洽计算（scf），就是多了一些控制离子步的参数，主要是这五个NSW，IBRION，POTIM，EIDFFG，PSTRESS。

（2）首先要确保初始结构合理，比如不能出现某些原子距离很近的情况，可以用原子的经验半径去判断，这个数据可查https://ptable.com/，一般两个原子距离不会小于0.8Å（除了一些特殊的分子，比如H2的键长0.74Å）；

（3）一定要查看result或者OUTCAR确保计算顺利完成。一方面，要确保每一个离子步中电子步都收敛，此时结构的演化才是正确的。如果偶尔几步离子步中出现电子步只迭代到了1E-3，始终没法达到预设的EDIFF；而后续的离子步中电子步又收敛了，此时的结果也是可信的。

同时，要确保最后压强趋于预设压强（1kbar量级的波动可以忽略），以及原子受力逐渐接近预设值EDIFFG。如果体系有磁性，那么结构优化过程中也需要开启磁性参数ISPIN，并且考虑设置初始磁矩；磁性部分的内容请看后面章节。

（4）如果你的relax计算顺利完成，可以读取最后一个离子步的能量，OSZICAR中的E0值，作为该结构优化后的总能。或者拿优化后的结构单独做一次静态计算以获得能量。可以自己测试对比二者结果的差异。

（5）有时候体系比较复杂时，比如大几十个原子甚至上百原子，结构优化的精度标准可以适当放宽，比如EDIFFG=-0.05；EDIFF=1E-5。或者采用两次优化的策略，第一次用低精度，k111+低截断能，目的是快速找到合理的局域能量最低构型附近的结构。第二次提高精度，读取第一次的CONTCAR进行relax计算。并且，在relax过程中，一定要查看压强和力，确保压强接近预设值（默认是0kbar，最终压强可以在1kbar量级附近抖动），以及力要小于EDIFFG；如果优化到最大离子步时，力仍然没有完全达到预设值，那需要把CONTCAR命名为POSCAR，其余INCAR，KPOINTS和POTCAR不变，新建文件夹，在里面继续优化。新建文件夹的目的是避免覆盖前一次relax计算，以便于后面查看数据。个人建议，养成新建文件夹计算的习惯，尽量不要覆盖原始数据，除非你已经很清楚地知道前面的数据对自己没有用。

（6）结构优化最常出现的报错就是当前EDIFF不够，此时result或者OUTCAR会提示你减小EDIFF或者复制CONTCAR为POSCAR继续优化，遵循程序的建议操作即可。

（7）结构优化要合理利用好ISIF参数， ISIF的取值有很多：

https://www.vasp.at/wiki/index.php/ISIF

通常情况下都需要ISIF=3的优化，相当于在某一压强下充分弛豫释放结构内部全部应力。但如果是研究超胞掺杂模型，或者是某一金属表面吸附官能团的问题，可以考虑做ISIF=2的固定晶胞的离子优化。ISIF=2和3是初学阶段用得较多的选择，其它设置可以根据需要采用。

3.3.5 结构优化计算不收敛怎么办

对于原子多的大体系，常出现的问题是结构优化不收敛。此时要分清楚两种情况：

（1）某一离子步中的静态计算不收敛，使得原子受力和压强不正常，最终结构演化到奇葩的情况。该情况，请参照“自洽计算不收敛”问题，先解决自洽计算不收敛。

（2）（静态计算收敛了）原子受力始终没法降低到EDIFFG以下。有时候查看relax过程中的原子受力，会发现达到最大离子步NSW时，原子受力仍然在0.1附近震荡，没法降低到0.01（假如EDIFFG=-0.01）。此时可能是INCAR中POTIM太大（默认0.5），体系一直在局域极小值附近震荡没法达到极小值，此时需要减小POTIM，比如改成0.2或者0.1，再复制CONTCAR为POSCAR，继续优化。或者提高EDIFF的精度，从1E-6提高到1E-8；或者修改IBRION算法，1和2都尝试，个人经验，2比1更好。或者网上搜索“结构优化不收敛”的相关经验。

3.3.6 补充内容

（1）另外一种寻找晶格常数的结构优化方法——做状态方程拟合

思路就是不断改变晶胞体积，做ISIF=2的计算，得到不同体积下的能量，拟合出一条关于体系能量-体积的曲线即E-V曲线，该曲线在常压体积附近是类似二次曲线，原则上要用Birch-Murnaghan方程拟合。通过拟合+解析函数求导，寻找该曲线最小值，即为平衡体积。一般来说，这种方法得到的晶格常数和ISIF=3直接优化得到的晶格常数差别不大，感兴趣的可以自行作为练习尝试。https://www.bigbrosci.com/2018/02/02/ex33/

（2）泛函类型对结构优化结果的影响

选用不同泛函优化得到的晶格常数会不一样，一般最常用的泛函是LDA和PBE（GGA下的最常用的泛函类型）。相对于实验值，LDA得到的晶格常数会偏小，PBE一般偏大。二者各有优点，很难直接说明谁更优。VASP提供PBE和LDA两套泛函下的POTCAR文件，请向组里师兄师姐寻求。或者你可以采用PBE下的POTCAR，然后在INCAR中设置GGA=CA，即可实现LDA泛函的计算；

The Journal of chemical physics, 2012, 136(13): 134704.

（3）高压下的结构优化？

前面是在0GPa下做的结构优化，如果我想得到10GPa下的的晶格常数该怎么办呢？只需要在INCAR中设置PSTRESS=100即可（100就是100kbar，10GPa=100kbar），其余参数不变。但注意，此时OSZICAR中的E0和F就包含了压强体积的能量项，即E0=U+pV（称为焓，Enthalpy，物理符号为H），比如这是bcc-Li在10GPa下的情况：

所以此时E0会比E更大一些，多的部分就是pV项。查看OUTCAR：

查看体系压强此时（结构的external pressure接近0了）：

做高压结构优化时，当前给定的初始POSCAR的体积和目标压强下对应的体积差别不能太大。比如，现有的结构是0GPa下优化得到的，如果想计算它在100GPa甚至300GPa下的焓或者晶格常数的话，直接设置PSTRESS=1000或者3000时做结构优化有时候会出错（往往发生在第二、三个离子步，可以自行查看压强随离子步变化寻找原因）。此时可以采用两种策略：方法一，分段优化。比如先做50GPa下的优化，再做100GPa下的优化，直到300GPa，通过施加较小的压强梯度，逐渐达到目标高压。方法二，对0GPa的结构进行晶格缩放，设置POSCAR中第二行参数，原先是1，可以设置为0.9，就相当于对POSCAR各个轴压缩0.9，此时POSCAR可能会更接近200GPa的晶格常数，设置PSTRESS=3000时就不太容易报错了，具体缩放系数依赖经验。

在做高压的结构优化时，由于relax过程中体积变化。而kmesh是基于原始poscar产生的。如果优化前后的晶格常数变大很大，可能会导致后面的晶格常数与kmesh不匹配（精度变化）。此时，最好基于优化成功的CONTCAR做二次relax或者静态计算以得到能量/焓。

（4）结构优化选原胞还是晶胞？

前面我们对金刚石Si做结构优化时，为了直观确认结构，取的是一个cubic的晶胞，里面有8个原子，这并非标准原胞。金刚石结构的标准原胞只有2个原子，借助vaspkit的602命令，它会读取当前路径下的POSCAR文件，并生成一个PRIMCELL.vasp文件，即为原胞。金刚石Si的原胞如下：

如果我们拿这个结构在0GPa下做结构优化（注意要重新生成KPOINTS，用0.025密度），得到的晶格常数是否会和前面用晶胞（8个原子）得到的结果一致？可以自行测试（此时应该比较平均单原子体积）。除以，此时k点密度要保持一致。如果你的KPOINTS采用的是Auto写法，那保持不变即可。如果是用vaspkit-102生成的类似7x7x7的写法，那需要重新对原胞生成相同密度下的Kmesh。

（5）晶格常数的温度依赖？

前面提到过，VASP在做DFT计算时，是不考虑温度对原子的激发的（热振动），所以我们优化得出的晶格常数是0K下的结果。那么如何得到有限温度下的晶格常数呢？目前主流有两种办法：

一是基于phonopy的QHA模块，做静态DFT+声子计算即可（请自行网上查阅phonopy的QHA模块）；

二是通过AIMD模拟，在NPT系综模拟结构在某一压强和温度下的动力学情况，并且模拟时间需要足够（通常得5-10ps以上）。在确认结构在MD过程中维持原有构型的前提下，对晶格常数做时间的统计平均。但VASP做NPT模拟时，体系很容易可能会变形（可能需要合理设置ISIF来控制形状），比如没法维持住cubic，三条轴夹角可能变化，对统计带来问题。此时可以考虑设置不同的晶格缩放系数，固定温度，在不同体积下做NVT模拟。每一个模拟达到5-10ps时，舍弃前几百步数据，做体系压强的系综平均，然后拟合P-V数据，得到该温度下所需压强下的体积。

不过对于大部分体系，300K和0K的晶格常数偏差都很小（可能不到0.5%），所以一般都是用DFT的基态0K结果就足以解释实验了。但对于某些高温高压（2000K）的情况，就必须考虑温度带来的效应了，并且对于金属体系，可能还需要考虑电子熵对力和压强的影响。

（6）结构优化可以自动帮助我们找到势能面全局最小值吗？

比如对于含有8个Si原子的体系，我随意放置这8个Si原子的位置，构建一个任意的胞，让VASP的结构优化算法自动帮我找到它在常压下最稳定的结构可以吗？不行，VASP的结构优化只能找到势能面局域最小值。搜索结构在势能面上的全局最小值事实上一直是一项挑战，由此衍生了“晶体（分子）结构预测”的研究方向。

Nature Reviews Materials, 2019, 4(5): 331-348.

Accounts of Chemical Research, 2022, 55(15): 2068-2076.

练习

现在你已经初步掌握结构优化计算的流程，你需要学会的技能是，给定一个初始结构（手动建模或者数据库下载），你可以自己生成合适的INCAR，KPOINTS和POTCAR并进行结构优化计算，通过查看OSZICAR文件和OUTCAR中原子受力以及压强信息，确保结构优化成功结束，并得到relaxed的结构（CONTCAR）以及相应的能量。建议做如下练习：

（1）构建一个晶格常数为4Å的面心立方（fcc）结构的Al单质（有四个原子），分别采用PBE泛函和LDA泛函对其进行结构优化（PSTRESS=0，常压），k点密度通过vaspkit-102生成，密度取0.025。并将relax后的两种晶格常数与实验对比，其中实验值自行查找ICSD或者Material Project数据库。如果不清楚这两个数据库，请自行b站搜索了解，或向师兄师姐寻求指导。

（2）将在步骤（1）中手动写好的含有4个原子的fcc-Al的结构转为fcc的标准原胞（含1个原子），通过vaspkit-602命令实现。然后做结构优化，KPOINTS用vaspkit重新生成，密度取0.025。将优化后的结构与晶胞优化得到的结构进行对比，比较单原子体积以及平均每原子能量。如果你的计算没有问题的话，两种结构得到的单原子体积和每原子能量在误差范围内应该是一样的。DFT的能量误差一般被认为是1meV以内（该误差源于ENCUT或者KPOINTS等计算精度）。

附：静态自洽计算和结构优化计算INCAR总结

静态自洽计算（SCF）和结构优化计算（relax）几乎是最常遇到的两个计算，因此非常有必要保存一个合适的INCAR。INCAR的写法因个人习惯而异，这里给出一个我个人常用的：

这里的参数在教程2.3.3和3.3中几乎都有说明。如果尚不清楚，请务必多浏览vaspwiki的参数介绍，一定要清楚这些常用参数的实际作用。

如果要做静态自洽计算，则需要设置NSW=0，此时#Ionic Relax部分的参数都不起作用。其余参数如LCHARG、LWAVE、ISPIN等根据任务需要进行设置。并且我这里加了NEDOS=2000以及LORBIT=11，是为了静态自洽计算完成后直接输出一个还不错的态密度数据DOSCAR，以及相应的分波态密度数据，有时候可以用来快速判断体系是否是金属、半导体。

如果需要做结构优化计算，只需要在静态自洽计算的基础上，开启NSW参数。有时候为了高精度，也可以设置PREC=Accurate，具体作用查官网。

注意，这里给出的INCAR中没有包含磁性、+U、SOC、IVDW等相关部分的参数。这些参数如何设置在后面的例子中会逐步介绍。因此这里的INCAR适用于那些没有磁性、不需要+U、不需要考虑VDW作用、不需要+U的体系的计算，且适用于常规的PBE或者LDA泛函。如果需要采用PBEsol或者SCAN等泛函，请自行查看vaspwiki的介绍。