VASP 计算的过程遇到的问题
01、第一原理计算的一些心得
(1)第一性原理其实是包括基于密度泛函的从头算和基于Hartree-Fock自洽计算的从头算,前者以电子密度作为基本变量(霍亨伯格-科洪定理),通过求解Kohn-Sham方程,迭代自洽得到体系的基态电子密度,然后求体系的基态性质;后者则通过自洽求解Hartree-Fock方程,获得体系的波函数,求基态性质;
评述:K-S方程的计算水平达到了H-F水平,同时还考虑了电子间的交换关联作用。 (2)关于DFT中密度泛函的Functional,其实是交换关联泛函 包括LDA,GGA,杂化泛函等等
一般LDA为局域密度近似,在空间某点用均匀电子气密度作为交换关联泛函的唯一变量,多数为参数化的CA-PZ方案;
GGA为广义梯度近似,不仅将电子密度作为交换关联泛函的变量,也考虑了密度的梯度为变量,包括PBE,PW,RPBE等方案,BLYP泛函也属于GGA; 此外还有一些杂化泛函,B3LYP等。 (3)关于赝势
在处理计算体系中原子的电子态时,有两种方法,一种是考虑所有电子,叫做全电子法,比如WIEN2K中的FLAPW方法(线性缀加平面波);此外还有一种方法是只考虑价电子,而把芯电子和原子核构成离子实放在一起考虑,即赝势法,一般赝势法是选取一个截断半径,截断半径以内,波函数变化较平滑,和真实的不同,截断半径以外则和真实情况相同,而且赝势法得到的能量本征值和全电子法应该相同。
赝势包括模守恒和超软,模守恒较硬,一般需要较大的截断能,超软势则可以用较小的截断能即可。另外,模守恒势的散射特性和全电子相同,因此一般红外,拉曼等光谱的计算需要用模守恒势。
赝势的测试标准应是赝势与全电子法计算结果的匹配度,而不是赝势与实验结果的匹配度,因为和实验结果的匹配可能是偶然的。 (4)关于收敛测试
(a)Ecut,也就是截断能,一般情况下,总能相对于不同Ecut做计算,当Ecut增大时总能变化不明显了即可;然而,在需要考虑体系应力时,还需对应力进行收敛测试,而且应力相对于Ecut的收敛要比总能更为苛刻,也就是某个截断能下总能已经收敛了,但应力未必收敛。
(b)K-point,即K网格,一般金属需要较大的K网格,采用超晶胞时可以选用相对较小的K网格,但实际上还是要经过测试。 (5)关于磁性
一般何时考虑自旋呢?举例子,例如BaTiO3中,Ba、Ti和O分别为+2,+4和-2价,离子全部为各个轨道满壳层的结构,就不必考虑自旋了;对于BaMnO3中,由于Mn+3价时d轨道还有电子,但未满,因此需考虑Mn的自旋,至于Ba和O则不必考虑。其实设定自旋就是给定一个原子磁矩的初始值,只在刚开始计算时作为初始值使用,具体的可参照磁性物理。
(6)关于几何优化
包括很多种了,比如晶格常数和原子位置同时优化,只优化原子位置,只优化晶格常数,还有晶格常数和原子位置分开优化等等。
在PRL一篇文章中见到过只优化原子位置,晶格常数用实验值的例子(PRL 100, 186402 (2008));也见到过晶格常数先优化,之后固定晶格常数优化原子位置的情况;更多的情况则是Full geometry optimization。 一般情况下,也有不优化几何结构直接计算电子结构的,但是对于缺陷形成能的计算则往往要优化。
(7)关于软件
软件大致分为基于平面波的软件,如CASTEP、PWSCF和ABINIT等等,计算量大概和体系原子数目的三次方相关;还有基于原子轨道线性组合的软件(LCAO),比如openmx,siesta,dmol等,计算量和体系原子数目相关,一般可模拟较多原子数目的体系。
VASP是使用赝势和平面波基组,进行从头量子力学分子动力学计算的软件包,它基于CASTEP 1989版开发。VAMP/VASP中的方法基于有限温度下的局域密度近似(用自由能作为变量)以及对每一MD步骤用有效矩阵对角方案和有效Pulay混 合求解瞬时电子基态。这些技术可以避免原始的Car-Parrinello方法存在的一切问题,而后者是基于电子、离子运动方程同时积分的方法。离子和电 子的相互作用超缓Vanderbilt赝势(US-PP)或投影扩充波(PAW)方法描述。两种技术都可以相当程度地减少过渡金属或第一行元素的每个原子 所必需的平面波数量。力与张量可以用VAMP/VASP很容易地计算,用于把原子衰减到其瞬时基态中。
02、VASP程序的亮点:
1. VASP使用PAW方法或超软赝势,因此基组尺寸非常小,描述体材料一般需要每原子不超过100个平面波,大多数情况下甚至每原子50个平面波就能得到可靠结果。
2. 在平面波程序中,某些部分代码的执行是三次标度。在VASP中,三次标度部分的前因子足可忽略,导致关于体系尺寸的高效标度。因此可以在实空间求解势的非局域贡献,并使正交化的次数最少。当体系具有大约2000个电子能带时,三次标度部分与其它部分可比,因此VASP可用于直到4000个价电子的体系。 3. VASP使用传统的自洽场循环计算电子基态。这一方案与数值方法组合会实现有效、稳定、快速的Kohn-Sham方程自洽求解方案。程序使用的迭代矩阵对角化方案(RMM-DISS和分块Davidson)可能是目前最快的方案。
4. VASP包含全功能的对称性代码,可以自动确定任意构型的对称性。
5. 对称性代码还用于设定Monkhorst-Pack特殊点,可以有效计算体材料和对称的团簇。Brillouin区的积分使用模糊方法或四面体方法。四面体方法可以用Blöchl校正去掉线性四面体方法的二次误差,实现更快的k点收敛速度。
03、VASP 5.2的新功能:
1. 大规模并行计算需要较少的内存。
2. 加入新的梯度校正泛函AM05和PBEsol;用标准PBE POTCAR文件提供新泛函;改善了单中心处理。
3. 离子位置和格矢中加入有限差分,从而得到二阶导,用于计算原子间力常数和声子(需要超晶胞近似),和弹性常数。计算中自动考虑对称性。 4. 离子位置和静电场中加入线性响应,从而得到二阶导,用于计算原子间力常数和声子(需要超晶胞近似),Born有效电荷张量,静态介电张量(电子和离子贡献),内应变张量,压电张量(电子和离子贡献)。线性响应只能用于局域和半局域泛函。
5. 精确的非局域交换和杂化泛函:Hartree-Fock方法;杂化泛函,特别是PBE0和HSE06;屏蔽交换;(实验性的)简单模型势GW-COHSEX,用于经验的屏蔽交换内核;(实验性的)
杂化泛函B3LYP。
6. 通过本征态求和计算含频介电张量:使用粒子无关近似,或通过GW的随机相近似。可用于局域,半局域,杂化泛函,屏蔽交换,和Hartree-Fock。
7. 完全含频GW,速度达到等离子极点模型:单发G0W0;在G和W中迭代本征矢直至自洽;(实验性的)迭代G(也可以选W)本征矢的自洽GW;(实验性的)对相关能使用RPA近似的GW总能量;用LDA计算G和W的顶点校正(局域场效应),仅能用于非自旋极化的情况;(实验性的)W的多体顶点校正,仅能用于非自旋极化的情况。
8. 实验性的功能:用TD-HF和TD-杂化泛函求解Cassida方程(仅能用于非自旋极化的Tamm-Dancoff近似);GW顶点的Bethe-Salpeter(仅能用于非自旋极化的Tamm-Dancoff近似)。
1、VASP能够进行哪些过程的计算?怎样设置?
我们平时最常用的研究方法是做单点能计算,结构优化、从头计算的分子动力学和电子结构相关性质的计算。
一般我们的研究可以按照这样的过程来进行
如果要研究一个体系的最优化构型问题可以首先进行结构弛豫优化,然后对优化后的结构进行性质计算或者单点能计算。
如果要研究一个体系的热力学变化过程可以首先进行分子动力学过程模拟,然后在某个温度或压强下进行性质计算或者单点能计算。
如果要研究一个体系的热力学结构变化可以首先在初始温度下进行NVT计算,然后进行分子动力学退火,然后在结束温度下进行性质计算研究。
2、什么是单点能计算(single point energy)?如何计算?
跟其它软件类似,VASP具有单点能计算的功能。也就是说,对一个给定的固定不变的结构(包括原子、分子、表面或体材料)能够计算其总能,即静态计算功能。
单点能计算需要的参数最少,最多只要在KPOINTS文件中设置一下合适的K点或者在INCAR文件中给定一个截断能ENCUT就可以了。还有一个参数就是电子步的收敛标准的设置EDIFF,默认值为EDIFF=1E-4,一般不需要修改这个值。
具体来说要计算单点能,只要在INCAR中设置IBRION=-1也就是让离子不移动就可以了。
3、什么是结构优化(structure optimization)?如何计算?
结构优化又叫结构弛豫(structure relax),是指通过对体系的坐标进行调整,使得其能量或内力达到最小的过程,与动力学退火不同,它是一种在0K下用原子间静力进行优化的方法。可以认为结构优化后的结构是相对稳定的基态结构,能够在实验之中获得的几率要大些(当然这只是理论计算的结果,必须由实验来验证)。
一般要做弛豫计算,需要设置弛豫收敛标准,也就是告诉系统收敛达成的判据(convergence break condition),当系统检测到能量变化减小到一个确定值时例如EDIFFG=1E-3时视为收敛中断计算,移动离子位置尝试进行下一步计算。EDIFFG这个值可以为负,例如EDIFFG=-0.02,这时的收敛标准是当系统发现所有离子间作用力都小于给定的数值,如0.02eV/A时视为收敛而中断。
弛豫计算主要有两种方式:准牛顿方法(quasi-Newton RMM-DIIS)和共轭梯度法(CG)两种。准牛顿方法计算速度较快,适合于初始结构与平衡结构(势能面上全局最小值)比较
接近的情况,而CG方法慢一些,找到全局最小的可能性也要大一些。选择方法为IBRION=1时为准牛顿方法而IBRION=2时为CG方法。
具体来说要做弛豫计算,设置IBRION=1或者2就可以了,其它参数根据需要来设置。NSW是进行弛豫的最大步数,例如设置NSW=100,当计算在100步之内达到收敛时计算自动中断,而100步内没有达到收敛的话系统将在第100步后强制中止(平常计算步数不会超过100步,超过100步可能是计算的体系出了问题)。参数通常可以从文献中发现,例如收敛标准EDIFFG等。
有的时候我们需要一些带限制条件的弛豫计算,例如冻结部分原子、限制自旋的计算等等。冻结部分原子可以在POSCAR文件中设置selective dynamic来实现。自旋多重度限制可以在INCAR中以NUPDOWN选项来设置。另外ISIF选项可以控制弛豫时的晶胞变化情况,例如晶胞的形状和体积等。
费米面附近能级电子分布的smearing是一种促进收敛的有效方法,可能产生物理意义不明确的分数占据态情况,不过问题不大。在INCAR文件中以ISMEAR来设置。一般来说K点只有一两个的时候采用ISMEAR=0,金属体材料用ISMEAR=1或2,半导体材料用ISMEAR=-5等等。不过有时电子步收敛速度依然很慢,还需要设置一些算法控制选项,例如设置ALGO=Very_Fast,减小真空层厚度,减少K点数目等。
弛豫是一种非常有效的分析计算手段,虽然是静力学计算但是往往获得一些动力学得不到的结果。
4、vasp的分子动力学模拟
vasp做分子动力学的好处,由于vasp是近些年开发的比较成熟的软件,在做电子scf速度方面有较好的优势。缺点:可选系综太少。
尽管如此,对于大多数有关分子动力学的任务还是可以胜任的。主要使用的系综是 NVT 和 NVE。
一般做分子动力学的时候都需要较多原子,一般都超过100个。
当原子数多的时候,k点实际就需要较少了。有的时候用一个k点就行,不过这都需要严格的测试。通常超过200个原子的时候,用一个k点,即Gamma点就可以了。
INCAR:
EDIFF 一般来说,用1E-4 或者1E-5都可以,这个参数只是对第一个离子步的自洽影响大一些,对于长时间的分子动力学的模拟,精度小一点也无所谓,但不能太小。 IBRION=0 分子动力学模拟
IALGO=48 一般用48,对于原子数较多,这个优化方式较好。 NSW=1000 多少个时间步长。
POTIM=3 时间步长,单位fs, 通常1到3. ISIF=2 计算外界的压力.
NBLOCK= 1 多少个时间步长,写一次CONTCAR,CHG和CHGCAR,PCDAT.
KBLOCK=50 NBLOCK*KBLOCK 个步长写一次XDATCAR.(个离子步写一次PCDAT.) ISMEAR=-1 费米迪拉克分布. SIGMA =0.05 单位:电子伏
NELMIN=8 一般用6到8, 最小的电子scf数.太少的话,收敛的不好. LREAL=A
APACO=10 径向分布函数距离, 单位是埃. NPACO=200 径向分布函数插的点数.
LCHARG=F 尽量不写电荷密度,否则CHG文件太大. TEBEG=300 初始温度.
TEEND=300 终态温度。 不设的话,等于TEBEG.
SMASS=-3 NVE ensemble;-1 用来做模拟退火。大于0 NVT 系综。正确:SMASS=1,2,3 是没有区别的。都是NVT ensemble。SMASS只要是大于0就是NVT系综。 CONTCAR是每个离子步之后都会写出来的,但是会用新的把老的覆盖 CHG是在每10个离子步写一次,不会覆盖 CHGCAR是在任务正常结束之后才写的。
5、收敛判据的选择
结构弛豫的判据一般有两中选择:能量和力。这两者是相关的,理想情况下,能量收敛到基态,力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大,力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的,能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。
到底是以能量为收敛判据,还是以力为收敛判据呢?关心能量的人,觉得以能量为判据就够了;关心力相关量的人,没有选择,只能用力作为收敛标准。对于超胞体系的结构优化,文献大部分采用Gamma点做单点优化。这个时候即使采用力为判据(EDIFFG=-0.02),在做静态自洽计算能量的时候,会发现,原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。这个时候,是不是该怀疑对超胞仅做Gamma点结构优化的合理性呢?是不是要提高K点密度再做结构优化呢。
在我看来,这取决于所研究的问题的复杂程度。我们的计算从原胞开始,到超胞,到掺杂结构,到吸附结构,到反应和解离。每一步都在增加复杂程度。结构优化终点与初始结构是有关的,如果遇到对初始结构敏感的优化,那就头疼了。而且,还要注意到,催化反应不仅与原子本身及其化学环境有关,还会与几何构型有关。气固催化反应过程是电子的传递过程,也是分子拆分与重新组合的过程。如果优化终点的构型不同,可能会导致化学反应的途径上的差异。仅从这一点来看,第一性原理计算的复杂性,结果上的合理性判断都不是手册上写的那么简单。
对于涉及构型敏感性的结构优化过程,我觉得,以力作为收敛判据更合适。而且需要在Gamma点优化的基础上再提高K点密度继续优化,直到静态自洽计算时力达到收敛标准的。
6、结构优化参数设置
结构优化,或者叫弛豫,是后续计算的基础。其收敛性受两个主要因素影响:初始结构的合理性和弛豫参数的设置 初始结构