计算材料学
0%

VASP教学

发表于 分类于
本文字数： 6.8k 阅读时长 ≈ 6 分钟

2020年10月末北京科技大学学生

大纲

  • 安装软件，登陆集群，熟悉常见的linux操作，熟悉vasp的提交格式，掌握slurm任务调度系统，成功提交计算并分析；

    HW1:提交任务计算Mg，Cu晶格常数并与实验和其他计算值做比较

  • 纯金属和第二相表面能的计算；

    HW2:分别计算Mg，Cu的(001)的slab能量及表面能，刘阳调节好MS，MS如何切表面。

  • 第二相表面吸附H原子，电子逸出功（功函数）；

    HW3:

    (1)计算Al最稳定表面的H原子吸附能（要考虑到不同的代表性吸附位点），计算参杂Zn原子后的吸附能并且与纯Al吸附比较;

    (2)提交Al不同表面(001),(110),(111)的功函数计算。

  • 表面析氢动力学；

  • 基体、析出相电极电位差，金属第二相的表面能和吸附；

  • 自由能计算（考虑熵的影响）

第一天

常用软件

1
2
3
4
5
6
Atom: 文本编辑软件
VESTA： 原子结构可视化软件
FileZilla： 远程传输软件
Putty： Windows登陆集群软件
Material Studio： 原子建模软件
P4 vasp：vasp结果后处理软件
SSH命令远程登陆集群
常用的linux命令
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
clear
ls
cd <文件夹名>
cd ..
Tab 补全文件名称
cd ~
pwd
rm <文件名/文件夹名>
mkdir
vim <文件名> - i进入编辑状态，esc退出编辑状态，:wq退出文本编辑
cat POSCAR
cp /lustre/home/acct-umjzhh/umjzhh-3/BJUT/wyw/get.py ./
cp ../ly/get.py ./
grep "text" OUTCAR
head/tail -n 5 <文件名>
mv <文件名> 文件地址。eg： mv test ..    mv ./* bulk
cp ../../bulk/get.py ./
VASP计算的输入文件

POSCAR：各个原子的位置信息，可以自己手动输入，也可以从MP库中下载

INCAR：告诉VASP算什么，怎么算。各个标签的含义可以从vasp manuel中查找

KPOINTS：包含计算的K点信息(晶胞的边长*k=25)

POTCAR：对应的是计算过程中每个原子的平面波基组，描述原子中原子核和电子的相关信息。

VASP计算的输入文件

CONTCAR：弛豫后的结构，弛豫是物理学用语，指的是在某一个渐变物理过程中，从某一个状态逐渐地恢复到平衡态的过程。

WAVECAR：一般很大，存放波函数信息

CHGCAR：存放电荷信息

OSZICAR：存放能量信息

OUTCAR：比较综合的输出文件，很多有很多信息。

1
2
3
4
5
6
7
8
9
10
11
12
Al4 #代表了体系的名称，第一行随便写。 I'm from China
1.0 #缩放系数
4.038930 0.000000 0.000000
0.000000 4.038930 0.000000
0.000000 0.000000 4.038930
Al #体系中元素的名称
4  #体系中元素的个数
direct #可以写成分数坐标和直角（Cartersian）坐标，这里面写成分数坐标。
0.000000 0.000000 0.000000 Al
0.000000 0.500000 0.500000 Al
0.500000 0.000000 0.500000 Al
0.500000 0.500000 0.000000 Al
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
EDIFF = 1e-5 #能量的收敛标准，10**-5～10**-6 eV/atom
EDIFFG = -0.015 #力的收敛标准，暂时不考.BULK计算可以不放，slab计算要放。
ENCUT = 480 #截断能，grep “ENMAX” POTCAR，POTCAR一般使用PBE。
IBRION = 2 #conjugate-gradient algorithm
NSW = 200 #离子收敛步数
NELM = 100 #电子收敛步数
ISIF = 3 #controls whether the stress tensor is calculated.
ISMEAR = 1 #金属设置1就可以，半导体绝缘体设置0
ISPIN = 1 #电子自旋，Fe，Co，Ni或者O原子需要设置为2，其余为1

SIGMA = 0.05
POTIM = 0.1
LREAL = Auto
PREC = Accurate
ALGO = Fast

python脚本自动生成KPOINTS文件和POTCAR文件

Material project

https://materialsproject.org/（排在第一的绿色的），下载POSCAR。

https://www.vasp.at/forum/

https://cms.mpi.univie.ac.at/vasp/vasp/vasp.html

虚拟环境

1
2
3
4
5
6
7
conda env list #查看虚拟环境个数及名称
conda activate atomate_wyw #进入虚拟环境
source deactivate atomate_wyw #退出虚拟环境
python get.py #获得KPOINTS文件和POTCAR文件
Ctrl+D 退出集群登陆

scancel+id 取消任务

提交任务

1
2
3
sbatch vasp_std.slurm
squeue #查看任务计算情况
cat CONTCAR

第二天

表面能的计算公式

$E_{surf}=\frac{E_{slab}-E_{bulk}}{2S}$

其中$E_{surf}$代表体系对应的表面能，$E_{slab}$,$E_{bulk}$ 分别代表slab和bulk(体相)的能量，S代表表面积。

结果读取
1
2
3
4
5
6
# 判断体系是否正常结束？收敛？
OUTCAR, "reached required accuracy - stopping structural energy minimisation"
- grep "reached required accuracy" OUTCAR

# 体系的能量在哪里读？
OSZICAR 中E0后面的 E0= -.14946477E+02 = -14.946477 eV
切表面的过程

1.将计算好的CONTCAR保存为POSCAR,并在VESTA中将.vasp文件转位MS支持的.cif文件。File-Export data-save as .alcif

第三天

之前计算了Al(001), (110), (111)表面的表面能，并且发现Al(111)表面是最稳定的。

graph LR
ISIF=3优化晶格常数-->MS切表面
MS切表面-->计算低指数表面的表面能
计算低指数表面的表面能-->选择的最稳定表面做H吸附

吸附能的计算公式

$E_{ads}=E_{slabH}-E_{slab}-\frac{nE_{H_2}}{2}$

n代表吸附原子的个数

找到Al(111)表面具有代表性的吸附位点

在表面添加吸附H原子

编辑CONTCAR文件

功函数计算流程

$W_{F}=E_{vac}-E_{fermi}$

费米能的计算：VASP 弛豫计算结束之后，通过命令，grep 'E-fermi' OUTCAR 即可提取出来。

静电势能的计算: 通过在 INCAR 中添加： LVHAR =.TRUE. 这个参数。

加入这个参数，计算结束后，VASP 会输出一个文件： LOCPOT 文件。我们可以通过脚本，或者程序对这个文件后处理来获取静电势能。

功函数计算流程总结：

graph LR
优化获取稳定的结构-->将CONTCAR复制成POSCAR
将CONTCAR复制成POSCAR-->更改INCAR参数
更改INCAR参数-->提交任务
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
EDIFF = 1e-4
EDIFFG = -0.015
ENCUT = 480
IBRION = 2
# 计算功函数时需要更改！
NSW = 0
NELM = 100
ISIF = 3
ISMEAR = 1
ISPIN = 1

SIGMA = 0.05
POTIM = 0.1
LREAL = Auto
PREC = Accurate
ALGO = Fast

# 需要增加的参数
LVHAR = .TRUE.
IDIPOL =.TRUE.
IDIPOL = 3

计算$H_2$分子能量

基体、析出相电极电位差

能斯特方程

在标准状态下，标准H电极的电位是0.

常见的离子化学势网站：https://job-stiftung.de/index.php?id=63

Pourbaix diagram, E-PH图， 布拜图定腐蚀产物

我们假设的腐蚀环境，离子浓度为$10^{-6}$或者$10^{-8}$,环境为中性

NEB的计算

1.准备结构优化后的NEB初始状态和结束状态，并将两个文件夹改名为“start”，“end”。

2.在平行于begin文件夹的目录下输入：

nebmake.pl begin/CONTCAR end/CONTCAR 5

这行命令的意思是：利用nebmake脚本，以begin文件夹下的CONTCAR为初始结构，end文件夹下的CONTCAR为末态结构，在这个过度态过程中插入5个点。 nebmake的脚本需要自行下载编译。

3.将begin或end文件夹下的OUTCAR复制到00和06中，KPOINTS和POTCATR复制到00-06中。INCAR需要加入四个VASP相关的计算参数。

1
2
3
4
5
6
IBRION 改为3

LCLIMB = .TRUE.
IMAGES = 5
SPRING = -5.0
ICHAIN = 0

4.提交任务

sbatch -N 5 vasp_neb.slurm

5.结果处理脚本:

nebresultsrev.pl

化合物的表面建立过程

1.找出某一晶体方向的所有表面

2.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
from pymatgen.core.surface import Slab, SlabGenerator, generate_all_slabs, Structure, Lattice, ReconstructionGenerator, get_symmetrically_distinct_miller_indices
from pymatgen.core.structure import Structure
from pymatgen.analysis.adsorption import *
from pymatgen.ext.matproj import MPRester
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.io.vasp.inputs import Poscar
from matplotlib import pyplot as plt
import os

mpr = MPRester('Bw7HdCARiXvzEWJK')
path = '/Users/wangyaowei/Desktop/codes'

# struct = Structure.from_file("Mg.vasp")
struct = mpr.get_structure_by_material_id('mp-1124')
struct = SpacegroupAnalyzer(struct).get_conventional_standard_structure()

formula = struct.composition.reduced_formula
print(formula)
print("输入的结构是:{}".format(formula))

os.mkdir(formula)
os.chdir(formula)
ls = get_symmetrically_distinct_miller_indices(struct, 1, return_hkil=False)
print(ls)

for i in ls:
    slab = SlabGenerator(struct, miller_index=i, min_slab_size=25,min_vacuum_size=15.0, lll_reduce=True, center_slab=True)
    for n, slabs in enumerate(slab.get_slabs(bonds=None, ftol=0.1, tol=0.1, max_broken_bonds=0, repair=False, symmetrize=True)):
        slabs.make_supercell([[2,0,0],[0,2,0],[0,0,1]])
        name = str(i).split(',')[0][1]+str(i).split(',')[1][1]+str(i).split(',')[2][1]
        open(formula+'_'+name +'_' + str(n+1) + '.vasp', 'w').write(str(Poscar(slabs)))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
from pymatgen.io.vasp import Locpot
import matplotlib.pyplot as plt
import numpy as np

locpot = Locpot.from_file("LOCPOT") 

z_locpot = locpot.get_average_along_axis(2)

fermi_level =  1.7968

import matplotlib.pyplot as plt
import numpy as np
import matplotlib.pyplot as plt 
from pylab import * #控制刻度线
from matplotlib import cm
from matplotlib.collections import LineCollection
cm = plt.cm.get_cmap('viridis')

plt.figure(figsize=(11.3,6))#图片大小设置
# plt.ylim(-2.54+2.534,-2.44+2.534)
label_font = {'family':'Arial','weight':'normal','size':22}#设置标签和图例字体
legend_font = {'family':'Arial','weight':'normal','size':12}
plt.ylabel('Electrostatic potential',label_font)
plt.xlabel('Distance along Z direction',label_font)
plt.yticks(fontproperties = 'Arial', size = 24)
plt.xticks(fontproperties = 'Arial', size = 24)
# plt.xticks(np.arange(1, 5.1, 1))#横坐标范围及间隔
matplotlib.rcParams['patch.linewidth'] = 1.5#图例边框线宽度
plt.rcParams['xtick.direction'] = 'in'#将x周的刻度线方向设置向内
plt.rcParams['ytick.direction'] = 'in'#将y轴的刻度方向设置向内
plt.tick_params(which='major',width=2,length=6)#设置刻度线宽度
bwith = 2#加粗图四周边框
ax=plt.gca()
ax.spines['bottom'].set_linewidth(bwith)    
ax.spines['top'].set_linewidth(bwith)
ax.spines['left'].set_linewidth(bwith)  
ax.spines['right'].set_linewidth(bwith) 

plt.plot(z_locpot, lw = 3)
ax.axvline(300,color='red',ls='dashed')
ax.axvline(352,color='red',ls='dashed')

vac = np.average(z_locpot[300:352])

plt.savefig('plotential.tif',dpi=200)

print(vac)

print("Al(001)表面的功函数为：{}".format(round(vac-fermi_level,3)))