跳转至

VASP 使用

约 3662 个字 199 行代码 预计阅读时间 15 分钟

介绍

  • VASP 全称：Vienna Ab-initio Simulation Package

多粒子体系的复杂性主要体现在交换关联能 \(E_{xc}\) 项上，对交换关联能 \(E_{xc}\) 的精确描述是求解 KS 方程的关键所在。

引入关于电子密度函数的近似泛函，如 LDA、GGA、杂化泛函等

LDA 忽略了非均匀效应，认为在具有同等电子密度的前提下，空间任意点的 \(E_{xc}\) 相同与均匀电子云的相同

GGA 对 LDA 中忽略的电子密度的非均匀效应进行了修正，即加入电子密度梯度的作用

Perdew-Burke-Ernzerhof (PBE) 形式的 generalized gradient approximation (GGA) 泛函（表达电子间的交换关联作用）

全电子势描述电子 - 原子核相互作用，如 EMTO（The Exact Muffin-Tin Orbitals）

投影缀加平面波赝势（PAW）方法（描述离子 - 电子相互作用）

PAW (Projected Augmented Wave) 投影缀加波，是基于密度泛函理论（DFT）开发的描述电子、原子核行为的全电子方法

APW 是增广平面波方法（Augmented Plane Wave），也是一种全电子方法，将电子分为软、硬两部分，前者用 PW 描述，后者用 LCAO 描述（Linear Combination of Atomic Orbitals，原子轨道的线性组合）

PBE 是泛函，描述电子的交换 - 相关能的拟合函数，以 Perdew-Burke-Ernzerhof 三位开发者的名字缩写命名

PW 是指平面波基组，用于展开波函数或者说原子、分子轨道

PAW 是独立于 PBE 的理论方法，但是我们常常会见到 POTCAR 中称为 PAW-PBE 赝势，把三个概念放在一起了，意思实际上是针对不同的的泛函利用 PAW 方法相应调参优化得到的一致性赝势文件

PAW (Projected Augmented Wave) 全电子理论计算方法

含各种类型计算的通用 INCAR 文件： GitHub - WMD-group/INCAR: A generic INCAR file for the density functional theory package VASP

VASP 计算问题小结 - chempeng

计算一致性：K 点密度、截断能、SMEAR 平面波数量也应一致（对于空位形成能计算）

DFT-D3：vdW 相互作用修正

Heyd–Scuseria–Ernzerhof 泛函 (HSE06)：更精确，处理电子和光学性质

METAGGA

SCAN (Strongly constrained and appropriately normed)

METAGGA - VASP Wiki

参考资料

VASP 中计算电子基态的算法

Algorithms used in VASP to calculate the electronic groundstate - Vaspwiki

k 点积分

K-point integration - Vaspwiki

phonon dispersion 计算

Computing the phonon dispersion - Vaspwiki

原子计算

Calculation of atoms - Vaspwiki

分子动力学计算

Molecular dynamics calculations - Vaspwiki

Molecular dynamics - Tutorial - Vaspwiki

GW 计算

Practical guide to GW calculations - Vaspwiki

GW approximation

内含 VASP 计算相关的案例

GitHub - hello-arun/tutorials: Tutorials of codes such as VASP, Quantum Espresso and Lammps

DFT磁性的计算(本例为用VASP计算FCC Ni的磁矩)

VASP教学 - 计算材料学

使用

工具

  • VASPMO：用于显示 VASP 计算的波函（或分子轨道）。它能够读取 VASP 的输出文件 PROCAR 和 CONTCAR，并产生 Gaussian 输出格式的输出文件，用于其它显示工具，如 Molekel、Chemcraft、Gabedit、Molden 和 JMol 等）读取，进而绘制和观看体系的分子轨道

算例

晶体高对称点 - 知乎

示例

# 能带结构、DOS 计算
ISTART   = 1
ICHARG   = 11

# 其余计算
ISTART   = 0
ICHARG   = 2

# 断点后续算；其他参数值不改变
ISTART   = 1
ICHARG   = 1

静态计算

# INCAR 重要参数
IBRION   = -1
NSW      = 0
ISIF     = 2

孤立原子能量

# INCAR 重要参数
ISPIN    = 2              # 开启自旋

KSPACING = 1.0            # 可确保 K 点密度为 1x1x1
KGAMMA   = .TRUE.

ISMEAR   = 0              # 使用 Gaussian smearing
SIGMA    = 0.001          # 展宽值取小；或 0.05
NELM     = 300            # 增加电子步

ISYM     = 0              # 关闭对称性

内聚能

  • 内聚能（cohesive energy/atomization energy）定义：孤立原子结合成固体/晶体时释放的能量，为正值（一般以这个为准）；部分定义为其逆过程，为负值

  • 分子动力学中部分经验势函数不考虑自能，孤立原子能量为 0，体系平均原子能量即为内聚能

  • DFT 中孤立原子能量不为 0

收敛性测试

弛豫/结构优化计算

  • 建议结构优化（尤其变胞优化）后再进行一次静态计算的原因（能量更准确）：Hartree 势积分格点的问题；因为积分格点是从最开始定死的，变胞优化完之后沿用的是优化前的格点结果自然会有点问题（来自 GPUMD 学习交流群）

  • 高精度结构优化：ISIF=3 优化 --> ISIF=2 优化（基本 1-2 步结束）--> scf 单点能计算

  • INCAR 参数示例：

Global Parameters
SYSTEM   = Relaxation
ISTART   = 0
ICHARG   = 2
ISPIN    = 1
LCHARG   = .FALSE.
LWAVE    = .FALSE.
PREC     = Accurate
LREAL    = .FALSE.
ENCUT    = 500

KSPACING = 0.15
KGAMMA   = .TRUE.

Electronic Relaxation
ALGO     = Normal
ISMEAR   = 1
SIGMA    = 0.05
EDIFF    = 1E-06
NELM     = 300
NELMIN   = 6

Ionic Relxation
NSW      = 150
IBRION   = 2
ISIF     = 3
EDIFFG   = -1E-02

NPAR     = 4

（平衡）晶格常数

  • EOS 拟合（扫描法）：在晶格常数附近取 N 个点，分别进行单点能计算，进行 EOS 拟合

  • 直接弛豫

用力弛豫时，若存在内变量的体系，如固溶体，溶质原子几乎不会待在在平衡位置（有可能），近邻的溶剂原子偏离其理想位置，因此需要对它们的位置进行弛豫，用扫描的方法，不做弛豫，得不到准确的位置，因此需要用到弛豫的方法；对于计算胞对称性不太好，用扫描法，一般适合只有一个变量变化，变量多，需要用弛豫的方法（hcp 结构）

表面能

点缺陷形成能

  • 空位、间隙原子形成能

  • 超胞尽可能是立方体形状

  • 要做原子数 n（不能太大）和空位形成能的收敛性测试

态密度、能带计算

  • 参考：

  • 计算流程：

    • 弛豫计算（或结构优化；初始构型很好，可忽略此步）
    • 静态自洽计算（可不用拷贝自洽计算生成的 WAVECAR）
    • 态密度计算：拷贝自洽计算生成的 CHGCAR，非自洽计算（ICHARG=11，增加 K 点密度；K*a=45 可满足要求）
    • 能带计算：拷贝自洽计算生成的 CHGCAR，非自洽计算（ICHARG=11，K-path）

  • 自洽与非自洽计算的区别：电子密度是否匹配；不是静态与弛豫计算的区别！

  • 开启自旋极化，DOS 会有上下两条线（自旋向上、向下；上下对称、不对称的含义是什么）

  • 建议 NEDOS 数值稍微取密一些（NEDOS=2000/3000 足够好）

  • DOS、能带计算耗时相对较少

  • DOS 计算相比弛豫对截断能没有那么敏感，截断能值可以设小一些（？）

  • 态密度相关输出文件：DOSCAR、PROCAR

ISMEAR   = -5             # DOS计算

ISMEAR   = 0              # 能带计算

Bader 电荷计算

# INCAR 参数设置
IBRION   = -1
NSW      = 0
LCHARG   = .TRUE.
LAECHG   = .TRUE.         # Bader 计算


# 后处理
chgsum.pl AECCAR0 AECCAR2     # 得到 CHGCAR_sum，包含了 VASP 中定义的原子总电荷密度
bader CHGCAR -ref CHGCAR_sum  # 得到 BCF.dat、ACF.dat、AVF.dat
cat ACF.dat                   # 查看 CHARGE 列

# BCF.dat
# TODO: 待确认
# 电荷极大值序号和坐标；体积内的bader电荷积分；距离极大值最近的原子和距离

# ACF.dat；包含了所有原子的价电子数信息
# TODO: 待确认
# 原子序号和坐标；价电荷数；到零通量面最小距离；原子体积）
# 示例
#         X           Y           Z       CHARGE      MIN DIST  ATOMIC VOL
--------------------------------------------------------------------------------
1    3.164973    3.164973   3.164973    8.129202     2.022066   50.790052
2    0.000000    0.000000   0.000000   13.223566     1.494570   30.667423
3    0.000000    0.000000   3.164973    7.549077     1.582486   57.390231
4    3.164973    0.000000   0.000000    7.549079     1.582486   57.390910
5    0.000000    3.164973   0.000000    7.549079     1.582486   57.390910
--------------------------------------------------------------------------------
VACUUM CHARGE:              0.0000
VACUUM VOLUME:              0.0000
 NUMBER OF ELECTRONS:        44.0000

ELF

  • 参考：VASP视频教程-电荷局域分析_哔哩哔哩_bilibili

  • 电子局域化函数 (electron localization function, ELF)

  • 在自洽计算中添加/修改以下参数；计算结束后，得到 ELFCAR 文件；使用 VESTA 进行可视化；可将 ELF 图与构型视图叠放在一起，对照效果更好

PREC     = Accurate
LELF     = .TRUE.

弹性常数计算

# INCAR 参数设置
ENCUT    = 520            # 尽量高一些；个人设置 700

IBRION   = 6
NSW      = 1
NFREE    = 2        # 4 或 2
POTIM    = 0.015
ISIF     = 3


# 单位换算
1 kBar = 0.1 GPa


# 获取弹性常数数据
# 方式 1
grep -A9 'TOTAL ELASTIC MODULI' OUTCAR
# 方式 2
vaspkit - 2 - 203


# OSZICAR 相关内容
Found     1 degrees of freedom:
Adding 6 more for cell shape distortion

# log/out 文件相关内容
Found     1 degrees of freedom:
Adding 6 more for cell shape distortion
Finite differences POTIM= 0.01500 DOF=   7

Finite differences progress:
 Degree of freedom:   1/  7
 Displacement:        2/  2
 Total:               2/ 14

# OUTCAR 相关内容
Found     1 degrees of freedom:


        Directions for atom     1:
        --------------------------

              1.0000000000000000      0.0000000000000000      0.0000000000000000
        Strain: 6 additional degrees of freedom
        ---------------------------------------
 Finite differences:
   Step               POTIM =   1.500000000000000E-002
   Degrees of freedom DOF   =            7


Finite differences progress:
 Degree of freedom:   1/  7
 Displacement:        2/  2
 Total:               2/ 14


# 相关 Warning（有无影响，影响是否大）
VERY BAD NEWS! internal error in subroutine IBZKPT:
Reciprocal lattice and k-lattice belong to different class of lattices. Often results are still useful...      96

功函数计算

  • VASP:功函数

  • 步骤：结构弛豫；静态计算（修改以下 INCAR 参数）

NSW      = 0
LVHAR    = .TRUE.

AIMD 计算

PSTRESS                     # Parinello-Rahman 压浴控制目标压强

# LANGEVIN_GAMMA、LANGEVIN_GAMMA_L 这两个额外 tag 必须设置
# PMASS 可选 tag
LANGEVIN_GAMMA              # 原子自由度的 langevin 摩擦系数
LANGEVIN_GAMMA_L            # 点阵自由度的 langevin 摩擦系数
PMASS                       # 点阵自由度的虚拟质量
  • NVT 参数设置
# INCAR 参数设置
ENCUT    = 400            # 不用很高
EDIFF    = 1E-5           # 不用很小；EDIFFG 不用于 AIMD
ISMEAR   = 0              # 不用设置成 1
SIGMA    = 0.05
LWAVE    = .FALSE.
LCHARG   = .FALSE.
LREAL    = Auto
PREC     = Normal
ALGO     = Fast

IBRION   = 0              # 开启 AIMD
MDALGO   = 2              # MDALGO、ISIF 控制系综和热浴
ISIF     = 2
SMASS    = 0              # 控制 AIMD 运行过程中的速度
ISYM     = 0              # 关闭对称性

POTIM   = 1               # 时间步长（单位 fs）
NSW     = 10000           # 跑多少步；10 ps

TEBEG   = 300             # 起始温度
TEEND   = 300             # 终止温度

NWRITE  = 1               # 长时 AIMD，建议取 1/0


# 数据获取
# 压强
grep 'external pressure' OUTCAR | awk '{print $4}'

# 运行步数，温度
grep 'T= ' OSZICAR | awk '{print $1 " " $3 " " }'

进行 NVT 模拟时，不会设置 PSTRESS，所以 Pullay stress=0

原子数越多，温度波动越小

外压 external pressure 数值 = Total in kB XX、YY、ZZ 分量之和的平均值 - Pullay stress

基于分子动力学模拟，可以通过对速度自关联函数（velocity autocorrelation function，VACF）进行傅里叶变换得到材料的振动态密度（vibrational density of states， VDOS）。VACF 是根据动力学模拟出来的轨迹文件和速度文件，求算系统在某一时刻的速度与另一时刻速度的关联程度的函数，直接看 VACF 并不能很直观的得到一些信息，而 VDOS 直接对应实验红外光谱，可以直观的对高温或高压下的振动变化情况等进行分析。

AIMD结合vaspkit计算振动态密度

其他

  • VASP 拉伸模拟
  • SOC：自旋轨道耦合、旋轨耦合
    • VASP INCAR 参数 LSORBIT = T 开启 SOC 时，需使用 vasp_ncl，使用 vasp_std 会出现以下报错
ERROR: non collinear calculations require that VASP is compiled
 without the flag -DNGXhalf and -DNGZhalf
ERROR: non collinear calculations require that VASP is compiled
 without the flag -DNGXhalf and -DNGZhalf

# pymatgen 解析 SOC 的 Vasprun 文件出现以下 Warning（实际不影响）
pymatgen/io/vasp/outputs.py:161: UserWarning: Float overflow (*******) encountered in vasprun
warnings.warn("Float overflow (*******) encountered in vasprun")