VASP 简明教程
密度泛函理论 (DFT) 主要流程
结构优化 -> 静态自洽 (SCF) -> 非自洽 (NSCF)-> ? (探索中...)
非自洽 (NSCF)：包含电子态密度、能带结构等。
INCAR 参数
| SYSTEM |
体系名称（不重要） |
|
计算方式相关参数 |
| ISTART |
运行方式，一般就选0，有WAVECAR的情况，默认是1 |
| ICHARG |
计算初始电子密度的方式。0：从初始轨道推导计算出电荷密度。1：读取已有的CHARGCAR文件。2：利用叠加原理直接拿原子电荷密度重叠作为电荷密度。11：电荷密度从始至终都是采用已有的CHARGCAR文件进行计算，不迭代更新，这是计算DOS和BAND时的常用办法。一般选2 |
| SYMPREC |
POSCAR文件中的原子位置的精度，一般选1e-4 |
| electronic optimization |
电子优化相关参数 |
| ENCUT |
截断能量，单位eV，ENCUT越大，计算精度越高，但计算量会越大。VASP可以直接从POTCAR中得到每个元素默认的截断能，并且取最大值作为整个计算ENCUT的默认值 |
| PREC |
选精度，一般选Accurate |
| ALGO |
计算电子优化的算法，一般选Normal，此外有Very_Fast和Fast |
| NELM |
自洽步数的最大值，默认是60，也可以选到300 |
| EDIFF |
结束自洽的条件，一般选1e-6 |
| ISMEAR |
smearing使用的方法。0：Gaussian smearing。-1：Fermi smearing。-2：从文件中读取。一般选0或1 |
| SIGMA |
ISMEAR对应的参数。有两种常用组合 0：0.05；1：0.20 |
| ionic relaxation |
离子弛豫相关参数 |
| EDIFFG |
结束离子优化的条件，一般选1e-4 |
| NSW |
最大离子弛豫步数，不做离子弛豫就设置为0，做离子弛豫一般设置为500 |
| ISIF |
确定是否计算应力张量以及允许在松弛和分子动力学运行中改变哪些主自由度,IBRION=0时选0，其他时候选2，这样会计算力和压力 |
| IBRION |
决定离子更新和移动的方式。-1：不移动。0：用于分子动力学。1：离子弛豫（RMM方法）。2：离子弛豫（CG算法）。3：离子弛豫（阻尼分子动力学算法）算离子弛豫一般选2。5，6：算HESS矩阵和声子频率（有限差分算法）。7，8：算HESS矩阵和声子频率（微扰理论） |
| POTIM |
设置分子动力学中的时间步长或离子弛豫中的步长。对于IBRION = 0，POTIM给出了所有从头算分子动力学运行中的时间步长（以 fs 为单位），因此必须提供它，否则 VASP 在启动后会立即崩溃，一般选不知道。对于IBRION =1、2 和 3，一般选0.10，这分别对应于使用准牛顿算法、共轭梯度算法和阻尼分子动力学的离子弛豫，POTIM标记用作步宽的缩放常数。准牛顿算法对这个参数的选择特别敏感。对于IBRION = 5 和 6，使用有限差分方法进行声子计算，其中POTIM是每个离子的位移宽度以计算 Hessian 矩阵。 |
| performance optimization |
加快计算速度 |
| KPAR |
例子里没有设置过这个参数，make 4 groups, each group working on one set of k-points |
| NCORE |
例子里没有设置过这个参数，one orbital handled by 4 cores |
| LREAL |
例子里没有设置过这个参数，real space projection; slightly less accurate but faster |
| k网格相关 |
k网格相关 |
| KSPACING = 0.2 |
K网格取点间距 |
| PSTRESS |
目标压力，单位是KBar=0.1GPa |
| 读文件相关 |
读文件相关 |
| LCHARG |
输出CHG和CHGCAR，一般是F，SCF时需要设置为T |
| LWAVE |
输出WAVECAR，一般是F，SCF时需要设置为T |
| NPAR |
并行用，一般选4 |
对INCAR参数解释的一些网页：https://zhuanlan.zhihu.com/p/559467814
结构优化 (弛豫Relaxation)
准备工作
结构文件POSCAR
INCAR
| Text Only |
|---|
| PREC = A # 精度，一般选Accurate
ISTART = 0 # 有WAVECAR的情况,默认是1,否则设0(0代表⼀个全新的计算)
ICHARG = 2 # ISTART=0时,默认是2,否则设0,计算DOS和BAND时设为11,代表从已有的CHARGCAR进行计算
ENCUT = 850 # Cutoff eV,一般550左右
EDIFF = 1e-5 # 结束自洽的条件,一般选1e-6,越小精度越高
EDIFFG = -0.001 # 结束离子优化的条件,一般选1e-4
IBRION = 2 # 默认-1(固定离子),决定离子更新和移动的方式,对于初始结构接近局域极小时是最好的选择IBRION=2:比较可靠,弛豫困难时,推荐使用IBRION=3:对非常糟糕的初始结构进行弛豫时可考虑使用
NSW = 500 # 最大离子弛豫步数,不做离子弛豫(即不做结构优化)就设置为0,做离子弛豫一般设置为500
ISIF = 2 # IBRION=0时选0, 否则选2
ISMEAR = 0 # 一般选0或1
SIGMA = 0.05 # ISMEAR对应的参数,有两种常用组合 0:0.05; 1:0.20
# Determine whether the orbitals (file WAVECAR), the charge densities (file CHGCAR and CHG) are written.
LCHARG = F # write WAVECAR
LWAVE = F # write CHGCAR
#NPAR = 4 # 并行用,一般选4,How many cores work on one orbital
#ISPIN = 2 # 是否考虑自旋极化,默认1(不考虑,non spin polarized),反之填2
ALGO = F # 计算电子优化的算法,一般选Normal,此外有Very_Fast和Fast
IVDW = 11 # vdW corrections, 默认0(不校正), 11(zero damping DFT-D3 method of Grimme)
|
KPOINTS
可用vaspkit生成。
| Text Only |
|---|
| KPT-Resolved Value to Generate K-Mesh: 0.030
0
Monkhorst-Pack
8 8 1
0.0 0.0 0.0
|
The description of each line is given as follows:
- Header (comment).
- Specifies the k mesh generation type. : automatic generation scheme.
- Gamma-centered (or Monkhorst-Pack) grid.
- Number of subdivisions in each direction.
- Optional shift of the mesh.
关于K点密度的选取：
第一是对称性，不同的晶体结构有不同的对称性，我们应该根据晶体结构选择合适的K点，这个在各种教材中已经给出了。不过手动选择K点只在能带计算的时候是必要的，在做优化和性质计算时我们通常选择自动撒点。
一般来说，K点的选取跟晶胞的边长有关，边长越长，所需的K点就越少。如果是块体的单胞，K点的分割要取密一些，比如我们的TiO2单胞，K点分割12 12 12就非常精确了。如果我们将这个晶胞延xyz方向各扩展一倍，这时只要6 6 6足够。越大的晶胞，所需的K点数就越小。如果是表面，那么在某个方向有大于10Å的真空层，那么这个方向上K点最多只须取2足矣。
三个方向K点数和晶格常数乘积近似一致，但一般乘积不要大于50，尤其是针对比较大的体系，否则会使得优化速度比较慢。
真空层一定不要加K点，用1即可，因为这个方向没有结构，加K点没有意义，比如二维体系在xy方向是周期性体系，z方向存在真空层，那么这个方向一定不要加K点，用1即可。
POTCAR
从VASP赝势文件中获得，如果结构单元包含多种原子，则需要将对应原子种类的POTCAR文件顺序连接。
- Data that was required for generating the pseudo potentials.
- Number of valence electrons.
- Atomic mass.
- Energy cut-off…
- If the cell contains different atomic species, the corresponding POTCAR files have to be concatenated, in the same order as the atomic species are given in the POSCAR file.
自洽 (SCF)
静态自洽计算的准备工作：
使用结构优化输出的CONTCAR拷贝成新的POSCAR。
静态自洽主要涉及的参数：
| Text Only |
|---|
| NSW=0 # 不再进行离子弛豫（固定）
LWAVE=TRUE # 需要在WAVECAR中输出波函数
LCHARGE=TRUE # 需要在CHG和CHGCAR中输出电荷密度
|
非自洽 (NSCF)
能带计算 (Band)
非自洽用于求解态密度(DOS)，能带(Band) 电子结构分析 或者光学等其他性质。
能带计算主要涉及的参数：
| Text Only |
|---|
| ISTART=1 # 接着计算
ICHAGE=11 # 读入自洽的CHGCAR，从而开始能带或者态密度的非自洽计算
NSW=0 # 不再进行离子弛豫
LWAVE=F # 不再需要输出波函数
LCHARGE=F # 不再需要输出电荷密度
# 如果需要元素投影的话就加上这个
LORBIT = 11 # or 10 for per atom
|
关于能带计算时的 KPOINTS：
https://mattermodeling.stackexchange.com/questions/4335/k-points-value-in-kpoints-file-for-the-vasp-band-calculation
在算正方形结构的时候，只需要 GAMMA->X->S->GAMMA 就行，长方形的时候 GAMMA->X->S->Y->GAMMA ，具体请参照布里渊区，下为C13F8结构的KPOINTS示例：
| Text Only |
|---|
| Special k-points for band structure
20 ! intersections # 一般撒20个点，少于20基本不太行
line-mode
reciprocal
0.0000000000 0.0000000000 0.0000000000 GAMMA
0.5000000000 0.0000000000 0.0000000000 X
0.5000000000 0.0000000000 0.0000000000 X
0.5000000000 0.5000000000 0.0000000000 S
0.5000000000 0.5000000000 0.0000000000 S
0.0000000000 0.0000000000 0.0000000000 GAMMA
|
可使用 Python pymatgen 库直接绘制 Band 图：
| Python |
|---|
| import matplotlib.pyplot as plt
from pymatgen.io.vasp.outputs import Vasprun
from pymatgen.electronic_structure.plotter import BSPlotter
# parse_projected_eigen=True 提取投影信息
bs_vasprun=Vasprun("./vasprun.xml", parse_projected_eigen=True)
bs_data=bs_vasprun.get_band_structure(line_mode=1)
plotter = BSPlotter(bs=bs_data)
plotter.get_plot()
plt.savefig("element_band.png")
|
DOS涉及的参数(基于上述参数后继续修改)：
| Text Only |
|---|
| # 带有Blöchl校正的四面体方法，使用以伽马为中心的k-mesh。计算块体材料中的总能量。很好地解释了电子DOS。
ISMEAR = -5
# 注释掉SIGMA
# SIGMA = 0.05
# write DOSCAR and lm decomposed PROCAR file
LORBIT = 11 # or 10 for per atom
# 加大态密度计算的撒点数，这样可以让态密度的峰值更平顺
NEDOS = 2000
|
在能带及态密度均计算完成后，可使用 Python pymatgen 库直接绘制 Band-DOS 图：
将 BSDOSPlotte 参数 bs_projection 与 dos_projection 的值均设置为 elements ，表示读取能带和态密度的投影到元素的信息：
| Python |
|---|
| # parse_projected_eigen=True 提取投影信息
dos_vasprun = Vasprun("./dos/vasprun.xml", parse_projected_eigen=True)
dos_data = dos_vasprun.complete_dos
bs_vasprun = Vasprun("./band/vasprun.xml", parse_projected_eigen=True)
bs_data = bs_vasprun.get_band_structure(line_mode=True)
plt_1 = BSDOSPlotter(bs_projection='elements', dos_projection='elements')
|
将参数 bs_projection 设置为 element ，dos_projection 的值设置为 orbitals，表示将能带投影到元素，态密度投影到轨道：
| Python |
|---|
| # parse_projected_eigen=True 提取投影信息
dos_vasprun = Vasprun("./dos/vasprun.xml", parse_projected_eigen=True)
dos_data = dos_vasprun.complete_dos
bs_vasprun = Vasprun("./band/vasprun.xml", parse_projected_eigen=True)
bs_data = bs_vasprun.get_band_structure(line_mode=True)
plt_1 = BSDOSPlotter(bs_projection='elements', dos_projection='orbitals')
|
以C13F8结构为例：
| Python |
|---|
| import matplotlib.pyplot as plt
from pymatgen.io.vasp.outputs import Vasprun
from pymatgen.electronic_structure.plotter import BSDOSPlotter,BSPlotter,BSPlotterProjected,DosPlotter
# parse_projected_eigen=True 提取投影信息
dos_vasprun = Vasprun("./dos/vasprun.xml", parse_projected_eigen=True)
dos_data = dos_vasprun.complete_dos
bs_vasprun = Vasprun("./band/vasprun.xml", parse_projected_eigen=True)
bs_data = bs_vasprun.get_band_structure(line_mode=True)
plt_1 = BSDOSPlotter(bs_projection="elements", dos_projection="elements", fig_size=(16,12), font="DejaVu Sans")
plt_1.get_plot(bs=bs_data, dos=dos_data)
plt.savefig('plt.png')
|
态密度计算 (DOS)
态密度计算所需要的文件：
关于DOS计算的KPOINTS：
撒点密度取值尽量大一点，一般保证 \(晶格常数*kpoint \ge 30\) 即可，但最好不要超过50，会影响计算效率。
| Text Only |
|---|
| KPT-Resolved Value to Generate K-Mesh: 0.030
0
Monkhorst-Pack
20 20 1
0.0 0.0 0.0
|
计算完DOS后，可使用 Python pymatgen 库直接绘制 DOS 图：
| Python |
|---|
| from pymatgen.io.vasp import Vasprun
from pymatgen.electronic_structure.plotter import DosPlotter
from pymatgen.electronic_structure.core import OrbitalType
import matplotlib.pyplot as plt
op = input("1.按元素投影 DOS\n2.按轨道投影 DOS\n please input option:\n")
if op == "1":
v = Vasprun('./vasprun.xml')
cdos = v.complete_dos
dos = cdos.get_element_dos()
plotter = DosPlotter()
plotter.add_dos('Total DOS', v.tdos)
plotter.add_dos_dict(dos)
# 按需更改轴区间
plotter.get_plot(xlim=[-5, 5], ylim=[-30, 30])
plt.savefig('element_dos.png')
else:
result = Vasprun('./vasprun.xml', parse_potcar_file=False)
complete_dos = result.complete_dos
# 需要按自己的结构更改
pdos_Ce = complete_dos.get_element_spd_dos('Ce')
pdos_Ni = complete_dos.get_element_spd_dos('Ni')
pdos_Mn = complete_dos.get_element_spd_dos('Mn')
# 需要按自己的结构更改
plotter = DosPlotter()
plotter.add_dos('Total DOS', result.tdos)
plotter.add_dos('Ce(f)', pdos_Ce[OrbitalType.f])
plotter.add_dos('Ni(d)', pdos_Ni[OrbitalType.d])
plotter.add_dos('Mn(d)', pdos_Mn[OrbitalType.d])
# 按需更改轴区间
plotter.get_plot(xlim=(-5, 2), ylim=(-150, 350))
plt.savefig('element_dos_spd.png')
|
以C13F8结构为例：
晶体结构力学性能相关计算
VASP5.0及更新版本在做完结构优化之后，可以直接计算弹性常数EC，INCAR示例如下：
| Text Only |
|---|
| System = EC
PREC = A
ISTART = 0
ICHARG = 2
ENCUT = 850
EDIFF = 1e-6
EDIFFG = -0.001
#NSW = 500
ISMEAR = 0
SIGMA = 0.05
LCHARG = F
LWAVE = F
#NPAR = 4
#ISPIN = 2
#ALGO = F
#IVDW = 11
#For EC
IBRION = 6
ISIF = 3
NFREE = 4
|
关于EC计算的KPOINTS：
撒点密度取值尽量大一点，一般保证 \(晶格常数*kpoint \approx 50\) 即可。
| Text Only |
|---|
| KPT-Resolved Value to Generate K-Mesh: 0.030
0
Monkhorst-Pack
10 10 1
0.0 0.0 0.0
|
计算完成后，查看OUTCAR文件：
| OUTCAR |
|---|
| ELASTIC MODULI (kBar)
Direction XX YY ZZ XY YZ ZX
--------------------------------------------------------------------------------
XX 625.7835 36.2769 92.3649 0.0000 0.0000 0.0000
YY 36.2797 625.7786 92.3701 0.0000 0.0000 0.0000
ZZ 92.3806 92.3808 1098.5611 0.0000 0.0000 0.0000
XY -0.0234 0.0335 -0.0043 61.3585 0.0000 0.0000
YZ -0.0063 -0.0055 -0.0358 0.0000 246.1519 0.0000
ZX -0.0022 0.0000 -0.0011 0.0000 0.0000 247.0661
--------------------------------------------------------------------------------
|
晶体结构稳定性判据
与 Materials Project 的结果进行比较，需注意单位转换 kBar -> GPa：
对于 Orthorhombic (斜方晶系) 结构，需要满足以下机械稳定性条件：
\[
\begin{cases} {{C_{1 1} > 0, \, C_{1 1} C_{2 2} > C_{1 2}^{\ 2},}} \\ {{C_{4 4} > 0, \, C_{5 5} > 0, \, C_{6 6} > 0,}} \\ {{C_{1 1} C_{2 2} C_{3 3}+\, 2 C_{1 2} C_{1 3} C_{2 3}-\, C_{1 1} C_{2 3}^{\ 2}}} \\ {{-\, \, C_{2 2} C_{1 3}^{\ 2}-\, C_{3 3} C_{1 2}^{\ 2} > 0}} \end{cases}
\]
弹性常数->Elastic moduli
在Voigt 近似条件下：
BV = [(C11 + C22 +C33) + 2(C12 + C13 + C23)]/9
GV = [(C11 + C22 +C33) − (C12 + C13 + C23) + 3(C44 + C55 + C66)]/15
而对于立方体系：
C11 = C22 = C33, C12 = C13 = C23, C44 = C55 = C66
因此：
BV = [C11 + 2C12]/3,
GV = [C11 − 2C12 + 3C44]/5
后台运行 VASP
只针对于无任务调度系统单机运行VASP的情况，如果服务器或集群已经安装如 Slurm 等的任务调度工具，请忽略！
集群提交任务方式请参考：集群教程
| Bash |
|---|
| # 首先写好当前主机运行 vasp 所需要的步骤至 shell(vasp.sh)
# 后台执行vasp.sh文件，将标准日志输出到output.log文件中，将错误日志也输出到output.log文件中
# 命令末尾的 & 代表后台运行
# 2>&1是一个整体，将错误内容重定向到标准输出中
nohup ./vasp.sh > output.log 2>&1 &
|
常见错误处理
SIGSEGV, segmentation fault occurred
只需将堆栈大小(stack size)调大，一般可以设置成256M或者没有限制。
| Bash |
|---|
| # 查看当前限制
ulimit -s
# 256M
ulimit -s 262140
# or unlimited
ulimit -s unlimited
|
Inconsistent Bravais lattice types found for crystalline and reciprocal lattice
是什么原因导致的不清楚...
解决办法：
| Text Only |
|---|
| # INCAR
SYMPREC = 1E-4
# or
SYMPREC = 1E-6
|
BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES (EXIT CODE:9)
这种报错一般是爆内存了，我是在能带计算的时候出现的，可能的原因：KPOINTS文件中选的高对称点路径太多了，撒点密度过高。我虚拟机给了20G内存依然不能正常计算，那就没办法了...
一些解决办法：
删掉后面的一些高对称点，调低KPOINTS第二行数值，会导致计算比较粗略。
http://bbs.keinsci.com/thread-45731-1-1.html
KPOINTS第二行的20表示每条能带上撒20个点，一般不需要这么多，调小这个数就可以避免爆内存的问题。虽然撒点数减少之后计算比较粗略，但可以在 ~/.vaspkit 配置文件里可以设置能带外推，另外，得到粗略的能带图后可以找出自己感兴趣的少数几条能带，用较高的分辨率重算。
OrginPro 有时候导入数据后，无法对数据列进行编辑
点击表格左上方的黄色小文件夹图标，点击“解锁导入的数据”即可。