干货 | 如何利用VASP计算任务?看看这10点分享

1、VASP 能够进行哪些过程的计算?怎样设置?

我们平时最常用的研究方法是做单点能计算,结构优化、从头计算的分子动力学和电子结构相关性质的计算。 一般我们的研究可以按照这样的过程来进行 如果要研究一个体系的最优化构型问题可以首先进行结构弛豫优化 然后对优化后的结构进 , 行性质计算或者单点能计算。如果要研究一个体系的热力学变化过程可以首先进行分子动力学过程模拟 然后在某个温度 , 或压强下进行性质计算或者单点能计算。 如果要研究一个体系的热力学结构变化可以首先在初始温度下进行 NVT 计算,然后进行分 子动力学退火,然后在结束温度下进行性质计算研究。

2、什么是单点能计算(single point energy)?如何计算?

跟其它软件类似,VASP 具有单点能计算的功能。也就是说,对一个给定的固定不变的结构(包括原子、分子、表面或体材料)能够计算其总能,即静态计算功能。 单点能计算需要的参数最少,最多只要在 KPOINTS 文件中设置一下合适的 K 点或者在 INCAR 文件中给定一个截断能 ENCUT 就可以了。还有一个参数就是电子步的收敛标准的 设置 EDIFF,默认值为 EDIFF=1E-4,一般不需要修改这个值。 具体来说要计算单点能,只要在 INCAR中设置 IBRION=-1 也就是让离子不移动就可以了。

3、什么是结构优化(structureoptimization)?如何计算?

结构优化又叫结构弛豫(structurerelax),是指通过对体系的坐标进行调整,使得其能量或内力达到最小的过程,与动力学退火不同,它是一种在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=AAPACO=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、如何进行结构优化参数设置?

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

初始结构包括原子堆积方式,和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。比如掺杂,表面吸附,空位等结构,初始原子的距离,角度等的设置需要有一定的 经验积累。DFT 计算短程强相互作用(相对于范德华力),如果初始距离设置过远(如超过 4 埃) ,则明显导致收敛很慢甚至得到不合理的结果。 比较好的设置方法可以参照键长。比如 CO 在 O 顶位的吸附,可以参照CO2 中 C-O 键 长来设置(如增长 20%)。也可以参照文献。记住一些常见键长,典型晶体中原子间距离等 参数,有助于提高初始结构设置的合理性。实在不行,可以先在小体系上测试,然后再放到 大体系中算。 弛豫参数弛豫参数对收敛速度影响很大,这一点在计算工作没有全部铺开时可能不会觉察到有 什么不妥,反正就给 NSW 设置个“无穷大”的数,最后总会有结果的。但是,时间是宝贵 的,恰当的设置 3 小时就收敛的结果,不恰当的设置可能要一个白天加一个黑夜。如果你赶 文章或者赶着毕业,你就知道这意味这什么。 结构优化分电子迭代和离子弛豫两个嵌套的过程。电子迭代自洽的速度, 有四个响很大 的因素:初始结构的合理性,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 步(默认的), 后面就会越来越快了。 总的说来,一般入门者,多看手册,多想多理解,多上机实践总结,比较容易提高到一 个熟练操作工的水平。 如果要想做到“精确打击” ,做到能在问题始发的时候就立刻采取有效措施来解决,就 需要回归基础理论和计算方法上来了。

 

7、如何用 VASP 计算铁磁、反铁磁和顺磁

顺磁,意味进行 non-spinpolarized 的计算,也就是 ISPIN=1。 铁磁,意味进行 spin-polarized 的计算,ISPIN=2,而且每个磁性原子的初始磁矩设置为一样的值,也就是磁性原子的 MAGMOM 设置为一样的值。对非磁性原子也可以设置成一 样的非零值(与磁性原子的一样)或零,最后收敛的结果,非磁性原子的 local 磁矩很小, 快接近 0,很小的情况,很可能意味着真的是非磁性原子也会被极化而出现很小的 local 磁 矩。 反铁磁,也意味着要进行 spin-polarized 的计算,ISPIN=2,这是需采用反铁磁的磁胞来 进行计算, 意味着此时计算所采用的晶胞不再是铁磁计算时的最小原胞。 比如对铁晶体的铁磁状态,你可以采用 bcc 的原胞来计算,但是在进行反铁磁的 Fe 计算,这是你需要采用 sc 的结构来计算,计算的晶胞中包括两个原子,你要设置一个原子的 MAGMOM 为正的,另一个原子的 MAGMOM 设置为负,但是它们的绝对值一样。因此在进行反铁磁的计算时,

 

应该确定好反铁磁的磁胞,以及磁序,要判断哪种磁序和磁胞是最可能的反铁磁状态,那只能是先做好各种可能的排列组合, 然后分别计算这些可能组合的情况, 最后比较它们的总能, 总能最低的就是可能的磁序。 同样也可以与它们同铁磁或顺磁的进行比较。 了解到该材料究竟是铁磁的、还是顺磁或反铁磁的。 亚铁磁,也意味要进行 spin-polarized 的计算,ISPIN=2,与反铁磁的计算类似,不同的 是原子正负磁矩的绝对值不是样大。非共线的磁性,那需采用专门的 non-collinear 的来进行 计算,除了要设置 ISPIN,MAGMOM 的设置还需要指定每个原子在 x,y,z 方向上的大小。 这种情况会复杂一些。举个例子来说,对于 Mn-Cu(001)c(2×2)这种体系,原胞里面有 2个 Mn 原子,那么你直 接让两个 Mn 原子的 MAGMOM 的绝对值一样,符号相反就可以了,再加上 ISPIN=2。这样就可以实现进行反铁磁的计算了

 

8、VASP在计算磁性的时候,oszicar 中得到的磁矩和 outcar 中得到各原子磁矩之和不一致在投稿的是否曾碰到有审稿人质疑,对于这个不一致你们一般是怎么解释?

 

答:OSZICAR 中得到的磁矩是 OUTCAR 中最后一步得到的总磁矩是相等的。总磁矩和各 原子的磁矩(RMT 球内的磁矩)之和之差就是间隙区的磁矩。因为有间隙区存在,不一致是正 常的。

 

9、磁性计算应该比较负责。你应该还使用别的程序计算过磁性,与vasp 结果比较是否一

致,对磁性计算采用的程序有什么推荐。 ps:由于曾使用 vasp 和 dmol 算过非周期体系磁性,结构对磁性影响非常大,因此使用这两个程序计算的磁性要一致很麻烦。还不敢确定到底是哪个程序可能不可靠。

 

答:如果算磁性,全电子的结果更精确,我的一些计算结果显示磁性原子对在最近邻的位置时,PAW 与 FPLAW 给出的能量差不一致,在长程时符合的很好。虽然并没有改变定性结论。 感觉 PAW 似乎不能很好地描述较强耦合。 我试图在找出原因, 主要使用 exciting 和 vasp 做比较。计算磁性推荐使用 FP-LAPW, FP-LMTO, FPLO 很吸引人(不过是商业的),后者是 O(N)算法。

10、VASP 中过渡态计算设置的一点体会

计算过渡态先要摆正心态,不急于下手。步骤如下:

 

(1)做模型,初态 IS 和终态 FS,分别结构优化到基态;

 

(2)线形插入 images: nebmake.pl POSCAR.IS POSCAR.FS N N 为 image 个数。

 

(3)nebmovie.pl,生成 movie.xyz。用 Xcrysden –xyz movie.xyz 反复观看动画,仔细检查过 程的合理性。这里要提醒,POSCAR.IS 和 POSCAR.FS 中原子坐标列表的顺序必须对应。

 

(4) INCAR,选 IOPT。 写 注意, 最好忘记 vasp 自带的 NEB, 而全部改用包含 vtstool 的 vasp. IBRION=3,POTIM=0 关闭 vasp 自带的 NEB 功能。

 

(5)过渡态计算第一个离子步最耗时,也最容易出问题,也是模型设计合理性检验的首要环节。所以可以选小一些的 ENCUT,可以不用考虑自旋(ISPIN=1),也不用考虑 DFT+U。 而且用最快最粗糙的算法(IOPT=3,其他默认)。

 

(6)带 vtstool 的 vasp-ClNEB(NEB)过渡态计算 ICHAIN=0 作为入口,这个也是默认的。 LCLIMB=TRUE 也是默认的。如果不要 climb image,可以设置 LCLIMB = False.

 

(7) 收敛判据 EDIFFG<0。 过渡态计算要以力为收敛判据, 而不是能量。 一般EDIFFG=-0.05 就可以接受, -0.02 或者-0.01更好。 但是作为开始的过渡态计算, 可以设置很宽的收敛条件, 如 EDIFFG=-1.

 

(8)初步过渡态收敛后,修改 INCAR 中的优化器(IOPT),并修改相应参数(参考 vtstool 官方论坛) ,EDIFFG 改小(如-0.05) ,然后运行 vfin.pl,这个脚本自动帮你准备在原来的基 础上继续运行新的过渡态计算(完成 cp CONTCAR POSCAR, 保留电荷密度和波函数的操 作) 。

 

(9)过渡态如何验算虚频呢?比如一个 6 层原子层的 slab 上表面吸附小分子。 slab 下部 3 层原子是固定的。 验算虚频的时 候,是不是还是固定下面三层原子,然后按照一般频率计算方法来算虚频?这样的话,可以移动的原子数在 20 数量级上,考虑三个自由度,及其组合,就有很多很多可能了。请问该 怎么设置这样的过渡态虚频计算呢?

原创文章,作者:菜菜欧尼酱,如若转载,请注明来源华算科技,注明出处:https://www.v-suan.com/index.php/2023/12/01/b4b76b2bc6/

(0)

相关推荐