cineb过渡态频率微信超出验证频率限制问题求助

VASP&资料搜集_minmin3683609_新浪博客
VASP&资料搜集
(1)这个范德华力是在吸附的时候用的到,其他情况,比如晶格常数优化是没用的。因为有些气体吸附在体系表面时(弱相互作用),算到的吸附能与实验值偏差大,原因是没有考虑长程的范德华作用力。最近在看这方面的文献,希望对你有帮助:THE
JOURNAL OF CHEMICAL PHYSICS 134, 11)
(2)应该是VASP5.0以上版本有VDW校正,使用PBE赝势,在INCAR文件中加入LVDW=.T.就可以了。
范德瓦耳斯力(van der Waals'
force)在化学中指分子之间非定向的、无饱和性的、较弱的相互作用力,根据荷兰物理学家约翰内斯·范德瓦耳斯命名。范德瓦耳斯力是一种电性引力,但它比化学键或氢键弱得多,通常其能量小于5kJ/mol。分子的大小和范德瓦耳斯力的大小成正比。
范德瓦耳斯力的主要来源有三种机制:
极性分子与极性分子之间的永久偶极矩相互作用,称为“取向力”。
极性分子对非极性分子有极化作用,使之产生诱导偶极矩,永久偶极矩与其诱导出的偶极矩相互作用,称为“诱导力”。
一对非极性分子本身由于电子的概率运动,可以相互配合产生一对方向相反的瞬时偶极矩,这一对瞬时偶极矩相互作用,称为“色散力”。这种机制是范德瓦耳斯力的主要来源,1930年由F.W.伦敦首先根据量子力学原理给出解释,因此也称为“伦敦力”。
范德瓦耳斯力的大小会影响物质尤其是分子晶体的熔点和沸点,通常分子的相对分子质量越大,范德瓦耳斯力越大。水(氧化氢)比硫化氢的相对分子质量小,因此范德瓦耳斯力比后者弱,但由于水分子间存在更强的氢键,熔沸点反而更高。壁虎能够在墙及各种表面上行走,便是因为脚上极细致的匙突(spatulae)和接触面产生的范德瓦耳斯力所致。
(转)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, 08))
;也见到过晶格常数先优化,之后固定晶格常数优化原子位置的情况;更多的情况 则是 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
-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、优化结果对初始结构和“优化路径”的依赖
原子吸附问题不大,但是小分子吸附,存在初始构型上的差异。slab 上水平放置,还是 垂直放置,可能导致收敛结果上的差异。根据
H-K 理论,理想情况下,优化得到的应该是 全局最小,但在数值计算的时候可能经常碰到不是全局最小的情况。实际操作中发现,多个
不同初始结构优化收敛后在能量和结构上存在一定差异。 为了加快收敛速度, 特别是对于表面-分子吸附结构, 初始放松约束, 比如
EDIFF=1E-3, EDIFFG=-0.3,NSW=30 可能是很好的设置。但是下面的情况应当慎重: EDIFF=1E-3;
EDIFFG=-0.1;!或者更小 NSW=500; !或者更大
电子步收敛约束较小,而离子步约束偏大,离子步数又很多,这种情况下,可能导致的 结果是结构弛豫到严重未知的区间。
再在这个基础上提高约束来优化,可能就是徒劳的了——结果不可逆转的偏向不正常的 区间。
好的做法,是对初始结构做比较松弛的约束,弛豫离子步 NSW 应该限制在一个较小的 数值内。EDIFF=1E-3 的话,EDIFFG
也最好是偏大一些,如-0.3 而不是-0.1. 这样可以在较 少的步数内达到初步收敛。
对于远离基态的初始结构,一开始在非常松弛的约束下跑若干离子步,时间上带来的好 处是很大的。对于 100 个原子的体系用 vasp 做
Gamma 点优化,如果一开始就是正常优化 (
EDIFF=1E-4,EDIFFG=-0.02)设置,开始十个离子步可能都要花上几个小时。如果这个
时候才发现输入文件有错误,那下午的时间就白费了,顺便带上晚上机器空转。 所以,我习惯的做法,是在初始几步优化后,会用 xcrysden
检查一下 XDATCAR 中的 数据,用 xdat2xyz.pl 生成
movie.xyz,然后看看弛豫过程是不是按照设想的那样。后续过程 跑完一个收敛过程,就再检查一下
movie.xyz。如此这般,才放心的展开后续计算。
8、目的导向的结构优化
结构优化到这个阶段,是高级的了。为了得到特定结构,或者为了验证某些猜想,需要设计 合理的初始结构,然后在这个基础上小心优化,比如
POTIM=0.1 跑几步看看,然后修改优 化参数。 我 遇 到 过 的 一 件 跟 结 构 优 化 关 系 很 大 的 算 例 是
CeO2 氧 空 位 结 构 电 子 局 域 的 问 题
(http://emuch.net/bbs/viewthread.php?tid=2954558)。按照一般方式(从优化好的 bulk
建 slab 模型,然后优化)得到一个 O 空位留下的两个电子均匀局域到 O 次外层三个 Ce 原子上, 得到空位形成能
2.34eV.经高人指点后, 调节空位附近 O 原子位置, 打破对称性后重新优化, 两个电子完美的局域到两个 Ce
原子上了。并且空位形成能降低到 2.0X eV。从这个例子可 以看到, 结构优化存在不少技巧的,
这些技巧建立在研究者对模拟对象的物理意义的理解上。 对物理图像的直观深入理解, 才能做好模型预设,
在此引导下才可能有目的的优化出不比寻 常的结果。
目前第一性原理理论中的交换关联泛函部分包含经验参数。 考虑这一点对优化结果的影响也 很有意思。比如有专家提到,DFT+U
参数对某些结构的收敛终态构型有影响。构型的变化 可能影响表面反应过程。基于这一点,一个好的计算研究可能就出来了。
真实过程总是复杂多变的。无论何种模拟,估计都可以找到一些试验现象来验证。但是到底 应该如何评判模拟结果,
如何从第一性原理研究中得出有意义的结论需要很好的洞察力。 这 样的模拟不见得就必须建立的试验的基础上,
完全凭空设计的模型有可能更能优美的解释本 质。
9、在单机上计算态密度好像不会出问题。我先谈一下我的看法:
第一个 WARNING,可以在 INCAR 文件中设置 NGX,NGY 和 NGZ 的值,设置的值要 足够大,就可以消除这个
warning。设置多大合适呢?这就要用到编译 vasp 时,同时也编译 得到的 make param 小程序,make param
可以帮助你预先检查你设置的文件是否正确,以及 某些参数的值是否合适。要得到合适的 NGX,NGY,NGZ 以及 NBANDS,先在
INCAR 中不设 置这些参数的值,然后运行 makeparam &param.inc,其中
param.inc 是包含了输出结果的文 件,在 param.inc 文件中你可以看到这些参数的值,以及计算大概需要多少的内存。然后把
param.inc 文件中的 NGX,NGY,NGZ 和 NBANDS 的值拷贝到 INCAR 文件中。
第二个是计算态密度时,我个人的做法是,一般把 KPOINTS 文件中的 k 点增多,然后 把 INCAR 文件中的
ISTART=1,ICHARG=11,当然还设置 RWIGS。最后把静止自洽计算得 到的 CHG 和 CHGCAR
文件拷贝到当前目录下。 从我在单机上的计算来看, 没有 WAVECAR 文件也是可以计算态密度的。我想你出现的这个问题,可能是你
cluster 上计算时,每个节 点上的 CHGCAR 和 WAVECAR 文件不一致造成的。 第三个是当 k 点数增加了,
会出现一个 WARING, 要把此 WARNING 消失掉, INCAR 在 文件中设置 NELMDL, 它的值小于等于默认值
(默认值好像是-5,你可以设为-6)。 没有 cluster
的系统用来计算,也没有这样的经历,我仅从在单机上的计算经验来谈,有错还请包涵
10、如何用 VASP 计算铁磁、反铁磁和顺磁
顺磁,意味进行 non-spin polarized 的计算,也就是 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(2x2)这种体系,原胞里面有 2 个 Mn
原子,那么你直 接让两个 Mn 原子的 MAGMOM 的绝对值一样,符号相反就可以了,再加上 ISPIN=2。这样
就可以实现进行反铁磁的计算了
11、vasp 在计算磁性的时候,oszicar 中得到的磁矩和 outcar 中得到各原子磁矩
之和不一致在投稿的是否曾碰到有审稿人质疑, 对于这个不一致你们一般是怎么 解释的了?
答:OSZICAR 中得到的磁矩是 OUTCAR 中最后一步得到的总磁矩是相等的。总磁矩和各 原子的磁矩(RMT
球内的磁矩)之和之差就是间隙区的磁矩。因为有间隙区存在,不一致是正 常的。
12、磁性计算应该比较负责。你应该还使用别的程序计算过磁性,与 vasp 结果比较是否一
致,对磁性计算采用的程序有什么推荐。 ps:由于曾使用 vasp 和 dmol
算过非周期体系磁性,结构对磁性影响非常大,因此使用这两 个程序计算的磁性要一致很麻烦。还不敢确定到底是哪个程序可能不可靠。
答:如果算磁性,全电子的结果更精确,我的一些计算结果显示磁性原子对在最近邻的位置 时,PAW 与 FPLAW
给出的能量差不一致,在长程时符合的很好。虽然并没有改变定性结 论。 感觉 PAW 似乎不能很好地描述较强耦合。 我试图在找出原因,
主要使用 exciting 和 vasp 做比较。计算磁性推荐使用 FP-LAPW, FP-LMTO, FPLO
很吸引人(不过是商业的),后者是 O(N)算法。
13、vasp 学习笔记 POTCAR 的建立
POTCAR 将要告诉 vasp 计算的系统中所包含的各种元素的赝势 pesudopotential,vasp 本身
就带有比较完善的赝势包, 我们需要做的就是选择我们需要具体哪种赝势, 然后把相应的文 件拷贝形成我们具体的 POTCAR 文件。我们以
GaAs 为例。 1)赝势的选择: vasp 的赝势文件放在目录 ~/vasp/potentials
下,可以看到该目录又包含五个子目录 pot pot_GGA potpaw potpaw_GGA potpaw_PBE
,其中每一个子目录对应一种赝势形式。 赝势按产生方法可以分为 PP (standard pesudopotential,其中大部分是
USPP, ultrasoft pesudopotential) 和 PAW (projector augmented wave
method)。按交换关联函数的不同又可以 有 LDA (local density approximation) 和 GGA
(generalized gradient approximation),其中 GGA 之下又可以再分为 PW91 和 PBE。
以上各个目录对应起来分别是 pot ==& PP, LDA ; pot_GGA
==& PP, GGA ; potpaw ==& PAW, LDA ;
potpaw_GGA ==& PAW, GGA, PW91 ; potpaw_PBE
==& PAW , GGA, PBE。选择某个目 录进去,
我们还会发现对应每种元素往往还会有多种赝势存在。 这是因为根据对截断能量的 选取不同还可以分为
Ga,Ga_s,Ga_h,或者根据半芯态的不同还可以分为 Ga,Ga_sv,Ga_pv 的
不同。 一般推荐选取 PAW_PBE。其中各个元素具体推荐哪种形式的赝势可以参考 vasp workshop 中有关赝势部分的
ppt。当然自己能测试之后在选择是最好不过的了,以后再聊。 2).POTCAR 的建立: 选好哪一种赝势之后, 进入对应的目录,
你会看到里边有这么几个文件, POTCAR.Z PSCTR.Z V_RHFIN.Z WS_FTP.LOG
。我们需要的是第一个。把它解压,如 zcat POTCAR.Z & Ga 。 对 As
元素我们也可以类似得到一个 As 文件。用 cp 命令或者 mv 命令把这两个文件都移 到我们的工作目录里。然后再用 cat
命令把这两个文件合并在一起,如 cat Ga As & POTCAR ,这样就得到了我们需要的
POTCAR。同理,有多个元素的 POTCAR 也可以这 样产生。这里需要注意的是,记住元素的排列顺序,以后在 POSCAR
里各个元素的排列就 是按着这里来的。 3).POTCAR 里的信息: 如果你想看 POTCAR 长什么样,可以用 vim POTCAR
命令,进去后可以用上下键移动光 标。想出来的时候,可以敲入:q!就可以。具体的 vim 的命令可以在网上查到。一般我会看 POTCAR
里的截断能量为多大,用 grep -in "enmax" POTCAR 。 据说 B3LYP 的赝势计算比较准,我在 MS
上面测试过,好像 DOS 和能带图的计算确实比 较准。不过不知道 vasp 有没有类似的赝势包。 hybrid functional
的计算,并不需要特定的 hybrid functional 的赝势。大部分就是基于 GGA-PBE
的赝势来做,也就是芯电子与价电子的交换关联作用,以及芯电子与芯电子的交 换关联作用还是基于 GGA-PBE
的,只是将价电子与价电子的交换关联作用通过 hybrid functional 交换关联来描述。
谢谢老师的解答。那具体操作是不是像网上写的那样,使用 GGA 的赝势,设置 GGA = B3, 然后更改 POTCAR 里面的
LEXCH =B3 就行了。我试过了,可以跑,不过结果没做详细的 分析。
14、VASP 中所有能量的物理意义及它们之间的区别,让你彻底 搞清楚 VASP 的所有能量
(一)首先我们应明白,固体的结合能就是固体的内能 E(结合)=U(内能) , 原因如下:
一般情况都把孤立原子的能量作为能量参考点。 前段时间有个同学问 VASP 中得出的绝 对能量是相对于什么的,其实就是相对孤立原子得。
(二)其次我们根据自由能与内能之间的关系 F=U-TS 而且我们都知道 VASP 的所有计算都是在绝对 0 度下的情况,T=0
代入上式,有 F=U。所以 结合就等于内能等于自由能。肯定有 Free energy TOTEN=energy without
entropy 恒成立... 这时候肯定有人会说不对啊,可以看 VASP 手册,候博的参考书作证,肯定不对得。
现在我告诉你确实它们二者确实有区别,区别在下面的情况 (1)当我们用 ISMEAR=-5
时,费米能这儿没有展宽,它算出来的就是完全在绝对 0 度的 能量。Free energy TOTEN=energy without
entropy 恒成立。 (2)有时为了在数学上处理的方便,为了更容易积分,我们也用 ISMEAR!=-5(!=是不
等于的意思)的方法,这个时候费米能这儿有一定的展宽。此时,我们容易想到,有展宽不 就是相当有一定的熵值吗?所以这个时候虽然算的是绝对
0 度的情况,但是有一定的熵值 (我们应明白, 这个熵值不是由一定的温度带来的, 而是数学处理的结果) 所以在 SMEAR! 。
=-5 的方法我们会发现 Free energy TOTEN 和 energy without entropy
有一定的差别。此时
energy without entropy 是 Free energy TOTEN 在 SIGMA 趋于 0 的极限。
注意:(1)有人在算单个原子的能量时会发现单个原子的能量虽然很小但并不是 0,但是按我
上面的推导,固体中的结合能是相对孤立体系的能量而来的,所以单个原子得到的 TOTEN 肯 定是 0 啊,原因在于我们的 POTCAR
不可能绝对合理,而且我们也知道计算单个原子的能量就 是为了检测赝势,单原子得到的 TOTEN 越小说明赝势越好。 但一般不会正好是
0.对这个说法 我还存在点疑问,写在了最后面。 (2)如果你注意的话,energy without entropy 与 Free
energy TOTEN 在 SIGMA 趋于 0 也不 是完全相等,但是也会发现它们之间的差别在 10E-3
左右,原因在于计算机求积分、求极限 不能像我们人一样达到任意的精度。
15、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 数量级上,考虑三个自由度,及其组合,就有很多很多可能了。请问该 怎么设置这样的过渡态虚频计算呢?
16、关于概念的问题做个讨论
(一)关于结合能。你说“结合能是定义为相距无穷远的原子结合形成一定结构的物质所放 出的能量”
你和我说的没区别,我说的是结合能是相对于“孤立原子做参考点的”,也就是它与周围任何
原子没有相互作用,和你所说的相距无穷远一回事,我这个好像没有任何错误。 ===============
这里你说的是没有错误,但是我觉得有必要先澄清一下。
(二)关于单点能。你说“它是第一性原理计算直接得到的能量,或者说是赝能,是一个空
间点阵平均每阵点上采用赝势计算所得到的能量,其中包含了结合能的贡献,但是更多的, 也包含了靠近芯区附近的电子在采用赝势近似下的能量,
这一部分能量既不是原子芯区附近 电子能量的真实反应,也不会影响化学键性质,不会对结合能有所贡献”。 我赞同你的大部分观点,
也提出你说的几点错误, 单点能准确的来说它包含了所有哈密顿的 量, 而且这儿的单点能不是你所说的“平均每个原子的能量”,
而是你计算的整个原胞的能量。 但是这个能量有一个参考点。你可以看候博得,也可以看我回的下一个贴子,至于“影响不
影响成键之类的内容”固体力学上已经说的很清楚了。 ============== 从你的回复中,
我可以知道你肯定没有学过晶体学或者空间群理论, 你应该看看晶体学国际
表中对于阵点的定义,阵点并不是每个原子,这里你的理解有问题,阵点是一个抽象点,一
个晶体中包含所有对称性的可以仅通过平移来构造整个晶体的结构所占据的位置就是一个
阵点,换句话说,一个阵点,就是一个满足平移对称性的原子集团,且该集团内部的位置满 足该晶体结构的全部对称性,而且它不仅仅是“原胞”
(三)你说“Free energy TOTEN 是体系总能,要减去阵点上分布的原子的能量再除以平均原
子数才是结合能(当然,这个和你的计算脚本的设计有关) ,而且这还没有考虑不加展宽时 没有被计算到的能带的因素”
我不赞同你后面说的几点。 Free energy TOTEN 从字面意思上我们也知道它的结果是自由能,
你可以说它是总能,因为根据我上面的推导,它们至少在数值是相等得。不在于你把它说成
什么,你就是把它说成总能,其实它还是等于结合能,等于自由能,等于内能。至于除不除
原子总数在于你想得到的是平均每个原子的还是总体系的,这在于个人。考虑不考虑展宽, 那要看 ISMEAR
等于几,做几个实例就会感觉到它考虑没考虑了。 =============== 这个能量确切的说应该是叫做考虑电子振动熵的体系总自由能,
当不考虑展宽的时候, 它是 等于总能的,如果你读过 Vasp 的代码,就知道 TOTEN 在 vasp 的计算中就是总能,这个和
结合能不是一个概念,还包含有非成键部分的贡献,至于内能的定义,如果你阅读过塞兹的 现代固体理论,或者 Pauling 的书,或者读过
Morse 当时提出 morse 势的那篇文献,就应该
知道,固体物理中所使用的内能,指的是离子实的动能和原子的“结合能”之和——这里结合 能之所以要打引号, 是因为按定义,
是要形成稳定结构或者亚稳态结构时才能称之为结合能, 内能定义并没有考虑粒子芯区附近电子能量的影响,正如你所说,是“相对于“孤立原子做参
考点的””,在 0K 下已经不考虑动能,因此就应是总能减去孤立原子的能量和才行,至于结 合能,则是稳定状态下结构的这个能量。
(四)你说“是否考虑展宽和结合能的定义没有关系” 我也没说和它的定义有什么关系啊, 但是由于数学处理带来的误差, 它对结合能的结果有一
定的影响啊。 (五)对单个原子的结合能的计算应该只计算 Gamma 点的能量,且用削除简并。 你说得第五点我不懂,
我算单个原子的能量时一直是按三个表面的构造方法来算的, 也没想 过什么简并。还望有高手给我帮助,第五点怎么理解。
============= 简单的操作是计算单个原子能量只考虑 Gamma 点,然后三边都设置在 10A 以上,且不相等
至于原因,应该去查量子力学的书,记得本科的时候,老师都会讲到的。 根据 F=U-TS,E(结合)=E(bulk)-nE(单个离散 。而
E(bulk)就是 U,通常我们取 E 单个离散)。 (结合) ( ) 单个离散 ( ) , 离散)为参考点。
(离散)为参考点。也就是把 E(离散)看为 0,这样推出来,在 T=0 时,F=E(结合) 它 (离散) ,这样推出来, (结合)
是包含你所说得那些,而且就像我在前面的贴子中说得, 得到的能量, 是包含你所说得那些,而且就像我在前面的贴子中说得,它就是总得
H 得到的能量,但是 说得 它有一个参考, 它有一个参考,而这个参考就是离散原子的能量
17、用 vasp 软件研究表面小分子的吸附解离遇到几个问题。
1 关于反应物和产物的结构,请问你们是如何构建的?(参考文献吗?)能
不能介绍一下相关的经验,个人认为好的开端是成功的一半,所以需要更多的经 验帮助。
答:.据说现在流行的过渡态算法多用 CINEB,vasp 自带的 NEB 也可以,但是镜像受
力优化和势垒精确度以及收敛速度略逊于 CINEB,lz 可以 google 到其官方网站了解。初始
和终了的结构是局部能量最低的构型(vasp 做弛豫收敛的构型) ,至于 lz 所关心的是从具体
结构之间的过渡的化,一是猜测,二是文献吧。基本出发点是能量趋向降低的构型。
2 关于结构的优化,能量收敛和力收敛的标准是否应该比默认值高,有利于 过渡态的搜寻? 答:DFT
理论计算能量可能更准确,对力的计算由于是对能量再求导,涉及数值算法
原因可能不准,所以收敛标准基本采用默认即可,有时适当提高(一个量级以内)也是可以
的。在合理的精度范围内讨论出合理的结果就达到目的了,这是师兄告我的。基本上收敛精 度高于默认值去跑的话,那就 god bless
3 关于结构虚频的处理,通过学习了解了消除虚频的方法--将对应原子的坐 标加到原坐标上。由于 vasp
的振动模式不能可视化,所以请问如何找到对应的 原子?另外关于虚频,大家是如何处理的呢。
答:.虚频没有关注过。只知道在鞍点位置才应该有一个虚频,属于鞍点具有的特性吧。
为何要消除呢?根据这个好像才能按速率理论算出反应率的吧, 这个我没有经验。 只关注过 势垒。
4 关于 images 的问题,最近利用脚本产生了 2、4 的 POSCAR 的文件,发 现只有 2 的能跑起来,我的服务器是
4 个核的。看别人发的贴:指定 4 个 cpu 同 时计算,应该是 1 个 cpu 一个 image 。为什么我设置 4 个
image 跑不起来呢?
答:image 数量越多,相互之间受力需要同时满足收敛准则,速度肯定会降低,但是这
样找到的能量最低路径也会与实际更加符合。没错, 每个 image 事实上是需要相当于单独计 算的构型的 cpu 的数目的,所以 4
个核的话,平均一个 image 一个 cpu,实在是不一定跑得 动。如果体系原子数较多,k
点网格较密,能量截断值再一高,估计经过很久才可能有结果。 5、在计算 dos 时,为什么会出现 WARNING:stress and
forces are not correct。这是那些 原因造成的,计算的是多铁物质。还有就是在优化时,TOTAL-FORCE
这项怎么全是 0 啊, 原子这个没有关系的 答:用优化后的 CONTCAR 算 band 或 dos
都会出现这个提示位置没有变化啊,这是怎 么回事啊? 6、我在做能带计算,发现自洽计算和非自洽计算中,OUTCAR 文件中都可以找到一个
费米能级的大小, 但是现在迷茫的是, 不知道在画能带时减去的费米能级的大小应该取哪个 结果文件中得值?是取自洽计算的结果文件中的
E-fermi ,还是取非自洽的结果文件中的 E-fermi 的大小?? 答: 算能带的那次计算, 指的是自洽 (求解 WAVECAR
和 CHGCAR 的那一步) 的那步计算, 还是非自洽(直接读取自洽的电荷和波函数,进行非自洽计算)的那次计算啊?我对这个说
法的标准不是很理解,希望您能再指导一下,我觉得画能带图时,费米能级的取值是否正确 是很关键的,取自洽计算得到的 E-fermi
热度(2)全文链接
(转)用VASP计算能量态密度(DOS)和能带
一、用VASP计算态密度
使用的软件:VASP
辅助分析计算的小程序:read_dos.f90
操作流程及主要步骤:
主要分两步:一、静态自洽计算;二、非自洽计算(以fcc结构的Al为例)
第一步静态自洽计算(得到自洽的电荷密度CHG、CHGCAR和E-fermi,提供给下一步非自洽计算用)
准备INCAR:即定义PREC,EDIFF,ENCUT,ISTART= 0(默认),IBRION =
2,ISMEAR(对于半导体和绝缘体:ISMEAR= 0,SIGMA = 0.05;ISMEAR= 0 表示采用Gaussian
smearing 方法; 对于金属:ISMEAR= 1,SIGMA = 0.2;ISMEAR= 1
表示采用1阶smearing方法)
SYSTEM=Al-fcc
PREC=Accurate
EDIFF= 1e-5
ENCUT= 250.0
ISTART= 0 ;IBRION = 2 #【默认NSW=0,此处定义IBRION=2无意义】
ISMEAR= 1 ; SIGMA = 0.2
准备好POSCAR文件,或以优化的晶格参数作为基础,把结构优化产生的CONTCAR拷贝成POSCAR
0.0 0.5 0.5
0.5 0.0 0.5
0.5 0.5 0.0
0.0 0.0 0.0
POTCAR直接从赝势库中得到
Automatic generation
MohkorstPack
0.0 0.0 0.0
记下OUTCAR中E-fermi的数值,下一步处理结果时会用到。
grep“E-fermi”OUTCAR
第二步非自洽计算
准备INCAR:即定义ENCUT,PREC,ISTART= 1(存在WAVECAR文件,所以取1),ICHARG =
11(表示从CHGCAR中读入电荷分布,并且在计算中保持不变),ISMEAR= -5 ( -5表示带Bloch修正的四面体方法;采用
tetrahedron method with Bl&ochl
corrections计算DOS更准确一些,在PWscf中也是如此) ,RWIGS (或LORBIT =
10,这时可不设RWIGS)
SYSTEM = Al-fcc
PREC = Accurate
ENCUT = 250 #【截断能不能大于上一步非自洽计算?最好与上一步计算保持一致】
ISTART = 1;ICHARG = 11
ISMEAR = -5
LORBIT = 10
# 【对于PAW势,可设置LORBIT = 10,此时可不用设置RWIGS参数;具体解释见附】
#【一般上述参数用来做非自洽计算、产生DOS已经足够;如NSW表示离子运动步数,默认为0(此时,IBRION默认为-1,表示
原子位置不移动),即做静态计算。自洽还是非自洽,如何判断?】
准备好KPOINTS文件,增加k点网格(增加k点可加大积分计算精度,因为是非自洽计算,k点数虽增大但仍可以很快)
Automatic generation
MohkorstPack
0.0 0.0 0.0
POSCAR、POTCAR同上一步自洽计算。
将上一步自洽计算得到的CHG、CHGCAR拷贝至同一目录下,提交运行。【ICHARG=11,ISMEAR=-5】
计算完后,得到包含了态密度值的DOSCAR文件,采用read_dos.f90小程序对态密度文件DOSCAR进行处理,得到总态密度dos_total.dat,及分波态密度*_dos.dat。
运行read_dos.f90时会提示:please input the
Ef,此时需要把第二步自洽计算得到的E-fermi的值给出。接着会提示:please input the number of
atoms in the POSCAR file,此时需要把POSCAR中的原子数给出。【操作./read_dos.x
输出文件的能量值是以费米能级作为能量参考零点。dos_total.dat的第一列数据是能量值,单位为eV;第二列数据是总态密度的值,单位State/eV.unit
cell;第三列数据是总态密度的积分值,也就是电子数,单位为electrons
(用origin等软件绘图时把dos_total.dat的最后一列删除即可) 。
*_dos.dat是相应原子【unit
cell中含有的每个原子都有各自对应的s、p、d分波态密度值】的分波态密度值,其中的第一列数据是能量值,单位为eV;第二、三、四列数据分别对应于s、p、d态的分波态密度值,单位为State/eV.atom。
#【如果将s、p、d态密度按每种元素求和,则态密度和的单位是State/eV.atom?如何将单位转换为State/eV.el.】
#【为什么自洽完了,要做非自洽计算的英文解释见后面附中的英文解释。提问:非自洽产生的DOSCAR
文件与上一步自洽计算产生的DOSCAR之间有什么区别?】
二、用VASP计算能带
使用的软件:VASP
辅助分析计算的小程序:read_band.f、get_kpoints.sh
操作流程及主要步骤:
计算材料的能带结构即色散曲线E(k),主要分两步:一、静态自洽计算(fcc);二、非自洽计算(以fcc结构的Al为例)
第一步静态自洽计算(同上)【因此,可单独建一个scf的文件夹,用以产生并存储电荷密度CHG、CHGCAR和E-fermi】
第二步非自洽计算(在固定电子密度的情况下,得到选取K点的能量本征值。)
准备INCAR:即定义ENCUT,PREC,ISTART= 1,ICHARG =
11(表示从CHGCAR中读入电荷分布,并且在计算中保持不变),并增加NBANDS(默认值为NELECT/2+NIONS/2,NELECT和NIONS分别为电子数和离子数,可以上一步自洽计算产生的OUTCAR文件中找到这两个参数,如grep
"NIONS" scf/OUTCAR,和 grep "NELECT"
scf/OUTCAR。)【NBANDS的设置,详见后面的英文注释】
SYSTEM = Al-fcc
PREC = Accurate
EDIFF = 1e-5
ENCUT = 250
ISTART = 1;ICHARG = 11
ISMEAR= 1;SIGMA =0.2
NBANDS = 12
准备KPOINTS文件,使用Line-mode,给出高对称性k点之间的分割点数。【分割越密,则路径积分越准确,VASP
5.x格式有变吗?】
k-points along high symmetry lines !注释行,无特别的意义
20 ! intersections,沿G-X特殊点之间产生10个k点
Line-mode ! 程序自动产生特殊k点间的k点
Rec ! 各k点相对于倒格子基失来写的 Al的元胞的高对称性点
0.5 0.0 0.5 ! X
0.0 0.0 0.0 ! G
0.0 0.0 0.0 ! G
0.5 0.5 0.5 ! L
0.5 0.5 0.5 ! L
0.5 0.25 0.75 ! W
0.5 0.25 0.75 ! W
0.375 0.375 0.75 ! K
0.375 0.375 0.75 ! K
0.0 0.0 0.0 ! G
#【对高对称性点的理解及注释见后附注部分。】
说明:通过指定Line-mode,VASP会自动在起点和终点之间插入指定的K点数,比如上面的文件就是指定VASP计算沿着X到Gamma点再到L、W、K回到Gamma点,每个方向上各取20个K点。
POSCAR、POTCAR同上一步自洽计算。
将上一步自洽计算得到的CHG、CHGCAR拷贝至同一目录下,提交运行,计算完后可以从OUTCAR文件或者EIGENVAL文件里得到需要的每个K点的能级信息【与运用WAVECAR得到的能带结构一样,CHGCAR是WAVECAR的积分。是的,波函数代表一切】
运行小程序get_kpoints.sh,得到kpoints.dat(即高对称点的位置;)。编译程序read_band.f,运行,系统提示:
please input thefermi energy:E-fermi from scf OUTCAR
将上一步得到的E-fermi值输入到屏幕,则输出bands.dat,并产生spoint.dat文件
(该文件给出了作图时高对称性点的横坐标位置),因此,可用origin等绘图软件作图了。【read_band.f程序的解释见附】
lpf的经验和心得:
对于第一步自洽计算:由于在能带计算时k点是一些在倒空间高对称线上的点,不能进行自洽计算【解释得好啊】,因此在进行能带计算前必须加一次自洽计算以得到精确的电荷密度值。
自洽计算得到的电荷密度文件CHGCAR是能带和DOS计算需要的输入文件。
对于非自洽计算:计算能带时ICHARG= 11,金属用ISMEAR=1;半导体或绝缘体,用ISMEAR=0 。
计算态密度:ICHARG = 11,ISMEAR=-5
【WAVECAR文件在计算DOS和BAND过程中有用么?是万能的,但要开计算band和dos的源程序调用谁了】
用origin绘图技巧,用Vi打开spoint.dat文件,得到作图时高对称点的横坐标位置,将横坐标端点分别设为第一个、最后一个高对称点对应的坐标;画竖线,双击,修改Objectproperties---Cordinates:Units选为Scale,设置开始x和终点x均为spoint.dat文件中的高对称点数值;y值即设为纵坐标范围。【实际上图形中的横坐标范围已经定死了,无法修改;其实也恰好与spoint.dat中点的范围一致。如何断定spoint.dat文件中数值具体对应哪个高对称点符号?】
赝势文件夹各个目录对应起来分别是pot ==& PP, LDA ; pot_GGA
==& PP, GGA ; potpaw ==& PAW, LDA
potpaw_GGA ==& PAW, GGA, PW91 ; potpaw_PBE
==&PAW , GGA, PBE。
选择某个目录进去,我们还会发现对应每种元素往往还会有多种赝势存在。这是因为根据对截断能量的选取不同还可以分为Ga,Ga_s,Ga_h,或者根据半芯态的不同还可以分为Ga,Ga_sv,Ga_pv的不同。
2:non-scf
ISTART=1;ICHARG=11时,PREC必须与做静态计算的设置一致(尤其是截断能),否则会报错
WARNING:dimensions on CHGCAR file are different
ERROR: chargedensity could not be read from file CHGCAR for
POSCAR中的那个缩放系数对于正格矢起作用如果原子坐标是用Car格式,它也起作用。如果原子坐标是Direct格式,则它不起做作用。【yexq评:缩放系数即对应PWscf中的alat,只影响晶格基矢量的大小,个人感觉不会影响原子的crystal位置】。
NSW是优化原子位置和晶胞个开启项,默认=0。如果没有设置NSW,即使ISIF=3,依然只做静态计算。
5:KPOINTS
采用Line-mode时,”0 0 0 !gama“的最后一个0和!之间需要有空格。【yexq评:真是经验总结,很仔细】
附录-yexq总结:
1. 为什么自洽计算完了,还要做非自洽计算,请见下面的英文解释
Calculating a DOS can be done in two ways:
The simple one is toperform a static (NSW=0, IBRION=-1)
selfconsistent calculationand to take the DOSCAR file from this
calculation.
However, the simple approach discussed above is not applicable
in all cases: A high quality DOS requires usually very fine
kmeshes. You should think at least in orders of 16x16x16 meshes for
small cells and even for large cells you might need something like
6x6x6- or 8x8x8-meshes.
The usual way,to do DOS or band-structure calculations in this
case is the following: the charge density and the effective
potential converge rapidly with increasing number of k-points.
So, as a first step one generates a high quality charge density
using a few k-points in a static selfconsistent run. The next step
is to perform a non-selfconsistent calculation using the CHGCAR
file from this selfconsistent run (i.e. ICHARG is set to 11, see
section 6.14) .
Mind, this is the only way to calculate the band structure,
because for band-structure calculationsthe supplied k-points form
usually no regular three-dimensional grid and therefore a
selfconsistent calculation gives pure nonsense !
【但是还是需要前面的自洽计算】
For ICHARG=11, all k-points can be treated independently, there
is no coupling between them, because the charge density and the
potential are kept fixed. Therefore there is also no need to treat
the k-points within one single run simultaneously.
【可新建1个pdos文件夹,非自洽计算可安排在该文件夹中进行】
2. NBANDS的设置问题
可见NELECT是所采用赝势中的价电子数吗?例如 对TiH2:取14/2 + 3/2 = 8.5,PdH:
取17/2+2/2=9.5
3. 对高对称性点选取的理解
采用与自洽计算所用的相同大小的胞(元胞或单胞),来选取高对称性点;
高对称性点并不一定形成闭合的曲线;
应对结构施加正确的symmetry后,选高对称性点。这也说明,能带计算不是在uniform的格子上进行,而是沿高对称性点划定的路径积分。所以,如果结构没施加对称性,则选取的高对称性点将与施加了对称性操作的“相同”结构不一样;同理按元胞找的对称性点与按单胞找的对称性点将不一样。
对于具有相同对称性的结构,其高对称性点都有一致的约定的取法,即与具体物质无关。例如同为fcc结构的不同物质,其高对称性点将是一样的。因此可以通过下述网站输入结构的空间群号寻找高对称性点:
http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-kv-list,或采用Material
studio 得到高对称性点。
4.能带计算完成后,绘制能带图时,如何确定图中高对称性点的位置。
spoint.dat文件中数值按大小排列,其对应的高对称性点的顺序可能与计算前KPOINTS中输入的高对称点的顺序一致,实际上高对称性点的crystal坐标(即KPOINTS中的三个数值)可在EIGENVAL文件中找到(或OUTCAR中找到)。
其实,read_band.f程序很简单,即将路径上的各k点crystal坐标求均方根,然后逐点累加,并记路每一k点对应的engenvalue;当遇到特殊k点(高对称性点)时,同时将该点坐标输出到spoint.dat文件中。
该程序跟据KPOINTS文件中特殊k点之间的取样间隔来判断某点是否是特殊k点,如是则同时输出其坐标到到spoint.dat文件中。因此,当特殊k点之间插入的一般k点数变化时,则应修改源文件中的if(mod((i-1),20).eq.0
.or.i.eq.nks语句,例如将20修改为所插入的一般k点数或实际分割的k点数。因此,该程序缺乏通用性,然而程序虽简单,但却高度体现了作者对能带结构的理解。
实际高对称性点的位置可以直接求均方根确定吗?(不能,没有加和,还应该加上路径上其他点坐标产生的加和)
2.部分参数的英文注释
(1)ISMEAR = -5
SIGMA = width of the smearing in eV
ISMEAR = 1
SIGMA = 0.2
ISMEAR determines how the partial occupancies fnk are set for
each wavefunction. For the finite temperature LDA SIGMA determines
the width of the smearing in eV.
-5 tetrahedron method with Bl&ochl
corrections
For the calculation of the total energy in bulk materials we
recommend the tetrahedron method with Bl&ochl
corrections (ISMEAR=-5). This method also gives a good account for
the electronic density of states (DOS).
The only drawback is that the methods is not variational with
respect to the partial occupancies. Therefore the calculated forces
and the stress tensor can be wrong by up to 5 to 10 % for
For the calculation of phonon frequencies based on forces we
recommend the method of Methfessel-Paxton
(ISMEAR&0). For semiconductors and insulators the
forces are correct, because partial occupancies do not vary and are
zero or one.
For the calculations of the DOS and very accurate total energy
calculations (no relaxation in metals) use the tetrahedron method
(ISMEAR=-5).
For semiconductors or insulators use the tetrahedron method
(ISMEAR=-5), if the cell is too large (or if you use only a single
or two k-points) use ISMEAR=0 in combination with a small
SIGMA=0.05.
For relaxations in metals always use ISMEAR=1 or ISMEAR=2 and an
appropriated SIGMA value (the entropy term should be less than 1
meV per atom). Mind: Avoid to use ISMEAR&0 for
semiconductors and insulators, since it might cause problems. For
metals a sensible value is usually SIGMA= 0.2 (which is the
(2)RWIGS or LORBIT
If RWIGS or LORBIT (Wigner Seitz radii, see section 6.326.33) is
set in the INCAR file, a lm- and site-projected DOS is calculated
and also written to the file DOSCAR. One set of data is written for
each ion, each set of data holds NDOS lines with the following
energy s-DOS p-DOS d-DOS
energy s-DOS(up) s-DOS(down) p-DOS(up) p-DOS(dwn) d-DOS(up)
d-DOS(dwn)
for the non spin-polarized and spin polarized case respectively.
As before the written densities are understood as the difference of
the integrated DOS between two pins.
(3) DOSCAR file
Mind: For relaxations, the DOSCAR is usually useless. If you
want to get an accurate DOS for the final configuration, first copy
CONTCAR to POSCAR and continue with one static (ISTART=1; NSW=0)
calculation.【静态自洽计算?】
(4)LORBIT
The default for LORBIT is .FALSE. (respectively 0).
This flag determines, together with an appropriate RWIGS (see
section 6.32), whether the PROCAR or PROOUT files (see section
5.21) are written.
The file PROCAR contains the spd- and site projected wave
function character of each band.
The wave function character is calculated, either by projecting
the wavefunctions onto spherical harmonics that are non zero within
spheres of a radius RWIGS around each ion (LORIT=1, 2), or using a
quick projection scheme relying that works only for the PAW method
(LORBIT=10,11,12, see below). If the LORBIT flag is not equal zero,
the site and l-projected density of states is also calculated.
If the projector augmented wave method is used, LORBIT can also
be set to 10, 11 or 12. This alternative setting selects a quick
method for the determination of the spd- and site projected wave
function character and does not require the specification of a
Wigner-Seitz radius in the INCAR file (the RWIGS line is neglected
in this case).The method works only for PAW POTCAR files and not
for ultrasoft or norm conserving pseudopotentials. 即对每个能带的波函数进行s
pd和site分解或投影;LORBIT = 10,不设RWIGS,针对PAW势的计算。
LORBIT = 10表示对波函数按原子site进行s 、p、d轨道分解或投影;并输出投影波函数到PROCAR or
PROOUT 文件中,详见上面英文注释部分:This alternative setting selects a quick
method for the determination of the spd- and site projected wave
function character。什么样的quick method?
d-band center相关的帖子
http://emuch.net/html/9537.html
http://emuch.net/html/4578.html
/forum.php?mod=viewthread&tid=164973
Nrskov等提出的 d-能带中心模型,即金属表面原子与吸附质相互作用随其 d-能带中心渊着d冤负移而减弱, 且吸附能与
着d成近似线性关系。
[转载]VASP电荷密度差
[http://emuch.net/html/9228.html
一般分为:差分电荷密度图(deformation charge density),二次差分电荷密度图(difference
charge density)
两者的区别有很多种说法,其中一种定义认为:
差分电荷定义为成键后的电荷密度与对应的点的原子电荷密度之差。通过差分电荷密度的计算和分析,可以清楚地得到在成键和成键电子耦合过程中的电荷移动以及成键极化方向等性质。
“二次”是指同一个体系化学成分或者几何构型改变之后电荷的重新分布。
如何做对电荷密度做差,目前有两个版本的公式:
个人认为公式(1)和差分电荷密度的定义比较吻合,但是这种方法会产生很多干扰信息,而文献中常用的是公式(2)。
如何用VASP做电荷密度差图?
① 对优化后的结构做静态自洽计算
要注意的是,AB、A和B要分别放在同样大小的空间格子中并保证A和B与AB中相应坐标不变,计算时也要保证三次自洽计算所采用的FFT
mesh 一致(NGXF,NGYF,NGZF)。
INCAR中几个注意的参数:
IBRION = -1;NSW = 0;
NGXF,NGYF,NGZF
② 做差
将第一步得到的三个CHGCAR文件依次命名为CHGCAR1、CHGCAR2、CHGCAR3,采用脚本程序按照公式②对三个文件做差,得到CHGdiff文件。
③ 视图
3D显示可以将得到的CHGdiff文件添加vasp后缀(CHGdiff.vasp),用VESTA读取,另外VESTA也可以做2D切面,缺点是没有标度。
2D的显示也可以用lev00去做切面,然后用origin作图,得到的2D面电荷密度差带有标度便于分析(具体做法可参阅lev00说明书)。
④ 以CO为例
图1左图是CO的电荷密度差黄色部分表示电荷密度增加,蓝色表示电荷密度减少,右图是Gaussian模拟的CO的HOMO和LUMO轨道,可以看出黄色部分对应HOMO,蓝色部分对应LUMO。
图2 是用lev00切面,origin做图得到的CO二维图示。
VASP表面计算
表面计算对精度和技巧要求相当高,所以在进行此类计算时需要非常细心,计算之前应仔细阅读和消化第8章中的基本知识。在下面几节中将介绍表面计算的典型步骤,需要注意的是,即使你严格按照给定步骤进行计算,依然可能会出现错误,得到没有任何物理意义的结果,这就需要你对可能出现的错误(结合第8章)进行仔细的分析,比如,FFT网格大小、K点数目是否足够,计算结果是否收敛,原子位置是否正确等,更一般的来说就是INCAR、POSCAR、KPOINTS、param.inc等输入文件是否正确。还需注意,前一步的计算误差可能导致后面计算中出现更大的错误,比如晶格常数1%的计算误差可能导致表面弛豫结果中3%的误差。因此,在初始计算中(包括体性质计算,确定FFT网格大小和K点的合适数目等)多花点时间是非常值得的。
1. 体性质计算
材料体性质计算的第一步,可以使用四面体方法(ISMEAR=-5或-4)、通过提高K点数目直到计算收敛到预期精度。这里所说的预期精度没有特定标准,一般说来如果表面能计算误差要控制在10
meV以内,可以提高K点数目直到总能收敛到1
meV左右即可。计算材料体性质的slab模型通常含有20-100个原子,因此为得到可靠的表面能结果需要进行非常精确的体能量计算,实际上越精确越好,推荐选择PREC=High。
第二步,将四面体方法改为有限温度方法(ISMEAR=0或N,N为正整数),有两种选择:
Gaussian方法(ISMEAR=0),较小的SIGMA值(SIGMA=0.1)。前面计算中多次用到这种方法,但更推荐第二种方法:
Methfessel-Paxton方法(ISMEAR=1),SIGMA值要尽量大,但要保证OUTCAR文件中的自由能和总能之差(比如entropy
T*S项)可忽略不计。entropy项可以很好的估计可能误差的大小,必须在预期误差范围内(如前例中是1 meV左右)
从现在开始,将无须改变ISMEAR、SIGMA和ENCUT值,而通过逐步提高K点数目来重复体性质计算,此时K点的收敛速率应该和四面体方法基本差不多。
选择一个用在表面计算中合适的K点数目和截断能,计算平衡晶格常数,注意尽量避免包络错误(Wrap around
errors,设置PREC=High)。现在得出的晶格常数就是将来表面计算中使用的晶格常数,而自由能是所有后续计算的参考值,此外还要计算entropy值(OUTCAR中的entropy*TS项)或者记下总能和physical能(σ-&0)。
2. 表面计算的FFT网格和K点数目
第一步需要确定合理的FFT网格大小,为从根本上杜绝包络错误,可以选择VASP的推荐参数(或者设置PREC=High)。开始尝试计算时,可选择5个原子层和5个真空层的slab模型,取一个合理的不太大的K点数目,在OUTCAR文件中将给出能严格避免包络错误的FFT网格大小值:
WARNING: wrap around error must be expected
Set NGX to 22
采用这个推荐网格会导致计算时间的延长,但至少应该可以做一次这样精确的计算。如果想减小计算时间,可以尝试使用3/4规则(PREC=Med),将其计算结果与精确收敛结果进行一下对比。
下一步需要确定合适的K点网格数目。体性质计算已经提供了一个参考,在表面计算中,通常含有一个长的晶格矢量和两个短的晶格矢量,在长晶格方向上,一个K点网格一般就足够了,因为在这个方向上色散是由真空层造成的。在短晶格方向上,网格划分数的收敛速度和体材料基本相同,提高K点数目知道自由能足够收敛即可。要再次提醒,避免包络错误!
可能用到的互相检查方法:
(1) 每个原子的entropy值应该和体性质计算相同,否则降低SIGMA值重新所有计算。
(2) 总的力偏移量应该足够小,否则原因是FFT网格不够大,应当相应地提高。
(3) 检查力随着K点网格和FFT网格的收敛速度。
3. 原子层和真空层数目
从现在开始,固定K点网格数目和包络错误大小(比如,保持严格杜绝包络错误的FFT网格大小和实际的FFT网格大小的比例为一定值),测试需要多少原子层和真空层,才能得到合理的表面能值,以及表面第一层(或第二层)原子收敛的作用力。
注意:在一次计算和下一次计算中修改的参数不要超过一个,比较k点数目和超胞大小都不同的两次计算结果几乎是不可能的。特别要注意FFT网格:如果增大了超胞尺寸而没有相应地提高FFT网格的尺寸,计算结果不但没有改进反而更糟糕,因为包络错误增大了
本文引用地址: /blog-479.html
● pro
热度(1)全文链接转载自:N.
[转载]VASP计算结构弛豫设置及收敛判据
1)收敛判据的选择
结构弛豫的判据一般有两种选择:能量和力。这两者是相关的,理想情况下,能量收敛到基态,力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大,力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的,能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。
到底是以能量为收敛判据,还是以力为收敛判据呢?关心能量的人,觉得以能量为判据就够了;关心力相关量的人,没有选择,只能用力作为收敛标准。对于超胞体系的结构优化,文献大部分采用Gamma点做单点优化。这个时候即使采用力为判据(EDIFFG=-0.02),在做静态自洽计算能量的时候,会发现,原本已经收敛得好好的力在不少敏感位置还是超过了结构优化时设置的标准。这个时候,是不是该怀疑对超胞仅做Gamma点结构优化的合理性呢?是不是要提高K点密度再做结构优化呢
在我看来,这取决于所研究的问题的复杂程度。我们的计算从原胞开始,到超胞,到掺杂结构,到吸附结构,到反应和解离。每一步都在增加复杂程度。结构优化终点与初始结构是有关的,如果遇到对初始结构敏感的优化,那就头疼了。而且,还要注意到,催化反应不仅与原子本身及其化学环境有关,还会与几何构型有关。气固催化反应过程是电子的传递过程,也是分子拆分与重新组合的过程。如果优化终点的构型不同,可能会导致化学反应的途径上的差异。仅从这一点来看,第一性原理计算的复杂性,结果上的合理性判断都不是手册上写的那么简单。
对于涉及构型敏感性的结构优化过程,我觉得,以力作为收敛判据更合适。而且需要在Gamma点优化的基础上再提高K点密度继续优化,直到静态自洽计算时力也是达到收敛标准的。
(2)结构优化参数设置
结构优化,或者叫弛豫,是后续计算的基础。其收敛性受两个主要因素影响:初始结构的合理性和弛豫参数的设置。
初始结构包括原子堆积方式,和自旋、磁性、电荷、偶极等具有明确物理意义的模型相关参数。比如掺杂,表面吸附,空位等结构,初始原子的距离,角度等的设置需要有一定的经验积累。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步(默认的),后面就会越来越快了。
总的说来,一般入门者,多看手册,多想多理解,多上机实践总结,比较容易提高到一个熟练操作工的水平。
如果要想做到“精确打击”,做到能在问题始发的时候就立刻采取有效措施来解决,就需要回归基础理论和计算方法上来了。
(3)优化结果对初始结构和“优化路径”的依赖
原子吸附问题不大,但是小分子吸附,存在初始构型上的差异。slab上水平放置,还是垂直放置,可能导致收敛结果上的差异。根据H-K理论,理想情况下,优化得到的应该是全局最小,但在数值计算的时候可能经常碰到不是全局最小的情况。实际操作中发现,多个不同初始结构优化收敛后在能量和结构上存在一定差异。
为了加快收敛速度,特别是对于表面-分子吸附结构,初始放松约束,比如EDIFF=1E-3,EDIFFG=-0.3,NSW=30可能是很好的设置。但是下面的情况应当慎重:
EDIFF=1E-3;
EDIFFG=-0.1;!或者更小
NSW=500;!或者更大
电子步收敛约束较小,而离子步约束偏大,离子步数又很多,这种情况下,可能导致的结果是结构弛豫到严重未知的区间。
再在这个基础上提高约束来优化,可能就是徒劳的了——结果不可逆转的偏向不正常的区间。
好的做法,是对初始结构做比较松弛的约束,弛豫离子步NSW应该限制在一个较小的数值内。EDIFF=1E-3的话,EDIFFG也最好是偏大一些,如-0.3而不是-0.1.
这样可以在较少的步数内达到初步收敛。
对于远离基态的初始结构,一开始在非常松弛的约束下跑若干离子步,时间上带来的好处是很大的。对于100个原子的体系用vasp做Gamma点优化,如果一开始就是正常优化(
EDIFF=1E-4,EDIFFG=-0.02)设置,开始十个离子步可能都要花上几个小时。如果这个时候才发现输入文件有错误,那下午的时间就白费了,顺便带上晚上机器空转。
所以,我习惯的做法,是在初始几步优化后,会用xcrysden
检查一下XDATCAR中的数据,用xdat2xyz.pl生成movie.xyz,然后看看弛豫过程是不是按照设想的那样。后续过程跑完一个收敛过程,就再检查一下movie.xyz。如此这般,才放心的展开后续计算。
本文来自: 小木虫论坛
http://emuch.net/bbs/viewthread.php?tid=2512497&fpage=1
● pro
热度(1)全文链接转载自:N.
电荷密度图、能带结构、态密度的分析
能带图的横坐标是在模型对称性基础上取的K点。为什么要取K点呢?因为晶体的周期性使得薛定谔方程的解也具有了周期性。按照对称性取K点,可以保证以最小的计算量获得最全的能量特征解。能带图横坐标是K点,其实就是倒格空间中的几何点。纵坐标是能量。那么能带图应该就是表示了研究体系中,各个具有对称性位置的点的能量。我们所得到的体系总能量,应该就是整个体系各个点能量的加和。
主要是从以下三个方面进行定性/定量的讨论:
1、电荷密度图(charge density);
2、能带结构(Energy Band Structure);
3、态密度(Density of States,简称DOS)。
电荷密度图是以图的形式出现在文章中,非常直观,因此对于一般的入门级研究人员来讲不会有任何的疑问。唯一需要注意的就是这种分析的种种衍生形式,比如差分电荷密图(def-ormation
charge density)和二次差分图(difference charge
density)等等,加自旋极化的工作还可能有自旋极化电荷密度图(spin-polarized charge
density)。所谓“差分”是指原子组成体系(团簇)之后电荷的重新分布,“二次”是指同一个体系化学成分或者几何构型改变之后电荷的重新分布,因此通过这种差分图可以很直观地看出体系中个原子的成键情况。通过电荷聚集(accumulation)/损失(depletion)的具体空间分布,看成键的极性强弱;通过某格点附近的电荷分布形状判断成键的轨道(这个主要是对d轨道的分析,对于s或者p轨道的形状分析我还没有见过)。分析总电荷密度图的方法类似,不过相对而言,这种图所携带的信息量较小。
成键前后电荷转移的电荷密度差。此时电荷密度差定义为:delta_RHO = RHO_sc - RHO_atom
其中 RHO_sc 为自洽的面电荷密度,而 RHO_atom
为相应的非自洽的面电荷密度,是由理想的原子周围电荷分布堆彻得到的,即为原子电荷密度的叠加(a superposition of
atomic charge densities)。需要特别注意的,应保持前后两次计算(自洽和非自洽)中的 FFT-mesh
一致。因为,只有维数一样,我们才能对两个RHO作相应的矩阵相减。
能带结构分析现在在各个领域的第一原理计算工作中用得非常普遍了。首先当然可以看出这个体系是金属、半导体还是绝缘体。对于本征半导体,还可以看出是直接能隙还是间接能隙:如果导带的最低点和价带的最高点在同一个k点处,则为直接能隙,否则为间接能隙。
1)因为目前的计算大多采用超单胞(supercell)的形式,在一个单胞里有几十个原子以及上百个电子,所以得到的能带图往往在远低于费米能级处非常平坦,也非常密集。原则上讲,这个区域的能带并不具备多大的解说/阅读价值。因此,不要被这种现象吓住,一般的工作中,我们主要关心的还是费米能级附近的能带形状。
能带的宽窄在能带的分析中占据很重要的位置。能带越宽,也即在能带图中的起伏越大,说明处于这个带中的电子有效质量越小、非局域(non-local)的程度越大、组成这条能带的原子轨道扩展性越强。如果形状近似于抛物线形状,一般而言会被冠以类sp带(sp-like
band)之名(此陈述有待考证—博主加)。反之,一条比较窄的能带表明对应于这条能带的本征态主要是由局域于某个格点的原子轨道组成,这条带上的电子局域性非常强,有效质量相对较大。
3)如果体系为掺杂的非本征半导体,注意与本征半导体的能带结构图进行对比,一般而言在能隙处会出现一条新的、比较窄的能带。这就是通常所谓的杂质态(doping
state),或者按照掺杂半导体的类型称为受主态或者施主态。
4) 关于自旋极化的能带,一般是画出两幅图:majority spin和minority
spin。经典的说,分别代表自旋向上和自旋向下的轨道所组成的能带结构。注意它们在费米能级处的差异。如果费米能级与majority
spin的能带图相交而处于minority spin的能隙中,则此体系具有明显的自旋极化现象,而该体系也可称之为半金属(half
metal)。如果majority spin与费米能级相交的能带主要由杂质原子轨道组成,可以此为出发点讨论杂质的磁性特征。
5)做界面问题时,衬底材料的能带图显得非常重要,各高对称点之间有可能出现不同的情况。具体地说,在某两点之间,费米能级与能带相交;而在另外的k的区间上,费米能级正好处在导带和价带之间。这样,衬底材料就呈现出各项异性:对于前者,呈现金属性,而对于后者,呈现绝缘性。因此,有的工作是通过某种材料的能带图而选择不同的面作为生长面。具体的分析应该结合试验结果给出。
原则上讲,态密度可以作为能带结构的一个可视化结果。很多分析和能带的分析结果可以一一对应,很多术语也和能带分析相通。但是因为它更直观,因此在结果讨论中用得比能带分析更广泛一些。简要总结分析要点如下:
在整个能量区间之内分布较为平均、没有局域尖峰的DOS,对应的是类sp带(此陈述有待考证—博主加),表明电子的非局域化性质很强。相反,对于一般的过渡金属而言,d轨道的DOS一般是一个很大的尖峰,说明d电子相对比较局域,相应的能带也比较窄。
2)从DOS图也可分析能隙特性:若费米能级处于DOS值为零的区间中,说明该体系是半导体或绝缘体;若有分波DOS跨过费米能级,则该体系是金属。此外,可以画出分波(PDOS)和局域(LDOS)两种态密度,更加细致的研究在各点处的分波成键情况。
3)从DOS图中还可引入“赝能隙”(pseudogap)的概念。也即在费米能级两侧分别有两个尖峰。而两个尖峰之间的DOS并不为零。赝能隙直接反映了该体系成键的共价性的强弱:越宽,说明共价性越强。如果分析的是局域态密度(LDOS),那么赝能隙反映的则是相邻两个原子成键的强弱:赝能隙越宽,说明两个原子成键越强。上述分析的理论基础可从紧束缚理论出发得到解释:实际上,可以认为赝能隙的宽度直接和Hamiltonian矩阵的非对角元相关,彼此间成单调递增的函数关系。
4) 对于自旋极化的体系,与能带分析类似,也应该将majority spin和minority
spin分别画出,若费米能级与majority的DOS相交而处于minority的DOS的能隙之中,可以说明该体系的自旋极化。
5)考虑LDOS,如果相邻原子的LDOS在同一个能量上同时出现了尖峰,则我们将其称之为杂化峰(hybridized
peak),这个概念直观地向我们展示了相邻原子之间的作用强弱。
由于金属的能带有可能穿越fermi能级,从而引起总能计算时的不连续变化。为了避免这种情况,需要引入分数的占据态smearing。
● pro
热度(1)全文链接转载自:N.
倒空间及Bloch定理的基本讨论
对于凝聚态方向的研究人员而言,由晶体平移不变性所引入的倒空间是非常基本的概念。但是倒空间的物理意义、扩展布里渊区和简约布里渊区的关系,以及由此而
来的能带等概念很难理解,这在一定程度上也限制了研究人员对于自己研究方向的深入理解与分析。
1、倒空间、倒格矢的引入
晶体最显著的特征就是其点阵的平移不变性:可看成是最小结构单元在三维无限扩展所形成的无限大体系。这种结构特点直接导致关于晶体的一个结论:无法唯一的确定体系的原点。假设两个人分别站在某两个相隔为一个实空间基矢Rm的格点上,那么这两个人的观察不会有丝毫的差别。因此,对于晶体而言,其体系的波函数应该满足周期性边界条件:
ψ(r)=ψ(r+{N}R)
其中R是三维实空间的基矢,{N}是一组(三个)整数。很显然,平面波解满足这个条件,而且为了满足周期性边界条件,k{N}R=2nπ
(*)。考虑最简单的情况,取n=1,那么对应于特定的R,可以唯一的确定一个矢量K。接下来在三个方向上分别改变n,那么由上式可以得到k的一个空间分布,注意到这个分布中的所有点都是K的整数倍,因此可以将K看作是另一个空间的基矢。这个空间我们称为倒空间。相应的,这个空间同样也存在着平移不变性。
在这里我们应当指出,如果k和K乘
上普朗克常数,都是电子动量的量纲,因此在倒空间进行的这种平移看起来在物理上并不等价,这也是让很多人困惑的一个地方。一般的固体理论教材上往往直接将
数学表达式给出,但对这个问题没有进行太多的解释。想澄清这个问题,必须和上文所说的“无法唯一的确定晶体体系的原点”联系在一起来理解。事实上,在倒空
间上远离原点的平移,等于将观察点在实空间里向原点移近。简而言之,倒空间中的平移,等于同时改变k和{N}以保证(*)式中右侧的n不变。从这个意义上说,K同样可以定义一个最}

我要回帖

更多关于 微信超出验证频率限制 的文章

更多推荐

版权声明:文章内容来源于网络,版权归原作者所有,如有侵权请点击这里与我们联系,我们将及时删除。

点击添加站长微信