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泛函的计算；
TABLE II. Lattice constants (A) for SLC34.
Solid Structurea Expt.b SOGGA PBEsol LSDA HSE
C A4 3.555 3.556 3.561 3.537 3.553
Si A4 5.422 5.430 5.442 5.410 5.444
Ge A4 5.644 5.672 5.692 5.634 5.701
SiC B3 4.348 4.374 4.380 4.355 4.372
BP B3 4.527 4.528 4.540 4.509 4.543
BAs B3 4.764 4.774 4.788 4.750 4.794
AIP B3 5.450 5.461 5.472 5.436 5.472
AlAs B3 5.649 5.670 5.682 5.639 5.691
AISb B3 6.126 6.131 6.146 6.079 6.146
GaNc B4 3.179 3.190 3.197 3.167 3.198The 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下的情况：
N E dE eps ncg
DAV: 1 -0.352027129522E+01 -0.84064E-02 -0.34789E-03 1104
DAV: 2 -0.352029216784E+01 -0.20873E-04 -0.89364E-04 1032
DAV: 3 -0.352027508702E+01 0.17081E-04 -0.28088E-04 1696
DAV: 0.352023590319E+01 0.39184E-04 -0.57721E-05 1648
DAV: 5 0.352023643290E+01 -0.52970E-06 -0.12432E-06 752
5 F=-.17909779E+01 E0= -.17907357E+01 E =-.499574E+00所以此时E0会比E更大一些，多的部分就是pV项。查看OUTCAR：
FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)
free energy TOTEN -3.52023643 BeV
energy without entropy= -3.51950978 energy(sigma->0) -3.51999422
enthalpy is TOTEN -1.79097795 ev P V= 1.72925848查看体系压强此时（结构的external pressure接近0了）：
Total 1.72346 1.72346 1.72346 0.00000 0.00000 -0.00000
in kB 99.66487 99.66487 99.66487 0.00000 0.00000 -0.00000
external pressure -0.34 kB Pullay stress 100.00 kB做高压结构优化时，当前给定的初始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的原胞如下：
Primitive Cell
1.000000
0.00000000000000 2.73536078266694 2.73536078266694
2.73536078266694 0.00000000000000 2.73536078266694
2.73536078266694 2.73536078266694 0.00000000000000
Si
2
DIRECT
0.7500000000000000 0.7500000000000000 0.7500000000000000 Si1
0.5000000000000000 0.5000000000000000 0.5000000000000000 Si2如果我们拿这个结构在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的写法因个人习惯而异，这里给出一个我个人常用的：
ISTART = 0
ICHARG = 2
LCHARG = .FALSE.
LWAVE = .FALSE.
ENCUT = 520
LREAL = Auto
#PREC = Accurate
ISPIN = 1
# Ionic Relaxation
NSW = 100
IBRION = 2
ISIF = 3
PSTRESS = 0
EDIFFG = -0.01
# SCF
ISMEAR = 0
SIGMA = 0.05
NEDOS = 2000
NELM = 60
EDIFF = 1E-06
LORBIT = 11
NCORE = 4
KPAR = 4这里的参数在教程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的介绍。