VASP 结构优化、静态自洽、非自洽计算

时间:2022-07-28 19:32:41

第一步:结构优化(迟豫)


结构优化也叫结构迟豫。是指对整个输入体系的坐标进行调整,得到一个相对稳定的基态结构。结构优化分原子迟豫和电子迭代两个嵌套的过程,每次计算中都进行原子迟豫和电子迭代计算(电子迭代嵌套在原子迟豫中),达到原子迟豫收敛标准时进行下一步计算,直到达到自动中断或者最大原子迟豫步数。(通常以前后两次总*能之差或 原子所受最大的力 作为原子迟豫的收敛标准,电子自洽迭代计算以总能作为收敛标准,默认为10-4,因电子迭代的嵌套特性,以及电子迭代和原子迟豫的收敛速度不同,最终以原子迟豫最大步数作为强制停止参数)


通过结构优化使体系达到相对稳定结构,因此可以得到体系的总能量(哪个能量??)、晶体结构参数(可以得到包括晶格常数、结合能等)

同时得到CONTCAR(离子迟豫时,每步移动后的体系的晶格参数,和POSCAR内容相同。因此结构优化完成后,要使用输出的优化CONTCAR拷贝成POSCAR进行下一步的静态自洽计算,对电子结构做进一步调整 


结构优化的过程分析:

结构优化中原子迟豫和电子迭代嵌套,在每步原子迟豫过程中,对原子体系施加很小的力(应力还是H-F力??),决定原子如何移动,如何优化原子的位置(准牛顿法或共轭梯度法规定原子位置),并在前后两个位置判断原子迟豫收敛的标准,若未收敛,即在原子迟豫的每步内进行电子自洽迭代计算(也判断电子迭代收敛标准和最大迭代步数。电子停止迭代,并进行下一步原子迟豫),最终要达到原子迟豫的标准停止迟豫 或者 达到原子迟豫允许的最大步数强行停止。

在此过程中,若参数以及收敛值设置恰当,即可在最大允许步数内使得原子迟豫达到收敛标准,从而得到体系相对稳定的基态结构。(原子迟豫的收敛和电子自洽迭代的收敛需要设定一个适当的值,从而保证达到最大步数时体系已经达到相对稳定状态,但是收敛值也不能太小,会增加计算次数。

??每一步原子迟豫,是通过判断收敛标准开始下一步还是位置达到??——>应该是位置,否则最终的原子迟豫停止无法判断


弛豫计算的两种方式:

迟豫的方式,即原子位置的优化算法,即是原子迟豫每步的原子位置,在每步内做电子迭代,比较两步的能量或力判断收敛。

准牛顿方法(quasi-Newton RMM-DIIS)和共轭梯度法(CG)两种。准牛顿方法计算速度较快,适合于初始结构与平衡结构(势能面上全局最小值)比较接近的情况,而CG方法慢一些,找到全局最小的可能性也要大一些。


主要涉及到的参数:

(1)原子迟豫相关参数:

EDIFFG(原子迟豫收敛标准);

NSW(原子迟豫最大步数,当计算在最大步数之内达到收敛时计算自动中断,而最大步数内没有收敛的话结构迟豫将被强制中止(平常计算步数不会超过100步,超过100步可能是计算的体系出了问题)。同时在每一步原子迟豫内,电子进行自洽计算;每一步内,精确计算原子所受的H-F力和应力。);

IBRION(决定原子如何移动和迟豫,即迟豫计算的方式:-1 表示原子不移动(NSW为0或-1时,默认取-1);1、2、3 表示以不同的算法优化原子位置。准牛顿方法(IBRION=1)计算速度较快,适合于初始结构与平衡结构(势能面上全局最小值)比较接近的情况,而CG方法(IBRION=2)慢一些,找到全局最小的可能性也要大一些);

POTIM(IBRION=1,2,3时表示原子每步移动的大小,默认为0.5。即是使用准牛顿法、共轭梯度法、最速下降法优化原子位置时需要设置的步长)

(2)电子自洽迭代相关参数:

NELM(允许电子自洽迭代的最大步数,默认60);

NELMIN(电子迭代的最小步数。默认为2,一般不设置);

EDIFF(电子自洽迭代循环的收敛标准,以总能衡量,默认为10-4

(3)结构优化相关参数:

ISIF(可以控制迟豫时晶胞的变化情况,例如原胞的形状和体积等,一般选2、3)


结构优化收敛性分析1:

结构优化,或者叫弛豫,是后续计算的基础。其收敛性受两个主要因素影响:初始结构的合理性和弛豫参数的设置。

(1)初始结构

初始结构包括原子堆积方式,和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。比如掺杂,表面吸附,空位等结构,初始原子的距离,角度等的设置需要有一定的经验积累。DFT计算短程强相互作用(相对于范德华力),如果初始距离设置过远(如超过4埃),则明显导致收敛很慢甚至得到不合理的结果。

比较好的设置方法可以参照键长。比如CO在O顶位的吸附,可以参照CO2中C-O键长来设置(如增长20%)。也可以参照文献。记住一些常见键长,典型晶体中原子间距离等参数,有助于提高初始结构设置的合理性。实在不行,可以先在小体系上测试,然后再放到大体系中算。

(2)弛豫参数

弛豫参数对收敛速度影响很大,这一点在计算工作没有全部铺开时可能不会觉察到有什么不妥,反正就给NSW设置个“无穷大”的数,最后总会有结果的。但是,时间是宝贵的,恰当的设置3小时就收敛的结果,不恰当的设置可能要一个白天加一个黑夜。如果你赶文章或者赶着毕业,你就知道这意味这什么。


结构优化的收敛性分析2:

结构优化分电子迭代和离子弛豫两个嵌套的过程。

电子迭代自洽的速度,有四个响很大的因素:初始结构的合理性,k点密度,是否考虑自旋和高斯展宽(SIGMA);

离子弛豫的收敛速度,有三个很大的影响因素:弛豫方法(IBRION),步长(POTIM)和收敛判据(EDIFFG)。

一般来说,针对理论催化的计算,初始结构都是不太合理的。因此一开始采用很粗糙的优化(EDIFF=0.001,EDIFFG= -0.2),很低的K点密度(Gamma),不考虑自旋就可以了,这样NSW<60的设置就比较好。其它参数可以默认。

经过第一轮优化,就可以进入下一步细致的优化了。就我的经验,EDIFF=1E-4,EDIFFG=-0.05,不考虑自旋,IBRION=2,其它默认,NSW=100;跑完后可以设置IBRION = 1,减小OPTIM(默认为0.5,可以设置0.2)继续优化。

优化的时候让它自己闷头跑是不对的,经常看看中间过程,根据情况调节优化参数是可以很好的提高优化速度。这个时候,提交两个以上的任务排队是好的方式,一个在调整的时候,下一个可以接着运行,不会因为停下当前任务导致机器空闲。

无论结构优化还是静态自洽,电子步的收敛也常常让新手头痛。如果电子步不能在40步内收敛,要么是参数设置的问题,要么是初始模型太糟糕(糟糕的不是一点点)。

静态自洽过程电子步不收敛一般是参数设置有问题。这个时候,改变迭代算法(ALGO),提高高斯展宽(SIGMA增加),设置自洽延迟(NELMDL)都是不错的方法。对于大体系比较难收敛的话,可以先调节AMIN,BMIX跑十多步,得到电荷密度和波函数,再重新计算。实在没办法了,可以先放任它跑40步,没有收敛的迹象的话,停下来,得到电荷密度和波函数后重新计算。一般都能在40步内收敛。

对于离子弛豫过程,不调节关系也不大。开始两个离子步电子迭代可能要跑满60步(默认的),后面就会越来越快了。

总的说来,一般入门者,多看手册,多想多理解,多上机实践总结,比较容易提高到一个熟练操作工的水平。

如果要想做到“精确打击”,做到能在问题始发的时候就立刻采取有效措施来解决,就需要回归基础理论和计算方法上来了。


迟豫收敛判据(强制停止)的选择:

结构弛豫的判据一般有两种选择:能量和力。这两者是相关的,理想情况下,能量收敛到基态,力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大,力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的,能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。

到底是以能量为收敛判据,还是以力为收敛判据呢?关心能量的人,觉得以能量为判据就够了;关心力相关量的人,用力作为收敛标准。

对于超胞体系的结构优化,文献大部分采用Gamma点做单点优化。这个时候即使采用力为判据(EDIFFG=-0.02),在做静态自洽计算能量的时候,会发现,原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。这个时候,是不是该怀疑对超胞仅做Gamma点结构优化的合理性呢?是不是要提高K点密度再做结构优化呢。

在我看来,这取决于所研究的问题的复杂程度。我们的计算从原胞开始,到超胞,到掺杂结构,到吸附结构,到反应和解离。每一步都在增加复杂程度。结构优化终点与初始结构是有关的,如果遇到对初始结构敏感的优化,那就头疼了。而且,还要注意到,催化反应不仅与原子本身及其化学环境有关,还会与几何构型有关。气固催化反应过程是电子的传递过程,也是分子拆分与重新组合的过程。如果优化终点的构型不同,可能会导致化学反应的途径上的差异。仅从这一点来看,第一性原理计算的复杂性,结果上的合理性判断都不是手册上写的那么简单。

对于涉及构型敏感性的结构优化过程,我觉得,以力作为收敛判据更合适。而且需要在Gamma点优化的基础上再提高K点密度继续优化,直到静态自洽计算时力达到收敛标准的。



促进收敛的方法:

费米面附近能级电子分布的smearing是一种促进收敛的有效方法,可能产生物理意义不明确的分数占据态情况,不过问题不大。在INCAR文件中以ISMEAR来设置。一般来说K点只有一两个的时候采用ISMEAR=0,金属体材料用ISMEAR=1或2,半导体材料用ISMEAR=-5等等。

不过有时电子步收敛速度依然很慢,还需要设置一些算法控制选项,例如设置ALGO=Very_Fast,减小真空层厚度,减少K点数目等。


带限制条件的迟豫计算:

有的时候我们需要一些带限制条件的弛豫计算,例如冻结部分原子、限制自旋的计算等等。冻结部分原子可以在POSCAR文件中设置selectivedynamic来实现。自旋多重度限制可以在INCAR中以NUPDOWN选项来设置。另外ISIF选项可以控制弛豫时的晶胞变化情况,例如晶胞的形状和体积等。


第二步:静态自洽计算


静态:原子位置保持不动,不再进行原子迟豫(???是否正确);
自洽:电子再进行自洽计算;
因此,静态自洽计算是在结构优化的基础上,在体系能量达到 较低、体系 较稳定 的情况下固定原子的位置坐标,再对体系中的电子进行调整,以达到体系的 最低能量

因此,静态自洽计算前需要提供 较稳定体系 的晶格结构信息(即结构优化完的CONTCAR),从而通过电子自洽计算 (通过自洽迭代求解薛定谔方程(微观中描述电子的状态,相当于宏观的运动方程))  完整地计算出 体系基态下费米能级(准确值,The fermi level location is accurate only in self-consistent calculation.)、 电子的波函数、电荷密度等信息,可以直接分析原子间的键合作用,也可以在非自洽之后进一步分析晶体的电子结构和材料的相关性质

静态自洽计算的准备工作:

使用结构优化输出的CONTCAR拷贝成新的POSCAR

静态自洽主要涉及的参数:

NSW=0  不再进行原子迟豫
LWAVE=.TRUE  (需要在WAVECAR中输出波函数)
LCHARGE=.TRUE    (需要在CHG和CHGCAR中输出电荷密度)


注释:第一步和第二部都含有电子迭代的自洽计算。第三步开始都不是非自洽计算。

第三步:非自洽计算


非自洽计算是在自洽基础上改变k点(重新生成k点)等参数,根据不同需要选取能量或势函数或电子密度作为初始值,进行非自洽迭代计算,可用于求解DOS,能带(电子结构分析)或者光学等其他性质

——>非自洽计算可以得到EIGENCAR文件,与syml文件通过脚本做能带和态密度性质。


能带的非自洽计算在侯博士的手册中有讲。下面只是一些重要的地方。

重新生成k点:

syml文件:需要静态自洽计算之后输出的OUTCAR中的基矢和倒基矢、费米能级(准确值);  取高对称点路径以及分割点数。
具体以以下格式书写:
6   
30  20  30  20  30
G 0.0    0.0    0.0
M 0.0    0.5    0.0
K -0.333333333333   0.6666666666667   0.0
G 0.0    0.0    0.0
A 0.0    0.0    0.5
L 0.0    0.5    0.5
0.000000000 -4.052756964  0.000000000     0.142458646 -0.246745613  0.000000000
3.509790486  2.026378482  0.000000000     0.284917292  0.000000000  0.000000000
0.000000000  0.000000000  5.572964774     0.000000000  0.000000000  0.179437703
-20.0  20.0
0.0

上面每一行意思:
#############comment
1 line: nhighk   !how many high symmetry k-points
2     : ndiv(i)  !number of grid points between two symmetry k-points
3     : labhk(i), phighk(i,j)   !label of high sysmmetry k-points, coordinates
......
9     : a(i,j),b(i,j)   !direct & reciprocal lattice vectors
......
12    : emin, emax   ! energy width of band plot
13    : efermi       ! fermi level
##############
\\\ G--M--K--G--A---L   NKPT=131
\\\ 1--31--51--81--101---131
--------------------------------------------
Hex
M:0.0    0.5    0.0
G:0.0    0.0    0.0
H:-1/3   2/3    0.5
K:-1/3   2/3    0.0
L:0.0    0.5    0.5
A:0.0    0.0    0.5

kpoints有几种生成方式:自动模式,line模式(用于能带计算),全手动模式,用SYML作为输入文件使用gk.f编译的程序生成。

不同的方式有不同的用途:
auto的方式,除了不能用在能带计算中,其他的都可以用到。
line模式,只用在能带的计算中。
syml的模式,也只是用在能带计算中。这个的功能等同于line 模式,只是自己手动产生而已。


涉及的参数设置:

ISTART=1  (接着计算)
ICHAGE=11   (读入自洽的CHGCAR,从而开始能带或者态密度的非自洽计算)
NSW=0  (不再进行原子迟豫)


VASP 结构优化、静态自洽、非自洽计算

注:结构优化和静态自洽计算都设置ISTART=0,因此都是从头开始新的计算,因此会刷新各个文件中的数据,所以拷贝文件要全。同理,非自洽计算设置ISTART=1,即进行接着算,也要注意文件内容。